Skip to content
Snippets Groups Projects

adding full_trackpluscalo option

Merged adding full_trackpluscalo option
trackpluscalo into master
Merged Miguel Arratia requested to merge trackpluscalo into master
Compare
and
1 file
+ 502
− 0
Compare changes
  • Side-by-side
  • Inline
benchmarks/clustering/options/full_trackpluscalo_reco.py 0 → 100644
+ 502
− 0
 
'''
 
An example option file to digitize/reconstruct/clustering calorimeter hits
 
'''
 
from Gaudi.Configuration import *
 
import os
 
import ROOT
 
 
from GaudiKernel.DataObjectHandleBase import DataObjectHandleBase
 
from Configurables import ApplicationMgr, EICDataSvc, PodioInput, PodioOutput, GeoSvc
 
from GaudiKernel.SystemOfUnits import MeV, GeV, mm, cm, mrad
 
 
detector_name = str(os.environ.get("JUGGLER_DETECTOR", "athena"))
 
detector_path = str(os.environ.get("JUGGLER_DETECTOR_PATH", "."))
 
compact_path = os.path.join(detector_path, detector_name)
 
 
# get sampling fractions from system environment variable, 1.0 by default
 
ci_ecal_sf = float(os.environ.get("CI_ECAL_SAMP_FRAC", 0.253))
 
cb_ecal_sf = float(os.environ.get("CB_ECAL_SAMP_FRAC", 0.01324))
 
cb_hcal_sf = float(os.environ.get("CB_HCAL_SAMP_FRAC", 0.038))
 
ci_hcal_sf = float(os.environ.get("CI_HCAL_SAMP_FRAC", 0.025))
 
ce_hcal_sf = float(os.environ.get("CE_HCAL_SAMP_FRAC", 0.025))
 
scifi_barrel_sf = float(os.environ.get("CB_EMCAL_SCFI_SAMP_FRAC",0.0938))
 
# input and output
 
input_sims = [f.strip() for f in str.split(os.environ["JUGGLER_SIM_FILE"], ",") if f.strip()]
 
output_rec = str(os.environ["JUGGLER_REC_FILE"])
 
n_events = int(os.environ["JUGGLER_N_EVENTS"])
 
 
# geometry service
 
geo_service = GeoSvc("GeoSvc", detectors=["{}.xml".format(compact_path)], OutputLevel=INFO)
 
# data service
 
podioevent = EICDataSvc("EventDataSvc", inputs=input_sims)
 
 
 
# juggler components
 
from Configurables import Jug__Base__InputCopier_dd4pod__Geant4ParticleCollection_dd4pod__Geant4ParticleCollection_ as MCCopier
 
from Configurables import Jug__Base__InputCopier_dd4pod__CalorimeterHitCollection_dd4pod__CalorimeterHitCollection_ as CalCopier
 
from Configurables import Jug__Base__InputCopier_dd4pod__TrackerHitCollection_dd4pod__TrackerHitCollection_ as TrkCopier
 
 
from Configurables import Jug__Digi__UFSDTrackerDigi as TrackerDigi
 
from Configurables import Jug__Reco__TrackerHitReconstruction as TrackerHitReconstruction
 
from Configurables import Jug__Reco__TrackerSourceLinker as TrackerSourceLinker
 
from Configurables import Jug__Reco__TrackerSourcesLinker as TrackerSourcesLinker
 
#from Configurables import Jug__Reco__TrackingHitsSourceLinker as TrackingHitsSourceLinker
 
from Configurables import Jug__Reco__TrackParamTruthInit as TrackParamTruthInit
 
from Configurables import Jug__Reco__TrackParamClusterInit as TrackParamClusterInit
 
from Configurables import Jug__Reco__TrackParamVertexClusterInit as TrackParamVertexClusterInit
 
 
from Configurables import Jug__Reco__TrackFindingAlgorithm as TrackFindingAlgorithm
 
from Configurables import Jug__Reco__ParticlesFromTrackFit as ParticlesFromTrackFit
 
 
 
 
 
 
from Configurables import Jug__Digi__CalorimeterHitDigi as CalHitDigi
 
from Configurables import Jug__Reco__CalorimeterHitReco as CalHitReco
 
from Configurables import Jug__Reco__CalorimeterHitsMerger as CalHitsMerger
 
