Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • EIC/benchmarks/reconstruction_benchmarks
1 result
Show changes
Commits on Source (10)
Hide whitespace changes
Inline Side-by-side
Showing
with 304 additions and 699 deletions
.gitignore
View file @ 22ed77c9
...@@ -41,3 +41,13 @@ __pycache__/ ...@@ -41,3 +41,13 @@ __pycache__/
calorimeters/test/ calorimeters/test/
*.d *.d
*.pcm *.pcm
# transient directories and files
setup
results
*.root
fieldmaps
eic-env.sh
sim_output
*.obj
*.hepmc
.gitlab-ci.yml
View file @ 22ed77c9
...@@ -3,6 +3,7 @@ image: eicweb.phy.anl.gov:4567/eic/juggler/juggler:$JUGGLER_TAG ...@@ -3,6 +3,7 @@ image: eicweb.phy.anl.gov:4567/eic/juggler/juggler:$JUGGLER_TAG
default: default:
before_script: before_script:
- source .local/bin/env.sh - source .local/bin/env.sh
interruptible: true
artifacts: artifacts:
expire_in: 3 days expire_in: 3 days
paths: paths:
...@@ -58,13 +59,14 @@ common:detector: ...@@ -58,13 +59,14 @@ common:detector:
- ln -s "${LOCAL_DATA_PATH}/sim_output" sim_output - ln -s "${LOCAL_DATA_PATH}/sim_output" sim_output
- ln -s "${LOCAL_DATA_PATH}/datasets/data" data - ln -s "${LOCAL_DATA_PATH}/datasets/data" data
- ls -lrtha - ls -lrtha
- bash bin/get_calibrations
interruptible: true interruptible: true
include: include:
- local: 'benchmarks/ecal/config.yml' - local: 'benchmarks/ecal/config.yml'
- local: 'benchmarks/tracking/config.yml' - local: 'benchmarks/tracking/config.yml'
- local: 'benchmarks/clustering/config.yml' - local: 'benchmarks/clustering/config.yml'
- local: 'benchmarks/rich/config.yml' # - local: 'benchmarks/rich/config.yml'
- local: 'benchmarks/imaging_ecal/config.yml' - local: 'benchmarks/imaging_ecal/config.yml'
- local: 'benchmarks/imaging_shower_ML/config.yml' - local: 'benchmarks/imaging_shower_ML/config.yml'
- local: 'benchmarks/full/config.yml' - local: 'benchmarks/full/config.yml'
......
benchmarks/B0/scripts/gen_far_forward_protons.cxx
View file @ 22ed77c9
...@@ -53,12 +53,12 @@ void gen_far_forward_protons(int n_events = 100, ...@@ -53,12 +53,12 @@ void gen_far_forward_protons(int n_events = 100,
//std::cout << std::sqrt(px*px + py*py + pz*pz) - p << " is zero? \n"; //std::cout << std::sqrt(px*px + py*py + pz*pz) - p << " is zero? \n";
// type 1 is final state // type 1 is final state
// pdgid 11 - proton 0.510 MeV/c^2 // pdgid 11 - proton 938.272 MeV/c^2
GenParticlePtr p3 = std::make_shared<GenParticle>( GenParticlePtr p3 = std::make_shared<GenParticle>(
FourVector( FourVector(
px, py, pz, px, py, pz,
sqrt(p*p + (0.000511 * 0.000511))), sqrt(p*p + (0.938272 * 0.938272))),
11, 1); 2212, 1);
GenVertexPtr v1 = std::make_shared<GenVertex>(); GenVertexPtr v1 = std::make_shared<GenVertex>();
v1->add_particle_in(p1); v1->add_particle_in(p1);
......
benchmarks/B0/scripts/hits_far_forward_protons.cxx
View file @ 22ed77c9
...@@ -31,7 +31,7 @@ auto p_track = [](std::vector<eic::TrackParametersData> const& in) { ...@@ -31,7 +31,7 @@ auto p_track = [](std::vector<eic::TrackParametersData> const& in) {
std::vector<float> pt (std::vector<dd4pod::Geant4ParticleData> const& in){ std::vector<float> pt (std::vector<dd4pod::Geant4ParticleData> const& in){
std::vector<float> result; std::vector<float> result;
for (size_t i = 0; i < in.size(); ++i) { for (size_t i = 0; i < in.size(); ++i) {
result.push_back(std::sqrt(in[i].psx * in[i].psx + in[i].psy * in[i].psy)); result.push_back(std::sqrt(in[i].ps.x * in[i].ps.x + in[i].ps.y * in[i].ps.y));
} }
return result; return result;
} }
...@@ -54,7 +54,7 @@ auto fourvec = [](ROOT::VecOps::RVec<dd4pod::Geant4ParticleData> const& in) { ...@@ -54,7 +54,7 @@ auto fourvec = [](ROOT::VecOps::RVec<dd4pod::Geant4ParticleData> const& in) {
std::vector<ROOT::Math::PxPyPzMVector> result; std::vector<ROOT::Math::PxPyPzMVector> result;
ROOT::Math::PxPyPzMVector lv; ROOT::Math::PxPyPzMVector lv;
for (size_t i = 0; i < in.size(); ++i) { for (size_t i = 0; i < in.size(); ++i) {
lv.SetCoordinates(in[i].psx, in[i].psy, in[i].psz, in[i].mass); lv.SetCoordinates(in[i].ps.x, in[i].ps.y, in[i].ps.z, in[i].mass);
result.push_back(lv); result.push_back(lv);
} }
return result; return result;
......
benchmarks/B0/scripts/rec_far_forward_protons.cxx
View file @ 22ed77c9
...@@ -29,7 +29,7 @@ auto p_track = [](std::vector<eic::TrackParametersData> const& in) { ...@@ -29,7 +29,7 @@ auto p_track = [](std::vector<eic::TrackParametersData> const& in) {
std::vector<float> pt (std::vector<dd4pod::Geant4ParticleData> const& in){ std::vector<float> pt (std::vector<dd4pod::Geant4ParticleData> const& in){
std::vector<float> result; std::vector<float> result;
for (size_t i = 0; i < in.size(); ++i) { for (size_t i = 0; i < in.size(); ++i) {
result.push_back(std::sqrt(in[i].psx * in[i].psx + in[i].psy * in[i].psy)); result.push_back(std::sqrt(in[i].ps.x * in[i].ps.x + in[i].ps.y * in[i].ps.y));
} }
return result; return result;
} }
...@@ -52,7 +52,7 @@ auto fourvec = [](ROOT::VecOps::RVec<dd4pod::Geant4ParticleData> const& in) { ...@@ -52,7 +52,7 @@ auto fourvec = [](ROOT::VecOps::RVec<dd4pod::Geant4ParticleData> const& in) {
std::vector<ROOT::Math::PxPyPzMVector> result; std::vector<ROOT::Math::PxPyPzMVector> result;
ROOT::Math::PxPyPzMVector lv; ROOT::Math::PxPyPzMVector lv;
for (size_t i = 0; i < in.size(); ++i) { for (size_t i = 0; i < in.size(); ++i) {
lv.SetCoordinates(in[i].psx, in[i].psy, in[i].psz, in[i].mass); lv.SetCoordinates(in[i].ps.x, in[i].ps.y, in[i].ps.z, in[i].mass);
result.push_back(lv); result.push_back(lv);
} }
return result; return result;
......
benchmarks/clustering/full_cal_clusters.sh
View file @ 22ed77c9
...@@ -131,7 +131,7 @@ fi ...@@ -131,7 +131,7 @@ fi
FULL_CAL_SCRIPT_DIR=benchmarks/clustering/scripts FULL_CAL_SCRIPT_DIR=benchmarks/clustering/scripts
python ${FULL_CAL_SCRIPT_DIR}/cluster_plots.py ${JUGGLER_SIM_FILE} ${JUGGLER_REC_FILE} -o results \ python ${FULL_CAL_SCRIPT_DIR}/cluster_plots.py ${JUGGLER_SIM_FILE} ${JUGGLER_REC_FILE} -o results \
--collections "EcalEndcapNClusters, EcalEndcapPClusters, EcalBarrelClusters, --collections "EcalEndcapNClusters, EcalEndcapPClusters, EcalBarrelClusters,
HcalElectronEndcapClusters, HcalHadronEndcapClusters, HcalBarrelClusters" HcalEndcapNClusters, HcalEndcapPClusters, HcalBarrelClusters"
root_filesize=$(stat --format=%s "${JUGGLER_REC_FILE}") root_filesize=$(stat --format=%s "${JUGGLER_REC_FILE}")
if [[ "${JUGGLER_N_EVENTS}" -lt "500" ]] ; then if [[ "${JUGGLER_N_EVENTS}" -lt "500" ]] ; then
......
benchmarks/clustering/options/full_cal_reco.py
View file @ 22ed77c9
...@@ -2,6 +2,7 @@ ...@@ -2,6 +2,7 @@
An example option file to digitize/reconstruct/clustering calorimeter hits An example option file to digitize/reconstruct/clustering calorimeter hits
''' '''
from Gaudi.Configuration import * from Gaudi.Configuration import *
import json
import os import os
import ROOT import ROOT
...@@ -12,13 +13,21 @@ detector_name = str(os.environ.get("JUGGLER_DETECTOR", "athena")) ...@@ -12,13 +13,21 @@ detector_name = str(os.environ.get("JUGGLER_DETECTOR", "athena"))
detector_path = str(os.environ.get("DETECTOR_PATH", ".")) detector_path = str(os.environ.get("DETECTOR_PATH", "."))
compact_path = os.path.join(detector_path, detector_name) compact_path = os.path.join(detector_path, detector_name)
# input arguments from calibration file
with open('config/emcal_barrel_calibration.json') as f:
calib_data = json.load(f)['electron']
print(calib_data)
cb_ecal_sf = float(calib_data['sampling_fraction_img'])
scifi_barrel_sf = float(calib_data['sampling_fraction_scfi'])
# get sampling fractions from system environment variable, 1.0 by default # get sampling fractions from system environment variable, 1.0 by default
ci_ecal_sf = float(os.environ.get("CI_ECAL_SAMP_FRAC", 0.253)) ci_ecal_sf = float(os.environ.get("CI_ECAL_SAMP_FRAC", 0.253))
cb_ecal_sf = float(os.environ.get("CB_ECAL_SAMP_FRAC", 0.01324))
cb_hcal_sf = float(os.environ.get("CB_HCAL_SAMP_FRAC", 0.038)) cb_hcal_sf = float(os.environ.get("CB_HCAL_SAMP_FRAC", 0.038))
ci_hcal_sf = float(os.environ.get("CI_HCAL_SAMP_FRAC", 0.025)) ci_hcal_sf = float(os.environ.get("CI_HCAL_SAMP_FRAC", 0.025))
ce_hcal_sf = float(os.environ.get("CE_HCAL_SAMP_FRAC", 0.025)) ce_hcal_sf = float(os.environ.get("CE_HCAL_SAMP_FRAC", 0.025))
scifi_barrel_sf = float(os.environ.get("CB_EMCAL_SCFI_SAMP_FRAC",0.0938))
# input and output # input and output
input_sims = [f.strip() for f in str.split(os.environ["JUGGLER_SIM_FILE"], ",") if f.strip()] input_sims = [f.strip() for f in str.split(os.environ["JUGGLER_SIM_FILE"], ",") if f.strip()]
output_rec = str(os.environ["JUGGLER_REC_FILE"]) output_rec = str(os.environ["JUGGLER_REC_FILE"])
...@@ -51,8 +60,8 @@ sim_coll = [ ...@@ -51,8 +60,8 @@ sim_coll = [
"EcalBarrelHits", "EcalBarrelHits",
"HcalBarrelHits", "HcalBarrelHits",
"EcalBarrelScFiHits", "EcalBarrelScFiHits",
"HcalHadronEndcapHits", "HcalEndcapPHits",
"HcalElectronEndcapHits", "HcalEndcapNHits",
] ]
# input and output # input and output
...@@ -76,7 +85,7 @@ ce_ecal_digi = CalHitDigi("ce_ecal_digi", ...@@ -76,7 +85,7 @@ ce_ecal_digi = CalHitDigi("ce_ecal_digi",
**ce_ecal_daq) **ce_ecal_daq)
ce_ecal_reco = CalHitReco("ce_ecal_reco", ce_ecal_reco = CalHitReco("ce_ecal_reco",
inputHitCollection="EcalEndcapNHitsDigi", inputHitCollection=ce_ecal_digi.outputHitCollection,
outputHitCollection="EcalEndcapNHitsReco", outputHitCollection="EcalEndcapNHitsReco",
thresholdFactor=4, # 4 sigma cut on pedestal sigma thresholdFactor=4, # 4 sigma cut on pedestal sigma
readoutClass="EcalEndcapNHits", readoutClass="EcalEndcapNHits",
...@@ -85,8 +94,8 @@ ce_ecal_reco = CalHitReco("ce_ecal_reco", ...@@ -85,8 +94,8 @@ ce_ecal_reco = CalHitReco("ce_ecal_reco",
ce_ecal_cl = IslandCluster("ce_ecal_cl", ce_ecal_cl = IslandCluster("ce_ecal_cl",
# OutputLevel=DEBUG, # OutputLevel=DEBUG,
inputHitCollection="EcalEndcapNHitsReco", inputHitCollection=ce_ecal_reco.outputHitCollection,
outputHitCollection="EcalEndcapNClusterHits", outputProtoClusterCollection="EcalEndcapNProtoClusters",
splitCluster=False, splitCluster=False,
minClusterHitEdep=1.0*MeV, # discard low energy hits minClusterHitEdep=1.0*MeV, # discard low energy hits
minClusterCenterEdep=30*MeV, minClusterCenterEdep=30*MeV,
...@@ -94,8 +103,10 @@ ce_ecal_cl = IslandCluster("ce_ecal_cl", ...@@ -94,8 +103,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 dimScaledLocalDistXY=[1.8, 1.8]) # dimension scaled dist is good for hybrid sectors with different module size
ce_ecal_clreco = RecoCoG("ce_ecal_clreco", ce_ecal_clreco = RecoCoG("ce_ecal_clreco",
inputHitCollection="EcalEndcapNClusterHits", inputHitCollection=ce_ecal_cl.inputHitCollection,
inputProtoClusterCollection=ce_ecal_cl.outputProtoClusterCollection,
outputClusterCollection="EcalEndcapNClusters", outputClusterCollection="EcalEndcapNClusters",
outputInfoCollection="EcalEndcapNClustersInfo",
samplingFraction=0.998, # this accounts for a small fraction of leakage samplingFraction=0.998, # this accounts for a small fraction of leakage
logWeightBase=4.6) logWeightBase=4.6)
...@@ -112,7 +123,7 @@ ci_ecal_digi = CalHitDigi("ci_ecal_digi", ...@@ -112,7 +123,7 @@ ci_ecal_digi = CalHitDigi("ci_ecal_digi",
**ci_ecal_daq) **ci_ecal_daq)
ci_ecal_reco = CalHitReco("ci_ecal_reco", ci_ecal_reco = CalHitReco("ci_ecal_reco",
inputHitCollection="EcalEndcapPHitsDigi", inputHitCollection=ci_ecal_digi.outputHitCollection,
outputHitCollection="EcalEndcapPHitsReco", outputHitCollection="EcalEndcapPHitsReco",
thresholdFactor=5.0, thresholdFactor=5.0,
**ci_ecal_daq) **ci_ecal_daq)
...@@ -120,7 +131,7 @@ ci_ecal_reco = CalHitReco("ci_ecal_reco", ...@@ -120,7 +131,7 @@ ci_ecal_reco = CalHitReco("ci_ecal_reco",
# merge hits in different layer (projection to local x-y plane) # merge hits in different layer (projection to local x-y plane)
ci_ecal_merger = CalHitsMerger("ci_ecal_merger", ci_ecal_merger = CalHitsMerger("ci_ecal_merger",
# OutputLevel=DEBUG, # OutputLevel=DEBUG,
inputHitCollection="EcalEndcapPHitsReco", inputHitCollection=ci_ecal_reco.outputHitCollection,
outputHitCollection="EcalEndcapPHitsRecoXY", outputHitCollection="EcalEndcapPHitsRecoXY",
fields=["layer", "slice"], fields=["layer", "slice"],
fieldRefNumbers=[1, 0], fieldRefNumbers=[1, 0],
...@@ -128,15 +139,17 @@ ci_ecal_merger = CalHitsMerger("ci_ecal_merger", ...@@ -128,15 +139,17 @@ ci_ecal_merger = CalHitsMerger("ci_ecal_merger",
ci_ecal_cl = IslandCluster("ci_ecal_cl", ci_ecal_cl = IslandCluster("ci_ecal_cl",
# OutputLevel=DEBUG, # OutputLevel=DEBUG,
inputHitCollection="EcalEndcapPHitsRecoXY", inputHitCollection=ci_ecal_merger.inputHitCollection,
outputHitCollection="EcalEndcapPClusterHits", outputProtoClusterCollection="EcalEndcapPProtoClusters",
splitCluster=False, splitCluster=False,
minClusterCenterEdep=10.*MeV, minClusterCenterEdep=10.*MeV,
localDistXY=[10*mm, 10*mm]) localDistXY=[10*mm, 10*mm])
ci_ecal_clreco = RecoCoG("ci_ecal_clreco", ci_ecal_clreco = RecoCoG("ci_ecal_clreco",
inputHitCollection="EcalEndcapPClusterHits", inputHitCollection=ci_ecal_cl.inputHitCollection,
inputProtoClusterCollection=ci_ecal_cl.outputProtoClusterCollection,
outputClusterCollection="EcalEndcapPClusters", outputClusterCollection="EcalEndcapPClusters",
outputInfoCollection="EcalEndcapPClustersInfo",
logWeightBase=6.2, logWeightBase=6.2,
samplingFraction=ci_ecal_sf) samplingFraction=ci_ecal_sf)
...@@ -154,7 +167,7 @@ cb_ecal_digi = CalHitDigi("cb_ecal_digi", ...@@ -154,7 +167,7 @@ cb_ecal_digi = CalHitDigi("cb_ecal_digi",
**cb_ecal_daq) **cb_ecal_daq)
cb_ecal_reco = ImCalPixelReco("cb_ecal_reco", cb_ecal_reco = ImCalPixelReco("cb_ecal_reco",
inputHitCollection="EcalBarrelHitsDigi", inputHitCollection=cb_ecal_digi.outputHitCollection,
outputHitCollection="EcalBarrelHitsReco", outputHitCollection="EcalBarrelHitsReco",
thresholdFactor=3, # about 20 keV thresholdFactor=3, # about 20 keV
readoutClass="EcalBarrelHits", # readout class readoutClass="EcalBarrelHits", # readout class
...@@ -163,8 +176,8 @@ cb_ecal_reco = ImCalPixelReco("cb_ecal_reco", ...@@ -163,8 +176,8 @@ cb_ecal_reco = ImCalPixelReco("cb_ecal_reco",
**cb_ecal_daq) **cb_ecal_daq)
cb_ecal_cl = ImagingCluster("cb_ecal_cl", cb_ecal_cl = ImagingCluster("cb_ecal_cl",
inputHitCollection="EcalBarrelHitsReco", inputHitCollection=cb_ecal_reco.outputHitCollection,
outputHitCollection="EcalBarrelClusterHits", outputProtoClusterCollection="EcalBarrelProtoClusters",
localDistXY=[2.*mm, 2*mm], # same layer localDistXY=[2.*mm, 2*mm], # same layer
layerDistEtaPhi=[10*mrad, 10*mrad], # adjacent layer layerDistEtaPhi=[10*mrad, 10*mrad], # adjacent layer
neighbourLayersRange=2, # id diff for adjacent layer neighbourLayersRange=2, # id diff for adjacent layer
...@@ -172,8 +185,10 @@ cb_ecal_cl = ImagingCluster("cb_ecal_cl", ...@@ -172,8 +185,10 @@ cb_ecal_cl = ImagingCluster("cb_ecal_cl",
cb_ecal_clreco = ImagingClusterReco("cb_ecal_clreco", cb_ecal_clreco = ImagingClusterReco("cb_ecal_clreco",
samplingFraction=cb_ecal_sf, samplingFraction=cb_ecal_sf,
inputHitCollection="EcalBarrelClusterHits", inputHitCollection=cb_ecal_cl.inputHitCollection,
inputProtoClusterCollection=cb_ecal_cl.outputProtoClusterCollection,
outputClusterCollection="EcalBarrelClusters", outputClusterCollection="EcalBarrelClusters",
outputInfoCollection="EcalBarrelClustersInfo",
outputLayerCollection="EcalBarrelLayers") outputLayerCollection="EcalBarrelLayers")
#Central ECAL SciFi #Central ECAL SciFi
...@@ -190,7 +205,7 @@ scfi_barrel_digi = CalHitDigi("scfi_barrel_digi", ...@@ -190,7 +205,7 @@ scfi_barrel_digi = CalHitDigi("scfi_barrel_digi",
**scfi_barrel_daq) **scfi_barrel_daq)
scfi_barrel_reco = CalHitReco("scfi_barrel_reco", scfi_barrel_reco = CalHitReco("scfi_barrel_reco",
inputHitCollection="EcalBarrelScFiHitsDigi", inputHitCollection=scfi_barrel_digi.outputHitCollection,
outputHitCollection="EcalBarrelScFiHitsReco", outputHitCollection="EcalBarrelScFiHitsReco",
thresholdFactor=5.0, thresholdFactor=5.0,
readoutClass="EcalBarrelScFiHits", readoutClass="EcalBarrelScFiHits",
...@@ -202,7 +217,7 @@ scfi_barrel_reco = CalHitReco("scfi_barrel_reco", ...@@ -202,7 +217,7 @@ scfi_barrel_reco = CalHitReco("scfi_barrel_reco",
# merge hits in different layer (projection to local x-y plane) # merge hits in different layer (projection to local x-y plane)
scfi_barrel_merger = CalHitsMerger("scfi_barrel_merger", scfi_barrel_merger = CalHitsMerger("scfi_barrel_merger",
# OutputLevel=DEBUG, # OutputLevel=DEBUG,
inputHitCollection="EcalBarrelScFiHitsReco", inputHitCollection=scfi_barrel_reco.outputHitCollection,
outputHitCollection="EcalBarrelScFiGridReco", outputHitCollection="EcalBarrelScFiGridReco",
fields=["fiber"], fields=["fiber"],
fieldRefNumbers=[1], fieldRefNumbers=[1],
...@@ -210,15 +225,17 @@ scfi_barrel_merger = CalHitsMerger("scfi_barrel_merger", ...@@ -210,15 +225,17 @@ scfi_barrel_merger = CalHitsMerger("scfi_barrel_merger",
scfi_barrel_cl = IslandCluster("scfi_barrel_cl", scfi_barrel_cl = IslandCluster("scfi_barrel_cl",
# OutputLevel=DEBUG, # OutputLevel=DEBUG,
inputHitCollection="EcalBarrelScFiGridReco", inputHitCollection=scfi_barrel_merger.outputHitCollection,
outputHitCollection="EcalBarrelScFiClusterHits", outputProtoClusterCollection="EcalBarrelScFiProtoClusters",
splitCluster=False, splitCluster=False,
minClusterCenterEdep=10.*MeV, minClusterCenterEdep=10.*MeV,
localDistXZ=[30*mm, 30*mm]) localDistXZ=[30*mm, 30*mm])
scfi_barrel_clreco = RecoCoG("scfi_barrel_clreco", scfi_barrel_clreco = RecoCoG("scfi_barrel_clreco",
inputHitCollection="EcalBarrelScFiClusterHits", inputHitCollection=scfi_barrel_cl.inputHitCollection,
inputProtoClusterCollection=scfi_barrel_cl.outputProtoClusterCollection,
outputClusterCollection="EcalBarrelScFiClusters", outputClusterCollection="EcalBarrelScFiClusters",
outputInfoCollection="EcalBarrelScFiClustersInfo",
logWeightBase=6.2, logWeightBase=6.2,
samplingFraction= scifi_barrel_sf) samplingFraction= scifi_barrel_sf)
...@@ -236,7 +253,7 @@ cb_hcal_digi = CalHitDigi("cb_hcal_digi", ...@@ -236,7 +253,7 @@ cb_hcal_digi = CalHitDigi("cb_hcal_digi",
**cb_hcal_daq) **cb_hcal_daq)
cb_hcal_reco = CalHitReco("cb_hcal_reco", cb_hcal_reco = CalHitReco("cb_hcal_reco",
inputHitCollection="HcalBarrelHitsDigi", inputHitCollection=cb_hcal_digi.outputHitCollection,
outputHitCollection="HcalBarrelHitsReco", outputHitCollection="HcalBarrelHitsReco",
thresholdFactor=5.0, thresholdFactor=5.0,
readoutClass="HcalBarrelHits", readoutClass="HcalBarrelHits",
...@@ -245,22 +262,24 @@ cb_hcal_reco = CalHitReco("cb_hcal_reco", ...@@ -245,22 +262,24 @@ cb_hcal_reco = CalHitReco("cb_hcal_reco",
**cb_hcal_daq) **cb_hcal_daq)
cb_hcal_merger = CalHitsMerger("cb_hcal_merger", cb_hcal_merger = CalHitsMerger("cb_hcal_merger",
inputHitCollection="HcalBarrelHitsReco", inputHitCollection=cb_hcal_reco.outputHitCollection,
outputHitCollection="HcalBarrelHitsRecoXY", outputHitCollection="HcalBarrelHitsRecoXY",
readoutClass="HcalBarrelHits", readoutClass="HcalBarrelHits",
fields=["layer", "slice"], fields=["layer", "slice"],
fieldRefNumbers=[1, 0]) fieldRefNumbers=[1, 0])
cb_hcal_cl = IslandCluster("cb_hcal_cl", cb_hcal_cl = IslandCluster("cb_hcal_cl",
inputHitCollection="HcalBarrelHitsRecoXY", inputHitCollection=cb_hcal_merger.outputHitCollection,
outputHitCollection="HcalBarrelClusterHits", outputProtoClusterCollection="HcalBarrelProtoClusters",
splitCluster=False, splitCluster=False,
minClusterCenterEdep=30.*MeV, minClusterCenterEdep=30.*MeV,
localDistXY=[15.*cm, 15.*cm]) localDistXY=[15.*cm, 15.*cm])
cb_hcal_clreco = RecoCoG("cb_hcal_clreco", cb_hcal_clreco = RecoCoG("cb_hcal_clreco",
inputHitCollection="HcalBarrelClusterHits", inputHitCollection=cb_hcal_cl.inputHitCollection,
inputProtoClusterCollection=cb_hcal_cl.outputProtoClusterCollection,
outputClusterCollection="HcalBarrelClusters", outputClusterCollection="HcalBarrelClusters",
outputInfoCollection="HcalBarrelClustersInfo",
logWeightBase=6.2, logWeightBase=6.2,
samplingFraction=cb_hcal_sf) samplingFraction=cb_hcal_sf)
...@@ -272,33 +291,35 @@ ci_hcal_daq = dict( ...@@ -272,33 +291,35 @@ ci_hcal_daq = dict(
pedestalMean=400, pedestalMean=400,
pedestalSigma=10) pedestalSigma=10)
ci_hcal_digi = CalHitDigi("ci_hcal_digi", ci_hcal_digi = CalHitDigi("ci_hcal_digi",
inputHitCollection="HcalHadronEndcapHits", inputHitCollection="HcalEndcapPHits",
outputHitCollection="HcalHadronEndcapHitsDigi", outputHitCollection="HcalEndcapPHitsDigi",
**ci_hcal_daq) **ci_hcal_daq)
ci_hcal_reco = CalHitReco("ci_hcal_reco", ci_hcal_reco = CalHitReco("ci_hcal_reco",
inputHitCollection="HcalHadronEndcapHitsDigi", inputHitCollection=ci_hcal_digi.outputHitCollection,
outputHitCollection="HcalHadronEndcapHitsReco", outputHitCollection="HcalEndcapPHitsReco",
thresholdFactor=5.0, thresholdFactor=5.0,
**ci_hcal_daq) **ci_hcal_daq)
ci_hcal_merger = CalHitsMerger("ci_hcal_merger", ci_hcal_merger = CalHitsMerger("ci_hcal_merger",
inputHitCollection="HcalHadronEndcapHitsReco", inputHitCollection=ci_hcal_reco.outputHitCollection,
outputHitCollection="HcalHadronEndcapHitsRecoXY", outputHitCollection="HcalEndcapPHitsRecoXY",
readoutClass="HcalHadronEndcapHits", readoutClass="HcalEndcapPHits",
fields=["layer", "slice"], fields=["layer", "slice"],
fieldRefNumbers=[1, 0]) fieldRefNumbers=[1, 0])
ci_hcal_cl = IslandCluster("ci_hcal_cl", ci_hcal_cl = IslandCluster("ci_hcal_cl",
inputHitCollection="HcalHadronEndcapHitsRecoXY", inputHitCollection=ci_hcal_merger.outputHitCollection,
outputHitCollection="HcalHadronEndcapClusterHits", outputProtoClusterCollection="HcalEndcapPProtoClusters",
splitCluster=False, splitCluster=False,
minClusterCenterEdep=30.*MeV, minClusterCenterEdep=30.*MeV,
localDistXY=[15.*cm, 15.*cm]) localDistXY=[15.*cm, 15.*cm])
ci_hcal_clreco = RecoCoG("ci_hcal_clreco", ci_hcal_clreco = RecoCoG("ci_hcal_clreco",
inputHitCollection="HcalHadronEndcapClusterHits", inputHitCollection=ci_hcal_cl.inputHitCollection,
outputClusterCollection="HcalHadronEndcapClusters", inputProtoClusterCollection=ci_hcal_cl.outputProtoClusterCollection,
outputClusterCollection="HcalEndcapPClusters",
outputInfoCollection="HcalEndcapPClustersInfo",
logWeightBase=6.2, logWeightBase=6.2,
samplingFraction=ci_hcal_sf) samplingFraction=ci_hcal_sf)
...@@ -310,33 +331,35 @@ ce_hcal_daq = dict( ...@@ -310,33 +331,35 @@ ce_hcal_daq = dict(
pedestalSigma=10) pedestalSigma=10)
ce_hcal_digi = CalHitDigi("ce_hcal_digi", ce_hcal_digi = CalHitDigi("ce_hcal_digi",
inputHitCollection="HcalElectronEndcapHits", inputHitCollection="HcalEndcapNHits",
outputHitCollection="HcalElectronEndcapHitsDigi", outputHitCollection="HcalEndcapNHitsDigi",
**ce_hcal_daq) **ce_hcal_daq)
ce_hcal_reco = CalHitReco("ce_hcal_reco", ce_hcal_reco = CalHitReco("ce_hcal_reco",
inputHitCollection="HcalElectronEndcapHitsDigi", inputHitCollection=ce_hcal_digi.outputHitCollection,
outputHitCollection="HcalElectronEndcapHitsReco", outputHitCollection="HcalEndcapNHitsReco",
thresholdFactor=5.0, thresholdFactor=5.0,
**ce_hcal_daq) **ce_hcal_daq)
ce_hcal_merger = CalHitsMerger("ce_hcal_merger", ce_hcal_merger = CalHitsMerger("ce_hcal_merger",
inputHitCollection="HcalElectronEndcapHitsReco", inputHitCollection=ce_hcal_reco.outputHitCollection,
outputHitCollection="HcalElectronEndcapHitsRecoXY", outputHitCollection="HcalEndcapNHitsRecoXY",
readoutClass="HcalElectronEndcapHits", readoutClass="HcalEndcapNHits",
fields=["layer", "slice"], fields=["layer", "slice"],
fieldRefNumbers=[1, 0]) fieldRefNumbers=[1, 0])
ce_hcal_cl = IslandCluster("ce_hcal_cl", ce_hcal_cl = IslandCluster("ce_hcal_cl",
inputHitCollection="HcalElectronEndcapHitsRecoXY", inputHitCollection=ce_hcal_merger.outputHitCollection,
outputHitCollection="HcalElectronEndcapClusterHits", outputProtoClusterCollection="HcalEndcapNProtoClusters",
splitCluster=False, splitCluster=False,
minClusterCenterEdep=30.*MeV, minClusterCenterEdep=30.*MeV,
localDistXY=[15.*cm, 15.*cm]) localDistXY=[15.*cm, 15.*cm])
ce_hcal_clreco = RecoCoG("ce_hcal_clreco", ce_hcal_clreco = RecoCoG("ce_hcal_clreco",
inputHitCollection="HcalElectronEndcapClusterHits", inputHitCollection=ce_hcal_cl.inputHitCollection,
outputClusterCollection="HcalElectronEndcapClusters", inputProtoClusterCollection=ce_hcal_cl.outputProtoClusterCollection,
outputClusterCollection="HcalEndcapNClusters",
outputInfoCollection="HcalEndcapNClustersInfo",
logWeightBase=6.2, logWeightBase=6.2,
samplingFraction=ce_hcal_sf) samplingFraction=ce_hcal_sf)
...@@ -345,8 +368,7 @@ podout.outputCommands = ['drop *', ...@@ -345,8 +368,7 @@ podout.outputCommands = ['drop *',
'keep mcparticles2', 'keep mcparticles2',
'keep *Digi', 'keep *Digi',
'keep *Reco*', 'keep *Reco*',
'keep *ClusterHits', 'keep *Cluster*',
'keep *Clusters',
'keep *Layers'] 'keep *Layers']
ApplicationMgr( ApplicationMgr(
......
benchmarks/clustering/options/full_trackpluscalo_reco.py deleted 100644 → 0
View file @ 1e66addd
'''
An example option file to digitize/reconstruct/clustering calorimeter hits
'''
from Gaudi.Configuration import *
import os
import ROOT
from Configurables import ApplicationMgr, EICDataSvc, PodioInput, PodioOutput, GeoSvc
from GaudiKernel.SystemOfUnits import MeV, GeV, mm, cm, mrad
detector_name = str(os.environ.get("JUGGLER_DETECTOR", "athena"))
detector_path = str(os.environ.get("DETECTOR_PATH", "."))
compact_path = os.path.join(detector_path, detector_name)
# get sampling fractions from system environment variable, 1.0 by default
ci_ecal_sf = float(os.environ.get("CI_ECAL_SAMP_FRAC", 0.253))
cb_ecal_sf = float(os.environ.get("CB_ECAL_SAMP_FRAC", 0.01324))
cb_hcal_sf = float(os.environ.get("CB_HCAL_SAMP_FRAC", 0.038))
ci_hcal_sf = float(os.environ.get("CI_HCAL_SAMP_FRAC", 0.025))
ce_hcal_sf = float(os.environ.get("CE_HCAL_SAMP_FRAC", 0.025))
scifi_barrel_sf = float(os.environ.get("CB_EMCAL_SCFI_SAMP_FRAC",0.0938))
# input and output
input_sims = [f.strip() for f in str.split(os.environ["JUGGLER_SIM_FILE"], ",") if f.strip()]
output_rec = str(os.environ["JUGGLER_REC_FILE"])
n_events = int(os.environ["JUGGLER_N_EVENTS"])
# geometry service
geo_service = GeoSvc("GeoSvc", detectors=["{}.xml".format(compact_path)], OutputLevel=INFO)
# data service
podioevent = EICDataSvc("EventDataSvc", inputs=input_sims)
# juggler components
from Configurables import Jug__Base__InputCopier_dd4pod__Geant4ParticleCollection_dd4pod__Geant4ParticleCollection_ as MCCopier
from Configurables import Jug__Base__InputCopier_dd4pod__CalorimeterHitCollection_dd4pod__CalorimeterHitCollection_ as CalCopier
from Configurables import Jug__Base__InputCopier_dd4pod__TrackerHitCollection_dd4pod__TrackerHitCollection_ as TrkCopier
from Configurables import Jug__Digi__UFSDTrackerDigi as TrackerDigi
from Configurables import Jug__Reco__TrackerHitReconstruction as TrackerHitReconstruction
from Configurables import Jug__Reco__TrackerSourceLinker as TrackerSourceLinker
from Configurables import Jug__Reco__TrackerSourcesLinker as TrackerSourcesLinker
#from Configurables import Jug__Reco__TrackingHitsSourceLinker as TrackingHitsSourceLinker
from Configurables import Jug__Reco__TrackParamTruthInit as TrackParamTruthInit
from Configurables import Jug__Reco__TrackParamClusterInit as TrackParamClusterInit
from Configurables import Jug__Reco__TrackParamVertexClusterInit as TrackParamVertexClusterInit
from Configurables import Jug__Reco__TrackFindingAlgorithm as TrackFindingAlgorithm
from Configurables import Jug__Reco__ParticlesFromTrackFit as ParticlesFromTrackFit
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__ImagingPixelReco as ImCalPixelReco
from Configurables import Jug__Reco__ImagingTopoCluster as ImagingCluster
from Configurables import Jug__Reco__ClusterRecoCoG as RecoCoG
from Configurables import Jug__Reco__ImagingClusterReco as ImagingClusterReco
# branches needed from simulation root file
sim_coll = [
"mcparticles",
"EcalEndcapNHits",
"EcalEndcapPHits",
"EcalBarrelHits",
"HcalBarrelHits",
"EcalBarrelScFiHits",
"HcalHadronEndcapHits",
"HcalElectronEndcapHits",
"TrackerEndcapHits",
"TrackerBarrelHits",
"VertexBarrelHits",
"VertexEndcapHits"
]
# input and output
podin = PodioInput("PodioReader", collections=sim_coll)
podout = PodioOutput("out", filename=output_rec)
# copier needed to get around input --> output copy bug. So truth (mcparticles) can be saved in output file
copier = MCCopier("MCCopier", inputCollection="mcparticles", outputCollection="mcparticles2")
trkcopier = TrkCopier("TrkCopier",
inputCollection="TrackerBarrelHits",
outputCollection="TrackerBarrelHits2")
##Tracking digi
trk_b_digi = TrackerDigi("trk_b_digi",
inputHitCollection="TrackerBarrelHits",
outputHitCollection="TrackerBarrelRawHits",
timeResolution=8)
trk_ec_digi = TrackerDigi("trk_ec_digi",
inputHitCollection="TrackerEndcapHits",
outputHitCollection="TrackerEndcapRawHits",
timeResolution=8)
vtx_b_digi = TrackerDigi("vtx_b_digi",
inputHitCollection="VertexBarrelHits",
outputHitCollection="VertexBarrelRawHits",
timeResolution=8)
vtx_ec_digi = TrackerDigi("vtx_ec_digi",
inputHitCollection="VertexEndcapHits",
outputHitCollection="VertexEndcapRawHits",
timeResolution=8)
# Tracker and vertex reconstruction
trk_b_reco = TrackerHitReconstruction("trk_b_reco",
inputHitCollection = trk_b_digi.outputHitCollection,
outputHitCollection="TrackerBarrelRecHits")
trk_ec_reco = TrackerHitReconstruction("trk_ec_reco",
inputHitCollection = trk_ec_digi.outputHitCollection,
outputHitCollection="TrackerEndcapRecHits")
vtx_b_reco = TrackerHitReconstruction("vtx_b_reco",
inputHitCollection = vtx_b_digi.outputHitCollection,
outputHitCollection="VertexBarrelRecHits")
vtx_ec_reco = TrackerHitReconstruction("vtx_ec_reco",
inputHitCollection = vtx_ec_digi.outputHitCollection,
outputHitCollection="VertexEndcapRecHits")
# Hit Source linker
sourcelinker = TrackerSourceLinker("sourcelinker",
inputHitCollection=trk_b_reco.outputHitCollection,
outputSourceLinks="BarrelTrackSourceLinks",
outputMeasurements="BarrelTrackMeasurements",
OutputLevel=DEBUG)
#trk_hits_srclnkr = TrackerSourcesLinker("trk_srcslnkr",
# ITrackerBarrelHits = vtx_b_reco.outputHitCollection,
# ITrackerEndcapHits = vtx_ec_reco.outputHitCollection,
# OTrackerBarrelHits = trk_b_reco.outputHitCollection,
# OTrackerEndcapHits = trk_ec_reco.outputHitCollection,
# outputSourceLinks="TrackerMeasurements",
# OutputLevel=DEBUG)
## Track param init
truth_trk_init = TrackParamTruthInit("truth_trk_init",
inputMCParticles="mcparticles",
outputInitialTrackParameters="InitTrackParams",
OutputLevel=DEBUG)
#clust_trk_init = TrackParamClusterInit("clust_trk_init",
# inputClusters="SimpleClusters",
# outputInitialTrackParameters="InitTrackParamsFromClusters",
# OutputLevel=DEBUG)
#vtxcluster_trk_init = TrackParamVertexClusterInit("vtxcluster_trk_init",
# inputVertexHits="VertexBarrelRecHits",
# inputClusters="SimpleClusters",
# outputInitialTrackParameters="InitTrackParamsFromVtxClusters",
# maxHitRadius=40.0*units.mm,
# OutputLevel=WARNING)
# Tracking algorithms
trk_find_alg = TrackFindingAlgorithm("trk_find_alg",
inputSourceLinks = sourcelinker.outputSourceLinks,
inputMeasurements = sourcelinker.outputMeasurements,
inputInitialTrackParameters= "InitTrackParams",#"InitTrackParamsFromClusters",
outputTrajectories="trajectories",
OutputLevel=DEBUG)
parts_from_fit = ParticlesFromTrackFit("parts_from_fit",
inputTrajectories="trajectories",
outputParticles="ReconstructedParticles",
outputTrackParameters="outputTrackParameters",
OutputLevel=DEBUG)
#trk_find_alg1 = TrackFindingAlgorithm("trk_find_alg1",
# inputSourceLinks = sourcelinker.outputSourceLinks,
# inputMeasurements = sourcelinker.outputMeasurements,
# inputInitialTrackParameters= "InitTrackParamsFromClusters",
# outputTrajectories="trajectories1",
# OutputLevel=DEBUG)
#parts_from_fit1 = ParticlesFromTrackFit("parts_from_fit1",
# inputTrajectories="trajectories1",
# outputParticles="ReconstructedParticles1",
# outputTrackParameters="outputTrackParameters1",
# OutputLevel=DEBUG)
# Crystal Endcap Ecal
ce_ecal_daq = dict(
dynamicRangeADC=5.*GeV,
capacityADC=32768,
pedestalMean=400,
pedestalSigma=3)
ce_ecal_digi = CalHitDigi("ce_ecal_digi",
inputHitCollection="EcalEndcapNHits",
outputHitCollection="EcalEndcapNHitsDigi",
energyResolutions=[0., 0.02, 0.],
**ce_ecal_daq)
ce_ecal_reco = CalHitReco("ce_ecal_reco",
inputHitCollection="EcalEndcapNHitsDigi",
outputHitCollection="EcalEndcapNHitsReco",
thresholdFactor=4, # 4 sigma cut on pedestal sigma
readoutClass="EcalEndcapNHits",
sectorField="sector",
**ce_ecal_daq)
ce_ecal_cl = IslandCluster("ce_ecal_cl",
# OutputLevel=DEBUG,
inputHitCollection="EcalEndcapNHitsReco",
outputHitCollection="EcalEndcapNClusterHits",
splitCluster=False,
minClusterHitEdep=1.0*MeV, # discard low energy hits
minClusterCenterEdep=30*MeV,
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",
outputClusterCollection="EcalEndcapNClusters",
samplingFraction=0.998, # this accounts for a small fraction of leakage
logWeightBase=4.6)
# Endcap Sampling Ecal
ci_ecal_daq = dict(
dynamicRangeADC=50.*MeV,
capacityADC=32768,
pedestalMean=400,
pedestalSigma=10)
ci_ecal_digi = CalHitDigi("ci_ecal_digi",
inputHitCollection="EcalEndcapPHits",
outputHitCollection="EcalEndcapPHitsDigi",
**ci_ecal_daq)
ci_ecal_reco = CalHitReco("ci_ecal_reco",
inputHitCollection="EcalEndcapPHitsDigi",
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",
# OutputLevel=DEBUG,
inputHitCollection="EcalEndcapPHitsReco",
outputHitCollection="EcalEndcapPHitsRecoXY",
fields=["layer", "slice"],
fieldRefNumbers=[1, 0],
readoutClass="EcalEndcapPHits")
ci_ecal_cl = IslandCluster("ci_ecal_cl",
# OutputLevel=DEBUG,
inputHitCollection="EcalEndcapPHitsRecoXY",
outputHitCollection="EcalEndcapPClusterHits",
splitCluster=False,
minClusterCenterEdep=10.*MeV,
localDistXY=[10*mm, 10*mm])
ci_ecal_clreco = RecoCoG("ci_ecal_clreco",
inputHitCollection="EcalEndcapPClusterHits",
outputClusterCollection="EcalEndcapPClusters",
logWeightBase=6.2,
samplingFraction=ci_ecal_sf)
# Central Barrel Ecal (Imaging Cal.)
cb_ecal_daq = dict(
dynamicRangeADC=3*MeV,
capacityADC=8192,
pedestalMean=400,
pedestalSigma=20) # about 6 keV
cb_ecal_digi = CalHitDigi("cb_ecal_digi",
inputHitCollection="EcalBarrelHits",
outputHitCollection="EcalBarrelHitsDigi",
energyResolutions=[0., 0.02, 0.], # 2% flat resolution
**cb_ecal_daq)
cb_ecal_reco = ImCalPixelReco("cb_ecal_reco",
inputHitCollection="EcalBarrelHitsDigi",
outputHitCollection="EcalBarrelHitsReco",
thresholdFactor=3, # about 20 keV
readoutClass="EcalBarrelHits", # readout class
layerField="layer", # field to get layer id
sectorField="module", # field to get sector id
**cb_ecal_daq)
cb_ecal_cl = ImagingCluster("cb_ecal_cl",
inputHitCollection="EcalBarrelHitsReco",
outputHitCollection="EcalBarrelClusterHits",
localDistXY=[2.*mm, 2*mm], # same layer
layerDistEtaPhi=[10*mrad, 10*mrad], # adjacent layer
neighbourLayersRange=2, # id diff for adjacent layer
sectorDist=3.*cm) # different sector
cb_ecal_clreco = ImagingClusterReco("cb_ecal_clreco",
samplingFraction=cb_ecal_sf,
inputHitCollection="EcalBarrelClusterHits",
outputClusterCollection="EcalBarrelClusters",
outputLayerCollection="EcalBarrelLayers")
#Central ECAL SciFi
# 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",
layerField="layer",
sectorField="module",
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= scifi_barrel_sf)
# Central Barrel Hcal
cb_hcal_daq = dict(
dynamicRangeADC=50.*MeV,
capacityADC=32768,
pedestalMean=400,
pedestalSigma=10)
cb_hcal_digi = CalHitDigi("cb_hcal_digi",
inputHitCollection="HcalBarrelHits",
outputHitCollection="HcalBarrelHitsDigi",
**cb_hcal_daq)
cb_hcal_reco = CalHitReco("cb_hcal_reco",
inputHitCollection="HcalBarrelHitsDigi",
outputHitCollection="HcalBarrelHitsReco",
thresholdFactor=5.0,
readoutClass="HcalBarrelHits",
layerField="layer",
sectorField="module",
**cb_hcal_daq)
cb_hcal_merger = CalHitsMerger("cb_hcal_merger",
inputHitCollection="HcalBarrelHitsReco",
outputHitCollection="HcalBarrelHitsRecoXY",
readoutClass="HcalBarrelHits",
fields=["layer", "slice"],
fieldRefNumbers=[1, 0])
cb_hcal_cl = IslandCluster("cb_hcal_cl",
inputHitCollection="HcalBarrelHitsRecoXY",
outputHitCollection="HcalBarrelClusterHits",
splitCluster=False,
minClusterCenterEdep=30.*MeV,
localDistXY=[15.*cm, 15.*cm])
cb_hcal_clreco = RecoCoG("cb_hcal_clreco",
inputHitCollection="HcalBarrelClusterHits",
outputClusterCollection="HcalBarrelClusters",
logWeightBase=6.2,
samplingFraction=cb_hcal_sf)
# Hcal Hadron Endcap
ci_hcal_daq = dict(
dynamicRangeADC=50.*MeV,
capacityADC=32768,
pedestalMean=400,
pedestalSigma=10)
ci_hcal_digi = CalHitDigi("ci_hcal_digi",
inputHitCollection="HcalHadronEndcapHits",
outputHitCollection="HcalHadronEndcapHitsDigi",
**ci_hcal_daq)
ci_hcal_reco = CalHitReco("ci_hcal_reco",
inputHitCollection="HcalHadronEndcapHitsDigi",
outputHitCollection="HcalHadronEndcapHitsReco",
thresholdFactor=5.0,
**ci_hcal_daq)
ci_hcal_merger = CalHitsMerger("ci_hcal_merger",
inputHitCollection="HcalHadronEndcapHitsReco",
outputHitCollection="HcalHadronEndcapHitsRecoXY",
readoutClass="HcalHadronEndcapHits",
fields=["layer", "slice"],
fieldRefNumbers=[1, 0])
ci_hcal_cl = IslandCluster("ci_hcal_cl",
inputHitCollection="HcalHadronEndcapHitsRecoXY",
outputHitCollection="HcalHadronEndcapClusterHits",
splitCluster=False,
minClusterCenterEdep=30.*MeV,
localDistXY=[15.*cm, 15.*cm])
ci_hcal_clreco = RecoCoG("ci_hcal_clreco",
inputHitCollection="HcalHadronEndcapClusterHits",
outputClusterCollection="HcalHadronEndcapClusters",
logWeightBase=6.2,
samplingFraction=ci_hcal_sf)
# Hcal Electron Endcap
ce_hcal_daq = dict(
dynamicRangeADC=50.*MeV,
capacityADC=32768,
pedestalMean=400,
pedestalSigma=10)
ce_hcal_digi = CalHitDigi("ce_hcal_digi",
inputHitCollection="HcalElectronEndcapHits",
outputHitCollection="HcalElectronEndcapHitsDigi",
**ce_hcal_daq)
ce_hcal_reco = CalHitReco("ce_hcal_reco",
inputHitCollection="HcalElectronEndcapHitsDigi",
outputHitCollection="HcalElectronEndcapHitsReco",
thresholdFactor=5.0,
**ce_hcal_daq)
ce_hcal_merger = CalHitsMerger("ce_hcal_merger",
inputHitCollection="HcalElectronEndcapHitsReco",
outputHitCollection="HcalElectronEndcapHitsRecoXY",
readoutClass="HcalElectronEndcapHits",
fields=["layer", "slice"],
fieldRefNumbers=[1, 0])
ce_hcal_cl = IslandCluster("ce_hcal_cl",
inputHitCollection="HcalElectronEndcapHitsRecoXY",
outputHitCollection="HcalElectronEndcapClusterHits",
splitCluster=False,
minClusterCenterEdep=30.*MeV,
localDistXY=[15.*cm, 15.*cm])
ce_hcal_clreco = RecoCoG("ce_hcal_clreco",
inputHitCollection="HcalElectronEndcapClusterHits",
outputClusterCollection="HcalElectronEndcapClusters",
logWeightBase=6.2,
samplingFraction=ce_hcal_sf)
podout.outputCommands = ['drop *',
'keep mcparticles2',
'keep *Digi',
'keep *Reco*',
'keep *ClusterHits',
'keep *Clusters',
'keep *Layers',
'keep *outputTrackParameters*',
'keep *Track*']
ApplicationMgr(
TopAlg = [podin, copier,trkcopier,
trk_b_digi, trk_ec_digi, vtx_b_digi, vtx_ec_digi,
trk_b_reco, trk_ec_reco, vtx_b_reco, vtx_ec_reco,
sourcelinker,
#clust_trk_init,
truth_trk_init,
trk_find_alg, parts_from_fit,
#trk_find_alg1, parts_from_fit1,
ce_ecal_digi, ce_ecal_reco, ce_ecal_cl, ce_ecal_clreco,
ci_ecal_digi, ci_ecal_reco, ci_ecal_merger, ci_ecal_cl, ci_ecal_clreco,
cb_ecal_digi, cb_ecal_reco, cb_ecal_cl, cb_ecal_clreco,
scfi_barrel_digi, scfi_barrel_reco, scfi_barrel_merger, scfi_barrel_cl, scfi_barrel_clreco,
cb_hcal_digi, cb_hcal_reco, cb_hcal_merger, cb_hcal_cl, cb_hcal_clreco,
ce_hcal_digi, ce_hcal_reco, ce_hcal_merger, ce_hcal_cl, ce_hcal_clreco,
ci_hcal_digi, ci_hcal_reco, ci_hcal_merger, ci_hcal_cl, ci_hcal_clreco,
podout],
EvtSel = 'NONE',
EvtMax = n_events,
ExtSvc = [podioevent],
OutputLevel=DEBUG
)
benchmarks/clustering/scripts/cluster_plots.py
View file @ 22ed77c9
...@@ -25,6 +25,7 @@ def load_root_macros(arg_macros): ...@@ -25,6 +25,7 @@ def load_root_macros(arg_macros):
def flatten_collection(rdf, collection, cols=None): def flatten_collection(rdf, collection, cols=None):
if not cols: 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: else:
cols = ['{}.{}'.format(collection, c) for c in cols] cols = ['{}.{}'.format(collection, c) for c in cols]
if not cols: if not cols:
...@@ -50,8 +51,9 @@ def flatten_collection(rdf, collection, cols=None): ...@@ -50,8 +51,9 @@ def flatten_collection(rdf, collection, cols=None):
def thrown_particles_figure(rdf, save, mcbranch="mcparticles"): def thrown_particles_figure(rdf, save, mcbranch="mcparticles"):
# define truth particle info # define truth particle info
dft = flatten_collection(rdf, mcbranch, ['genStatus', 'pdgID', 'psx', 'psy', 'psz', 'mass']) dft = flatten_collection(rdf, mcbranch, ['genStatus', 'pdgID', 'ps.x', 'ps.y', 'ps.z', 'mass'])
dft.rename(columns={c: c.replace(mcbranch + '.', '') for c in dft.columns}, inplace=True) 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 # select thrown particles
dft = dft[dft['genStatus'] == 1] dft = dft[dft['genStatus'] == 1]
...@@ -77,7 +79,7 @@ def thrown_particles_figure(rdf, save, mcbranch="mcparticles"): ...@@ -77,7 +79,7 @@ def thrown_particles_figure(rdf, save, mcbranch="mcparticles"):
# calculate kinematics # calculate kinematics
get_4vecs = np.vectorize(lambda x, y, z, m: ROOT.Math.PxPyPzMVector(x, y, z, m)) get_4vecs = np.vectorize(lambda x, y, z, m: ROOT.Math.PxPyPzMVector(x, y, z, m))
fourvecs = get_4vecs(*dft[['psx', 'psy', 'psz', 'mass']].values.T) fourvecs = get_4vecs(*dft[['ps.x', 'ps.y', 'ps.z', 'mass']].values.T)
dft.loc[:, 'p'] = [v.P() for v in fourvecs] dft.loc[:, 'p'] = [v.P() for v in fourvecs]
dft.loc[:, 'eta'] = [v.Eta() for v in fourvecs] dft.loc[:, 'eta'] = [v.Eta() for v in fourvecs]
...@@ -152,17 +154,16 @@ if __name__ == '__main__': ...@@ -152,17 +154,16 @@ if __name__ == '__main__':
('eta', '$\eta$', 50), ('eta', '$\eta$', 50),
] ]
for coll in colls: for coll in colls:
print(coll)
df = flatten_collection(rdf_rec, coll) df = flatten_collection(rdf_rec, coll)
if not df.shape[0]: if not df.shape[0]:
print('{} do not have valid entries, skip it'.format(coll)) print('{} do not have valid entries, skip it'.format(coll))
continue continue
df.rename(columns={c: c.replace(coll + '.', '') for c in df.columns}, inplace=True) df.rename(columns={c: c.replace(coll + '.', '') for c in df.columns}, inplace=True)
# calculate eta df.rename(columns={c: c.replace(coll + 'Info.', '') for c in df.columns}, inplace=True)
if 'eta' not in df.columns:
df.loc[:, 'eta'] = -np.log(np.tan(df['polar.theta'].values/2.))
# print(df[['eta', 'polar.theta', 'position.x', 'position.y', 'position.z']]) # print(df[['eta', 'polar.theta', 'position.x', 'position.y', 'position.z']])
fig, axs = plt.subplots(2, 2, figsize=(12, 8), dpi=160) 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]), 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') bins=np.arange(0, 10), align='mid', ec='k')
axs[0][0].set_xlabel('Number of Clusters', fontsize=16) axs[0][0].set_xlabel('Number of Clusters', fontsize=16)
......
benchmarks/ecal/options/barrel.py
View file @ 22ed77c9
...@@ -2,6 +2,7 @@ ...@@ -2,6 +2,7 @@
An example script to digitize/reconstruct/clustering endcap ecal hits An example script to digitize/reconstruct/clustering endcap ecal hits
''' '''
from Gaudi.Configuration import * from Gaudi.Configuration import *
import json
import os import os
import ROOT import ROOT
...@@ -16,7 +17,12 @@ compact_path = str(os.environ.get("JUGGLER_COMPACT_PATH", "{}.xml".format(os.pat ...@@ -16,7 +17,12 @@ compact_path = str(os.environ.get("JUGGLER_COMPACT_PATH", "{}.xml".format(os.pat
input_sims = [f.strip() for f in str.split(os.environ["JUGGLER_SIM_FILE"], ",") if f.strip()] input_sims = [f.strip() for f in str.split(os.environ["JUGGLER_SIM_FILE"], ",") if f.strip()]
output_rec = str(os.environ["JUGGLER_REC_FILE"]) output_rec = str(os.environ["JUGGLER_REC_FILE"])
n_events = int(os.environ["JUGGLER_N_EVENTS"]) n_events = int(os.environ["JUGGLER_N_EVENTS"])
cb_ecal_sf = float(os.environ.get("CB_ECAL_SAMP_FRAC", 0.01324))
# input arguments from calibration file
with open('config/emcal_barrel_calibration.json') as f:
calib_data = json.load(f)['electron']
cb_ecal_sf = float(calib_data['sampling_fraction_img'])
# geometry service # geometry service
geo_service = GeoSvc("GeoSvc", detectors=[compact_path], OutputLevel=INFO) geo_service = GeoSvc("GeoSvc", detectors=[compact_path], OutputLevel=INFO)
...@@ -58,7 +64,7 @@ cb_ecal_digi = CalHitDigi("cb_ecal_digi", ...@@ -58,7 +64,7 @@ cb_ecal_digi = CalHitDigi("cb_ecal_digi",
**cb_ecal_daq) **cb_ecal_daq)
cb_ecal_reco = ImCalPixelReco("cb_ecal_reco", cb_ecal_reco = ImCalPixelReco("cb_ecal_reco",
inputHitCollection="EcalBarrelHitsDigi", inputHitCollection=cb_ecal_digi.outputHitCollection,
outputHitCollection="EcalBarrelHitsReco", outputHitCollection="EcalBarrelHitsReco",
thresholdFactor=3, # about 20 keV thresholdFactor=3, # about 20 keV
readoutClass="EcalBarrelHits", # readout class readoutClass="EcalBarrelHits", # readout class
...@@ -67,8 +73,8 @@ cb_ecal_reco = ImCalPixelReco("cb_ecal_reco", ...@@ -67,8 +73,8 @@ cb_ecal_reco = ImCalPixelReco("cb_ecal_reco",
**cb_ecal_daq) **cb_ecal_daq)
cb_ecal_cl = ImagingCluster("cb_ecal_cl", cb_ecal_cl = ImagingCluster("cb_ecal_cl",
inputHitCollection="EcalBarrelHitsReco", inputHitCollection=cb_ecal_reco.outputHitCollection,
outputHitCollection="EcalBarrelClusterHits", outputProtoClusterCollection="EcalBarrelProtoClusters",
localDistXY=[2.*mm, 2*mm], # same layer localDistXY=[2.*mm, 2*mm], # same layer
layerDistEtaPhi=[10*mrad, 10*mrad], # adjacent layer layerDistEtaPhi=[10*mrad, 10*mrad], # adjacent layer
neighbourLayersRange=2, # id diff for adjacent layer neighbourLayersRange=2, # id diff for adjacent layer
...@@ -76,16 +82,17 @@ cb_ecal_cl = ImagingCluster("cb_ecal_cl", ...@@ -76,16 +82,17 @@ cb_ecal_cl = ImagingCluster("cb_ecal_cl",
cb_ecal_clreco = ImagingClusterReco("cb_ecal_clreco", cb_ecal_clreco = ImagingClusterReco("cb_ecal_clreco",
samplingFraction=cb_ecal_sf, samplingFraction=cb_ecal_sf,
inputHitCollection="EcalBarrelClusterHits", inputHitCollection=cb_ecal_cl.inputHitCollection,
inputProtoClusterCollection=cb_ecal_cl.outputProtoClusterCollection,
outputClusterCollection="EcalBarrelClusters", outputClusterCollection="EcalBarrelClusters",
outputLayerCollection="EcalBarrelLayers") outputLayerCollection="EcalBarrelLayers",
outputInfoCollection="EcalBarrelClustersInfo")
podout.outputCommands = ['drop *', podout.outputCommands = ['drop *',
'keep mcparticles2', 'keep mcparticles2',
'keep *HitsReco', 'keep *HitsReco',
'keep *HitsDigi', 'keep *HitsDigi',
'keep *ClusterHits', 'keep *Cluster*']
'keep *Clusters']
ApplicationMgr( ApplicationMgr(
TopAlg = [podin, copier, TopAlg = [podin, copier,
......
benchmarks/ecal/options/endcap_e.py
View file @ 22ed77c9
...@@ -57,7 +57,7 @@ ce_ecal_digi = CalHitDigi("ce_ecal_digi", ...@@ -57,7 +57,7 @@ ce_ecal_digi = CalHitDigi("ce_ecal_digi",
**ce_ecal_daq) **ce_ecal_daq)
ce_ecal_reco = CalHitReco("ce_ecal_reco", ce_ecal_reco = CalHitReco("ce_ecal_reco",
inputHitCollection="EcalEndcapNHitsDigi", inputHitCollection=ce_ecal_digi.outputHitCollection,
outputHitCollection="EcalEndcapNHitsReco", outputHitCollection="EcalEndcapNHitsReco",
thresholdFactor=4, # 4 sigma cut on pedestal sigma thresholdFactor=4, # 4 sigma cut on pedestal sigma
readoutClass="EcalEndcapNHits", readoutClass="EcalEndcapNHits",
...@@ -66,8 +66,8 @@ ce_ecal_reco = CalHitReco("ce_ecal_reco", ...@@ -66,8 +66,8 @@ ce_ecal_reco = CalHitReco("ce_ecal_reco",
ce_ecal_cl = IslandCluster("ce_ecal_cl", ce_ecal_cl = IslandCluster("ce_ecal_cl",
# OutputLevel=DEBUG, # OutputLevel=DEBUG,
inputHitCollection="EcalEndcapNHitsReco", inputHitCollection=ce_ecal_reco.outputHitCollection,
outputHitCollection="EcalEndcapNClusterHits", outputProtoClusterCollection="EcalEndcapNProtoClusters",
splitCluster=False, splitCluster=False,
minClusterHitEdep=1.0*MeV, # discard low energy hits minClusterHitEdep=1.0*MeV, # discard low energy hits
minClusterCenterEdep=30*MeV, minClusterCenterEdep=30*MeV,
...@@ -75,8 +75,10 @@ ce_ecal_cl = IslandCluster("ce_ecal_cl", ...@@ -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 dimScaledLocalDistXY=[1.8, 1.8]) # hybrid calorimeter with different dimensions, using a scaled dist
ce_ecal_clreco = RecoCoG("ce_ecal_clreco", ce_ecal_clreco = RecoCoG("ce_ecal_clreco",
inputHitCollection="EcalEndcapNClusterHits", inputHitCollection=ce_ecal_cl.inputHitCollection,
inputProtoClusterCollection=ce_ecal_cl.outputProtoClusterCollection,
outputClusterCollection="EcalEndcapNClusters", outputClusterCollection="EcalEndcapNClusters",
outputInfoCollection="EcalEndcapNClustersInfo",
samplingFraction=0.998, # this accounts for a small fraction of leakage samplingFraction=0.998, # this accounts for a small fraction of leakage
logWeightBase=4.6) logWeightBase=4.6)
...@@ -84,8 +86,7 @@ podout.outputCommands = ['drop *', ...@@ -84,8 +86,7 @@ podout.outputCommands = ['drop *',
'keep mcparticles2', 'keep mcparticles2',
'keep *HitsReco', 'keep *HitsReco',
'keep *HitsDigi', 'keep *HitsDigi',
'keep *ClusterHits', 'keep *Cluster*']
'keep *Clusters']
ApplicationMgr( ApplicationMgr(
TopAlg = [podin, copier, TopAlg = [podin, copier,
......
benchmarks/ecal/options/endcap_i.py
View file @ 22ed77c9
...@@ -59,7 +59,7 @@ ci_ecal_digi = CalHitDigi("ci_ecal_digi", ...@@ -59,7 +59,7 @@ ci_ecal_digi = CalHitDigi("ci_ecal_digi",
**ci_ecal_daq) **ci_ecal_daq)
ci_ecal_reco = CalHitReco("ci_ecal_reco", ci_ecal_reco = CalHitReco("ci_ecal_reco",
inputHitCollection="EcalEndcapPHitsDigi", inputHitCollection=ci_ecal_digi.outputHitCollection,
outputHitCollection="EcalEndcapPHitsReco", outputHitCollection="EcalEndcapPHitsReco",
thresholdFactor=5.0, thresholdFactor=5.0,
**ci_ecal_daq) **ci_ecal_daq)
...@@ -67,7 +67,7 @@ ci_ecal_reco = CalHitReco("ci_ecal_reco", ...@@ -67,7 +67,7 @@ ci_ecal_reco = CalHitReco("ci_ecal_reco",
# merge hits in different layer (projection to local x-y plane) # merge hits in different layer (projection to local x-y plane)
ci_ecal_merger = CalHitsMerger("ci_ecal_merger", ci_ecal_merger = CalHitsMerger("ci_ecal_merger",
# OutputLevel=DEBUG, # OutputLevel=DEBUG,
inputHitCollection="EcalEndcapPHitsReco", inputHitCollection=ci_ecal_reco.outputHitCollection,
outputHitCollection="EcalEndcapPHitsRecoXY", outputHitCollection="EcalEndcapPHitsRecoXY",
fields=["layer", "slice"], fields=["layer", "slice"],
fieldRefNumbers=[1, 0], fieldRefNumbers=[1, 0],
...@@ -75,15 +75,17 @@ ci_ecal_merger = CalHitsMerger("ci_ecal_merger", ...@@ -75,15 +75,17 @@ ci_ecal_merger = CalHitsMerger("ci_ecal_merger",
ci_ecal_cl = IslandCluster("ci_ecal_cl", ci_ecal_cl = IslandCluster("ci_ecal_cl",
# OutputLevel=DEBUG, # OutputLevel=DEBUG,
inputHitCollection="EcalEndcapPHitsRecoXY", inputHitCollection=ci_ecal_merger.outputHitCollection,
outputHitCollection="EcalEndcapPClusterHits", outputProtoClusterCollection="EcalEndcapPProtoClusters",
splitCluster=False, splitCluster=False,
minClusterCenterEdep=10.*MeV, minClusterCenterEdep=10.*MeV,
localDistXY=[10*mm, 10*mm]) localDistXY=[10*mm, 10*mm])
ci_ecal_clreco = RecoCoG("ci_ecal_clreco", ci_ecal_clreco = RecoCoG("ci_ecal_clreco",
inputHitCollection="EcalEndcapPClusterHits", inputHitCollection=ci_ecal_cl.inputHitCollection,
inputProtoClusterCollection=ci_ecal_cl.outputProtoClusterCollection,
outputClusterCollection="EcalEndcapPClusters", outputClusterCollection="EcalEndcapPClusters",
outputInfoCollection="EcalEndcapPClustersInfo",
logWeightBase=6.2, logWeightBase=6.2,
samplingFraction=ci_ecal_sf) samplingFraction=ci_ecal_sf)
...@@ -91,8 +93,7 @@ podout.outputCommands = ['drop *', ...@@ -91,8 +93,7 @@ podout.outputCommands = ['drop *',
'keep mcparticles2', 'keep mcparticles2',
'keep *HitsReco*', 'keep *HitsReco*',
'keep *HitsDigi', 'keep *HitsDigi',
'keep *ClusterHits', 'keep *Cluster*']
'keep *Clusters']
ApplicationMgr( ApplicationMgr(
TopAlg = [podin, copier, TopAlg = [podin, copier,
......
benchmarks/ecal/run_emcal_benchmarks.py
View file @ 22ed77c9
...@@ -15,19 +15,19 @@ import argparse ...@@ -15,19 +15,19 @@ import argparse
default_type = { default_type = {
'endcap_e': [ 'endcap_e': [
['endcap_e.py'], ['endcap_e.py'],
['draw_cluters.py'], ['draw_clusters.py'],
['EcalEndcapNClusters']], ['EcalEndcapNClusters']],
'endcap_i': [ 'endcap_i': [
['endcap_i.py'], ['endcap_i.py'],
['draw_cluters.py'], ['draw_clusters.py'],
['EcalEndcapPClusters']], ['EcalEndcapPClusters']],
'barrel': [ 'barrel': [
['barrel.py'], ['barrel.py'],
['draw_cluters.py'], ['draw_clusters.py'],
['EcalBarrelClusters']], ['EcalBarrelClusters']],
# 'all': [ # 'all': [
# ['all_ecal.py'], # ['all_ecal.py'],
# ['draw_cluters.py'], # ['draw_clusters.py'],
# ['EcalEndcapNClusters', 'EcalEndcapPClusters', 'EcalBarrelClusters']], # ['EcalEndcapNClusters', 'EcalEndcapPClusters', 'EcalBarrelClusters']],
} }
......
benchmarks/ecal/scripts/draw_cluters.py → benchmarks/ecal/scripts/draw_clusters.py
View file @ 22ed77c9
...@@ -25,7 +25,8 @@ def load_root_macros(arg_macros): ...@@ -25,7 +25,8 @@ def load_root_macros(arg_macros):
# read from RDataFrame and flatten a given collection, return pandas dataframe # read from RDataFrame and flatten a given collection, return pandas dataframe
def flatten_collection(rdf, collection, cols=None): def flatten_collection(rdf, collection, cols=None):
if not cols: 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: else:
cols = ['{}.{}'.format(collection, c) for c in cols] cols = ['{}.{}'.format(collection, c) for c in cols]
if not cols: if not cols:
...@@ -92,10 +93,10 @@ if __name__ == '__main__': ...@@ -92,10 +93,10 @@ if __name__ == '__main__':
] ]
for coll in [c.strip() for c in args.coll.split(',')]: for coll in [c.strip() for c in args.coll.split(',')]:
df = flatten_collection(rdf_rec, coll) df = flatten_collection(rdf_rec, coll)
# calculate eta # Get eta from info table
df.loc[:, '{}.eta'.format(coll)] = -np.log(np.tan(df['{}.polar.theta'.format(coll)].values/2.)) df[coll + '.eta'] = df[coll + 'Info.eta']
fig, axs = plt.subplots(2, 2, figsize=(12, 8), dpi=160) 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]), 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') bins=np.arange(0, 10), align='mid', ec='k')
axs[0][0].set_xlabel('Number of Clusters', fontsize=16) axs[0][0].set_xlabel('Number of Clusters', fontsize=16)
......
benchmarks/full/options/full_reconstruction.py
View file @ 22ed77c9
...@@ -4,6 +4,8 @@ from Configurables import ApplicationMgr, EICDataSvc, PodioOutput, GeoSvc ...@@ -4,6 +4,8 @@ from Configurables import ApplicationMgr, EICDataSvc, PodioOutput, GeoSvc
from GaudiKernel import SystemOfUnits as units from GaudiKernel import SystemOfUnits as units
from GaudiKernel.SystemOfUnits import MeV, GeV, mm, cm, mrad from GaudiKernel.SystemOfUnits import MeV, GeV, mm, cm, mrad
import json
detector_name = "athena" detector_name = "athena"
if "JUGGLER_DETECTOR" in os.environ : if "JUGGLER_DETECTOR" in os.environ :
detector_name = str(os.environ["JUGGLER_DETECTOR"]) detector_name = str(os.environ["JUGGLER_DETECTOR"])
...@@ -18,13 +20,20 @@ compact_path = os.path.join(detector_path, detector_name) ...@@ -18,13 +20,20 @@ compact_path = os.path.join(detector_path, detector_name)
qe_data = [(1.0, 0.25), (7.5, 0.25),] qe_data = [(1.0, 0.25), (7.5, 0.25),]
# CAL reconstruction # CAL reconstruction
# get sampling fractions from system environment variable, 1.0 by default # get sampling fractions from system environment variable
ci_ecal_sf = float(os.environ.get("CI_ECAL_SAMP_FRAC", 0.253)) ci_ecal_sf = float(os.environ.get("CI_ECAL_SAMP_FRAC", 0.253))
cb_ecal_sf = float(os.environ.get("CB_ECAL_SAMP_FRAC", 0.01324))
cb_hcal_sf = float(os.environ.get("CB_HCAL_SAMP_FRAC", 0.038)) cb_hcal_sf = float(os.environ.get("CB_HCAL_SAMP_FRAC", 0.038))
ci_hcal_sf = float(os.environ.get("CI_HCAL_SAMP_FRAC", 0.025)) ci_hcal_sf = float(os.environ.get("CI_HCAL_SAMP_FRAC", 0.025))
ce_hcal_sf = float(os.environ.get("CE_HCAL_SAMP_FRAC", 0.025)) ce_hcal_sf = float(os.environ.get("CE_HCAL_SAMP_FRAC", 0.025))
scifi_barrel_sf = float(os.environ.get("CB_EMCAL_SCFI_SAMP_FRAC",0.0938))
# input arguments from calibration file
with open('config/emcal_barrel_calibration.json') as f:
calib_data = json.load(f)['electron']
print(calib_data)
cb_ecal_sf = float(calib_data['sampling_fraction_img'])
scifi_barrel_sf = float(calib_data['sampling_fraction_scfi'])
# input and output # input and output
input_sims = [f.strip() for f in str.split(os.environ["JUGGLER_SIM_FILE"], ",") if f.strip()] input_sims = [f.strip() for f in str.split(os.environ["JUGGLER_SIM_FILE"], ",") if f.strip()]
...@@ -86,13 +95,14 @@ sim_coll = [ ...@@ -86,13 +95,14 @@ sim_coll = [
"EcalBarrelHits", "EcalBarrelHits",
"EcalBarrelScFiHits", "EcalBarrelScFiHits",
"HcalBarrelHits", "HcalBarrelHits",
"HcalHadronEndcapHits", "HcalEndcapPHits",
"HcalElectronEndcapHits", "HcalEndcapNHits",
"TrackerEndcapHits", "TrackerEndcapHits",
"TrackerBarrelHits", "TrackerBarrelHits",
"GEMTrackerEndcapHits",
"VertexBarrelHits", "VertexBarrelHits",
"VertexEndcapHits", "VertexEndcapHits",
"ForwardRICHHits" "DRICHHits"
] ]
# input and output # input and output
...@@ -107,8 +117,8 @@ trkcopier = TrkCopier("TrkCopier", ...@@ -107,8 +117,8 @@ trkcopier = TrkCopier("TrkCopier",
inputCollection="TrackerBarrelHits", inputCollection="TrackerBarrelHits",
outputCollection="TrackerBarrelHits2") outputCollection="TrackerBarrelHits2")
pmtcopier = PMTCopier("PMTCopier", pmtcopier = PMTCopier("PMTCopier",
inputCollection="ForwardRICHHits", inputCollection="DRICHHits",
outputCollection="ForwardRICHHits2") outputCollection="DRICHHits2")
# Dummy reconstruction # Dummy reconstruction
dummy = MC2DummyParticle("MC2Dummy", dummy = MC2DummyParticle("MC2Dummy",
...@@ -122,13 +132,14 @@ ecal_digi = EMCalorimeterDigi("ecal_digi", ...@@ -122,13 +132,14 @@ ecal_digi = EMCalorimeterDigi("ecal_digi",
outputHitCollection="EcalBarrelHitsSimpleDigi") outputHitCollection="EcalBarrelHitsSimpleDigi")
ecal_reco = EMCalReconstruction("ecal_reco", ecal_reco = EMCalReconstruction("ecal_reco",
inputHitCollection="EcalBarrelHitsSimpleDigi", inputHitCollection=ecal_digi.outputHitCollection,
outputHitCollection="EcalBarrelHitsSimpleReco", outputHitCollection="EcalBarrelHitsSimpleReco",
minModuleEdep=0.0*units.MeV) minModuleEdep=0.0*units.MeV)
simple_cluster = SimpleClustering("simple_cluster", simple_cluster = SimpleClustering("simple_cluster",
inputHitCollection="EcalBarrelHitsSimpleReco", inputHitCollection=ecal_reco.outputHitCollection,
outputClusters="SimpleClusters", outputClusterCollection="SimpleClusters",
outputProtoClusterCollection="SimpleProtoClusters",
minModuleEdep=1.0*units.MeV, minModuleEdep=1.0*units.MeV,
maxDistance=50.0*units.cm) maxDistance=50.0*units.cm)
...@@ -146,7 +157,7 @@ ce_ecal_digi = CalHitDigi("ce_ecal_digi", ...@@ -146,7 +157,7 @@ ce_ecal_digi = CalHitDigi("ce_ecal_digi",
**ce_ecal_daq) **ce_ecal_daq)
ce_ecal_reco = CalHitReco("ce_ecal_reco", ce_ecal_reco = CalHitReco("ce_ecal_reco",
inputHitCollection="EcalEndcapNHitsDigi", inputHitCollection=ce_ecal_digi.outputHitCollection,
outputHitCollection="EcalEndcapNHitsReco", outputHitCollection="EcalEndcapNHitsReco",
thresholdFactor=4, # 4 sigma cut on pedestal sigma thresholdFactor=4, # 4 sigma cut on pedestal sigma
readoutClass="EcalEndcapNHits", readoutClass="EcalEndcapNHits",
...@@ -154,8 +165,8 @@ ce_ecal_reco = CalHitReco("ce_ecal_reco", ...@@ -154,8 +165,8 @@ ce_ecal_reco = CalHitReco("ce_ecal_reco",
**ce_ecal_daq) **ce_ecal_daq)
ce_ecal_cl = IslandCluster("ce_ecal_cl", ce_ecal_cl = IslandCluster("ce_ecal_cl",
inputHitCollection="EcalEndcapNHitsReco", inputHitCollection=ce_ecal_reco.outputHitCollection,
outputHitCollection="EcalEndcapNClusterHits", outputProtoClusterCollection="EcalEndcapNProtoClusters",
splitCluster=False, splitCluster=False,
minClusterHitEdep=1.0*units.MeV, # discard low energy hits minClusterHitEdep=1.0*units.MeV, # discard low energy hits
minClusterCenterEdep=30*units.MeV, minClusterCenterEdep=30*units.MeV,
...@@ -163,7 +174,8 @@ ce_ecal_cl = IslandCluster("ce_ecal_cl", ...@@ -163,7 +174,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 dimScaledLocalDistXY=[1.8, 1.8]) # dimension scaled dist is good for hybrid sectors with different module size
ce_ecal_clreco = RecoCoG("ce_ecal_clreco", ce_ecal_clreco = RecoCoG("ce_ecal_clreco",
inputHitCollection="EcalEndcapNClusterHits", inputHitCollection=ce_ecal_cl.inputHitCollection,
inputProtoClusterCollection=ce_ecal_cl.outputProtoClusterCollection,
outputClusterCollection="EcalEndcapNClusters", outputClusterCollection="EcalEndcapNClusters",
samplingFraction=0.998, # this accounts for a small fraction of leakage samplingFraction=0.998, # this accounts for a small fraction of leakage
logWeightBase=4.6) logWeightBase=4.6)
...@@ -181,29 +193,31 @@ ci_ecal_digi = CalHitDigi("ci_ecal_digi", ...@@ -181,29 +193,31 @@ ci_ecal_digi = CalHitDigi("ci_ecal_digi",
**ci_ecal_daq) **ci_ecal_daq)
ci_ecal_reco = CalHitReco("ci_ecal_reco", ci_ecal_reco = CalHitReco("ci_ecal_reco",
inputHitCollection="EcalEndcapPHitsDigi", inputHitCollection=ci_ecal_digi.outputHitCollection,
outputHitCollection="EcalEndcapPHitsReco", outputHitCollection="EcalEndcapPHitsReco",
thresholdFactor=5.0, thresholdFactor=5.0,
**ci_ecal_daq) **ci_ecal_daq)
# merge hits in different layer (projection to local x-y plane) # merge hits in different layer (projection to local x-y plane)
ci_ecal_merger = CalHitsMerger("ci_ecal_merger", ci_ecal_merger = CalHitsMerger("ci_ecal_merger",
inputHitCollection="EcalEndcapPHitsReco", inputHitCollection=ci_ecal_reco.outputHitCollection,
outputHitCollection="EcalEndcapPHitsRecoXY", outputHitCollection="EcalEndcapPHitsRecoXY",
fields=["layer", "slice"], fields=["layer", "slice"],
fieldRefNumbers=[1, 0], fieldRefNumbers=[1, 0],
readoutClass="EcalEndcapPHits") readoutClass="EcalEndcapPHits")
ci_ecal_cl = IslandCluster("ci_ecal_cl", ci_ecal_cl = IslandCluster("ci_ecal_cl",
inputHitCollection="EcalEndcapPHitsRecoXY", inputHitCollection=ci_ecal_merger.outputHitCollection,
outputHitCollection="EcalEndcapPClusterHits", outputProtoClusterCollection="EcalEndcapPProtoClusters",
splitCluster=False, splitCluster=False,
minClusterCenterEdep=10.*units.MeV, minClusterCenterEdep=10.*units.MeV,
localDistXY=[10*units.mm, 10*units.mm]) localDistXY=[10*units.mm, 10*units.mm])
ci_ecal_clreco = RecoCoG("ci_ecal_clreco", ci_ecal_clreco = RecoCoG("ci_ecal_clreco",
inputHitCollection="EcalEndcapPClusterHits", inputHitCollection=ci_ecal_cl.inputHitCollection,
inputProtoClusterCollection=ci_ecal_cl.outputProtoClusterCollection,
outputClusterCollection="EcalEndcapPClusters", outputClusterCollection="EcalEndcapPClusters",
outputInfoCollection="EcalEndcapPClustersInfo",
logWeightBase=6.2, logWeightBase=6.2,
samplingFraction=ci_ecal_sf) samplingFraction=ci_ecal_sf)
...@@ -221,7 +235,7 @@ cb_ecal_digi = CalHitDigi("cb_ecal_digi", ...@@ -221,7 +235,7 @@ cb_ecal_digi = CalHitDigi("cb_ecal_digi",
**cb_ecal_daq) **cb_ecal_daq)
cb_ecal_reco = ImCalPixelReco("cb_ecal_reco", cb_ecal_reco = ImCalPixelReco("cb_ecal_reco",
inputHitCollection="EcalBarrelHitsDigi", inputHitCollection=cb_ecal_digi.outputHitCollection,
outputHitCollection="EcalBarrelHitsReco", outputHitCollection="EcalBarrelHitsReco",
thresholdFactor=3, # about 20 keV thresholdFactor=3, # about 20 keV
readoutClass="EcalBarrelHits", # readout class readoutClass="EcalBarrelHits", # readout class
...@@ -230,8 +244,8 @@ cb_ecal_reco = ImCalPixelReco("cb_ecal_reco", ...@@ -230,8 +244,8 @@ cb_ecal_reco = ImCalPixelReco("cb_ecal_reco",
**cb_ecal_daq) **cb_ecal_daq)
cb_ecal_cl = ImagingCluster("cb_ecal_cl", cb_ecal_cl = ImagingCluster("cb_ecal_cl",
inputHitCollection="EcalBarrelHitsReco", inputHitCollection=cb_ecal_reco.outputHitCollection,
outputHitCollection="EcalBarrelClusterHits", outputProtoClusterCollection="EcalBarrelProtoClusters",
localDistXY=[2.*units.mm, 2*units.mm], # same layer localDistXY=[2.*units.mm, 2*units.mm], # same layer
layerDistEtaPhi=[10*units.mrad, 10*units.mrad], # adjacent layer layerDistEtaPhi=[10*units.mrad, 10*units.mrad], # adjacent layer
neighbourLayersRange=2, # id diff for adjacent layer neighbourLayersRange=2, # id diff for adjacent layer
...@@ -239,8 +253,10 @@ cb_ecal_cl = ImagingCluster("cb_ecal_cl", ...@@ -239,8 +253,10 @@ cb_ecal_cl = ImagingCluster("cb_ecal_cl",
cb_ecal_clreco = ImagingClusterReco("cb_ecal_clreco", cb_ecal_clreco = ImagingClusterReco("cb_ecal_clreco",
samplingFraction=cb_ecal_sf, samplingFraction=cb_ecal_sf,
inputHitCollection="EcalBarrelClusterHits", inputHitCollection=cb_ecal_cl.inputHitCollection,
inputProtoClusterCollection=cb_ecal_cl.outputProtoClusterCollection,
outputClusterCollection="EcalBarrelClusters", outputClusterCollection="EcalBarrelClusters",
outputInfoCollection="EcalBarrelClustersInfo",
outputLayerCollection="EcalBarrelLayers") outputLayerCollection="EcalBarrelLayers")
# Central ECAL SciFi # Central ECAL SciFi
...@@ -256,7 +272,7 @@ scfi_barrel_digi = CalHitDigi("scfi_barrel_digi", ...@@ -256,7 +272,7 @@ scfi_barrel_digi = CalHitDigi("scfi_barrel_digi",
**scfi_barrel_daq) **scfi_barrel_daq)
scfi_barrel_reco = CalHitReco("scfi_barrel_reco", scfi_barrel_reco = CalHitReco("scfi_barrel_reco",
inputHitCollection="EcalBarrelScFiHitsDigi", inputHitCollection=scfi_barrel_digi.outputHitCollection,
outputHitCollection="EcalBarrelScFiHitsReco", outputHitCollection="EcalBarrelScFiHitsReco",
thresholdFactor=5.0, thresholdFactor=5.0,
readoutClass="EcalBarrelScFiHits", readoutClass="EcalBarrelScFiHits",
...@@ -267,22 +283,24 @@ scfi_barrel_reco = CalHitReco("scfi_barrel_reco", ...@@ -267,22 +283,24 @@ scfi_barrel_reco = CalHitReco("scfi_barrel_reco",
# merge hits in different layer (projection to local x-y plane) # merge hits in different layer (projection to local x-y plane)
scfi_barrel_merger = CalHitsMerger("scfi_barrel_merger", scfi_barrel_merger = CalHitsMerger("scfi_barrel_merger",
inputHitCollection="EcalBarrelScFiHitsReco", inputHitCollection=scfi_barrel_reco.outputHitCollection,
outputHitCollection="EcalBarrelScFiGridReco", outputHitCollection="EcalBarrelScFiGridReco",
fields=["fiber"], fields=["fiber"],
fieldRefNumbers=[1], fieldRefNumbers=[1],
readoutClass="EcalBarrelScFiHits") readoutClass="EcalBarrelScFiHits")
scfi_barrel_cl = IslandCluster("scfi_barrel_cl", scfi_barrel_cl = IslandCluster("scfi_barrel_cl",
inputHitCollection="EcalBarrelScFiGridReco", inputHitCollection=scfi_barrel_merger.outputHitCollection,
outputHitCollection="EcalBarrelScFiClusterHits", outputProtoClusterCollection="EcalBarrelScFiProtoClusters",
splitCluster=False, splitCluster=False,
minClusterCenterEdep=10.*MeV, minClusterCenterEdep=10.*MeV,
localDistXZ=[30*mm, 30*mm]) localDistXZ=[30*mm, 30*mm])
scfi_barrel_clreco = RecoCoG("scfi_barrel_clreco", scfi_barrel_clreco = RecoCoG("scfi_barrel_clreco",
inputHitCollection="EcalBarrelScFiClusterHits", inputHitCollection=scfi_barrel_cl.inputHitCollection,
inputProtoClusterCollection=scfi_barrel_cl.outputProtoClusterCollection,
outputClusterCollection="EcalBarrelScFiClusters", outputClusterCollection="EcalBarrelScFiClusters",
outputInfoCollection="EcalBarrelScFiClustersInfo",
logWeightBase=6.2, logWeightBase=6.2,
samplingFraction= scifi_barrel_sf) samplingFraction= scifi_barrel_sf)
...@@ -299,7 +317,7 @@ cb_hcal_digi = CalHitDigi("cb_hcal_digi", ...@@ -299,7 +317,7 @@ cb_hcal_digi = CalHitDigi("cb_hcal_digi",
**cb_hcal_daq) **cb_hcal_daq)
cb_hcal_reco = CalHitReco("cb_hcal_reco", cb_hcal_reco = CalHitReco("cb_hcal_reco",
inputHitCollection="HcalBarrelHitsDigi", inputHitCollection=cb_hcal_digi.outputHitCollection,
outputHitCollection="HcalBarrelHitsReco", outputHitCollection="HcalBarrelHitsReco",
thresholdFactor=5.0, thresholdFactor=5.0,
readoutClass="HcalBarrelHits", readoutClass="HcalBarrelHits",
...@@ -308,22 +326,24 @@ cb_hcal_reco = CalHitReco("cb_hcal_reco", ...@@ -308,22 +326,24 @@ cb_hcal_reco = CalHitReco("cb_hcal_reco",
**cb_hcal_daq) **cb_hcal_daq)
cb_hcal_merger = CalHitsMerger("cb_hcal_merger", cb_hcal_merger = CalHitsMerger("cb_hcal_merger",
inputHitCollection="HcalBarrelHitsReco", inputHitCollection=cb_hcal_reco.outputHitCollection,
outputHitCollection="HcalBarrelHitsRecoXY", outputHitCollection="HcalBarrelHitsRecoXY",
readoutClass="HcalBarrelHits", readoutClass="HcalBarrelHits",
fields=["layer", "slice"], fields=["layer", "slice"],
fieldRefNumbers=[1, 0]) fieldRefNumbers=[1, 0])
cb_hcal_cl = IslandCluster("cb_hcal_cl", cb_hcal_cl = IslandCluster("cb_hcal_cl",
inputHitCollection="HcalBarrelHitsRecoXY", inputHitCollection=cb_hcal_merger.outputHitCollection,
outputHitCollection="HcalBarrelClusterHits", outputProtoClusterCollection="HcalBarrelProtoClusters",
splitCluster=False, splitCluster=False,
minClusterCenterEdep=30.*units.MeV, minClusterCenterEdep=30.*units.MeV,
localDistXY=[15.*units.cm, 15.*units.cm]) localDistXY=[15.*units.cm, 15.*units.cm])
cb_hcal_clreco = RecoCoG("cb_hcal_clreco", cb_hcal_clreco = RecoCoG("cb_hcal_clreco",
inputHitCollection="HcalBarrelClusterHits", inputHitCollection=cb_hcal_cl.inputHitCollection,
inputProtoClusterCollection=cb_hcal_cl.outputProtoClusterCollection,
outputClusterCollection="HcalBarrelClusters", outputClusterCollection="HcalBarrelClusters",
outputInfoCollection="HcalBarrelClustersInfo",
logWeightBase=6.2, logWeightBase=6.2,
samplingFraction=cb_hcal_sf) samplingFraction=cb_hcal_sf)
...@@ -335,33 +355,35 @@ ci_hcal_daq = dict( ...@@ -335,33 +355,35 @@ ci_hcal_daq = dict(
pedestalMean=400, pedestalMean=400,
pedestalSigma=10) pedestalSigma=10)
ci_hcal_digi = CalHitDigi("ci_hcal_digi", ci_hcal_digi = CalHitDigi("ci_hcal_digi",
inputHitCollection="HcalHadronEndcapHits", inputHitCollection="HcalEndcapPHits",
outputHitCollection="HcalHadronEndcapHitsDigi", outputHitCollection="HcalEndcapPHitsDigi",
**ci_hcal_daq) **ci_hcal_daq)
ci_hcal_reco = CalHitReco("ci_hcal_reco", ci_hcal_reco = CalHitReco("ci_hcal_reco",
inputHitCollection="HcalHadronEndcapHitsDigi", inputHitCollection=ci_hcal_digi.outputHitCollection,
outputHitCollection="HcalHadronEndcapHitsReco", outputHitCollection="HcalEndcapPHitsReco",
thresholdFactor=5.0, thresholdFactor=5.0,
**ci_hcal_daq) **ci_hcal_daq)
ci_hcal_merger = CalHitsMerger("ci_hcal_merger", ci_hcal_merger = CalHitsMerger("ci_hcal_merger",
inputHitCollection="HcalHadronEndcapHitsReco", inputHitCollection=ci_hcal_reco.outputHitCollection,
outputHitCollection="HcalHadronEndcapHitsRecoXY", outputHitCollection="HcalEndcapPHitsRecoXY",
readoutClass="HcalHadronEndcapHits", readoutClass="HcalEndcapPHits",
fields=["layer", "slice"], fields=["layer", "slice"],
fieldRefNumbers=[1, 0]) fieldRefNumbers=[1, 0])
ci_hcal_cl = IslandCluster("ci_hcal_cl", ci_hcal_cl = IslandCluster("ci_hcal_cl",
inputHitCollection="HcalHadronEndcapHitsRecoXY", inputHitCollection=ci_hcal_merger.outputHitCollection,
outputHitCollection="HcalHadronEndcapClusterHits", outputProtoClusterCollection="HcalEndcapPProtoClusters",
splitCluster=False, splitCluster=False,
minClusterCenterEdep=30.*units.MeV, minClusterCenterEdep=30.*units.MeV,
localDistXY=[15.*units.cm, 15.*units.cm]) localDistXY=[15.*units.cm, 15.*units.cm])
ci_hcal_clreco = RecoCoG("ci_hcal_clreco", ci_hcal_clreco = RecoCoG("ci_hcal_clreco",
inputHitCollection="HcalHadronEndcapClusterHits", inputHitCollection=ci_hcal_cl.inputHitCollection,
outputClusterCollection="HcalHadronEndcapClusters", inputProtoClusterCollection=ci_hcal_cl.outputProtoClusterCollection,
outputClusterCollection="HcalEndcapPClusters",
outputInfoCollection="HcalEndcapPClustersInfo",
logWeightBase=6.2, logWeightBase=6.2,
samplingFraction=ci_hcal_sf) samplingFraction=ci_hcal_sf)
...@@ -373,33 +395,35 @@ ce_hcal_daq = dict( ...@@ -373,33 +395,35 @@ ce_hcal_daq = dict(
pedestalSigma=10) pedestalSigma=10)
ce_hcal_digi = CalHitDigi("ce_hcal_digi", ce_hcal_digi = CalHitDigi("ce_hcal_digi",
inputHitCollection="HcalElectronEndcapHits", inputHitCollection="HcalEndcapNHits",
outputHitCollection="HcalElectronEndcapHitsDigi", outputHitCollection="HcalEndcapNHitsDigi",
**ce_hcal_daq) **ce_hcal_daq)
ce_hcal_reco = CalHitReco("ce_hcal_reco", ce_hcal_reco = CalHitReco("ce_hcal_reco",
inputHitCollection="HcalElectronEndcapHitsDigi", inputHitCollection=ce_hcal_digi.outputHitCollection,
outputHitCollection="HcalElectronEndcapHitsReco", outputHitCollection="HcalEndcapNHitsReco",
thresholdFactor=5.0, thresholdFactor=5.0,
**ce_hcal_daq) **ce_hcal_daq)
ce_hcal_merger = CalHitsMerger("ce_hcal_merger", ce_hcal_merger = CalHitsMerger("ce_hcal_merger",
inputHitCollection="HcalElectronEndcapHitsReco", inputHitCollection=ce_hcal_reco.outputHitCollection,
outputHitCollection="HcalElectronEndcapHitsRecoXY", outputHitCollection="HcalEndcapNHitsRecoXY",
readoutClass="HcalElectronEndcapHits", readoutClass="HcalEndcapNHits",
fields=["layer", "slice"], fields=["layer", "slice"],
fieldRefNumbers=[1, 0]) fieldRefNumbers=[1, 0])
ce_hcal_cl = IslandCluster("ce_hcal_cl", ce_hcal_cl = IslandCluster("ce_hcal_cl",
inputHitCollection="HcalElectronEndcapHitsRecoXY", inputHitCollection=ce_hcal_merger.outputHitCollection,
outputHitCollection="HcalElectronEndcapClusterHits", outputProtoClusterCollection="HcalEndcapNProtoClusters",
splitCluster=False, splitCluster=False,
minClusterCenterEdep=30.*units.MeV, minClusterCenterEdep=30.*units.MeV,
localDistXY=[15.*units.cm, 15.*units.cm]) localDistXY=[15.*units.cm, 15.*units.cm])
ce_hcal_clreco = RecoCoG("ce_hcal_clreco", ce_hcal_clreco = RecoCoG("ce_hcal_clreco",
inputHitCollection="HcalElectronEndcapClusterHits", inputHitCollection=ce_hcal_cl.inputHitCollection,
outputClusterCollection="HcalElectronEndcapClusters", inputProtoClusterCollection=ce_hcal_cl.outputProtoClusterCollection,
outputClusterCollection="HcalEndcapNClusters",
outputInfoCollection="HcalEndcapNClustersInfo",
logWeightBase=6.2, logWeightBase=6.2,
samplingFraction=ce_hcal_sf) samplingFraction=ce_hcal_sf)
...@@ -424,6 +448,11 @@ vtx_ec_digi = TrackerDigi("vtx_ec_digi", ...@@ -424,6 +448,11 @@ vtx_ec_digi = TrackerDigi("vtx_ec_digi",
outputHitCollection="VertexEndcapRawHits", outputHitCollection="VertexEndcapRawHits",
timeResolution=8) timeResolution=8)
gem_ec_digi = TrackerDigi("gem_ec_digi",
inputHitCollection="GEMTrackerEndcapHits",
outputHitCollection="GEMTrackerEndcapRawHits",
timeResolution=10)
# Tracker and vertex reconstruction # Tracker and vertex reconstruction
trk_b_reco = TrackerHitReconstruction("trk_b_reco", trk_b_reco = TrackerHitReconstruction("trk_b_reco",
inputHitCollection = trk_b_digi.outputHitCollection, inputHitCollection = trk_b_digi.outputHitCollection,
...@@ -441,6 +470,11 @@ vtx_ec_reco = TrackerHitReconstruction("vtx_ec_reco", ...@@ -441,6 +470,11 @@ vtx_ec_reco = TrackerHitReconstruction("vtx_ec_reco",
inputHitCollection = vtx_ec_digi.outputHitCollection, inputHitCollection = vtx_ec_digi.outputHitCollection,
outputHitCollection="VertexEndcapRecHits") outputHitCollection="VertexEndcapRecHits")
gem_ec_reco = TrackerHitReconstruction("gem_ec_digi",
inputHitCollection=gem_ec_digi.outputHitCollection,
outputHitCollection="GEMTrackerEndcapRecHits",
timeResolution=10)
# Hit Source linker # Hit Source linker
sourcelinker = TrackerSourcesLinker("trk_srcslnkr", sourcelinker = TrackerSourcesLinker("trk_srcslnkr",
inputHitCollections = ["VertexBarrelRecHits", "TrackerBarrelRecHits"], inputHitCollections = ["VertexBarrelRecHits", "TrackerBarrelRecHits"],
...@@ -504,13 +538,13 @@ parts_from_fit = ParticlesFromTrackFit("parts_from_fit", ...@@ -504,13 +538,13 @@ parts_from_fit = ParticlesFromTrackFit("parts_from_fit",
pmtdigi = PhotoMultiplierDigi("pmtdigi", pmtdigi = PhotoMultiplierDigi("pmtdigi",
inputHitCollection="ForwardRICHHits", inputHitCollection="DRICHHits",
outputHitCollection="ForwardRICHHitsDigi", outputHitCollection="DRICHHitsDigi",
quantumEfficiency=[(a*units.eV, b) for a, b in qe_data]) quantumEfficiency=[(a*units.eV, b) for a, b in qe_data])
pmtreco = PhotoMultiplierReco("pmtreco", pmtreco = PhotoMultiplierReco("pmtreco",
inputHitCollection="ForwardRICHHitsDigi", inputHitCollection="DRICHHitsDigi",
outputHitCollection="ForwardRICHHitsReco") outputHitCollection="DRICHHitsReco")
# FIXME # FIXME
#richcluster = PhotoRingClusters("richcluster", #richcluster = PhotoRingClusters("richcluster",
......
benchmarks/imaging_ecal/options/hybrid_cluster.py
View file @ 22ed77c9
import json
import os import os
import ROOT import ROOT
from Gaudi.Configuration import * from Gaudi.Configuration import *
from GaudiKernel.SystemOfUnits import MeV, mm, cm, mrad, rad, ns from GaudiKernel.SystemOfUnits import MeV, mm, cm, mrad, rad, ns
from Configurables import ApplicationMgr, EICDataSvc, PodioOutput, GeoSvc from Configurables import ApplicationMgr, EICDataSvc, PodioOutput, GeoSvc
from Configurables import PodioInput from Configurables import PodioInput
...@@ -16,12 +16,17 @@ from Configurables import Jug__Reco__ImagingPixelReco as ImagingPixelReco ...@@ -16,12 +16,17 @@ from Configurables import Jug__Reco__ImagingPixelReco as ImagingPixelReco
from Configurables import Jug__Reco__ImagingTopoCluster as ImagingTopoCluster from Configurables import Jug__Reco__ImagingTopoCluster as ImagingTopoCluster
from Configurables import Jug__Reco__ImagingClusterReco as ImagingClusterReco from Configurables import Jug__Reco__ImagingClusterReco as ImagingClusterReco
# input arguments from calibration file
with open('config/emcal_barrel_calibration.json') as f:
calib_data = json.load(f)['electron']
print(calib_data)
# input arguments through environment variables
kwargs = dict() kwargs = dict()
kwargs['img_sf'] = float(os.environ.get('CB_EMCAL_IMG_SF', '0.01324')) kwargs['img_sf'] = float(calib_data['sampling_fraction_img'])
kwargs['scfi_sf'] = float(os.environ.get('CB_EMCAL_SCFI_SF', '0.0938')) kwargs['scfi_sf'] = float(calib_data['sampling_fraction_scfi'])
# raise error if not defined
# input arguments through environment variables
kwargs['input'] = os.environ['CB_EMCAL_SIM_FILE'] kwargs['input'] = os.environ['CB_EMCAL_SIM_FILE']
kwargs['output'] = os.environ['CB_EMCAL_REC_FILE'] kwargs['output'] = os.environ['CB_EMCAL_REC_FILE']
kwargs['compact'] = os.environ['CB_EMCAL_COMPACT_PATH'] kwargs['compact'] = os.environ['CB_EMCAL_COMPACT_PATH']
...@@ -65,7 +70,7 @@ imcaldigi = CalHitDigi('imcal_digi', ...@@ -65,7 +70,7 @@ imcaldigi = CalHitDigi('imcal_digi',
**imcaldaq) **imcaldaq)
imcalreco = ImagingPixelReco('imcal_reco', imcalreco = ImagingPixelReco('imcal_reco',
# OutputLevel=DEBUG, # OutputLevel=DEBUG,
inputHitCollection='DigiEcalBarrelHits', inputHitCollection=imcaldigi.outputHitCollection,
outputHitCollection='RecoEcalBarrelHits', outputHitCollection='RecoEcalBarrelHits',
readoutClass='EcalBarrelHits', readoutClass='EcalBarrelHits',
layerField='layer', layerField='layer',
...@@ -74,17 +79,19 @@ imcalreco = ImagingPixelReco('imcal_reco', ...@@ -74,17 +79,19 @@ imcalreco = ImagingPixelReco('imcal_reco',
imcalcluster = ImagingTopoCluster('imcal_cluster', imcalcluster = ImagingTopoCluster('imcal_cluster',
# OutputLevel=DEBUG, # OutputLevel=DEBUG,
inputHitCollection='RecoEcalBarrelHits', inputHitCollection=imcalreco.outputHitCollection,
outputHitCollection='EcalBarrelClusterHits', outputProtoClusterCollection='EcalBarrelProtoClusters',
localDistXY=[2.*mm, 2*mm], localDistXY=[2.*mm, 2*mm],
layerDistEtaPhi=[10*mrad, 10*mrad], layerDistEtaPhi=[10*mrad, 10*mrad],
neighbourLayersRange=2, neighbourLayersRange=2,
sectorDist=3.*cm) sectorDist=3.*cm)
clusterreco = ImagingClusterReco('imcal_clreco', clusterreco = ImagingClusterReco('imcal_clreco',
# OutputLevel=DEBUG, # OutputLevel=DEBUG,
inputHitCollection='EcalBarrelClusterHits', inputHitCollection=imcalcluster.inputHitCollection,
inputProtoClusterCollection=imcalcluster.outputProtoClusterCollection,
outputLayerCollection='EcalBarrelClustersLayers', outputLayerCollection='EcalBarrelClustersLayers',
outputClusterCollection='EcalBarrelClusters', outputClusterCollection='EcalBarrelClusters',
outputInfoCollection='EcalBarrelClustersInfo',
samplingFraction=kwargs['img_sf']) samplingFraction=kwargs['img_sf'])
# scfi layers # scfi layers
...@@ -100,7 +107,7 @@ scfi_barrel_digi = CalHitDigi("scfi_barrel_digi", ...@@ -100,7 +107,7 @@ scfi_barrel_digi = CalHitDigi("scfi_barrel_digi",
**scfi_barrel_daq) **scfi_barrel_daq)
scfi_barrel_reco = CalHitReco("scfi_barrel_reco", scfi_barrel_reco = CalHitReco("scfi_barrel_reco",
inputHitCollection="EcalBarrelScFiHitsDigi", inputHitCollection=scfi_barrel_digi.outputHitCollection,
outputHitCollection="EcalBarrelScFiHitsReco", outputHitCollection="EcalBarrelScFiHitsReco",
thresholdFactor=5.0, thresholdFactor=5.0,
readoutClass="EcalBarrelScFiHits", readoutClass="EcalBarrelScFiHits",
...@@ -112,7 +119,7 @@ scfi_barrel_reco = CalHitReco("scfi_barrel_reco", ...@@ -112,7 +119,7 @@ scfi_barrel_reco = CalHitReco("scfi_barrel_reco",
# merge hits in different layer (projection to local x-y plane) # merge hits in different layer (projection to local x-y plane)
scfi_barrel_merger = CalHitsMerger("scfi_barrel_merger", scfi_barrel_merger = CalHitsMerger("scfi_barrel_merger",
# OutputLevel=DEBUG, # OutputLevel=DEBUG,
inputHitCollection="EcalBarrelScFiHitsReco", inputHitCollection=scfi_barrel_reco.outputHitCollection,
outputHitCollection="EcalBarrelScFiGridReco", outputHitCollection="EcalBarrelScFiGridReco",
fields=["fiber"], fields=["fiber"],
fieldRefNumbers=[1], fieldRefNumbers=[1],
...@@ -120,16 +127,18 @@ scfi_barrel_merger = CalHitsMerger("scfi_barrel_merger", ...@@ -120,16 +127,18 @@ scfi_barrel_merger = CalHitsMerger("scfi_barrel_merger",
scfi_barrel_cl = IslandCluster("scfi_barrel_cl", scfi_barrel_cl = IslandCluster("scfi_barrel_cl",
# OutputLevel=DEBUG, # OutputLevel=DEBUG,
inputHitCollection="EcalBarrelScFiGridReco", inputHitCollection=scfi_barrel_merger.outputHitCollection,
outputHitCollection="EcalBarrelScFiClusterHits", outputProtoClusterCollection="EcalBarrelScFiProtoClusters",
splitCluster=False, splitCluster=False,
minClusterCenterEdep=10.*MeV, minClusterCenterEdep=10.*MeV,
localDistXZ=[30*mm, 30*mm], localDistXZ=[30*mm, 30*mm],
sectorDist=5.*cm) sectorDist=5.*cm)
scfi_barrel_clreco = RecoCoG("scfi_barrel_clreco", scfi_barrel_clreco = RecoCoG("scfi_barrel_clreco",
inputHitCollection="EcalBarrelScFiClusterHits", inputHitCollection=scfi_barrel_cl.inputHitCollection,
inputProtoClusterCollection=scfi_barrel_cl.outputProtoClusterCollection,
outputClusterCollection="EcalBarrelScFiClusters", outputClusterCollection="EcalBarrelScFiClusters",
outputInfoCollection="EcalBarrelScFiClustersInfo",
logWeightBase=6.2, logWeightBase=6.2,
samplingFraction=kwargs['scfi_sf']) samplingFraction=kwargs['scfi_sf'])
......
benchmarks/imaging_ecal/options/imaging_2dcluster.py
View file @ 22ed77c9
import json
import os import os
import ROOT import ROOT
from Gaudi.Configuration import * from Gaudi.Configuration import *
...@@ -11,10 +12,14 @@ from Configurables import Jug__Reco__CalorimeterHitsEtaPhiProjector as CalHitsPr ...@@ -11,10 +12,14 @@ from Configurables import Jug__Reco__CalorimeterHitsEtaPhiProjector as CalHitsPr
from Configurables import Jug__Reco__CalorimeterIslandCluster as IslandCluster from Configurables import Jug__Reco__CalorimeterIslandCluster as IslandCluster
from Configurables import Jug__Reco__ClusterRecoCoG as RecoCoG from Configurables import Jug__Reco__ClusterRecoCoG as RecoCoG
# input arguments from calibration file
with open('config/emcal_barrel_calibration.json') as f:
calib_data = json.load(f)['electron']
# input arguments through environment variables
kwargs = dict() kwargs = dict()
kwargs['sf'] = float(os.environ.get('CB_EMCAL_SCFI, SAMP_FRAC', '0.0134')) kwargs['img_sf'] = float(calib_data['sampling_fraction_img'])
# input arguments through environment variables
kwargs['input'] = os.environ.get('CB_EMCAL_SIM_FILE', '../topside/barrel_pion0_5GeV.root') 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['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['compact'] = os.environ.get('CB_EMCAL_COMPACT_PATH', '../topside/test.xml')
...@@ -26,8 +31,6 @@ if kwargs['nev'] < 1: ...@@ -26,8 +31,6 @@ if kwargs['nev'] < 1:
print(kwargs) 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) geo_service = GeoSvc('GeoSvc', detectors=kwargs['compact'].split(','), OutputLevel=INFO)
podioevent = EICDataSvc('EventDataSvc', inputs=kwargs['input'].split(','), OutputLevel=DEBUG) podioevent = EICDataSvc('EventDataSvc', inputs=kwargs['input'].split(','), OutputLevel=DEBUG)
...@@ -64,7 +67,7 @@ imcal_barrel_merger = CalHitsProj('imcal_barrel_merger', ...@@ -64,7 +67,7 @@ imcal_barrel_merger = CalHitsProj('imcal_barrel_merger',
imcal_barrel_cl = IslandCluster('imcal_barrel_cl', imcal_barrel_cl = IslandCluster('imcal_barrel_cl',
# OutputLevel=DEBUG, # OutputLevel=DEBUG,
inputHitCollection=imcal_barrel_merger.outputHitCollection, inputHitCollection=imcal_barrel_merger.outputHitCollection,
outputHitCollection='EcalBarrelClusterHits', outputProtoClusterCollection='EcalBarrelProtoClusters',
minClusterCenterEdep=1.*MeV, minClusterCenterEdep=1.*MeV,
minClusterHitEdep=0.1*MeV, minClusterHitEdep=0.1*MeV,
splitCluster=True, splitCluster=True,
...@@ -72,10 +75,12 @@ imcal_barrel_cl = IslandCluster('imcal_barrel_cl', ...@@ -72,10 +75,12 @@ imcal_barrel_cl = IslandCluster('imcal_barrel_cl',
imcal_barrel_clreco = RecoCoG('imcal_barrel_clreco', imcal_barrel_clreco = RecoCoG('imcal_barrel_clreco',
# OutputLevel=DEBUG, # OutputLevel=DEBUG,
inputHitCollection=imcal_barrel_cl.outputHitCollection, inputHitCollection=imcal_barrel_cl.inputHitCollection,
inputProtoClusterCollection = imcal_barrel_cl.outputProtoClusterCollection,
outputClusterCollection='EcalBarrelClusters', outputClusterCollection='EcalBarrelClusters',
outputInfoCollection='EcalBarrelClustersInfo',
logWeightBase=6.2, logWeightBase=6.2,
samplingFraction=kwargs['sf']) samplingFraction=kwargs['img_sf'])
podout.outputCommands = ['keep *'] podout.outputCommands = ['keep *']
......
benchmarks/imaging_ecal/options/imaging_topocluster.py
View file @ 22ed77c9
import json
import os import os
import ROOT import ROOT
from Gaudi.Configuration import * from Gaudi.Configuration import *
...@@ -14,9 +15,14 @@ from Configurables import Jug__Reco__ImagingTopoCluster as ImagingTopoCluster ...@@ -14,9 +15,14 @@ from Configurables import Jug__Reco__ImagingTopoCluster as ImagingTopoCluster
from Configurables import Jug__Reco__ImagingClusterReco as ImagingClusterReco from Configurables import Jug__Reco__ImagingClusterReco as ImagingClusterReco
# input arguments through environment variables # input arguments from calibration file
with open('config/emcal_barrel_calibration.json') as f:
calib_data = json.load(f)['electron']
kwargs = dict() kwargs = dict()
kwargs['sf'] = float(os.environ.get('CB_EMCAL_SAMP_FRAC', '1.0')) kwargs['sf'] = float(calib_data['sampling_fraction_img'])
# input arguments through environment variables
kwargs['input'] = os.environ.get('CB_EMCAL_SIM_FILE', '../topside/barrel_pion0_5GeV.root') 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['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['compact'] = os.environ.get('CB_EMCAL_COMPACT_PATH', '../topside/test.xml')
...@@ -28,8 +34,6 @@ if kwargs['nev'] < 1: ...@@ -28,8 +34,6 @@ if kwargs['nev'] < 1:
print(kwargs) 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) geo_service = GeoSvc("GeoSvc", detectors=kwargs['compact'].split(','), OutputLevel=INFO)
podioevent = EICDataSvc("EventDataSvc", inputs=kwargs['input'].split(','), OutputLevel=DEBUG) podioevent = EICDataSvc("EventDataSvc", inputs=kwargs['input'].split(','), OutputLevel=DEBUG)
out = PodioOutput("out", filename=kwargs['output']) out = PodioOutput("out", filename=kwargs['output'])
...@@ -59,7 +63,7 @@ imcaldigi = CalorimeterHitDigi("imcal_digi", ...@@ -59,7 +63,7 @@ imcaldigi = CalorimeterHitDigi("imcal_digi",
**daq_setting) **daq_setting)
imcalreco = ImagingPixelReco("imcal_reco", imcalreco = ImagingPixelReco("imcal_reco",
# OutputLevel=DEBUG, # OutputLevel=DEBUG,
inputHitCollection="DigiEcalBarrelHits", inputHitCollection=imcaldigi.outputHitCollection,
outputHitCollection="RecoEcalBarrelHits", outputHitCollection="RecoEcalBarrelHits",
readoutClass="EcalBarrelHits", readoutClass="EcalBarrelHits",
layerField="layer", layerField="layer",
...@@ -68,8 +72,8 @@ imcalreco = ImagingPixelReco("imcal_reco", ...@@ -68,8 +72,8 @@ imcalreco = ImagingPixelReco("imcal_reco",
imcalcluster = ImagingTopoCluster("imcal_cluster", imcalcluster = ImagingTopoCluster("imcal_cluster",
# OutputLevel=DEBUG, # OutputLevel=DEBUG,
inputHitCollection="RecoEcalBarrelHits", inputHitCollection=imcalreco.outputHitCollection,
outputHitCollection="EcalBarrelClusterHits", outputProtoClusterCollection="EcalBarrelProtoClusters",
localDistXY=[2.*units.mm, 2*units.mm], localDistXY=[2.*units.mm, 2*units.mm],
layerDistEtaPhi=[10*units.mrad, 10*units.mrad], layerDistEtaPhi=[10*units.mrad, 10*units.mrad],
neighbourLayersRange=2, neighbourLayersRange=2,
...@@ -77,9 +81,11 @@ imcalcluster = ImagingTopoCluster("imcal_cluster", ...@@ -77,9 +81,11 @@ imcalcluster = ImagingTopoCluster("imcal_cluster",
minClusterNhits=5) minClusterNhits=5)
clusterreco = ImagingClusterReco("imcal_clreco", clusterreco = ImagingClusterReco("imcal_clreco",
# OutputLevel=DEBUG, # OutputLevel=DEBUG,
inputHitCollection="EcalBarrelClusterHits", inputHitCollection=imcalcluster.inputHitCollection,
inputProtoClusterCollection=imcalcluster.outputProtoClusterCollection,
outputLayerCollection="EcalBarrelClustersLayers", outputLayerCollection="EcalBarrelClustersLayers",
outputClusterCollection="EcalBarrelClusters", outputClusterCollection="EcalBarrelClusters",
outputInfoCollection="EcalBarrelClustersInfo",
samplingFraction=sf) samplingFraction=sf)
out.outputCommands = ["keep *"] out.outputCommands = ["keep *"]
......
benchmarks/imaging_ecal/options/scfi_cluster.py
View file @ 22ed77c9
import json
import os import os
import ROOT import ROOT
from Gaudi.Configuration import * from Gaudi.Configuration import *
...@@ -11,10 +12,14 @@ from Configurables import Jug__Reco__CalorimeterHitsMerger as CalHitsMerger ...@@ -11,10 +12,14 @@ from Configurables import Jug__Reco__CalorimeterHitsMerger as CalHitsMerger
from Configurables import Jug__Reco__CalorimeterIslandCluster as IslandCluster from Configurables import Jug__Reco__CalorimeterIslandCluster as IslandCluster
from Configurables import Jug__Reco__ClusterRecoCoG as RecoCoG from Configurables import Jug__Reco__ClusterRecoCoG as RecoCoG
# input arguments from calibration file
with open('config/emcal_barrel_calibration.json') as f:
calib_data = json.load(f)['electron']
# input arguments through environment variables
kwargs = dict() kwargs = dict()
kwargs['sf'] = float(os.environ.get('CB_EMCAL_SCFI, SAMP_FRAC', '0.0938')) kwargs['sf'] = float(calib_data['sampling_fraction_scfi'])
# input arguments through environment variables
kwargs['input'] = os.environ.get('CB_EMCAL_SIM_FILE', '../topside/barrel_pion0_5GeV.root') 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['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['compact'] = os.environ.get('CB_EMCAL_COMPACT_PATH', '../topside/test.xml')
...@@ -47,7 +52,7 @@ scfi_barrel_digi = CalHitDigi("scfi_barrel_digi", ...@@ -47,7 +52,7 @@ scfi_barrel_digi = CalHitDigi("scfi_barrel_digi",
**scfi_barrel_daq) **scfi_barrel_daq)
scfi_barrel_reco = CalHitReco("scfi_barrel_reco", scfi_barrel_reco = CalHitReco("scfi_barrel_reco",
inputHitCollection="EcalBarrelScFiHitsDigi", inputHitCollection=scfi_barrel_digi.outputHitCollection,
outputHitCollection="EcalBarrelScFiHitsReco", outputHitCollection="EcalBarrelScFiHitsReco",
thresholdFactor=5.0, thresholdFactor=5.0,
readoutClass="EcalBarrelScFiHits", readoutClass="EcalBarrelScFiHits",
...@@ -59,7 +64,7 @@ scfi_barrel_reco = CalHitReco("scfi_barrel_reco", ...@@ -59,7 +64,7 @@ scfi_barrel_reco = CalHitReco("scfi_barrel_reco",
# merge hits in different layer (projection to local x-y plane) # merge hits in different layer (projection to local x-y plane)
scfi_barrel_merger = CalHitsMerger("scfi_barrel_merger", scfi_barrel_merger = CalHitsMerger("scfi_barrel_merger",
# OutputLevel=DEBUG, # OutputLevel=DEBUG,
inputHitCollection="EcalBarrelScFiHitsReco", inputHitCollection=scfi_barrel_reco.outputHitCollection,
outputHitCollection="EcalBarrelScFiGridReco", outputHitCollection="EcalBarrelScFiGridReco",
fields=["fiber"], fields=["fiber"],
fieldRefNumbers=[1], fieldRefNumbers=[1],
...@@ -67,15 +72,17 @@ scfi_barrel_merger = CalHitsMerger("scfi_barrel_merger", ...@@ -67,15 +72,17 @@ scfi_barrel_merger = CalHitsMerger("scfi_barrel_merger",
scfi_barrel_cl = IslandCluster("scfi_barrel_cl", scfi_barrel_cl = IslandCluster("scfi_barrel_cl",
# OutputLevel=DEBUG, # OutputLevel=DEBUG,
inputHitCollection="EcalBarrelScFiGridReco", inputHitCollection=scfi_barrel_reco.outputHitCollection,
outputHitCollection="EcalBarrelScFiClusterHits", outputProtoClusterCollection="EcalBarrelScFiProtoClusters",
splitCluster=False, splitCluster=False,
minClusterCenterEdep=10.*MeV, minClusterCenterEdep=10.*MeV,
localDistXZ=[30*mm, 30*mm]) localDistXZ=[30*mm, 30*mm])
scfi_barrel_clreco = RecoCoG("scfi_barrel_clreco", scfi_barrel_clreco = RecoCoG("scfi_barrel_clreco",
inputHitCollection="EcalBarrelScFiClusterHits", inputHitCollection=scfi_barrel_cl.inputHitCollection,
inputProtoClusterCollection=scfi_barrel.outputProtoClusterCollection,
outputClusterCollection="EcalBarrelScFiClusters", outputClusterCollection="EcalBarrelScFiClusters",
outputInfoCollection="EcalBarrelScFiClustersInfo",
logWeightBase=6.2, logWeightBase=6.2,
samplingFraction=kwargs['sf']) samplingFraction=kwargs['sf'])
......
benchmarks/imaging_ecal/scripts/draw_cluster.py
View file @ 22ed77c9
...@@ -26,8 +26,8 @@ import sys ...@@ -26,8 +26,8 @@ import sys
# note z and x axes are switched # note z and x axes are switched
def draw_hits3d(axis, data, cmap, units=('mm', 'mm', 'mm', 'MeV'), fontsize=24, **kwargs): def draw_hits3d(axis, data, cmap, units=('mm', 'mm', 'mm', 'MeV'), fontsize=24, **kwargs):
# normalize energy to get colors # normalize energy to get colors
x, y, z, edep = np.transpose(data) x, y, z, energy = np.transpose(data)
cvals = edep - min(edep) / (max(edep) - min(edep)) cvals = energy - min(energy) / (max(energy) - min(energy))
cvals[np.isnan(cvals)] = 1.0 cvals[np.isnan(cvals)] = 1.0
colors = cmap(cvals) colors = cmap(cvals)
...@@ -37,7 +37,7 @@ def draw_hits3d(axis, data, cmap, units=('mm', 'mm', 'mm', 'MeV'), fontsize=24, ...@@ -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_zlabel('x ({})'.format(units[2]), fontsize=fontsize + 2, labelpad=fontsize)
axis.set_ylabel('y ({})'.format(units[1]), 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) 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) ax=axis, shrink=0.85)
cb.ax.tick_params(labelsize=fontsize) cb.ax.tick_params(labelsize=fontsize)
cb.ax.get_yaxis().labelpad = fontsize cb.ax.get_yaxis().labelpad = fontsize
...@@ -158,7 +158,7 @@ if __name__ == '__main__': ...@@ -158,7 +158,7 @@ if __name__ == '__main__':
# filter out outliers # filter out outliers
dfc = df[(df['eta'] <= df['eta'].quantile(0.95)) & (df['eta'] >= df['eta'].quantile(0.05)) & 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))] (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) vec = vec/np.linalg.norm(vec)
# particle line from (0, 0, 0) to the inner Ecal surface # particle line from (0, 0, 0) to the inner Ecal surface
...@@ -174,7 +174,7 @@ if __name__ == '__main__': ...@@ -174,7 +174,7 @@ if __name__ == '__main__':
# TODO need to implement motion of particles in a field # TODO need to implement motion of particles in a field
# ax.plot(*pline[[2, 1]], '--', zs=pline[0], color='green') # ax.plot(*pline[[2, 1]], '--', zs=pline[0], color='green')
# draw hits # 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, length, rstride=10, cstride=10, color='royalblue')
draw_cylinder3d(ax, rmin + thickness, length, rstride=10, cstride=10, color='forestgreen') draw_cylinder3d(ax, rmin + thickness, length, rstride=10, cstride=10, color='forestgreen')
ax.set_zlim(-(rmin + thickness), rmin + thickness) ax.set_zlim(-(rmin + thickness), rmin + thickness)
...@@ -191,7 +191,7 @@ if __name__ == '__main__': ...@@ -191,7 +191,7 @@ if __name__ == '__main__':
# TODO need to implement motion of particles in a field # TODO need to implement motion of particles in a field
# ax.plot(*pline[[2, 1]], '--', zs=pline[0], color='green') # ax.plot(*pline[[2, 1]], '--', zs=pline[0], color='green')
# draw hits # 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, # 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)) # scat_kw=dict(color='k', s=50.0), line_kw=dict(linestyle=':', color='k', lw=3))
# view range # view range
...@@ -213,7 +213,7 @@ if __name__ == '__main__': ...@@ -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. 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]}) 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)), 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')) cmap=cmap, cmin=0., pc_kw=dict(alpha=0.8, edgecolor='k'))
......