diff --git a/Examples/options/tut_detsim_SDT.py b/Examples/options/tut_detsim_SDT.py index 027d581388a4aa46983f6d34629364e4c694ce80..ca5498b772ef7afdf10146dc2a551cd4c751a277 100644 --- a/Examples/options/tut_detsim_SDT.py +++ b/Examples/options/tut_detsim_SDT.py @@ -121,6 +121,24 @@ detsimalg.RootDetElem = "WorldDetElemTool" from Configurables import AnExampleDetElemTool example_dettool = AnExampleDetElemTool("AnExampleDetElemTool") +from Configurables import CalorimeterSensDetTool +from Configurables import DriftChamberSensDetTool + +calo_sensdettool = CalorimeterSensDetTool("CalorimeterSensDetTool") +driftchamber_sensdettool = DriftChamberSensDetTool("DriftChamberSensDetTool") + +# dedxoption = "DummyDedxSimTool" +dedxoption = "BetheBlochEquationDedxSimTool" + +driftchamber_sensdettool.DedxSimTool = dedxoption + +from Configurables import DummyDedxSimTool +from Configurables import BetheBlochEquationDedxSimTool + +if dedxoption == "DummyDedxSimTool": + dedx_simtool = DummyDedxSimTool("DummyDedxSimTool") +elif dedxoption == "BetheBlochEquationDedxSimTool": + dedx_simtool = BetheBlochEquationDedxSimTool("BetheBlochEquationDedxSimTool") ############################################################################## # POD I/O diff --git a/Examples/options/tut_detsim_pan_matrix.py b/Examples/options/tut_detsim_pan_matrix.py index 1d5bd4561d3f5bd4318455d7c6a9a54df88c5aa2..ab3b6cbb1499b30cca09ce964ecfadbbb6205032 100644 --- a/Examples/options/tut_detsim_pan_matrix.py +++ b/Examples/options/tut_detsim_pan_matrix.py @@ -3,7 +3,6 @@ import os print(os.environ["DD4HEP_LIBRARY_PATH"]) import sys -# sys.exit(0) from Gaudi.Configuration import * @@ -12,20 +11,13 @@ from Gaudi.Configuration import * ############################################################################## from Configurables import RndmGenSvc, HepRndm__Engine_CLHEP__RanluxEngine_ -# rndmengine = HepRndm__Engine_CLHEP__RanluxEngine_() # The default engine in Gaudi rndmengine = HepRndm__Engine_CLHEP__HepJamesRandom_() # The default engine in Geant4 rndmengine.SetSingleton = True rndmengine.Seeds = [42] -# rndmgensvc = RndmGenSvc("RndmGenSvc") -# rndmgensvc.Engine = rndmengine.name() - - ############################################################################## # Event Data Svc ############################################################################## -#from Configurables import CEPCDataSvc -#dsvc = CEPCDataSvc("EventDataSvc") from Configurables import K4DataSvc dsvc = K4DataSvc("EventDataSvc") @@ -34,7 +26,6 @@ dsvc = K4DataSvc("EventDataSvc") # Geometry Svc ############################################################################## -# geometry_option = "CepC_v4-onlyTracker.xml" geometry_option = "CepC_v4-onlyVXD.xml" if not os.getenv("DETCEPCV4ROOT"): @@ -49,7 +40,7 @@ if not os.path.exists(geometry_path): from Configurables import GeoSvc geosvc = GeoSvc("GeoSvc") #geosvc.compact = geometry_path -geosvc.compact = "../Detector/DetEcalMatrix/compact/det.xml" +geosvc.compact = "./Detector/DetEcalMatrix/compact/det.xml" ############################################################################## # Physics Generator @@ -63,18 +54,16 @@ from Configurables import GenPrinter gun = GtGunTool("GtGunTool") gun.Particles = ["gamma","gamma"] -gun.EnergyMins= [10, 10] # GeV -gun.EnergyMaxs= [10, 10] # GeV +gun.EnergyMins= [5 , 10] # GeV +gun.EnergyMaxs= [5 , 10] # GeV gun.ThetaMins = [90, 90] # degree gun.ThetaMaxs = [90, 90] # degree -gun.PhiMins = [0, 1 ] # degree -gun.PhiMaxs = [0, 1 ] # degree +gun.PhiMins = [0, 4 ] # degree +gun.PhiMaxs = [0, 4 ] # degree -stdheprdr = StdHepRdr("StdHepRdr") -#stdheprdr.Input = "/cefs/data/stdhep/CEPC250/2fermions/E250.Pbhabha.e0.p0.whizard195/bhabha.e0.p0.00001.stdhep" -#stdheprdr.Input = "/cefs/data/stdhep/CEPC250/2fermions/E250.Pbhabha.e0.p0.whizard195/bhabha.e0.p0.00001.stdhep" -stdheprdr.Input = "/cefs/data/stdhep/CEPC250/higgs/E250.Pbbh.whizard195/E250.Pbbh_X.e0.p0.whizard195/Pbbh_X.e0.p0.00001.stdhep" +#stdheprdr = StdHepRdr("StdHepRdr") +#stdheprdr.Input = "/cefs/data/stdhep/CEPC250/higgs/E250.Pbbh.whizard195/E250.Pbbh_X.e0.p0.whizard195/Pbbh_X.e0.p0.00001.stdhep" # lciordr = SLCIORdr("SLCIORdr") # lciordr.Input = "/cefs/data/stdhep/lcio250/signal/Higgs/E250.Pbbh.whizard195/E250.Pbbh_X.e0.p0.whizard195/Pbbh_X.e0.p0.00001.slcio" @@ -95,31 +84,21 @@ genalg.GenTools = ["GtGunTool"] # Detector Simulation ############################################################################## from Configurables import DetSimSvc - detsimsvc = DetSimSvc("DetSimSvc") - -# from Configurables import ExampleAnaElemTool -# example_anatool = ExampleAnaElemTool("ExampleAnaElemTool") - from Configurables import DetSimAlg - detsimalg = DetSimAlg("DetSimAlg") - -# detsimalg.VisMacs = ["vis.mac"] - +detsimalg.VisMacs = ["Examples/options/vis.mac"] detsimalg.RunCmds = [ # "/tracking/verbose 1", ] detsimalg.AnaElems = [ - # example_anatool.name() - # "ExampleAnaElemTool" "Edm4hepWriterAnaElemTool" ] detsimalg.RootDetElem = "WorldDetElemTool" - from Configurables import AnExampleDetElemTool example_dettool = AnExampleDetElemTool("AnExampleDetElemTool") - +############################################################################## +# Detector digitization ############################################################################## from Configurables import CaloDigiAlg example_CaloDigiAlg = CaloDigiAlg("CaloDigiAlg") @@ -130,39 +109,23 @@ example_CaloDigiAlg.CaloAssociationCollection = "RecoCaloAssociation_ECALBarr ############################################################################## from Configurables import GearSvc gearSvc = GearSvc("GearSvc") -#gearSvc.GearXMLFile = "/junofs/users/wxfang/CEPC/CEPCOFF/doSim/fullDet/GearOutput.xml" -gearSvc.GearXMLFile = "../Detector/DetCEPCv4/compact/FullDetGear.xml" +gearSvc.GearXMLFile = "./Detector/DetCEPCv4/compact/FullDetGear.xml" +############################################################################## +# Pandora ############################################################################## -#from Configurables import PandoraPFAlg from Configurables import PandoraMatrixAlg -#pandoralg = PandoraPFAlg("PandoraPFAlg") pandoralg = PandoraMatrixAlg("PandoraMatrixAlg") -## KEEP same with lcioinput name for the ReadXXX ########### -pandoralg.ReadMCParticle = "MCParticle" -pandoralg.ReadECALBarrel = "ECALBarrel" -pandoralg.ReadECALEndcap = "ECALEndcap" -pandoralg.ReadECALOther = "ECALOther" -pandoralg.ReadHCALBarrel = "HCALBarrel" -pandoralg.ReadHCALEndcap = "HCALEndcap" -pandoralg.ReadHCALOther = "HCALOther" -pandoralg.ReadMUON = "MUON" -pandoralg.ReadLCAL = "LCAL" -pandoralg.ReadLHCAL = "LHCAL" -pandoralg.ReadBCAL = "BCAL" -pandoralg.ReadKinkVertices = "KinkVertices" -pandoralg.ReadProngVertices = "ProngVertices" -pandoralg.ReadSplitVertices = "SplitVertices" -pandoralg.ReadV0Vertices = "V0Vertices" -pandoralg.ReadTracks = "MarlinTrkTracks" -pandoralg.MCRecoCaloAssociation = "RecoCaloAssociation_ECALBarrel" +pandoralg.collections = [ + "MCParticle:MCParticle", + "CalorimeterHit:ECALBarrel", + "MCRecoCaloAssociation:RecoCaloAssociation_ECALBarrel" + ] pandoralg.WriteClusterCollection = "PandoraClusters" pandoralg.WriteReconstructedParticleCollection = "PandoraPFOs" pandoralg.WriteVertexCollection = "PandoraPFANewStartVertices" -pandoralg.AnaOutput = "/cefs/higgs/wxfang/cepc/Pandora/Ana/gamma/Ana_gamma_Matrix_Rel_10GeV_test.root" - -pandoralg.PandoraSettingsDefault_xml = "/junofs/users/wxfang/MyGit/MarlinPandora/scripts/PandoraSettingsDefault_wx.xml" -#### Do not chage the collection name, only add or delete ############### +pandoralg.AnaOutput = "AnaMatrix.root" +pandoralg.PandoraSettingsDefault_xml = "./Reconstruction/PFA/Pandora/PandoraSettingsDefault.xml" pandoralg.TrackCollections = ["MarlinTrkTracks"] pandoralg.ECalCaloHitCollections= ["ECALBarrel", "ECALEndcap", "ECALOther"] pandoralg.HCalCaloHitCollections= ["HCALBarrel", "HCALEndcap", "HCALOther"] @@ -170,32 +133,31 @@ pandoralg.LCalCaloHitCollections= ["LCAL"] pandoralg.LHCalCaloHitCollections= ["LHCAL"] pandoralg.MuonCaloHitCollections= ["MUON"] pandoralg.MCParticleCollections = ["MCParticle"] -pandoralg.RelCaloHitCollections = ["RecoCaloAssociation_ECALBarrel", "RecoCaloAssociation_ECALEndcap", "RecoCaloAssociation_ECALOther", "RecoCaloAssociation_HCALBarrel", "RecoCaloAssociation_HCALEndcap", "RecoCaloAssociation_HCALOther", "RecoCaloAssociation_LCAL", "RecoCaloAssociation_LHCAL", "RecoCaloAssociation_MUON"] +pandoralg.RelCaloHitCollections = ["RecoCaloAssociation_ECALBarrel"] pandoralg.RelTrackCollections = ["MarlinTrkTracksMCTruthLink"] pandoralg.KinkVertexCollections = ["KinkVertices"] pandoralg.ProngVertexCollections= ["ProngVertices"] pandoralg.SplitVertexCollections= ["SplitVertices"] pandoralg.V0VertexCollections = ["V0Vertices"] -pandoralg.ECalToMipCalibration = 112 #1000MeV/8.918 #CEPC :160.0 +pandoralg.ECalToMipCalibration = 112 #1000MeV/8.918 pandoralg.HCalToMipCalibration = 34.8 -pandoralg.ECalMipThreshold = 0.225# 8.918*0.225=2.00655 #CEPC 0.5 +pandoralg.ECalMipThreshold = 0.225# 8.918*0.225=2.00655 pandoralg.HCalMipThreshold = 0.3 -pandoralg.ECalToEMGeVCalibration= 1.# BGO, to be tuned # CEPC: 0.9 for G2CD Digi, 1.007 for NewLDCaloDigi +pandoralg.ECalToEMGeVCalibration= 1.# BGO, to be tuned pandoralg.HCalToEMGeVCalibration= 1.007 -pandoralg.ECalToHadGeVCalibrationBarrel= 1.12 #very small effect +pandoralg.ECalToHadGeVCalibrationBarrel= 1.12 pandoralg.ECalToHadGeVCalibrationEndCap= 1.12 pandoralg.HCalToHadGeVCalibration= 1.07 pandoralg.MuonToMipCalibration= 10.0 pandoralg.DigitalMuonHits= 0 pandoralg.MaxHCalHitHadronicEnergy = 1.0 pandoralg.UseOldTrackStateCalculation= 0 -pandoralg.AbsorberRadLengthECal= 0.08945 #BG0: 1/11.18 mm , CEPC: 0.2854 = 1/3.504 mm -pandoralg.AbsorberIntLengthECal= 0.00448 #BG0: 1/223.2 mm , CEPC: 0.0101 = 1/99.46 mm +pandoralg.AbsorberRadLengthECal= 0.08945 #BG0: 1/11.18 mm +pandoralg.AbsorberIntLengthECal= 0.00448 #BG0: 1/223.2 mm pandoralg.AbsorberRadLengthHCal= 0.0569 pandoralg.AbsorberIntLengthHCal= 0.006 pandoralg.AbsorberRadLengthOther= 0.0569 pandoralg.AbsorberIntLengthOther= 0.006 - ############################################################################## # write PODIO file @@ -207,13 +169,10 @@ write.outputCommands = ["keep *"] # ApplicationMgr from Configurables import ApplicationMgr ApplicationMgr( - TopAlg = [genalg, detsimalg, example_CaloDigiAlg, pandoralg], - #TopAlg = [genalg, detsimalg], - #TopAlg = [genalg, detsimalg, write], - #TopAlg = [read, pandoralg], + TopAlg = [genalg, detsimalg], + #TopAlg = [genalg, detsimalg, example_CaloDigiAlg, pandoralg], EvtSel = 'NONE', - EvtMax = 10, - #ExtSvc = [rndmengine, dsvc, geosvc, gearSvc], + EvtMax = 50, ExtSvc = [rndmengine, dsvc, geosvc, gearSvc,detsimsvc], OutputLevel=INFO ) diff --git a/Examples/options/tut_detsim_pandora.py b/Examples/options/tut_detsim_pandora.py index e1ac4a7fa9053fac9cd929609d125fac1ca5f727..be552eaaf41ce2bfc393043ff924ef7ce067b84a 100644 --- a/Examples/options/tut_detsim_pandora.py +++ b/Examples/options/tut_detsim_pandora.py @@ -3,7 +3,6 @@ import os print(os.environ["DD4HEP_LIBRARY_PATH"]) import sys -# sys.exit(0) from Gaudi.Configuration import * @@ -17,15 +16,9 @@ rndmengine = HepRndm__Engine_CLHEP__HepJamesRandom_() # The default engine in Ge rndmengine.SetSingleton = True rndmengine.Seeds = [42] -# rndmgensvc = RndmGenSvc("RndmGenSvc") -# rndmgensvc.Engine = rndmengine.name() - - ############################################################################## # Event Data Svc ############################################################################## -#from Configurables import CEPCDataSvc -#dsvc = CEPCDataSvc("EventDataSvc") from Configurables import K4DataSvc dsvc = K4DataSvc("EventDataSvc") @@ -34,8 +27,7 @@ dsvc = K4DataSvc("EventDataSvc") # Geometry Svc ############################################################################## -# geometry_option = "CepC_v4-onlyTracker.xml" -geometry_option = "CepC_v4-onlyVXD.xml" +geometry_option = "CepC_v4-onlyECAL.xml" if not os.getenv("DETCEPCV4ROOT"): print("Can't find the geometry. Please setup envvar DETCEPCV4ROOT." ) @@ -70,10 +62,8 @@ gun.PhiMins = [0, 1 ] # degree gun.PhiMaxs = [0, 1 ] # degree -stdheprdr = StdHepRdr("StdHepRdr") -#stdheprdr.Input = "/cefs/data/stdhep/CEPC250/2fermions/E250.Pbhabha.e0.p0.whizard195/bhabha.e0.p0.00001.stdhep" -#stdheprdr.Input = "/cefs/data/stdhep/CEPC250/2fermions/E250.Pbhabha.e0.p0.whizard195/bhabha.e0.p0.00001.stdhep" -stdheprdr.Input = "/cefs/data/stdhep/CEPC250/higgs/E250.Pbbh.whizard195/E250.Pbbh_X.e0.p0.whizard195/Pbbh_X.e0.p0.00001.stdhep" +#stdheprdr = StdHepRdr("StdHepRdr") +#stdheprdr.Input = "/cefs/data/stdhep/CEPC250/higgs/E250.Pbbh.whizard195/E250.Pbbh_X.e0.p0.whizard195/Pbbh_X.e0.p0.00001.stdhep" # lciordr = SLCIORdr("SLCIORdr") # lciordr.Input = "/cefs/data/stdhep/lcio250/signal/Higgs/E250.Pbbh.whizard195/E250.Pbbh_X.e0.p0.whizard195/Pbbh_X.e0.p0.00001.slcio" @@ -97,9 +87,6 @@ from Configurables import DetSimSvc detsimsvc = DetSimSvc("DetSimSvc") -# from Configurables import ExampleAnaElemTool -# example_anatool = ExampleAnaElemTool("ExampleAnaElemTool") - from Configurables import DetSimAlg detsimalg = DetSimAlg("DetSimAlg") @@ -110,8 +97,6 @@ detsimalg.RunCmds = [ # "/tracking/verbose 1", ] detsimalg.AnaElems = [ - # example_anatool.name() - # "ExampleAnaElemTool" "Edm4hepWriterAnaElemTool" ] detsimalg.RootDetElem = "WorldDetElemTool" diff --git a/Reconstruction/PFA/Pandora/GaudiPandora/include/PandoraPFAlg.h b/Reconstruction/PFA/Pandora/GaudiPandora/include/PandoraPFAlg.h index e357bd34b55532325758d9f45ac9ea0dbfb34868..8188fb965dab87b4a2cdefcae3be7c24d7bccea2 100644 --- a/Reconstruction/PFA/Pandora/GaudiPandora/include/PandoraPFAlg.h +++ b/Reconstruction/PFA/Pandora/GaudiPandora/include/PandoraPFAlg.h @@ -128,7 +128,7 @@ protected: - Gaudi::Property< std::string > m_PandoraSettingsXmlFile { this, "PandoraSettingsDefault_xml", "/junofs/users/wxfang/MyGit/MarlinPandora/scripts/PandoraSettingsDefault_wx.xml" }; + Gaudi::Property< std::string > m_PandoraSettingsXmlFile { this, "PandoraSettingsDefault_xml", "PandoraSettingsDefault_wx.xml" }; Gaudi::Property<int> m_NEventsToSkip { this, "NEventsToSkip", 0 }; Gaudi::Property< std::vector<std::string> > m_TrackCollections{ this, "TrackCollections", {"MarlinTrkTracks"} }; diff --git a/Reconstruction/PFA/Pandora/GaudiPandora/src/CaloHitCreator.cpp b/Reconstruction/PFA/Pandora/GaudiPandora/src/CaloHitCreator.cpp index 12120d3f56e3e8b2e8de9a03ddd0ceeb90dab51d..e1c24b1e182b76ab503e007a0ec2d402424fe117 100644 --- a/Reconstruction/PFA/Pandora/GaudiPandora/src/CaloHitCreator.cpp +++ b/Reconstruction/PFA/Pandora/GaudiPandora/src/CaloHitCreator.cpp @@ -193,7 +193,7 @@ pandora::StatusCode CaloHitCreator::CreateECalCaloHits(const CollectionMaps& col } //caloHitParameters.m_mipEquivalentEnergy = pCaloHit->getEnergy() * eCalToMip * absorberCorrection; - caloHitParameters.m_mipEquivalentEnergy = pCaloHit->getEnergy() * eCalToMip;//FIXME. is absorberCorrection it needed for digi input + caloHitParameters.m_mipEquivalentEnergy = pCaloHit->getEnergy() * eCalToMip;//is absorberCorrection it needed for digi input, also the m_mipEquivalentEnergy seems is not used in alg if (caloHitParameters.m_mipEquivalentEnergy.Get() < eCalMipThreshold) continue; @@ -266,7 +266,7 @@ pandora::StatusCode CaloHitCreator::CreateHCalCaloHits(const CollectionMaps& col PandoraApi::CaloHit::Parameters caloHitParameters; caloHitParameters.m_hitType = pandora::HCAL; - caloHitParameters.m_isDigital = false; + caloHitParameters.m_isDigital = false;// if it is DHCAL or AHCAL caloHitParameters.m_layer = cellIdDecoder(pCaloHit)[layerCoding.c_str()]; caloHitParameters.m_isInOuterSamplingLayer = (this->GetNLayersFromEdge(pCaloHit) <= m_settings.m_nOuterSamplingLayers); this->GetCommonCaloHitProperties(pCaloHit, caloHitParameters); diff --git a/Reconstruction/PFA/Pandora/MatrixPandora/include/CaloHitCreator.h b/Reconstruction/PFA/Pandora/MatrixPandora/include/CaloHitCreator.h index c7714205c558c0df913d3f28b2198204aec59f58..7dfb684f8d155291f5d1e1e3a02154da7f2ad1dc 100644 --- a/Reconstruction/PFA/Pandora/MatrixPandora/include/CaloHitCreator.h +++ b/Reconstruction/PFA/Pandora/MatrixPandora/include/CaloHitCreator.h @@ -1,5 +1,4 @@ /** - * @file MarlinPandora/include/CaloHitCreator.h * * @brief Header file for the calo hit creator class. * @@ -116,9 +115,7 @@ public: /** * @brief Create calo hits * - * @param pLCEvent the lcio event */ - //pandora::StatusCode CreateCaloHits(const LCEvent *const pLCEvent); pandora::StatusCode CreateCaloHits(const CollectionMaps& collectionMaps); /** @@ -137,44 +134,36 @@ private: /** * @brief Create ecal calo hits * - * @param pLCEvent the lcio event */ - //pandora::StatusCode CreateECalCaloHits(const EVENT::LCEvent *const pLCEvent); pandora::StatusCode CreateECalCaloHits(const CollectionMaps& collectionMaps); /** * @brief Create hcal calo hits * - * @param pLCEvent the lcio event */ pandora::StatusCode CreateHCalCaloHits(const CollectionMaps& collectionMaps); /** * @brief Create muon calo hits * - * @param pLCEvent the lcio event */ pandora::StatusCode CreateMuonCaloHits(const CollectionMaps& collectionMaps); /** * @brief Create lcal calo hits * - * @param pLCEvent the lcio event */ pandora::StatusCode CreateLCalCaloHits(const CollectionMaps& collectionMaps); /** * @brief Create lhcal calo hits * - * @param pLCEvent the lcio event */ pandora::StatusCode CreateLHCalCaloHits(const CollectionMaps& collectionMaps); /** * @brief Get common calo hit properties: position, parent address, input energy and time * - * @param pCaloHit the lcio calorimeter hit - * @param caloHitParameters the calo hit parameters to populate */ void GetCommonCaloHitProperties(const edm4hep::CalorimeterHit *const pCaloHit, PandoraApi::CaloHit::Parameters &caloHitParameters) const; diff --git a/Reconstruction/PFA/Pandora/MatrixPandora/include/GeometryCreator.h b/Reconstruction/PFA/Pandora/MatrixPandora/include/GeometryCreator.h index 0a0c1ae2d6953e7d1b066e4107fda52c55396ff1..6bd9c4ba2f9c3cfc48591533377d296d7df62d17 100644 --- a/Reconstruction/PFA/Pandora/MatrixPandora/include/GeometryCreator.h +++ b/Reconstruction/PFA/Pandora/MatrixPandora/include/GeometryCreator.h @@ -1,5 +1,4 @@ /** - * @file MarlinPandora/include/GeometryCreator.h * * @brief Header file for the geometry creator class. * diff --git a/Reconstruction/PFA/Pandora/MatrixPandora/include/MCParticleCreator.h b/Reconstruction/PFA/Pandora/MatrixPandora/include/MCParticleCreator.h index 4ad9f27553b03b7c795d6fc457fa6ecddaa46463..3c1ea57efa4ebd5d6015830af7feec301deac09f 100644 --- a/Reconstruction/PFA/Pandora/MatrixPandora/include/MCParticleCreator.h +++ b/Reconstruction/PFA/Pandora/MatrixPandora/include/MCParticleCreator.h @@ -1,5 +1,4 @@ /** - * @file MarlinPandora/include/MCParticleCreator.h * * @brief Header file for the mc particle creator class. * diff --git a/Reconstruction/PFA/Pandora/MatrixPandora/include/PandoraMatrixAlg.h b/Reconstruction/PFA/Pandora/MatrixPandora/include/PandoraMatrixAlg.h index 63210760fdcfe3b7b4e08697ce340639dee6edf3..41c0e809b64bc248134aac25c4479300fa82482c 100644 --- a/Reconstruction/PFA/Pandora/MatrixPandora/include/PandoraMatrixAlg.h +++ b/Reconstruction/PFA/Pandora/MatrixPandora/include/PandoraMatrixAlg.h @@ -49,7 +49,6 @@ */ namespace pandora {class Pandora;} -class IEventSeeder; class CollectionMaps { @@ -59,10 +58,6 @@ public: */ CollectionMaps(); void clear(); - std::map<std::string, const edm4hep::MCParticleCollection*> CollectionMap_MC; - std::map<std::string, const edm4hep::CalorimeterHitCollection*> CollectionMap_CaloHit; - std::map<std::string, const edm4hep::VertexCollection*> CollectionMap_Vertex; - std::map<std::string, const edm4hep::TrackCollection*> CollectionMap_Track; std::map<std::string, std::vector<edm4hep::MCParticle> > collectionMap_MC; std::map<std::string, std::vector<edm4hep::CalorimeterHit> > collectionMap_CaloHit; @@ -76,7 +71,6 @@ public: class PandoraMatrixAlg : public GaudiAlgorithm { - //friend class AlgFactory<PandoraMatrixAlg>;//gives error in 97 version public: @@ -95,7 +89,6 @@ public: */ virtual StatusCode finalize() ; - //void FinaliseSteeringParameters(); void FinaliseSteeringParameters(ISvcLocator* svcloc); pandora::StatusCode RegisterUserComponents() const; void Reset(); @@ -129,20 +122,17 @@ public: */ const pandora::Pandora *GetPandora() const; StatusCode updateMap(); - StatusCode updateMap(CollectionMaps & tmp_map); StatusCode Ana(); StatusCode CreateMCRecoParticleAssociation(); - //StatusCode Create_MC(); protected: typedef std::vector<float> FloatVec; int _nEvt ; - IEventSeeder * _SEEDER; - Gaudi::Property< std::string > m_PandoraSettingsXmlFile { this, "PandoraSettingsDefault_xml", "/junofs/users/wxfang/MyGit/MarlinPandora/scripts/PandoraSettingsDefault_wx.xml" }; + Gaudi::Property< std::string > m_PandoraSettingsXmlFile { this, "PandoraSettingsDefault_xml", "PandoraSettingsDefault_wx.xml" }; Gaudi::Property<int> m_NEventsToSkip { this, "NEventsToSkip", 0 }; Gaudi::Property< std::vector<std::string> > m_TrackCollections{ this, "TrackCollections", {"MarlinTrkTracks"} }; @@ -296,26 +286,12 @@ protected: Gaudi::Property< std::string > m_AnaOutput{ this, "AnaOutput", "/junofs/users/wxfang/MyGit/CEPCSW/Reconstruction/PFA/Pandora/GaudiPandora/Ana.root" }; //###################### + std::map< std::string, std::string > m_collections; + Gaudi::Property<std::vector<std::string>> m_readCols{this, "collections", {}, "Places of collections to read"}; + //the map of collection name to its corresponding DataHandle + std::map<std::string, DataObjectHandleBase*> m_dataHandles; DataHandle<edm4hep::MCParticleCollection> m_mcParCol_r {"MCParticle", Gaudi::DataHandle::Reader, this}; - DataHandle<edm4hep::CalorimeterHitCollection> m_ECALBarrel_r{"ECALBarrel", Gaudi::DataHandle::Reader, this}; - DataHandle<edm4hep::CalorimeterHitCollection> m_ECALEndcap_r{"ECALEndcap", Gaudi::DataHandle::Reader, this}; - DataHandle<edm4hep::CalorimeterHitCollection> m_ECALOther_r {"ECALOther", Gaudi::DataHandle::Reader, this}; - DataHandle<edm4hep::CalorimeterHitCollection> m_HCALBarrel_r{"HCALBarrel", Gaudi::DataHandle::Reader, this}; - DataHandle<edm4hep::CalorimeterHitCollection> m_HCALEndcap_r{"HCALEndcap", Gaudi::DataHandle::Reader, this}; - DataHandle<edm4hep::CalorimeterHitCollection> m_HCALOther_r {"HCALOther", Gaudi::DataHandle::Reader, this}; - DataHandle<edm4hep::CalorimeterHitCollection> m_MUON_r {"MUON", Gaudi::DataHandle::Reader, this}; - DataHandle<edm4hep::CalorimeterHitCollection> m_LCAL_r {"LCAL", Gaudi::DataHandle::Reader, this}; - DataHandle<edm4hep::CalorimeterHitCollection> m_LHCAL_r {"LHCAL", Gaudi::DataHandle::Reader, this}; - DataHandle<edm4hep::CalorimeterHitCollection> m_BCAL_r {"BCAL", Gaudi::DataHandle::Reader, this}; - - DataHandle<edm4hep::VertexCollection> m_KinkVertices_r {"KinkVertices",Gaudi::DataHandle::Reader, this}; - DataHandle<edm4hep::VertexCollection> m_ProngVertices_r {"ProngVertices",Gaudi::DataHandle::Reader, this}; - DataHandle<edm4hep::VertexCollection> m_SplitVertices_r {"SplitVertices",Gaudi::DataHandle::Reader, this}; - DataHandle<edm4hep::VertexCollection> m_V0Vertices_r {"V0Vertices",Gaudi::DataHandle::Reader, this}; - DataHandle<edm4hep::TrackCollection> m_MarlinTrkTracks_r {"MarlinTrkTracks",Gaudi::DataHandle::Reader, this}; - DataHandle<edm4hep::MCRecoCaloAssociationCollection> m_MCRecoCaloAssociation_r {"MCRecoCaloAssociation",Gaudi::DataHandle::Reader, this}; - DataHandle<edm4hep::MCRecoTrackerAssociationCollection> m_MCRecoTrackerAssociation_r {"MCRecoTrackerAssociation",Gaudi::DataHandle::Reader, this}; DataHandle<edm4hep::ClusterCollection> m_ClusterCollection_w {"PandoraClusters",Gaudi::DataHandle::Writer, this}; DataHandle<edm4hep::ReconstructedParticleCollection> m_ReconstructedParticleCollection_w {"PandoraPFOs" ,Gaudi::DataHandle::Writer, this}; diff --git a/Reconstruction/PFA/Pandora/MatrixPandora/include/PfoCreator.h b/Reconstruction/PFA/Pandora/MatrixPandora/include/PfoCreator.h index b9679a48933ff9aaf90550d1a9e6c170c1f0c91f..2cc8d6327e0a99f3782164346c992770d9fe0dd4 100644 --- a/Reconstruction/PFA/Pandora/MatrixPandora/include/PfoCreator.h +++ b/Reconstruction/PFA/Pandora/MatrixPandora/include/PfoCreator.h @@ -1,5 +1,4 @@ /** - * @file MarlinPandora/include/PfoCreator.h * * @brief Header file for the pfo creator class. * diff --git a/Reconstruction/PFA/Pandora/MatrixPandora/include/TrackCreator.h b/Reconstruction/PFA/Pandora/MatrixPandora/include/TrackCreator.h index 1e1588c2629db7a4f7728b9edb62b86f5cf3fd0e..b5daa829d8b6ef94a5d14d7a01c9a93558fde5e3 100644 --- a/Reconstruction/PFA/Pandora/MatrixPandora/include/TrackCreator.h +++ b/Reconstruction/PFA/Pandora/MatrixPandora/include/TrackCreator.h @@ -1,5 +1,4 @@ /** - * @file MarlinPandora/include/TrackCreator.h * * @brief Header file for the track creator class. * @@ -27,9 +26,7 @@ namespace gear { class GearMgr; } class CollectionMaps; typedef std::vector<const edm4hep::Track *> TrackVector; -//typedef std::set<const edm4hep::Track *> TrackList; typedef std::set<unsigned int> TrackList; -//typedef std::map<edm4hep::Track *, int> TrackToPidMap; typedef std::map<edm4hep::ConstTrack, int> TrackToPidMap; /* inline LCCollectionVec *newTrkCol(const std::string &name, LCEvent *evt , bool isSubset) @@ -127,17 +124,14 @@ public: /** * @brief Create associations between tracks, V0s, kinks, etc * - * @param pLCEvent the lcio event */ - //pandora::StatusCode CreateTrackAssociations(const EVENT::LCEvent *const pLCEvent); pandora::StatusCode CreateTrackAssociations(const CollectionMaps& collectionMaps); + const edm4hep::Track* GetTrackAddress(const CollectionMaps& collectionMaps, const edm4hep::ConstTrack& pTrack ); /** * @brief Create tracks, insert user code here * - * @param pLCEvent the lcio event */ - //pandora::StatusCode CreateTracks(EVENT::LCEvent *pLCEvent); pandora::StatusCode CreateTracks(const CollectionMaps& collectionMaps); /** @@ -156,140 +150,102 @@ private: /** * @brief Extract kink information from specified lcio collections * - * @param pLCEvent the lcio event */ - //pandora::StatusCode ExtractKinks(const EVENT::LCEvent *const pLCEvent); pandora::StatusCode ExtractKinks(const CollectionMaps& collectionMaps); /** * @brief Extract prong and split information from specified lcio collections * - * @param pLCEvent the lcio event */ - //pandora::StatusCode ExtractProngsAndSplits(const EVENT::LCEvent *const pLCEvent); + pandora::StatusCode ExtractProngsAndSplits(const CollectionMaps& collectionMaps); /** * @brief Extract v0 information from specified lcio collections * - * @param pLCEvent the lcio event */ - //pandora::StatusCode ExtractV0s(const EVENT::LCEvent *const pLCEvent); + pandora::StatusCode ExtractV0s(const CollectionMaps& collectionMaps); /** * @brief Whether the track vertex conflicts with previously provided relationship information * - * @param trackVec the vector of tracks associated with the vertex */ - //bool IsConflictingRelationship(const EVENT::TrackVec &trackVec) const; - //bool IsConflictingRelationship(edm4hep::ConstTrack &trackVec) const; bool IsConflictingRelationship(const edm4hep::ConstReconstructedParticle &Particle) const; /** * @brief Whether a track is a v0 track * - * @param pTrack the lcio track * * @return boolean */ - //bool IsV0(const EVENT::Track *const pTrack) const; bool IsV0(unsigned int pTrack_id) const; /** * @brief Whether a track is a parent track * - * @param pTrack the lcio track * * @return boolean */ - //bool IsParent(const EVENT::Track *const pTrack) const; bool IsParent(unsigned int pTrack_id) const; /** * @brief Whether a track is a daughter track * - * @param pTrack the lcio track * * @return boolean */ - //bool IsDaughter(const EVENT::Track *const pTrack) const; bool IsDaughter(unsigned int pTrack_id) const; /** * @brief Copy track states stored in lcio tracks to pandora track parameters * - * @param pTrack the lcio track - * @param trackParameters the track parameters */ - //void GetTrackStates(const EVENT::Track *const pTrack, PandoraApi::Track::Parameters &trackParameters) const; void GetTrackStates(const edm4hep::Track *const pTrack, PandoraApi::Track::Parameters &trackParameters) const; /** * @brief Copy track state from lcio track state instance to pandora input track state * - * @param pTrackState the lcio track state instance - * @param inputTrackState the pandora input track state */ - //void CopyTrackState(const TrackState *const pTrackState, pandora::InputTrackState &inputTrackState) const; void CopyTrackState(const edm4hep::TrackState & pTrackState, pandora::InputTrackState &inputTrackState) const; /** * @brief Obtain track time when it reaches ECAL * - * @param pTrack the lcio track */ - //float CalculateTrackTimeAtCalorimeter(const EVENT::Track *const pTrack) const; float CalculateTrackTimeAtCalorimeter(const edm4hep::Track *const pTrack) const; /** * @brief Decide whether track reaches the ecal surface * - * @param pTrack the lcio track - * @param trackParameters the track parameters */ - //void TrackReachesECAL(const EVENT::Track *const pTrack, PandoraApi::Track::Parameters &trackParameters) const; void TrackReachesECAL(const edm4hep::Track *const pTrack, PandoraApi::Track::Parameters &trackParameters) const; - //void TrackReachesECAL(const edm4hep::Track& pTrack, PandoraApi::Track::Parameters &trackParameters) const; /** * @brief Determine whether a track can be used to form a pfo under the following conditions: * 1) if the track proves to be associated with a cluster, OR * 2) if the track proves to have no cluster associations * - * @param pTrack the lcio track - * @param trackParameters the track parameters */ - //void DefineTrackPfoUsage(const EVENT::Track *const pTrack, PandoraApi::Track::Parameters &trackParameters) const; void DefineTrackPfoUsage(const edm4hep::Track *const pTrack, PandoraApi::Track::Parameters &trackParameters) const; /** * @brief Whether track passes the quality cuts required in order to be used to form a pfo * - * @param pTrack the lcio track - * @param trackParameters the track parameters * * @return boolean */ - //bool PassesQualityCuts(const EVENT::Track *const pTrack, const PandoraApi::Track::Parameters &trackParameters) const; bool PassesQualityCuts(const edm4hep::Track *const pTrack, const PandoraApi::Track::Parameters &trackParameters) const; /** * @brief Get number of hits in TPC of a track * - * @param pTrack the lcio track - * - * @return number of hits in TPC of a track */ - //int GetNTpcHits(const EVENT::Track *const pTrack) const; int GetNTpcHits(const edm4hep::Track *const pTrack) const; /** * @brief Get number of hits in FTD of a track * - * @param pTrack the lcio track * - * @return number of hits in FTDof a track */ - //int GetNFtdHits(const EVENT::Track *const pTrack) const; int GetNFtdHits(const edm4hep::Track *const pTrack) const; const Settings m_settings; ///< The track creator settings @@ -345,28 +301,22 @@ inline void TrackCreator::Reset() //------------------------------------------------------------------------------------------------------------------------------------------ -//inline bool TrackCreator::IsV0(const Track *const pTrack) const inline bool TrackCreator::IsV0(unsigned int pTrack_id) const // should check here, if id is correct one to do this { - //return (m_v0TrackList.end() != m_v0TrackList.find(pTrack)); return (m_v0TrackList.end() != m_v0TrackList.find(pTrack_id)); } //------------------------------------------------------------------------------------------------------------------------------------------ -//inline bool TrackCreator::IsParent(const Track *const pTrack) const inline bool TrackCreator::IsParent(unsigned int pTrack_id) const { - //return (m_parentTrackList.end() != m_parentTrackList.find(pTrack)); return (m_parentTrackList.end() != m_parentTrackList.find(pTrack_id)); } //------------------------------------------------------------------------------------------------------------------------------------------ -//inline bool TrackCreator::IsDaughter(const Track *const pTrack) const inline bool TrackCreator::IsDaughter(unsigned int pTrack_id) const { - //return (m_daughterTrackList.end() != m_daughterTrackList.find(pTrack)); return (m_daughterTrackList.end() != m_daughterTrackList.find(pTrack_id)); } diff --git a/Reconstruction/PFA/Pandora/MatrixPandora/src/CaloHitCreator.cpp b/Reconstruction/PFA/Pandora/MatrixPandora/src/CaloHitCreator.cpp index 3196cde227b8c14c647aa2e042c89b36ba91ad9e..a6e812806716d76aaadbcda877c1e38959ca6ca0 100644 --- a/Reconstruction/PFA/Pandora/MatrixPandora/src/CaloHitCreator.cpp +++ b/Reconstruction/PFA/Pandora/MatrixPandora/src/CaloHitCreator.cpp @@ -1,5 +1,4 @@ /** - * @file MarlinPandora/src/CaloHitCreator.cc * * @brief Implementation of the calo hit creator class. * @@ -69,12 +68,6 @@ CaloHitCreator::CaloHitCreator(const Settings &settings, const pandora::Pandora dd4hep::rec::LayeredCalorimeterData* Data = detElement.extension<dd4hep::rec::LayeredCalorimeterData>() ; if(!Data) throw "Failed to get LayeredCalorimeterData ..."; m_cellIDConverter = new dd4hep::rec::CellIDPositionConverter(*m_dd4hep); - /* - m_compact = "/junofs/users/wxfang/MyGit/lcg97/UsingK4FWCore/0514/CEPCSW/Detector/DetEcalMatrix/compact/det.xml"; - dd4hep::Detector& description = dd4hep::Detector::getInstance(); - description.fromCompact( m_compact ); - m_cellIDConverter = new dd4hep::rec::CellIDPositionConverter(description); - */ m_eCalBarrelOuterZ = Data->extent[3]; m_eCalBarrelInnerPhi0 = Data->inner_phi0; @@ -202,14 +195,7 @@ pandora::StatusCode CaloHitCreator::CreateECalCaloHits(const CollectionMaps& col PandoraApi::CaloHit::Parameters caloHitParameters; caloHitParameters.m_hitType = pandora::ECAL; - caloHitParameters.m_isDigital = false;//FIXME, maybe this means the MIP thershold cut haven't be applied yet ? - /* - dd4hep::Position position = m_cellIDConverter->position(pCaloHit->getCellID()); - long id_sys,id_x,id_y,id_z ; - GetCoding(pCaloHit, id_sys, id_x, id_y, id_z); - std::cout << "ECAL id =" << pCaloHit->getCellID()<<",sys="<<id_sys<<",x="<<id_x<<",y="<<id_y<<",z="<<id_z << std::endl; - std::cout << "ECAL id =" << pCaloHit->getCellID()<<",x="<<position.x()<<",y="<<position.y()<<",z="<<position.z() << std::endl; - */ + caloHitParameters.m_isDigital = false; caloHitParameters.m_isInOuterSamplingLayer = false; this->GetCommonCaloHitProperties(pCaloHit, caloHitParameters); @@ -224,7 +210,7 @@ pandora::StatusCode CaloHitCreator::CreateECalCaloHits(const CollectionMaps& col } else { // will not be used for ECAL Matrix - caloHitParameters.m_layer = cellIdDecoder(pCaloHit)[layerCoding.c_str()] + 1;//FIXME, should use + 1? because the decoded layer is start from 0. + caloHitParameters.m_layer = cellIdDecoder(pCaloHit)[layerCoding.c_str()] + 1; this->GetEndCapCaloHitProperties(pCaloHit, endcapLayerLayout, caloHitParameters, absorberCorrection); caloHitParameters.m_hadronicEnergy = eCalToHadGeVEndCap * pCaloHit->getEnergy(); } diff --git a/Reconstruction/PFA/Pandora/MatrixPandora/src/GeometryCreator.cpp b/Reconstruction/PFA/Pandora/MatrixPandora/src/GeometryCreator.cpp index 3ab6cc805b00cc2f2a74c5baecdeba10c9cb1907..343b6045398c032c9b71c3c64e18081437a5f2b4 100644 --- a/Reconstruction/PFA/Pandora/MatrixPandora/src/GeometryCreator.cpp +++ b/Reconstruction/PFA/Pandora/MatrixPandora/src/GeometryCreator.cpp @@ -1,5 +1,4 @@ /** - * @file MarlinPandora/src/GeometryCreator.cc * * @brief Implementation of the geometry creator class. * diff --git a/Reconstruction/PFA/Pandora/MatrixPandora/src/MCParticleCreator.cpp b/Reconstruction/PFA/Pandora/MatrixPandora/src/MCParticleCreator.cpp index 21a2a741bcc28617dba766d43600ee1f818ded8e..c05190ebebc7efc39813bcaf22820743406afa1d 100644 --- a/Reconstruction/PFA/Pandora/MatrixPandora/src/MCParticleCreator.cpp +++ b/Reconstruction/PFA/Pandora/MatrixPandora/src/MCParticleCreator.cpp @@ -1,5 +1,4 @@ /** - * @file MarlinPandora/src/MCParticleCreator.cc * * @brief Implementation of the mc particle creator class. * diff --git a/Reconstruction/PFA/Pandora/MatrixPandora/src/PandoraMatrixAlg.cpp b/Reconstruction/PFA/Pandora/MatrixPandora/src/PandoraMatrixAlg.cpp index 566f61e5b116ba626db23dba934b495e3f2e8b8d..d59c1a063d9ae4734ef1944d42db35930e594edf 100644 --- a/Reconstruction/PFA/Pandora/MatrixPandora/src/PandoraMatrixAlg.cpp +++ b/Reconstruction/PFA/Pandora/MatrixPandora/src/PandoraMatrixAlg.cpp @@ -41,24 +41,6 @@ PandoraMatrixAlg::PandoraMatrixAlg(const std::string& name, ISvcLocator* svcLoc) { m_CollectionMaps = new CollectionMaps(); - declareProperty("ReadMCParticle" , m_mcParCol_r, "Handle of the MCParticle input collection" ); - declareProperty("ReadECALBarrel" , m_ECALBarrel_r, "Handle of the ECALBarrel input collection" ); - declareProperty("ReadECALEndcap" , m_ECALEndcap_r, "Handle of the ECALEndcap input collection" ); - declareProperty("ReadECALOther" , m_ECALOther_r, "Handle of the ECALOther input collection" ); - declareProperty("ReadHCALBarrel" , m_HCALBarrel_r, "Handle of the HCALBarrel input collection" ); - declareProperty("ReadHCALEndcap" , m_HCALEndcap_r, "Handle of the HCALEndcap input collection" ); - declareProperty("ReadHCALOther" , m_HCALOther_r, "Handle of the HCALOther input collection" ); - declareProperty("ReadMUON" , m_MUON_r, "Handle of the MUON input collection" ); - declareProperty("ReadLCAL" , m_LCAL_r, "Handle of the LCAL input collection" ); - declareProperty("ReadLHCAL" , m_LHCAL_r, "Handle of the LHCAL input collection" ); - declareProperty("ReadBCAL" , m_BCAL_r, "Handle of the BCAL input collection" ); - declareProperty("ReadKinkVertices" , m_KinkVertices_r, "Handle of the KinkVertices input collection" ); - declareProperty("ReadProngVertices" , m_ProngVertices_r, "Handle of the ProngVertices input collection" ); - declareProperty("ReadSplitVertices" , m_SplitVertices_r, "Handle of the SplitVertices input collection" ); - declareProperty("ReadV0Vertices" , m_V0Vertices_r, "Handle of the V0Vertices input collection" ); - declareProperty("ReadTracks" , m_MarlinTrkTracks_r, "Handle of the Tracks input collection" ); - declareProperty("MCRecoCaloAssociation" , m_MCRecoCaloAssociation_r, "Handle of the MCRecoCaloAssociation input collection" ); - declareProperty("MCRecoTrackerAssociation" , m_MCRecoTrackerAssociation_r, "Handle of the MCRecoTrackerAssociation input collection" ); declareProperty("WriteClusterCollection" , m_ClusterCollection_w, "Handle of the ClusterCollection output collection" ); declareProperty("WriteReconstructedParticleCollection", m_ReconstructedParticleCollection_w, "Handle of the ReconstructedParticleCollection output collection" ); declareProperty("WriteVertexCollection" , m_VertexCollection_w, "Handle of the VertexCollection output collection" ); @@ -138,6 +120,41 @@ StatusCode PandoraMatrixAlg::initialize() m_tree->Branch("m_hits_z", &m_hits_z); m_tree->Branch("m_hits_E", &m_hits_E); + for ( const auto& col : m_readCols ) { + auto seperater = col.find(':'); + std::string colType = col.substr(0, seperater); + std::string colName = col.substr(seperater+1); + m_collections[colName] = colType; + + if ( colType == "MCParticle" ) { + m_dataHandles[colName] = + new DataHandle<edm4hep::MCParticleCollection>(colName, Gaudi::DataHandle::Reader, this); + } + else if ( colType == "Track" ) { + m_dataHandles[colName] = + new DataHandle<edm4hep::TrackCollection>(colName, Gaudi::DataHandle::Reader, this); + } + else if ( colType == "CalorimeterHit" ) { + m_dataHandles[colName] = + new DataHandle<edm4hep::CalorimeterHitCollection>(colName, Gaudi::DataHandle::Reader, this); + } + else if ( colType == "Vertex" ) { + m_dataHandles[colName] = + new DataHandle<edm4hep::VertexCollection>(colName, Gaudi::DataHandle::Reader, this); + } + else if ( colType == "MCRecoTrackerAssociation" ) { + m_dataHandles[colName] = + new DataHandle<edm4hep::MCRecoTrackerAssociationCollection>(colName, Gaudi::DataHandle::Reader, this); + } + else if ( colType == "MCRecoCaloAssociation" ) { + m_dataHandles[colName] = + new DataHandle<edm4hep::MCRecoCaloAssociationCollection>(colName, Gaudi::DataHandle::Reader, this); + } + else { + error() << "invalid collection type: " << colType << endmsg; + return StatusCode::FAILURE; + } + } // XML file m_settings.m_pandoraSettingsXmlFile = m_PandoraSettingsXmlFile ; @@ -320,9 +337,9 @@ StatusCode PandoraMatrixAlg::execute() PANDORA_THROW_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, m_pMCParticleCreator->CreateMCParticles(*m_CollectionMaps)); PANDORA_THROW_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, m_pCaloHitCreator->CreateCaloHits(*m_CollectionMaps)); PANDORA_THROW_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, m_pMCParticleCreator->CreateCaloHitToMCParticleRelationships(*m_CollectionMaps, m_pCaloHitCreator->GetCalorimeterHitVector() )); - //PANDORA_THROW_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, m_pTrackCreator->CreateTrackAssociations(*m_CollectionMaps)); - //PANDORA_THROW_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, m_pTrackCreator->CreateTracks(*m_CollectionMaps)); - //PANDORA_THROW_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, m_pMCParticleCreator->CreateTrackToMCParticleRelationships(*m_CollectionMaps, m_pTrackCreator->GetTrackVector() )); + PANDORA_THROW_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, m_pTrackCreator->CreateTrackAssociations(*m_CollectionMaps)); + PANDORA_THROW_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, m_pTrackCreator->CreateTracks(*m_CollectionMaps)); + PANDORA_THROW_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, m_pMCParticleCreator->CreateTrackToMCParticleRelationships(*m_CollectionMaps, m_pTrackCreator->GetTrackVector() )); PANDORA_THROW_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, PandoraApi::ProcessEvent(*m_pPandora)); PANDORA_THROW_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, m_pPfoCreator->CreateParticleFlowObjects(*m_CollectionMaps, m_ClusterCollection_w, m_ReconstructedParticleCollection_w, m_VertexCollection_w)); @@ -439,10 +456,6 @@ CollectionMaps::CollectionMaps() } void CollectionMaps::clear() { -CollectionMap_MC.clear(); -CollectionMap_CaloHit.clear(); -CollectionMap_Vertex.clear(); -CollectionMap_Track.clear(); collectionMap_MC.clear(); collectionMap_CaloHit.clear(); collectionMap_Vertex.clear(); @@ -451,373 +464,99 @@ collectionMap_CaloRel.clear(); collectionMap_TrkRel.clear(); } + + StatusCode PandoraMatrixAlg::updateMap() { - const edm4hep::MCParticleCollection* MCParticle = nullptr; - const edm4hep::CalorimeterHitCollection* ECALBarrel = nullptr; - const edm4hep::CalorimeterHitCollection* ECALEndcap = nullptr; - const edm4hep::CalorimeterHitCollection* ECALOther = nullptr; - const edm4hep::CalorimeterHitCollection* HCALBarrel = nullptr; - const edm4hep::CalorimeterHitCollection* HCALEndcap = nullptr; - const edm4hep::CalorimeterHitCollection* HCALOther = nullptr; - const edm4hep::CalorimeterHitCollection* MUON = nullptr; - const edm4hep::CalorimeterHitCollection* LCAL = nullptr; - const edm4hep::CalorimeterHitCollection* LHCAL = nullptr; - const edm4hep::CalorimeterHitCollection* BCAL = nullptr; - const edm4hep::VertexCollection* KinkVertices = nullptr; - const edm4hep::VertexCollection* ProngVertices = nullptr; - const edm4hep::VertexCollection* SplitVertices = nullptr; - const edm4hep::VertexCollection* V0Vertices = nullptr; - const edm4hep::TrackCollection* MarlinTrkTracks = nullptr; - const edm4hep::MCRecoCaloAssociationCollection* mcRecoCaloAssociation = nullptr; - const edm4hep::MCRecoTrackerAssociationCollection* mcRecoTrackerAssociation = nullptr; - StatusCode sc = StatusCode::SUCCESS; - sc = getCol(m_mcParCol_r , MCParticle ); - sc = getCol(m_ECALBarrel_r, ECALBarrel ); - sc = getCol(m_ECALEndcap_r, ECALEndcap ); - sc = getCol(m_ECALOther_r , ECALOther ); - sc = getCol(m_HCALBarrel_r, HCALBarrel ); - sc = getCol(m_HCALEndcap_r, HCALEndcap ); - sc = getCol(m_HCALOther_r , HCALOther ); - sc = getCol(m_MUON_r , MUON ); - sc = getCol(m_LCAL_r , LCAL ); - sc = getCol(m_LHCAL_r , LHCAL ); - sc = getCol(m_BCAL_r , BCAL ); - sc = getCol(m_KinkVertices_r , KinkVertices ); - sc = getCol(m_ProngVertices_r , ProngVertices); - sc = getCol(m_SplitVertices_r , SplitVertices); - sc = getCol(m_V0Vertices_r , V0Vertices ); - sc = getCol(m_MarlinTrkTracks_r , MarlinTrkTracks ); - sc = getCol(m_MCRecoCaloAssociation_r , mcRecoCaloAssociation ); - sc = getCol(m_MCRecoTrackerAssociation_r , mcRecoTrackerAssociation); - - if (NULL != MCParticle ) - { - std::vector<edm4hep::MCParticle> v_mc; - m_CollectionMaps->CollectionMap_MC ["MCParticle"] = MCParticle ; - m_CollectionMaps->collectionMap_MC ["MCParticle"] = v_mc; - for(unsigned int i=0 ; i< MCParticle->size(); i++) m_CollectionMaps->collectionMap_MC ["MCParticle"].push_back(MCParticle->at(i)); - } - if (NULL != ECALBarrel ) - { - std::vector<edm4hep::CalorimeterHit> v_cal; - m_CollectionMaps->CollectionMap_CaloHit["ECALBarrel"] = ECALBarrel ; - m_CollectionMaps->collectionMap_CaloHit["ECALBarrel"] = v_cal ; - for(unsigned int i=0 ; i< ECALBarrel->size(); i++) m_CollectionMaps->collectionMap_CaloHit ["ECALBarrel"].push_back(ECALBarrel->at(i)); - } - if (NULL != ECALEndcap ) - { - std::vector<edm4hep::CalorimeterHit> v_cal; - m_CollectionMaps->CollectionMap_CaloHit["ECALEndcap"] = ECALEndcap ; - m_CollectionMaps->collectionMap_CaloHit["ECALEndcap"] = v_cal ; - for(unsigned int i=0 ; i< ECALEndcap->size(); i++) m_CollectionMaps->collectionMap_CaloHit ["ECALEndcap"].push_back(ECALEndcap->at(i)); - } - if (NULL != ECALOther ) - { - std::vector<edm4hep::CalorimeterHit> v_cal; - m_CollectionMaps->CollectionMap_CaloHit["ECALOther"] = ECALOther ; - m_CollectionMaps->collectionMap_CaloHit["ECALOther"] = v_cal ; - for(unsigned int i=0 ; i< ECALOther->size(); i++) m_CollectionMaps->collectionMap_CaloHit ["ECALOther"].push_back(ECALOther->at(i)); - } - if (NULL != HCALBarrel ) - { - std::vector<edm4hep::CalorimeterHit> v_cal; - m_CollectionMaps->CollectionMap_CaloHit["HCALBarrel"] = HCALBarrel ; - m_CollectionMaps->collectionMap_CaloHit["HCALBarrel"] = v_cal ; - for(unsigned int i=0 ; i< HCALBarrel->size(); i++) m_CollectionMaps->collectionMap_CaloHit ["HCALBarrel"].push_back(HCALBarrel->at(i)); - } - if (NULL != HCALEndcap ) - { - std::vector<edm4hep::CalorimeterHit> v_cal; - m_CollectionMaps->CollectionMap_CaloHit["HCALEndcap"] = HCALEndcap ; - m_CollectionMaps->collectionMap_CaloHit["HCALEndcap"] = v_cal ; - for(unsigned int i=0 ; i< HCALEndcap->size(); i++) m_CollectionMaps->collectionMap_CaloHit ["HCALEndcap"].push_back(HCALEndcap->at(i)); - } - if (NULL != HCALOther ) - { - std::vector<edm4hep::CalorimeterHit> v_cal; - m_CollectionMaps->CollectionMap_CaloHit["HCALOther"] = HCALOther ; - m_CollectionMaps->collectionMap_CaloHit["HCALOther"] = v_cal ; - for(unsigned int i=0 ; i< HCALOther->size(); i++) m_CollectionMaps->collectionMap_CaloHit ["HCALOther"].push_back(HCALOther->at(i)); - } - if (NULL != MUON ) - { - std::vector<edm4hep::CalorimeterHit> v_cal; - m_CollectionMaps->CollectionMap_CaloHit["MUON"] = MUON ; - m_CollectionMaps->collectionMap_CaloHit["MUON"] = v_cal ; - for(unsigned int i=0 ; i< MUON->size(); i++) m_CollectionMaps->collectionMap_CaloHit ["MUON"].push_back(MUON->at(i)); - } - if (NULL != LCAL ) - { - std::vector<edm4hep::CalorimeterHit> v_cal; - m_CollectionMaps->CollectionMap_CaloHit["LCAL"] = LCAL ; - m_CollectionMaps->collectionMap_CaloHit["LCAL"] = v_cal ; - for(unsigned int i=0 ; i< LCAL->size(); i++) m_CollectionMaps->collectionMap_CaloHit ["LCAL"].push_back(LCAL->at(i)); - } - if (NULL != LHCAL ) - { - std::vector<edm4hep::CalorimeterHit> v_cal; - m_CollectionMaps->CollectionMap_CaloHit["LHCAL"] = LHCAL ; - m_CollectionMaps->collectionMap_CaloHit["LHCAL"] = v_cal ; - for(unsigned int i=0 ; i< LHCAL->size(); i++) m_CollectionMaps->collectionMap_CaloHit ["LHCAL"].push_back(LHCAL->at(i)); - } - if (NULL != BCAL ) - { - std::vector<edm4hep::CalorimeterHit> v_cal; - m_CollectionMaps->CollectionMap_CaloHit["BCAL"] = BCAL ; - m_CollectionMaps->collectionMap_CaloHit["BCAL"] = v_cal ; - for(unsigned int i=0 ; i< BCAL->size(); i++) m_CollectionMaps->collectionMap_CaloHit ["BCAL"].push_back(BCAL->at(i)); - } - if (NULL != KinkVertices ) - { - std::vector<edm4hep::Vertex> v_cal; - m_CollectionMaps->CollectionMap_Vertex["KinkVertices"] = KinkVertices ; - m_CollectionMaps->collectionMap_Vertex["KinkVertices"] = v_cal ; - for(unsigned int i=0 ; i< KinkVertices->size(); i++) m_CollectionMaps->collectionMap_Vertex ["KinkVertices"].push_back(KinkVertices->at(i)); - } - if (NULL != ProngVertices ) - { - std::vector<edm4hep::Vertex> v_cal; - m_CollectionMaps->CollectionMap_Vertex["ProngVertices"] = ProngVertices ; - m_CollectionMaps->collectionMap_Vertex["ProngVertices"] = v_cal ; - for(unsigned int i=0 ; i< ProngVertices->size(); i++) m_CollectionMaps->collectionMap_Vertex ["ProngVertices"].push_back(ProngVertices->at(i)); - } - if (NULL != SplitVertices ) - { - std::vector<edm4hep::Vertex> v_cal; - m_CollectionMaps->CollectionMap_Vertex["SplitVertices"] = SplitVertices ; - m_CollectionMaps->collectionMap_Vertex["SplitVertices"] = v_cal ; - for(unsigned int i=0 ; i< SplitVertices->size(); i++) m_CollectionMaps->collectionMap_Vertex ["SplitVertices"].push_back(SplitVertices->at(i)); - } - if (NULL != V0Vertices ) - { - std::vector<edm4hep::Vertex> v_cal; - m_CollectionMaps->CollectionMap_Vertex["V0Vertices"] = V0Vertices ; - m_CollectionMaps->collectionMap_Vertex["V0Vertices"] = v_cal ; - for(unsigned int i=0 ; i< V0Vertices->size(); i++) m_CollectionMaps->collectionMap_Vertex ["V0Vertices"].push_back(V0Vertices->at(i)); - } - if (NULL != MarlinTrkTracks ) - { - std::vector<edm4hep::Track> v_cal; - m_CollectionMaps->CollectionMap_Track["MarlinTrkTracks"] = MarlinTrkTracks ; - m_CollectionMaps->collectionMap_Track["MarlinTrkTracks"] = v_cal ; - for(unsigned int i=0 ; i< MarlinTrkTracks->size(); i++) m_CollectionMaps->collectionMap_Track ["MarlinTrkTracks"].push_back(MarlinTrkTracks->at(i)); - } - if (NULL != mcRecoCaloAssociation ) - { - std::vector<edm4hep::MCRecoCaloAssociation> v_cal; - m_CollectionMaps->collectionMap_CaloRel["RecoCaloAssociation_ECALBarrel"] = v_cal ; - for(unsigned int i=0 ; i< mcRecoCaloAssociation->size(); i++) m_CollectionMaps->collectionMap_CaloRel ["RecoCaloAssociation_ECALBarrel"].push_back(mcRecoCaloAssociation->at(i)); - } - - else - { - if (NULL != MCParticle ) - { - for(unsigned int i=0 ; i< MCParticle->size(); i++) - { - if(MCParticle->at(i).parents_size()==0) - { - std::cout<<"create recoCaloAssociation by hand now"<<std::endl; - for(std::map<std::string, std::vector<edm4hep::CalorimeterHit> >::iterator iter = m_CollectionMaps->collectionMap_CaloHit.begin(); iter != m_CollectionMaps->collectionMap_CaloHit.end(); iter++) - { - std::string prefix = "RecoCaloAssociation_"; - std::string key = prefix + iter->first; - std::cout<<"create for "<<key<<std::endl; - std::vector<edm4hep::MCRecoCaloAssociation> v_cal; - m_CollectionMaps->collectionMap_CaloRel[key] = v_cal ; - for(std::vector<edm4hep::CalorimeterHit>::iterator it=iter->second.begin(); it != iter->second.end(); it ++) - { - edm4hep::SimCalorimeterHit sim_hit( it->getCellID(), it->getEnergy(), it->getPosition() ); - edm4hep::CaloHitContribution conb ( MCParticle->at(i).getPDG(), it->getEnergy(), 0, it->getPosition() ); - conb.setParticle( MCParticle->at(i) ); - sim_hit.addToContributions(conb); - edm4hep::MCRecoCaloAssociation calo_association; - calo_association.setRec(*it); - calo_association.setSim(sim_hit); - m_CollectionMaps->collectionMap_CaloRel[key].push_back(calo_association); - } - } - break; - } + for(auto &v : m_dataHandles){ + try{ + if(m_collections[v.first]=="MCParticle"){ + auto handle = dynamic_cast<DataHandle<edm4hep::MCParticleCollection>*> (v.second); + auto po = handle->get(); + if(po != NULL){ + std::vector<edm4hep::MCParticle> v_mc; + m_CollectionMaps->collectionMap_MC [v.first] = v_mc; + for(unsigned int i=0 ; i< po->size(); i++) m_CollectionMaps->collectionMap_MC [v.first].push_back(po->at(i)); + std::cout<<"saved col name="<<v.first<<std::endl; + } + else{ + std::cout<<"don't find col name="<<v.first<<std::endl; } } + else if(m_collections[v.first]=="CalorimeterHit"){ + auto handle = dynamic_cast<DataHandle<edm4hep::CalorimeterHitCollection>*> (v.second); + auto po = handle->get(); + if(po != NULL){ + std::vector<edm4hep::CalorimeterHit> v_cal; + m_CollectionMaps->collectionMap_CaloHit[v.first] = v_cal ; + for(unsigned int i=0 ; i< po->size(); i++) m_CollectionMaps->collectionMap_CaloHit [v.first].push_back(po->at(i)); + std::cout<<"saved col name="<<v.first<<std::endl; + } + else{ + std::cout<<"don't find col name="<<v.first<<std::endl; + } + } + else if(m_collections[v.first]=="Track"){ + auto handle = dynamic_cast<DataHandle<edm4hep::TrackCollection>*> (v.second); + auto po = handle->get(); + if(po != NULL){ + std::vector<edm4hep::Track> v_cal; + m_CollectionMaps->collectionMap_Track[v.first] = v_cal ; + for(unsigned int i=0 ; i< po->size(); i++) m_CollectionMaps->collectionMap_Track [v.first].push_back(po->at(i)); + std::cout<<"saved col name="<<v.first<<std::endl; + } + else{ + std::cout<<"don't find col name="<<v.first<<std::endl; + } + } + else if(m_collections[v.first]=="Vertex"){ + auto handle = dynamic_cast<DataHandle<edm4hep::VertexCollection>*> (v.second); + auto po = handle->get(); + if(po != NULL){ + std::vector<edm4hep::Vertex> v_cal; + m_CollectionMaps->collectionMap_Vertex[v.first] = v_cal ; + for(unsigned int i=0 ; i< po->size(); i++) m_CollectionMaps->collectionMap_Vertex [v.first].push_back(po->at(i)); + std::cout<<"saved col name="<<v.first<<std::endl; + } + else{ + std::cout<<"don't find col name="<<v.first<<std::endl; + } + } + else if(m_collections[v.first]=="MCRecoCaloAssociation"){ + auto handle = dynamic_cast<DataHandle<edm4hep::MCRecoCaloAssociationCollection>*> (v.second); + auto po = handle->get(); + if(po != NULL){ + std::vector<edm4hep::MCRecoCaloAssociation> v_cal; + m_CollectionMaps->collectionMap_CaloRel[v.first] = v_cal ; + for(unsigned int i=0 ; i< po->size(); i++) m_CollectionMaps->collectionMap_CaloRel [v.first].push_back(po->at(i)); + std::cout<<"saved col name="<<v.first<<std::endl; + } + else{ + std::cout<<"don't find col name="<<v.first<<std::endl; + } + } + else if(m_collections[v.first]=="MCRecoTrackerAssociation"){ + auto handle = dynamic_cast<DataHandle<edm4hep::MCRecoTrackerAssociationCollection>*> (v.second); + auto po = handle->get(); + if(po != NULL){ + std::vector<edm4hep::MCRecoTrackerAssociation> v_cal; + m_CollectionMaps->collectionMap_TrkRel[v.first] = v_cal ; + for(unsigned int i=0 ; i< po->size(); i++) m_CollectionMaps->collectionMap_TrkRel [v.first].push_back(po->at(i)); + std::cout<<"saved col name="<<v.first<<std::endl; + } + else{ + std::cout<<"don't find col name="<<v.first<<std::endl; + } + } + else{ + std::cout<<"wrong type name for col :"<<v.first<<std::endl; + } + }//try + catch(...){ + std::cout<<"don't find "<<v.first<<"in event"<<std::endl; + std::cout<<"don't find col name="<<v.first<<",with type="<<m_collections[v.first]<<" in this event"<<std::endl; } - - if (NULL != mcRecoTrackerAssociation ) - { - std::vector<edm4hep::MCRecoTrackerAssociation> v_cal; - m_CollectionMaps->collectionMap_TrkRel["RecoTrackerAssociation"] = v_cal ; - for(unsigned int i=0 ; i< mcRecoTrackerAssociation->size(); i++) m_CollectionMaps->collectionMap_TrkRel ["RecoTrackerAssociation"].push_back(mcRecoTrackerAssociation->at(i)); - } - return StatusCode::SUCCESS; -} - - -StatusCode PandoraMatrixAlg::updateMap(CollectionMaps & tmp_map) -{ - const edm4hep::MCParticleCollection* MCParticle = nullptr; - const edm4hep::CalorimeterHitCollection* ECALBarrel = nullptr; - const edm4hep::CalorimeterHitCollection* ECALEndcap = nullptr; - const edm4hep::CalorimeterHitCollection* ECALOther = nullptr; - const edm4hep::CalorimeterHitCollection* HCALBarrel = nullptr; - const edm4hep::CalorimeterHitCollection* HCALEndcap = nullptr; - const edm4hep::CalorimeterHitCollection* HCALOther = nullptr; - const edm4hep::CalorimeterHitCollection* MUON = nullptr; - const edm4hep::CalorimeterHitCollection* LCAL = nullptr; - const edm4hep::CalorimeterHitCollection* LHCAL = nullptr; - const edm4hep::CalorimeterHitCollection* BCAL = nullptr; - const edm4hep::VertexCollection* KinkVertices = nullptr; - const edm4hep::VertexCollection* ProngVertices = nullptr; - const edm4hep::VertexCollection* SplitVertices = nullptr; - const edm4hep::VertexCollection* V0Vertices = nullptr; - const edm4hep::TrackCollection* MarlinTrkTracks = nullptr; - const edm4hep::MCRecoCaloAssociationCollection* mcRecoCaloAssociation = nullptr; - const edm4hep::MCRecoTrackerAssociationCollection* mcRecoTrackerAssociation = nullptr; - StatusCode sc = StatusCode::SUCCESS; - sc = getCol(m_mcParCol_r , MCParticle ); - sc = getCol(m_ECALBarrel_r, ECALBarrel ); - sc = getCol(m_ECALEndcap_r, ECALEndcap ); - sc = getCol(m_ECALOther_r , ECALOther ); - sc = getCol(m_HCALBarrel_r, HCALBarrel ); - sc = getCol(m_HCALEndcap_r, HCALEndcap ); - sc = getCol(m_HCALOther_r , HCALOther ); - sc = getCol(m_MUON_r , MUON ); - sc = getCol(m_LCAL_r , LCAL ); - sc = getCol(m_LHCAL_r , LHCAL ); - sc = getCol(m_BCAL_r , BCAL ); - sc = getCol(m_KinkVertices_r , KinkVertices ); - sc = getCol(m_ProngVertices_r , ProngVertices); - sc = getCol(m_SplitVertices_r , SplitVertices); - sc = getCol(m_V0Vertices_r , V0Vertices ); - sc = getCol(m_MarlinTrkTracks_r , MarlinTrkTracks ); - sc = getCol(m_MCRecoCaloAssociation_r , mcRecoCaloAssociation ); - sc = getCol(m_MCRecoTrackerAssociation_r , mcRecoTrackerAssociation ); - - if (NULL != MCParticle ) - { - std::vector<edm4hep::MCParticle> v_mc; - tmp_map.CollectionMap_MC ["MCParticle"] = MCParticle ; - tmp_map.collectionMap_MC ["MCParticle"] = v_mc; - for(unsigned int i=0 ; i< MCParticle->size(); i++) tmp_map.collectionMap_MC ["MCParticle"].push_back(MCParticle->at(i)); - } - if (NULL != ECALBarrel ) - { - std::vector<edm4hep::CalorimeterHit> v_cal; - tmp_map.CollectionMap_CaloHit["ECALBarrel"] = ECALBarrel ; - tmp_map.collectionMap_CaloHit["ECALBarrel"] = v_cal ; - for(unsigned int i=0 ; i< ECALBarrel->size(); i++) tmp_map.collectionMap_CaloHit ["ECALBarrel"].push_back(ECALBarrel->at(i)); - } - if (NULL != ECALEndcap ) - { - std::vector<edm4hep::CalorimeterHit> v_cal; - tmp_map.CollectionMap_CaloHit["ECALEndcap"] = ECALEndcap ; - tmp_map.collectionMap_CaloHit["ECALEndcap"] = v_cal ; - for(unsigned int i=0 ; i< ECALEndcap->size(); i++) tmp_map.collectionMap_CaloHit ["ECALEndcap"].push_back(ECALEndcap->at(i)); - } - if (NULL != ECALOther ) - { - std::vector<edm4hep::CalorimeterHit> v_cal; - tmp_map.CollectionMap_CaloHit["ECALOther"] = ECALOther ; - tmp_map.collectionMap_CaloHit["ECALOther"] = v_cal ; - for(unsigned int i=0 ; i< ECALOther->size(); i++) tmp_map.collectionMap_CaloHit ["ECALOther"].push_back(ECALOther->at(i)); - } - if (NULL != HCALBarrel ) - { - std::vector<edm4hep::CalorimeterHit> v_cal; - tmp_map.CollectionMap_CaloHit["HCALBarrel"] = HCALBarrel ; - tmp_map.collectionMap_CaloHit["HCALBarrel"] = v_cal ; - for(unsigned int i=0 ; i< HCALBarrel->size(); i++) tmp_map.collectionMap_CaloHit ["HCALBarrel"].push_back(HCALBarrel->at(i)); - } - if (NULL != HCALEndcap ) - { - std::vector<edm4hep::CalorimeterHit> v_cal; - tmp_map.CollectionMap_CaloHit["HCALEndcap"] = HCALEndcap ; - tmp_map.collectionMap_CaloHit["HCALEndcap"] = v_cal ; - for(unsigned int i=0 ; i< HCALEndcap->size(); i++) tmp_map.collectionMap_CaloHit ["HCALEndcap"].push_back(HCALEndcap->at(i)); - } - if (NULL != HCALOther ) - { - std::vector<edm4hep::CalorimeterHit> v_cal; - tmp_map.CollectionMap_CaloHit["HCALOther"] = HCALOther ; - tmp_map.collectionMap_CaloHit["HCALOther"] = v_cal ; - for(unsigned int i=0 ; i< HCALOther->size(); i++) tmp_map.collectionMap_CaloHit ["HCALOther"].push_back(HCALOther->at(i)); - } - if (NULL != MUON ) - { - std::vector<edm4hep::CalorimeterHit> v_cal; - tmp_map.CollectionMap_CaloHit["MUON"] = MUON ; - tmp_map.collectionMap_CaloHit["MUON"] = v_cal ; - for(unsigned int i=0 ; i< MUON->size(); i++) tmp_map.collectionMap_CaloHit ["MUON"].push_back(MUON->at(i)); - } - if (NULL != LCAL ) - { - std::vector<edm4hep::CalorimeterHit> v_cal; - tmp_map.CollectionMap_CaloHit["LCAL"] = LCAL ; - tmp_map.collectionMap_CaloHit["LCAL"] = v_cal ; - for(unsigned int i=0 ; i< LCAL->size(); i++) tmp_map.collectionMap_CaloHit ["LCAL"].push_back(LCAL->at(i)); - } - if (NULL != LHCAL ) - { - std::vector<edm4hep::CalorimeterHit> v_cal; - tmp_map.CollectionMap_CaloHit["LHCAL"] = LHCAL ; - tmp_map.collectionMap_CaloHit["LHCAL"] = v_cal ; - for(unsigned int i=0 ; i< LHCAL->size(); i++) tmp_map.collectionMap_CaloHit ["LHCAL"].push_back(LHCAL->at(i)); - } - if (NULL != BCAL ) - { - std::vector<edm4hep::CalorimeterHit> v_cal; - tmp_map.CollectionMap_CaloHit["BCAL"] = BCAL ; - tmp_map.collectionMap_CaloHit["BCAL"] = v_cal ; - for(unsigned int i=0 ; i< BCAL->size(); i++) tmp_map.collectionMap_CaloHit ["BCAL"].push_back(BCAL->at(i)); - } - if (NULL != KinkVertices ) - { - std::vector<edm4hep::Vertex> v_cal; - tmp_map.CollectionMap_Vertex["KinkVertices"] = KinkVertices ; - tmp_map.collectionMap_Vertex["KinkVertices"] = v_cal ; - for(unsigned int i=0 ; i< KinkVertices->size(); i++) tmp_map.collectionMap_Vertex ["KinkVertices"].push_back(KinkVertices->at(i)); - } - if (NULL != ProngVertices ) - { - std::vector<edm4hep::Vertex> v_cal; - tmp_map.CollectionMap_Vertex["ProngVertices"] = ProngVertices ; - tmp_map.collectionMap_Vertex["ProngVertices"] = v_cal ; - for(unsigned int i=0 ; i< ProngVertices->size(); i++) tmp_map.collectionMap_Vertex ["ProngVertices"].push_back(ProngVertices->at(i)); - } - if (NULL != SplitVertices ) - { - std::vector<edm4hep::Vertex> v_cal; - tmp_map.CollectionMap_Vertex["SplitVertices"] = SplitVertices ; - tmp_map.collectionMap_Vertex["SplitVertices"] = v_cal ; - for(unsigned int i=0 ; i< SplitVertices->size(); i++) tmp_map.collectionMap_Vertex ["SplitVertices"].push_back(SplitVertices->at(i)); - } - if (NULL != V0Vertices ) - { - std::vector<edm4hep::Vertex> v_cal; - tmp_map.CollectionMap_Vertex["V0Vertices"] = V0Vertices ; - tmp_map.collectionMap_Vertex["V0Vertices"] = v_cal ; - for(unsigned int i=0 ; i< V0Vertices->size(); i++) tmp_map.collectionMap_Vertex ["V0Vertices"].push_back(V0Vertices->at(i)); - } - if (NULL != MarlinTrkTracks ) - { - std::vector<edm4hep::Track> v_cal; - tmp_map.CollectionMap_Track["MarlinTrkTracks"] = MarlinTrkTracks ; - tmp_map.collectionMap_Track["MarlinTrkTracks"] = v_cal ; - for(unsigned int i=0 ; i< MarlinTrkTracks->size(); i++) tmp_map.collectionMap_Track ["MarlinTrkTracks"].push_back(MarlinTrkTracks->at(i)); - } - if (NULL != mcRecoCaloAssociation ) - { - std::vector<edm4hep::MCRecoCaloAssociation> v_cal; - tmp_map.collectionMap_CaloRel["RecoCaloAssociation"] = v_cal ; - for(unsigned int i=0 ; i< mcRecoCaloAssociation->size(); i++) tmp_map.collectionMap_CaloRel ["RecoCaloAssociation"].push_back(mcRecoCaloAssociation->at(i)); - } - if (NULL != mcRecoTrackerAssociation ) - { - std::vector<edm4hep::MCRecoTrackerAssociation> v_cal; - tmp_map.collectionMap_TrkRel["RecoTrackerAssociation"] = v_cal ; - for(unsigned int i=0 ; i< mcRecoTrackerAssociation->size(); i++) tmp_map.collectionMap_TrkRel ["RecoTrackerAssociation"].push_back(mcRecoTrackerAssociation->at(i)); - } + } return StatusCode::SUCCESS; } @@ -887,19 +626,6 @@ StatusCode PandoraMatrixAlg::Ana() if(hasEm && hasEp) m_hasConversion=1; } } - //sanity check calo info - const edm4hep::CalorimeterHitCollection* CaloCol = nullptr; - sc = getCol(m_ECALBarrel_r, CaloCol); - if (NULL != CaloCol ) - { - for(unsigned int i=0 ; i< CaloCol->size(); i++) - { - m_hits_x.push_back(CaloCol->at(i).getPosition()[0]); - m_hits_y.push_back(CaloCol->at(i).getPosition()[1]); - m_hits_z.push_back(CaloCol->at(i).getPosition()[2]); - m_hits_E.push_back(CaloCol->at(i).getEnergy() ); - } - } m_tree->Fill(); return StatusCode::SUCCESS; diff --git a/Reconstruction/PFA/Pandora/MatrixPandora/src/PfoCreator.cpp b/Reconstruction/PFA/Pandora/MatrixPandora/src/PfoCreator.cpp index b197f2dbb907b3c98597c5e7acf9b700208b155e..c20c597faa4f59a440718ed4e874576038b6f294 100644 --- a/Reconstruction/PFA/Pandora/MatrixPandora/src/PfoCreator.cpp +++ b/Reconstruction/PFA/Pandora/MatrixPandora/src/PfoCreator.cpp @@ -1,12 +1,10 @@ /** - * @file MarlinPandora/src/PfoCreator.cc * * @brief Implementation of the pfo creator class. * * $Log: $ */ -//#include "CalorimeterHitType.h" #include "Api/PandoraApi.h" @@ -46,27 +44,20 @@ pandora::StatusCode PfoCreator::CreateParticleFlowObjects(CollectionMaps& collec PANDORA_RETURN_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, PandoraApi::GetCurrentPfoList(*m_pPandora, pPandoraPfoList)); - //IMPL::LCFlagImpl lcFlagImpl(pClusterCollection->getFlag()); - //lcFlagImpl.setBit(LCIO::CLBIT_HITS); - //pClusterCollection->setFlag(lcFlagImpl.getFlag()); pandora::StringVector subDetectorNames; this->InitialiseSubDetectorNames(subDetectorNames); - //pClusterCollection->parameters().setValues("ClusterSubdetectorNames", subDetectorNames); - // Create lcio "reconstructed particles" from the pandora "particle flow objects" std::cout<<"pPandoraPfoList size="<<pPandoraPfoList->size()<<std::endl; for (pandora::PfoList::const_iterator pIter = pPandoraPfoList->begin(), pIterEnd = pPandoraPfoList->end(); pIter != pIterEnd; ++pIter) { const pandora::ParticleFlowObject *const pPandoraPfo(*pIter); - //IMPL::ReconstructedParticleImpl *const pReconstructedParticle(new ReconstructedParticleImpl()); edm4hep::ReconstructedParticle pReconstructedParticle0 = pReconstructedParticleCollection->create(); edm4hep::ReconstructedParticle* pReconstructedParticle = &pReconstructedParticle0; const bool hasTrack(!pPandoraPfo->GetTrackList().empty()); const pandora::ClusterList &clusterList(pPandoraPfo->GetClusterList()); - //std::cout<<"ClusterList size="<<clusterList.size()<<std::endl; float clustersTotalEnergy(0.f); pandora::CartesianVector referencePoint(0.f, 0.f, 0.f), clustersWeightedPosition(0.f, 0.f, 0.f); for (pandora::ClusterList::const_iterator cIter = clusterList.begin(), cIterEnd = clusterList.end(); cIter != cIterEnd; ++cIter) @@ -77,7 +68,6 @@ pandora::StatusCode PfoCreator::CreateParticleFlowObjects(CollectionMaps& collec pandoraCaloHitList.insert(pandoraCaloHitList.end(), pPandoraCluster->GetIsolatedCaloHitList().begin(), pPandoraCluster->GetIsolatedCaloHitList().end()); pandora::FloatVector hitE, hitX, hitY, hitZ; - //IMPL::ClusterImpl *const p_Cluster(new ClusterImpl()); edm4hep::Cluster p_Cluster0 = pClusterCollection->create(); edm4hep::Cluster* p_Cluster = &p_Cluster0; this->SetClusterSubDetectorEnergies(subDetectorNames, p_Cluster, pandoraCaloHitList, hitE, hitX, hitY, hitZ); @@ -95,7 +85,6 @@ pandora::StatusCode PfoCreator::CreateParticleFlowObjects(CollectionMaps& collec clustersTotalEnergy += clusterCorrectEnergy; } - //pClusterCollection->addElement(p_Cluster); edm4hep::ConstCluster p_ClusterCon = *p_Cluster; pReconstructedParticle->addToClusters(p_ClusterCon); } @@ -123,7 +112,6 @@ pandora::StatusCode PfoCreator::CreateParticleFlowObjects(CollectionMaps& collec edm4hep::Vertex pStartVertex0 = pStartVertexCollection->create(); edm4hep::Vertex* pStartVertex = &pStartVertex0; - //pStartVertex->setAlgorithmType(m_settings.m_startVertexAlgName.c_str()); pStartVertex->setAlgorithmType(0); const float ref_value[3] = {referencePoint.GetX(),referencePoint.GetY(),referencePoint.GetZ()}; pStartVertex->setPosition(edm4hep::Vector3f(ref_value)); @@ -260,7 +248,6 @@ pandora::StatusCode PfoCreator::CalculateTrackBasedReferencePoint(const pandora: const float z0(pPandoraTrack->GetZ0()); pandora::CartesianVector intersectionPoint(0.f, 0.f, 0.f); - //intersectionPoint.SetValues(pLcioTrack->getD0() * std::cos(pLcioTrack->getPhi()), pLcioTrack->getD0() * std::sin(pLcioTrack->getPhi()), z0); if(pLcioTrack.trackStates_size()==0) throw "zero trackStates size find"; intersectionPoint.SetValues(pLcioTrack.getTrackStates(0).D0 * std::cos(pLcioTrack.getTrackStates(0).phi), pLcioTrack.getTrackStates(0).D0 * std::sin(pLcioTrack.getTrackStates(0).phi), z0); const float trackMomentumAtDca((pPandoraTrack->GetMomentumAtDca()).GetMagnitude()); diff --git a/Reconstruction/PFA/Pandora/MatrixPandora/src/TrackCreator.cpp b/Reconstruction/PFA/Pandora/MatrixPandora/src/TrackCreator.cpp index 1a3c0f15a7107e3671d246b80d69918f82c952de..2a5c8f0084f15dcb3f22a22e97a67081db215844 100644 --- a/Reconstruction/PFA/Pandora/MatrixPandora/src/TrackCreator.cpp +++ b/Reconstruction/PFA/Pandora/MatrixPandora/src/TrackCreator.cpp @@ -1,5 +1,4 @@ /** - * @file MarlinPandora/src/TrackCreator.cc * * @brief Implementation of the track creator class. * @@ -141,67 +140,44 @@ TrackCreator::~TrackCreator() pandora::StatusCode TrackCreator::CreateTrackAssociations(const CollectionMaps& collectionMaps) { -// Don't use it now, because the Vertex.getAssociatedParticle() doesn't work for LCIO to plcio transfer -// PANDORA_RETURN_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, this->ExtractKinks(collectionMaps)); -// PANDORA_RETURN_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, this->ExtractProngsAndSplits(const CollectionMaps& collectionMaps)); -// PANDORA_RETURN_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, this->ExtractV0s(const CollectionMaps& collectionMaps)); + PANDORA_RETURN_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, this->ExtractKinks(collectionMaps)); + PANDORA_RETURN_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, this->ExtractProngsAndSplits(collectionMaps)); + PANDORA_RETURN_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, this->ExtractV0s(collectionMaps)); return pandora::STATUS_CODE_SUCCESS; } //------------------------------------------------------------------------------------------------------------------------------------------ -/* pandora::StatusCode TrackCreator::ExtractKinks(const CollectionMaps& collectionMaps) { std::cout<<"start TrackCreator::ExtractKinks:"<<std::endl; for (StringVector::const_iterator iter = m_settings.m_kinkVertexCollections.begin(), iterEnd = m_settings.m_kinkVertexCollections.end(); iter != iterEnd; ++iter) { - if(collectionMaps.CollectionMap_Vertex.find(*iter) == collectionMaps.CollectionMap_Vertex.end()) { std::cout<<"not find "<<(*iter)<<std::endl; continue;} + if(collectionMaps.collectionMap_Vertex.find(*iter) == collectionMaps.collectionMap_Vertex.end()) { std::cout<<"not find "<<(*iter)<<std::endl; continue;} try { - const edm4hep::VertexCollection *pKinkCollection = (collectionMaps.CollectionMap_Vertex.find(*iter))->second; + const std::vector<edm4hep::Vertex>& pKinkCollection = (collectionMaps.collectionMap_Vertex.find(*iter))->second; - for (int i = 0, iMax = pKinkCollection->size(); i < iMax; ++i) + for (int i = 0, iMax = pKinkCollection.size(); i < iMax; ++i) { try { - const edm4hep::Vertex pVertex0 = pKinkCollection->at(i); + const edm4hep::Vertex pVertex0 = pKinkCollection.at(i); const edm4hep::Vertex* pVertex = &(pVertex0); if (NULL == pVertex) throw ("Collection type mismatch"); - std::cout<<"pVertex0 getChi2="<<pVertex0.getChi2()<<std::endl; - std::cout<<"pVertex getChi2="<<pVertex->getChi2()<<std::endl; - std::cout<<"Hi 0 "<<std::endl; - const plcio::ConstReconstructedParticle pReconstructedParticle0 = pVertex0.getAssociatedParticle(); - std::cout<<"Hi 1"<<std::endl; - //EVENT::ReconstructedParticle *pReconstructedParticle = pVertex->getAssociatedParticle(); - const plcio::ConstReconstructedParticle pReconstructedParticle = pVertex->getAssociatedParticle(); - std::cout<<"Hi 2:"<<&pReconstructedParticle<<std::endl; - //const EVENT::TrackVec &trackVec(pReconstructedParticle->getTracks()); - //plcio::ConstTrack trackVec = pReconstructedParticle.getTracks(); - - std::cout<<"pReconstructedParticle0 en="<<pReconstructedParticle0.getEnergy()<<std::endl; - std::cout<<"pReconstructedParticle en="<<pReconstructedParticle.getEnergy()<<std::endl; - //if (this->IsConflictingRelationship(trackVec))continue; + const edm4hep::ConstReconstructedParticle pReconstructedParticle = pVertex->getAssociatedParticle(); if (this->IsConflictingRelationship(pReconstructedParticle))continue; - //const int vertexPdgCode(pReconstructedParticle->getType()); const int vertexPdgCode(pReconstructedParticle.getType()); // Extract the kink vertex information - //for (unsigned int iTrack = 0, nTracks = trackVec.size(); iTrack < nTracks; ++iTrack) - //for (unsigned int iTrack = 0, nTracks = trackVec.tracks_size(); iTrack < nTracks; ++iTrack) for (unsigned int iTrack = 0, nTracks = pReconstructedParticle.tracks_size(); iTrack < nTracks; ++iTrack) { - //EVENT::Track *pTrack = trackVec[iTrack]; - //plcio::ConstTrack pTrack = trackVec.getTracks(iTrack); - plcio::ConstTrack pTrack = pReconstructedParticle.getTracks(iTrack); - //(0 == iTrack) ? m_parentTrackList.insert(pTrack) : m_daughterTrackList.insert(pTrack); + edm4hep::ConstTrack pTrack = pReconstructedParticle.getTracks(iTrack); (0 == iTrack) ? m_parentTrackList.insert(pTrack.id()) : m_daughterTrackList.insert(pTrack.id()); - //std::cout << "KinkTrack " << iTrack << ", nHits " << pTrack.getTrackerHits().size() << std::endl; - std::cout << "KinkTrack " << iTrack << ", nHits " << pTrack.trackerHits_size() << std::endl; int trackPdgCode = pandora::UNKNOWN_PARTICLE_TYPE; @@ -230,7 +206,6 @@ pandora::StatusCode TrackCreator::ExtractKinks(const CollectionMaps& collectionM trackPdgCode = pandora::PI_PLUS; break; default : - //(pTrack->getOmega() > 0) ? trackPdgCode = pandora::PI_PLUS : trackPdgCode = pandora::PI_MINUS; (pTrack.getTrackStates(0).omega > 0) ? trackPdgCode = pandora::PI_PLUS : trackPdgCode = pandora::PI_MINUS; break; } @@ -246,9 +221,7 @@ pandora::StatusCode TrackCreator::ExtractKinks(const CollectionMaps& collectionM { for (unsigned int jTrack = iTrack + 1; jTrack < nTracks; ++jTrack) { - //PANDORA_RETURN_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, PandoraApi::SetTrackParentDaughterRelationship(*m_pPandora, pTrack, trackVec[jTrack])); - //PANDORA_RETURN_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, PandoraApi::SetTrackParentDaughterRelationship(*m_pPandora, (int*)(pTrack.id()), (int*)(trackVec.getTracks(jTrack).id()) ) ); - PANDORA_RETURN_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, PandoraApi::SetTrackParentDaughterRelationship(*m_pPandora, (int*)(pTrack.id()), (int*)(pReconstructedParticle.getTracks(jTrack).id()) ) ); + PANDORA_RETURN_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, PandoraApi::SetTrackParentDaughterRelationship(*m_pPandora, GetTrackAddress(collectionMaps, pTrack), GetTrackAddress(collectionMaps, pReconstructedParticle.getTracks(jTrack) ) ) ); } } @@ -257,9 +230,7 @@ pandora::StatusCode TrackCreator::ExtractKinks(const CollectionMaps& collectionM { for (unsigned int jTrack = iTrack + 1; jTrack < nTracks; ++jTrack) { - //PANDORA_RETURN_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, PandoraApi::SetTrackSiblingRelationship(*m_pPandora, pTrack, trackVec[jTrack])); - //PANDORA_RETURN_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, PandoraApi::SetTrackSiblingRelationship(*m_pPandora, (int*)(pTrack.id()), (int*)(trackVec.getTracks(jTrack).id()) )); - PANDORA_RETURN_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, PandoraApi::SetTrackSiblingRelationship(*m_pPandora, (int*)(pTrack.id()), (int*)(pReconstructedParticle.getTracks(jTrack).id()) )); + PANDORA_RETURN_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, PandoraApi::SetTrackSiblingRelationship(*m_pPandora, GetTrackAddress(collectionMaps, pTrack), GetTrackAddress(collectionMaps, pReconstructedParticle.getTracks(jTrack) ) ) ); } } } @@ -278,50 +249,44 @@ pandora::StatusCode TrackCreator::ExtractKinks(const CollectionMaps& collectionM return pandora::STATUS_CODE_SUCCESS; } -*/ //------------------------------------------------------------------------------------------------------------------------------------------ -/* -pandora::StatusCode TrackCreator::ExtractProngsAndSplits(const EVENT::LCEvent *const pLCEvent) + +pandora::StatusCode TrackCreator::ExtractProngsAndSplits(const CollectionMaps& collectionMaps) { - for (StringVector::const_iterator iter = m_settings.m_prongSplitVertexCollections.begin(), iterEnd = m_settings.m_prongSplitVertexCollections.end(); - iter != iterEnd; ++iter) + std::cout<<"start TrackCreator::ExtractProngsAndSplits:"<<std::endl; + for (StringVector::const_iterator iter = m_settings.m_prongSplitVertexCollections.begin(), iterEnd = m_settings.m_prongSplitVertexCollections.end(); iter != iterEnd; ++iter) { + if(collectionMaps.collectionMap_Vertex.find(*iter) == collectionMaps.collectionMap_Vertex.end()) { std::cout<<"not find "<<(*iter)<<std::endl; continue;} try { - const EVENT::LCCollection *pProngOrSplitCollection = pLCEvent->getCollection(*iter); + const std::vector<edm4hep::Vertex>& pProngOrSplitCollection = (collectionMaps.collectionMap_Vertex.find(*iter))->second; - for (int i = 0, iMax = pProngOrSplitCollection->getNumberOfElements(); i < iMax; ++i) + for (int i = 0, iMax = pProngOrSplitCollection.size(); i < iMax; ++i) { try { - EVENT::Vertex *pVertex = dynamic_cast<Vertex*>(pProngOrSplitCollection->getElementAt(i)); - - if (NULL == pVertex) - throw EVENT::Exception("Collection type mismatch"); + const edm4hep::Vertex pVertex0 = pProngOrSplitCollection.at(i); + const edm4hep::Vertex* pVertex = &(pVertex0); - EVENT::ReconstructedParticle *pReconstructedParticle = pVertex->getAssociatedParticle(); - const EVENT::TrackVec &trackVec(pReconstructedParticle->getTracks()); + if (NULL == pVertex) throw ("Collection type mismatch"); + const edm4hep::ConstReconstructedParticle pReconstructedParticle = pVertex->getAssociatedParticle(); - if (this->IsConflictingRelationship(trackVec)) - continue; + if (this->IsConflictingRelationship(pReconstructedParticle))continue; // Extract the prong/split vertex information - for (unsigned int iTrack = 0, nTracks = trackVec.size(); iTrack < nTracks; ++iTrack) + for (unsigned int iTrack = 0, nTracks = pReconstructedParticle.tracks_size(); iTrack < nTracks; ++iTrack) { - EVENT::Track *pTrack = trackVec[iTrack]; - (0 == iTrack) ? m_parentTrackList.insert(pTrack) : m_daughterTrackList.insert(pTrack); - streamlog_out(DEBUG) << "Prong or Split Track " << iTrack << ", nHits " << pTrack->getTrackerHits().size() << std::endl; + edm4hep::ConstTrack pTrack = pReconstructedParticle.getTracks(iTrack); + (0 == iTrack) ? m_parentTrackList.insert(pTrack.id()) : m_daughterTrackList.insert(pTrack.id()); - if (0 == m_settings.m_shouldFormTrackRelationships) - continue; + if (0 == m_settings.m_shouldFormTrackRelationships) continue; // Make track parent-daughter relationships if (0 == iTrack) { for (unsigned int jTrack = iTrack + 1; jTrack < nTracks; ++jTrack) { - PANDORA_RETURN_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, PandoraApi::SetTrackParentDaughterRelationship(*m_pPandora, - pTrack, trackVec[jTrack])); + PANDORA_RETURN_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, PandoraApi::SetTrackParentDaughterRelationship(*m_pPandora, GetTrackAddress(collectionMaps, pTrack), GetTrackAddress(collectionMaps, pReconstructedParticle.getTracks(jTrack) ) ) ); } } @@ -330,21 +295,20 @@ pandora::StatusCode TrackCreator::ExtractProngsAndSplits(const EVENT::LCEvent *c { for (unsigned int jTrack = iTrack + 1; jTrack < nTracks; ++jTrack) { - PANDORA_RETURN_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, PandoraApi::SetTrackSiblingRelationship(*m_pPandora, - pTrack, trackVec[jTrack])); + PANDORA_RETURN_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, PandoraApi::SetTrackSiblingRelationship(*m_pPandora, GetTrackAddress(collectionMaps, pTrack), GetTrackAddress(collectionMaps, pReconstructedParticle.getTracks(jTrack) ) ) ); } } } } - catch (EVENT::Exception &exception) + catch (...) { - streamlog_out(WARNING) << "Failed to extract prong/split vertex: " << exception.what() << std::endl; + std::cout << "Failed to extract prong/split vertex" <<std::endl; } } } - catch (EVENT::Exception &exception) + catch (...) { - streamlog_out(DEBUG5) << "Failed to extract prong/split vertex collection: " << *iter << ", " << exception.what() << std::endl; + std::cout<< "Failed to extract prong/split vertex collection: " << *iter << std::endl; } } @@ -353,102 +317,92 @@ pandora::StatusCode TrackCreator::ExtractProngsAndSplits(const EVENT::LCEvent *c //------------------------------------------------------------------------------------------------------------------------------------------ -pandora::StatusCode TrackCreator::ExtractV0s(const EVENT::LCEvent *const pLCEvent) +pandora::StatusCode TrackCreator::ExtractV0s(const CollectionMaps& collectionMaps) { - for (StringVector::const_iterator iter = m_settings.m_v0VertexCollections.begin(), iterEnd = m_settings.m_v0VertexCollections.end(); - iter != iterEnd; ++iter) + std::cout<<"start TrackCreator::ExtractV0s:"<<std::endl; + for (StringVector::const_iterator iter = m_settings.m_v0VertexCollections.begin(), iterEnd = m_settings.m_v0VertexCollections.end(); iter != iterEnd; ++iter) { + if(collectionMaps.collectionMap_Vertex.find(*iter) == collectionMaps.collectionMap_Vertex.end()) { std::cout<<"not find "<<(*iter)<<std::endl; continue;} try { - const EVENT::LCCollection *pV0Collection = pLCEvent->getCollection(*iter); + const std::vector<edm4hep::Vertex>& pV0Collection = (collectionMaps.collectionMap_Vertex.find(*iter))->second; - for (int i = 0, iMax = pV0Collection->getNumberOfElements(); i < iMax; ++i) + for (int i = 0, iMax = pV0Collection.size(); i < iMax; ++i) { try { - EVENT::Vertex *pVertex = dynamic_cast<Vertex*>(pV0Collection->getElementAt(i)); + const edm4hep::Vertex pVertex0 = pV0Collection.at(i); + const edm4hep::Vertex* pVertex = &(pVertex0); - if (NULL == pVertex) - throw EVENT::Exception("Collection type mismatch"); + if (NULL == pVertex) throw ("Collection type mismatch"); - EVENT::ReconstructedParticle *pReconstructedParticle = pVertex->getAssociatedParticle(); - const EVENT::TrackVec &trackVec(pReconstructedParticle->getTracks()); + const edm4hep::ConstReconstructedParticle pReconstructedParticle = pVertex->getAssociatedParticle(); - if (this->IsConflictingRelationship(trackVec)) - continue; + if (this->IsConflictingRelationship(pReconstructedParticle))continue; // Extract the v0 vertex information - const int vertexPdgCode(pReconstructedParticle->getType()); + const int vertexPdgCode(pReconstructedParticle.getType()); - for (unsigned int iTrack = 0, nTracks = trackVec.size(); iTrack < nTracks; ++iTrack) + for (unsigned int iTrack = 0, nTracks = pReconstructedParticle.tracks_size(); iTrack < nTracks; ++iTrack) { - EVENT::Track *pTrack = trackVec[iTrack]; - m_v0TrackList.insert(pTrack); - streamlog_out(DEBUG) << "V0Track " << iTrack << ", nHits " << pTrack->getTrackerHits().size() << std::endl; + edm4hep::ConstTrack pTrack = pReconstructedParticle.getTracks(iTrack); + m_v0TrackList.insert(pTrack.id()); int trackPdgCode = pandora::UNKNOWN_PARTICLE_TYPE; switch (vertexPdgCode) { case pandora::PHOTON : - (pTrack->getOmega() > 0) ? trackPdgCode = pandora::E_PLUS : trackPdgCode = pandora::E_MINUS; + (pTrack.getTrackStates(0).omega > 0) ? trackPdgCode = pandora::E_PLUS : trackPdgCode = pandora::E_MINUS; break; case pandora::LAMBDA : - (pTrack->getOmega() > 0) ? trackPdgCode = pandora::PROTON : trackPdgCode = pandora::PI_MINUS; + (pTrack.getTrackStates(0).omega > 0) ? trackPdgCode = pandora::PROTON : trackPdgCode = pandora::PI_MINUS; break; case pandora::LAMBDA_BAR : - (pTrack->getOmega() > 0) ? trackPdgCode = pandora::PI_PLUS : trackPdgCode = pandora::PROTON_BAR; + (pTrack.getTrackStates(0).omega > 0) ? trackPdgCode = pandora::PI_PLUS : trackPdgCode = pandora::PROTON_BAR; break; case pandora::K_SHORT : - (pTrack->getOmega() > 0) ? trackPdgCode = pandora::PI_PLUS : trackPdgCode = pandora::PI_MINUS; + (pTrack.getTrackStates(0).omega > 0) ? trackPdgCode = pandora::PI_PLUS : trackPdgCode = pandora::PI_MINUS; break; default : - (pTrack->getOmega() > 0) ? trackPdgCode = pandora::PI_PLUS : trackPdgCode = pandora::PI_MINUS; + (pTrack.getTrackStates(0).omega > 0) ? trackPdgCode = pandora::PI_PLUS : trackPdgCode = pandora::PI_MINUS; break; } m_trackToPidMap.insert(TrackToPidMap::value_type(pTrack, trackPdgCode)); - if (0 == m_settings.m_shouldFormTrackRelationships) - continue; + if (0 == m_settings.m_shouldFormTrackRelationships) continue; // Make track sibling relationships for (unsigned int jTrack = iTrack + 1; jTrack < nTracks; ++jTrack) { - PANDORA_RETURN_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, PandoraApi::SetTrackSiblingRelationship(*m_pPandora, - pTrack, trackVec[jTrack])); + PANDORA_RETURN_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, PandoraApi::SetTrackSiblingRelationship(*m_pPandora, GetTrackAddress(collectionMaps, pTrack), GetTrackAddress(collectionMaps, pReconstructedParticle.getTracks(jTrack) ) ) ); } } } - catch (EVENT::Exception &exception) + catch (...) { - streamlog_out(WARNING) << "Failed to extract v0 vertex: " << exception.what() << std::endl; + std::cout<< "Failed to extract v0 vertex" << std::endl; } } } - catch (EVENT::Exception &exception) + catch (...) { - streamlog_out(DEBUG5) << "Failed to extract v0 vertex collection: " << *iter << ", " << exception.what() << std::endl; + std::cout << "Failed to extract v0 vertex collection: " << *iter << std::endl; } } return pandora::STATUS_CODE_SUCCESS; } -*/ //------------------------------------------------------------------------------------------------------------------------------------------ bool TrackCreator::IsConflictingRelationship(const edm4hep::ConstReconstructedParticle &Particle) const { - //for (unsigned int iTrack = 0, nTracks = trackVec.size(); iTrack < nTracks; ++iTrack) - //for (unsigned int iTrack = 0, nTracks = trackVec.tracks_size(); iTrack < nTracks; ++iTrack) - std::cout<<"Particle en="<<Particle.getEnergy()<<std::endl; for (unsigned int iTrack = 0, nTracks = Particle.tracks_size(); iTrack < nTracks; ++iTrack) { - //EVENT::Track *pTrack = trackVec[iTrack]; edm4hep::ConstTrack pTrack = Particle.getTracks(iTrack) ; unsigned int pTrack_id = pTrack.id() ; - //if (this->IsDaughter(pTrack) || this->IsParent(pTrack) || this->IsV0(pTrack)) if (this->IsDaughter(pTrack_id) || this->IsParent(pTrack_id) || this->IsV0(pTrack_id)) return true; } @@ -456,6 +410,20 @@ bool TrackCreator::IsConflictingRelationship(const edm4hep::ConstReconstructedPa return false; } +const edm4hep::Track* TrackCreator::GetTrackAddress(const CollectionMaps& collectionMaps, const edm4hep::ConstTrack& pTrack ) +{ + for (StringVector::const_iterator iter = m_settings.m_trackCollections.begin(), iterEnd = m_settings.m_trackCollections.end(); iter != iterEnd; ++iter) + { + if(collectionMaps.collectionMap_Track.find(*iter) == collectionMaps.collectionMap_Track.end()) { std::cout<<"not find "<<(*iter)<<std::endl; continue;} + const std::vector<edm4hep::Track>& pTrackCollection = (collectionMaps.collectionMap_Track.find(*iter))->second; + for (int i = 0, iMax = pTrackCollection.size(); i < iMax; ++i) + { + const edm4hep::Track& pTrack0 = pTrackCollection.at(i); + if (pTrack.id() == pTrack0.id()) return (const edm4hep::Track*) (&pTrack0); + } + } + return NULL; +} //------------------------------------------------------------------------------------------------------------------------------------------ pandora::StatusCode TrackCreator::CreateTracks(const CollectionMaps& collectionMaps) @@ -480,8 +448,6 @@ pandora::StatusCode TrackCreator::CreateTracks(const CollectionMaps& collectionM if (NULL == pTrack) throw ("Collection type mismatch"); int minTrackHits = m_settings.m_minTrackHits; - //const float tanLambda(std::fabs(pTrack->getTanLambda())); - //std::cout<<"track states size="<<pTrack->trackStates_size()<<std::endl; const float tanLambda(std::fabs(pTrack->getTrackStates(0).tanLambda)); if (tanLambda > m_tanLambdaFtd) @@ -507,9 +473,7 @@ pandora::StatusCode TrackCreator::CreateTracks(const CollectionMaps& collectionM // Proceed to create the pandora track PandoraApi::Track::Parameters trackParameters; - //trackParameters.m_d0 = pTrack->getD0(); trackParameters.m_d0 = pTrack->getTrackStates(0).D0; - //trackParameters.m_z0 = pTrack->getZ0(); trackParameters.m_z0 = pTrack->getTrackStates(0).Z0; trackParameters.m_pParentAddress = pTrack; // By default, assume tracks are charged pions @@ -559,31 +523,21 @@ pandora::StatusCode TrackCreator::CreateTracks(const CollectionMaps& collectionM void TrackCreator::GetTrackStates(const edm4hep::Track *const pTrack, PandoraApi::Track::Parameters &trackParameters) const { - //const TrackState *pTrackState = pTrack->getTrackState(TrackState::AtIP); edm4hep::TrackState pTrackState = pTrack->getTrackStates(1); // ref /cvmfs/cepcsw.ihep.ac.cn/prototype/LCIO/include/EVENT/TrackState.h - //if (!pTrackState) throw pandora::StatusCodeException(pandora::STATUS_CODE_NOT_INITIALIZED); - //const double pt(m_bField * 2.99792e-4 / std::fabs(pTrackState->getOmega())); const double pt(m_bField * 2.99792e-4 / std::fabs(pTrackState.omega)); - //trackParameters.m_momentumAtDca = pandora::CartesianVector(std::cos(pTrackState->getPhi()), std::sin(pTrackState->getPhi()), pTrackState->getTanLambda()) * pt; trackParameters.m_momentumAtDca = pandora::CartesianVector(std::cos(pTrackState.phi), std::sin(pTrackState.phi), pTrackState.tanLambda) * pt; - //this->CopyTrackState(pTrack->getTrackState(TrackState::AtFirstHit), trackParameters.m_trackStateAtStart); this->CopyTrackState(pTrack->getTrackStates(2), trackParameters.m_trackStateAtStart); - //fg: curling TPC tracks have pointers to track segments stored -> need to get track states from last segment! - //const EVENT::Track *pEndTrack = (pTrack->getTracks().empty() ? pTrack : pTrack->getTracks().back()); auto pEndTrack = (pTrack->tracks_size() ==0 ) ? *pTrack : pTrack->getTracks(pTrack->tracks_size()-1); - //std::cout<<"GetTrackStates 1.6"<<", end track trackStates_size()="<<pEndTrack.trackStates_size()<<std::endl; - //this->CopyTrackState(pEndTrack->getTrackState(TrackState::AtLastHit), trackParameters.m_trackStateAtEnd); - //this->CopyTrackState(pEndTrack->getTrackState(TrackState::AtCalorimeter), trackParameters.m_trackStateAtCalorimeter); this->CopyTrackState(pEndTrack.getTrackStates(3), trackParameters.m_trackStateAtEnd); - //this->CopyTrackState(pEndTrack.getTrackStates(4), trackParameters.m_trackStateAtCalorimeter); //FIXME ? LCIO input only has 4 states, so 4 can't be used. - this->CopyTrackState(pEndTrack.getTrackStates(3), trackParameters.m_trackStateAtCalorimeter); + if( pEndTrack.trackStates_size()<5) this->CopyTrackState(pEndTrack.getTrackStates(3), trackParameters.m_trackStateAtCalorimeter); + else this->CopyTrackState(pEndTrack.getTrackStates(4), trackParameters.m_trackStateAtCalorimeter); trackParameters.m_isProjectedToEndCap = ((std::fabs(trackParameters.m_trackStateAtCalorimeter.Get().GetPosition().GetZ()) < m_eCalEndCapInnerZ) ? false : true); @@ -599,7 +553,6 @@ void TrackCreator::GetTrackStates(const edm4hep::Track *const pTrack, PandoraApi float TrackCreator::CalculateTrackTimeAtCalorimeter(const edm4hep::Track *const pTrack) const { - //const pandora::Helix helix(pTrack->getPhi(), pTrack->getD0(), pTrack->getZ0(), pTrack->getOmega(), pTrack->getTanLambda(), m_bField); const pandora::Helix helix(pTrack->getTrackStates(0).phi, pTrack->getTrackStates(0).D0, pTrack->getTrackStates(0).Z0, pTrack->getTrackStates(0).omega, pTrack->getTrackStates(0).tanLambda, m_bField); const pandora::CartesianVector &referencePoint(helix.GetReferencePoint()); @@ -655,15 +608,10 @@ float TrackCreator::CalculateTrackTimeAtCalorimeter(const edm4hep::Track *const void TrackCreator::CopyTrackState(const edm4hep::TrackState & pTrackState, pandora::InputTrackState &inputTrackState) const { - //if (!pTrackState) throw pandora::StatusCodeException(pandora::STATUS_CODE_NOT_INITIALIZED); - //const double pt(m_bField * 2.99792e-4 / std::fabs(pTrackState->getOmega())); const double pt(m_bField * 2.99792e-4 / std::fabs(pTrackState.omega)); - //const double px(pt * std::cos(pTrackState->getPhi())); const double px(pt * std::cos(pTrackState.phi)); - //const double py(pt * std::sin(pTrackState->getPhi())); const double py(pt * std::sin(pTrackState.phi)); - //const double pz(pt * pTrackState->getTanLambda()); const double pz(pt * pTrackState.tanLambda); const double xs(pTrackState.referencePoint[0]); @@ -713,7 +661,6 @@ void TrackCreator::TrackReachesECAL(const edm4hep::Track *const pTrack, PandoraA (std::fabs(z) - m_settings.m_reachesECalFtdZMaxDistance < m_ftdZPositions[j]) && (std::fabs(z) + m_settings.m_reachesECalFtdZMaxDistance > m_ftdZPositions[j])) { - //if (static_cast<int>(j) > maxOccupiedFtdLayer) maxOccupiedFtdLayer = static_cast<int>(j); if ( j > maxOccupiedFtdLayer) maxOccupiedFtdLayer = j; break; } @@ -765,19 +712,14 @@ void TrackCreator::DefineTrackPfoUsage(const edm4hep::Track *const pTrack, Pando bool canFormPfo(false); bool canFormClusterlessPfo(false); - //if (trackParameters.m_reachesCalorimeter.Get() && !this->IsParent(pTrack)) if (trackParameters.m_reachesCalorimeter.Get() && !this->IsParent(pTrack->id())) { - //const float d0(std::fabs(pTrack->getD0())), z0(std::fabs(pTrack->getZ0())); const float d0(std::fabs(pTrack->getTrackStates(0).D0)), z0(std::fabs(pTrack->getTrackStates(0).Z0)); - //EVENT::TrackerHitVec trackerHitvec(pTrack->getTrackerHits()); float rInner(std::numeric_limits<float>::max()), zMin(std::numeric_limits<float>::max()); - //for (EVENT::TrackerHitVec::const_iterator iter = trackerHitvec.begin(), iterEnd = trackerHitvec.end(); iter != iterEnd; ++iter) for (std::vector<edm4hep::ConstTrackerHit>::const_iterator iter = pTrack->trackerHits_begin(), iterEnd = pTrack->trackerHits_end(); iter != iterEnd; ++iter) { - //const double *pPosition((*iter)->getPosition()); const edm4hep::Vector3d pPosition = (*iter).getPosition(); const float x(pPosition[0]), y(pPosition[1]), absoluteZ(std::fabs(pPosition[2])); const float r(std::sqrt(x * x + y * y)); @@ -855,7 +797,6 @@ bool TrackCreator::PassesQualityCuts(const edm4hep::Track *const pTrack, const P if (trackParameters.m_trackStateAtCalorimeter.Get().GetPosition().GetMagnitude() < m_settings.m_minTrackECalDistanceFromIp) return false; - //if (std::fabs(pTrack->getOmega()) < std::numeric_limits<float>::epsilon()) if (std::fabs(pTrack->getTrackStates(0).omega) < std::numeric_limits<float>::epsilon()) { std::cout<<"ERROR Track has Omega = 0 " << std::endl; @@ -864,7 +805,6 @@ bool TrackCreator::PassesQualityCuts(const edm4hep::Track *const pTrack, const P // Check momentum uncertainty is reasonable to use track const pandora::CartesianVector &momentumAtDca(trackParameters.m_momentumAtDca.Get()); - //const float sigmaPOverP(std::sqrt(pTrack->getCovMatrix()[5]) / std::fabs(pTrack->getOmega())); const float sigmaPOverP(std::sqrt(pTrack->getTrackStates(0).covMatrix[5]) / std::fabs(pTrack->getTrackStates(0).omega)); if (sigmaPOverP > m_settings.m_maxTrackSigmaPOverP) @@ -930,12 +870,6 @@ bool TrackCreator::PassesQualityCuts(const edm4hep::Track *const pTrack, const P int TrackCreator::GetNTpcHits(const edm4hep::Track *const pTrack) const { - // ATTN - //fg: hit numbers are now given in different order wrt LOI: - // trk->subdetectorHitNumbers()[ 2 * ILDDetID::TPC - 1 ] = hitsInFit ; - // trk->subdetectorHitNumbers()[ 2 * ILDDetID::TPC - 2 ] = hitCount ; - // ---- use hitsInFit : - //return pTrack->getSubdetectorHitNumbers()[ 2 * lcio::ILDDetID::TPC - 1 ]; return pTrack->getSubDetectorHitNumbers(2 * lcio::ILDDetID::TPC - 1);// still use LCIO code now } @@ -943,12 +877,6 @@ int TrackCreator::GetNTpcHits(const edm4hep::Track *const pTrack) const int TrackCreator::GetNFtdHits(const edm4hep::Track *const pTrack) const { - // ATTN - //fg: hit numbers are now given in different order wrt LOI: - // trk->subdetectorHitNumbers()[ 2 * ILDDetID::TPC - 1 ] = hitsInFit ; - // trk->subdetectorHitNumbers()[ 2 * ILDDetID::TPC - 2 ] = hitCount ; - // ---- use hitsInFit : - //return pTrack->getSubdetectorHitNumbers()[ 2 * lcio::ILDDetID::FTD - 1 ]; return pTrack->getSubDetectorHitNumbers( 2 * lcio::ILDDetID::FTD - 1 ); } diff --git a/Reconstruction/PFA/Pandora/PandoraSettingsDefault.xml b/Reconstruction/PFA/Pandora/PandoraSettingsDefault.xml index be26c3643a7017db4f3158da6d6ef338fd30a92e..81a41e39c6bf686f35ea5664624c8e2ffde2efd9 100644 --- a/Reconstruction/PFA/Pandora/PandoraSettingsDefault.xml +++ b/Reconstruction/PFA/Pandora/PandoraSettingsDefault.xml @@ -69,7 +69,8 @@ <ClusterListName>PhotonClusters</ClusterListName> <ReplaceCurrentClusterList>false</ReplaceCurrentClusterList> <ShouldMakePdfHistograms>false</ShouldMakePdfHistograms> - <HistogramFile>/junofs/users/wxfang/MyGit/tmp/fork_update_pandora/CEPCSW/Reconstruction/PFA/Pandora/PandoraLikelihoodData9EBin.xml</HistogramFile> + <HistogramFile>./Reconstruction/PFA/Pandora/PandoraLikelihoodData9EBin.xml</HistogramFile> + <!--HistogramFile>/junofs/users/wxfang/MyGit/tmp/fork_update_pandora/CEPCSW/Reconstruction/PFA/Pandora/PandoraLikelihoodData9EBin.xml</HistogramFile--> </algorithm> <!-- Clustering parent algorithm runs a daughter clustering algorithm --> diff --git a/Service/DedxSvc/CMakeLists.txt b/Service/DedxSvc/CMakeLists.txt deleted file mode 100644 index 0fce402f937a8b428db7e04bd3b82b83b963ac4b..0000000000000000000000000000000000000000 --- a/Service/DedxSvc/CMakeLists.txt +++ /dev/null @@ -1,13 +0,0 @@ -gaudi_subdir(DedxSvc v0r0) - -set(DedxSvc_srcs - src/*.cpp -) -find_package(Geant4 REQUIRED ui_all vis_all) -include(${Geant4_USE_FILE}) -gaudi_install_headers(DedxSvc) - -gaudi_add_module(DedxSvc ${DedxSvc_srcs} - INCLUDE_DIRS GaudiKernel - LINK_LIBRARIES GaudiKernel -) diff --git a/Service/DedxSvc/DedxSvc/IDedxSvc.h b/Service/DedxSvc/DedxSvc/IDedxSvc.h deleted file mode 100644 index 14b76620c9cd503bd72d724eb46acb8d9b9fc8bd..0000000000000000000000000000000000000000 --- a/Service/DedxSvc/DedxSvc/IDedxSvc.h +++ /dev/null @@ -1,18 +0,0 @@ -#ifndef I_Dedx_SVC_H -#define I_Dedx_SVC_H - -#include "GaudiKernel/IService.h" - -#include "G4Step.hh" - -class IDedxSvc: virtual public IService { -public: - DeclareInterfaceID(IDedxSvc, 0, 1); // major/minor version - - virtual ~IDedxSvc() = default; - virtual float pred(const G4Step* aStep)=0 ; - -}; - - -#endif diff --git a/Service/DedxSvc/src/DedxSvc.cpp b/Service/DedxSvc/src/DedxSvc.cpp deleted file mode 100644 index 9a62508dd0497e95fc7cfa5dc4cd3d4504354b8b..0000000000000000000000000000000000000000 --- a/Service/DedxSvc/src/DedxSvc.cpp +++ /dev/null @@ -1,45 +0,0 @@ -#include "DedxSvc.h" - -//https://folk.uib.no/ruv004/ -DECLARE_COMPONENT(DedxSvc) - -DedxSvc::DedxSvc(const std::string& name, ISvcLocator* svc) - : base_class(name, svc) -{ -} - -DedxSvc::~DedxSvc() -{ -} - -float DedxSvc::pred(const G4Step* aStep) -{ - G4Track* gTrack = aStep->GetTrack() ; - G4int z = gTrack->GetDefinition()->GetPDGCharge(); - if (z == 0) return 0; - G4double M = gTrack->GetDefinition()->GetPDGMass();//MeV - M = pow(10,6)*M; //to eV - G4double gammabeta=aStep->GetPreStepPoint()->GetBeta() * aStep->GetPreStepPoint()->GetGamma(); - if(gammabeta<0.01)return 0;//too low momentum - float beta = gammabeta/sqrt(1.0+pow(gammabeta,2)); - float gamma = gammabeta/beta; - float Tmax = 2*m_me*pow(gammabeta,2)/(1+(2*gamma*m_me/M)+pow(m_me/M,2)); - float dedx = m_K*pow(z,2)*m_material_Z*(0.5*log(2*m_me*pow(gammabeta,2)*Tmax/pow(m_I,2))-pow(beta,2))/(m_material_A*pow(beta,2)); - dedx = dedx*m_scale;// the material density can be absorbed in scale - dedx = dedx*(1+((*m_distribution)(m_generator))); - return dedx*m_material_density; // MeV / cm -} - -StatusCode DedxSvc::initialize() -{ - m_distribution = new std::normal_distribution<double>(0, m_resolution); - m_me = 0.511*pow(10,6);//0.511 MeV to eV - m_K = 0.307075;//const - m_I = m_material_Z*10; // Approximate - return StatusCode::SUCCESS; -} - -StatusCode DedxSvc::finalize() -{ - return StatusCode::SUCCESS; -} diff --git a/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.cpp b/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.cpp index 9ef16425454ad9a520659b07536d7ebf344018cc..954664bf21a5897f86285f1380299bf16068df1a 100644 --- a/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.cpp +++ b/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.cpp @@ -173,7 +173,7 @@ Edm4hepWriterAnaElemTool::EndOfEventAction(const G4Event* anEvent) { } dd4hep::sim::Geant4TrackerHit* trk_hit = dynamic_cast<dd4hep::sim::Geant4TrackerHit*>(h); - if (trk_hit) { + if (trk_hit && tracker_col_ptr) { ++n_trk_hit; // auto edm_trk_hit = trackercols->create(); auto edm_trk_hit = tracker_col_ptr->create(); @@ -192,10 +192,11 @@ Edm4hepWriterAnaElemTool::EndOfEventAction(const G4Event* anEvent) { trk_hit->momentum.y()/CLHEP::GeV, trk_hit->momentum.z()/CLHEP::GeV}; edm_trk_hit.setMomentum(edm4hep::Vector3f(mom)); + } dd4hep::sim::Geant4CalorimeterHit* cal_hit = dynamic_cast<dd4hep::sim::Geant4CalorimeterHit*>(h); - if (cal_hit) { + if (cal_hit && calo_col_ptr) { ++n_cal_hit; auto edm_calo_hit = calo_col_ptr->create(); edm_calo_hit.setCellID(cal_hit->cellID); diff --git a/Simulation/DetSimDedx/CMakeLists.txt b/Simulation/DetSimDedx/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..824c35b37314dc7b007227c1082b032ea631336b --- /dev/null +++ b/Simulation/DetSimDedx/CMakeLists.txt @@ -0,0 +1,24 @@ +gaudi_subdir(DetSimDedx v0r0) + +gaudi_depends_on_subdirs( + FWCore + Simulation/DetSimInterface +) + +find_package(Geant4 REQUIRED ui_all vis_all) +include(${Geant4_USE_FILE}) +find_package(DD4hep COMPONENTS DDG4 REQUIRED) + +set(DetSimDedx_srcs + src/DummyDedxSimTool.cpp + src/BetheBlochEquationDedxSimTool.cpp +) + +gaudi_add_module(DetSimDedx ${DetSimDedx_srcs} + INCLUDE_DIRS + LINK_LIBRARIES + DD4hep + ${DD4hep_COMPONENT_LIBRARIES} + GaudiKernel +) + diff --git a/Simulation/DetSimDedx/src/BetheBlochEquationDedxSimTool.cpp b/Simulation/DetSimDedx/src/BetheBlochEquationDedxSimTool.cpp new file mode 100644 index 0000000000000000000000000000000000000000..6aecca2bb3bf4651e1c5c06f70fc7bfd2d8ea6bc --- /dev/null +++ b/Simulation/DetSimDedx/src/BetheBlochEquationDedxSimTool.cpp @@ -0,0 +1,52 @@ +#include "BetheBlochEquationDedxSimTool.h" +#include "G4Step.hh" +#include "G4SystemOfUnits.hh" + +// https://folk.uib.no/ruv004/ +DECLARE_COMPONENT(BetheBlochEquationDedxSimTool) + + +double BetheBlochEquationDedxSimTool::dedx(const G4Step* aStep) +{ + G4Track* gTrack = aStep->GetTrack() ; + G4int z = gTrack->GetDefinition()->GetPDGCharge(); + if (z == 0) return 0; + + G4Material* material = gTrack->GetMaterial(); + G4double material_density = material->GetDensity() / (CLHEP::g/CLHEP::cm3); // conert from G4 unit. + G4double material_Z = material->GetZ(); + G4double material_A = material->GetA(); + + m_I = material_Z*10; // Approximate + + G4double M = gTrack->GetDefinition()->GetPDGMass();//MeV + M = pow(10,6)*M; //to eV + G4double gammabeta=aStep->GetPreStepPoint()->GetBeta() * aStep->GetPreStepPoint()->GetGamma(); + if(gammabeta<0.01)return 0;//too low momentum + float beta = gammabeta/sqrt(1.0+pow(gammabeta,2)); + float gamma = gammabeta/beta; + float Tmax = 2*m_me*pow(gammabeta,2)/(1+(2*gamma*m_me/M)+pow(m_me/M,2)); + float dedx = m_K*pow(z,2)*material_Z*(0.5*log(2*m_me*pow(gammabeta,2)*Tmax/pow(m_I,2))-pow(beta,2))/(material_A*pow(beta,2)); + dedx = dedx*m_scale;// the material density can be absorbed in scale + dedx = dedx*(1+((*m_distribution)(m_generator))); + return dedx*material_density; // MeV / cm +} + +StatusCode BetheBlochEquationDedxSimTool::initialize() +{ + m_distribution = new std::normal_distribution<double>(0, m_resolution); + m_me = 0.511*pow(10,6);//0.511 MeV to eV + m_K = 0.307075;//const + + + info() << "Initialize BetheBlochEquationDedxSimTool with following parameters" << endmsg; + info() << "-> m_me: " << m_me << endmsg; + info() << "-> m_K: " << m_K << endmsg; + + return StatusCode::SUCCESS; +} + +StatusCode BetheBlochEquationDedxSimTool::finalize() +{ + return StatusCode::SUCCESS; +} diff --git a/Service/DedxSvc/src/DedxSvc.h b/Simulation/DetSimDedx/src/BetheBlochEquationDedxSimTool.h similarity index 71% rename from Service/DedxSvc/src/DedxSvc.h rename to Simulation/DetSimDedx/src/BetheBlochEquationDedxSimTool.h index 6f1f85eb53dbabf872cbbe56aee9a44e4656a639..6e80ad7972d1fcac831f709cde20f77f8bf97f5a 100644 --- a/Service/DedxSvc/src/DedxSvc.h +++ b/Simulation/DetSimDedx/src/BetheBlochEquationDedxSimTool.h @@ -1,21 +1,17 @@ -#ifndef Dedx_SVC_H -#define Dedx_SVC_H +#ifndef BetheBlochEquationDedxSimTool_h +#define BetheBlochEquationDedxSimTool_h -#include "DedxSvc/IDedxSvc.h" -#include <GaudiKernel/Service.h> -#include "G4Step.hh" +#include "DetSimInterface/IDedxSimTool.h" +#include <GaudiKernel/AlgTool.h> #include <random> -class DedxSvc : public extends<Service, IDedxSvc> -{ +class BetheBlochEquationDedxSimTool: public extends<AlgTool, IDedxSimTool> { public: - DedxSvc(const std::string& name, ISvcLocator* svc); - ~DedxSvc(); - + using extends::extends; StatusCode initialize() override; StatusCode finalize() override; - float pred(const G4Step* aStep) override; + double dedx(const G4Step* aStep) override; private: @@ -27,6 +23,7 @@ class DedxSvc : public extends<Service, IDedxSvc> float m_me;// Here me is the electron rest mass float m_K; // K was set as a constant. float m_I; // Mean excitation energy + std::default_random_engine m_generator; std::normal_distribution<double>* m_distribution; }; diff --git a/Simulation/DetSimDedx/src/DummyDedxSimTool.cpp b/Simulation/DetSimDedx/src/DummyDedxSimTool.cpp new file mode 100644 index 0000000000000000000000000000000000000000..291b8c205eafd630a9e45f9be64a99cbd0827c87 --- /dev/null +++ b/Simulation/DetSimDedx/src/DummyDedxSimTool.cpp @@ -0,0 +1,24 @@ +#include "DummyDedxSimTool.h" + +#include "G4Step.hh" + +DECLARE_COMPONENT(DummyDedxSimTool); + +StatusCode DummyDedxSimTool::initialize() { + StatusCode sc; + + return sc; +} + +StatusCode DummyDedxSimTool::finalize() { + StatusCode sc; + + return sc; +} + +double DummyDedxSimTool::dedx(const G4Step* aStep) { + double result = aStep->GetTotalEnergyDeposit(); + + + return result; +} diff --git a/Simulation/DetSimDedx/src/DummyDedxSimTool.h b/Simulation/DetSimDedx/src/DummyDedxSimTool.h new file mode 100644 index 0000000000000000000000000000000000000000..ce7f3e3e4a49430a2f9736e03c4538c2cb87878f --- /dev/null +++ b/Simulation/DetSimDedx/src/DummyDedxSimTool.h @@ -0,0 +1,21 @@ +#ifndef DummyDedxSimTool_h +#define DummyDedxSimTool_h + +#include "GaudiKernel/AlgTool.h" +#include "DetSimInterface/IDedxSimTool.h" + +class DummyDedxSimTool: public extends<AlgTool, IDedxSimTool> { + +public: + using extends::extends; + + /// Overriding initialize and finalize + StatusCode initialize() override; + StatusCode finalize() override; + + /// Overriding dedx tool + double dedx(const G4Step* aStep) override; + +}; + +#endif diff --git a/Simulation/DetSimGeom/src/AnExampleDetElemTool.cpp b/Simulation/DetSimGeom/src/AnExampleDetElemTool.cpp index 87ad61b35a7fb05c9bd18999bcc3ea5e40d66c9f..e78170bc5eaf9534a29126630458285166421772 100644 --- a/Simulation/DetSimGeom/src/AnExampleDetElemTool.cpp +++ b/Simulation/DetSimGeom/src/AnExampleDetElemTool.cpp @@ -110,6 +110,17 @@ AnExampleDetElemTool::ConstructSDandField() { } } else if (typ=="tracker") { + // if drift chamber + if (nam == "DriftChamber") { + m_driftchamber_sdtool = ToolHandle<ISensDetTool>("DriftChamberSensDetTool"); + if (m_driftchamber_sdtool) { + info() << "Find the DriftChamberSensDetTool" << endmsg; + g4sd = m_driftchamber_sdtool->createSD(nam); + } else { + warning() << "DriftChamberSensDetTool is not found. " << endmsg; + } + } + } } @@ -186,6 +197,9 @@ AnExampleDetElemTool::initialize() { return StatusCode::FAILURE; } + m_calo_sdtool = ToolHandle<ISensDetTool>("CalorimeterSensDetTool"); + m_driftchamber_sdtool = ToolHandle<ISensDetTool>("DriftChamberSensDetTool"); + return sc; } diff --git a/Simulation/DetSimGeom/src/AnExampleDetElemTool.h b/Simulation/DetSimGeom/src/AnExampleDetElemTool.h index b0b22861b892c0817876e9c325853e843acf7dc4..ba1990715d10ca9f07eb5a2791314f8ce0d090ee 100644 --- a/Simulation/DetSimGeom/src/AnExampleDetElemTool.h +++ b/Simulation/DetSimGeom/src/AnExampleDetElemTool.h @@ -33,6 +33,7 @@ private: SmartIF<IGeoSvc> m_geosvc; ToolHandle<ISensDetTool> m_calo_sdtool; + ToolHandle<ISensDetTool> m_driftchamber_sdtool; }; #endif diff --git a/Simulation/DetSimInterface/DetSimInterface/IDedxSimTool.h b/Simulation/DetSimInterface/DetSimInterface/IDedxSimTool.h new file mode 100644 index 0000000000000000000000000000000000000000..f47ceb0211c3a61553adca623e443f25e46c5038 --- /dev/null +++ b/Simulation/DetSimInterface/DetSimInterface/IDedxSimTool.h @@ -0,0 +1,29 @@ +#ifndef IDedxSimTool_h +#define IDedxSimTool_h + +/* + * Description: + * IDedxSimTool is used to give a dE/dx value during simulation. + * + * The interface: + * * dedx: predict the dE/dx according to Geant4 Step + * + * Author: Tao Lin <lintao@ihep.ac.cn> + */ + + +#include "GaudiKernel/IAlgTool.h" + +class G4Step; + +class IDedxSimTool: virtual public IAlgTool { +public: + + DeclareInterfaceID(IDedxSimTool, 0, 1); + virtual ~IDedxSimTool() {} + + virtual double dedx(const G4Step* aStep) = 0; + +}; + +#endif diff --git a/Simulation/DetSimSD/CMakeLists.txt b/Simulation/DetSimSD/CMakeLists.txt index 10b1a15a3f0d9ed6cf91a635ea6d53eca49daba3..d917ecef05a6bb5253e4a1e1529472ab939e184b 100644 --- a/Simulation/DetSimSD/CMakeLists.txt +++ b/Simulation/DetSimSD/CMakeLists.txt @@ -16,6 +16,9 @@ set(DetSimSD_srcs src/DDG4SensitiveDetector.cpp src/CaloSensitiveDetector.cpp + + src/DriftChamberSensDetTool.cpp + src/DriftChamberSensitiveDetector.cpp ) gaudi_add_module(DetSimSD ${DetSimSD_srcs} diff --git a/Simulation/DetSimSD/src/DriftChamberSensDetTool.cpp b/Simulation/DetSimSD/src/DriftChamberSensDetTool.cpp new file mode 100644 index 0000000000000000000000000000000000000000..bf0532165fcf820a10563af25f9ff172fbdacc8e --- /dev/null +++ b/Simulation/DetSimSD/src/DriftChamberSensDetTool.cpp @@ -0,0 +1,48 @@ +#include "DriftChamberSensDetTool.h" + +#include "G4VSensitiveDetector.hh" + +#include "DD4hep/Detector.h" + +#include "DriftChamberSensitiveDetector.h" + +DECLARE_COMPONENT(DriftChamberSensDetTool); + +StatusCode DriftChamberSensDetTool::initialize() { + StatusCode sc; + + m_geosvc = service<IGeoSvc>("GeoSvc"); + if (!m_geosvc) { + error() << "Failed to find GeoSvc." << endmsg; + return StatusCode::FAILURE; + } + + m_dedx_simtool = ToolHandle<IDedxSimTool>(m_dedx_sim_option.value()); + if (!m_dedx_simtool) { + error() << "Failed to find dedx simtoo." << endmsg; + return StatusCode::FAILURE; + } + + return sc; +} + +StatusCode DriftChamberSensDetTool::finalize() { + StatusCode sc; + + return sc; +} + +G4VSensitiveDetector* +DriftChamberSensDetTool::createSD(const std::string& name) { + dd4hep::Detector* dd4hep_geo = m_geosvc->lcdd(); + + G4VSensitiveDetector* sd = nullptr; + + if (name == "DriftChamber") { + DriftChamberSensitiveDetector* dcsd = new DriftChamberSensitiveDetector(name, *dd4hep_geo); + dcsd->setDedxSimTool(m_dedx_simtool); + } + + + return sd; +} diff --git a/Simulation/DetSimSD/src/DriftChamberSensDetTool.h b/Simulation/DetSimSD/src/DriftChamberSensDetTool.h new file mode 100644 index 0000000000000000000000000000000000000000..ae727ed09fbcb58c5ea918cf3e86813bd2625a2f --- /dev/null +++ b/Simulation/DetSimSD/src/DriftChamberSensDetTool.h @@ -0,0 +1,41 @@ +#ifndef DriftChamberSensDetTool_h +#define DriftChamberSensDetTool_h + +/* + * DriftChamberSensDetTool is used to create Drift Chamber SD. + * + * It will use DedxSimTool to give the dE/dx value. + * + * -- 17 Sept 2020, Tao Lin <lintao@ihep.ac.cn> + */ + +#include "GaudiKernel/AlgTool.h" +#include "GaudiKernel/ToolHandle.h" +#include "DetSimInterface/ISensDetTool.h" +#include "DetSimInterface/IDedxSimTool.h" +#include "DetInterface/IGeoSvc.h" + + +class DriftChamberSensDetTool: public extends<AlgTool, ISensDetTool> { + +public: + + using extends::extends; + + /// Overriding initialize and finalize + StatusCode initialize() override; + StatusCode finalize() override; + + /// Override ISensDetTool + virtual G4VSensitiveDetector* createSD(const std::string& name) override; + +private: + + // in order to initialize SD, we need to get the lcdd() + SmartIF<IGeoSvc> m_geosvc; + ToolHandle<IDedxSimTool> m_dedx_simtool; + Gaudi::Property<std::string> m_dedx_sim_option{this, "DedxSimTool"}; + +}; + +#endif diff --git a/Simulation/DetSimSD/src/DriftChamberSensitiveDetector.cpp b/Simulation/DetSimSD/src/DriftChamberSensitiveDetector.cpp new file mode 100644 index 0000000000000000000000000000000000000000..1ba9f7bff59a19d93519bf36368daee33035a7e9 --- /dev/null +++ b/Simulation/DetSimSD/src/DriftChamberSensitiveDetector.cpp @@ -0,0 +1,73 @@ +#include "DriftChamberSensitiveDetector.h" + +#include "G4SDManager.hh" + +DriftChamberSensitiveDetector::DriftChamberSensitiveDetector(const std::string& name, + dd4hep::Detector& description) + : DDG4SensitiveDetector(name, description), + m_hc(nullptr) { + + const std::string& coll_name = m_sensitive.hitsCollection(); + + collectionName.insert(coll_name); +} + +bool DriftChamberSensitiveDetector::setDedxSimTool(ToolHandle<IDedxSimTool> simtool) { + m_dedx_simtool = simtool; + + return true; +} + +void +DriftChamberSensitiveDetector::Initialize(G4HCofThisEvent* HCE) { + + const std::string& coll_name = collectionName[0]; + m_hc = new HitCollection(GetName(), coll_name); + + int HCID = -1; + if(HCID<0) HCID = G4SDManager::GetSDMpointer()->GetCollectionID(m_hc); + HCE->AddHitsCollection( HCID, m_hc ); + +} + +G4bool +DriftChamberSensitiveDetector::ProcessHits(G4Step* step, G4TouchableHistory*) { + // Refer to: DDG4/legacy/Geant4TrackerSD.cpp (note: there's bug in momentum calculation) + // DDCore/include/DD4hep/Objects.h (mean_direction and mean_length) + + dd4hep::sim::Geant4StepHandler h(step); + + dd4hep::Position prePos = h.prePos(); + dd4hep::Position postPos = h.postPos(); + dd4hep::Position direction = postPos - prePos; + dd4hep::Position position = mean_direction(prePos,postPos); // (pre+post)/2 + double hit_len = direction.R(); + + HitContribution contrib = dd4hep::sim::Geant4Hit::extractContribution(step); + // Now, invokes the dE/dx simulator + double dedx = 0.0; + + double de = hit_len * dedx; + // contrib.deposit = de; // if need the de from dedx simulator + + // create a new hit + TrackerHit* hit = new TrackerHit( + h.track->GetTrackID(), + h.track->GetDefinition()->GetPDGEncoding(), + de, // not the Geant4's deposit energy. from dE/dx simulator + h.track->GetGlobalTime() + ); + hit->cellID = getCellID(step); + hit->energyDeposit = de; // FIXME: also use the dedx + hit->position = position; + hit->momentum = (h.preMom() + h.postMom() )/2; + hit->length = hit_len; + m_hc->insert(hit); + + return true; +} + +void +DriftChamberSensitiveDetector::EndOfEvent(G4HCofThisEvent* HCE) { + +} diff --git a/Simulation/DetSimSD/src/DriftChamberSensitiveDetector.h b/Simulation/DetSimSD/src/DriftChamberSensitiveDetector.h new file mode 100644 index 0000000000000000000000000000000000000000..638e146bcad688dc5a0fdefb4f1178394c6c09ad --- /dev/null +++ b/Simulation/DetSimSD/src/DriftChamberSensitiveDetector.h @@ -0,0 +1,40 @@ +#ifndef DriftChamberSensitiveDetector_h +#define DriftChamberSensitiveDetector_h + +/* + * DriftChamberSensitiveDetector is used in Drift Chamber with dE/dx simulator. + * + * 19 Sept. 2020, Tao Lin <lintao@ihep.ac.cn> + */ + +#include "DetSimSD/DDG4SensitiveDetector.h" +#include "DetSimInterface/IDedxSimTool.h" +#include "GaudiKernel/ToolHandle.h" + +class DriftChamberSensitiveDetector: public DDG4SensitiveDetector { +public: + typedef dd4hep::sim::Geant4TrackerHit TrackerHit; + typedef G4THitsCollection<TrackerHit> TrackerHitCollection; + +public: + DriftChamberSensitiveDetector(const std::string& name, dd4hep::Detector& description); + + bool setDedxSimTool(ToolHandle<IDedxSimTool>); + +public: + // Geant4 interface + + virtual void Initialize(G4HCofThisEvent* HCE); + virtual G4bool ProcessHits(G4Step* step,G4TouchableHistory* history); + virtual void EndOfEvent(G4HCofThisEvent* HCE); + +protected: + + HitCollection* m_hc; + + // this is passed from SensDetTool + ToolHandle<IDedxSimTool> m_dedx_simtool; + +}; + +#endif diff --git a/docs/quickstart.md b/docs/quickstart.md index 5f2de3f2b366bab6e3dabfab255ff65e68089735..7f31084bdcdffd68d76c3763933591d37dafdf6c 100644 --- a/docs/quickstart.md +++ b/docs/quickstart.md @@ -7,12 +7,28 @@ Start the container in lxslc7 (OS: CentOS7): $ /cvmfs/container.ihep.ac.cn/bin/hep_container shell SL6 ``` -## Get code using git +## Manage code using git +Fork the CEPCSW into your own repo. For an example: +* https://github.com/*USERNAME*/CEPCSW + +Get the source code from your own repo: ``` -$ git clone git@github.com:cepc/CEPCSW.git +$ git clone git@github.com:USERNAME/CEPCSW.git $ cd CEPCSW ``` + +Add the upstream repo: +``` +$ git remote add cepc https://github.com/cepc/CEPCSW.git +``` + +Sync and merge the upstream repo: +``` +$ git remote update +$ git merge cepc/master +``` + ## Setup and build ``` diff --git a/docs/simulation-tutorial/index.html b/docs/simulation-tutorial/index.html new file mode 100644 index 0000000000000000000000000000000000000000..6a59e836c53728605595a379e217b7fb7ffc8462 --- /dev/null +++ b/docs/simulation-tutorial/index.html @@ -0,0 +1,415 @@ + +<!doctype html> +<html lang="en"> + <head> + <meta charset="utf-8"> + <meta name="apple-mobile-web-app-capable" content="yes"> + <meta name="apple-mobile-web-app-status-bar-style" content="black"> + <meta name="mobile-web-app-capable" content="yes"> + + + <meta name="description" content=" # Tutorial on CEPCSW simulation Tao Lin IHEP 17 Sept. 2020 [CEPCSW Tutorial, 2020, IHEP](https://"> + + <title>Tutorial on CEPCSW simulation - CodiMD</title> + <link rel="icon" type="image/png" href="https://jupyter.ihep.ac.cn/favicon.png"> + <link rel="apple-touch-icon" href="https://jupyter.ihep.ac.cn/apple-touch-icon.png"> + + + <link rel="stylesheet" href="https://cdnjs.cloudflare.com/ajax/libs/fork-awesome/1.1.7/css/fork-awesome.min.css" integrity="sha256-gsmEoJAws/Kd3CjuOQzLie5Q3yshhvmo7YNtBG7aaEY=" crossorigin="anonymous" /> + <link rel="stylesheet" href="https://cdnjs.cloudflare.com/ajax/libs/ionicons/2.0.1/css/ionicons.min.css" integrity="sha256-3iu9jgsy9TpTwXKb7bNQzqWekRX7pPK+2OLj3R922fo=" crossorigin="anonymous" /> + <link rel="stylesheet" href="https://cdn.jsdelivr.net/npm/reveal.js@3.9.2/css/reveal.min.css" integrity="sha256-h2NhWerL2k7KAzo6YqYMo1T5B6+QT2Bb/CprRV2aW20=" crossorigin="anonymous" /> + <link rel="stylesheet" href="https://cdn.jsdelivr.net/npm/@hackmd/emojify.js@2.1.0/dist/css/basic/emojify.min.css" integrity="sha256-UOrvMOsSDSrW6szVLe8ZDZezBxh5IoIfgTwdNDgTjiU=" crossorigin="anonymous" /> + <link rel="stylesheet" href="https://cdnjs.cloudflare.com/ajax/libs/leaflet/1.6.0/leaflet.css" integrity="sha256-SHMGCYmST46SoyGgo4YR/9AlK1vf3ff84Aq9yK4hdqM=" crossorigin="anonymous" /> + <link href="https://jupyter.ihep.ac.cn/build/font.css" rel="stylesheet"><link href="https://jupyter.ihep.ac.cn/build/slide-styles.css" rel="stylesheet"><link href="https://jupyter.ihep.ac.cn/build/slide.css" rel="stylesheet"> + <!-- HTML5 shim and Respond.js for IE8 support of HTML5 elements and media queries --> +<!-- WARNING: Respond.js doesn't work if you view the page via file:// --> +<!--[if lt IE 9]> + <script src="https://cdnjs.cloudflare.com/ajax/libs/html5shiv/3.7.3/html5shiv.min.js" integrity="sha256-3Jy/GbSLrg0o9y5Z5n1uw0qxZECH7C6OQpVBgNFYa0g=" crossorigin="anonymous"></script> + <script src="https://cdnjs.cloudflare.com/ajax/libs/respond.js/1.4.2/respond.min.js" integrity="sha256-g6iAfvZp+nDQ2TdTR/VVKJf3bGro4ub5fvWSWVRi2NE=" crossorigin="anonymous"></script> + <script src="https://cdnjs.cloudflare.com/ajax/libs/es5-shim/4.5.9/es5-shim.min.js" integrity="sha256-8E4Is26QH0bD52WoQpcB+R/tcWQtpzlCojrybUd7Mxo=" crossorigin="anonymous"></script> +<![endif]--> + + + + <!-- For reveal.js theme --> + + <link rel="stylesheet" href="https://jupyter.ihep.ac.cn/build/reveal.js/css/theme/black.css" id="theme"> + + <!-- For syntax highlighting --> + <link rel="stylesheet" href="https://jupyter.ihep.ac.cn/build/reveal.js/lib/css/zenburn.css"> + <!-- For overwrite reveal.js --> + <link rel="stylesheet" href="https://jupyter.ihep.ac.cn/css/slide.css"> + + <!-- Printing and PDF exports --> + <script nonce="840f3704-7c5d-407c-bee3-d2af62757acb"> + var link = document.createElement( 'link' ); + link.rel = 'stylesheet'; + link.type = 'text/css'; + link.href = 'https://jupyter.ihep.ac.cn/build/reveal.js/' + (window.location.search.match( /print-pdf/gi ) ? 'css/print/pdf.css' : 'css/print/paper.css'); + document.getElementsByTagName( 'head' )[0].appendChild( link ); + </script> + </head> + <body> + <div class="container"> + <div class="reveal"> + <div class="slides"> + +# Tutorial on CEPCSW simulation + +Tao Lin +IHEP +17 Sept. 2020 +[CEPCSW Tutorial, 2020, IHEP](https://indico.ihep.ac.cn/event/12341/other-view?view=standard) + +--- + +## What will learn in this Tutorial? + +* As Users + * Run a simple simulation job in CEPCSW + * Understand and Customize the simulation + * Analyze the simulation output +* As Developers + * Understand the simulation framework + * Learn basics on Geant4 simulation + * Simulate with different detector options +* Note: the emoji :writing_hand: means exercises + +--- + +### CEPCSW Simulation Framework + + + +* The simulation chain is driven by Gaudi. +* Detector description is from DD4hep. +* Event Data Model is in EDM4hep format. +* Detector response is done by Geant4. + +---- + +#### Code is on GitHub +* Detector description: [See Detector](https://github.com/cepc/CEPCSW/tree/master/Detector) +* Event generator interface: [See Generator](https://github.com/cepc/CEPCSW/tree/master/Generator) +* Detector simulation: [See Simulation](https://github.com/cepc/CEPCSW/tree/master/Simulation) + * [DetSimInterface: Gaudi Tool interface](https://github.com/cepc/CEPCSW/tree/master/Simulation/DetSimInterface) + * [DetSimCore: integrate Gaudi and Geant4](https://github.com/cepc/CEPCSW/tree/master/Simulation/DetSimCore) + * [DetSimGeom: integrate with DD4hep](https://github.com/cepc/CEPCSW/tree/master/Simulation/DetSimGeom) + * [DetSimAna: collect data from Geant4](https://github.com/cepc/CEPCSW/tree/master/Simulation/DetSimAna) + * [DetSimSD: detector response](https://github.com/cepc/CEPCSW/tree/master/Simulation/DetSimSD) +* Job options: [See Examples/options](https://github.com/cepc/CEPCSW/tree/master/Examples/options) + +--- + +### Run simulation in CEPCSW + +* The simulation is run by following command: +```bash +$ ./run.sh Examples/options/tut_detsim.py +``` +* The job option: `Examples/options/tut_detsim.py`. <!-- .element: style="color:yellow; " --> +* :writing_hand: copy the job option into your current directory. Edit your job option in the later exercises. +```bash +$ cp Examples/options/tut_detsim.py my_detsim.py +``` + +--- + +### What's inside the job option? +* Random Number Service + * Use `Seeds` option to control the random number sequences. +* Event Data Service and PODIO writer +* Geometry Service + * Different detector options could be loaded here via `compact` option. +* Physics generator algorithm +* Detector simulation algorithm + + +--- + +#### Save detector response into ROOT file +```python +from Configurables import PodioOutput +out = PodioOutput("outputalg") +out.filename = "test-detsim10-seed42.root" +out.outputCommands = ["keep *"] +``` +* The EDM4hep format is used in the detector response. +* All the collections created in simulation will be saved. +* :writing_hand: modify the output file name. + +--- + +#### Control how many events to be simulated +```python +from Configurables import ApplicationMgr +ApplicationMgr( TopAlg = [genalg, detsimalg, out], + EvtSel = 'NONE', + EvtMax = 10, + ExtSvc = [rndmengine, dsvc, geosvc], +) +``` +* :writing_hand: modify the `EvtMax` property and check the entries in the output. + + +--- + +#### Random Number +```python +from Configurables import RndmGenSvc, HepRndm__Engine_CLHEP__RanluxEngine_ + +# rndmengine = HepRndm__Engine_CLHEP__RanluxEngine_() # The default engine in Gaudi +rndmengine = HepRndm__Engine_CLHEP__HepJamesRandom_() # The default engine in Geant4 +rndmengine.SetSingleton = True +rndmengine.Seeds = [42] + +``` + +* Seed is used to initialize the state of the random number engine. +* If two job set the same seed, the outputs will be same. +* :writing_hand: modify the seed and see the difference. + +--- + +#### Geometry / Detector Description +```python +geometry_option = "CepC_v4-onlyVXD.xml" + +geometry_path = os.path.join(os.getenv("DETCEPCV4ROOT"), + "compact", geometry_option) + +from Configurables import GeoSvc +geosvc = GeoSvc("GeoSvc") +geosvc.compact = geometry_path +``` +* The compact file is in XML format, which describes the detector. +* :writing_hand: change the geometry path and run simulation again. +``` +geometry_path = "Detector/DetEcalMatrix/compact/det.xml" +``` + +--- + +#### Customize primary particles +##### Particle Gun +```python +gun = GtGunTool("GtGunTool") +gun.Particles = ["pi+", "pi-"] +gun.EnergyMins = [100., 100] # GeV +gun.EnergyMaxs = [100., 100] # GeV + +gun.ThetaMins = [0, 0] # deg +gun.ThetaMaxs = [180., 180] # deg + +gun.PhiMins = [0., 0.] # deg +gun.PhiMaxs = [360., 360.] # deg +``` + +* Particle name can be found in [`$ROOTSYS/etc/pdg_table.txt`](https://github.com/root-project/root/blob/master/etc/pdg_table.txt) +* :writing_hand: change the particles, energies and directions. + +---- + +##### Event Generators +```python +stdheprdr = StdHepRdr("StdHepRdr") +stdheprdr.Input = "/cefs/data/stdhep/CEPC250/2fermions/E250.Pbhabha.e0.p0.whizard195/bhabha.e0.p0.00001.stdhep" + +genalg = GenAlgo("GenAlgo") +genalg.GenTools = ["StdHepRdr"] + +``` +* There are several readers to read the output of event generators in different formats + * StdHep: `StdHepRdr`, lcio: `SLCIORdr`, HepMC: `HepMCRdr`. +* :writing_hand: use the different readers to load different samples. +* [The existing samples could be found here.](http://cepcsoft.ihep.ac.cn/guides/Generation/docs/ExistingSamples/) + +--- + +#### Customize Geant4 using [built-in commands](http://geant4-userdoc.web.cern.ch/geant4-userdoc/UsersGuides/ForApplicationDeveloper/html/Control/AllResources/Control/UIcommands/_.html) +##### :writing_hand: Turn on the verbose during tracking +```python +detsimalg.RunCmds = [ + "/tracking/verbose 1", +] +# Or +detsimalg.RunMacs = [ + "run.mac", +] +``` +```bash +# Below is the content of run.mac +/tracking/verbose 1 +``` + +* Each step will be print out. Remeber to redirect the output to a file. +```bash +$ ./run.sh my_detsim.py >& mylog +``` + +---- + +##### Geant4 tracking output +``` +********************************************************************************************************* +* G4Track Information: Particle = gamma, Track ID = 1, Parent ID = 0 +********************************************************************************************************* + +Step# X(mm) Y(mm) Z(mm) KinE(MeV) dE(MeV) StepLeng TrackLeng NextVolume ProcName + 0 0 0 0 1e+05 0 0 0 pWorld initStep + 1 1.03e+03 1.8e+04 1.1e-12 0 0 1.8e+04 1.8e+04 pWorld conv + +********************************************************************************************************* +* G4Track Information: Particle = e+, Track ID = 4, Parent ID = 1 +********************************************************************************************************* + +Step# X(mm) Y(mm) Z(mm) KinE(MeV) dE(MeV) StepLeng TrackLeng NextVolume ProcName + 0 1.03e+03 1.8e+04 1.1e-12 2.91e+04 0 0 0 pWorld initStep + 1 1.04e+03 1.8e+04 -0.00523 2.91e+04 0.0117 54.9 54.9 pWorld eIoni + 2 1.04e+03 1.81e+04 -0.00912 2.91e+04 0.00604 38.7 93.6 pWorld eIoni + 3 1.06e+03 1.84e+04 -0.0404 2.91e+04 0.0634 321 415 pWorld eIoni + 4 1.06e+03 1.84e+04 -0.0416 2.91e+04 0.00234 13.1 428 pWorld eIoni + +``` +<!-- .element: style="color:yellow; font-size:small" --> + +* From this output, you can see the current track and its stepping information. +* Particle name, current track ID, parent track ID +* Step position, deposit energy + +---- + +##### :writing_hand: Visualize using Geant4 + +* Enable following command in your job option: +```python +detsimalg.VisMacs = ["Examples/options/vis.mac"] +``` +* :warning: If your X Server supports the G4 OpenGL, the detector will be shown. + * Try [Xming X Server](https://sourceforge.net/projects/xming/) in Windows. +* :notebook: [Visualization in Geant4 Documentation](http://geant4-userdoc.web.cern.ch/geant4-userdoc/UsersGuides/ForApplicationDeveloper/html/Visualization/visualization.html) + * G4 UI commands during visualization +``` +/vis/scene/add/axes 0 0 0 3 m +/vis/scene/add/magneticField +``` + +---- + + +##### Snapshot: The Qt based Geant4 visualization + + + +:writing_hand: Play with Geant4 Visualization +<!-- .element: style="color:yellow" --> + +--- + +### Analyze the simulation output +* :writing_hand: Modify the geometry option and run the simulation +```python +geometry_option = "CepC_v4-onlyECAL.xml" +``` +* :writing_hand: Plot the `EcalBarrelCollection` in ROOT +``` +root [] events->Draw("EcalBarrelCollection.position.y:EcalBarrelCollection.position.x") +root [] events->Draw("EcalBarrelCollection.position.y:EcalBarrelCollection.position.x", "Entry$==0") +``` +<!-- .element: style="color:yellow; font-size:small" --> + +---- + +#### See the branches in the `events` tree +``` +root [] events->Print() +*............................................................................* +*Br 146 :EcalBarrelCollection.cellID : * +* | ULong64_t cellID[EcalBarrelCollection_] * +*Entries : 10 : Total Size= 1398779 bytes File Size = 420214 * +*Baskets : 4 : Basket Size= 32000 bytes Compression= 3.33 * +*............................................................................* +*Br 147 :EcalBarrelCollection.energy : Float_t energy[EcalBarrelCollection_]* +*Entries : 10 : Total Size= 699855 bytes File Size = 163107 * +*Baskets : 3 : Basket Size= 32000 bytes Compression= 4.29 * +*............................................................................* +*Br 148 :EcalBarrelCollection.position.x : Float_t x[EcalBarrelCollection_] * +*Entries : 10 : Total Size= 699865 bytes File Size = 466951 * +*Baskets : 3 : Basket Size= 32000 bytes Compression= 1.50 * +*............................................................................* +*Br 149 :EcalBarrelCollection.position.y : Float_t y[EcalBarrelCollection_] * +*Entries : 10 : Total Size= 699865 bytes File Size = 469393 * +*Baskets : 3 : Basket Size= 32000 bytes Compression= 1.49 * +*............................................................................* +*Br 150 :EcalBarrelCollection.position.z : Float_t z[EcalBarrelCollection_] * +*Entries : 10 : Total Size= 699865 bytes File Size = 560229 * +*Baskets : 3 : Basket Size= 32000 bytes Compression= 1.25 * +*............................................................................* +``` +<!-- .element: style="color:yellow; font-size:small" --> + + +--- + + +# Thank you for your attention + +* [Create issue](https://github.com/cepc/CEPCSW/issues): Report a bug + +* [Pull Request](https://github.com/cepc/CEPCSW/pulls): Fix a bug or Implement a feature + +Your contributions are welcome! +</div> + </div> + + <div id="meta">{"type":"slide","slideOptions":{"transition":"slide","controls":true,"showSlideNumber":"all","navigationMode":"linear","center":false}}</div> + + <div class="footer"> + <div class="unselectable hidden-print gray-font"> + <small> + <span> + + <span class="ui-lastchangeuser"> <i class="ui-user-icon small" style="background-image: url(https://www.gravatar.com/avatar/a696972a0dc29e70843f619695acb1e9?s=80&d=identicon);" data-toggle="tooltip" data-placement="right" title="lintao@ihep.ac.cn"></i></span> + + <span class="text-uppercase ui-status-lastchange"></span> + <span class="ui-lastchange text-uppercase" data-createtime="Mon Sep 14 2020 21:17:56 GMT+0800 (China Standard Time)" data-updatetime="Wed Sep 16 2020 23:50:34 GMT+0800 (China Standard Time)"></span> + </span> + <span class="pull-right">230 views <a href="#" class="ui-edit" title="Edit this note"><i class="fa fa-fw fa-pencil"></i></a><a class="ui-print" href="" title="Open print view"><i class="fa fa-fw fa-print"></i></a></span> + <br> + + </small> + </div> + + </div> + </div> + + <script src="https://jupyter.ihep.ac.cn/js/mathjax-config-extra.js"></script> + + <script src="https://cdn.jsdelivr.net/npm/reveal.js@3.9.2/js/reveal.min.js" integrity="sha256-1fq1NvUmkMIWOBgIEzGFr0UUNuwWmOa29YqMkXnYlH4=" crossorigin="anonymous"></script> + <script src="https://cdnjs.cloudflare.com/ajax/libs/jquery/3.4.1/jquery.min.js" integrity="sha256-CSXorXvZcTkaix6Yvo6HppcZGetbYMGWSFlBw8HfCJo=" crossorigin="anonymous"></script> + <script src="https://cdnjs.cloudflare.com/ajax/libs/velocity/1.5.2/velocity.min.js" integrity="sha256-1HqoI76JGKA17K0C0s9K8L/iy8PAC43KVLt1hRD/Ojc=" crossorigin="anonymous" defer></script> + <script src="https://cdnjs.cloudflare.com/ajax/libs/jquery-mousewheel/3.1.13/jquery.mousewheel.min.js" integrity="sha256-jnOjDTXIPqall8M0MyTSt98JetJuZ7Yu+1Jm7hLTF7U=" crossorigin="anonymous" defer></script> + <script src="https://cdnjs.cloudflare.com/ajax/libs/js-yaml/3.13.1/js-yaml.min.js" integrity="sha256-ry6nlLQ1JmRl5RUPQ4eSuaSp/rGNPvl144WHHs7BiNE=" crossorigin="anonymous" defer></script> + <script src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.0/MathJax.js" integrity="sha256-yYfngbEKv4RENfGDvNUqJTqGFcKf31NJEe9OTnnMH3Y=" crossorigin="anonymous" defer></script> + <script src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.0/config/TeX-AMS-MML_HTMLorMML.js" integrity="sha256-immzXfCGLhnx3Zfi9F/dUcqxEM8K3o3oTFy9Bh6HCwg=" crossorigin="anonymous" defer></script> + <script src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.0/config/Safe.js" integrity="sha256-0ygBUDksNDXZS4vm5HMNH1a33KUu6QT1cdNTN+ZLF+4=" crossorigin="anonymous" defer></script> + <script src="https://cdnjs.cloudflare.com/ajax/libs/moment.js/2.24.0/moment-with-locales.min.js" integrity="sha256-AdQN98MVZs44Eq2yTwtoKufhnU+uZ7v2kXnD5vqzZVo=" crossorigin="anonymous" defer></script> + <script src="https://cdnjs.cloudflare.com/ajax/libs/mermaid/8.4.8/mermaid.min.js" integrity="sha256-lyWCDMnMeZiXRi7Zl54sZGKYmgQs4izcT7+tKc+KUBk=" crossorigin="anonymous" defer></script> + <script src="https://cdn.jsdelivr.net/npm/@hackmd/emojify.js@2.1.0/dist/js/emojify-browser.min.js" integrity="sha256-swgfXtqk2bC98KzPoE8tXRz5tmrzpjJWhhXlhYo/wRA=" crossorigin="anonymous" defer></script> + <script src="https://cdnjs.cloudflare.com/ajax/libs/handlebars.js/4.1.2/handlebars.min.js" integrity="sha256-ngJY93C4H39YbmrWhnLzSyiepRuQDVKDNCWO2iyMzFw=" crossorigin="anonymous" defer></script> + <script src="https://cdnjs.cloudflare.com/ajax/libs/highlight.js/9.15.10/highlight.min.js" integrity="sha256-1zu+3BnLYV9LdiY85uXMzii3bdrkelyp37e0ZyTAQh0=" crossorigin="anonymous" defer></script> + <script src="https://cdnjs.cloudflare.com/ajax/libs/gist-embed/2.6.0/gist-embed.min.js" integrity="sha256-KyF2D6xPIJUW5sUDSs93vWyZm+1RzIpKCexxElmxl8g=" crossorigin="anonymous" defer></script> + <script src="https://cdnjs.cloudflare.com/ajax/libs/viz.js/2.1.2/viz.js" integrity="sha256-8RHyK+AFzq9iXwbFo2unqidwPbwHU5FFWe3RwkcVtuU=" crossorigin="anonymous" defer></script> + <script src="https://cdnjs.cloudflare.com/ajax/libs/viz.js/2.1.2/full.render.js" integrity="sha256-Ogqs510LFnekr9o7OLdpelaaAmNss9egQRTyzCqV2NQ=" crossorigin="anonymous" defer></script> + <script src="https://cdnjs.cloudflare.com/ajax/libs/abcjs/3.1.1/abcjs_basic-min.js" integrity="sha256-Sq1r2XXWXQoShQKsS0Wrf5r7fRkErd9Fat9vHYeU68s=" crossorigin="anonymous" defer></script> + <script src="https://cdnjs.cloudflare.com/ajax/libs/vega/5.9.1/vega.min.js" integrity="sha256-xVmd2OiOTh73s2iPfGy1DNyu/lCKvaDto452MU1O+xs=" crossorigin="anonymous" defer></script> + <script src="https://cdnjs.cloudflare.com/ajax/libs/vega-lite/4.4.0/vega-lite.min.js" integrity="sha256-ollz/GSuG0/f7aV4v8LGDYxPs4G2DwEk9+hALicqp9I=" crossorigin="anonymous" defer></script> + <script src="https://cdnjs.cloudflare.com/ajax/libs/vega-embed/6.2.2/vega-embed.min.js" integrity="sha256-AW13lGYqQzWT9PymwqUEJqQHaz9ntM5m5jQVkvtzja4=" crossorigin="anonymous" defer></script> + <script src="https://cdnjs.cloudflare.com/ajax/libs/leaflet/1.6.0/leaflet.js" integrity="sha256-fNoRrwkP2GuYPbNSJmMJOCyfRB2DhPQe0rGTgzRsyso=" crossorigin="anonymous" defer></script> + <script src="https://jupyter.ihep.ac.cn/config"></script><script src="https://jupyter.ihep.ac.cn/build/slide.857c50b01c7ad0375da9.js" defer="defer"></script> + + </body> +</html> + + + diff --git a/docs/tutorial.md b/docs/tutorial.md index 6dac5d48a290820cdcb0e28b776234b842625b38..ed53bf9d1b95e1ab587c27d3db3a0b5325153ffc 100644 --- a/docs/tutorial.md +++ b/docs/tutorial.md @@ -1,3 +1,4 @@ # Tutorial * [New CEPCSW Tutorial and detector study, IHEP, 17-18 Sept 2020](https://indico.ihep.ac.cn/event/12341/) + * [Tutorial on CEPCSW simulation](https://cepc.github.io/CEPCSW/simulation-tutorial/index.html#/)