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 402bc11e93861ffb7a1b8ad7b63855cd9947dea0..3f23faa1d75c975fb769cacc2b22e2fa784045b3 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 336c35d4d8fbe201c43ab704aa1595c5a59e2719..d5a890acbba2213adc8d0e923a6b4dc5cdb85804 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)