from Configurables import Jug__Reco__CalorimeterIslandCluster as IslandCluster
 
 
from Configurables import Jug__Reco__ImagingPixelReco as ImCalPixelReco
 
from Configurables import Jug__Reco__ImagingTopoCluster as ImagingCluster
 
 
from Configurables import Jug__Reco__ClusterRecoCoG as RecoCoG
 
from Configurables import Jug__Reco__ImagingClusterReco as ImagingClusterReco
 
 
# branches needed from simulation root file
 
sim_coll = [
 
"mcparticles",
 
"EcalEndcapNHits",
 
"EcalEndcapPHits",
 
"EcalBarrelHits",
 
"HcalBarrelHits",
 
"EcalBarrelScFiHits",
 
"HcalHadronEndcapHits",
 
"HcalElectronEndcapHits",
 
"TrackerEndcapHits",
 
"TrackerBarrelHits",
 
"VertexBarrelHits",
 
"VertexEndcapHits"
 
]
 
 
# input and output
 
podin = PodioInput("PodioReader", collections=sim_coll)
 
podout = PodioOutput("out", filename=output_rec)
 
# copier needed to get around input --> output copy bug. So truth (mcparticles) can be saved in output file
 
copier = MCCopier("MCCopier", inputCollection="mcparticles", outputCollection="mcparticles2")
 
trkcopier = TrkCopier("TrkCopier",
 
inputCollection="TrackerBarrelHits",
 
outputCollection="TrackerBarrelHits2")
 
 
 
 
##Tracking digi
 
trk_b_digi = TrackerDigi("trk_b_digi",
 
inputHitCollection="TrackerBarrelHits",
 
outputHitCollection="TrackerBarrelRawHits",
 
timeResolution=8)
 
trk_ec_digi = TrackerDigi("trk_ec_digi",
 
inputHitCollection="TrackerEndcapHits",
 
outputHitCollection="TrackerEndcapRawHits",
 
timeResolution=8)
 
 
vtx_b_digi = TrackerDigi("vtx_b_digi",
 
inputHitCollection="VertexBarrelHits",
 
outputHitCollection="VertexBarrelRawHits",
 
timeResolution=8)
 
 
vtx_ec_digi = TrackerDigi("vtx_ec_digi",
 
inputHitCollection="VertexEndcapHits",
 
outputHitCollection="VertexEndcapRawHits",
 
timeResolution=8)
 
 
# Tracker and vertex reconstruction
 
trk_b_reco = TrackerHitReconstruction("trk_b_reco",
 
inputHitCollection = trk_b_digi.outputHitCollection,
 
outputHitCollection="TrackerBarrelRecHits")
 
 
trk_ec_reco = TrackerHitReconstruction("trk_ec_reco",
 
inputHitCollection = trk_ec_digi.outputHitCollection,
 
outputHitCollection="TrackerEndcapRecHits")
 
 
vtx_b_reco = TrackerHitReconstruction("vtx_b_reco",
 
inputHitCollection = vtx_b_digi.outputHitCollection,
 
outputHitCollection="VertexBarrelRecHits")
 
 
vtx_ec_reco = TrackerHitReconstruction("vtx_ec_reco",
 
inputHitCollection = vtx_ec_digi.outputHitCollection,
 
outputHitCollection="VertexEndcapRecHits")
 
 
# Hit Source linker
 
sourcelinker = TrackerSourceLinker("sourcelinker",
 
inputHitCollection=trk_b_reco.outputHitCollection,
 
outputSourceLinks="BarrelTrackSourceLinks",
 
outputMeasurements="BarrelTrackMeasurements",
 
OutputLevel=DEBUG)
 
 
#trk_hits_srclnkr = TrackerSourcesLinker("trk_srcslnkr",
 
# ITrackerBarrelHits = vtx_b_reco.outputHitCollection,
 
# ITrackerEndcapHits = vtx_ec_reco.outputHitCollection,
 
# OTrackerBarrelHits = trk_b_reco.outputHitCollection,
 
# OTrackerEndcapHits = trk_ec_reco.outputHitCollection,
 
# outputSourceLinks="TrackerMeasurements",
 
# OutputLevel=DEBUG)
 
 
## Track param init
 
