From 4692d2e1571df947ef7165a7a7b2f2be327baec8 Mon Sep 17 00:00:00 2001 From: Chao Peng <cpeng@anl.gov> Date: Wed, 14 Jul 2021 00:11:55 +0000 Subject: [PATCH] Add ScFi Clusters --- benchmarks/clustering/full_cal_clusters.sh | 10 +- .../clustering/options/full_cal_reco.py | 3 +- .../clustering/scripts/cluster_plots.py | 16 +- benchmarks/ecal/options/endcap_e.py | 3 +- benchmarks/ecal/scripts/draw_cluters.py | 5 +- .../imaging_ecal/options/hybrid_cluster.py | 153 ++++++++++++++++++ .../imaging_ecal/options/scfi_cluster.py | 91 +++++++++++ benchmarks/imaging_ecal/run_emcal_barrel.sh | 13 +- .../imaging_ecal/scripts/draw_cluster.py | 8 +- .../scripts/draw_cluster_layers.py | 8 +- benchmarks/imaging_ecal/scripts/utils.py | 2 + 11 files changed, 279 insertions(+), 33 deletions(-) create mode 100644 benchmarks/imaging_ecal/options/hybrid_cluster.py create mode 100644 benchmarks/imaging_ecal/options/scfi_cluster.py diff --git a/benchmarks/clustering/full_cal_clusters.sh b/benchmarks/clustering/full_cal_clusters.sh index ba97c36e..f4dce6fc 100644 --- a/benchmarks/clustering/full_cal_clusters.sh +++ b/benchmarks/clustering/full_cal_clusters.sh @@ -118,11 +118,8 @@ rootls -t "${JUGGLER_SIM_FILE}" # Directory for results mkdir -p results -FULL_CAL_OPTION_DIR=benchmarks/clustering/options -FULL_CAL_SCRIPT_DIR=benchmarks/clustering/scripts - # Run Juggler -# Digitization & Clustering +FULL_CAL_OPTION_DIR=benchmarks/clustering/options xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv \ gaudirun.py ${FULL_CAL_OPTION_DIR}/full_cal_reco.py if [[ "$?" -ne "0" ]] ; then @@ -131,7 +128,10 @@ if [[ "$?" -ne "0" ]] ; then fi # Run analysis scripts -python ${FULL_CAL_SCRIPT_DIR}/cluster_plots.py ${JUGGLER_SIM_FILE} ${JUGGLER_REC_FILE} -o results +FULL_CAL_SCRIPT_DIR=benchmarks/clustering/scripts +python ${FULL_CAL_SCRIPT_DIR}/cluster_plots.py ${JUGGLER_SIM_FILE} ${JUGGLER_REC_FILE} -o results \ + --collections "EcalEndcapNClusters, EcalEndcapPClusters, EcalBarrelClusters, + HcalElectronEndcapClusters, HcalHadronEndcapClusters, HcalBarrelClusters" root_filesize=$(stat --format=%s "${JUGGLER_REC_FILE}") if [[ "${JUGGLER_N_EVENTS}" -lt "500" ]] ; then diff --git a/benchmarks/clustering/options/full_cal_reco.py b/benchmarks/clustering/options/full_cal_reco.py index 830726e4..2769bdeb 100644 --- a/benchmarks/clustering/options/full_cal_reco.py +++ b/benchmarks/clustering/options/full_cal_reco.py @@ -90,7 +90,8 @@ ce_ecal_cl = IslandCluster("ce_ecal_cl", splitCluster=False, minClusterHitEdep=1.0*MeV, # discard low energy hits minClusterCenterEdep=30*MeV, - dimScaledDist=1.8) # dimension scaled dist is good for hybrid sectors with different module size + 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", diff --git a/benchmarks/clustering/scripts/cluster_plots.py b/benchmarks/clustering/scripts/cluster_plots.py index 8802a94f..8d4ec745 100644 --- a/benchmarks/clustering/scripts/cluster_plots.py +++ b/benchmarks/clustering/scripts/cluster_plots.py @@ -103,6 +103,8 @@ if __name__ == '__main__': help='Macro files to be loaded by root, separated by \",\".') parser.add_argument('--mc-branch', dest='mc', default='mcparticles', help='Branch name of MC particles truth info.') + parser.add_argument('--collections', dest='coll', default='', + help='Collection name of clusters to plot') args = parser.parse_args() # multi-threading for RDataFrame @@ -124,14 +126,7 @@ if __name__ == '__main__': rdf_rec = ROOT.RDataFrame('events', args.rec_file) # cluster collections - colls = [ - 'EcalEndcapNClusters', - 'EcalEndcapPClusters', - 'EcalBarrelClusters', - 'HcalElectronEndcapClusters', - 'HcalHadronEndcapClusters', - 'HcalBarrelClusters', - ] + colls = [c.strip() for c in args.coll.split(',')] # a function to extract eta ROOT.gInterpreter.Declare(''' template<typename T> @@ -164,7 +159,8 @@ if __name__ == '__main__': axs[0][0].set_xlabel('Number of Clusters', fontsize=16) for ax, (branch, xlabel, bins) in zip(axs.flat[1:], plots): - ax.hist(df[branch].values, bins=bins, align='mid', ec='k') + ax.hist(df[branch].values, weights=np.repeat(1./float(df.shape[0]), df.shape[0]), + bins=bins, align='mid', ec='k') ax.set_xlabel(xlabel, fontsize=16) # universal axis settings @@ -172,6 +168,6 @@ if __name__ == '__main__': ax.tick_params(labelsize=16, direction='in') ax.set_axisbelow(True) ax.grid(linestyle=':') - fig.text(0.05, 0.5, 'Normalized Counts', rotation=90, va='center', fontsize=18) + fig.text(0.06, 0.5, 'Normalized Counts', rotation=90, va='center', fontsize=18) fig.savefig(os.path.join(args.outdir, '{}.png'.format(coll))) diff --git a/benchmarks/ecal/options/endcap_e.py b/benchmarks/ecal/options/endcap_e.py index 9dee7fb8..b70f0d80 100644 --- a/benchmarks/ecal/options/endcap_e.py +++ b/benchmarks/ecal/options/endcap_e.py @@ -72,7 +72,8 @@ ce_ecal_cl = IslandCluster("ce_ecal_cl", splitCluster=False, minClusterHitEdep=1.0*MeV, # discard low energy hits minClusterCenterEdep=30*MeV, - dimScaledDist=1.8) # hybrid calorimeter with different dimensions, using a scaled dist + 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", inputHitCollection="EcalEndcapNClusterHits", diff --git a/benchmarks/ecal/scripts/draw_cluters.py b/benchmarks/ecal/scripts/draw_cluters.py index d46e1a2d..a346731d 100644 --- a/benchmarks/ecal/scripts/draw_cluters.py +++ b/benchmarks/ecal/scripts/draw_cluters.py @@ -93,7 +93,8 @@ if __name__ == '__main__': axs[0][0].set_xlabel('Number of Clusters', fontsize=16) for ax, (branch, xlabel, bins) in zip(axs.flat[1:], plots): - ax.hist(df['{}.{}'.format(coll, branch)].values, bins=bins, align='mid', ec='k') + ax.hist(df['{}.{}'.format(coll, branch)].values, weights=np.repeat(1./float(df.shape[0]), df.shape[0]), + bins=bins, align='mid', ec='k') ax.set_xlabel(xlabel, fontsize=16) # universal axis settings @@ -101,7 +102,7 @@ if __name__ == '__main__': ax.tick_params(labelsize=16, direction='in') ax.set_axisbelow(True) ax.grid(linestyle=':') - fig.text(0.05, 0.5, 'Normalized Counts', rotation=90, va='center', fontsize=18) + fig.text(0.06, 0.5, 'Normalized Counts', rotation=90, va='center', fontsize=18) fig.savefig(os.path.join(args.outdir, '{}.png'.format(coll))) diff --git a/benchmarks/imaging_ecal/options/hybrid_cluster.py b/benchmarks/imaging_ecal/options/hybrid_cluster.py new file mode 100644 index 00000000..75705fc0 --- /dev/null +++ b/benchmarks/imaging_ecal/options/hybrid_cluster.py @@ -0,0 +1,153 @@ +import os +import ROOT +from Gaudi.Configuration import * +from GaudiKernel.SystemOfUnits import MeV, mm, cm, mrad, rad, ns + +from GaudiKernel.DataObjectHandleBase import DataObjectHandleBase +from Configurables import ApplicationMgr, EICDataSvc, PodioOutput, GeoSvc + +from Configurables import PodioInput +from Configurables import Jug__Base__InputCopier_dd4pod__Geant4ParticleCollection_dd4pod__Geant4ParticleCollection_ as MCCopier +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 +from Configurables import Jug__Reco__ImagingPixelReco as ImagingPixelReco +from Configurables import Jug__Reco__ImagingTopoCluster as ImagingTopoCluster +from Configurables import Jug__Reco__ImagingClusterReco as ImagingClusterReco + + +# input arguments through environment variables +kwargs = dict() +kwargs['img_sf'] = float(os.environ.get('CB_EMCAL_IMG_SF', '0.01324')) +kwargs['scfi_sf'] = float(os.environ.get('CB_EMCAL_SCFI_SF', '0.0938')) +# raise error if not defined +kwargs['input'] = os.environ['CB_EMCAL_SIM_FILE'] +kwargs['output'] = os.environ['CB_EMCAL_REC_FILE'] +kwargs['compact'] = os.environ['CB_EMCAL_COMPACT_PATH'] +kwargs['nev'] = int(os.environ['CB_EMCAL_NUMEV']) + +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 +geo_service = GeoSvc('GeoSvc', detectors=kwargs['compact'].split(','), OutputLevel=INFO) +podioevent = EICDataSvc('EventDataSvc', inputs=kwargs['input'].split(','), OutputLevel=DEBUG) + +sim_colls = [ + 'mcparticles', + 'EcalBarrelHits', + 'EcalBarrelScFiHits', +] +podin = PodioInput('PodioReader', collections=sim_colls, OutputLevel=DEBUG) +podout = PodioOutput('podout', filename=kwargs['output']) + +mccopier = MCCopier('MCCopier', + # OutputLevel=DEBUG, + inputCollection='mcparticles', + outputCollection='mcparticles2') + +# use the same daq_setting for digi/reco pair +imcaldaq = dict( + dynamicRangeADC=3*MeV, + capacityADC=32767, + pedestalMean=400, + pedestalSigma=50) # 50/32767*3 MeV ~ 5 keV + +imcaldigi = CalHitDigi('imcal_digi', + # OutputLevel=DEBUG, + inputHitCollection='EcalBarrelHits', + outputHitCollection='DigiEcalBarrelHits', + energyResolutions=[0., 0.02, 0.], + **imcaldaq) +imcalreco = ImagingPixelReco('imcal_reco', + # OutputLevel=DEBUG, + inputHitCollection='DigiEcalBarrelHits', + outputHitCollection='RecoEcalBarrelHits', + readoutClass='EcalBarrelHits', + layerField='layer', + sectorField='module', + **imcaldaq) + +imcalcluster = ImagingTopoCluster('imcal_cluster', + # OutputLevel=DEBUG, + inputHitCollection='RecoEcalBarrelHits', + outputHitCollection='EcalBarrelClusterHits', + localDistXY=[2.*mm, 2*mm], + layerDistEtaPhi=[10*mrad, 10*mrad], + neighbourLayersRange=2, + sectorDist=3.*cm) +clusterreco = ImagingClusterReco('imcal_clreco', + # OutputLevel=DEBUG, + inputHitCollection='EcalBarrelClusterHits', + outputLayerCollection='EcalBarrelClustersLayers', + outputClusterCollection='EcalBarrelClusters', + samplingFraction=kwargs['img_sf']) + +# scfi layers +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", + 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=kwargs['scfi_sf']) + +# TODO: merge two types of clusters + +podout.outputCommands = ['drop *', + 'keep mcparticles2', + 'keep *Reco*', + 'keep *Digi*', + 'keep *Cluster*', + 'keep *Layer*', + ] + +ApplicationMgr( + TopAlg=[podin, mccopier, + imcaldigi, imcalreco, imcalcluster, clusterreco, + 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 +) diff --git a/benchmarks/imaging_ecal/options/scfi_cluster.py b/benchmarks/imaging_ecal/options/scfi_cluster.py new file mode 100644 index 00000000..6775a0aa --- /dev/null +++ b/benchmarks/imaging_ecal/options/scfi_cluster.py @@ -0,0 +1,91 @@ +import os +import ROOT +from Gaudi.Configuration import * +from GaudiKernel.SystemOfUnits import mm, MeV, rad, ns +from GaudiKernel.DataObjectHandleBase import DataObjectHandleBase +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 + + +# input arguments through environment variables +kwargs = dict() +kwargs['sf'] = float(os.environ.get('CB_EMCAL_SCFI, SAMP_FRAC', '0.0938')) +kwargs['input'] = os.environ.get('CB_EMCAL_SIM_FILE', '../topside/barrel_pion0_5GeV.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) + +podin = PodioInput("PodioReader", collections=["mcparticles", "EcalBarrelScFiHits"], OutputLevel=DEBUG) +podout = PodioOutput("out", filename=kwargs['output']) + +# 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", + 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=kwargs['sf']) + + +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 +) + diff --git a/benchmarks/imaging_ecal/run_emcal_barrel.sh b/benchmarks/imaging_ecal/run_emcal_barrel.sh index 8584bef2..73bae02c 100644 --- a/benchmarks/imaging_ecal/run_emcal_barrel.sh +++ b/benchmarks/imaging_ecal/run_emcal_barrel.sh @@ -103,23 +103,24 @@ rootls -t "${CB_EMCAL_SIM_FILE}" # Directory for plots mkdir -p results -# CB_EMCAL_OPTION_DIR=${JUGGLER_INSTALL_PREFIX}/Examples/options -# CB_EMCAL_SCRIPT_DIR=${JUGGLER_INSTALL_PREFIX}/Examples/scripts/imaging_calorimeter CB_EMCAL_OPTION_DIR=benchmarks/imaging_ecal/options -CB_EMCAL_SCRIPT_DIR=benchmarks/imaging_ecal/scripts - # Run Juggler xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv \ - gaudirun.py ${CB_EMCAL_OPTION_DIR}/imaging_topocluster.py + gaudirun.py ${CB_EMCAL_OPTION_DIR}/hybrid_cluster.py # gaudirun.py ${CB_EMCAL_OPTION_DIR}/imaging_topocluster.py if [[ "$?" -ne "0" ]] ; then echo "ERROR running juggler" exit 1 fi +# Plot clusters first +FULL_CAL_SCRIPT_DIR=benchmarks/clustering/scripts +python ${FULL_CAL_SCRIPT_DIR}/cluster_plots.py ${CB_EMCAL_SIM_FILE} ${CB_EMCAL_REC_FILE} -o results \ + --collections "EcalBarrelClusters, EcalBarrelScFiClusters" + # check required python modules python -m pip install -r benchmarks/imaging_ecal/requirements.txt - +CB_EMCAL_SCRIPT_DIR=benchmarks/imaging_ecal/scripts IFS=',' read -ra ADDR <<< "$CB_EMCAL_IEV" for iev in "${ADDR[@]}"; do if [[ $iev -ge "${CB_EMCAL_NUMEV}" ]] ; then diff --git a/benchmarks/imaging_ecal/scripts/draw_cluster.py b/benchmarks/imaging_ecal/scripts/draw_cluster.py index d03db476..88f338d0 100644 --- a/benchmarks/imaging_ecal/scripts/draw_cluster.py +++ b/benchmarks/imaging_ecal/scripts/draw_cluster.py @@ -114,10 +114,10 @@ if __name__ == '__main__': # we can read these values from xml file desc = compact_constants(args.compact, [ - 'EcalBarrel_rmin', - 'EcalBarrel_ReadoutLayerThickness', - 'EcalBarrel_ReadoutLayerNumber', - 'EcalBarrel_length', +# 'EcalBarrel_rmin', +# 'EcalBarrel_ReadoutLayerThickness', +# 'EcalBarrel_ReadoutLayerNumber', +# 'EcalBarrel_length', ]) if not len(desc): # or define Ecal shapes diff --git a/benchmarks/imaging_ecal/scripts/draw_cluster_layers.py b/benchmarks/imaging_ecal/scripts/draw_cluster_layers.py index 25f102e0..7baafebb 100644 --- a/benchmarks/imaging_ecal/scripts/draw_cluster_layers.py +++ b/benchmarks/imaging_ecal/scripts/draw_cluster_layers.py @@ -68,10 +68,10 @@ if __name__ == '__main__': # we can read these values from xml file desc = compact_constants(args.compact, [ - 'EcalBarrel_rmin', - 'EcalBarrel_ReadoutLayerThickness', - 'EcalBarrel_ReadoutLayerNumber', - 'EcalBarrel_length', +# 'EcalBarrel_rmin', +# 'EcalBarrel_ReadoutLayerThickness', +# 'EcalBarrel_ReadoutLayerNumber', +# 'EcalBarrel_length', ]) if not len(desc): # or define Ecal shapes diff --git a/benchmarks/imaging_ecal/scripts/utils.py b/benchmarks/imaging_ecal/scripts/utils.py index 8c9bfadf..da9954ac 100644 --- a/benchmarks/imaging_ecal/scripts/utils.py +++ b/benchmarks/imaging_ecal/scripts/utils.py @@ -157,6 +157,8 @@ def compact_constants(path, names): if not os.path.exists(path): print('Cannot find compact file \"{}\".'.format(path)) return [] + if not names: + return [] kernel = DDG4.Kernel() description = kernel.detectorDescription() kernel.loadGeometry("file:{}".format(path)) -- GitLab