From 203880ac3e5cb98df19cfd81e65717f0b0c9026b Mon Sep 17 00:00:00 2001 From: Sylvester Joosten <sjoosten@anl.gov> Date: Sat, 7 Aug 2021 21:35:39 +0000 Subject: [PATCH] Update eicd --- .gitignore | 10 +++ .../full_trackpluscalo_reco.py | 0 .../clustering/options/full_cal_reco.py | 83 ++++++++++-------- .../clustering/scripts/cluster_plots.py | 9 +- benchmarks/ecal/options/barrel.py | 15 ++-- benchmarks/ecal/options/endcap_e.py | 13 +-- benchmarks/ecal/options/endcap_i.py | 15 ++-- benchmarks/ecal/run_emcal_benchmarks.py | 8 +- .../{draw_cluters.py => draw_clusters.py} | 9 +- .../full/options/full_reconstruction.py | 86 +++++++++++-------- .../imaging_ecal/options/hybrid_cluster.py | 22 +++-- .../imaging_ecal/options/imaging_2dcluster.py | 6 +- .../options/imaging_topocluster.py | 10 ++- .../imaging_ecal/options/scfi_cluster.py | 12 +-- .../imaging_ecal/scripts/draw_cluster.py | 14 +-- .../scripts/draw_cluster_layers.py | 6 +- .../imaging_ecal/scripts/energy_profile.py | 22 ++--- benchmarks/imaging_ecal/scripts/utils.py | 18 ++-- .../options/imaging_ml_data.py | 4 +- .../options/imagingcluster_seeded_tracking.py | 10 ++- options/tracker_reconstruction.py | 5 +- 21 files changed, 218 insertions(+), 159 deletions(-) rename benchmarks/clustering/options/{ => deprecated}/full_trackpluscalo_reco.py (100%) rename benchmarks/ecal/scripts/{draw_cluters.py => draw_clusters.py} (92%) diff --git a/.gitignore b/.gitignore index ec2bfe82..8fdeae71 100644 --- a/.gitignore +++ b/.gitignore @@ -41,3 +41,13 @@ __pycache__/ calorimeters/test/ *.d *.pcm + +# transient directories and files +setup +results +*.root +fieldmaps +eic-env.sh +sim_output +*.obj +*.hepmc diff --git a/benchmarks/clustering/options/full_trackpluscalo_reco.py b/benchmarks/clustering/options/deprecated/full_trackpluscalo_reco.py similarity index 100% rename from benchmarks/clustering/options/full_trackpluscalo_reco.py rename to benchmarks/clustering/options/deprecated/full_trackpluscalo_reco.py diff --git a/benchmarks/clustering/options/full_cal_reco.py b/benchmarks/clustering/options/full_cal_reco.py index 402bc11e..3f23faa1 100644 --- a/benchmarks/clustering/options/full_cal_reco.py +++ b/benchmarks/clustering/options/full_cal_reco.py @@ -76,7 +76,7 @@ ce_ecal_digi = CalHitDigi("ce_ecal_digi", **ce_ecal_daq) ce_ecal_reco = CalHitReco("ce_ecal_reco", - inputHitCollection="EcalEndcapNHitsDigi", + inputHitCollection=ce_ecal_digi.outputHitCollection, outputHitCollection="EcalEndcapNHitsReco", thresholdFactor=4, # 4 sigma cut on pedestal sigma readoutClass="EcalEndcapNHits", @@ -85,8 +85,8 @@ ce_ecal_reco = CalHitReco("ce_ecal_reco", ce_ecal_cl = IslandCluster("ce_ecal_cl", # OutputLevel=DEBUG, - inputHitCollection="EcalEndcapNHitsReco", - outputHitCollection="EcalEndcapNClusterHits", + inputHitCollection=ce_ecal_reco.outputHitCollection, + outputProtoClusterCollection="EcalEndcapNProtoClusters", splitCluster=False, minClusterHitEdep=1.0*MeV, # discard low energy hits minClusterCenterEdep=30*MeV, @@ -94,8 +94,10 @@ ce_ecal_cl = IslandCluster("ce_ecal_cl", 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", + inputHitCollection=ce_ecal_cl.inputHitCollection, + inputProtoClusterCollection=ce_ecal_cl.outputProtoClusterCollection, outputClusterCollection="EcalEndcapNClusters", + outputInfoCollection="EcalEndcapNClustersInfo", samplingFraction=0.998, # this accounts for a small fraction of leakage logWeightBase=4.6) @@ -112,7 +114,7 @@ ci_ecal_digi = CalHitDigi("ci_ecal_digi", **ci_ecal_daq) ci_ecal_reco = CalHitReco("ci_ecal_reco", - inputHitCollection="EcalEndcapPHitsDigi", + inputHitCollection=ci_ecal_digi.outputHitCollection, outputHitCollection="EcalEndcapPHitsReco", thresholdFactor=5.0, **ci_ecal_daq) @@ -120,7 +122,7 @@ ci_ecal_reco = CalHitReco("ci_ecal_reco", # merge hits in different layer (projection to local x-y plane) ci_ecal_merger = CalHitsMerger("ci_ecal_merger", # OutputLevel=DEBUG, - inputHitCollection="EcalEndcapPHitsReco", + inputHitCollection=ci_ecal_reco.outputHitCollection, outputHitCollection="EcalEndcapPHitsRecoXY", fields=["layer", "slice"], fieldRefNumbers=[1, 0], @@ -128,15 +130,17 @@ ci_ecal_merger = CalHitsMerger("ci_ecal_merger", ci_ecal_cl = IslandCluster("ci_ecal_cl", # OutputLevel=DEBUG, - inputHitCollection="EcalEndcapPHitsRecoXY", - outputHitCollection="EcalEndcapPClusterHits", + inputHitCollection=ci_ecal_merger.inputHitCollection, + outputProtoClusterCollection="EcalEndcapPProtoClusters", splitCluster=False, minClusterCenterEdep=10.*MeV, localDistXY=[10*mm, 10*mm]) ci_ecal_clreco = RecoCoG("ci_ecal_clreco", - inputHitCollection="EcalEndcapPClusterHits", + inputHitCollection=ci_ecal_cl.inputHitCollection, + inputProtoClusterCollection=ci_ecal_cl.outputProtoClusterCollection, outputClusterCollection="EcalEndcapPClusters", + outputInfoCollection="EcalEndcapPClustersInfo", logWeightBase=6.2, samplingFraction=ci_ecal_sf) @@ -154,7 +158,7 @@ cb_ecal_digi = CalHitDigi("cb_ecal_digi", **cb_ecal_daq) cb_ecal_reco = ImCalPixelReco("cb_ecal_reco", - inputHitCollection="EcalBarrelHitsDigi", + inputHitCollection=cb_ecal_digi.outputHitCollection, outputHitCollection="EcalBarrelHitsReco", thresholdFactor=3, # about 20 keV readoutClass="EcalBarrelHits", # readout class @@ -163,8 +167,8 @@ cb_ecal_reco = ImCalPixelReco("cb_ecal_reco", **cb_ecal_daq) cb_ecal_cl = ImagingCluster("cb_ecal_cl", - inputHitCollection="EcalBarrelHitsReco", - outputHitCollection="EcalBarrelClusterHits", + inputHitCollection=cb_ecal_reco.outputHitCollection, + outputProtoClusterCollection="EcalBarrelProtoClusters", localDistXY=[2.*mm, 2*mm], # same layer layerDistEtaPhi=[10*mrad, 10*mrad], # adjacent layer neighbourLayersRange=2, # id diff for adjacent layer @@ -172,8 +176,10 @@ cb_ecal_cl = ImagingCluster("cb_ecal_cl", cb_ecal_clreco = ImagingClusterReco("cb_ecal_clreco", samplingFraction=cb_ecal_sf, - inputHitCollection="EcalBarrelClusterHits", + inputHitCollection=cb_ecal_cl.inputHitCollection, + inputProtoClusterCollection=cb_ecal_cl.outputProtoClusterCollection, outputClusterCollection="EcalBarrelClusters", + outputInfoCollection="EcalBarrelClustersInfo", outputLayerCollection="EcalBarrelLayers") #Central ECAL SciFi @@ -190,7 +196,7 @@ scfi_barrel_digi = CalHitDigi("scfi_barrel_digi", **scfi_barrel_daq) scfi_barrel_reco = CalHitReco("scfi_barrel_reco", - inputHitCollection="EcalBarrelScFiHitsDigi", + inputHitCollection=scfi_barrel_digi.outputHitCollection, outputHitCollection="EcalBarrelScFiHitsReco", thresholdFactor=5.0, readoutClass="EcalBarrelScFiHits", @@ -202,7 +208,7 @@ scfi_barrel_reco = CalHitReco("scfi_barrel_reco", # merge hits in different layer (projection to local x-y plane) scfi_barrel_merger = CalHitsMerger("scfi_barrel_merger", # OutputLevel=DEBUG, - inputHitCollection="EcalBarrelScFiHitsReco", + inputHitCollection=scfi_barrel_reco.outputHitCollection, outputHitCollection="EcalBarrelScFiGridReco", fields=["fiber"], fieldRefNumbers=[1], @@ -210,15 +216,17 @@ scfi_barrel_merger = CalHitsMerger("scfi_barrel_merger", scfi_barrel_cl = IslandCluster("scfi_barrel_cl", # OutputLevel=DEBUG, - inputHitCollection="EcalBarrelScFiGridReco", - outputHitCollection="EcalBarrelScFiClusterHits", + inputHitCollection=scfi_barrel_merger.outputHitCollection, + outputProtoClusterCollection="EcalBarrelScFiProtoClusters", splitCluster=False, minClusterCenterEdep=10.*MeV, localDistXZ=[30*mm, 30*mm]) scfi_barrel_clreco = RecoCoG("scfi_barrel_clreco", - inputHitCollection="EcalBarrelScFiClusterHits", + inputHitCollection=scfi_barrel_cl.inputHitCollection, + inputProtoClusterCollection=scfi_barrel_cl.outputProtoClusterCollection, outputClusterCollection="EcalBarrelScFiClusters", + outputInfoCollection="EcalBarrelScFiClustersInfo", logWeightBase=6.2, samplingFraction= scifi_barrel_sf) @@ -236,7 +244,7 @@ cb_hcal_digi = CalHitDigi("cb_hcal_digi", **cb_hcal_daq) cb_hcal_reco = CalHitReco("cb_hcal_reco", - inputHitCollection="HcalBarrelHitsDigi", + inputHitCollection=cb_hcal_digi.outputHitCollection, outputHitCollection="HcalBarrelHitsReco", thresholdFactor=5.0, readoutClass="HcalBarrelHits", @@ -245,22 +253,24 @@ cb_hcal_reco = CalHitReco("cb_hcal_reco", **cb_hcal_daq) cb_hcal_merger = CalHitsMerger("cb_hcal_merger", - inputHitCollection="HcalBarrelHitsReco", + inputHitCollection=cb_hcal_reco.outputHitCollection, outputHitCollection="HcalBarrelHitsRecoXY", readoutClass="HcalBarrelHits", fields=["layer", "slice"], fieldRefNumbers=[1, 0]) cb_hcal_cl = IslandCluster("cb_hcal_cl", - inputHitCollection="HcalBarrelHitsRecoXY", - outputHitCollection="HcalBarrelClusterHits", + inputHitCollection=cb_hcal_merger.outputHitCollection, + outputProtoClusterCollection="HcalBarrelProtoClusters", splitCluster=False, minClusterCenterEdep=30.*MeV, localDistXY=[15.*cm, 15.*cm]) cb_hcal_clreco = RecoCoG("cb_hcal_clreco", - inputHitCollection="HcalBarrelClusterHits", + inputHitCollection=cb_hcal_cl.inputHitCollection, + inputProtoClusterCollection=cb_hcal_cl.outputProtoClusterCollection, outputClusterCollection="HcalBarrelClusters", + outputInfoCollection="HcalBarrelClustersInfo", logWeightBase=6.2, samplingFraction=cb_hcal_sf) @@ -277,28 +287,30 @@ ci_hcal_digi = CalHitDigi("ci_hcal_digi", **ci_hcal_daq) ci_hcal_reco = CalHitReco("ci_hcal_reco", - inputHitCollection="HcalHadronEndcapHitsDigi", + inputHitCollection=ci_hcal_digi.outputHitCollection, outputHitCollection="HcalHadronEndcapHitsReco", thresholdFactor=5.0, **ci_hcal_daq) ci_hcal_merger = CalHitsMerger("ci_hcal_merger", - inputHitCollection="HcalHadronEndcapHitsReco", + inputHitCollection=ci_hcal_reco.outputHitCollection, outputHitCollection="HcalHadronEndcapHitsRecoXY", readoutClass="HcalHadronEndcapHits", fields=["layer", "slice"], fieldRefNumbers=[1, 0]) ci_hcal_cl = IslandCluster("ci_hcal_cl", - inputHitCollection="HcalHadronEndcapHitsRecoXY", - outputHitCollection="HcalHadronEndcapClusterHits", + inputHitCollection=ci_hcal_merger.outputHitCollection, + outputProtoClusterCollection="HcalHadronEndcapProtoClusters", splitCluster=False, minClusterCenterEdep=30.*MeV, localDistXY=[15.*cm, 15.*cm]) ci_hcal_clreco = RecoCoG("ci_hcal_clreco", - inputHitCollection="HcalHadronEndcapClusterHits", + inputHitCollection=ci_hcal_cl.inputHitCollection, + inputProtoClusterCollection=ci_hcal_cl.outputProtoClusterCollection, outputClusterCollection="HcalHadronEndcapClusters", + outputInfoCollection="HcalHadronEndcapClustersInfo", logWeightBase=6.2, samplingFraction=ci_hcal_sf) @@ -315,28 +327,30 @@ ce_hcal_digi = CalHitDigi("ce_hcal_digi", **ce_hcal_daq) ce_hcal_reco = CalHitReco("ce_hcal_reco", - inputHitCollection="HcalElectronEndcapHitsDigi", + inputHitCollection=ce_hcal_digi.outputHitCollection, outputHitCollection="HcalElectronEndcapHitsReco", thresholdFactor=5.0, **ce_hcal_daq) ce_hcal_merger = CalHitsMerger("ce_hcal_merger", - inputHitCollection="HcalElectronEndcapHitsReco", + inputHitCollection=ce_hcal_reco.outputHitCollection, outputHitCollection="HcalElectronEndcapHitsRecoXY", readoutClass="HcalElectronEndcapHits", fields=["layer", "slice"], fieldRefNumbers=[1, 0]) ce_hcal_cl = IslandCluster("ce_hcal_cl", - inputHitCollection="HcalElectronEndcapHitsRecoXY", - outputHitCollection="HcalElectronEndcapClusterHits", + inputHitCollection=ce_hcal_merger.outputHitCollection, + outputProtoClusterCollection="HcalElectronEndcapProtoClusters", splitCluster=False, minClusterCenterEdep=30.*MeV, localDistXY=[15.*cm, 15.*cm]) ce_hcal_clreco = RecoCoG("ce_hcal_clreco", - inputHitCollection="HcalElectronEndcapClusterHits", + inputHitCollection=ce_hcal_cl.inputHitCollection, + inputProtoClusterCollection=ce_hcal_cl.outputProtoClusterCollection, outputClusterCollection="HcalElectronEndcapClusters", + outputInfoCollection="HcalElectronEndcapClustersInfo", logWeightBase=6.2, samplingFraction=ce_hcal_sf) @@ -345,8 +359,7 @@ podout.outputCommands = ['drop *', 'keep mcparticles2', 'keep *Digi', 'keep *Reco*', - 'keep *ClusterHits', - 'keep *Clusters', + 'keep *Cluster*', 'keep *Layers'] ApplicationMgr( diff --git a/benchmarks/clustering/scripts/cluster_plots.py b/benchmarks/clustering/scripts/cluster_plots.py index 336c35d4..d5a890ac 100644 --- a/benchmarks/clustering/scripts/cluster_plots.py +++ b/benchmarks/clustering/scripts/cluster_plots.py @@ -25,6 +25,7 @@ def load_root_macros(arg_macros): def flatten_collection(rdf, collection, cols=None): if not cols: cols = [str(c) for c in rdf.GetColumnNames() if str(c).startswith('{}.'.format(collection))] + cols += [str(c) for c in rdf.GetColumnNames() if str(c).startswith('{}Info.'.format(collection))] else: cols = ['{}.{}'.format(collection, c) for c in cols] if not cols: @@ -52,6 +53,7 @@ def thrown_particles_figure(rdf, save, mcbranch="mcparticles"): # define truth particle info dft = flatten_collection(rdf, mcbranch, ['genStatus', 'pdgID', 'psx', 'psy', 'psz', 'mass']) dft.rename(columns={c: c.replace(mcbranch + '.', '') for c in dft.columns}, inplace=True) + dft.rename(columns={c: c.replace(mcbranch + 'Info.', '') for c in dft.columns}, inplace=True) # select thrown particles dft = dft[dft['genStatus'] == 1] @@ -152,17 +154,16 @@ if __name__ == '__main__': ('eta', '$\eta$', 50), ] for coll in colls: + print(coll) df = flatten_collection(rdf_rec, coll) if not df.shape[0]: print('{} do not have valid entries, skip it'.format(coll)) continue df.rename(columns={c: c.replace(coll + '.', '') for c in df.columns}, inplace=True) - # calculate eta - if 'eta' not in df.columns: - df.loc[:, 'eta'] = -np.log(np.tan(df['polar.theta'].values/2.)) + df.rename(columns={c: c.replace(coll + 'Info.', '') for c in df.columns}, inplace=True) # print(df[['eta', 'polar.theta', 'position.x', 'position.y', 'position.z']]) fig, axs = plt.subplots(2, 2, figsize=(12, 8), dpi=160) - ncl = df.groupby('event')['clusterID'].nunique().values + ncl = df.groupby('event')['ID.value'].nunique().values axs[0][0].hist(ncl, weights=np.repeat(1./float(ncl.shape[0]), ncl.shape[0]), bins=np.arange(0, 10), align='mid', ec='k') axs[0][0].set_xlabel('Number of Clusters', fontsize=16) diff --git a/benchmarks/ecal/options/barrel.py b/benchmarks/ecal/options/barrel.py index e2df5b3b..870f6bf0 100644 --- a/benchmarks/ecal/options/barrel.py +++ b/benchmarks/ecal/options/barrel.py @@ -58,7 +58,7 @@ cb_ecal_digi = CalHitDigi("cb_ecal_digi", **cb_ecal_daq) cb_ecal_reco = ImCalPixelReco("cb_ecal_reco", - inputHitCollection="EcalBarrelHitsDigi", + inputHitCollection=cb_ecal_digi.outputHitCollection, outputHitCollection="EcalBarrelHitsReco", thresholdFactor=3, # about 20 keV readoutClass="EcalBarrelHits", # readout class @@ -67,8 +67,8 @@ cb_ecal_reco = ImCalPixelReco("cb_ecal_reco", **cb_ecal_daq) cb_ecal_cl = ImagingCluster("cb_ecal_cl", - inputHitCollection="EcalBarrelHitsReco", - outputHitCollection="EcalBarrelClusterHits", + inputHitCollection=cb_ecal_reco.outputHitCollection, + outputProtoClusterCollection="EcalBarrelProtoClusters", localDistXY=[2.*mm, 2*mm], # same layer layerDistEtaPhi=[10*mrad, 10*mrad], # adjacent layer neighbourLayersRange=2, # id diff for adjacent layer @@ -76,16 +76,17 @@ cb_ecal_cl = ImagingCluster("cb_ecal_cl", cb_ecal_clreco = ImagingClusterReco("cb_ecal_clreco", samplingFraction=cb_ecal_sf, - inputHitCollection="EcalBarrelClusterHits", + inputHitCollection=cb_ecal_cl.inputHitCollection, + inputProtoClusterCollection=cb_ecal_cl.outputProtoClusterCollection, outputClusterCollection="EcalBarrelClusters", - outputLayerCollection="EcalBarrelLayers") + outputLayerCollection="EcalBarrelLayers", + outputInfoCollection="EcalBarrelClustersInfo") podout.outputCommands = ['drop *', 'keep mcparticles2', 'keep *HitsReco', 'keep *HitsDigi', - 'keep *ClusterHits', - 'keep *Clusters'] + 'keep *Cluster*'] ApplicationMgr( TopAlg = [podin, copier, diff --git a/benchmarks/ecal/options/endcap_e.py b/benchmarks/ecal/options/endcap_e.py index fc907180..36b9376f 100644 --- a/benchmarks/ecal/options/endcap_e.py +++ b/benchmarks/ecal/options/endcap_e.py @@ -57,7 +57,7 @@ ce_ecal_digi = CalHitDigi("ce_ecal_digi", **ce_ecal_daq) ce_ecal_reco = CalHitReco("ce_ecal_reco", - inputHitCollection="EcalEndcapNHitsDigi", + inputHitCollection=ce_ecal_digi.outputHitCollection, outputHitCollection="EcalEndcapNHitsReco", thresholdFactor=4, # 4 sigma cut on pedestal sigma readoutClass="EcalEndcapNHits", @@ -66,8 +66,8 @@ ce_ecal_reco = CalHitReco("ce_ecal_reco", ce_ecal_cl = IslandCluster("ce_ecal_cl", # OutputLevel=DEBUG, - inputHitCollection="EcalEndcapNHitsReco", - outputHitCollection="EcalEndcapNClusterHits", + inputHitCollection=ce_ecal_reco.outputHitCollection, + outputProtoClusterCollection="EcalEndcapNProtoClusters", splitCluster=False, minClusterHitEdep=1.0*MeV, # discard low energy hits minClusterCenterEdep=30*MeV, @@ -75,8 +75,10 @@ ce_ecal_cl = IslandCluster("ce_ecal_cl", dimScaledLocalDistXY=[1.8, 1.8]) # hybrid calorimeter with different dimensions, using a scaled dist ce_ecal_clreco = RecoCoG("ce_ecal_clreco", - inputHitCollection="EcalEndcapNClusterHits", + inputHitCollection=ce_ecal_cl.inputHitCollection, + inputProtoClusterCollection=ce_ecal_cl.outputProtoClusterCollection, outputClusterCollection="EcalEndcapNClusters", + outputInfoCollection="EcalEndcapNClustersInfo", samplingFraction=0.998, # this accounts for a small fraction of leakage logWeightBase=4.6) @@ -84,8 +86,7 @@ podout.outputCommands = ['drop *', 'keep mcparticles2', 'keep *HitsReco', 'keep *HitsDigi', - 'keep *ClusterHits', - 'keep *Clusters'] + 'keep *Cluster*'] ApplicationMgr( TopAlg = [podin, copier, diff --git a/benchmarks/ecal/options/endcap_i.py b/benchmarks/ecal/options/endcap_i.py index 48b6408e..43a2b158 100644 --- a/benchmarks/ecal/options/endcap_i.py +++ b/benchmarks/ecal/options/endcap_i.py @@ -59,7 +59,7 @@ ci_ecal_digi = CalHitDigi("ci_ecal_digi", **ci_ecal_daq) ci_ecal_reco = CalHitReco("ci_ecal_reco", - inputHitCollection="EcalEndcapPHitsDigi", + inputHitCollection=ci_ecal_digi.outputHitCollection, outputHitCollection="EcalEndcapPHitsReco", thresholdFactor=5.0, **ci_ecal_daq) @@ -67,7 +67,7 @@ ci_ecal_reco = CalHitReco("ci_ecal_reco", # merge hits in different layer (projection to local x-y plane) ci_ecal_merger = CalHitsMerger("ci_ecal_merger", # OutputLevel=DEBUG, - inputHitCollection="EcalEndcapPHitsReco", + inputHitCollection=ci_ecal_reco.outputHitCollection, outputHitCollection="EcalEndcapPHitsRecoXY", fields=["layer", "slice"], fieldRefNumbers=[1, 0], @@ -75,15 +75,17 @@ ci_ecal_merger = CalHitsMerger("ci_ecal_merger", ci_ecal_cl = IslandCluster("ci_ecal_cl", # OutputLevel=DEBUG, - inputHitCollection="EcalEndcapPHitsRecoXY", - outputHitCollection="EcalEndcapPClusterHits", + inputHitCollection=ci_ecal_merger.outputHitCollection, + outputProtoClusterCollection="EcalEndcapPProtoClusters", splitCluster=False, minClusterCenterEdep=10.*MeV, localDistXY=[10*mm, 10*mm]) ci_ecal_clreco = RecoCoG("ci_ecal_clreco", - inputHitCollection="EcalEndcapPClusterHits", + inputHitCollection=ci_ecal_cl.inputHitCollection, + inputProtoClusterCollection=ci_ecal_cl.outputProtoClusterCollection, outputClusterCollection="EcalEndcapPClusters", + outputInfoCollection="EcalEndcapPClustersInfo", logWeightBase=6.2, samplingFraction=ci_ecal_sf) @@ -91,8 +93,7 @@ podout.outputCommands = ['drop *', 'keep mcparticles2', 'keep *HitsReco*', 'keep *HitsDigi', - 'keep *ClusterHits', - 'keep *Clusters'] + 'keep *Cluster*'] ApplicationMgr( TopAlg = [podin, copier, diff --git a/benchmarks/ecal/run_emcal_benchmarks.py b/benchmarks/ecal/run_emcal_benchmarks.py index f3ecf610..2380e404 100755 --- a/benchmarks/ecal/run_emcal_benchmarks.py +++ b/benchmarks/ecal/run_emcal_benchmarks.py @@ -15,19 +15,19 @@ import argparse default_type = { 'endcap_e': [ ['endcap_e.py'], - ['draw_cluters.py'], + ['draw_clusters.py'], ['EcalEndcapNClusters']], 'endcap_i': [ ['endcap_i.py'], - ['draw_cluters.py'], + ['draw_clusters.py'], ['EcalEndcapPClusters']], 'barrel': [ ['barrel.py'], - ['draw_cluters.py'], + ['draw_clusters.py'], ['EcalBarrelClusters']], # 'all': [ # ['all_ecal.py'], -# ['draw_cluters.py'], +# ['draw_clusters.py'], # ['EcalEndcapNClusters', 'EcalEndcapPClusters', 'EcalBarrelClusters']], } diff --git a/benchmarks/ecal/scripts/draw_cluters.py b/benchmarks/ecal/scripts/draw_clusters.py similarity index 92% rename from benchmarks/ecal/scripts/draw_cluters.py rename to benchmarks/ecal/scripts/draw_clusters.py index fd03f2dc..f4c81896 100644 --- a/benchmarks/ecal/scripts/draw_cluters.py +++ b/benchmarks/ecal/scripts/draw_clusters.py @@ -25,7 +25,8 @@ def load_root_macros(arg_macros): # read from RDataFrame and flatten a given collection, return pandas dataframe def flatten_collection(rdf, collection, cols=None): if not cols: - cols = [str(c) for c in rdf.GetColumnNames() if str(c).startswith('{}.'.format(collection))] + cols = [str(c) for c in rdf.GetColumnNames() if str(c).startswith('{}.'.format(collection))] + cols += [str(c) for c in rdf.GetColumnNames() if str(c).startswith('{}Info.'.format(collection))] else: cols = ['{}.{}'.format(collection, c) for c in cols] if not cols: @@ -92,10 +93,10 @@ if __name__ == '__main__': ] for coll in [c.strip() for c in args.coll.split(',')]: df = flatten_collection(rdf_rec, coll) - # calculate eta - df.loc[:, '{}.eta'.format(coll)] = -np.log(np.tan(df['{}.polar.theta'.format(coll)].values/2.)) + # Get eta from info table + df[coll + '.eta'] = df[coll + 'Info.eta'] fig, axs = plt.subplots(2, 2, figsize=(12, 8), dpi=160) - ncl = df.groupby('event')['{}.clusterID'.format(coll)].nunique().values + ncl = df.groupby('event')['{}.ID.value'.format(coll)].nunique().values axs[0][0].hist(ncl, weights=np.repeat(1./float(ncl.shape[0]), ncl.shape[0]), bins=np.arange(0, 10), align='mid', ec='k') axs[0][0].set_xlabel('Number of Clusters', fontsize=16) diff --git a/benchmarks/full/options/full_reconstruction.py b/benchmarks/full/options/full_reconstruction.py index ed05a1da..7ff92cfd 100644 --- a/benchmarks/full/options/full_reconstruction.py +++ b/benchmarks/full/options/full_reconstruction.py @@ -122,13 +122,14 @@ ecal_digi = EMCalorimeterDigi("ecal_digi", outputHitCollection="EcalBarrelHitsSimpleDigi") ecal_reco = EMCalReconstruction("ecal_reco", - inputHitCollection="EcalBarrelHitsSimpleDigi", + inputHitCollection=ecal_digi.outputHitCollection, outputHitCollection="EcalBarrelHitsSimpleReco", minModuleEdep=0.0*units.MeV) simple_cluster = SimpleClustering("simple_cluster", - inputHitCollection="EcalBarrelHitsSimpleReco", - outputClusters="SimpleClusters", + inputHitCollection=ecal_reco.outputHitCollection, + outputClusterCollection="SimpleClusters", + outputProtoClusterCollection="SimpleProtoClusters", minModuleEdep=1.0*units.MeV, maxDistance=50.0*units.cm) @@ -146,7 +147,7 @@ ce_ecal_digi = CalHitDigi("ce_ecal_digi", **ce_ecal_daq) ce_ecal_reco = CalHitReco("ce_ecal_reco", - inputHitCollection="EcalEndcapNHitsDigi", + inputHitCollection=ce_ecal_digi.outputHitCollection, outputHitCollection="EcalEndcapNHitsReco", thresholdFactor=4, # 4 sigma cut on pedestal sigma readoutClass="EcalEndcapNHits", @@ -154,8 +155,8 @@ ce_ecal_reco = CalHitReco("ce_ecal_reco", **ce_ecal_daq) ce_ecal_cl = IslandCluster("ce_ecal_cl", - inputHitCollection="EcalEndcapNHitsReco", - outputHitCollection="EcalEndcapNClusterHits", + inputHitCollection=ce_ecal_reco.outputHitCollection, + outputProtoClusterCollection="EcalEndcapNProtoClusters", splitCluster=False, minClusterHitEdep=1.0*units.MeV, # discard low energy hits minClusterCenterEdep=30*units.MeV, @@ -163,7 +164,8 @@ ce_ecal_cl = IslandCluster("ce_ecal_cl", 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", + inputHitCollection=ce_ecal_cl.inputHitCollection, + inputProtoClusterCollection=ce_ecal_cl.outputProtoClusterCollection, outputClusterCollection="EcalEndcapNClusters", samplingFraction=0.998, # this accounts for a small fraction of leakage logWeightBase=4.6) @@ -181,29 +183,31 @@ ci_ecal_digi = CalHitDigi("ci_ecal_digi", **ci_ecal_daq) ci_ecal_reco = CalHitReco("ci_ecal_reco", - inputHitCollection="EcalEndcapPHitsDigi", + inputHitCollection=ci_ecal_digi.outputHitCollection, outputHitCollection="EcalEndcapPHitsReco", thresholdFactor=5.0, **ci_ecal_daq) # merge hits in different layer (projection to local x-y plane) ci_ecal_merger = CalHitsMerger("ci_ecal_merger", - inputHitCollection="EcalEndcapPHitsReco", + inputHitCollection=ci_ecal_reco.outputHitCollection, outputHitCollection="EcalEndcapPHitsRecoXY", fields=["layer", "slice"], fieldRefNumbers=[1, 0], readoutClass="EcalEndcapPHits") ci_ecal_cl = IslandCluster("ci_ecal_cl", - inputHitCollection="EcalEndcapPHitsRecoXY", - outputHitCollection="EcalEndcapPClusterHits", + inputHitCollection=ci_ecal_merger.outputHitCollection, + outputProtoClusterCollection="EcalEndcapPProtoClusters", splitCluster=False, minClusterCenterEdep=10.*units.MeV, localDistXY=[10*units.mm, 10*units.mm]) ci_ecal_clreco = RecoCoG("ci_ecal_clreco", - inputHitCollection="EcalEndcapPClusterHits", + inputHitCollection=ci_ecal_cl.inputHitCollection, + inputProtoClusterCollection=ci_ecal_cl.outputProtoClusterCollection, outputClusterCollection="EcalEndcapPClusters", + outputInfoCollection="EcalEndcapPClustersInfo", logWeightBase=6.2, samplingFraction=ci_ecal_sf) @@ -221,7 +225,7 @@ cb_ecal_digi = CalHitDigi("cb_ecal_digi", **cb_ecal_daq) cb_ecal_reco = ImCalPixelReco("cb_ecal_reco", - inputHitCollection="EcalBarrelHitsDigi", + inputHitCollection=cb_ecal_digi.outputHitCollection, outputHitCollection="EcalBarrelHitsReco", thresholdFactor=3, # about 20 keV readoutClass="EcalBarrelHits", # readout class @@ -230,8 +234,8 @@ cb_ecal_reco = ImCalPixelReco("cb_ecal_reco", **cb_ecal_daq) cb_ecal_cl = ImagingCluster("cb_ecal_cl", - inputHitCollection="EcalBarrelHitsReco", - outputHitCollection="EcalBarrelClusterHits", + inputHitCollection=cb_ecal_reco.outputHitCollection, + outputProtoClusterCollection="EcalBarrelProtoClusters", localDistXY=[2.*units.mm, 2*units.mm], # same layer layerDistEtaPhi=[10*units.mrad, 10*units.mrad], # adjacent layer neighbourLayersRange=2, # id diff for adjacent layer @@ -239,8 +243,10 @@ cb_ecal_cl = ImagingCluster("cb_ecal_cl", cb_ecal_clreco = ImagingClusterReco("cb_ecal_clreco", samplingFraction=cb_ecal_sf, - inputHitCollection="EcalBarrelClusterHits", + inputHitCollection=cb_ecal_cl.inputHitCollection, + inputProtoClusterCollection=cb_ecal_cl.outputProtoClusterCollection, outputClusterCollection="EcalBarrelClusters", + outputInfoCollection="EcalBarrelClustersInfo", outputLayerCollection="EcalBarrelLayers") # Central ECAL SciFi @@ -256,7 +262,7 @@ scfi_barrel_digi = CalHitDigi("scfi_barrel_digi", **scfi_barrel_daq) scfi_barrel_reco = CalHitReco("scfi_barrel_reco", - inputHitCollection="EcalBarrelScFiHitsDigi", + inputHitCollection=scfi_barrel_digi.outputHitCollection, outputHitCollection="EcalBarrelScFiHitsReco", thresholdFactor=5.0, readoutClass="EcalBarrelScFiHits", @@ -267,22 +273,24 @@ scfi_barrel_reco = CalHitReco("scfi_barrel_reco", # merge hits in different layer (projection to local x-y plane) scfi_barrel_merger = CalHitsMerger("scfi_barrel_merger", - inputHitCollection="EcalBarrelScFiHitsReco", + inputHitCollection=scfi_barrel_reco.outputHitCollection, outputHitCollection="EcalBarrelScFiGridReco", fields=["fiber"], fieldRefNumbers=[1], readoutClass="EcalBarrelScFiHits") scfi_barrel_cl = IslandCluster("scfi_barrel_cl", - inputHitCollection="EcalBarrelScFiGridReco", - outputHitCollection="EcalBarrelScFiClusterHits", + inputHitCollection=scfi_barrel_merger.outputHitCollection, + outputProtoClusterCollection="EcalBarrelScFiProtoClusters", splitCluster=False, minClusterCenterEdep=10.*MeV, localDistXZ=[30*mm, 30*mm]) scfi_barrel_clreco = RecoCoG("scfi_barrel_clreco", - inputHitCollection="EcalBarrelScFiClusterHits", + inputHitCollection=scfi_barrel_cl.inputHitCollection, + inputProtoClusterCollection=scfi_barrel_cl.outputProtoClusterCollection, outputClusterCollection="EcalBarrelScFiClusters", + outputInfoCollection="EcalBarrelScFiClustersInfo", logWeightBase=6.2, samplingFraction= scifi_barrel_sf) @@ -299,7 +307,7 @@ cb_hcal_digi = CalHitDigi("cb_hcal_digi", **cb_hcal_daq) cb_hcal_reco = CalHitReco("cb_hcal_reco", - inputHitCollection="HcalBarrelHitsDigi", + inputHitCollection=cb_hcal_digi.outputHitCollection, outputHitCollection="HcalBarrelHitsReco", thresholdFactor=5.0, readoutClass="HcalBarrelHits", @@ -308,22 +316,24 @@ cb_hcal_reco = CalHitReco("cb_hcal_reco", **cb_hcal_daq) cb_hcal_merger = CalHitsMerger("cb_hcal_merger", - inputHitCollection="HcalBarrelHitsReco", + inputHitCollection=cb_hcal_reco.outputHitCollection, outputHitCollection="HcalBarrelHitsRecoXY", readoutClass="HcalBarrelHits", fields=["layer", "slice"], fieldRefNumbers=[1, 0]) cb_hcal_cl = IslandCluster("cb_hcal_cl", - inputHitCollection="HcalBarrelHitsRecoXY", - outputHitCollection="HcalBarrelClusterHits", + inputHitCollection=cb_hcal_merger.outputHitCollection, + outputProtoClusterCollection="HcalBarrelProtoClusters", splitCluster=False, minClusterCenterEdep=30.*units.MeV, localDistXY=[15.*units.cm, 15.*units.cm]) cb_hcal_clreco = RecoCoG("cb_hcal_clreco", - inputHitCollection="HcalBarrelClusterHits", + inputHitCollection=cb_hcal_cl.inputHitCollection, + inputProtoClusterCollection=cb_hcal_cl.outputProtoClusterCollection, outputClusterCollection="HcalBarrelClusters", + outputInfoCollection="HcalBarrelClustersInfo", logWeightBase=6.2, samplingFraction=cb_hcal_sf) @@ -340,28 +350,30 @@ ci_hcal_digi = CalHitDigi("ci_hcal_digi", **ci_hcal_daq) ci_hcal_reco = CalHitReco("ci_hcal_reco", - inputHitCollection="HcalHadronEndcapHitsDigi", + inputHitCollection=ci_hcal_digi.outputHitCollection, outputHitCollection="HcalHadronEndcapHitsReco", thresholdFactor=5.0, **ci_hcal_daq) ci_hcal_merger = CalHitsMerger("ci_hcal_merger", - inputHitCollection="HcalHadronEndcapHitsReco", + inputHitCollection=ci_hcal_reco.outputHitCollection, outputHitCollection="HcalHadronEndcapHitsRecoXY", readoutClass="HcalHadronEndcapHits", fields=["layer", "slice"], fieldRefNumbers=[1, 0]) ci_hcal_cl = IslandCluster("ci_hcal_cl", - inputHitCollection="HcalHadronEndcapHitsRecoXY", - outputHitCollection="HcalHadronEndcapClusterHits", + inputHitCollection=ci_hcal_merger.outputHitCollection, + outputProtoClusterCollection="HcalHadronEndcapProtoClusters", splitCluster=False, minClusterCenterEdep=30.*units.MeV, localDistXY=[15.*units.cm, 15.*units.cm]) ci_hcal_clreco = RecoCoG("ci_hcal_clreco", - inputHitCollection="HcalHadronEndcapClusterHits", + inputHitCollection=ci_hcal_cl.inputHitCollection, + inputProtoClusterCollection=ci_hcal_cl.outputProtoClusterCollection, outputClusterCollection="HcalHadronEndcapClusters", + outputInfoCollection="HcalHadronEndcapClustersInfo", logWeightBase=6.2, samplingFraction=ci_hcal_sf) @@ -378,28 +390,30 @@ ce_hcal_digi = CalHitDigi("ce_hcal_digi", **ce_hcal_daq) ce_hcal_reco = CalHitReco("ce_hcal_reco", - inputHitCollection="HcalElectronEndcapHitsDigi", + inputHitCollection=ce_hcal_digi.outputHitCollection, outputHitCollection="HcalElectronEndcapHitsReco", thresholdFactor=5.0, **ce_hcal_daq) ce_hcal_merger = CalHitsMerger("ce_hcal_merger", - inputHitCollection="HcalElectronEndcapHitsReco", + inputHitCollection=ce_hcal_reco.outputHitCollection, outputHitCollection="HcalElectronEndcapHitsRecoXY", readoutClass="HcalElectronEndcapHits", fields=["layer", "slice"], fieldRefNumbers=[1, 0]) ce_hcal_cl = IslandCluster("ce_hcal_cl", - inputHitCollection="HcalElectronEndcapHitsRecoXY", - outputHitCollection="HcalElectronEndcapClusterHits", + inputHitCollection=ce_hcal_merger.outputHitCollection, + outputProtoClusterCollection="HcalElectronEndcapProtoClusters", splitCluster=False, minClusterCenterEdep=30.*units.MeV, localDistXY=[15.*units.cm, 15.*units.cm]) ce_hcal_clreco = RecoCoG("ce_hcal_clreco", - inputHitCollection="HcalElectronEndcapClusterHits", + inputHitCollection=ce_hcal_cl.inputHitCollection, + inputProtoClusterCollection=ce_hcal_cl.outputProtoClusterCollection, outputClusterCollection="HcalElectronEndcapClusters", + outputInfoCollection="HcalElectronEndcapClustersInfo", logWeightBase=6.2, samplingFraction=ce_hcal_sf) diff --git a/benchmarks/imaging_ecal/options/hybrid_cluster.py b/benchmarks/imaging_ecal/options/hybrid_cluster.py index e89cef98..5b2e1ead 100644 --- a/benchmarks/imaging_ecal/options/hybrid_cluster.py +++ b/benchmarks/imaging_ecal/options/hybrid_cluster.py @@ -65,7 +65,7 @@ imcaldigi = CalHitDigi('imcal_digi', **imcaldaq) imcalreco = ImagingPixelReco('imcal_reco', # OutputLevel=DEBUG, - inputHitCollection='DigiEcalBarrelHits', + inputHitCollection=imcaldigi.outputHitCollection, outputHitCollection='RecoEcalBarrelHits', readoutClass='EcalBarrelHits', layerField='layer', @@ -74,17 +74,19 @@ imcalreco = ImagingPixelReco('imcal_reco', imcalcluster = ImagingTopoCluster('imcal_cluster', # OutputLevel=DEBUG, - inputHitCollection='RecoEcalBarrelHits', - outputHitCollection='EcalBarrelClusterHits', + inputHitCollection=imcalreco.outputHitCollection, + outputProtoClusterCollection='EcalBarrelProtoClusters', localDistXY=[2.*mm, 2*mm], layerDistEtaPhi=[10*mrad, 10*mrad], neighbourLayersRange=2, sectorDist=3.*cm) clusterreco = ImagingClusterReco('imcal_clreco', # OutputLevel=DEBUG, - inputHitCollection='EcalBarrelClusterHits', + inputHitCollection=imcalcluster.inputHitCollection, + inputProtoClusterCollection=imcalcluster.outputProtoClusterCollection, outputLayerCollection='EcalBarrelClustersLayers', outputClusterCollection='EcalBarrelClusters', + outputInfoCollection='EcalBarrelClustersInfo', samplingFraction=kwargs['img_sf']) # scfi layers @@ -100,7 +102,7 @@ scfi_barrel_digi = CalHitDigi("scfi_barrel_digi", **scfi_barrel_daq) scfi_barrel_reco = CalHitReco("scfi_barrel_reco", - inputHitCollection="EcalBarrelScFiHitsDigi", + inputHitCollection=scfi_barrel_digi.outputHitCollection, outputHitCollection="EcalBarrelScFiHitsReco", thresholdFactor=5.0, readoutClass="EcalBarrelScFiHits", @@ -112,7 +114,7 @@ scfi_barrel_reco = CalHitReco("scfi_barrel_reco", # merge hits in different layer (projection to local x-y plane) scfi_barrel_merger = CalHitsMerger("scfi_barrel_merger", # OutputLevel=DEBUG, - inputHitCollection="EcalBarrelScFiHitsReco", + inputHitCollection=scfi_barrel_reco.outputHitCollection, outputHitCollection="EcalBarrelScFiGridReco", fields=["fiber"], fieldRefNumbers=[1], @@ -120,16 +122,18 @@ scfi_barrel_merger = CalHitsMerger("scfi_barrel_merger", scfi_barrel_cl = IslandCluster("scfi_barrel_cl", # OutputLevel=DEBUG, - inputHitCollection="EcalBarrelScFiGridReco", - outputHitCollection="EcalBarrelScFiClusterHits", + inputHitCollection=scfi_barrel_merger.outputHitCollection, + outputProtoClusterCollection="EcalBarrelScFiProtoClusters", splitCluster=False, minClusterCenterEdep=10.*MeV, localDistXZ=[30*mm, 30*mm], sectorDist=5.*cm) scfi_barrel_clreco = RecoCoG("scfi_barrel_clreco", - inputHitCollection="EcalBarrelScFiClusterHits", + inputHitCollection=scfi_barrel_cl.inputHitCollection, + inputProtoClusterCollection=scfi_barrel_cl.outputProtoClusterCollection, outputClusterCollection="EcalBarrelScFiClusters", + outputInfoCollection="EcalBarrelScFiClustersInfo", logWeightBase=6.2, samplingFraction=kwargs['scfi_sf']) diff --git a/benchmarks/imaging_ecal/options/imaging_2dcluster.py b/benchmarks/imaging_ecal/options/imaging_2dcluster.py index e58e962f..9baee078 100644 --- a/benchmarks/imaging_ecal/options/imaging_2dcluster.py +++ b/benchmarks/imaging_ecal/options/imaging_2dcluster.py @@ -64,7 +64,7 @@ imcal_barrel_merger = CalHitsProj('imcal_barrel_merger', imcal_barrel_cl = IslandCluster('imcal_barrel_cl', # OutputLevel=DEBUG, inputHitCollection=imcal_barrel_merger.outputHitCollection, - outputHitCollection='EcalBarrelClusterHits', + outputProtoClusterCollection='EcalBarrelProtoClusters', minClusterCenterEdep=1.*MeV, minClusterHitEdep=0.1*MeV, splitCluster=True, @@ -72,8 +72,10 @@ imcal_barrel_cl = IslandCluster('imcal_barrel_cl', imcal_barrel_clreco = RecoCoG('imcal_barrel_clreco', # OutputLevel=DEBUG, - inputHitCollection=imcal_barrel_cl.outputHitCollection, + inputHitCollection=imcal_barrel_cl.inputHitCollection, + inputProtoClusterCollection = imcal_barrel_cl.outputProtoClusterCollection, outputClusterCollection='EcalBarrelClusters', + outputInfoCollection='EcalBarrelClustersInfo', logWeightBase=6.2, samplingFraction=kwargs['sf']) diff --git a/benchmarks/imaging_ecal/options/imaging_topocluster.py b/benchmarks/imaging_ecal/options/imaging_topocluster.py index 572544ef..9c715fcc 100644 --- a/benchmarks/imaging_ecal/options/imaging_topocluster.py +++ b/benchmarks/imaging_ecal/options/imaging_topocluster.py @@ -59,7 +59,7 @@ imcaldigi = CalorimeterHitDigi("imcal_digi", **daq_setting) imcalreco = ImagingPixelReco("imcal_reco", # OutputLevel=DEBUG, - inputHitCollection="DigiEcalBarrelHits", + inputHitCollection=imcaldigi.outputHitCollection, outputHitCollection="RecoEcalBarrelHits", readoutClass="EcalBarrelHits", layerField="layer", @@ -68,8 +68,8 @@ imcalreco = ImagingPixelReco("imcal_reco", imcalcluster = ImagingTopoCluster("imcal_cluster", # OutputLevel=DEBUG, - inputHitCollection="RecoEcalBarrelHits", - outputHitCollection="EcalBarrelClusterHits", + inputHitCollection=imcalreco.outputHitCollection, + outputProtoClusterCollection="EcalBarrelProtoClusters", localDistXY=[2.*units.mm, 2*units.mm], layerDistEtaPhi=[10*units.mrad, 10*units.mrad], neighbourLayersRange=2, @@ -77,9 +77,11 @@ imcalcluster = ImagingTopoCluster("imcal_cluster", minClusterNhits=5) clusterreco = ImagingClusterReco("imcal_clreco", # OutputLevel=DEBUG, - inputHitCollection="EcalBarrelClusterHits", + inputHitCollection=imcalcluster.inputHitCollection, + inputProtoClusterCollection=imcalcluster.outputProtoClusterCollection, outputLayerCollection="EcalBarrelClustersLayers", outputClusterCollection="EcalBarrelClusters", + outputInfoCollection="EcalBarrelClustersInfo", samplingFraction=sf) out.outputCommands = ["keep *"] diff --git a/benchmarks/imaging_ecal/options/scfi_cluster.py b/benchmarks/imaging_ecal/options/scfi_cluster.py index 873c3b89..3926326b 100644 --- a/benchmarks/imaging_ecal/options/scfi_cluster.py +++ b/benchmarks/imaging_ecal/options/scfi_cluster.py @@ -47,7 +47,7 @@ scfi_barrel_digi = CalHitDigi("scfi_barrel_digi", **scfi_barrel_daq) scfi_barrel_reco = CalHitReco("scfi_barrel_reco", - inputHitCollection="EcalBarrelScFiHitsDigi", + inputHitCollection=scfi_barrel_digi.outputHitCollection, outputHitCollection="EcalBarrelScFiHitsReco", thresholdFactor=5.0, readoutClass="EcalBarrelScFiHits", @@ -59,7 +59,7 @@ scfi_barrel_reco = CalHitReco("scfi_barrel_reco", # merge hits in different layer (projection to local x-y plane) scfi_barrel_merger = CalHitsMerger("scfi_barrel_merger", # OutputLevel=DEBUG, - inputHitCollection="EcalBarrelScFiHitsReco", + inputHitCollection=scfi_barrel_reco.outputHitCollection, outputHitCollection="EcalBarrelScFiGridReco", fields=["fiber"], fieldRefNumbers=[1], @@ -67,15 +67,17 @@ scfi_barrel_merger = CalHitsMerger("scfi_barrel_merger", scfi_barrel_cl = IslandCluster("scfi_barrel_cl", # OutputLevel=DEBUG, - inputHitCollection="EcalBarrelScFiGridReco", - outputHitCollection="EcalBarrelScFiClusterHits", + inputHitCollection=scfi_barrel_reco.outputHitCollection, + outputProtoClusterCollection="EcalBarrelScFiProtoClusters", splitCluster=False, minClusterCenterEdep=10.*MeV, localDistXZ=[30*mm, 30*mm]) scfi_barrel_clreco = RecoCoG("scfi_barrel_clreco", - inputHitCollection="EcalBarrelScFiClusterHits", + inputHitCollection=scfi_barrel_cl.inputHitCollection, + inputProtoClusterCollection=scfi_barrel.outputProtoClusterCollection, outputClusterCollection="EcalBarrelScFiClusters", + outputInfoCollection="EcalBarrelScFiClustersInfo", logWeightBase=6.2, samplingFraction=kwargs['sf']) diff --git a/benchmarks/imaging_ecal/scripts/draw_cluster.py b/benchmarks/imaging_ecal/scripts/draw_cluster.py index 88f338d0..3c896f2c 100644 --- a/benchmarks/imaging_ecal/scripts/draw_cluster.py +++ b/benchmarks/imaging_ecal/scripts/draw_cluster.py @@ -26,8 +26,8 @@ import sys # note z and x axes are switched def draw_hits3d(axis, data, cmap, units=('mm', 'mm', 'mm', 'MeV'), fontsize=24, **kwargs): # normalize energy to get colors - x, y, z, edep = np.transpose(data) - cvals = edep - min(edep) / (max(edep) - min(edep)) + x, y, z, energy = np.transpose(data) + cvals = energy - min(energy) / (max(energy) - min(energy)) cvals[np.isnan(cvals)] = 1.0 colors = cmap(cvals) @@ -37,7 +37,7 @@ def draw_hits3d(axis, data, cmap, units=('mm', 'mm', 'mm', 'MeV'), fontsize=24, axis.set_zlabel('x ({})'.format(units[2]), fontsize=fontsize + 2, labelpad=fontsize) axis.set_ylabel('y ({})'.format(units[1]), fontsize=fontsize + 2, labelpad=fontsize) axis.set_xlabel('z ({})'.format(units[0]), fontsize=fontsize + 2, labelpad=fontsize) - cb = plt.colorbar(cm.ScalarMappable(norm=matplotlib.colors.Normalize(vmin=min(edep), vmax=max(edep)), cmap=cmap), + cb = plt.colorbar(cm.ScalarMappable(norm=matplotlib.colors.Normalize(vmin=min(energy), vmax=max(energy)), cmap=cmap), ax=axis, shrink=0.85) cb.ax.tick_params(labelsize=fontsize) cb.ax.get_yaxis().labelpad = fontsize @@ -158,7 +158,7 @@ if __name__ == '__main__': # filter out outliers dfc = df[(df['eta'] <= df['eta'].quantile(0.95)) & (df['eta'] >= df['eta'].quantile(0.05)) & (df['phi'] <= df['phi'].quantile(0.95)) & (df['phi'] >= df['phi'].quantile(0.05))] - vec = np.average(dfc[['x', 'y', 'z']].values, axis=0, weights=dfc['edep'].values) + vec = np.average(dfc[['x', 'y', 'z']].values, axis=0, weights=dfc['energy'].values) vec = vec/np.linalg.norm(vec) # particle line from (0, 0, 0) to the inner Ecal surface @@ -174,7 +174,7 @@ if __name__ == '__main__': # TODO need to implement motion of particles in a field # ax.plot(*pline[[2, 1]], '--', zs=pline[0], color='green') # draw hits - draw_hits3d(ax, df[['x', 'y', 'z', 'edep']].values, cmap, s=5.0) + draw_hits3d(ax, df[['x', 'y', 'z', 'energy']].values, cmap, s=5.0) draw_cylinder3d(ax, rmin, length, rstride=10, cstride=10, color='royalblue') draw_cylinder3d(ax, rmin + thickness, length, rstride=10, cstride=10, color='forestgreen') ax.set_zlim(-(rmin + thickness), rmin + thickness) @@ -191,7 +191,7 @@ if __name__ == '__main__': # TODO need to implement motion of particles in a field # ax.plot(*pline[[2, 1]], '--', zs=pline[0], color='green') # draw hits - draw_hits3d(ax, df[['x', 'y', 'z', 'edep']].values, cmap, s=20.0) + draw_hits3d(ax, df[['x', 'y', 'z', 'energy']].values, cmap, s=20.0) # draw_track_fit(ax, df, stop_layer=args.stop, # scat_kw=dict(color='k', s=50.0), line_kw=dict(linestyle=':', color='k', lw=3)) # view range @@ -213,7 +213,7 @@ if __name__ == '__main__': eta_rg = np.resize(-np.log(np.tan(vecp[0]/1000./2.)), 2) + np.asarray([-args.topo_range, args.topo_range])/1000. fig, axs = plt.subplots(1, 2, figsize=(13, 12), dpi=160, gridspec_kw={'wspace':0., 'width_ratios': [12, 1]}) - ax, sm = draw_heatmap(axs[0], df['eta'].values, df['phi'].values, weights=df['edep'].values, + ax, sm = draw_heatmap(axs[0], df['eta'].values, df['phi'].values, weights=df['energy'].values, bins=(np.arange(*eta_rg, step=args.topo_size/1000.), np.arange(*phi_rg, step=args.topo_size)), cmap=cmap, cmin=0., pc_kw=dict(alpha=0.8, edgecolor='k')) diff --git a/benchmarks/imaging_ecal/scripts/draw_cluster_layers.py b/benchmarks/imaging_ecal/scripts/draw_cluster_layers.py index 7baafebb..13a8eb1b 100644 --- a/benchmarks/imaging_ecal/scripts/draw_cluster_layers.py +++ b/benchmarks/imaging_ecal/scripts/draw_cluster_layers.py @@ -110,7 +110,7 @@ if __name__ == '__main__': else: dfc = df[(df['eta'] <= df['eta'].quantile(0.95)) & (df['eta'] >= df['eta'].quantile(0.05)) & (df['phi'] <= df['phi'].quantile(0.95)) & (df['phi'] >= df['phi'].quantile(0.05))] - vec = np.average(dfc[['x', 'y', 'z']].values, axis=0, weights=dfc['edep'].values) + vec = np.average(dfc[['x', 'y', 'z']].values, axis=0, weights=dfc['energy'].values) vec = vec/np.linalg.norm(vec) # particle line from (0, 0, 0) to the inner Ecal surface @@ -133,7 +133,7 @@ if __name__ == '__main__': for i in np.arange(1, df['layer'].max() + 1, dtype=int): data = df[df['layer'] == i] fig, axs = plt.subplots(1, 2, figsize=(17, 16), dpi=160, gridspec_kw={'wspace':0., 'width_ratios': [16, 1]}) - ax, sm = draw_heatmap(axs[0], data['eta'].values, data['phi'].values, weights=data['edep'].values, + ax, sm = draw_heatmap(axs[0], data['eta'].values, data['phi'].values, weights=data['energy'].values, bins=(np.arange(*eta_rg, step=args.topo_size/1000.), np.arange(*phi_rg, step=args.topo_size)), cmap=cmap, vmin=0.0, vmax=1.0, pc_kw=dict(alpha=0.8, edgecolor='k')) ax.set_ylabel(r'$\phi$ (mrad)', fontsize=28) @@ -170,7 +170,7 @@ if __name__ == '__main__': data = df[df['layer'] == i] fig, axs = plt.subplots(1, 2, figsize=(17, 16), dpi=160, gridspec_kw={'wspace':0., 'width_ratios': [16, 1]}) ax = axs[0] - colors = cmap(data['edep']/1.0) + colors = cmap(data['energy']/1.0) ax.scatter(data['eta'].values, data['phi'].values, c=colors, marker='s', s=15.0) ax.set_xlim(*eta_rg) ax.set_ylim(*phi_rg) diff --git a/benchmarks/imaging_ecal/scripts/energy_profile.py b/benchmarks/imaging_ecal/scripts/energy_profile.py index fff925ce..a8406f0c 100644 --- a/benchmarks/imaging_ecal/scripts/energy_profile.py +++ b/benchmarks/imaging_ecal/scripts/energy_profile.py @@ -20,10 +20,10 @@ import sys from utils import * -def find_start_layer(grp, min_edep=0.5): - ids, edeps = grp.sort_values('layer')[['layer', 'edep']].values.T - edeps = np.cumsum(edeps) - return min(ids[edeps > min_edep]) if sum(edeps > min_edep) > 0 else -1 +def find_start_layer(grp, min_energy=0.5): + ids, energys = grp.sort_values('layer')[['layer', 'energy']].values.T + energys = np.cumsum(energys) + return min(ids[energys > min_energy]) if sum(energys > min_energy) > 0 else -1 # execute this script @@ -44,9 +44,9 @@ if __name__ == '__main__': load_root_macros(args.macros) # prepare data dfe = get_layers_data(args.file, branch=args.branch) - dfe.loc[:, 'total_edep'] = dfe.groupby(['event', 'cluster'])['edep'].transform('sum') - # dfe.loc[:, 'efrac'] = dfe['edep']/dfe['total_edep'] - dfe.loc[:, 'efrac'] = dfe['edep']/(args.energy*1000.) + dfe.loc[:, 'total_energy'] = dfe.groupby(['event', 'cluster'])['energy'].transform('sum') + # dfe.loc[:, 'efrac'] = dfe['energy']/dfe['total_energy'] + dfe.loc[:, 'efrac'] = dfe['energy']/(args.energy*1000.) # dfe = dfe[dfe['cluster'] == 0] os.makedirs(args.outdir, exist_ok=True) @@ -54,7 +54,7 @@ if __name__ == '__main__': slayer = dfe.groupby('event').apply(lambda x: find_start_layer(x, 1.0)).astype(int) dfe = dfe.merge(slayer.to_frame(name='slayer'), on='event') # dfe.loc[:, 'eff_layer'] = dfe['layer'] - dfe['slayer'] - # prof = dfe[dfe['eff_layer'] > 0].groupby('eff_layer')['edep'].describe() + # prof = dfe[dfe['eff_layer'] > 0].groupby('eff_layer')['energy'].describe() prof = dfe.groupby('layer')['efrac'].describe() # print(prof['mean']) @@ -73,7 +73,7 @@ if __name__ == '__main__': ax.set_ylabel('Normalized Counts', fontsize=26) ax.text(*bpos, 'Mininum Edep\n' + '{:.1f} MeV'.format(1.0), transform=ax.transAxes, fontsize=26, verticalalignment='top', bbox=bprops) - fig.savefig(os.path.join(args.outdir, 'edep_start_{}_{}.png'.format(args.type, int(args.energy)))) + fig.savefig(os.path.join(args.outdir, 'energy_start_{}_{}.png'.format(args.type, int(args.energy)))) fig, ax = plt.subplots(figsize=(16, 9), dpi=160) @@ -90,7 +90,7 @@ if __name__ == '__main__': ax.set_ylabel('Energy Deposit Percentage', fontsize=26) fig.savefig(os.path.join(args.outdir, 'efrac_{}_{}.png'.format(args.type, int(args.energy)))) - # edep fraction by layers + # energy fraction by layers layers = np.asarray([ [1, 5, 8,], [10, 15, 20], @@ -111,7 +111,7 @@ if __name__ == '__main__': ax.grid(linestyle=':', which='both') # ax.set_xlabel('Energy Deposit (MeV)', fontsize=26) # ax.set_ylabel('Normalized Counts', fontsize=26) - mean, std = data.describe().loc[['mean', 'std'], 'edep'].values + mean, std = data.describe().loc[['mean', 'std'], 'energy'].values label = 'Layer {}'.format(layer) # + '\n' + r'$\mu = {:.3f}$ MeV'.format(mean) # + '\n' + r'$\sigma = {:.3f}$ MeV'.format(std) diff --git a/benchmarks/imaging_ecal/scripts/utils.py b/benchmarks/imaging_ecal/scripts/utils.py index da9954ac..d3847618 100644 --- a/benchmarks/imaging_ecal/scripts/utils.py +++ b/benchmarks/imaging_ecal/scripts/utils.py @@ -98,10 +98,12 @@ def get_hits_data(path, evnums=None, branch='RecoEcalBarrelHits'): events.GetEntry(iev) for hit in getattr(events, branch): - dbuf[idb] = (iev, hit.clusterID, hit.layerID, hit.position.x, hit.position.y, hit.position.z, hit.edep*1000.) + dbuf[idb] = (iev, hit.ID.value, hit.layer, hit.position.x, hit.position.y, + hit.position.z, hit.energy*1000.) idb += 1 - return pd.DataFrame(data=dbuf[:idb], columns=['event', 'cluster', 'layer', 'x', 'y', 'z', 'edep']) + return pd.DataFrame(data=dbuf[:idb], columns=['event', 'cluster', 'layer', 'x', 'y', + 'z', 'energy']) # read layers data from root file @@ -122,11 +124,13 @@ def get_layers_data(path, evnums=None, branch="EcalBarrelClustersLayers"): events.GetEntry(iev) for layer in getattr(events, branch): - dbuf[idb] = (iev, layer.clusterID, layer.layerID, - layer.position.x, layer.position.y, layer.position.z, layer.edep*1000.) + dbuf[idb] = (iev, layer.clusterID.value, layer.layer, + layer.position.x, layer.position.y, layer.position.z, + layer.energy*1000.) idb += 1 - return pd.DataFrame(data=dbuf[:idb], columns=['event', 'cluster', 'layer', 'x', 'y', 'z', 'edep']) + return pd.DataFrame(data=dbuf[:idb], columns=['event', 'cluster', 'layer', 'x', 'y', + 'z', 'energy']) # read clusters data from root file @@ -147,10 +151,10 @@ def get_clusters_data(path, evnums=None, branch='EcalBarrelClustersReco'): events.GetEntry(iev) for k, cl in enumerate(getattr(events, branch)): - dbuf[idb] = (iev, k, cl.nhits, cl.edep*1000., cl.cl_theta, cl.cl_phi) + dbuf[idb] = (iev, k, cl.nhits, cl.energy*1000., cl.cl_theta, cl.cl_phi) idb += 1 - return pd.DataFrame(data=dbuf[:idb], columns=['event', 'cluster', 'nhits', 'edep', 'cl_theta', 'cl_phi']) + return pd.DataFrame(data=dbuf[:idb], columns=['event', 'cluster', 'nhits', 'energy', 'cl_theta', 'cl_phi']) def compact_constants(path, names): diff --git a/benchmarks/imaging_shower_ML/options/imaging_ml_data.py b/benchmarks/imaging_shower_ML/options/imaging_ml_data.py index 1bdc3f2f..7cf23956 100644 --- a/benchmarks/imaging_shower_ML/options/imaging_ml_data.py +++ b/benchmarks/imaging_shower_ML/options/imaging_ml_data.py @@ -46,7 +46,7 @@ imcaldigi = CalorimeterHitDigi('imcal_digi', dynamicRangeADC=3*units.MeV, pedestalSigma=40) imcalreco = ImagingPixelReco('imcal_reco', - inputHitCollection='EcalBarrelHitsDigi', + inputHitCollection=imcaldigi.outputHitCollection, outputHitCollection='EcalBarrelHitsReco', dynamicRangeADC=3.*units.MeV, pedestalSigma=40, @@ -56,7 +56,7 @@ imcalreco = ImagingPixelReco('imcal_reco', sectorField='module') imcaldata = ImagingPixelMerger('imcal_merger', # OutputLevel=DEBUG, - inputHitCollection='EcalBarrelHitsReco', + inputHitCollection=imcalreco.outputHitCollection, outputHitCollection='EcalBarrelHitsML', etaSize=0.001, phiSize=0.001, diff --git a/benchmarks/tracking/options/imagingcluster_seeded_tracking.py b/benchmarks/tracking/options/imagingcluster_seeded_tracking.py index f5a47fc0..9f691670 100644 --- a/benchmarks/tracking/options/imagingcluster_seeded_tracking.py +++ b/benchmarks/tracking/options/imagingcluster_seeded_tracking.py @@ -110,7 +110,7 @@ cb_ecal_digi = CalHitDigi("cb_ecal_digi", **cb_ecal_daq) cb_ecal_reco = ImCalPixelReco("cb_ecal_reco", - inputHitCollection="EcalBarrelHitsDigi", + inputHitCollection=cb_ecal_digi.outputHitCollection, outputHitCollection="EcalBarrelHitsReco", thresholdFactor=3, # about 20 keV readoutClass="EcalBarrelHits", # readout class @@ -119,8 +119,8 @@ cb_ecal_reco = ImCalPixelReco("cb_ecal_reco", **cb_ecal_daq) cb_ecal_cl = ImagingCluster("cb_ecal_cl", - inputHitCollection="EcalBarrelHitsReco", - outputHitCollection="EcalBarrelClusterHits", + inputHitCollection=cb_ecal_reco.outputHitCollection, + outputProtoClusterCollection="EcalBarrelProtoClusters", localDistXY=[2.*mm, 2*mm], # same layer layerDistEtaPhi=[10*mrad, 10*mrad], # adjacent layer neighbourLayersRange=2, # id diff for adjacent layer @@ -128,8 +128,10 @@ cb_ecal_cl = ImagingCluster("cb_ecal_cl", cb_ecal_clreco = ImagingClusterReco("cb_ecal_clreco", samplingFraction=cb_ecal_sf, - inputHitCollection="EcalBarrelClusterHits", + inputHitCollection=cb_ecal_cl.inputHitCollection, + inputProtoClusterCollection=cb_ecal_cl.outputProtoClusterCollection, outputClusterCollection="EcalBarrelClusters", + outputInfoCollection="EcalBarrelClustersInfo", outputLayerCollection="EcalBarrelLayers") ## Track param init diff --git a/options/tracker_reconstruction.py b/options/tracker_reconstruction.py index a9f0347d..38dcf0f7 100644 --- a/options/tracker_reconstruction.py +++ b/options/tracker_reconstruction.py @@ -86,13 +86,14 @@ ufsd_digi2 = UFSDTrackerDigi("ufsd_digi2", ecal_reco = EMCalReconstruction("ecal_reco", - inputHitCollection="RawEcalBarrelHits", + inputHitCollection=ecal_digi.outputHitCollection, outputHitCollection="RecEcalBarrelHits", minModuleEdep=0.0*units.MeV) #OutputLevel=DEBUG) simple_cluster = SimpleClustering("simple_cluster", - inputHitCollection="RecEcalBarrelHits", + inputHitCollection=ecal_reco.outputHitCollection, + outputProtoClusters="SimpleProtoClusters", outputClusters="SimpleClusters", minModuleEdep=1.0*units.MeV, maxDistance=50.0*units.cm) -- GitLab