truth_trk_init = TrackParamTruthInit("truth_trk_init",
 
inputMCParticles="mcparticles",
 
outputInitialTrackParameters="InitTrackParams",
 
OutputLevel=DEBUG)
 
 
#clust_trk_init = TrackParamClusterInit("clust_trk_init",
 
# inputClusters="SimpleClusters",
 
# outputInitialTrackParameters="InitTrackParamsFromClusters",
 
# OutputLevel=DEBUG)
 
 
#vtxcluster_trk_init = TrackParamVertexClusterInit("vtxcluster_trk_init",
 
# inputVertexHits="VertexBarrelRecHits",
 
# inputClusters="SimpleClusters",
 
# outputInitialTrackParameters="InitTrackParamsFromVtxClusters",
 
# maxHitRadius=40.0*units.mm,
 
# OutputLevel=WARNING)
 
 
# Tracking algorithms
 
trk_find_alg = TrackFindingAlgorithm("trk_find_alg",
 
inputSourceLinks = sourcelinker.outputSourceLinks,
 
inputMeasurements = sourcelinker.outputMeasurements,
 
inputInitialTrackParameters= "InitTrackParams",#"InitTrackParamsFromClusters",
 
outputTrajectories="trajectories",
 
OutputLevel=DEBUG)
 
 
parts_from_fit = ParticlesFromTrackFit("parts_from_fit",
 
inputTrajectories="trajectories",
 
outputParticles="ReconstructedParticles",
 
outputTrackParameters="outputTrackParameters",
 
OutputLevel=DEBUG)
 
 
#trk_find_alg1 = TrackFindingAlgorithm("trk_find_alg1",
 
# inputSourceLinks = sourcelinker.outputSourceLinks,
 
# inputMeasurements = sourcelinker.outputMeasurements,
 
# inputInitialTrackParameters= "InitTrackParamsFromClusters",
 
# outputTrajectories="trajectories1",
 
# OutputLevel=DEBUG)
 
 
#parts_from_fit1 = ParticlesFromTrackFit("parts_from_fit1",
 
# inputTrajectories="trajectories1",
 
# outputParticles="ReconstructedParticles1",
 
# outputTrackParameters="outputTrackParameters1",
 
# OutputLevel=DEBUG)
 
 
 
 
# Crystal Endcap Ecal
 
ce_ecal_daq = dict(
 
dynamicRangeADC=5.*GeV,
 
capacityADC=32768,
 
pedestalMean=400,
 
pedestalSigma=3)
 
 
ce_ecal_digi = CalHitDigi("ce_ecal_digi",
 
inputHitCollection="EcalEndcapNHits",
 
outputHitCollection="EcalEndcapNHitsDigi",
 
energyResolutions=[0., 0.02, 0.],
 
**ce_ecal_daq)
 
 
ce_ecal_reco = CalHitReco("ce_ecal_reco",
 
inputHitCollection="EcalEndcapNHitsDigi",
 
outputHitCollection="EcalEndcapNHitsReco",
 
thresholdFactor=4, # 4 sigma cut on pedestal sigma
 
readoutClass="EcalEndcapNHits",
 
sectorField="sector",
 
**ce_ecal_daq)
 
 
ce_ecal_cl = IslandCluster("ce_ecal_cl",
 
# OutputLevel=DEBUG,
 
inputHitCollection="EcalEndcapNHitsReco",
 
outputHitCollection="EcalEndcapNClusterHits",
 
splitCluster=False,
 
minClusterHitEdep=1.0*MeV, # discard low energy hits
 
minClusterCenterEdep=30*MeV,
 
sectorDist=5.0*cm,
 
dimScaledLocalDistXY=[1.8, 1.8]) # dimension scaled dist is good for hybrid sectors with different module size
 
 
ce_ecal_clreco = RecoCoG("ce_ecal_clreco",
 
inputHitCollection="EcalEndcapNClusterHits",
 
outputClusterCollection="EcalEndcapNClusters",
 
samplingFraction=0.998, # this accounts for a small fraction of leakage
 
logWeightBase=4.6)
 
 
# Endcap Sampling Ecal
 
ci_ecal_daq = dict(
 
dynamicRangeADC=50.*MeV,
 
capacityADC=32768,
 
pedestalMean=400,
 
pedestalSigma=10)
 
 
ci_ecal_digi = CalHitDigi("ci_ecal_digi",
 
inputHitCollection="EcalEndcapPHits",
 
outputHitCollection="EcalEndcapPHitsDigi",
 
**ci_ecal_daq)
 
 
ci_ecal_reco = CalHitReco("ci_ecal_reco",
 
inputHitCollection="EcalEndcapPHitsDigi",
 
outputHitCollection="EcalEndcapPHitsReco",
 
thresholdFactor=5.0,
 
**ci_ecal_daq)
 
 
# merge hits in different layer (projection to local x-y plane)
 
ci_ecal_merger = CalHitsMerger("ci_ecal_merger",
 
# OutputLevel=DEBUG,
 
inputHitCollection="EcalEndcapPHitsReco",
 
outputHitCollection="EcalEndcapPHitsRecoXY",
 
fields=["layer", "slice"],
 
fieldRefNumbers=[1, 0],
 
readoutClass="EcalEndcapPHits")
 
 
ci_ecal_cl = IslandCluster("ci_ecal_cl",
 
# OutputLevel=DEBUG,
 
inputHitCollection="EcalEndcapPHitsRecoXY",
 
outputHitCollection="EcalEndcapPClusterHits",
 
splitCluster=False,
 
minClusterCenterEdep=10.*MeV,
 
localDistXY=[10*mm, 10*mm])
 
 
ci_ecal_clreco = RecoCoG("ci_ecal_clreco",
 
inputHitCollection="EcalEndcapPClusterHits",
 
outputClusterCollection="EcalEndcapPClusters",
 
logWeightBase=6.2,
 
samplingFraction=ci_ecal_sf)
 
 
# Central Barrel Ecal (Imaging Cal.)
 
cb_ecal_daq = dict(
 
dynamicRangeADC=3*MeV,
 
capacityADC=8192,
 
pedestalMean=400,
 
pedestalSigma=20) # about 6 keV
 
 
cb_ecal_digi = CalHitDigi("cb_ecal_digi",
 
inputHitCollection="EcalBarrelHits",
 
outputHitCollection="EcalBarrelHitsDigi",
 
energyResolutions=[0., 0.02, 0.], # 2% flat resolution
 
**cb_ecal_daq)
 
 
cb_ecal_reco = ImCalPixelReco("cb_ecal_reco",
 
inputHitCollection="EcalBarrelHitsDigi",
 
outputHitCollection="EcalBarrelHitsReco",
 
thresholdFactor=3, # about 20 keV
 
readoutClass="EcalBarrelHits", # readout class
 
layerField="layer", # field to get layer id
 
sectorField="module", # field to get sector id
 
**cb_ecal_daq)
 
 
cb_ecal_cl = ImagingCluster("cb_ecal_cl",
 
inputHitCollection="EcalBarrelHitsReco",
 
outputHitCollection="EcalBarrelClusterHits",
 
localDistXY=[2.*mm, 2*mm], # same layer
 
layerDistEtaPhi=[10*mrad, 10*mrad], # adjacent layer
 
neighbourLayersRange=2, # id diff for adjacent layer
 
sectorDist=3.*cm) # different sector
 
 
cb_ecal_clreco = ImagingClusterReco("cb_ecal_clreco",
 
samplingFraction=cb_ecal_sf,
 
inputHitCollection="EcalBarrelClusterHits",
 
outputClusterCollection="EcalBarrelClusters",
 
outputLayerCollection="EcalBarrelLayers")
 
 
#Central ECAL SciFi
 
# use the same daq_setting for digi/reco pair
 
scfi_barrel_daq = dict(
 
dynamicRangeADC=50.*MeV,
 
capacityADC=32768,
 
pedestalMean=400,
 
pedestalSigma=10)
 
 
scfi_barrel_digi = CalHitDigi("scfi_barrel_digi",
 
inputHitCollection="EcalBarrelScFiHits",
 
outputHitCollection="EcalBarrelScFiHitsDigi",
 
**scfi_barrel_daq)
 
 
scfi_barrel_reco = CalHitReco("scfi_barrel_reco",
 
inputHitCollection="EcalBarrelScFiHitsDigi",
 
outputHitCollection="EcalBarrelScFiHitsReco",
 
thresholdFactor=5.0,
 
readoutClass="EcalBarrelScFiHits",
 
layerField="layer",
 
sectorField="module",
 
localDetFields=["system", "module"], # use local coordinates in each module (stave)
 
**scfi_barrel_daq)
 
 
# merge hits in different layer (projection to local x-y plane)
 
scfi_barrel_merger = CalHitsMerger("scfi_barrel_merger",
 
# OutputLevel=DEBUG,
 
inputHitCollection="EcalBarrelScFiHitsReco",
 
outputHitCollection="EcalBarrelScFiGridReco",
 
fields=["fiber"],
 
fieldRefNumbers=[1],
 
readoutClass="EcalBarrelScFiHits")
 
 
scfi_barrel_cl = IslandCluster("scfi_barrel_cl",
 
# OutputLevel=DEBUG,
 
inputHitCollection="EcalBarrelScFiGridReco",
 
outputHitCollection="EcalBarrelScFiClusterHits",
 
splitCluster=False,
 
minClusterCenterEdep=10.*MeV,
 
localDistXZ=[30*mm, 30*mm])
 
 
scfi_barrel_clreco = RecoCoG("scfi_barrel_clreco",
 
inputHitCollection="EcalBarrelScFiClusterHits",
 
outputClusterCollection="EcalBarrelScFiClusters",
 
logWeightBase=6.2,
 
samplingFraction= scifi_barrel_sf)
 
 
 
# Central Barrel Hcal
 
cb_hcal_daq = dict(
 
dynamicRangeADC=50.*MeV,
 
capacityADC=32768,
 
pedestalMean=400,
 
pedestalSigma=10)
 
 
cb_hcal_digi = CalHitDigi("cb_hcal_digi",
 
inputHitCollection="HcalBarrelHits",
 
outputHitCollection="HcalBarrelHitsDigi",
 
**cb_hcal_daq)
 
 
cb_hcal_reco = CalHitReco("cb_hcal_reco",
 
inputHitCollection="HcalBarrelHitsDigi",
 
outputHitCollection="HcalBarrelHitsReco",
 
thresholdFactor=5.0,
 
readoutClass="HcalBarrelHits",
 
layerField="layer",
 
sectorField="module",
 
**cb_hcal_daq)
 
 
cb_hcal_merger = CalHitsMerger("cb_hcal_merger",
 
inputHitCollection="HcalBarrelHitsReco",
 
outputHitCollection="HcalBarrelHitsRecoXY",
 
readoutClass="HcalBarrelHits",
 
fields=["layer", "slice"],
 
fieldRefNumbers=[1, 0])
 
 
cb_hcal_cl = IslandCluster("cb_hcal_cl",
 
inputHitCollection="HcalBarrelHitsRecoXY",
 
outputHitCollection="HcalBarrelClusterHits",
 
splitCluster=False,
 
minClusterCenterEdep=30.*MeV,
 
localDistXY=[15.*cm, 15.*cm])
 
 
cb_hcal_clreco = RecoCoG("cb_hcal_clreco",
 
inputHitCollection="HcalBarrelClusterHits",
 
outputClusterCollection="HcalBarrelClusters",
 
logWeightBase=6.2,
 
samplingFraction=cb_hcal_sf)
 
 
 
# Hcal Hadron Endcap
 
ci_hcal_daq = dict(
 
dynamicRangeADC=50.*MeV,
 
capacityADC=32768,
 
pedestalMean=400,
 
pedestalSigma=10)
 
ci_hcal_digi = CalHitDigi("ci_hcal_digi",
 
inputHitCollection="HcalHadronEndcapHits",
 
outputHitCollection="HcalHadronEndcapHitsDigi",
 
**ci_hcal_daq)
 
 
ci_hcal_reco = CalHitReco("ci_hcal_reco",
 
inputHitCollection="HcalHadronEndcapHitsDigi",
 
outputHitCollection="HcalHadronEndcapHitsReco",
 
thresholdFactor=5.0,
 
**ci_hcal_daq)
 
 
ci_hcal_merger = CalHitsMerger("ci_hcal_merger",
 
inputHitCollection="HcalHadronEndcapHitsReco",
 
outputHitCollection="HcalHadronEndcapHitsRecoXY",
 
readoutClass="HcalHadronEndcapHits",
 
fields=["layer", "slice"],
 
fieldRefNumbers=[1, 0])
 
 
ci_hcal_cl = IslandCluster("ci_hcal_cl",
 
inputHitCollection="HcalHadronEndcapHitsRecoXY",
 
outputHitCollection="HcalHadronEndcapClusterHits",
 
splitCluster=False,
 
minClusterCenterEdep=30.*MeV,
 
localDistXY=[15.*cm, 15.*cm])
 
 
ci_hcal_clreco = RecoCoG("ci_hcal_clreco",
 
inputHitCollection="HcalHadronEndcapClusterHits",
 
outputClusterCollection="HcalHadronEndcapClusters",
 
logWeightBase=6.2,
 
samplingFraction=ci_hcal_sf)
 
 
# Hcal Electron Endcap
 
