Updated cluster script

- added energy resolution.

clustering/scripts/barrel_clusters.cxx

@@ -49,6 +49,20 @@ auto eta = [](ROOT::VecOps::RVec<dd4pod::Geant4ParticleData> const& in) {
   return result;
 };
 
+auto delta_E_over_E = [](std::vector<ROOT::Math::PxPyPzMVector> const& thrown, const std::vector<eic::ClusterData>& clusters) {
+  std::vector<double> result;
+  double best = 1000000.0;
+  for (const auto& p : thrown) {
+    for (const auto& c : clusters) {
+      double dE = (p.E() - c.energy)/p.E();
+      result.push_back(dE );
+      if( dE < best) {
+        best = dE;
+      }
+    }
+  }
+  return result;
+};
 auto delta_E = [](std::vector<ROOT::Math::PxPyPzMVector> const& thrown, const std::vector<eic::ClusterData>& clusters) {
   std::vector<double> result;
   double best = 1000000.0;
@@ -77,6 +91,7 @@ void barrel_clusters(const char* in_fname = "topside/rec_barrel_clusters.root")
                 .Define("thrownTheta", theta, {"thrownP"})
                 .Define("thrownMomentum", momentum, {"thrownP"})
                 .Define("delta_E", delta_E, {"thrownP","SimpleClusters"})
+                .Define("delta_E_over_E", delta_E_over_E, {"thrownP","SimpleClusters"})
                 .Define("nclusters", "SimpleClusters.size()")
                 .Define("Ecluster",
                         [](const std::vector<eic::ClusterData>& in) {
@@ -99,7 +114,9 @@ void barrel_clusters(const char* in_fname = "topside/rec_barrel_clusters.root")
   auto h_Ecluster1        = d1.Histo1D({"h_Ecluster1", "One cluster events;  cluster E [GeV]", 100, 0, 30}, "Ecluster");
   auto h_momentum_thrown1 = d1.Histo1D({"h_momentum_thrown", "; E [GeV]", 100, 0.0, 30.0},"thrownMomentum");
 
-  auto h_delta_E        = d0.Histo1D({"h_delta_E", ";  #Delta E [GeV]", 100, -10, 10}, "delta_E");
+  auto h_delta_E = d0.Histo1D({"h_delta_E", ";  #Delta E [GeV]", 100, -3, 3}, "delta_E");
+  auto h_delta_E_over_E =
+      d0.Histo1D({"h_delta_E_over_E", ";  #frac{E_{thrown}-E_{cluster}}{E_{thrown}} ", 100, -1, 1}, "delta_E_over_E");
 
   //auto h_Ecluster2 = d1.Filter("thrownMomentum > 4").Histo1D({"h_Ecluster1", "One cluster events;  cluster E [GeV]",100, 0,30},"Ecluster");
   //auto h_momentum_thrown2 = d1.Filter("thrownMomentum > 4").Histo1D({"h_momentum_thrown", "; E [GeV]", 100, 0.0, 30.0},"thrownMomentum");
@@ -128,6 +145,9 @@ void barrel_clusters(const char* in_fname = "topside/rec_barrel_clusters.root")
   h_delta_E->DrawCopy();
   c->SaveAs("results/clustering/barrel_clusters_delta_E.png");
 
+  h_delta_E_over_E->DrawCopy();
+  c->SaveAs("results/clustering/barrel_clusters_delta_E_over_E.png");
+
   std::cout <<  (*c_nclusters1) << " / " << (*c_thrown) << " = single cluster events/all \n";
 
   //c->SaveAs("results/crystal_cal_electrons_Ecluster.png");

tracking/central_electrons.sh

@@ -28,14 +28,17 @@ mkdir ${JUGGLER_DETECTOR}/build
 pushd ${JUGGLER_DETECTOR}/build
 cmake ../. -DCMAKE_INSTALL_PREFIX=/usr/local && make -j30 install
 popd
-#
+
 ## generate the input events
-## note datasets is now only used to develop datasets.
-##git clone https://eicweb.phy.anl.gov/EIC/datasets.git datasets
 root -b -q "tracking/scripts/gen_central_electrons.cxx(${JUGGLER_N_EVENTS}, \"${JUGGLER_FILE_NAME_TAG}.hepmc\")"
+if [[ "$?" -ne "0" ]] ; then
+  echo "ERROR running script"
+  exit 1
+fi
+
 #
 pushd ${JUGGLER_DETECTOR}
-#ls -l
+
 ## run geant4 simulations
 npsim --runType batch \
       --part.minimalKineticEnergy 1000*GeV  \
@@ -47,6 +50,11 @@ npsim --runType batch \
 
 # Need to figure out how to pass file name to juggler from the commandline
 xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv gaudirun.py ../tracking/options/tracker_reconstruction.py
+
+if [[ "$?" -ne "0" ]] ; then
+  echo "ERROR running juggler"
+  exit 1
+fi
 ls -l
 popd
 

tracking/options/tracker_reconstruction.py

@@ -2,6 +2,7 @@ from Gaudi.Configuration import *
 
 from GaudiKernel.DataObjectHandleBase import DataObjectHandleBase
 from Configurables import ApplicationMgr, EICDataSvc, PodioOutput, GeoSvc
+from GaudiKernel import SystemOfUnits as units
 
 detector_name = "topside"
 if "JUGGLER_DETECTOR" in os.environ :
@@ -16,32 +17,84 @@ geo_service  = GeoSvc("GeoSvc", detectors=["{}.xml".format(detector_name)])
 podioevent   = EICDataSvc("EventDataSvc", inputs=[input_sim_file], OutputLevel=DEBUG)
 
 from Configurables import PodioInput
+from Configurables import Jug__Base__InputCopier_dd4pod__Geant4ParticleCollection_dd4pod__Geant4ParticleCollection_ as MCCopier
+from Configurables import Jug__Base__InputCopier_dd4pod__CalorimeterHitCollection_dd4pod__CalorimeterHitCollection_ as CalCopier
+from Configurables import Jug__Base__InputCopier_dd4pod__TrackerHitCollection_dd4pod__TrackerHitCollection_ as TrkCopier
+
 from Configurables import Jug__Digi__ExampleCaloDigi as ExampleCaloDigi
 from Configurables import Jug__Digi__UFSDTrackerDigi as UFSDTrackerDigi
+from Configurables import Jug__Digi__EMCalorimeterDigi as EMCalorimeterDigi
+
 from Configurables import Jug__Reco__TrackerHitReconstruction as TrackerHitReconstruction
 from Configurables import Jug__Reco__TrackerSourceLinker as TrackerSourceLinker
 from Configurables import Jug__Reco__TrackParamTruthInit as TrackParamTruthInit
+from Configurables import Jug__Reco__TrackParamClusterInit as TrackParamClusterInit
 from Configurables import Jug__Reco__TrackFindingAlgorithm as TrackFindingAlgorithm
 from Configurables import Jug__Reco__ParticlesFromTrackFit as ParticlesFromTrackFit
-from Configurables import Jug__Base__InputCopier_dd4pod__Geant4ParticleCollection_dd4pod__Geant4ParticleCollection_ as MCCopier
-from Configurables import Jug__Base__InputCopier_dd4pod__CalorimeterHitCollection_dd4pod__CalorimeterHitCollection_ as CalCopier
-from Configurables import Jug__Base__InputCopier_dd4pod__TrackerHitCollection_dd4pod__TrackerHitCollection_ as TrkCopier
+from Configurables import Jug__Reco__EMCalReconstruction as EMCalReconstruction
+
+from Configurables import Jug__Reco__SimpleClustering as SimpleClustering
+
 
 
 podioinput = PodioInput("PodioReader", 
-                        collections=["mcparticles","SiTrackerBarrelHits"])#, OutputLevel=DEBUG)
+                        collections=["mcparticles","SiTrackerBarrelHits","EcalBarrelHits"])#, OutputLevel=DEBUG)
 
 ## copiers to get around input --> output copy bug. Note the "2" appended to the output collection.
 copier = MCCopier("MCCopier", inputCollection="mcparticles", outputCollection="mcparticles2",OutputLevel=DEBUG) 
 trkcopier = TrkCopier("TrkCopier", inputCollection="SiTrackerBarrelHits", outputCollection="SiTrackerBarrelHits2",OutputLevel=DEBUG) 
 #caldigi = ExampleCaloDigi(inputHitCollection="FAEC_ShHits",outputHitCollection="RawFAECShowerHits")
 #ufsd_digi = UFSDTrackerDigi("ufsd_digi", inputHitCollection="SiVertexBarrelHits",outputHitCollection="VertexRawHits",timeResolution=8)
-ufsd_digi = UFSDTrackerDigi("ufsd_digi", inputHitCollection="SiTrackerBarrelHits",outputHitCollection="SiTrackerBarrelRawHits",timeResolution=8)
-trackpartruth = TrackParamTruthInit("trk_par_init",inputMCParticles="mcparticles",outputInitialTrackParameters="InitTrackParams",OutputLevel=DEBUG)
-trackerhit = TrackerHitReconstruction("trk_hit_reco",inputHitCollection="SiTrackerBarrelRawHits",outputHitCollection="TrackerBarrelRecHits",OutputLevel=DEBUG)
-sourcelinker = TrackerSourceLinker("trk_srclinker",inputHitCollection="TrackerBarrelRecHits",outputSourceLinks="BarrelTrackSourceLinks",OutputLevel=DEBUG)
-trk_find_alg = TrackFindingAlgorithm("trk_find_alg",inputSourceLinks="BarrelTrackSourceLinks",inputInitialTrackParameters= "InitTrackParams", outputTrajectories="trajectories",OutputLevel=DEBUG)
-parts_from_fit = ParticlesFromTrackFit("parts_from_fit",inputTrajectories="trajectories",outputParticles="ReconstructedParticles",OutputLevel=DEBUG)
+
+ecal_digi = EMCalorimeterDigi("ecal_digi", 
+        inputHitCollection="EcalBarrelHits", 
+        outputHitCollection="RawEcalBarrelHits")
+
+ecal_reco = EMCalReconstruction("ecal_reco", 
+        inputHitCollection="RawEcalBarrelHits", 
+        outputHitCollection="RecEcalBarrelHits",
+        minModuleEdep=1.0*units.MeV,OutputLevel=DEBUG)
+
+simple_cluster = SimpleClustering("simple_cluster", 
+        inputHitCollection="RecEcalBarrelHits", 
+        outputClusters="SimpleClusters",
+        minModuleEdep=1.0*units.MeV,OutputLevel=DEBUG)
+
+cluster_init = TrackParamClusterInit("trk_clust_init",
+        inputClusters="SimpleClusters",
+        outputInitialTrackParameters="InitTrackParamsFromClusters",
+        OutputLevel=DEBUG)
+
+ufsd_digi = UFSDTrackerDigi("ufsd_digi", 
+        inputHitCollection="SiTrackerBarrelHits",
+        outputHitCollection="SiTrackerBarrelRawHits",
+        timeResolution=8)
+
+trackpartruth = TrackParamTruthInit("trk_par_init",
+        inputMCParticles="mcparticles",
+        outputInitialTrackParameters="InitTrackParams",
+        OutputLevel=DEBUG)
+
+trackerhit = TrackerHitReconstruction("trk_hit_reco",
+        inputHitCollection="SiTrackerBarrelRawHits",
+        outputHitCollection="TrackerBarrelRecHits",
+        OutputLevel=DEBUG)
+
+sourcelinker = TrackerSourceLinker("trk_srclinker",
+        inputHitCollection="TrackerBarrelRecHits",
+        outputSourceLinks="BarrelTrackSourceLinks",
+        OutputLevel=DEBUG)
+
+trk_find_alg = TrackFindingAlgorithm("trk_find_alg",
+        inputSourceLinks="BarrelTrackSourceLinks",
+        inputInitialTrackParameters= "InitTrackParams",#"InitTrackParamsFromClusters", 
+        outputTrajectories="trajectories",
+        OutputLevel=DEBUG)
+
+parts_from_fit = ParticlesFromTrackFit("parts_from_fit",
+        inputTrajectories="trajectories",
+        outputParticles="ReconstructedParticles",
+        OutputLevel=DEBUG)
 
 #types = []
 ## this printout is useful to check that the type information is passed to python correctly
@@ -73,8 +126,11 @@ out.outputCommands = ["keep *",
 ApplicationMgr(
     TopAlg = [podioinput, 
               #copier, trkcopier,
+              ecal_digi,ecal_reco, 
+              simple_cluster,
               ufsd_digi, trackerhit, 
-              sourcelinker, trackpartruth, 
+              sourcelinker, 
+              cluster_init, trackpartruth, 
               trk_find_alg, parts_from_fit,
               out
               ],

tracking/scripts/gen_central_electrons.cxx

@@ -17,8 +17,8 @@ using namespace HepMC3;
 void gen_central_electrons(int n_events = 100, 
                      const char* out_fname = "central_electrons.hepmc")
 {
-  double cos_theta_min = std::cos( 70.0*(M_PI/180.0));
-  double cos_theta_max = std::cos(110.0*(M_PI/180.0));
+  double cos_theta_min = std::cos( 50.0*(M_PI/180.0));
+  double cos_theta_max = std::cos(130.0*(M_PI/180.0));
 
   WriterAscii hepmc_output(out_fname);
   int events_parsed = 0;
@@ -39,7 +39,7 @@ void gen_central_electrons(int n_events = 100,
         FourVector(0.0, 0.0, 0.0, 0.938), 2212, 4);
 
     // Define momentum
-    Double_t p     = r1->Uniform(1.0, 10.0);
+    Double_t p     = r1->Uniform(9.0, 10.0);
     Double_t phi   = r1->Uniform(0.0, 2.0 * M_PI);
     Double_t costh = r1->Uniform(cos_theta_min, cos_theta_max);
     Double_t th    = std::acos(costh);