diff --git a/Detector/DetCRD/scripts/TDR_o1_v01/calodigi.py b/Detector/DetCRD/scripts/TDR_o1_v01/calodigi.py index 46e93c7cb8f1c15bcabe7a92dccfcdb5b311f971..b90d25ae67adfccf5589f6ba4d0a5808adff96d1 100644 --- a/Detector/DetCRD/scripts/TDR_o1_v01/calodigi.py +++ b/Detector/DetCRD/scripts/TDR_o1_v01/calodigi.py @@ -3,7 +3,7 @@ import os from Gaudi.Configuration import * from Configurables import k4DataSvc -dsvc = k4DataSvc("EventDataSvc", input="rec00.root") +dsvc = k4DataSvc("EventDataSvc", input="rec_v01.root") from Configurables import RndmGenSvc, HepRndm__Engine_CLHEP__RanluxEngine_ seed = [12340] diff --git a/Detector/DetCRD/scripts/TDR_o1_v01/rec.py b/Detector/DetCRD/scripts/TDR_o1_v01/rec.py index 362421ef3cd9fac5eac067ce2a38173bd33f061d..172e8e938d849cdd4745d17470e2bd7c32a99266 100644 --- a/Detector/DetCRD/scripts/TDR_o1_v01/rec.py +++ b/Detector/DetCRD/scripts/TDR_o1_v01/rec.py @@ -47,7 +47,9 @@ inp.collections = [ # "HCALEndcapsParticleAssoCol", "MCParticle", "CompleteTracks", - "CompleteTracksParticleAssociation"] + "CompleteTracksParticleAssociation", + "RecTofCollection", + "DndxTracks"] ########################################## ######### Reconstruction ################ @@ -140,6 +142,10 @@ CyberPFAlg.AlgParValues = [ ["BarCol","Cluster1DCol","HalfClusterCol"],#1 ["EcalCluster", "SimpleHCALCluster", "outputPFO"], #16 ["1.26","4.", "1.", "4."] ]#17 +from Configurables import FinalPIDAlg +pid = FinalPIDAlg("FinalPIDAlg") +pid.OutputPFOName = "CyberPFOPID" + from Configurables import GenMatch genmatch = GenMatch("GenMatch") @@ -159,7 +165,7 @@ out.outputCommands = ["keep *"] from Configurables import ApplicationMgr ApplicationMgr( - TopAlg=[inp, CyberPFAlg, out ], + TopAlg=[inp, CyberPFAlg, pid, out ], EvtSel="NONE", EvtMax=5, ExtSvc=[podioevent, geomsvc], diff --git a/Detector/DetCRD/scripts/TDR_o1_v01/sim.py b/Detector/DetCRD/scripts/TDR_o1_v01/sim.py index 575297e7e073d4264a4ad772109ebe77c14f2bce..ff650d725fb30ce53b5c19bdf2aad87e5142e07d 100644 --- a/Detector/DetCRD/scripts/TDR_o1_v01/sim.py +++ b/Detector/DetCRD/scripts/TDR_o1_v01/sim.py @@ -90,7 +90,7 @@ evtseeder = MarlinEvtSeeder("EventSeeder") # output from Configurables import PodioOutput out = PodioOutput("outputalg") -out.filename = "sim00.root" +out.filename = "sim_v01.root" out.outputCommands = ["keep *"] # ApplicationMgr diff --git a/Detector/DetCRD/scripts/TDR_o1_v01/tracking.py b/Detector/DetCRD/scripts/TDR_o1_v01/tracking.py index a5f9324f580f7227a234858d914f6936b6985053..b5a94d6084f79d124b6a63ee872c37236b74f1c9 100644 --- a/Detector/DetCRD/scripts/TDR_o1_v01/tracking.py +++ b/Detector/DetCRD/scripts/TDR_o1_v01/tracking.py @@ -3,7 +3,7 @@ import os from Gaudi.Configuration import * from Configurables import k4DataSvc -dsvc = k4DataSvc("EventDataSvc", input="sim00.root") +dsvc = k4DataSvc("EventDataSvc", input="sim_v01.root") from Configurables import RndmGenSvc, HepRndm__Engine_CLHEP__RanluxEngine_ seed = [12340] @@ -253,6 +253,10 @@ from Configurables import TPCDndxAlg tpc_dndx = TPCDndxAlg("TPCDndxAlg") tpc_dndx.Method = "Simple" +from Configurables import TofRecAlg +tof = TofRecAlg("TofRecAlg") +#tof.OutputLevel = DEBUG + from Configurables import TrackParticleRelationAlg tpr = TrackParticleRelationAlg("Track2Particle") tpr.MCParticleCollection = "MCParticle" @@ -272,13 +276,13 @@ tmt.MuonDetTanTheta = 1.2 # muon det barrel/endcap separation tan(theta) # output from Configurables import PodioOutput out = PodioOutput("outputalg") -out.filename = "rec00.root" +out.filename = "rec_v01.root" out.outputCommands = ["keep *"] # ApplicationMgr from Configurables import ApplicationMgr mgr = ApplicationMgr( - TopAlg = [podioinput, digiVXD, digiSIT, digiOTKB, digiFTD, digiTPC, digiMuon, tracking, forward, subset, clupatra, full, tpr, tpc_dndx, tmt, out], + TopAlg = [podioinput, digiVXD, digiSIT, digiOTKB, digiFTD, digiTPC, digiMuon, tracking, forward, subset, clupatra, full, tpr, tpc_dndx, tof, tmt, out], EvtSel = 'NONE', EvtMax = 50, ExtSvc = [rndmengine, rndmgensvc, dsvc, evtseeder, geosvc, gearsvc, tracksystemsvc, pidsvc], diff --git a/Reconstruction/ParticleID/CMakeLists.txt b/Reconstruction/ParticleID/CMakeLists.txt index 4c03899bc8eef5dbc0dc61450f05e5823eb82d41..7d825c595cf255bc691189248ae5bf972a9210ee 100644 --- a/Reconstruction/ParticleID/CMakeLists.txt +++ b/Reconstruction/ParticleID/CMakeLists.txt @@ -5,6 +5,7 @@ gaudi_add_module(ParticleID SOURCES src/TPCDndxAlg.cpp src/DCHDndxAlg.cpp src/TofRecAlg.cpp + src/FinalPIDAlg.cpp LINK SimplePIDSvc Gaudi::GaudiAlgLib Gaudi::GaudiKernel diff --git a/Reconstruction/ParticleID/script/rec.py b/Reconstruction/ParticleID/script/rec.py new file mode 100644 index 0000000000000000000000000000000000000000..2765f73eac038f2269c8555f9ba9d56d6ad2a55c --- /dev/null +++ b/Reconstruction/ParticleID/script/rec.py @@ -0,0 +1,168 @@ +import os, sys +from Gaudi.Configuration import * + +############## GeomSvc ################# +geometry_option = "TDR_o1_v01/TDR_o1_v01.xml" + +if not os.getenv("DETCRDROOT"): + print("Can't find the geometry. Please setup envvar DETCRDROOT." ) + sys.exit(-1) + +geometry_path = os.path.join(os.getenv("DETCRDROOT"), "compact", geometry_option) +if not os.path.exists(geometry_path): + print("Can't find the compact geometry file: %s"%geometry_path) + sys.exit(-1) + +from Configurables import DetGeomSvc +geomsvc = DetGeomSvc("GeomSvc") +geomsvc.compact = geometry_path +####################################### + +########### k4DataSvc #################### +from Configurables import k4DataSvc +podioevent = k4DataSvc("EventDataSvc", input="Tracking_TDR_o1_v01.root") +########################################## + +########## CEPCSWData ################# +cepcswdatatop ="/cvmfs/cepcsw.ihep.ac.cn/prototype/releases/data/latest" +####################################### + +########## CrystalEcalEnergyCorrectionSvc ######## +from Configurables import CrystalEcalEnergyCorrectionSvc +crystalecalcorr = CrystalEcalEnergyCorrectionSvc("CrystalEcalEnergyCorrectionSvc") +crystalecalcorr.CorrectionFile = os.path.join(cepcswdatatop, "CEPCSWData/offline-data/Service/CrystalEcalSvc/data/CrackRegionEnergyCorrection.root") +################################################## + +########## Podio Input ################### +from Configurables import PodioInput +inp = PodioInput("InputReader") +inp.collections = [ + "ECALBarrel", + "ECALBarrelParticleAssoCol", +# "ECALEndcaps", +# "ECALEndcapsParticleAssoCol", + "HCALBarrel", + "HCALBarrelParticleAssoCol", +# "HCALEndcaps", +# "HCALEndcapsParticleAssoCol", + "MCParticle", + "CompleteTracks", + "CompleteTracksParticleAssociation", + "RecTofCollection", + "DndxTracks"] +########################################## + +######### Reconstruction ################ +from Configurables import CyberPFAlg +CyberPFAlg = CyberPFAlg("CyberPFAlg") +##----Global parameters---- +CyberPFAlg.Seed = 1024 +CyberPFAlg.BField = 3. +CyberPFAlg.Debug = 0 +CyberPFAlg.SkipEvt = 0 +CyberPFAlg.WriteAna = 1 +CyberPFAlg.AnaFileName = "RecAnaTuple_TDR_o1_v01.root" +CyberPFAlg.UseMCPTrack = 0 +CyberPFAlg.UseTruthMatchTrack = 0 +CyberPFAlg.DoCleanTrack = 1 +CyberPFAlg.TrackIDFile = "/cvmfs/cepcsw.ihep.ac.cn/prototype/releases/data/latest/CEPCSWData/offline-data/Reconstruction/CyberPFA_trackID/TrkID_BDT_BDTG.weights.xml" +CyberPFAlg.TrackIDMethod = "BDTG" +CyberPFAlg.EcalChargedCalib = 1.26 +CyberPFAlg.HcalChargedCalib = 4.0 +CyberPFAlg.EcalNeutralCalib = 1.0 +CyberPFAlg.HcalNeutralCalib = 4.0 +##----Readin collections---- +CyberPFAlg.MCParticleCollection = "MCParticle" +CyberPFAlg.TrackCollections = ["CompleteTracks"] +CyberPFAlg.MCRecoTrackParticleAssociationCollection = "CompleteTracksParticleAssociation" +#CyberPFAlg.ECalCaloHitCollections = ["ECALBarrel","ECALEndcaps"] +#CyberPFAlg.ECalReadOutNames = ["EcalBarrelCollection","EcalEndcapsCollection"] +#CyberPFAlg.ECalMCPAssociationName = ["ECALBarrelParticleAssoCol", "ECALEndcapsParticleAssoCol"] +#CyberPFAlg.HCalCaloHitCollections = ["HCALBarrel", "HCALEndcaps"] +#CyberPFAlg.HCalReadOutNames = ["HcalBarrelCollection", "HcalEndcapsCollection"] +#CyberPFAlg.HCalMCPAssociationName = ["HCALBarrelParticleAssoCol", "HCALEndcapsParticleAssoCol"] +CyberPFAlg.ECalCaloHitCollections = ["ECALBarrel"] +CyberPFAlg.ECalReadOutNames = ["EcalBarrelCollection"] +CyberPFAlg.ECalMCPAssociationName = ["ECALBarrelParticleAssoCol"] +CyberPFAlg.HCalCaloHitCollections = ["HCALBarrel"] +CyberPFAlg.HCalReadOutNames = ["HcalBarrelCollection"] +CyberPFAlg.HCalMCPAssociationName = ["HCALBarrelParticleAssoCol"] + +##--- Output collections --- +CyberPFAlg.OutputPFO = "outputPFO"; +CyberPFAlg.RecoPFOCollection = "CyberPFO" + +#----Algorithms---- +CyberPFAlg.AlgList = ["GlobalClusteringAlg", #1 + "LocalMaxFindingAlg", #2 + "TrackMatchingAlg", #3 + "HoughClusteringAlg", #4 + "ConeClustering2DAlg", #5 + "AxisMergingAlg", #6 + "EnergySplittingAlg", #9 + "EnergyTimeMatchingAlg", #11 + "HcalClusteringAlg", #12 + "TruthClusteringAlg", #15 + "TrackClusterConnectingAlg", #16 + "PFOReclusteringAlg" ] #17 +CyberPFAlg.AlgParNames = [ ["InputECALBars","OutputECAL1DClusters","OutputECALHalfClusters"],#1 + ["OutputLocalMaxName"],#2 + ["ReadinLocalMaxName","OutputLongiClusName"],#3 + ["ReadinLocalMaxName","LeftLocalMaxName","OutputLongiClusName"],#4 + ["ReadinLocalMaxName", "OutputLongiClusName"], #5 + ["OutputAxisName"], #6 + ["ReadinAxisName", "OutputClusName", "OutputTowerName"], #9 + ["ReadinHFClusterName", "ReadinTowerName","OutputClusterName"], #11 + ["InputHCALHits", "OutputHCALClusters"], #12 + ["DoECALClustering","DoHCALClustering","InputHCALHits","OutputHCALClusters"], #15 + ["ReadinECALClusterName", "ReadinHCALClusterName", "OutputCombPFO"], #16 + ["ECALChargedCalib", "HCALChargedCalib", "ECALNeutralCalib", "HCALNeutralCalib"] ]#17 +CyberPFAlg.AlgParTypes = [ ["string","string","string"],#1 + ["string"],#2 + ["string","string"],#3 + ["string","string","string"],#4 + ["string","string"], #5 + ["string"], #6 + ["string","string","string"], #9 + ["string","string","string"], #11 + ["string", "string"], #12 + ["bool","bool","string","string"], #15 + ["string","string","string"], #16 + ["double","double", "double","double"] ]#17 +CyberPFAlg.AlgParValues = [ ["BarCol","Cluster1DCol","HalfClusterCol"],#1 + ["AllLocalMax"],#2 + ["AllLocalMax","TrackAxis"],#3 + ["AllLocalMax","LeftLocalMax","HoughAxis"],#4 + ["LeftLocalMax","ConeAxis"], #5 + ["MergedAxis"], #6 + ["MergedAxis","ESHalfCluster","ESTower"], #9 + ["ESHalfCluster","ESTower","EcalCluster"], #11 + ["HCALBarrel", "SimpleHCALCluster"], #12 + ["0","1","HCALBarrel","HCALCluster"], #15 + ["EcalCluster", "SimpleHCALCluster", "outputPFO"], #16 + ["1.26","4.", "1.", "4."] ]#17 + + +from Configurables import FinalPIDAlg +pid = FinalPIDAlg("FinalPIDAlg") +pid.OutputPFOName = "CyberPFOPID" + +############################################################################## +# POD I/O +############################################################################## +from Configurables import PodioOutput +out = PodioOutput("outputalg") +out.filename = "Rec_TDR_o1_v01.root" +out.outputCommands = ["keep *"] + + +######################################## + +from Configurables import ApplicationMgr +ApplicationMgr( + TopAlg=[inp, CyberPFAlg, pid, out ], + EvtSel="NONE", + EvtMax=10, + ExtSvc=[podioevent, geomsvc], + #OutputLevel=DEBUG +) diff --git a/Reconstruction/ParticleID/script/tracking.py b/Reconstruction/ParticleID/script/tracking.py new file mode 100644 index 0000000000000000000000000000000000000000..dfb749c291d30fa6ffd1f4b83a021c78fa3cf018 --- /dev/null +++ b/Reconstruction/ParticleID/script/tracking.py @@ -0,0 +1,306 @@ +#!/usr/bin/env python +import os +from Gaudi.Configuration import * + +from Configurables import k4DataSvc +dsvc = k4DataSvc("EventDataSvc", input="CaloDigi_TDR_o1_v01.root") + +from Configurables import RndmGenSvc, HepRndm__Engine_CLHEP__RanluxEngine_ +seed = [1024] +# rndmengine = HepRndm__Engine_CLHEP__RanluxEngine_() # The default engine in Gaudi +rndmengine = HepRndm__Engine_CLHEP__HepJamesRandom_("RndmGenSvc.Engine") # The default engine in Geant4 +rndmengine.SetSingleton = True +rndmengine.Seeds = seed + +rndmgensvc = RndmGenSvc("RndmGenSvc") +rndmgensvc.Engine = rndmengine.name() + +geometry_option = "TDR_o1_v01/TDR_o1_v01.xml" + +if not os.getenv("DETCRDROOT"): + print("Can't find the geometry. Please setup envvar DETCRDROOT." ) + sys.exit(-1) + +geometry_path = os.path.join(os.getenv("DETCRDROOT"), "compact", geometry_option) +if not os.path.exists(geometry_path): + print("Can't find the compact geometry file: %s"%geometry_path) + sys.exit(-1) + +from Configurables import DetGeomSvc +geosvc = DetGeomSvc("GeomSvc") +geosvc.compact = geometry_path + +from Configurables import MarlinEvtSeeder +evtseeder = MarlinEvtSeeder("EventSeeder") + +from Configurables import GearSvc +gearsvc = GearSvc("GearSvc") + +from Configurables import TrackSystemSvc +tracksystemsvc = TrackSystemSvc("TrackSystemSvc") + +from Configurables import SimplePIDSvc +pidsvc = SimplePIDSvc("SimplePIDSvc") +cepcswdatatop = "/cvmfs/cepcsw.ihep.ac.cn/prototype/releases/data/latest" +pidsvc.ParFile = os.path.join(cepcswdatatop, "CEPCSWData/offline-data/Service/SimplePIDSvc/data/dNdx_TPC.root") + +from Configurables import PodioInput +podioinput = PodioInput("PodioReader", collections=[ +# "EventHeader", + "MCParticle", + "VXDCollection", + "SITCollection", + "TPCCollection", + "OTKBarrelCollection", + "FTDCollection", + "MuonBarrelCollection", + "MuonEndcapCollection" + ]) + +# digitization +vxdhitname = "VXDTrackerHits" +sithitname = "SITTrackerHits" +gashitname = "TPCTrackerHits" +sethitname = "OTKBarrelTrackerHits" +setspname = "OTKBarrelSpacePoints" +ftdhitname = "FTDTrackerHits" +ftdspname = "FTDSpacePoints" +from Configurables import SmearDigiTool +vxdtool = SmearDigiTool("VXD") +vxdtool.ResolutionU = [0.005] +vxdtool.ResolutionV = [0.005] +vxdtool.UsePlanarTag = True +vxdtool.ParameterizeResolution = False +vxdtool.ParametersU = [5.60959e-03, 5.74913e-03, 7.03433e-03, 1.99516, -663.952, 3.752e-03, 0, -0.0704734, 0.0454867e-03, 1.07359] +vxdtool.ParametersV = [5.60959e-03, 5.74913e-03, 7.03433e-03, 1.99516, -663.952, 3.752e-03, 0, -0.0704734, 0.0454867e-03, 1.07359] +#vxdtool.OutputLevel = DEBUG + +from Configurables import SiTrackerDigiAlg +digiVXD = SiTrackerDigiAlg("VXDDigi") +digiVXD.SimTrackHitCollection = "VXDCollection" +digiVXD.TrackerHitCollection = vxdhitname +digiVXD.TrackerHitAssociationCollection = "VXDTrackerHitAssociation" +digiVXD.DigiTool = "SmearDigiTool/VXD" +#digiVXD.OutputLevel = DEBUG + +from Configurables import PlanarDigiAlg +digiSIT = PlanarDigiAlg("SITDigi") +digiSIT.IsStrip = False +digiSIT.SimTrackHitCollection = "SITCollection" +digiSIT.TrackerHitCollection = sithitname +digiSIT.TrackerHitAssociationCollection = "SITTrackerHitAssociation" +digiSIT.ResolutionU = [0.0098] +digiSIT.ResolutionV = [0.0433] +digiSIT.UsePlanarTag = True +digiSIT.ParameterizeResolution = False +digiSIT.ParametersU = [2.29655e-03, 0.965899e-03, 0.584699e-03, 17.0856, 84.566, 12.4695e-03, -0.0643059, 0.168662, 1.87998e-03, 0.514452] +digiSIT.ParametersV = [1.44629e-02, 2.20108e-03, 1.03044e-02, 4.39195e+00, 3.29641e+00, 1.55167e+18, -5.41954e+01, 5.72986e+00, -6.80699e-03, 5.04095e-01] +#digiSIT.OutputLevel = DEBUG + +digiSET = PlanarDigiAlg("SETDigi") +digiSET.IsStrip = False +digiSET.SimTrackHitCollection = "OTKBarrelCollection" +digiSET.TrackerHitCollection = sethitname +digiSET.TrackerHitAssociationCollection = "OTKBarrelTrackerHitAssociation" +digiSET.ResolutionU = [0.010] +digiSET.ResolutionV = [1.000] +digiSET.UsePlanarTag = True +digiSET.ParameterizeResolution = False +digiSET.ParametersU = [2.29655e-03, 0.965899e-03, 0.584699e-03, 17.0856, 84.566, 12.4695e-03, -0.0643059, 0.168662, 1.87998e-03, 0.514452] +digiSET.ParametersV = [1.44629e-02, 2.20108e-03, 1.03044e-02, 4.39195e+00, 3.29641e+00, 1.55167e+18, -5.41954e+01, 5.72986e+00, -6.80699e-03, 5.04095e-01] +#digiSET.OutputLevel = DEBUG + +digiFTD = PlanarDigiAlg("FTDDigi") +digiFTD.IsStrip = False +digiFTD.SimTrackHitCollection = "FTDCollection" +digiFTD.TrackerHitCollection = ftdhitname +digiFTD.TrackerHitAssociationCollection = "FTDTrackerHitAssociation" +digiFTD.ResolutionU = [0.0072] +digiFTD.ResolutionV = [0.086] +digiFTD.UsePlanarTag = True +digiFTD.ParameterizeResolution = False +digiFTD.ParametersU = [2.29655e-03, 0.965899e-03, 0.584699e-03, 17.0856, 84.566, 12.4695e-03, -0.0643059, 0.168662, 1.87998e-03, 0.514452] +digiFTD.ParametersV = [1.44629e-02, 2.20108e-03, 1.03044e-02, 4.39195e+00, 3.29641e+00, 1.55167e+18, -5.41954e+01, 5.72986e+00, -6.80699e-03, 5.04095e-01] +#digiFTD.OutputLevel = DEBUG + +from Configurables import TPCDigiAlg +digiTPC = TPCDigiAlg("TPCDigi") +digiTPC.TPCCollection = "TPCCollection" +digiTPC.TPCLowPtCollection = "TPCLowPtCollection" +digiTPC.TPCTrackerHitsCol = gashitname +#default value, modify them according to future Garfield simulation results +#digiTPC.PixelClustering = True +#digiTPC.PointResolutionRPhi = 0.144 +#digiTPC.DiffusionCoeffRPhi = 0.0323 +#digiTPC.PointResolutionZ = 0.4 +#digiTPC.DiffusionCoeffZ = 0.23 +#digiTPC.N_eff = 30 +#digiTPC.OutputLevel = DEBUG + +## Muon Detector ## +from Configurables import MuonDigiAlg +digiMuon = MuonDigiAlg("MuonDigiAlg") +digiMuon.MuonBarrelHitsCollection = "MuonBarrelCollection" +digiMuon.MuonEndcapHitsCollection = "MuonEndcapCollection" +digiMuon.MuonBarrelTrackerHits = "MuonBarrelTrackerHits" +digiMuon.MuonEndcapTrackerHits = "MuonEndcapTrackerHits" +digiMuon.WriteNtuple = 0 +digiMuon.OutFileName = "Digi_MUON.root" + +# tracking +from Configurables import KalTestTool +# Close multiple scattering and smooth, used by clupatra +kt010 = KalTestTool("KalTest010") +kt010.MSOn = False +kt010.Smooth = False +#kt010.OutputLevel = DEBUG + +# Open multiple scattering, energy loss and smooth (default) +kt111 = KalTestTool("KalTest111") +#kt111.OutputLevel = DEBUG + +# Close smooth +kt110 = KalTestTool("KalTest110") +kt110.Smooth = False +#kt110.OutputLevel = DEBUG + +from Configurables import SiliconTrackingAlg +tracking = SiliconTrackingAlg("SiliconTracking") +tracking.LayerCombinationsFTD = [] +tracking.HeaderCol = "EventHeader" +tracking.VTXHitCollection = vxdhitname +tracking.SITHitCollection = sithitname +tracking.FTDPixelHitCollection = ftdhitname +tracking.FTDSpacePointCollection = ftdspname +tracking.SITRawHitCollection = "NotNeedForPixelSIT" +tracking.FTDRawHitCollection = ftdhitname +tracking.UseSIT = True +tracking.SmoothOn = False +tracking.NDivisionsInTheta = 10 +tracking.NDivisionsInPhi = 60 +tracking.NDivisionsInPhiFTD = 16 +tracking.MinDistCutAttach = 50 +# for p=1GeV, theta=10degree, Chi2FitCut = 1500, HelixMaxChi2 = 1000000, Chi2WZ = 0.02 +tracking.Chi2FitCut = 200 +tracking.MaxChi2PerHit = 200 +tracking.HelixMaxChi2 = 50000 +tracking.Chi2WZTriplet = 0.1 +tracking.Chi2WZQuartet = 0.1 +tracking.Chi2WZSeptet = 0.1 +#tracking.FitterTool = "KalTestTool/KalTest111" +#tracking.OutputLevel = DEBUG + +from Configurables import ForwardTrackingAlg +forward = ForwardTrackingAlg("ForwardTracking") +forward.FTDPixelHitCollection = ftdhitname +forward.FTDSpacePointCollection = ftdspname +forward.FTDRawHitCollection = ftdhitname +forward.Chi2ProbCut = 0.0 +forward.HitsPerTrackMin = 3 +forward.BestSubsetFinder = "SubsetSimple" +forward.Criteria = ["Crit2_DeltaPhi","Crit2_StraightTrackRatio","Crit3_3DAngle","Crit3_ChangeRZRatio","Crit3_IPCircleDist","Crit4_3DAngleChange","Crit4_DistToExtrapolation", + "Crit2_DeltaRho","Crit2_RZRatio","Crit3_PT"] +forward.CriteriaMin = [0, 0.9, 0, 0.995, 0, 0.8, 0, 20, 1.002, 0.1, 0, 0.99, 0, 0.999, 0, 0.99, 0] +forward.CriteriaMax = [30, 1.02, 10, 1.015, 20, 1.3, 1.0, 150, 1.08, 99999999, 0.8, 1.01, 0.35, 1.001, 1.5, 1.01, 0.05] +#forward.FitterTool = "KalTestTool/KalTest110" +#forward.OutputLevel = DEBUG + +from Configurables import TrackSubsetAlg +subset = TrackSubsetAlg("TrackSubset") +subset.TrackInputCollections = ["ForwardTracks", "SiTracks"] +subset.RawTrackerHitCollections = [vxdhitname, sithitname, ftdhitname, ftdspname] +subset.TrackSubsetCollection = "SubsetTracks" +#subset.FitterTool = "KalTestTool/KalTest111" +#subset.OutputLevel = DEBUG + +from Configurables import ClupatraAlg +clupatra = ClupatraAlg("Clupatra") +clupatra.TPCHitCollection = gashitname +#clupatra.OutputLevel = DEBUG + +from Configurables import FullLDCTrackingAlg +full = FullLDCTrackingAlg("FullTracking") +full.VTXTrackerHits = vxdhitname +full.SITTrackerHits = sithitname +full.TPCTrackerHits = gashitname +full.SETTrackerHits = sethitname +full.FTDPixelTrackerHits = ftdhitname +full.FTDSpacePoints = ftdspname +full.SITRawHits = "NotNeedForPixelSIT" +full.SETRawHits = "NotNeedForPixelSET" +full.FTDRawHits = ftdhitname +full.TPCTracks = "ClupatraTracks" # add standalone TPC track +full.SiTracks = "SubsetTracks" +full.OutputTracks = "CompleteTracks" # default name +#full.VTXHitToTrackDistance = 5. +full.FTDHitToTrackDistance = 5. +full.SITHitToTrackDistance = 3. +full.SETHitToTrackDistance = 5. +full.MinChi2ProbForSiliconTracks = 0 +full.MaxChi2PerHit = 200 +full.ForceSiTPCMerging = True +full.ForceTPCSegmentsMerging = True +#full.OutputLevel = DEBUG + +from Configurables import TPCDndxAlg +tpc_dndx = TPCDndxAlg("TPCDndxAlg") +tpc_dndx.Method = "Simple" + +from Configurables import TofRecAlg +tof = TofRecAlg("TofRecAlg") +#tof.OutputLevel = DEBUG + +from Configurables import TrackParticleRelationAlg +tpr = TrackParticleRelationAlg("Track2Particle") +tpr.MCParticleCollection = "MCParticle" +tpr.TrackList = ["CompleteTracks", "ClupatraTracks"] +tpr.TrackerAssociationList = ["VXDTrackerHitAssociation", "SITTrackerHitAssociation", "OTKBarrelTrackerHitAssociation", "FTDTrackerHitAssociation", "TPCTrackerHitAss"] +#tpr.OutputLevel = DEBUG + +from Configurables import TrueMuonTagAlg +tmt = TrueMuonTagAlg("TrueMuonTag") +tmt.MCParticleCollection = "MCParticle" +tmt.TrackList = ["CompleteTracks"] +tmt.TrackerAssociationList = ["VXDTrackerHitAssociation", "SITTrackerHitAssociation", "OTKBarrelTrackerHitAssociation", "FTDTrackerHitAssociation", "TPCTrackerHitAss"] +tmt.MuonTagEfficiency = 0.95 # muon true tag efficiency, default is 1.0 (100%) +tmt.MuonDetTanTheta = 1.2 # muon det barrel/endcap separation tan(theta) +#tmt.OutputLevel = DEBUG + +#from Configurables import ReadDigiAlg +#readtrk = ReadDigiAlg("ReadDigiAlg") +#readtrk.SiTracks = "SubsetTracks" +#readtrk.TPCTracks = "ClupatraTracks" +#readtrk.FullTracks = "CompleteTracks" +#readtrk.TPCTracksAssociation = "ClupatraTracksParticleAssociation" +#readtrk.FullTracksAssociation = "CompleteTracksParticleAssociation" +#readtrk.OutFileName = "TrackAnaTuple_mu.root" + +# output +from Configurables import PodioOutput +out = PodioOutput("outputalg") +out.filename = "Tracking_TDR_o1_v01.root" +out.outputCommands = ["drop *", + "keep ECALBarrel", + "keep ECALBarrelParticleAssoCol", + "keep HCALBarrel", + "keep HCALBarrelParticleAssoCol", + "keep ECALEndcaps", + "keep HCALEndcaps", + "keep ECALEndcapsParticleAssoCol", + "keep HCALEndcapsParticleAssoCol", + "keep MCParticle", + "keep CompleteTracks", + "keep CompleteTracksParticleAssociation", + "keep RecTofCollection", + "keep DndxTracks" ] + +# ApplicationMgr +from Configurables import ApplicationMgr +mgr = ApplicationMgr( + TopAlg = [podioinput, digiVXD, digiSIT, digiSET, digiFTD, digiTPC, digiMuon, tracking, forward, subset, clupatra, full, tpr, tpc_dndx, tof, tmt, out], + EvtSel = 'NONE', + EvtMax = 10, + ExtSvc = [rndmengine, rndmgensvc, dsvc, evtseeder, geosvc, gearsvc, tracksystemsvc, pidsvc], + #HistogramPersistency = 'ROOT', + OutputLevel = INFO +) diff --git a/Reconstruction/ParticleID/src/FinalPIDAlg.cpp b/Reconstruction/ParticleID/src/FinalPIDAlg.cpp new file mode 100644 index 0000000000000000000000000000000000000000..b1dc72f126b6590f278a7679e4620fecc36cff77 --- /dev/null +++ b/Reconstruction/ParticleID/src/FinalPIDAlg.cpp @@ -0,0 +1,156 @@ +#include "FinalPIDAlg.h" + +#include "GaudiKernel/DataObject.h" +#include "GaudiKernel/MsgStream.h" +#include "GaudiKernel/SmartDataPtr.h" +#include "DetInterface/IGeomSvc.h" + +using namespace edm4hep; + +DECLARE_COMPONENT( FinalPIDAlg ) + +//------------------------------------------------------------------------------ +FinalPIDAlg::FinalPIDAlg( const std::string& name, ISvcLocator* pSvcLocator ) + : Algorithm( name, pSvcLocator ) { + // input + declareProperty("PFOCollection", m_inPFOCol, "Reconstructed particles to be PIDed"); + declareProperty("TOFPIDCollection", m_inTofCol, "TOF PID information"); + declareProperty("TPCPIDCollection", m_inDqdxCol, "TPC PID information"); + // output + declareProperty("OutputPFOName", m_outPFOCol, "Reconstructed particles with PID information"); + // PID method + declareProperty("Method", m_method = "TPC+TOF", "PID method: TPC, TPC+TOF, TPC+TOF+ECAL"); +} + +//------------------------------------------------------------------------------ +StatusCode FinalPIDAlg::initialize(){ + debug() << "Initializing..." << endmsg; + if (m_method == "TPC" || m_method == "TPC+TOF") { //|| m_method == "TPC+TOF+ECAL" ) { // add ECAL later + debug() << "PID method: " << m_method << endmsg; + } else { + error() << "Unknown PID method: " << m_method << endmsg; + return StatusCode::FAILURE; + } + _nEvt = 0; + return StatusCode::SUCCESS; +} + +//------------------------------------------------------------------------------ +StatusCode FinalPIDAlg::execute(){ + + const edm4hep::ReconstructedParticleCollection* inpfocol = nullptr; + const edm4hep::RecTofCollection* tofcol = nullptr; + const edm4hep::RecDqdxCollection* dqdxcol = nullptr; + _hasTOF = true; + _hasTPC = true; + + //auto outpfocol = m_outPFOCol.createAndPut(); + edm4hep::ReconstructedParticleCollection* outpfocol = m_outPFOCol.createAndPut(); + + try { + inpfocol = m_inPFOCol.get(); + } + catch ( GaudiException &e ) { + debug() << " Reconstructed particle " << m_inPFOCol.fullKey() << " is unavailable in event " << _nEvt << endmsg; + _nEvt++; + return StatusCode::SUCCESS; + } + try { + tofcol = m_inTofCol.get(); + } + catch ( GaudiException &e ) { + _hasTOF = false; + } + try { + dqdxcol = m_inDqdxCol.get(); + } + catch ( GaudiException &e ) { + _hasTPC = false; + } + debug()<<" has TOF : "<<_hasTOF<<" has TPC : "<<_hasTPC<<" TOF size : "<<tofcol->size()<<" dqdx size : "<<dqdxcol->size()<<endmsg; + + if ( ( !_hasTPC && !_hasTOF ) || ( tofcol->size() == 0 && dqdxcol->size() == 0 ) ){ + + debug() << "TPC and TOF PID information are not available, skip event " << _nEvt << endmsg; + + for ( auto pfo : *inpfocol ){ + auto outpfo = pfo.clone(); + outpfocol->push_back(outpfo); + } + + _nEvt++; + return StatusCode::SUCCESS; + } + + for ( auto pfo : *inpfocol ){ + auto outpfo = pfo.clone(); + outpfocol->push_back(outpfo); + FillTPCTOFPID(dqdxcol, tofcol, outpfo); + debug()<<" cyber pfo pid : " << outpfo.getType() << " cyber pfo charge : " << outpfo.getCharge() << " cyber pfo energy "<< outpfo.getEnergy() << endmsg; + } + + _nEvt++; + + return StatusCode::SUCCESS; +}// end execute + +//------------------------------------------------------------------------------ +StatusCode FinalPIDAlg::finalize(){ + debug() << "Finalizing..." << endmsg; + return StatusCode::SUCCESS; +} + + +/* +void FinalPIDAlg::FillTPCPID(const edm4hep::ReconstructedParticleCollection* pfocol, const edm4hep::RecDqdxCollection* dqdxcol, edm4hep::ParticleIDCollection* pidcol){ + for (auto pfo : *pfocol){ + for (auto trk : pfo.getTracks()){ + for (auto dqdx : *dqdxcol){ + if (dqdx.getTrack() == trk){ + std::array<double, 5> chi2s = dqdx.getHypotheses(); + int minchi2idx = std::distance(chi2s.begin(), std::min_element(chi2s.begin(), chi2s.end())); + auto pid = pidcol->create(); + pid.setPDG( pfo.getCharge() * PDGIDs.at(minchi2idx) ); + pid.setLikelihood( chi2s[minchi2idx] ); + //pfo.setParticleID(pid); + + } + } + } + } +} +*/ + + +void FinalPIDAlg::FillTPCTOFPID(const edm4hep::RecDqdxCollection* dqdxcol, const edm4hep::RecTofCollection* tofcol, edm4hep::MutableReconstructedParticle& pfo){ + + for (auto trk : pfo.getTracks()){ + + std::array<double, 5> chi2s; chi2s.fill(0); + + for (auto dqdx : *dqdxcol){ + if (dqdx.getTrack() == trk){ + for (int i = 0; i < 5; i++){ + chi2s[i] = dqdx.getHypotheses(i).chi2; + } + } + } + + for (auto tof : *tofcol){ + if (tof.getTrack() == trk){ + double toft = tof.getTime(); + std::array<float, 5> tofexpts = tof.getTimeExp(); + double tofexpterr = tof.getSigma(); + std::array<double, 5> tofchi2s; + for (int i = 0; i < 5; i++){ + tofchi2s[i] = std::pow( (tofexpts[i] - toft) / tofexpterr, 2); + chi2s[i] += tofchi2s[i]; + } + } + } + + int minchi2idx = std::distance(chi2s.begin(), std::min_element(chi2s.begin(), chi2s.end())); + pfo.setType( pfo.getCharge() * PDGIDs.at(minchi2idx) ); + + } +} diff --git a/Reconstruction/ParticleID/src/FinalPIDAlg.h b/Reconstruction/ParticleID/src/FinalPIDAlg.h new file mode 100644 index 0000000000000000000000000000000000000000..5c5591368c7104737f240b1ba87a7afef34ada95 --- /dev/null +++ b/Reconstruction/ParticleID/src/FinalPIDAlg.h @@ -0,0 +1,48 @@ +#ifndef FinalPIDAlg_h +#define FinalPIDAlg_h 1 + +#include "k4FWCore/DataHandle.h" +#include "GaudiKernel/Algorithm.h" + +#include "edm4hep/ReconstructedParticleCollection.h" +#include "edm4cepc/RecTof.h" +#include "edm4cepc/RecTofCollection.h" +#include "edm4hep/RecDqdx.h" +#include "edm4hep/RecDqdxCollection.h" +#include "edm4hep/ParticleIDCollection.h" + +class FinalPIDAlg : public Algorithm { + public: + // Constructor of this form must be provided + FinalPIDAlg( const std::string& name, ISvcLocator* pSvcLocator ); + + // Three mandatory member functions of any algorithm + StatusCode initialize() override; + StatusCode execute() override; + StatusCode finalize() override; + + private: + DataHandle<edm4hep::ReconstructedParticleCollection> m_inPFOCol{"CyberPFO", Gaudi::DataHandle::Reader, this}; + DataHandle<edm4hep::RecTofCollection> m_inTofCol{"RecTofCollection", Gaudi::DataHandle::Reader, this}; + DataHandle<edm4hep::RecDqdxCollection> m_inDqdxCol{"DndxTracks", Gaudi::DataHandle::Reader, this}; + //DataHandle<edm4hep::ParticleIDCollection> m_PIDCol{"finalPID", Gaudi::DataHandle::Writer, this}; + DataHandle<edm4hep::ReconstructedParticleCollection> m_outPFOCol{"CyberPFOPID", Gaudi::DataHandle::Writer, this}; + Gaudi::Property<std::string> m_method{this, "Method", "TPC+TOF"}; + + //void FillTPCPID(const edm4hep::ReconstructedParticleCollection* pfocol, const edm4hep::RecDqdxCollection* dqdxcol, edm4hep::ParticleIDCollection* pidcol); + void FillTPCTOFPID(const edm4hep::RecDqdxCollection* dqdxcol, const edm4hep::RecTofCollection* tofcol, edm4hep::MutableReconstructedParticle& pfo); + + int _nEvt; + bool _hasTPC; + bool _hasTOF; + + const std::map<int, int> PDGIDs = { + {0, -11}, + {1, -13}, + {2, 211}, + {3, 321}, + {4, 2212}, + }; + +}; +#endif diff --git a/Reconstruction/RecPFACyber/script/rec.py b/Reconstruction/RecPFACyber/script/rec.py index e61b08432e13c26924be70ca2031a7d36b671523..3e2257cc2507f4764fc21242a015b9783573e54f 100644 --- a/Reconstruction/RecPFACyber/script/rec.py +++ b/Reconstruction/RecPFACyber/script/rec.py @@ -47,7 +47,9 @@ inp.collections = [ # "HCALEndcapsParticleAssoCol", "MCParticle", "CompleteTracks", - "CompleteTracksParticleAssociation"] + "CompleteTracksParticleAssociation", + "RecTofCollection", + "DndxTracks"] ########################################## ######### Reconstruction ################ @@ -140,6 +142,9 @@ CyberPFAlg.AlgParValues = [ ["BarCol","Cluster1DCol","HalfClusterCol"],#1 ["EcalCluster", "SimpleHCALCluster", "outputPFO"], #16 ["1.26","4.", "1.", "4."] ]#17 +from Configurables import FinalPIDAlg +pid = FinalPIDAlg("FinalPIDAlg") +pid.OutputPFOName = "CyberPFOPID" ############################################################################## # POD I/O @@ -154,7 +159,7 @@ out.outputCommands = ["keep *"] from Configurables import ApplicationMgr ApplicationMgr( - TopAlg=[inp, CyberPFAlg, out ], + TopAlg=[inp, CyberPFAlg, pid, out ], EvtSel="NONE", EvtMax=10, ExtSvc=[podioevent, geomsvc], diff --git a/Reconstruction/RecPFACyber/script/tracking.py b/Reconstruction/RecPFACyber/script/tracking.py index 76c38cb87933e5a3062a7ff65593585214431583..dfb749c291d30fa6ffd1f4b83a021c78fa3cf018 100644 --- a/Reconstruction/RecPFACyber/script/tracking.py +++ b/Reconstruction/RecPFACyber/script/tracking.py @@ -246,6 +246,10 @@ from Configurables import TPCDndxAlg tpc_dndx = TPCDndxAlg("TPCDndxAlg") tpc_dndx.Method = "Simple" +from Configurables import TofRecAlg +tof = TofRecAlg("TofRecAlg") +#tof.OutputLevel = DEBUG + from Configurables import TrackParticleRelationAlg tpr = TrackParticleRelationAlg("Track2Particle") tpr.MCParticleCollection = "MCParticle" @@ -286,12 +290,14 @@ out.outputCommands = ["drop *", "keep HCALEndcapsParticleAssoCol", "keep MCParticle", "keep CompleteTracks", - "keep CompleteTracksParticleAssociation" ] + "keep CompleteTracksParticleAssociation", + "keep RecTofCollection", + "keep DndxTracks" ] # ApplicationMgr from Configurables import ApplicationMgr mgr = ApplicationMgr( - TopAlg = [podioinput, digiVXD, digiSIT, digiSET, digiFTD, digiTPC, digiMuon, tracking, forward, subset, clupatra, full, tpr, tpc_dndx, tmt, out], + TopAlg = [podioinput, digiVXD, digiSIT, digiSET, digiFTD, digiTPC, digiMuon, tracking, forward, subset, clupatra, full, tpr, tpc_dndx, tof, tmt, out], EvtSel = 'NONE', EvtMax = 10, ExtSvc = [rndmengine, rndmgensvc, dsvc, evtseeder, geosvc, gearsvc, tracksystemsvc, pidsvc],