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 (11)
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 282 additions and 736 deletions
.gitlab-ci.yml
View file @ 0272a2a8
......@@ -3,6 +3,7 @@ image: eicweb.phy.anl.gov:4567/eic/juggler/juggler:$JUGGLER_TAG
default:
before_script:
- source .local/bin/env.sh
interruptible: true
artifacts:
expire_in: 3 days
paths:
......@@ -58,11 +59,12 @@ common:detector:
- ln -s "${LOCAL_DATA_PATH}/sim_output" sim_output
- ln -s "${LOCAL_DATA_PATH}/datasets/data" data
- ls -lrtha
- curl -sL https://eicweb.phy.anl.gov/EIC/benchmarks/detector_benchmarks/-/jobs/artifacts/master/raw/results/emcal_barrel_calibration.json?job=collect_results:barrel_ecal --output config/emcal_barrel_calibration.json
- bash bin/get_calibrations
interruptible: true
include:
- local: 'benchmarks/ecal/config.yml'
- local: 'benchmarks/track_finding/config.yml'
- local: 'benchmarks/tracking/config.yml'
- local: 'benchmarks/clustering/config.yml'
# - local: 'benchmarks/rich/config.yml'
......
benchmarks/B0/scripts/gen_far_forward_protons.cxx
View file @ 0272a2a8
......@@ -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";
// 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>(
FourVector(
px, py, pz,
sqrt(p*p + (0.000511 * 0.000511))),
11, 1);
sqrt(p*p + (0.938272 * 0.938272))),
2212, 1);
GenVertexPtr v1 = std::make_shared<GenVertex>();
v1->add_particle_in(p1);
......
benchmarks/B0/scripts/hits_far_forward_protons.cxx
View file @ 0272a2a8
......@@ -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> result;
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;
}
......@@ -54,7 +54,7 @@ auto fourvec = [](ROOT::VecOps::RVec<dd4pod::Geant4ParticleData> const& in) {
std::vector<ROOT::Math::PxPyPzMVector> result;
ROOT::Math::PxPyPzMVector lv;
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);
}
return result;
......
benchmarks/B0/scripts/rec_far_forward_protons.cxx
View file @ 0272a2a8
......@@ -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> result;
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;
}
......@@ -52,7 +52,7 @@ auto fourvec = [](ROOT::VecOps::RVec<dd4pod::Geant4ParticleData> const& in) {
std::vector<ROOT::Math::PxPyPzMVector> result;
ROOT::Math::PxPyPzMVector lv;
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);
}
return result;
......
benchmarks/clustering/full_cal_clusters.sh
View file @ 0272a2a8
......@@ -130,8 +130,8 @@ fi
# Run analysis scripts
FULL_CAL_SCRIPT_DIR=benchmarks/clustering/scripts
python ${FULL_CAL_SCRIPT_DIR}/cluster_plots.py ${JUGGLER_SIM_FILE} ${JUGGLER_REC_FILE} -o results \
--collections "EcalEndcapNClusters, EcalEndcapPClusters, EcalBarrelClusters,
HcalElectronEndcapClusters, HcalHadronEndcapClusters, HcalBarrelClusters"
--collections "EcalEndcapNClusters, EcalEndcapPClusters, EcalBarrelImagingClusters,
HcalEndcapNClusters, HcalEndcapPClusters, HcalBarrelClusters"
root_filesize=$(stat --format=%s "${JUGGLER_REC_FILE}")
if [[ "${JUGGLER_N_EVENTS}" -lt "500" ]] ; then
......
benchmarks/clustering/options/deprecated/full_trackpluscalo_reco.py deleted 100644 → 0
View file @ e0a3cd35
'''
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/options/full_cal_reco.py
View file @ 0272a2a8
......@@ -60,8 +60,8 @@ sim_coll = [
"EcalBarrelHits",
"HcalBarrelHits",
"EcalBarrelScFiHits",
"HcalHadronEndcapHits",
"HcalElectronEndcapHits",
"HcalEndcapPHits",
"HcalEndcapNHits",
]
# input and output
......@@ -162,13 +162,13 @@ cb_ecal_daq = dict(
cb_ecal_digi = CalHitDigi("cb_ecal_digi",
inputHitCollection="EcalBarrelHits",
outputHitCollection="EcalBarrelHitsDigi",
outputHitCollection="EcalBarrelImagingHitsDigi",
energyResolutions=[0., 0.02, 0.], # 2% flat resolution
**cb_ecal_daq)
cb_ecal_reco = ImCalPixelReco("cb_ecal_reco",
inputHitCollection=cb_ecal_digi.outputHitCollection,
outputHitCollection="EcalBarrelHitsReco",
outputHitCollection="EcalBarrelImagingHitsReco",
thresholdFactor=3, # about 20 keV
readoutClass="EcalBarrelHits", # readout class
layerField="layer", # field to get layer id
......@@ -177,7 +177,7 @@ cb_ecal_reco = ImCalPixelReco("cb_ecal_reco",
cb_ecal_cl = ImagingCluster("cb_ecal_cl",
inputHitCollection=cb_ecal_reco.outputHitCollection,
outputProtoClusterCollection="EcalBarrelProtoClusters",
outputProtoClusterCollection="EcalBarrelImagingProtoClusters",
localDistXY=[2.*mm, 2*mm], # same layer
layerDistEtaPhi=[10*mrad, 10*mrad], # adjacent layer
neighbourLayersRange=2, # id diff for adjacent layer
......@@ -187,9 +187,9 @@ cb_ecal_clreco = ImagingClusterReco("cb_ecal_clreco",
samplingFraction=cb_ecal_sf,
inputHitCollection=cb_ecal_cl.inputHitCollection,
inputProtoClusterCollection=cb_ecal_cl.outputProtoClusterCollection,
outputClusterCollection="EcalBarrelClusters",
outputInfoCollection="EcalBarrelClustersInfo",
outputLayerCollection="EcalBarrelLayers")
outputClusterCollection="EcalBarrelImagingClusters",
outputInfoCollection="EcalBarrelImagingClustersInfo",
outputLayerCollection="EcalBarrelImagingLayers")
#Central ECAL SciFi
# use the same daq_setting for digi/reco pair
......@@ -291,26 +291,26 @@ ci_hcal_daq = dict(
pedestalMean=400,
pedestalSigma=10)
ci_hcal_digi = CalHitDigi("ci_hcal_digi",
inputHitCollection="HcalHadronEndcapHits",
outputHitCollection="HcalHadronEndcapHitsDigi",
inputHitCollection="HcalEndcapPHits",
outputHitCollection="HcalEndcapPHitsDigi",
**ci_hcal_daq)
ci_hcal_reco = CalHitReco("ci_hcal_reco",
inputHitCollection=ci_hcal_digi.outputHitCollection,
outputHitCollection="HcalHadronEndcapHitsReco",
outputHitCollection="HcalEndcapPHitsReco",
thresholdFactor=5.0,
**ci_hcal_daq)
ci_hcal_merger = CalHitsMerger("ci_hcal_merger",
inputHitCollection=ci_hcal_reco.outputHitCollection,
outputHitCollection="HcalHadronEndcapHitsRecoXY",
readoutClass="HcalHadronEndcapHits",
outputHitCollection="HcalEndcapPHitsRecoXY",
readoutClass="HcalEndcapPHits",
fields=["layer", "slice"],
fieldRefNumbers=[1, 0])
ci_hcal_cl = IslandCluster("ci_hcal_cl",
inputHitCollection=ci_hcal_merger.outputHitCollection,
outputProtoClusterCollection="HcalHadronEndcapProtoClusters",
outputProtoClusterCollection="HcalEndcapPProtoClusters",
splitCluster=False,
minClusterCenterEdep=30.*MeV,
localDistXY=[15.*cm, 15.*cm])
......@@ -318,8 +318,8 @@ ci_hcal_cl = IslandCluster("ci_hcal_cl",
ci_hcal_clreco = RecoCoG("ci_hcal_clreco",
inputHitCollection=ci_hcal_cl.inputHitCollection,
inputProtoClusterCollection=ci_hcal_cl.outputProtoClusterCollection,
outputClusterCollection="HcalHadronEndcapClusters",
outputInfoCollection="HcalHadronEndcapClustersInfo",
outputClusterCollection="HcalEndcapPClusters",
outputInfoCollection="HcalEndcapPClustersInfo",
logWeightBase=6.2,
samplingFraction=ci_hcal_sf)
......@@ -331,26 +331,26 @@ ce_hcal_daq = dict(
pedestalSigma=10)
ce_hcal_digi = CalHitDigi("ce_hcal_digi",
inputHitCollection="HcalElectronEndcapHits",
outputHitCollection="HcalElectronEndcapHitsDigi",
inputHitCollection="HcalEndcapNHits",
outputHitCollection="HcalEndcapNHitsDigi",
**ce_hcal_daq)
ce_hcal_reco = CalHitReco("ce_hcal_reco",
inputHitCollection=ce_hcal_digi.outputHitCollection,
outputHitCollection="HcalElectronEndcapHitsReco",
outputHitCollection="HcalEndcapNHitsReco",
thresholdFactor=5.0,
**ce_hcal_daq)
ce_hcal_merger = CalHitsMerger("ce_hcal_merger",
inputHitCollection=ce_hcal_reco.outputHitCollection,
outputHitCollection="HcalElectronEndcapHitsRecoXY",
readoutClass="HcalElectronEndcapHits",
outputHitCollection="HcalEndcapNHitsRecoXY",
readoutClass="HcalEndcapNHits",
fields=["layer", "slice"],
fieldRefNumbers=[1, 0])
ce_hcal_cl = IslandCluster("ce_hcal_cl",
inputHitCollection=ce_hcal_merger.outputHitCollection,
outputProtoClusterCollection="HcalElectronEndcapProtoClusters",
outputProtoClusterCollection="HcalEndcapNProtoClusters",
splitCluster=False,
minClusterCenterEdep=30.*MeV,
localDistXY=[15.*cm, 15.*cm])
......@@ -358,8 +358,8 @@ ce_hcal_cl = IslandCluster("ce_hcal_cl",
ce_hcal_clreco = RecoCoG("ce_hcal_clreco",
inputHitCollection=ce_hcal_cl.inputHitCollection,
inputProtoClusterCollection=ce_hcal_cl.outputProtoClusterCollection,
outputClusterCollection="HcalElectronEndcapClusters",
outputInfoCollection="HcalElectronEndcapClustersInfo",
outputClusterCollection="HcalEndcapNClusters",
outputInfoCollection="HcalEndcapNClustersInfo",
logWeightBase=6.2,
samplingFraction=ce_hcal_sf)
......
benchmarks/clustering/scripts/cluster_plots.py
View file @ 0272a2a8
......@@ -51,7 +51,7 @@ def flatten_collection(rdf, collection, cols=None):
def thrown_particles_figure(rdf, save, mcbranch="mcparticles"):
# 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 + 'Info.', '') for c in dft.columns}, inplace=True)
# select thrown particles
......@@ -79,7 +79,7 @@ def thrown_particles_figure(rdf, save, mcbranch="mcparticles"):
# calculate kinematics
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[:, 'eta'] = [v.Eta() for v in fourvecs]
......
benchmarks/ecal/options/barrel.py
View file @ 0272a2a8
......@@ -59,13 +59,13 @@ cb_ecal_daq = dict(
cb_ecal_digi = CalHitDigi("cb_ecal_digi",
inputHitCollection="EcalBarrelHits",
outputHitCollection="EcalBarrelHitsDigi",
outputHitCollection="EcalBarrelImagingHitsDigi",
energyResolutions=[0., 0.02, 0.], # 2% flat resolution
**cb_ecal_daq)
cb_ecal_reco = ImCalPixelReco("cb_ecal_reco",
inputHitCollection=cb_ecal_digi.outputHitCollection,
outputHitCollection="EcalBarrelHitsReco",
outputHitCollection="EcalBarrelImagingHitsReco",
thresholdFactor=3, # about 20 keV
readoutClass="EcalBarrelHits", # readout class
layerField="layer", # field to get layer id
......@@ -74,7 +74,7 @@ cb_ecal_reco = ImCalPixelReco("cb_ecal_reco",
cb_ecal_cl = ImagingCluster("cb_ecal_cl",
inputHitCollection=cb_ecal_reco.outputHitCollection,
outputProtoClusterCollection="EcalBarrelProtoClusters",
outputProtoClusterCollection="EcalBarrelImagingProtoClusters",
localDistXY=[2.*mm, 2*mm], # same layer
layerDistEtaPhi=[10*mrad, 10*mrad], # adjacent layer
neighbourLayersRange=2, # id diff for adjacent layer
......@@ -84,9 +84,9 @@ cb_ecal_clreco = ImagingClusterReco("cb_ecal_clreco",
samplingFraction=cb_ecal_sf,
inputHitCollection=cb_ecal_cl.inputHitCollection,
inputProtoClusterCollection=cb_ecal_cl.outputProtoClusterCollection,
outputClusterCollection="EcalBarrelClusters",
outputLayerCollection="EcalBarrelLayers",
outputInfoCollection="EcalBarrelClustersInfo")
outputClusterCollection="EcalBarrelImagingClusters",
outputLayerCollection="EcalBarrelImagingLayers",
outputInfoCollection="EcalBarrelImagingClustersInfo")
podout.outputCommands = ['drop *',
'keep mcparticles2',
......
benchmarks/ecal/run_emcal_benchmarks.py
View file @ 0272a2a8
......@@ -24,11 +24,11 @@ default_type = {
'barrel': [
['barrel.py'],
['draw_clusters.py'],
['EcalBarrelClusters']],
['EcalBarrelImagingClusters']],
# 'all': [
# ['all_ecal.py'],
# ['draw_clusters.py'],
# ['EcalEndcapNClusters', 'EcalEndcapPClusters', 'EcalBarrelClusters']],
# ['EcalEndcapNClusters', 'EcalEndcapPClusters', 'EcalBarrelImagingClusters']],
}
default_compact = os.path.join(os.environ.get('DETECTOR_PATH', os.environ.get('DETECTOR_PATH', '')),
......
benchmarks/full/full_reconstruction.sh
View file @ 0272a2a8
......@@ -164,13 +164,14 @@ OPTION_DIR=benchmarks/full/options
SCRIPT_DIR=benchmarks/full/scripts
# Run Juggler
xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv \
gaudirun.py ${OPTION_DIR}/full_reconstruction.py
gaudirun.py ${OPTION_DIR}/full_reconstruction.py
if [[ "$?" -ne "0" ]] ; then
echo "ERROR running reconstruction (juggler)"
exit 1
fi
rootls -t "${JUGGLER_REC_FILE}"
# Run analysis scripts
#python ${SCRIPT_DIR}/full_reconstruction.py ${JUGGLER_SIM_FILE} ${JUGGLER_REC_FILE} -o results
......
benchmarks/imaging_ecal/options/hybrid_cluster.py
View file @ 0272a2a8
......@@ -65,13 +65,13 @@ imcaldaq = dict(
imcaldigi = CalHitDigi('imcal_digi',
# OutputLevel=DEBUG,
inputHitCollection='EcalBarrelHits',
outputHitCollection='DigiEcalBarrelHits',
outputHitCollection='DigiEcalBarrelImagingHits',
energyResolutions=[0., 0.02, 0.],
**imcaldaq)
imcalreco = ImagingPixelReco('imcal_reco',
# OutputLevel=DEBUG,
inputHitCollection=imcaldigi.outputHitCollection,
outputHitCollection='RecoEcalBarrelHits',
outputHitCollection='RecoEcalBarrelImagingHits',
readoutClass='EcalBarrelHits',
layerField='layer',
sectorField='module',
......@@ -80,7 +80,7 @@ imcalreco = ImagingPixelReco('imcal_reco',
imcalcluster = ImagingTopoCluster('imcal_cluster',
# OutputLevel=DEBUG,
inputHitCollection=imcalreco.outputHitCollection,
outputProtoClusterCollection='EcalBarrelProtoClusters',
outputProtoClusterCollection='EcalBarrelImagingProtoClusters',
localDistXY=[2.*mm, 2*mm],
layerDistEtaPhi=[10*mrad, 10*mrad],
neighbourLayersRange=2,
......@@ -89,9 +89,9 @@ clusterreco = ImagingClusterReco('imcal_clreco',
# OutputLevel=DEBUG,
inputHitCollection=imcalcluster.inputHitCollection,
inputProtoClusterCollection=imcalcluster.outputProtoClusterCollection,
outputLayerCollection='EcalBarrelClustersLayers',
outputClusterCollection='EcalBarrelClusters',
outputInfoCollection='EcalBarrelClustersInfo',
outputLayerCollection='EcalBarrelImagingClustersLayers',
outputClusterCollection='EcalBarrelImagingClusters',
outputInfoCollection='EcalBarrelImagingClustersInfo',
samplingFraction=kwargs['img_sf'])
# scfi layers
......
benchmarks/imaging_ecal/options/imaging_2dcluster.py
View file @ 0272a2a8
......@@ -6,6 +6,7 @@ from GaudiKernel.SystemOfUnits import mm, MeV, rad, ns
from Configurables import ApplicationMgr, EICDataSvc, PodioOutput, GeoSvc
from Configurables import PodioInput
from Configurables import Jug__Base__InputCopier_dd4pod__Geant4ParticleCollection_dd4pod__Geant4ParticleCollection_ as MCCopier
from Configurables import Jug__Digi__CalorimeterHitDigi as CalHitDigi
from Configurables import Jug__Reco__CalorimeterHitReco as CalHitReco
from Configurables import Jug__Reco__CalorimeterHitsEtaPhiProjector as CalHitsProj
......@@ -17,7 +18,7 @@ with open('config/emcal_barrel_calibration.json') as f:
calib_data = json.load(f)['electron']
kwargs = dict()
kwargs['sf'] = float(calib_data['sampling_fraction_scfi'])
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')
......@@ -37,6 +38,12 @@ podioevent = EICDataSvc('EventDataSvc', inputs=kwargs['input'].split(','), Outpu
podin = PodioInput('PodioReader', collections=['mcparticles', 'EcalBarrelHits'], OutputLevel=DEBUG)
podout = PodioOutput('out', filename=kwargs['output'])
mccopier = MCCopier('MCCopier',
# OutputLevel=DEBUG,
inputCollection='mcparticles',
outputCollection='mcparticles2')
# use the same daq_setting for digi/reco pair
imcal_barrel_daq = dict(
dynamicRangeADC=3.*MeV,
......@@ -46,12 +53,12 @@ imcal_barrel_daq = dict(
imcal_barrel_digi = CalHitDigi('imcal_barrel_digi',
inputHitCollection='EcalBarrelHits',
outputHitCollection='EcalBarrelHitsDigi',
outputHitCollection='EcalBarrelImagingHitsDigi',
**imcal_barrel_daq)
imcal_barrel_reco = CalHitReco('imcal_barrel_reco',
inputHitCollection=imcal_barrel_digi.outputHitCollection,
outputHitCollection='EcalBarrelHitsReco',
outputHitCollection='EcalBarrelImagingHitsReco',
thresholdFactor=5.0,
readoutClass='EcalBarrelHits',
localDetFields=['system', 'module'], # use local coordinates in each module (stave)
......@@ -61,13 +68,13 @@ imcal_barrel_reco = CalHitReco('imcal_barrel_reco',
imcal_barrel_merger = CalHitsProj('imcal_barrel_merger',
# OutputLevel=DEBUG,
inputHitCollection=imcal_barrel_reco.outputHitCollection,
outputHitCollection='EcalBarrelHitsMerge',
outputHitCollection='EcalBarrelImagingHitsMerge',
gridSizes=[0.004, 0.004*rad])
imcal_barrel_cl = IslandCluster('imcal_barrel_cl',
# OutputLevel=DEBUG,
inputHitCollection=imcal_barrel_merger.outputHitCollection,
outputProtoClusterCollection='EcalBarrelProtoClusters',
outputProtoClusterCollection='EcalBarrelImagingProtoClusters',
minClusterCenterEdep=1.*MeV,
minClusterHitEdep=0.1*MeV,
splitCluster=True,
......@@ -77,16 +84,16 @@ imcal_barrel_clreco = RecoCoG('imcal_barrel_clreco',
# OutputLevel=DEBUG,
inputHitCollection=imcal_barrel_cl.inputHitCollection,
inputProtoClusterCollection = imcal_barrel_cl.outputProtoClusterCollection,
outputClusterCollection='EcalBarrelClusters',
outputInfoCollection='EcalBarrelClustersInfo',
outputClusterCollection='EcalBarrelImagingClusters',
outputInfoCollection='EcalBarrelImagingClustersInfo',
logWeightBase=6.2,
samplingFraction=kwargs['sf'])
samplingFraction=kwargs['img_sf'])
podout.outputCommands = ['keep *']
ApplicationMgr(
TopAlg=[podin,
TopAlg=[podin, mccopier,
imcal_barrel_digi, imcal_barrel_reco,
imcal_barrel_merger, imcal_barrel_cl, imcal_barrel_clreco,
podout],
......
benchmarks/imaging_ecal/options/imaging_topocluster.py
View file @ 0272a2a8
......@@ -58,13 +58,13 @@ daq_setting = dict(
imcaldigi = CalorimeterHitDigi("imcal_digi",
# OutputLevel=DEBUG,
inputHitCollection="EcalBarrelHits",
outputHitCollection="DigiEcalBarrelHits",
outputHitCollection="DigiEcalBarrelImagingHits",
energyResolutions=[0., 0.02, 0.],
**daq_setting)
imcalreco = ImagingPixelReco("imcal_reco",
# OutputLevel=DEBUG,
inputHitCollection=imcaldigi.outputHitCollection,
outputHitCollection="RecoEcalBarrelHits",
outputHitCollection="RecoEcalBarrelImagingHits",
readoutClass="EcalBarrelHits",
layerField="layer",
sectorField="module",
......@@ -73,7 +73,7 @@ imcalreco = ImagingPixelReco("imcal_reco",
imcalcluster = ImagingTopoCluster("imcal_cluster",
# OutputLevel=DEBUG,
inputHitCollection=imcalreco.outputHitCollection,
outputProtoClusterCollection="EcalBarrelProtoClusters",
outputProtoClusterCollection="EcalBarrelImagingProtoClusters",
localDistXY=[2.*units.mm, 2*units.mm],
layerDistEtaPhi=[10*units.mrad, 10*units.mrad],
neighbourLayersRange=2,
......@@ -83,9 +83,9 @@ clusterreco = ImagingClusterReco("imcal_clreco",
# OutputLevel=DEBUG,
inputHitCollection=imcalcluster.inputHitCollection,
inputProtoClusterCollection=imcalcluster.outputProtoClusterCollection,
outputLayerCollection="EcalBarrelClustersLayers",
outputClusterCollection="EcalBarrelClusters",
outputInfoCollection="EcalBarrelClustersInfo",
outputLayerCollection="EcalBarrelImagingClustersLayers",
outputClusterCollection="EcalBarrelImagingClusters",
outputInfoCollection="EcalBarrelImagingClustersInfo",
samplingFraction=sf)
out.outputCommands = ["keep *"]
......
benchmarks/imaging_ecal/run_emcal_barrel.sh
View file @ 0272a2a8
......@@ -87,7 +87,7 @@ ls -lh ${CB_EMCAL_GEN_FILE}
# Run geant4 simulations
npsim --runType batch \
-v WARNING \
--part.minimalKineticEnergy "1*TeV" \
--part.minimalKineticEnergy "0.5*MeV" \
--numberOfEvents ${CB_EMCAL_NUMEV} \
--compactFile ${CB_EMCAL_COMPACT_PATH} \
--inputFiles ${CB_EMCAL_GEN_FILE} \
......@@ -116,7 +116,7 @@ fi
# Plot clusters first
FULL_CAL_SCRIPT_DIR=benchmarks/clustering/scripts
python ${FULL_CAL_SCRIPT_DIR}/cluster_plots.py ${CB_EMCAL_SIM_FILE} ${CB_EMCAL_REC_FILE} -o results \
--collections "EcalBarrelClusters, EcalBarrelScFiClusters"
--collections "EcalBarrelImagingClusters, EcalBarrelScFiClusters"
# check required python modules
python -m pip install -r benchmarks/imaging_ecal/requirements.txt
......
benchmarks/imaging_ecal/run_imcal_pion0.sh
View file @ 0272a2a8
......@@ -83,7 +83,7 @@ ls -lh ${CB_EMCAL_GEN_FILE}
# Run geant4 simulations
npsim --runType batch \
-v WARNING \
--part.minimalKineticEnergy "1*TeV" \
--part.minimalKineticEnergy "0.5*MeV" \
--numberOfEvents ${CB_EMCAL_NUMEV} \
--compactFile ${CB_EMCAL_COMPACT_PATH} \
--inputFiles ${CB_EMCAL_GEN_FILE} \
......@@ -111,7 +111,7 @@ fi
# Plot clusters first
FULL_CAL_SCRIPT_DIR=benchmarks/clustering/scripts
python ${FULL_CAL_SCRIPT_DIR}/cluster_plots.py ${CB_EMCAL_SIM_FILE} ${CB_EMCAL_REC_FILE} -o results \
--collections "EcalBarrelClusters"
--collections "EcalBarrelImagingClusters"
root_filesize=$(stat --format=%s "${CB_EMCAL_REC_FILE}")
if [[ "${CB_EMCAL_NUMEV}" -lt "500" ]] ; then
......
benchmarks/imaging_ecal/scripts/draw_cluster.py
View file @ 0272a2a8
......@@ -5,6 +5,10 @@
Author: Chao Peng (ANL)
Date: 04/30/2021
Added true decaying particles on eta-phi plane projection plot
Author: Jihee Kim (ANL)
Data: 08/06/2021
'''
import os
......@@ -101,7 +105,7 @@ if __name__ == '__main__':
parser.add_argument('--compact', type=str, default='', dest='compact', help='compact file')
parser.add_argument('-m', '--macros', type=str, default='rootlogon.C', dest='macros',
help='root macros to load (accept multiple paths separated by \",\")')
parser.add_argument('-b', '--branch-name', type=str, default='RecoEcalBarrelHits', dest='branch',
parser.add_argument('-b', '--branch-name', type=str, default='RecoEcalBarrelImagingHits', dest='branch',
help='branch name in the root file (outputLayerCollection from ImagingClusterReco)')
parser.add_argument('--topo-size', type=float, default=2.0, dest='topo_size',
help='bin size for projection plot (mrad)')
......@@ -143,9 +147,26 @@ if __name__ == '__main__':
df['theta'] = np.arccos(df['z'].values/df['r'].values)*1000.
df['eta'] = -np.log(np.tan(df['theta'].values/1000./2.))
# Read all mc particles
dfallmcp = get_all_mcp(args.file, args.iev, 'mcparticles2')
pdgbase = ROOT.TDatabasePDG()
# Select decaying particles
dftemp = dfallmcp[dfallmcp['g4Parent'] == 1.0]
if len(dftemp) > 0:
dfdecaymcp = dftemp.copy()
for iptl in [0, len(dfdecaymcp) - 1]:
infoptl = pdgbase.GetParticle(int(dfdecaymcp['pid'].iloc[iptl]))
print("{} Decaying particle = {}, pdgcode = {}, charge = {}, mass = {}"\
.format(iptl, infoptl.GetName(), infoptl.PdgCode(), infoptl.Charge(), infoptl.Mass()))
# Calculate geometric variables of decaying particles
dfdecaymcp['r'] = np.sqrt(dfdecaymcp['vex'].values**2 + dfdecaymcp['vey'].values**2 + dfdecaymcp['vez'].values**2)
dfdecaymcp['phi'] = np.arctan2(dfdecaymcp['vey'].values, dfdecaymcp['vex'].values)*1000.
dfdecaymcp['theta'] = np.arccos(dfdecaymcp['vez'].values/dfdecaymcp['r'].values)*1000.
dfdecaymcp['eta'] = -np.log(np.tan(dfdecaymcp['theta'].values/1000./2.))
# truth
dfmcp = get_mcp_simple(args.file, args.iev, 'mcparticles2').iloc[0]
pdgbase = ROOT.TDatabasePDG()
#pdgbase = ROOT.TDatabasePDG()
inpart = pdgbase.GetParticle(int(dfmcp['pid']))
print("Incoming particle = {}, pdgcode = {}, charge = {}, mass = {}"\
.format(inpart.GetName(), inpart.PdgCode(), inpart.Charge(), inpart.Mass()))
......@@ -216,6 +237,9 @@ if __name__ == '__main__':
ax, sm = draw_heatmap(axs[0], df['eta'].values, df['phi'].values, weights=df['energy'].values,
bins=(np.arange(*eta_rg, step=args.topo_size/1000.), np.arange(*phi_rg, step=args.topo_size)),
cmap=cmap, cmin=0., pc_kw=dict(alpha=0.8, edgecolor='k'))
# draw true decaying particle position
if len(dftemp) > 0:
ax.scatter(dfdecaymcp['eta'].values, dfdecaymcp['phi'].values, marker='x', color='red', s=22**2, linewidth=5.0)
ax.set_ylabel(r'$\phi$ (mrad)', fontsize=32)
ax.set_xlabel(r'$\eta$', fontsize=32)
......
benchmarks/imaging_ecal/scripts/draw_cluster_layers.py
View file @ 0272a2a8
......@@ -57,7 +57,7 @@ if __name__ == '__main__':
parser.add_argument('--compact', type=str, default='', dest='compact', help='compact file')
parser.add_argument('-m', '--macros', type=str, default='rootlogon.C', dest='macros',
help='root macros to load (accept multiple paths separated by \",\")')
parser.add_argument('-b', '--branch-name', type=str, default='RecoEcalBarrelHits', dest='branch',
parser.add_argument('-b', '--branch-name', type=str, default='RecoEcalBarrelImagingHits', dest='branch',
help='branch name in the root file (outputLayerCollection from ImagingClusterReco)')
parser.add_argument('--topo-size', type=float, default=2.0, dest='topo_size',
help='bin size for projection plot (mrad)')
......
benchmarks/imaging_ecal/scripts/energy_profile.py
View file @ 0272a2a8
......@@ -35,7 +35,7 @@ if __name__ == '__main__':
parser.add_argument('-e', '--energy', type=float, default=5., dest='energy', help='incident particle energy (GeV)')
parser.add_argument('-s', '--save', type=str, default='', dest='save', help='path to save profile')
parser.add_argument('-c', '--color', type=str, default='royalblue', dest='color', help='colors for bar plots')
parser.add_argument('-b', '--branch-name', type=str, default='EcalBarrelClustersLayers', dest='branch',
parser.add_argument('-b', '--branch-name', type=str, default='EcalBarrelImagingClustersLayers', dest='branch',
help='branch name in the root file (outputLayerCollection from ImagingClusterReco)')
parser.add_argument('-m', '--macros', type=str, default='rootlogon.C', dest='macros',
help='root macros to load (accept multiple paths separated by \",\")')
......