ce_hcal_daq = dict(
 
dynamicRangeADC=50.*MeV,
 
capacityADC=32768,
 
pedestalMean=400,
 
pedestalSigma=10)
 
 
ce_hcal_digi = CalHitDigi("ce_hcal_digi",
 
inputHitCollection="HcalElectronEndcapHits",
 
outputHitCollection="HcalElectronEndcapHitsDigi",
 
**ce_hcal_daq)
 
 
ce_hcal_reco = CalHitReco("ce_hcal_reco",
 
inputHitCollection="HcalElectronEndcapHitsDigi",
 
outputHitCollection="HcalElectronEndcapHitsReco",
 
thresholdFactor=5.0,
 
**ce_hcal_daq)
 
 
ce_hcal_merger = CalHitsMerger("ce_hcal_merger",
 
inputHitCollection="HcalElectronEndcapHitsReco",
 
outputHitCollection="HcalElectronEndcapHitsRecoXY",
 
readoutClass="HcalElectronEndcapHits",
 
fields=["layer", "slice"],
 
fieldRefNumbers=[1, 0])
 
 
ce_hcal_cl = IslandCluster("ce_hcal_cl",
 
inputHitCollection="HcalElectronEndcapHitsRecoXY",
 
outputHitCollection="HcalElectronEndcapClusterHits",
 
splitCluster=False,
 
minClusterCenterEdep=30.*MeV,
 
localDistXY=[15.*cm, 15.*cm])
 
 
ce_hcal_clreco = RecoCoG("ce_hcal_clreco",
 
inputHitCollection="HcalElectronEndcapClusterHits",
 
outputClusterCollection="HcalElectronEndcapClusters",
 
logWeightBase=6.2,
 
samplingFraction=ce_hcal_sf)
 
 
 
podout.outputCommands = ['drop *',
 
'keep mcparticles2',
 
'keep *Digi',
 
'keep *Reco*',
 
'keep *ClusterHits',
 
'keep *Clusters',
 
'keep *Layers',
 
'keep *outputTrackParameters*',
 
'keep *Track*']
 
 
ApplicationMgr(
 
TopAlg = [podin, copier,trkcopier,
 
trk_b_digi, trk_ec_digi, vtx_b_digi, vtx_ec_digi,
 
trk_b_reco, trk_ec_reco, vtx_b_reco, vtx_ec_reco,
 
sourcelinker,
 
#clust_trk_init,
 
truth_trk_init,
 
trk_find_alg, parts_from_fit,
 
#trk_find_alg1, parts_from_fit1,
 
ce_ecal_digi, ce_ecal_reco, ce_ecal_cl, ce_ecal_clreco,
 
ci_ecal_digi, ci_ecal_reco, ci_ecal_merger, ci_ecal_cl, ci_ecal_clreco,
 
cb_ecal_digi, cb_ecal_reco, cb_ecal_cl, cb_ecal_clreco,
 
scfi_barrel_digi, scfi_barrel_reco, scfi_barrel_merger, scfi_barrel_cl, scfi_barrel_clreco,
 
cb_hcal_digi, cb_hcal_reco, cb_hcal_merger, cb_hcal_cl, cb_hcal_clreco,
 
ce_hcal_digi, ce_hcal_reco, ce_hcal_merger, ce_hcal_cl, ce_hcal_clreco,
 
ci_hcal_digi, ci_hcal_reco, ci_hcal_merger, ci_hcal_cl, ci_hcal_clreco,
 
podout],
 
EvtSel = 'NONE',
 
EvtMax = n_events,
 
ExtSvc = [podioevent],
 
OutputLevel=DEBUG
 
)