Skip to content
Snippets Groups Projects
Commit fd33b147 authored by Wouter Deconinck's avatar Wouter Deconinck
Browse files

chore: black on ecal options before merge into vgawas branch

parent bd7b3826
No related branches found
No related tags found
1 merge request!268chore: black on ecal options before merge into vgawas branch
'''
"""
An example script to digitize/reconstruct/clustering endcap ecal hits
'''
"""
from Gaudi.Configuration import *
import json
import os
......@@ -13,23 +13,30 @@ detector_name = str(os.environ.get("DETECTOR", "epic"))
detector_config = str(os.environ.get("DETECTOR_CONFIG", detector_name))
detector_version = str(os.environ.get("DETECTOR_VERSION", "main"))
detector_path = str(os.environ.get("DETECTOR_PATH", "."))
compact_path = str(os.environ.get("JUGGLER_COMPACT_PATH", "{}.xml".format(os.path.join(detector_path, detector_name))))
compact_path = str(
os.environ.get(
"JUGGLER_COMPACT_PATH",
"{}.xml".format(os.path.join(detector_path, detector_name)),
)
)
# Detector features that affect reconstruction
has_ecal_barrel_scfi = False
if 'epic' in detector_name and 'imaging' in detector_config:
if "epic" in detector_name and "imaging" in detector_config:
has_ecal_barrel_scfi = True
# input and output
input_sims = [f.strip() for f in str.split(os.environ["JUGGLER_SIM_FILE"], ",") if f.strip()]
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"])
# input arguments from calibration file
with open(f'{detector_path}/calibrations/emcal_barrel_calibration.json') as f:
calib_data = json.load(f)['electron']
with open(f"{detector_path}/calibrations/emcal_barrel_calibration.json") as f:
calib_data = json.load(f)["electron"]
cb_ecal_sf = float(calib_data['sampling_fraction_img'])
cb_ecal_sf = float(calib_data["sampling_fraction_img"])
# geometry service
geo_service = GeoSvc("GeoSvc", detectors=[compact_path], OutputLevel=INFO)
......@@ -65,96 +72,110 @@ algorithms.append(podout)
if has_ecal_barrel_scfi:
# Central Barrel Ecal (Imaging Cal.)
cb_ecal_daq = dict(
dynamicRangeADC=3*MeV,
capacityADC=8192,
pedestalMean=400,
pedestalSigma=20) # about 6 keV
dynamicRangeADC=3 * MeV, capacityADC=8192, pedestalMean=400, pedestalSigma=20
) # about 6 keV
cb_ecal_digi = CalHitDigi("cb_ecal_digi",
cb_ecal_digi = CalHitDigi(
"cb_ecal_digi",
inputHitCollection="EcalBarrelHits",
outputHitCollection="EcalBarrelImagingHitsDigi",
energyResolutions=[0., 0.02, 0.], # 2% flat resolution
**cb_ecal_daq)
energyResolutions=[0.0, 0.02, 0.0], # 2% flat resolution
**cb_ecal_daq,
)
algorithms.append(cb_ecal_digi)
cb_ecal_reco = ImCalPixelReco("cb_ecal_reco",
cb_ecal_reco = ImCalPixelReco(
"cb_ecal_reco",
inputHitCollection=cb_ecal_digi.outputHitCollection,
outputHitCollection="EcalBarrelImagingHitsReco",
samplingFraction=cb_ecal_sf,
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)
layerField="layer", # field to get layer id
sectorField="module", # field to get sector id
**cb_ecal_daq,
)
algorithms.append(cb_ecal_reco)
cb_ecal_cl = ImagingCluster("cb_ecal_cl",
cb_ecal_cl = ImagingCluster(
"cb_ecal_cl",
inputHitCollection=cb_ecal_reco.outputHitCollection,
outputProtoClusterCollection="EcalBarrelImagingProtoClusters",
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
localDistXY=[2.0 * mm, 2 * mm], # same layer
layerDistEtaPhi=[10 * mrad, 10 * mrad], # adjacent layer
neighbourLayersRange=2, # id diff for adjacent layer
sectorDist=3.0 * cm,
) # different sector
algorithms.append(cb_ecal_cl)
cb_ecal_clreco = ImagingClusterReco("cb_ecal_clreco",
cb_ecal_clreco = ImagingClusterReco(
"cb_ecal_clreco",
inputProtoClusters=cb_ecal_cl.outputProtoClusterCollection,
outputClusters="EcalBarrelImagingClusters",
outputLayers="EcalBarrelImagingLayers",
mcHits="EcalBarrelHits")
mcHits="EcalBarrelHits",
)
algorithms.append(cb_ecal_clreco)
else:
# SciGlass calorimeter
sciglass_ecal_daq = dict(
dynamicRangeADC=5.*GeV,
capacityADC=32768,
pedestalMean=400,
pedestalSigma=3)
dynamicRangeADC=5.0 * GeV, capacityADC=32768, pedestalMean=400, pedestalSigma=3
)
sciglass_ecal_digi = CalHitDigi("sciglass_ecal_digi",
sciglass_ecal_digi = CalHitDigi(
"sciglass_ecal_digi",
inputHitCollection="EcalBarrelHits",
outputHitCollection="EcalBarrelHitsDigi",
energyResolutions=[0., 0.02, 0.], # 2% flat resolution
**sciglass_ecal_daq)
energyResolutions=[0.0, 0.02, 0.0], # 2% flat resolution
**sciglass_ecal_daq,
)
algorithms.append(sciglass_ecal_digi)
sciglass_ecal_reco = CalHitReco("sciglass_ecal_reco",
sciglass_ecal_reco = CalHitReco(
"sciglass_ecal_reco",
inputHitCollection=sciglass_ecal_digi.outputHitCollection,
outputHitCollection="EcalBarrelHitsReco",
thresholdFactor=3, # about 20 keV
readoutClass="EcalBarrelHits", # readout class
sectorField="sector", # field to get sector id
samplingFraction=0.998, # this accounts for a small fraction of leakage
**sciglass_ecal_daq)
sectorField="sector", # field to get sector id
samplingFraction=0.998, # this accounts for a small fraction of leakage
**sciglass_ecal_daq,
)
algorithms.append(sciglass_ecal_reco)
sciglass_ecal_cl = IslandCluster("sciglass_ecal_cl",
sciglass_ecal_cl = IslandCluster(
"sciglass_ecal_cl",
inputHitCollection=sciglass_ecal_reco.outputHitCollection,
outputProtoClusterCollection="EcalBarrelProtoClusters",
splitCluster=False,
minClusterHitEdep=1.0*MeV, # discard low energy hits
minClusterCenterEdep=30*MeV,
sectorDist=5.0*cm)
minClusterHitEdep=1.0 * MeV, # discard low energy hits
minClusterCenterEdep=30 * MeV,
sectorDist=5.0 * cm,
)
algorithms.append(sciglass_ecal_cl)
sciglass_ecal_clreco = ImagingClusterReco("sciglass_ecal_clreco",
sciglass_ecal_clreco = ImagingClusterReco(
"sciglass_ecal_clreco",
inputProtoClusters=sciglass_ecal_cl.outputProtoClusterCollection,
mcHits="EcalBarrelHits",
outputClusters="EcalBarrelClusters",
outputLayers="EcalBarrelLayers")
outputLayers="EcalBarrelLayers",
)
algorithms.append(sciglass_ecal_clreco)
podout.outputCommands = ['drop *',
'keep MCParticles',
'keep *HitsReco',
'keep *HitsDigi',
'keep *Cluster*']
podout.outputCommands = [
"drop *",
"keep MCParticles",
"keep *HitsReco",
"keep *HitsDigi",
"keep *Cluster*",
]
ApplicationMgr(
TopAlg = algorithms,
EvtSel = 'NONE',
EvtMax = n_events,
ExtSvc = [podioevent],
OutputLevel=DEBUG
TopAlg=algorithms,
EvtSel="NONE",
EvtMax=n_events,
ExtSvc=[podioevent],
OutputLevel=DEBUG,
)
'''
"""
An example script to digitize/reconstruct/clustering endcap ecal hits
'''
"""
from Gaudi.Configuration import *
import os
import ROOT
......@@ -10,10 +10,17 @@ from GaudiKernel.SystemOfUnits import MeV, GeV, mm, cm, mrad
detector_name = str(os.environ.get("DETECTOR_CONFIG", "epic"))
detector_path = str(os.environ.get("DETECTOR_PATH", "."))
compact_path = str(os.environ.get("JUGGLER_COMPACT_PATH", "{}.xml".format(os.path.join(detector_path, detector_name))))
compact_path = str(
os.environ.get(
"JUGGLER_COMPACT_PATH",
"{}.xml".format(os.path.join(detector_path, detector_name)),
)
)
# input and output
input_sims = [f.strip() for f in str.split(os.environ["JUGGLER_SIM_FILE"], ",") if f.strip()]
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"])
......@@ -43,53 +50,59 @@ podout = PodioOutput("out", filename=output_rec)
# 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=ce_ecal_digi.outputHitCollection,
outputHitCollection="EcalEndcapNHitsReco",
samplingFraction=0.998, # this accounts for a small fraction of leakage
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=ce_ecal_reco.outputHitCollection,
outputProtoClusterCollection="EcalEndcapNProtoClusters",
splitCluster=False,
minClusterHitEdep=1.0*MeV, # discard low energy hits
minClusterCenterEdep=30*MeV,
sectorDist=5.0*cm,
dimScaledLocalDistXY=[1.8, 1.8]) # hybrid calorimeter with different dimensions, using a scaled dist
ce_ecal_clreco = RecoCoG("ce_ecal_clreco",
inputProtoClusterCollection=ce_ecal_cl.outputProtoClusterCollection,
outputClusterCollection="EcalEndcapNClusters",
logWeightBase=4.6)
podout.outputCommands = ['drop *',
'keep MCParticles',
'keep *HitsReco',
'keep *HitsDigi',
'keep *Cluster*']
dynamicRangeADC=5.0 * GeV, capacityADC=32768, pedestalMean=400, pedestalSigma=3
)
ce_ecal_digi = CalHitDigi(
"ce_ecal_digi",
inputHitCollection="EcalEndcapNHits",
outputHitCollection="EcalEndcapNHitsDigi",
energyResolutions=[0.0, 0.02, 0.0],
**ce_ecal_daq
)
ce_ecal_reco = CalHitReco(
"ce_ecal_reco",
inputHitCollection=ce_ecal_digi.outputHitCollection,
outputHitCollection="EcalEndcapNHitsReco",
samplingFraction=0.998, # this accounts for a small fraction of leakage
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=ce_ecal_reco.outputHitCollection,
outputProtoClusterCollection="EcalEndcapNProtoClusters",
splitCluster=False,
minClusterHitEdep=1.0 * MeV, # discard low energy hits
minClusterCenterEdep=30 * MeV,
sectorDist=5.0 * cm,
dimScaledLocalDistXY=[1.8, 1.8],
) # hybrid calorimeter with different dimensions, using a scaled dist
ce_ecal_clreco = RecoCoG(
"ce_ecal_clreco",
inputProtoClusterCollection=ce_ecal_cl.outputProtoClusterCollection,
outputClusterCollection="EcalEndcapNClusters",
logWeightBase=4.6,
)
podout.outputCommands = [
"drop *",
"keep MCParticles",
"keep *HitsReco",
"keep *HitsDigi",
"keep *Cluster*",
]
ApplicationMgr(
TopAlg = [podin,
ce_ecal_digi, ce_ecal_reco, ce_ecal_cl, ce_ecal_clreco,
podout],
EvtSel = 'NONE',
EvtMax = n_events,
ExtSvc = [podioevent],
OutputLevel=DEBUG
TopAlg=[podin, ce_ecal_digi, ce_ecal_reco, ce_ecal_cl, ce_ecal_clreco, podout],
EvtSel="NONE",
EvtMax=n_events,
ExtSvc=[podioevent],
OutputLevel=DEBUG,
)
'''
"""
An example script to digitize/reconstruct/clustering endcap ecal hits
'''
"""
from Gaudi.Configuration import *
import os
import ROOT
......@@ -11,11 +10,18 @@ from GaudiKernel.SystemOfUnits import MeV, GeV, mm, cm, mrad
detector_name = str(os.environ.get("DETECTOR_CONFIG", "epic"))
detector_path = str(os.environ.get("DETECTOR_PATH", "."))
compact_path = str(os.environ.get("JUGGLER_COMPACT_PATH", "{}.xml".format(os.path.join(detector_path, detector_name))))
compact_path = str(
os.environ.get(
"JUGGLER_COMPACT_PATH",
"{}.xml".format(os.path.join(detector_path, detector_name)),
)
)
ci_ecal_sf = float(os.environ.get("CI_ECAL_SAMP_FRAC", 0.253))
# input and output
input_sims = [f.strip() for f in str.split(os.environ["JUGGLER_SIM_FILE"], ",") if f.strip()]
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"])
......@@ -46,59 +52,75 @@ podout = PodioOutput("out", filename=output_rec)
# 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=ci_ecal_digi.outputHitCollection,
outputHitCollection="EcalEndcapPHitsReco",
samplingFraction=ci_ecal_sf,
thresholdFactor=5.0,
**ci_ecal_daq)
dynamicRangeADC=50.0 * 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=ci_ecal_digi.outputHitCollection,
outputHitCollection="EcalEndcapPHitsReco",
samplingFraction=ci_ecal_sf,
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=ci_ecal_reco.outputHitCollection,
outputHitCollection="EcalEndcapPHitsRecoXY",
fields=["fiber_x", "fiber_y"],
fieldRefNumbers=[1, 1],
# fields=["layer", "slice"],
# fieldRefNumbers=[1, 0],
readoutClass="EcalEndcapPHits")
ci_ecal_cl = IslandCluster("ci_ecal_cl",
# OutputLevel=DEBUG,
inputHitCollection=ci_ecal_merger.outputHitCollection,
outputProtoClusterCollection="EcalEndcapPProtoClusters",
splitCluster=False,
minClusterCenterEdep=10.*MeV,
localDistXY=[10*mm, 10*mm])
ci_ecal_clreco = RecoCoG("ci_ecal_clreco",
inputProtoClusterCollection=ci_ecal_cl.outputProtoClusterCollection,
outputClusterCollection="EcalEndcapPClusters",
logWeightBase=6.2)
podout.outputCommands = ['drop *',
'keep MCParticles',
'keep *HitsReco*',
'keep *HitsDigi',
'keep *Cluster*']
ci_ecal_merger = CalHitsMerger(
"ci_ecal_merger",
# OutputLevel=DEBUG,
inputHitCollection=ci_ecal_reco.outputHitCollection,
outputHitCollection="EcalEndcapPHitsRecoXY",
fields=["fiber_x", "fiber_y"],
fieldRefNumbers=[1, 1],
# fields=["layer", "slice"],
# fieldRefNumbers=[1, 0],
readoutClass="EcalEndcapPHits",
)
ci_ecal_cl = IslandCluster(
"ci_ecal_cl",
# OutputLevel=DEBUG,
inputHitCollection=ci_ecal_merger.outputHitCollection,
outputProtoClusterCollection="EcalEndcapPProtoClusters",
splitCluster=False,
minClusterCenterEdep=10.0 * MeV,
localDistXY=[10 * mm, 10 * mm],
)
ci_ecal_clreco = RecoCoG(
"ci_ecal_clreco",
inputProtoClusterCollection=ci_ecal_cl.outputProtoClusterCollection,
outputClusterCollection="EcalEndcapPClusters",
logWeightBase=6.2,
)
podout.outputCommands = [
"drop *",
"keep MCParticles",
"keep *HitsReco*",
"keep *HitsDigi",
"keep *Cluster*",
]
ApplicationMgr(
TopAlg = [podin,
ci_ecal_digi, ci_ecal_reco, ci_ecal_merger, ci_ecal_cl, ci_ecal_clreco,
podout],
EvtSel = 'NONE',
EvtMax = n_events,
ExtSvc = [podioevent],
OutputLevel=DEBUG
TopAlg=[
podin,
ci_ecal_digi,
ci_ecal_reco,
ci_ecal_merger,
ci_ecal_cl,
ci_ecal_clreco,
podout,
],
EvtSel="NONE",
EvtMax=n_events,
ExtSvc=[podioevent],
OutputLevel=DEBUG,
)
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment