Update ECal Barrel ML: combine Imaging and ScFi layers

benchmarks/imaging_shower_ML/options/imaging_ml_data.py

@@ -3,13 +3,17 @@ import ROOT
 from Gaudi.Configuration import *
 from GaudiKernel import SystemOfUnits as units
 
-from Configurables import ApplicationMgr, EICDataSvc, PodioOutput, GeoSvc
+from Configurables import ApplicationMgr, EICDataSvc, PodioInput, PodioOutput, GeoSvc
+from GaudiKernel.SystemOfUnits import MeV, GeV, mm, cm, mrad
 
-from Configurables import PodioInput
 from Configurables import Jug__Base__InputCopier_dd4pod__Geant4ParticleCollection_dd4pod__Geant4ParticleCollection_ as MCCopier
-from Configurables import Jug__Digi__CalorimeterHitDigi as CalorimeterHitDigi
-from Configurables import Jug__Reco__ImagingPixelReco as ImagingPixelReco
-from Configurables import Jug__Reco__ImagingPixelMerger as ImagingPixelMerger
+from Configurables import Jug__Digi__CalorimeterHitDigi as CalHitDigi
+from Configurables import Jug__Digi__CalorimeterBirksCorr as CalBirksCorr
+from Configurables import Jug__Reco__CalorimeterHitReco as CalHitReco
+from Configurables import Jug__Reco__CalorimeterHitsMerger as CalHitsMerger
+from Configurables import Jug__Reco__ImagingPixelMerger as MLDataMerger
+from Configurables import Jug__Reco__ImagingPixelDataSorter as MLDataSorter
+from Configurables import Jug__Reco__ImagingPixelDataCombiner as MLDataCombiner
 
 
 # input arguments through environment variables
@@ -18,13 +22,13 @@ kwargs['input'] = os.environ.get('JUGGLER_SIM_FILE', '')
 kwargs['output'] = os.environ.get('JUGGLER_REC_FILE', '')
 kwargs['compact'] = os.environ.get('JUGGLER_COMPACT_PATH', '')
 kwargs['nev'] = int(os.environ.get('JUGGLER_N_EVENTS', -1))
-kwargs['nlayers'] = int(os.environ.get('IMCAL_ML_N_LAYERS', 20))
-kwargs['nhits'] = int(os.environ.get('IMCAL_ML_N_HITS', 20))
+kwargs['combine'] = os.environ.get('IMCAL_ML_COMBINE', 'concatenate')
+kwargs['img_nlayers'] = int(os.environ.get('IMCAL_ML_IMG_NLAYERS', 9))
+kwargs['nhits'] = int(os.environ.get('IMCAL_ML_NHITS', 20))
 
 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('JUGGLER_SAMP_FRAC', '1.0'))
@@ -32,44 +36,124 @@ sf = float(os.environ.get('JUGGLER_SAMP_FRAC', '1.0'))
 geo_service  = GeoSvc('GeoSvc', detectors=[f.strip() for f in kwargs['compact'].split(',')])
 podev = EICDataSvc('EventDataSvc', inputs=[f.strip() for f in kwargs['input'].split(',')])
 
-podin = PodioInput('PodioReader', collections=['mcparticles', 'EcalBarrelHits'])
+podin = PodioInput('PodioReader', collections=['mcparticles', 'EcalBarrelHits', 'EcalBarrelScFiHits'])
 podout = PodioOutput('out', filename=kwargs['output'])
 
 copier = MCCopier('MCCopier',
-        OutputLevel=DEBUG,
+        OutputLevel=WARNING,
         inputCollection='mcparticles',
         outputCollection='mcparticles2')
-imcaldigi = CalorimeterHitDigi('imcal_digi',
+
+# Central Barrel Ecal (Imaging Cal.)
+becal_img_daq = dict(
+        dynamicRangeADC=3*MeV,
+        capacityADC=8192,
+        pedestalMean=400,
+        pedestalSigma=20)   # about 6 keV
+
+becal_img_digi = CalHitDigi('becal_img_digi',
         inputHitCollection='EcalBarrelHits',
-        outputHitCollection='EcalBarrelHitsDigi',
-        energyResolutions=[0., 0.02, 0.],
-        dynamicRangeADC=3*units.MeV,
-        pedestalSigma=40)
-imcalreco = ImagingPixelReco('imcal_reco',
-        inputHitCollection=imcaldigi.outputHitCollection,
-        outputHitCollection='EcalBarrelHitsReco',
-        dynamicRangeADC=3.*units.MeV,
-        pedestalSigma=40,
-        thresholdFactor=3.5,
-        readoutClass='EcalBarrelHits',
-        layerField='layer',
-        sectorField='module')
-imcaldata = ImagingPixelMerger('imcal_merger',
-#        OutputLevel=DEBUG,
-        inputHitCollection=imcalreco.outputHitCollection,
-        outputHitCollection='EcalBarrelHitsML',
+        outputHitCollection='EcalBarrelImagingHitsDigi',
+        energyResolutions=[0., 0.02, 0.],   # 2% flat resolution
+        **becal_img_daq)
+
+becal_img_reco = CalHitReco('becal_img_reco',
+        inputHitCollection=becal_img_digi.outputHitCollection,
+        outputHitCollection='EcalBarrelImagingHitsReco',
+        thresholdFactor=3,  # about 20 keV
+        readoutClass='EcalBarrelHits',  # readout class
+        layerField='layer',             # field to get layer id
+        sectorField='module',           # field to get sector id
+        **becal_img_daq)
+
+becal_img_merger = MLDataMerger('becal_img_merger',
+        inputHitCollection=becal_img_reco.outputHitCollection,
+        outputHitCollection='EcalBarrelImagingHitsSeg',
         etaSize=0.001,
-        phiSize=0.001,
-        numberOfHits=kwargs['nhits'],
-        numberOfLayers=kwargs['nlayers'])
+        phiSize=0.001)
+
+becal_img_sorter = MLDataSorter('becal_img_sorter',
+        inputHitCollection=becal_img_merger.outputHitCollection,
+        outputHitCollection='EcalBarrelImagingHitsML',
+        numberOfLayers=kwargs['img_nlayers'],
+        numberOfHits=kwargs['nhits'])
+
+
+#Central ECAL SciFi
+# use the same daq_setting for digi/reco pair
+becal_scfi_daq = dict(
+        dynamicRangeADC=50.*MeV,
+        capacityADC=32768,
+        pedestalMean=400,
+        pedestalSigma=10)
+# becal_scfi_daq = dict(
+#         dynamicRangeADC=50.*MeV,
+#         capacityADC=2**16,
+#         pedestalMean=100,
+#         pedestalSigma=0)
+
+becal_scfi_birks = CalBirksCorr('becal_scfi_birks',
+        inputHitCollection='EcalBarrelScFiHits',
+        outputHitCollection='EcalBarrelScFiHitsCorr',
+        birksConstant=0.126*mm/MeV)
+
+becal_scfi_digi = CalHitDigi('becal_scfi_digi',
+        inputHitCollection=becal_scfi_birks.outputHitCollection,
+        # inputHitCollection='EcalBarrelScFiHits',
+        outputHitCollection='EcalBarrelScFiHitsDigi',
+        **becal_scfi_daq)
+
+becal_scfi_reco = CalHitReco('becal_scfi_reco',
+        inputHitCollection=becal_scfi_digi.outputHitCollection,
+        outputHitCollection='EcalBarrelScFiHitsReco',
+        thresholdFactor=5.0,
+        readoutClass='EcalBarrelScFiHits',
+        layerField='layer',
+        sectorField='module',
+        localDetFields=['system', 'module'], # use local coordinates in each module (stave)
+        **becal_scfi_daq)
+
+# merge hits in different layer (projection to local x-y plane)
+becal_scfi_merger = CalHitsMerger('becal_scfi_merger',
+        # OutputLevel=DEBUG,
+        inputHitCollection=becal_scfi_reco.outputHitCollection,
+        outputHitCollection='EcalBarrelScFiGridReco',
+        fields=['fiber'],
+        fieldRefNumbers=[1],
+        readoutClass='EcalBarrelScFiHits')
+
+becal_scfi_sorter = MLDataSorter('becal_scfi_sorter',
+        inputHitCollection=becal_scfi_merger.outputHitCollection,
+        outputHitCollection='EcalBarrelScFiHitsML',
+        numberOfLayers=20,
+        numberOfHits=kwargs['nhits'])
+
+# combine layers
+becal_combiner = MLDataCombiner('becal_combiner',
+        inputHitCollection1=becal_img_sorter.outputHitCollection,
+        inputHitCollection2=becal_scfi_sorter.outputHitCollection,
+        outputHitCollection='EcalBarrelHitsCombinedML',
+        layerIncrement=100,
+        rule=kwargs['combine'])
 
-podout.outputCommands = ['keep *']
+podout.outputCommands = [
+#        'keep *',
+        'drop *',
+        'keep mcparticles2',
+        'keep EcalBarrel*ML',
+        'keep *Corr',
+]
 
 ApplicationMgr(
-    TopAlg=[podin, copier, imcaldigi, imcalreco, imcaldata, podout],
+    TopAlg=[podin, copier,
+            becal_img_digi, becal_img_reco, becal_img_merger, becal_img_sorter,
+            becal_scfi_birks,
+            becal_scfi_digi, becal_scfi_reco, becal_scfi_merger, becal_scfi_sorter,
+            becal_combiner,
+            podout],
     EvtSel='NONE',
     EvtMax=kwargs['nev'],
     ExtSvc=[podev],
-    OutputLevel=DEBUG
+    OutputLevel=WARNING
 )