adjust scripts for crystal-glass endcap ecal

benchmarks/clustering/options/full_cal_clusters.py

@@ -48,8 +48,8 @@ from Configurables import Jug__Reco__ImagingClusterReco as ImagingClusterReco
 # branches needed from simulation root file
 digi_coll = [
 #    "mcparticles2",
-    "CrystalEcalHitsReco",
-    "EcalEndcapHitsReco",
+    "EcalEndcapNHitsReco",
+    "EcalEndcapPHitsReco",
     "EcalBarrelHitsReco",
     "HcalBarrelHitsReco",
     "HcalHadronEndcapHitsReco",
@@ -66,38 +66,38 @@ podout = PodioOutput("out", filename=output_rec)
 # Crystal Endcap Ecal
 ce_ecal_cl = IslandCluster("ce_ecal_cl",
         # OutputLevel=DEBUG,
-        inputHitCollection="CrystalEcalHitsReco",
-        outputClusterCollection="CrystalEcalClusters",
-        splitHitCollection="CrystalEcalHitsSplit",
+        inputHitCollection="EcalEndcapNHitsReco",
+        outputClusterCollection="EcalEndcapNClusters",
+        splitHitCollection="EcalEndcapNHitsSplit",
         splitCluster=False,
         minClusterCenterEdep=30*MeV,
         groupRanges=[2.2*cm, 2.2*cm])
 
 ce_ecal_clreco = RecoCoG("ce_ecal_clreco",
-        clusterCollection="CrystalEcalClusters",
+        clusterCollection="EcalEndcapNClusters",
         samplingFraction=0.998,      # this accounts for a small fraction of leakage
         logWeightBase=4.6)
 
 
 # Endcap Sampling Ecal
 # merge hits in different layer (projection to local x-y plane)
-ce_ecal2_merger = CalHitsMerger("ce_ecal2_merger",
-        inputHitCollection="EcalEndcapHitsReco",
-        outputHitCollection="EcalEndcapHitsRecoXY",
+ci_ecal_merger = CalHitsMerger("ci_ecal_merger",
+        inputHitCollection="EcalEndcapPHitsReco",
+        outputHitCollection="EcalEndcapPHitsRecoXY",
         fields=["layer", "slice"],
         fieldRefNumbers=[1, 0],
-        readoutClass="EcalEndcapHits")
+        readoutClass="EcalEndcapPHits")
 
-ce_ecal2_cl = IslandCluster("ce_ecal2_cl",
-        inputHitCollection="EcalEndcapHitsRecoXY",
-        outputClusterCollection="EcalEndcapClusters",
-        splitHitCollection="EcalEndcapHitsSplit",
+ci_ecal_cl = IslandCluster("ci_ecal_cl",
+        inputHitCollection="EcalEndcapPHitsRecoXY",
+        outputClusterCollection="EcalEndcapPClusters",
+        splitHitCollection="EcalEndcapPHitsSplit",
         splitCluster=False,
         minClusterCenterEdep=30.*MeV,
         groupRanges=[5*mm, 5*mm])
 
-ce_ecal2_clreco = RecoCoG("ce_ecal2_clreco",
-        clusterCollection="EcalEndcapClusters",
+ci_ecal_clreco = RecoCoG("ci_ecal_clreco",
+        clusterCollection="EcalEndcapPClusters",
         logWeightBase=6.2,
         samplingFraction=ce_ecal_sf)
 
@@ -187,7 +187,7 @@ podout.outputCommands = ["drop *", "keep mcparticles", "keep *Clusters", "keep *
 ApplicationMgr(
     TopAlg = [podin, # copier,
               ce_ecal_cl, ce_ecal_clreco,
-              ce_ecal2_merger, ce_ecal2_cl, ce_ecal2_clreco,
+              ci_ecal_merger, ci_ecal_cl, ci_ecal_clreco,
               cb_ecal_cl, cb_ecal_clreco,
               cb_hcal_merger, cb_hcal_cl, cb_hcal_clreco,
               ce_hcal_merger, ce_hcal_cl, ce_hcal_clreco,

benchmarks/clustering/options/full_cal_digi.py

@@ -40,8 +40,8 @@ from Configurables import Jug__Reco__ImagingClusterReco as ImagingClusterReco
 # branches needed from simulation root file
 sim_coll = [
     "mcparticles",
-    "CrystalEcalHits",
-    "EcalEndcapHits",
+    "EcalEndcapNHits",
+    "EcalEndcapPHits",
     "EcalBarrelHits",
     "HcalBarrelHits",
     "HcalHadronEndcapHits",
@@ -57,8 +57,8 @@ copier = MCCopier("MCCopier", inputCollection="mcparticles", outputCollection="m
 
 # Crystal Endcap Ecal
 ce_ecal_digi = CalHitDigi("ce_ecal_digi",
-        inputHitCollection="CrystalEcalHits",
-        outputHitCollection="CrystalEcalHitsDigi",
+        inputHitCollection="EcalEndcapNHits",
+        outputHitCollection="EcalEndcapNHitsDigi",
         energyResolutions=[0., 0.02, 0.],
         dynamicRangeADC=5.*GeV,     # digi settings must match with reco
         capacityADC=32768,
@@ -66,27 +66,29 @@ ce_ecal_digi = CalHitDigi("ce_ecal_digi",
         pedestalSigma=3)
 
 ce_ecal_reco = CalHitReco("ce_ecal_reco",
-        inputHitCollection="CrystalEcalHitsDigi",
-        outputHitCollection="CrystalEcalHitsReco",
+        inputHitCollection="EcalEndcapNHitsDigi",
+        outputHitCollection="EcalEndcapNHitsReco",
         dynamicRangeADC=5.*GeV,
         capacityADC=32768,
         pedestalMean=400,
         pedestalSigma=3,
         thresholdFactor=4,          # 4 sigma
-        minimumHitEdep=1.0*MeV)     # discard low energy hits
+        minimumHitEdep=1.0*MeV,     # discard low energy hits
+        readoutClass="EcalEndcapNHits",
+        sectorField="sector")
 
 # Endcap Sampling Ecal
-ce_ecal2_digi = CalHitDigi("ce_ecal2_digi",
-         inputHitCollection="EcalEndcapHits",
-         outputHitCollection="EcalEndcapHitsDigi",
+ci_ecal_digi = CalHitDigi("ci_ecal_digi",
+         inputHitCollection="EcalEndcapPHits",
+         outputHitCollection="EcalEndcapPHitsDigi",
          dynamicRangeADC=50.*MeV,
          capacityADC=32768,
          pedestalMean=400,
          pedestalSigma=10)
 
-ce_ecal2_reco = CalHitReco("ce_ecal2_reco",
-        inputHitCollection="EcalEndcapHitsDigi",
-        outputHitCollection="EcalEndcapHitsReco",
+ci_ecal_reco = CalHitReco("ci_ecal_reco",
+        inputHitCollection="EcalEndcapPHitsDigi",
+        outputHitCollection="EcalEndcapPHitsReco",
         dynamicRangeADC=50.*MeV,
         capacityADC=32768,
         pedestalMean=400,
@@ -130,7 +132,10 @@ cb_hcal_reco = CalHitReco("cb_hcal_reco",
         capacityADC=32768,
         pedestalMean=400,
         pedestalSigma=10,
-        thresholdFactor=5.0)
+        thresholdFactor=5.0,
+        readoutClass="HcalBarrelHits",
+        layerField="layer",
+        sectorField="module")
 
 # Hcal Hadron Endcap
 ci_hcal_digi = CalHitDigi("ci_hcal_digi",
@@ -174,7 +179,7 @@ podout.outputCommands = ['drop *', 'keep mcparticles2', 'keep *Reco', 'keep *Dig
 ApplicationMgr(
     TopAlg = [podin, copier,
               ce_ecal_digi, ce_ecal_reco,
-              ce_ecal2_digi, ce_ecal2_reco,
+              ci_ecal_digi, ci_ecal_reco,
               cb_ecal_digi, cb_ecal_reco,
               cb_hcal_digi, cb_hcal_reco,
               ce_hcal_digi, ce_hcal_reco,

benchmarks/clustering/scripts/cluster_plots.py

@@ -91,8 +91,8 @@ def general_clusters_figure(df, collection, save, min_nhits=3):
             data=np.vstack(list(data.values())).T)
     dfp.loc[:, 'evn'] = evns
     # select the max. energy cluster for each event
-    dfp = dfp.loc[dfp.groupby('evn')['edep'].idxmax()]
-    dfp = dfp.loc[dfp['nhits'] >= min_nhits]
+    dfc = dfp.loc[dfp.groupby('evn')['edep'].idxmax()]
+    dfc = dfc.loc[dfc['nhits'] >= min_nhits]
     # figure
     fig, axs = plt.subplots(2, 2, figsize=(16, 12), dpi=120)
     labels = [
@@ -101,7 +101,7 @@ def general_clusters_figure(df, collection, save, min_nhits=3):
         (r'$\theta$ (rad)', 'Counts'),
         (r'$\phi$ (rad)', 'Counts'),
     ]
-    for ax, label, vals in zip(axs.flat, labels, dfp[['nhits', 'edep', 'theta', 'phi']].values.T):
+    for ax, label, vals in zip(axs.flat, labels, dfc[['nhits', 'edep', 'theta', 'phi']].values.T):
         ax.hist(vals, bins=50, ec='k')
         ax.tick_params(labelsize=22, direction='in', which='both')
         ax.grid(linestyle=':', which='both')
@@ -112,6 +112,7 @@ def general_clusters_figure(df, collection, save, min_nhits=3):
     fig.text(0.5, 0.95, collection, ha='center', fontsize=24)
     fig.savefig(save)
     plt.close(fig)
+    return dfp
 
 
 if __name__ == '__main__':
@@ -143,12 +144,18 @@ if __name__ == '__main__':
     rdf_rec = ROOT.RDataFrame('events', args.rec_file)
 
     thrown_particles_figure(rdf_sim, save=os.path.join(args.outdir, 'thrown_particles.png'), mcbranch=args.mc)
-    general_clusters_figure(rdf_rec, collection='CrystalEcalClusters', min_nhits=5,
-            save=os.path.join(args.outdir, 'crystal_ecal_clusters.png'))
-    general_clusters_figure(rdf_rec, collection='EcalEndcapClusters', min_nhits=10,
-            save=os.path.join(args.outdir, 'ecal_endcap_clusters.png'))
+
+    general_clusters_figure(rdf_rec, collection='EcalEndcapNClusters', min_nhits=10,
+            save=os.path.join(args.outdir, 'ecal_electron_endcap_clusters.png'))
+    general_clusters_figure(rdf_rec, collection='EcalEndcapPClusters', min_nhits=5,
+            save=os.path.join(args.outdir, 'ecal_hadron_endcap_clusters.png'))
     general_clusters_figure(rdf_rec, collection='EcalBarrelClusters',
             save=os.path.join(args.outdir, 'ecal_barrel_clusters.png'))
+
+    general_clusters_figure(rdf_rec, collection='HcalElectronEndcapClusters', min_nhits=5,
+            save=os.path.join(args.outdir, 'hcal_electron_endcap_clusters.png'))
+    general_clusters_figure(rdf_rec, collection='HcalHadronEndcapClusters', min_nhits=10,
+            save=os.path.join(args.outdir, 'hcal_hadron_endcap_clusters.png'))
     general_clusters_figure(rdf_rec, collection='HcalBarrelClusters',
             save=os.path.join(args.outdir, 'hcal_barrel_clusters.png'))
 

benchmarks/ecal/options/crystal_calorimeter_reco.py

@@ -45,14 +45,14 @@ from Configurables import Jug__Reco__CrystalEndcapsReco as CrystalEndcapsReco
 from Configurables import Jug__Reco__CalorimeterIslandCluster as IslandCluster
 from Configurables import Jug__Reco__ClusterRecoCoG as RecoCoG
 
-podioinput = PodioInput("PodioReader", collections=["mcparticles","CrystalEcalHits"], OutputLevel=DEBUG)
+podioinput = PodioInput("PodioReader", collections=["mcparticles","EcalEndcapNHits"], OutputLevel=DEBUG)
 
 ## copiers to get around input --> output copy bug. Note the "2" appended to the output collection.
 copier    = MCCopier("MCCopier", inputCollection="mcparticles", outputCollection="mcparticles2",OutputLevel=INFO) 
-calcopier = CalCopier("CalCopier", inputCollection="CrystalEcalHits", outputCollection="CrystalEcalHits2",OutputLevel=INFO) 
+calcopier = CalCopier("CalCopier", inputCollection="EcalEndcapNHits", outputCollection="EcalEndcapNHits2",OutputLevel=INFO) 
 
 emcaldigi = CrystalEndcapsDigi("ecal_digi", 
-        inputHitCollection="CrystalEcalHits", 
+        inputHitCollection="EcalEndcapNHits", 
         outputHitCollection="RawDigiEcalHits")
 
 emcalreco = CrystalEndcapsReco("ecal_reco", 
@@ -69,7 +69,7 @@ emcalcluster = IslandCluster("emcal_cluster",
 clusterreco = RecoCoG("cluster_reco", 
         clusterCollection="EcalClusters", 
         logWeightBase=4.2, 
-        moduleDimZName="CrystalBox_z_length")
+        moduleDimZName="CrystalModule_sz")
 
 out = PodioOutput("out", filename=output_rec_file)
 

benchmarks/ecal/options/example_crystal.py

@@ -29,7 +29,7 @@ emcalreco = CrystalEndcapsReco("ecal_reco", inputHitCollection="RawDigiEcalHits"
                                minModuleEdep=1.0*units.MeV)
 emcalcluster = IslandCluster("emcal_cluster", inputHitCollection="RecoEcalHits", outputClusterCollection="EcalClusters",
                              minClusterCenterEdep=30*units.MeV, groupRange=2.0)
-clusterreco = RecoCoG("cluster_reco", clusterCollection="EcalClusters", logWeightBase=4.2, moduleDimZName="CrystalBox_z_length")
+clusterreco = RecoCoG("cluster_reco", clusterCollection="EcalClusters", logWeightBase=4.2, moduleDimZName="CrystalModule_sz")
 
 out = PodioOutput("out", filename=output_rec_file)
 

benchmarks/ecal/options/full_em_calorimeter_reco.py

@@ -42,7 +42,7 @@ from Configurables import Jug__Reco__CrystalEndcapsReco as CrystalEndcapsReco
 from Configurables import Jug__Reco__CalorimeterIslandCluster as IslandCluster
 from Configurables import Jug__Reco__ClusterRecoCoG as RecoCoG
 
-podioinput = PodioInput("PodioReader", collections=["mcparticles","CrystalEcalHits","EcalBarrelHits","EcalEndcapHits"], OutputLevel=DEBUG)
+podioinput = PodioInput("PodioReader", collections=["mcparticles","EcalEndcapNHits","EcalBarrelHits","EcalEndcapPHits"], OutputLevel=DEBUG)
 
 ## copiers to get around input --> output copy bug. Note the "2" appended to the output collection.
 copier = MCCopier("MCCopier", 
@@ -51,8 +51,8 @@ copier = MCCopier("MCCopier",
         OutputLevel=DEBUG) 
 
 calcopier = CalCopier("CalCopier", 
-        inputCollection="CrystalEcalHits", 
-        outputCollection="CrystalEcalHits2",
+        inputCollection="EcalEndcapNHits", 
+        outputCollection="EcalEndcapNHits2",
         OutputLevel=DEBUG)
 
 embarrelcopier = CalCopier("CalBarrelCopier", 
@@ -71,13 +71,13 @@ embarreldigi = EcalTungstenSamplingDigi("ecal_barrel_digi",
         OutputLevel=DEBUG)
 
 emendcapdigi = EcalTungstenSamplingDigi("ec_endcap_digi", 
-        inputHitCollection="EcalEndcapHits", 
-        outputHitCollection="RawEcalEndcapHits",
+        inputHitCollection="EcalEndcapPHits", 
+        outputHitCollection="RawEcalEndcapPHits",
         energyResolution=0.07,
         OutputLevel=DEBUG)
 
 emcaldigi = CrystalEndcapsDigi("ecal_digi", 
-        inputHitCollection="CrystalEcalHits", 
+        inputHitCollection="EcalEndcapNHits", 
         outputHitCollection="RawDigiEcalHits")
 
 emcalreco = CrystalEndcapsReco("ecal_reco", 
@@ -94,7 +94,7 @@ emcalcluster = IslandCluster("emcal_cluster",
 clusterreco = RecoCoG("cluster_reco", 
         clusterCollection="EcalClusters", 
         logWeightBase=4.2, 
-        moduleDimZName="CrystalBox_z_length")
+        moduleDimZName="CrystalModule_sz")
 
 out = PodioOutput("out", filename=output_rec_file)