Skip to content
Snippets Groups Projects
Commit 4692d2e1 authored by Chao Peng's avatar Chao Peng
Browse files

Add ScFi Clusters

parent 28364bf4
No related branches found
No related tags found
1 merge request!125Add ScFi Clusters
Hide whitespace changes
Inline Side-by-side
Showing
with 279 additions and 33 deletions
benchmarks/clustering/full_cal_clusters.sh
+ 5
− 5
View file @ 4692d2e1
......@@ -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
......
benchmarks/clustering/options/full_cal_reco.py
+ 2
− 1
View file @ 4692d2e1
......@@ -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",
......
benchmarks/clustering/scripts/cluster_plots.py
+ 6
− 10
View file @ 4692d2e1
......@@ -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)))
benchmarks/ecal/options/endcap_e.py
+ 2
− 1
View file @ 4692d2e1
......@@ -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",
......
benchmarks/ecal/scripts/draw_cluters.py
+ 3
− 2
View file @ 4692d2e1
......@@ -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)))
benchmarks/imaging_ecal/options/hybrid_cluster.py 0 → 100644
+ 153
− 0
View file @ 4692d2e1
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
)
benchmarks/imaging_ecal/options/scfi_cluster.py 0 → 100644
+ 91
− 0
View file @ 4692d2e1
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
)
benchmarks/imaging_ecal/run_emcal_barrel.sh
+ 7
− 6
View file @ 4692d2e1
......@@ -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
......
benchmarks/imaging_ecal/scripts/draw_cluster.py
+ 4
− 4
View file @ 4692d2e1
......@@ -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
......
benchmarks/imaging_ecal/scripts/draw_cluster_layers.py
+ 4
− 4
View file @ 4692d2e1
......@@ -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
......
benchmarks/imaging_ecal/scripts/utils.py
+ 2
− 0
View file @ 4692d2e1
......@@ -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))
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment