scfi_cluster.py

import json
import os
import ROOT
from Gaudi.Configuration import *
from GaudiKernel.SystemOfUnits import mm, MeV, rad, ns
from Configurables import ApplicationMgr, EICDataSvc, PodioOutput, GeoSvc

from Configurables import PodioInput
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__ClusterRecoCoG as RecoCoG

detector_name = str(os.environ.get("DETECTOR", "athena"))
detector_config = str(os.environ.get("DETECTOR_CONFIG", detector_name))
detector_version = str(os.environ.get("DETECTOR_VERSION", "master"))

detector_path = str(os.environ.get("DETECTOR_PATH", "."))

# Detector features that affect reconstruction
has_ecal_barrel_scfi = False
if 'athena' in detector_name:
    has_ecal_barrel_scfi = True
if 'ecce' in detector_name and 'imaging' in detector_config:
    has_ecal_barrel_scfi = True
if 'epic' in detector_name and 'imaging' in detector_config:
    has_ecal_barrel_scfi = True

# input arguments from calibration file
with open(f'{detector_path}/calibrations/emcal_barrel_calibration.json') as f:
    calib_data = json.load(f)['electron']

kwargs = dict()
kwargs['sf'] = float(calib_data['sampling_fraction_scfi'])

# input arguments through environment variables
kwargs['input'] = os.environ.get('CB_EMCAL_SIM_FILE', '../topside/barrel_pion0_5GeV.edm4hep.root')
kwargs['output'] = os.environ.get('CB_EMCAL_REC_FILE', 'barrel_pion0_5GeV_cluster.root')
kwargs['compact'] = os.environ.get('CB_EMCAL_COMPACT_PATH', '../topside/test.xml')
kwargs['nev'] = int(os.environ.get('CB_EMCAL_NUMEV', 100))

if kwargs['nev'] < 1:
    f = ROOT.TFile(kwargs['input'])
    kwargs['nev'] = f.events.GetEntries()

print(kwargs)

# get sampling fraction from system environment variable, 1.0 by default
sf = float(os.environ.get('CB_EMCAL_SAMP_FRAC', '1.0'))
geo_service = GeoSvc("GeoSvc", detectors=kwargs['compact'].split(','), OutputLevel=INFO)
podioevent = EICDataSvc("EventDataSvc", inputs=kwargs['input'].split(','), OutputLevel=DEBUG)

sim_coll = [
    "MCParticles",
]
if has_ecal_barrel_scfi:
    sim_coll += [
      "EcalBarrelHits",
      "EcalBarrelHitsContributions",
    ]

algorithms = []

podin = PodioInput("PodioReader", collections=sim_coll, OutputLevel=DEBUG)
algorithms.append(podin)

podout = PodioOutput("out", filename=kwargs['output'])
algorithms.append(podout)

if has_ecal_barrel_scfi:
    # 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)
    algorithms.append(scfi_barrel_digi)

    scfi_barrel_reco = CalHitReco("scfi_barrel_reco",
        inputHitCollection=scfi_barrel_digi.outputHitCollection,
        outputHitCollection="EcalBarrelScFiHitsReco",
        thresholdFactor=5.0,
        readoutClass="EcalBarrelScFiHits",
        layerField="layer",
        sectorField="module",
        localDetFields=["system", "module"], # use local coordinates in each module (stave)
        **scfi_barrel_daq)
    algorithms.append(scfi_barrel_reco)

    # merge hits in different layer (projection to local x-y plane)
    scfi_barrel_merger = CalHitsMerger("scfi_barrel_merger",
        # OutputLevel=DEBUG,
        inputHitCollection=scfi_barrel_reco.outputHitCollection,
        outputHitCollection="EcalBarrelScFiGridReco",
        fields=["fiber"],
        fieldRefNumbers=[1],
        readoutClass="EcalBarrelScFiHits")
    algorithms.append(scfi_barrel_merger)

    scfi_barrel_cl = IslandCluster("scfi_barrel_cl",
        # OutputLevel=DEBUG,
        inputHitCollection=scfi_barrel_reco.outputHitCollection,
        outputProtoClusterCollection="EcalBarrelScFiProtoClusters",
        splitCluster=False,
        minClusterCenterEdep=10.*MeV,
        localDistXZ=[30*mm, 30*mm])
    algorithms.append(scfi_barrel_cl)

    scfi_barrel_clreco = RecoCoG("scfi_barrel_clreco",
       inputProtoClusterCollection=scfi_barrel.outputProtoClusterCollection,
       outputClusterCollection="EcalBarrelScFiClusters",
       mcHits="EcalBarrelScFiHits",
       logWeightBase=6.2,
       samplingFraction=kwargs['sf'])
    algorithms.append(scfi_barrel_clreco)

podout.outputCommands = ["keep *"]

ApplicationMgr(
    TopAlg=[podin, scfi_barrel_digi, scfi_barrel_reco, scfi_barrel_merger, scfi_barrel_cl, scfi_barrel_clreco, podout],
    EvtSel='NONE',
    EvtMax=kwargs['nev'],
    ExtSvc=[podioevent],
    OutputLevel=DEBUG
)