diff --git a/clustering/scripts/barrel_clusters.cxx b/clustering/scripts/barrel_clusters.cxx index 2cfb8783cff975d8d369e3904a79bbf9ac7dd578..b756a970021a7627adb1e887112728bfbdb46b68 100644 --- a/clustering/scripts/barrel_clusters.cxx +++ b/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"); diff --git a/tracking/central_electrons.sh b/tracking/central_electrons.sh index 8fb2dadef829eeab71fc42ff91ae8fb307b28f56..376135cb8907357cb9fa1aee6c09489159e40316 100644 --- a/tracking/central_electrons.sh +++ b/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 diff --git a/tracking/options/tracker_reconstruction.py b/tracking/options/tracker_reconstruction.py index 68988a17d46d3e4c5408133cb5662df3f5984591..12260e4f30756ae23cd6b7d10f96908bbe2df46b 100644 --- a/tracking/options/tracker_reconstruction.py +++ b/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 ], diff --git a/tracking/scripts/gen_central_electrons.cxx b/tracking/scripts/gen_central_electrons.cxx index 06867dfaea7507522bf41a0ffc3ca87d413fb301..6629d85cd1337c4a4794ae4f3f1a1fca9203bcc8 100644 --- a/tracking/scripts/gen_central_electrons.cxx +++ b/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);