Draft: Resolve "Update pi0 reconstruction"

Closes #68

benchmarks/imaging_ecal/options/imaging_topocluster.py

@@ -2,8 +2,7 @@ import json
 import os
 import ROOT
 from Gaudi.Configuration import *
-from GaudiKernel import SystemOfUnits as units
-
+from GaudiKernel.SystemOfUnits import mm, MeV, rad, ns
 from Configurables import ApplicationMgr, EICDataSvc, PodioOutput, GeoSvc
 
 from Configurables import PodioInput
@@ -20,7 +19,7 @@ with open('config/emcal_barrel_calibration.json') as f:
     calib_data = json.load(f)['electron']
 
 kwargs = dict()
-kwargs['sf'] = float(calib_data['sampling_fraction_img'])
+kwargs['img_sf'] = float(calib_data['sampling_fraction_img'])
 
 # input arguments through environment variables
 kwargs['input'] = os.environ.get('CB_EMCAL_SIM_FILE', '../topside/barrel_pion0_5GeV.root')
@@ -41,52 +40,49 @@ out = PodioOutput("out", filename=kwargs['output'])
 podioinput = PodioInput("PodioReader", collections=["mcparticles", "EcalBarrelHits"], OutputLevel=DEBUG)
 
 copier = MCCopier("MCCopier",
-        # OutputLevel=DEBUG,
         inputCollection="mcparticles",
         outputCollection="mcparticles2")
 calcopier = CalCopier("CalCopier",
-        # OutputLevel=DEBUG,
         inputCollection="EcalBarrelHits",
         outputCollection="EcalBarrelHits2")
 
 # use the same daq_setting for digi/reco pair
 daq_setting = dict(
-        dynamicRangeADC=3*units.MeV,
+        dynamicRangeADC=3*MeV,
         capacityADC=32767,
         pedestalMean=400,
         pedestalSigma=50)   # 50/32767*3 MeV ~ 5 keV
 imcaldigi = CalorimeterHitDigi("imcal_digi",
-        # OutputLevel=DEBUG,
         inputHitCollection="EcalBarrelHits",
         outputHitCollection="DigiEcalBarrelImagingHits",
         energyResolutions=[0., 0.02, 0.],
         **daq_setting)
 imcalreco = ImagingPixelReco("imcal_reco",
-        # OutputLevel=DEBUG,
         inputHitCollection=imcaldigi.outputHitCollection,
         outputHitCollection="RecoEcalBarrelImagingHits",
         readoutClass="EcalBarrelHits",
         layerField="layer",
         sectorField="module",
+        thresholdFactor=5.0,
         **daq_setting)
 
 imcalcluster = ImagingTopoCluster("imcal_cluster",
-        # OutputLevel=DEBUG,
         inputHitCollection=imcalreco.outputHitCollection,
         outputProtoClusterCollection="EcalBarrelImagingProtoClusters",
-        localDistXY=[2.*units.mm, 2*units.mm],
-        layerDistEtaPhi=[10*units.mrad, 10*units.mrad],
+        localDistXY=[2.*mm, 2*mm],
+        layerDistEtaPhi=[0.01*rad, 0.01*rad],
         neighbourLayersRange=2,
-        sectorDist=3.*units.cm,
-        minClusterNhits=5)
+        sectorDist=30.*mm,
+        minClusterNhits=15,
+        minClusterEdep=0.5*MeV)
+
 clusterreco = ImagingClusterReco("imcal_clreco",
-        # OutputLevel=DEBUG,
         inputHitCollection=imcalcluster.inputHitCollection,
         inputProtoClusterCollection=imcalcluster.outputProtoClusterCollection,
         outputLayerCollection="EcalBarrelImagingClustersLayers",
         outputClusterCollection="EcalBarrelImagingClusters",
         outputInfoCollection="EcalBarrelImagingClustersInfo",
-        samplingFraction=sf)
+        samplingFraction=kwargs['img_sf'])
 
 out.outputCommands = ["keep *"]