File indexing completed on 2026-09-22 08:06:01
0001
0002
0003
0004
0005
0006
0007 #include "lfhcal_studiesProcessor.h"
0008
0009 #include <DD4hep/Detector.h>
0010 #include <DD4hep/IDDescriptor.h>
0011 #include <DD4hep/Readout.h>
0012 #include <JANA/JApplication.h>
0013 #include <JANA/JApplicationFwd.h>
0014 #include <JANA/JEvent.h>
0015 #include <JANA/Services/JGlobalRootLock.h>
0016 #include <RtypesCore.h>
0017 #include <TMath.h>
0018 #include <edm4eic/CalorimeterHitCollection.h>
0019 #include <edm4eic/ClusterCollection.h>
0020 #include <edm4hep/CaloHitContributionCollection.h>
0021 #include <edm4hep/MCParticleCollection.h>
0022 #include <edm4hep/SimCalorimeterHitCollection.h>
0023 #include <edm4hep/Vector3d.h>
0024 #include <edm4hep/Vector3f.h>
0025 #include <fmt/format.h>
0026 #include <podio/RelationRange.h>
0027 #include <algorithm>
0028 #include <cmath>
0029 #include <cstdint>
0030 #include <gsl/pointers>
0031 #include <iostream>
0032 #include <limits>
0033 #include <stdexcept>
0034 #include <vector>
0035
0036 #include "clusterizer_MA.h"
0037 #include "services/geometry/dd4hep/DD4hep_service.h"
0038 #include "services/log/Log_service.h"
0039 #include "services/rootfile/RootFile_service.h"
0040
0041
0042
0043
0044 void lfhcal_studiesProcessor::Init() {
0045 std::string plugin_name = ("lfhcal_studies");
0046
0047
0048
0049
0050 auto* app = GetApplication();
0051
0052 m_log = app->GetService<Log_service>()->logger(plugin_name);
0053
0054
0055 auto dd4hep_service = app->GetService<DD4hep_service>();
0056
0057
0058 auto root_file_service = app->GetService<RootFile_service>();
0059
0060
0061 auto globalRootLock = app->GetService<JGlobalRootLock>();
0062 globalRootLock->acquire_write_lock();
0063 auto* file = root_file_service->GetHistFile();
0064 globalRootLock->release_lock();
0065
0066
0067
0068
0069 m_dir_main = file->mkdir(plugin_name.c_str());
0070
0071
0072
0073
0074 hMCEnergyVsEta = new TH2D("hMCEnergyVsEta", "; E (GeV); #eta", 1500, 0., 150., 500, 0, 5);
0075 hMCEnergyVsEta->SetDirectory(m_dir_main);
0076
0077
0078
0079
0080 hClusterEcalib_E_eta =
0081 new TH3D("hClusterEcalib_E_eta", "; E_{MC} (GeV); E_{rec,rec hit}/E_{MC}; #eta", 1500, 0.,
0082 150.0, 200, 0., 2.0, 50, 0, 5);
0083 hClusterNCells_E_eta = new TH3D("hClusterNCells_E_eta", "; E_{MC} (GeV); N_{cells}; #eta", 1500,
0084 0., 150.0, 500, -0.5, 499.5, 50, 0, 5);
0085 hClusterEcalib_E_phi =
0086 new TH3D("hClusterEcalib_E_phi", "; E_{MC} (GeV); E_{rec,rec hit}/E_{MC}; #varphi (rad)",
0087 1500, 0., 150.0, 200, 0., 2.0, 360, -TMath::Pi(), TMath::Pi());
0088 hPosCaloHitsXY =
0089 new TH2D("hPosCaloHitsXY", "; X (cm); Y (cm)", 400, -400., 400., 400, -400., 400.);
0090 hPosCaloHitsZX =
0091 new TH2D("hPosCaloHitsZX", "; Z (cm); X (cm)", 200, 300., 500., 400, -400., 400.);
0092 hPosCaloHitsZY =
0093 new TH2D("hPosCaloHitsZY", "; Z (cm); Y (cm)", 200, 300., 500., 400, -400., 400.);
0094 hClusterEcalib_E_eta->SetDirectory(m_dir_main);
0095 hClusterNCells_E_eta->SetDirectory(m_dir_main);
0096 hClusterEcalib_E_phi->SetDirectory(m_dir_main);
0097 hPosCaloHitsXY->SetDirectory(m_dir_main);
0098 hPosCaloHitsZX->SetDirectory(m_dir_main);
0099 hPosCaloHitsZY->SetDirectory(m_dir_main);
0100
0101
0102
0103
0104 hClusterESimcalib_E_eta =
0105 new TH3D("hClusterESimcalib_E_eta", "; E_{MC} (GeV); E_{rec,sim hit}/E_{MC}; #eta", 1500, 0.,
0106 150.0, 200, 0., 2.0, 50, 0, 5);
0107 hClusterSimNCells_E_eta =
0108 new TH3D("hClusterSimNCells_E_eta", "; E_{MC} (GeV); N_{cells, sim}; #eta", 1500, 0., 150.0,
0109 500, -0.5, 499.5, 50, 0, 5);
0110 hClusterESimcalib_E_phi =
0111 new TH3D("hClusterESimcalib_E_phi", "; E_{MC} (GeV); E_{rec,sim hit}/E_{MC}; #varphi (rad)",
0112 1500, 0., 150.0, 200, 0., 2.0, 360, -TMath::Pi(), TMath::Pi());
0113 hCellESim_layerX = new TH2D("hCellESim_layerX", "; #cell ID X; E_{rec,sim hit} (GeV)", 240, -0.5,
0114 239.5, 5000, 0, 1);
0115 hCellESim_layerY = new TH2D("hCellESim_layerY", "; #cell ID Y; E_{rec,sim hit} (GeV)", 240, -0.5,
0116 239.5, 5000, 0, 1);
0117 hCellESim_layerZ = new TH2D("hCellESim_layerZ", "; #cell ID Z; E_{rec,sim hit} (GeV)", 70, -0.5,
0118 69.5, 5000, 0, 1);
0119 hCellTSim_layerZ = new TH2D("hCellTSim_layerZ", "; #cell ID Z; t_{rec,sim hit} (GeV)", 70, -0.5,
0120 69.5, 5000, 0, 1000);
0121 hPosCaloSimHitsXY =
0122 new TH2D("hPosCaloSimHitsXY", "; X (cm); Y (cm)", 400, -400., 400., 400, -400., 400.);
0123 hPosCaloSimHitsZX =
0124 new TH2D("hPosCaloSimHitsZX", "; Z (cm); X (cm)", 200, 300., 500., 400, -400., 400.);
0125 hPosCaloSimHitsZY =
0126 new TH2D("hPosCaloSimHitsZY", "; Z (cm); Y (cm)", 200, 300., 500., 400, -400., 400.);
0127 hClusterESimcalib_E_eta->SetDirectory(m_dir_main);
0128 hClusterSimNCells_E_eta->SetDirectory(m_dir_main);
0129 hClusterESimcalib_E_phi->SetDirectory(m_dir_main);
0130 hCellESim_layerX->SetDirectory(m_dir_main);
0131 hCellESim_layerY->SetDirectory(m_dir_main);
0132 hCellESim_layerZ->SetDirectory(m_dir_main);
0133 hCellTSim_layerZ->SetDirectory(m_dir_main);
0134 hPosCaloSimHitsXY->SetDirectory(m_dir_main);
0135 hPosCaloSimHitsZX->SetDirectory(m_dir_main);
0136 hPosCaloSimHitsZY->SetDirectory(m_dir_main);
0137
0138
0139
0140
0141 hRecClusterEcalib_E_eta =
0142 new TH3D("hRecClusterEcalib_E_eta", "; E_{MC} (GeV); E_{rec,rec clus}/E_{MC}; #eta", 1500, 0.,
0143 150.0, 200, 0., 2.0, 50, 0, 5);
0144 hRecNClusters_E_eta = new TH3D("hRecNClusters_E_eta", "; E_{MC} (GeV); N_{rec cl.}; #eta", 1500,
0145 0., 150.0, 10, -0.5, 9.5, 50, 0, 5);
0146
0147 hRecClusterEcalib_Ehigh_eta =
0148 new TH3D("hRecClusterEcalib_Ehigh_eta", "; E_{MC} (GeV); E_{rec,rec clus high.}/E_{MC}; #eta",
0149 1500, 0., 150.0, 200, 0., 2.0, 50, 0, 5);
0150 hRecClusterNCells_Ehigh_eta =
0151 new TH3D("hRecClusterNCells_Ehigh_eta", "; E_{MC} (GeV); N_{cells, rec cl., high.}; #eta",
0152 1500, 0., 150.0, 500, -0.5, 499.5, 50, 0, 5);
0153 hRecClusterEcalib_E_eta->SetDirectory(m_dir_main);
0154 hRecNClusters_E_eta->SetDirectory(m_dir_main);
0155 hRecClusterEcalib_Ehigh_eta->SetDirectory(m_dir_main);
0156 hRecClusterNCells_Ehigh_eta->SetDirectory(m_dir_main);
0157
0158
0159
0160
0161 hRecFClusterEcalib_E_eta =
0162 new TH3D("hRecFClusterEcalib_E_eta", "; E_{MC} (GeV); E_{rec,island clus}/E_{MC}; #eta", 1500,
0163 0., 150.0, 200, 0., 2.0, 50, 0, 5);
0164 hRecFNClusters_E_eta = new TH3D("hRecFNClusters_E_eta", "; E_{MC} (GeV); N_{rec f. cl.}; #eta",
0165 1500, 0., 150.0, 10, -0.5, 9.5, 50, 0, 5);
0166
0167 hRecFClusterEcalib_Ehigh_eta = new TH3D("hRecFClusterEcalib_Ehigh_eta",
0168 "; E_{MC} (GeV); E_{rec,island clus high.}/E_{MC}; #eta",
0169 1500, 0., 150.0, 200, 0., 2.0, 50, 0, 5);
0170 hRecFClusterNCells_Ehigh_eta =
0171 new TH3D("hRecFClusterNCells_Ehigh_eta", "; E_{MC} (GeV); N_{cells, rec f. cl., high.}; #eta",
0172 1500, 0., 150.0, 500, -0.5, 499.5, 50, 0, 5);
0173 hRecFClusterEcalib_E_eta->SetDirectory(m_dir_main);
0174 hRecFNClusters_E_eta->SetDirectory(m_dir_main);
0175 hRecFClusterEcalib_Ehigh_eta->SetDirectory(m_dir_main);
0176 hRecFClusterNCells_Ehigh_eta->SetDirectory(m_dir_main);
0177
0178
0179
0180
0181 hRecFEmClusterEcalib_E_eta = new TH3D("hRecFEmClusterEcalib_E_eta",
0182 "; E_{MC} (GeV); E_{Ecal, rec,island clus}/E_{MC}; #eta",
0183 1500, 0., 150.0, 200, 0., 2.0, 50, 0, 5);
0184 hRecFEmNClusters_E_eta =
0185 new TH3D("hRecFEmNClusters_E_eta", "; E_{MC} (GeV); N_{Ecal, rec f. cl.}; #eta", 1500, 0.,
0186 150.0, 10, -0.5, 9.5, 50, 0, 5);
0187
0188 hRecFEmClusterEcalib_Ehigh_eta =
0189 new TH3D("hRecFEmClusterEcalib_Ehigh_eta",
0190 "; E_{MC} (GeV); E_{Ecal, rec,island clus high.}/E_{MC}; #eta", 1500, 0., 150.0, 200,
0191 0., 2.0, 50, 0, 5);
0192 hRecFEmClusterEcalib_E_eta->SetDirectory(m_dir_main);
0193 hRecFEmNClusters_E_eta->SetDirectory(m_dir_main);
0194 hRecFEmClusterEcalib_Ehigh_eta->SetDirectory(m_dir_main);
0195
0196
0197
0198
0199 hSamplingFractionEta = new TH2D("hSamplingFractionEta", "; #eta; f", 400, 1., 5., 500, 0., 0.2);
0200 hSamplingFractionEta->SetDirectory(m_dir_main);
0201
0202
0203
0204
0205 if (enableTree) {
0206 event_tree = new TTree("event_tree", "event_tree");
0207 event_tree->SetDirectory(m_dir_main);
0208
0209 t_lFHCal_towers_cellE = new float[maxNTowers];
0210 t_lFHCal_towers_cellT = new float[maxNTowers];
0211 t_lFHCal_towers_cellIDx = new short[maxNTowers];
0212 t_lFHCal_towers_cellIDy = new short[maxNTowers];
0213 t_lFHCal_towers_cellIDz = new short[maxNTowers];
0214 t_lFHCal_towers_clusterIDA = new short[maxNTowers];
0215 t_lFHCal_towers_clusterIDB = new short[maxNTowers];
0216 t_lFHCal_towers_cellTrueID = new int[maxNTowers];
0217
0218
0219 event_tree->Branch("tower_LFHCAL_N", &t_lFHCal_towers_N, "tower_LFHCAL_N/I");
0220 event_tree->Branch("tower_LFHCAL_E", t_lFHCal_towers_cellE, "tower_LFHCAL_E[tower_LFHCAL_N]/F");
0221 event_tree->Branch("tower_LFHCAL_T", t_lFHCal_towers_cellT, "tower_LFHCAL_T[tower_LFHCAL_N]/F");
0222 event_tree->Branch("tower_LFHCAL_ix", t_lFHCal_towers_cellIDx,
0223 "tower_LFHCAL_ix[tower_LFHCAL_N]/S");
0224 event_tree->Branch("tower_LFHCAL_iy", t_lFHCal_towers_cellIDy,
0225 "tower_LFHCAL_iy[tower_LFHCAL_N]/S");
0226 event_tree->Branch("tower_LFHCAL_iz", t_lFHCal_towers_cellIDz,
0227 "tower_LFHCAL_iz[tower_LFHCAL_N]/S");
0228 event_tree->Branch("tower_LFHCAL_clusIDA", t_lFHCal_towers_clusterIDA,
0229 "tower_LFHCAL_clusIDA[tower_LFHCAL_N]/S");
0230 event_tree->Branch("tower_LFHCAL_clusIDB", t_lFHCal_towers_clusterIDB,
0231 "tower_LFHCAL_clusIDB[tower_LFHCAL_N]/S");
0232 event_tree->Branch("tower_LFHCAL_trueID", t_lFHCal_towers_cellTrueID,
0233 "tower_LFHCAL_trueID[tower_LFHCAL_N]/I");
0234 }
0235
0236
0237
0238
0239 if (enableTreeCluster) {
0240 cluster_tree = new TTree("cluster_tree", "cluster_tree");
0241 cluster_tree->SetDirectory(m_dir_main);
0242
0243 t_mc_E = new float[maxNMC];
0244 t_mc_Phi = new float[maxNMC];
0245 t_mc_Eta = new float[maxNMC];
0246 t_lFHCal_cluster_E = new float[maxNCluster];
0247 t_lFHCal_cluster_NCells = new int[maxNCluster];
0248 t_lFHCal_cluster_Phi = new float[maxNCluster];
0249 t_lFHCal_cluster_Eta = new float[maxNCluster];
0250 t_fEMC_cluster_E = new float[maxNCluster];
0251 t_fEMC_cluster_NCells = new int[maxNCluster];
0252 t_fEMC_cluster_Eta = new float[maxNCluster];
0253 t_fEMC_cluster_Phi = new float[maxNCluster];
0254
0255
0256 cluster_tree->Branch("mc_N", &t_mc_N, "mc_N/I");
0257 cluster_tree->Branch("mc_E", t_mc_E, "mc_E[mc_N]/F");
0258 cluster_tree->Branch("mc_Phi", t_mc_Phi, "mc_Phi[mc_N]/F");
0259 cluster_tree->Branch("mc_Eta", t_mc_Eta, "mc_Eta[mc_N]/F");
0260
0261 cluster_tree->Branch("cluster_LFHCAL_N", &t_lFHCal_clusters_N, "cluster_LFHCAL_N/I");
0262 cluster_tree->Branch("cluster_LFHCAL_E", t_lFHCal_cluster_E,
0263 "cluster_LFHCAL_E[cluster_LFHCAL_N]/F");
0264 cluster_tree->Branch("cluster_LFHCAL_Ncells", t_lFHCal_cluster_NCells,
0265 "cluster_LFHCAL_Ncells[cluster_LFHCAL_N]/I");
0266 cluster_tree->Branch("cluster_LFHCAL_Eta", t_lFHCal_cluster_Eta,
0267 "cluster_LFHCAL_Eta[cluster_LFHCAL_N]/F");
0268 cluster_tree->Branch("cluster_LFHCAL_Phi", t_lFHCal_cluster_Phi,
0269 "cluster_LFHCAL_Phi[cluster_LFHCAL_N]/F");
0270
0271 cluster_tree->Branch("cluster_FECAL_N", &t_fEMC_clusters_N, "cluster_FECAL_N/I");
0272 cluster_tree->Branch("cluster_FECAL_E", t_fEMC_cluster_E, "cluster_FECAL_E[cluster_FECAL_N]/F");
0273 cluster_tree->Branch("cluster_FECAL_Ncells", t_fEMC_cluster_NCells,
0274 "cluster_FECAL_Ncells[cluster_FECAL_N]/I");
0275 cluster_tree->Branch("cluster_FECAL_Eta", t_fEMC_cluster_Eta,
0276 "cluster_FECAL_Eta[cluster_FECAL_N]/F");
0277 cluster_tree->Branch("cluster_FECAL_Phi", t_fEMC_cluster_Phi,
0278 "cluster_FECAL_Phi[cluster_FECAL_N]/F");
0279 }
0280
0281 std::cout << __PRETTY_FUNCTION__ << " " << __LINE__ << std::endl;
0282 auto detector = dd4hep_service->detector();
0283 std::cout << "--------------------------\nID specification:\n";
0284 try {
0285 m_decoder = detector->readout("LFHCALHits").idSpec().decoder();
0286 std::cout << "1st: " << m_decoder << std::endl;
0287 iLx = m_decoder->index("towerx");
0288 iLy = m_decoder->index("towery");
0289 iLz = m_decoder->index("layerz");
0290 iPassive = m_decoder->index("passive");
0291 std::cout << "full list: "
0292 << " " << m_decoder->fieldDescription() << std::endl;
0293 } catch (...) {
0294 std::cout << "2nd: " << m_decoder << std::endl;
0295 m_log->error("readoutClass not in the output");
0296 throw std::runtime_error("readoutClass not in the output.");
0297 }
0298 }
0299
0300
0301
0302
0303 void lfhcal_studiesProcessor::Process(const std::shared_ptr<const JEvent>& event) {
0304
0305
0306
0307
0308
0309 const auto& mcParticles = *(event->GetCollection<edm4hep::MCParticle>("MCParticles"));
0310 double mceta = 0;
0311 double mcphi = 0;
0312 double mcp = 0;
0313 double mcenergy = 0;
0314 int iMC = 0;
0315 for (auto mcparticle : mcParticles) {
0316 if (mcparticle.getGeneratorStatus() != 1) {
0317 continue;
0318 }
0319 const auto& mom = mcparticle.getMomentum();
0320
0321 mcenergy = mcparticle.getEnergy();
0322
0323 mceta = -log(tan(atan2(sqrt(mom.x * mom.x + mom.y * mom.y), mom.z) / 2.));
0324
0325 mcphi = atan2(mom.y, mom.x);
0326
0327 mcp = sqrt(mom.x * mom.x + mom.y * mom.y + mom.z * mom.z);
0328 m_log->trace("MC particle:{} \t {} \t {} \t totmom: {} phi {} eta {}", mom.x, mom.y, mom.z, mcp,
0329 mcphi, mceta);
0330 hMCEnergyVsEta->Fill(mcp, mceta);
0331
0332 if (enableTreeCluster) {
0333 if (iMC < maxNMC) {
0334 t_mc_E[iMC] = mcenergy;
0335 t_mc_Phi[iMC] = mcphi;
0336 t_mc_Eta[iMC] = mceta;
0337 }
0338 }
0339 iMC++;
0340 }
0341 if (enableTreeCluster) {
0342 t_mc_N = iMC;
0343 }
0344
0345
0346
0347 std::vector<towersStrct> input_tower_sim;
0348 int nCaloHitsSim = 0;
0349 float sumActiveCaloEnergy = 0;
0350 float sumPassiveCaloEnergy = 0;
0351 const auto& simHits = *(event->GetCollection<edm4hep::SimCalorimeterHit>(nameSimHits));
0352 for (const auto caloHit : simHits) {
0353 float x = caloHit.getPosition().x / 10.;
0354 float y = caloHit.getPosition().y / 10.;
0355 float z = caloHit.getPosition().z / 10.;
0356 uint64_t cellID = caloHit.getCellID();
0357 float energy = caloHit.getEnergy();
0358 double time = std::numeric_limits<double>::max();
0359 for (const auto& c : caloHit.getContributions()) {
0360 time = std::min<double>(c.getTime(), time);
0361 }
0362
0363 auto detector_module_x = m_decoder->get(cellID, 1);
0364 auto detector_module_y = m_decoder->get(cellID, 2);
0365 auto detector_passive = m_decoder->get(cellID, iPassive);
0366 auto detector_layer_x = m_decoder->get(cellID, iLx);
0367 auto detector_layer_y = m_decoder->get(cellID, iLy);
0368 long detector_layer_rz = -1;
0369 if (isLFHCal) {
0370 detector_layer_rz = m_decoder->get(cellID, 7);
0371 }
0372 auto detector_layer_z = m_decoder->get(cellID, iLz);
0373 if (detector_passive == 0) {
0374 sumActiveCaloEnergy += energy;
0375 } else {
0376 sumPassiveCaloEnergy += energy;
0377 }
0378
0379 if (detector_passive > 0) {
0380 continue;
0381 }
0382
0383 long cellIDx = -1;
0384 long cellIDy = -1;
0385 long cellIDz = -1;
0386 if (isLFHCal) {
0387 cellIDx = 54LL * 2 - detector_module_x * 2 + detector_layer_x;
0388 cellIDy = 54LL * 2 - detector_module_y * 2 + detector_layer_y;
0389 cellIDz = detector_layer_rz * 10 + detector_layer_z;
0390 }
0391 nCaloHitsSim++;
0392
0393 hPosCaloSimHitsXY->Fill(x, y);
0394 hPosCaloSimHitsZX->Fill(z, x);
0395 hPosCaloSimHitsZY->Fill(z, y);
0396
0397 hCellESim_layerZ->Fill(cellIDz, energy);
0398 hCellESim_layerX->Fill(cellIDx, energy);
0399 hCellESim_layerY->Fill(cellIDy, energy);
0400 hCellTSim_layerZ->Fill(cellIDz, time);
0401
0402
0403 bool found = false;
0404 for (auto& tower : input_tower_sim) {
0405 if (tower.cellID == static_cast<decltype(tower.cellID)>(cellID)) {
0406 tower.energy += energy;
0407 found = true;
0408 break;
0409 }
0410 }
0411 if (!found) {
0412 towersStrct tempstructT;
0413 tempstructT.energy = energy;
0414 tempstructT.time = time;
0415 tempstructT.posx = x;
0416 tempstructT.posy = y;
0417 tempstructT.posz = z;
0418 tempstructT.cellID = cellID;
0419 tempstructT.cellIDx = cellIDx;
0420 tempstructT.cellIDy = cellIDy;
0421 tempstructT.cellIDz = cellIDz;
0422 tempstructT.tower_trueID = 0;
0423 input_tower_sim.push_back(tempstructT);
0424 }
0425 }
0426
0427
0428
0429
0430 const auto& recHits = *(event->GetCollection<edm4eic::CalorimeterHit>(nameRecHits));
0431 int nCaloHitsRec = 0;
0432 std::vector<towersStrct> input_tower_rec;
0433 std::vector<towersStrct> input_tower_recSav;
0434
0435 for (const auto caloHit : recHits) {
0436 float x = caloHit.getPosition().x / 10.;
0437 float y = caloHit.getPosition().y / 10.;
0438 float z = caloHit.getPosition().z / 10.;
0439 uint64_t cellID = caloHit.getCellID();
0440 float energy = caloHit.getEnergy();
0441 float time = caloHit.getTime();
0442
0443 auto detector_module_x = m_decoder->get(cellID, 1);
0444 auto detector_module_y = m_decoder->get(cellID, 2);
0445 auto detector_passive = m_decoder->get(cellID, iPassive);
0446 auto detector_layer_x = m_decoder->get(cellID, iLx);
0447 auto detector_layer_y = m_decoder->get(cellID, iLy);
0448 int detector_layer_rz = -1;
0449 if (isLFHCal) {
0450 detector_layer_rz = m_decoder->get(cellID, 7);
0451 }
0452
0453 if (detector_passive > 0) {
0454 continue;
0455 }
0456
0457
0458 long cellIDx = -1;
0459 long cellIDy = -1;
0460 if (isLFHCal) {
0461 cellIDx = 54LL * 2 - detector_module_x * 2 + detector_layer_x;
0462 cellIDy = 54LL * 2 - detector_module_y * 2 + detector_layer_y;
0463 }
0464
0465 hPosCaloHitsXY->Fill(x, y);
0466 hPosCaloHitsZX->Fill(z, x);
0467 hPosCaloHitsZY->Fill(z, y);
0468
0469 nCaloHitsRec++;
0470
0471
0472 bool found = false;
0473 for (auto& tower : input_tower_rec) {
0474 if (tower.cellID == static_cast<decltype(tower.cellID)>(cellID)) {
0475 tower.energy += energy;
0476 found = true;
0477 break;
0478 }
0479 }
0480 if (!found) {
0481 towersStrct tempstructT;
0482 tempstructT.energy = energy;
0483 tempstructT.time = time;
0484 tempstructT.posx = x;
0485 tempstructT.posy = y;
0486 tempstructT.posz = z;
0487 tempstructT.cellID = cellID;
0488 tempstructT.cellIDx = cellIDx;
0489 tempstructT.cellIDy = cellIDy;
0490 if (isLFHCal) {
0491 tempstructT.cellIDz = detector_layer_rz;
0492 }
0493 tempstructT.tower_trueID = 0;
0494 input_tower_rec.push_back(tempstructT);
0495 input_tower_recSav.push_back(tempstructT);
0496 }
0497 }
0498 m_log->trace("LFHCal mod: nCaloHits sim {}\t rec {}", nCaloHitsSim, nCaloHitsRec);
0499
0500 if (nCaloHitsRec > 0) {
0501 nEventsWithCaloHits++;
0502 }
0503
0504
0505
0506
0507 hSamplingFractionEta->Fill(mceta,
0508 sumActiveCaloEnergy / (sumActiveCaloEnergy + sumPassiveCaloEnergy));
0509 std::ranges::sort(input_tower_rec, &acompare);
0510 std::ranges::sort(input_tower_recSav, &acompare);
0511 std::ranges::sort(input_tower_sim, &acompare);
0512
0513
0514
0515
0516
0517
0518 double tot_energyRecHit = 0;
0519 for (auto& tower : input_tower_rec) {
0520 tower.energy = tower.energy;
0521 tot_energyRecHit += tower.energy;
0522 }
0523
0524 double samplingFractionFe = 0.037;
0525 double samplingFractionW = 0.019;
0526 int minCellIDzDiffSamp = 5;
0527
0528 double tot_energySimHit = 0;
0529 for (auto& tower : input_tower_sim) {
0530 if (tower.cellIDz < minCellIDzDiffSamp) {
0531 tower.energy = tower.energy / samplingFractionW;
0532 } else {
0533 tower.energy = tower.energy / samplingFractionFe;
0534 }
0535 tot_energySimHit += tower.energy;
0536 }
0537 m_log->trace("Mc E: {} \t eta: {} \t sim E rec: {}\t rec E rec: {}", mcenergy, mceta,
0538 tot_energySimHit, tot_energyRecHit);
0539
0540
0541
0542
0543
0544 hClusterNCells_E_eta->Fill(mcenergy, nCaloHitsRec, mceta);
0545 hClusterEcalib_E_eta->Fill(mcenergy, tot_energyRecHit / mcenergy, mceta);
0546 hClusterEcalib_E_phi->Fill(mcenergy, tot_energyRecHit / mcenergy, mcphi);
0547
0548 hClusterSimNCells_E_eta->Fill(mcenergy, nCaloHitsSim, mceta);
0549 hClusterESimcalib_E_eta->Fill(mcenergy, tot_energySimHit / mcenergy, mceta);
0550 hClusterESimcalib_E_phi->Fill(mcenergy, tot_energySimHit / mcenergy, mcphi);
0551
0552
0553
0554
0555 int removedCells = 0;
0556 float minAggE = 0.001;
0557 float seedE = 0.100;
0558
0559 if (!input_tower_rec.empty()) {
0560
0561
0562 while (input_tower_rec.at(input_tower_rec.size() - 1).energy < minAggE) {
0563 input_tower_rec.pop_back();
0564 removedCells++;
0565 }
0566 m_log->trace("removed {} with E < {} GeV", removedCells, minAggE);
0567
0568 int nclusters = 0;
0569
0570 std::vector<clustersStrct> clusters_calo;
0571
0572 std::vector<towersStrct> cluster_towers;
0573 while (!input_tower_rec.empty()) {
0574 cluster_towers.clear();
0575 clustersStrct tempstructC;
0576
0577 if (input_tower_rec.at(0).energy > seedE) {
0578 m_log->trace("seed: {}\t {} \t {}", input_tower_rec.at(0).energy,
0579 input_tower_rec.at(0).cellIDx, input_tower_rec.at(0).cellIDy,
0580 input_tower_rec.at(0).cellIDz);
0581 tempstructC = findMACluster(seedE, minAggE, input_tower_rec, cluster_towers);
0582
0583
0584 float* showershape_eta_phi =
0585 CalculateM02andWeightedPosition(cluster_towers, tempstructC.cluster_E, 4.5);
0586
0587 tempstructC.cluster_M02 = showershape_eta_phi[0];
0588 tempstructC.cluster_M20 = showershape_eta_phi[1];
0589 tempstructC.cluster_Eta = showershape_eta_phi[2];
0590 tempstructC.cluster_Phi = showershape_eta_phi[3];
0591 tempstructC.cluster_X = showershape_eta_phi[4];
0592 tempstructC.cluster_Y = showershape_eta_phi[5];
0593 tempstructC.cluster_Z = showershape_eta_phi[6];
0594
0595 tempstructC.cluster_towers = cluster_towers;
0596 m_log->trace("---------> \t {} \tcluster with E = {} \tEta: {} \tPhi: {} \tX: {} \tY: {} "
0597 "\tZ: {} \tntowers: {} \ttrueID: {}",
0598 nclusters, tempstructC.cluster_E, tempstructC.cluster_Eta,
0599 tempstructC.cluster_Phi, tempstructC.cluster_X, tempstructC.cluster_Y,
0600 tempstructC.cluster_Z, tempstructC.cluster_NTowers,
0601 tempstructC.cluster_trueID);
0602 clusters_calo.push_back(tempstructC);
0603
0604 clusters_calo.push_back(tempstructC);
0605
0606 nclusters++;
0607 } else {
0608 m_log->trace("remaining: {} largest: {} \t {} \t {} \t {}", (int)input_tower_rec.size(),
0609 input_tower_rec.at(0).energy, input_tower_rec.at(0).cellIDx,
0610 input_tower_rec.at(0).cellIDy, input_tower_rec.at(0).cellIDz);
0611 input_tower_rec.clear();
0612 }
0613 }
0614
0615
0616
0617
0618 std::ranges::sort(clusters_calo, &acompareCl);
0619 m_log->info("-----> found {} clusters", clusters_calo.size());
0620 hRecNClusters_E_eta->Fill(mcenergy, clusters_calo.size(), mceta);
0621 int iCl = 0;
0622 for (const auto& cluster : clusters_calo) {
0623 if (iCl < maxNCluster && enableTreeCluster) {
0624 t_lFHCal_cluster_E[iCl] = (float)cluster.cluster_E;
0625 t_lFHCal_cluster_NCells[iCl] = (int)cluster.cluster_NTowers;
0626 t_lFHCal_cluster_Eta[iCl] = (float)cluster.cluster_Eta;
0627 t_lFHCal_cluster_Phi[iCl] = (float)cluster.cluster_Phi;
0628 }
0629 hRecClusterEcalib_E_eta->Fill(mcenergy, cluster.cluster_E / mcenergy, mceta);
0630 for (const auto& cluster_tower : cluster.cluster_towers) {
0631 int pSav = 0;
0632 while (cluster_tower.cellID != input_tower_recSav.at(pSav).cellID &&
0633 pSav < (int)input_tower_recSav.size()) {
0634 pSav++;
0635 }
0636 if (cluster_tower.cellID == input_tower_recSav.at(pSav).cellID) {
0637 input_tower_recSav.at(pSav).tower_clusterIDA = iCl;
0638 }
0639 }
0640
0641 if (iCl == 0) {
0642 hRecClusterEcalib_Ehigh_eta->Fill(mcenergy, cluster.cluster_E / mcenergy, mceta);
0643 hRecClusterNCells_Ehigh_eta->Fill(mcenergy, cluster.cluster_NTowers, mceta);
0644 }
0645 iCl++;
0646 m_log->trace("MA cluster {}:\t {} \t {}", iCl, cluster.cluster_E, cluster.cluster_NTowers);
0647 }
0648 if (iCl < maxNCluster && enableTreeCluster) {
0649 t_lFHCal_clusters_N = iCl;
0650 }
0651
0652 clusters_calo.clear();
0653 } else {
0654 hRecNClusters_E_eta->Fill(mcenergy, 0., mceta);
0655 if (enableTreeCluster) {
0656 t_lFHCal_clusters_N = 0;
0657 }
0658 }
0659
0660
0661
0662
0663 int iClF = 0;
0664 float highestEFr = 0;
0665 int iClFHigh = 0;
0666
0667 const auto& lfhcalClustersF = *(event->GetCollection<edm4eic::Cluster>(nameClusters));
0668 for (const auto cluster : lfhcalClustersF) {
0669 if (cluster.getEnergy() > highestEFr) {
0670 iClFHigh = iClF;
0671 highestEFr = cluster.getEnergy();
0672 }
0673 hRecFClusterEcalib_E_eta->Fill(mcenergy, cluster.getEnergy() / mcenergy, mceta);
0674 m_log->trace("Island cluster {}:\t {} \t {}", iClF, cluster.getEnergy(), cluster.getNhits());
0675 iClF++;
0676 }
0677 hRecFNClusters_E_eta->Fill(mcenergy, iClF, mceta);
0678
0679 iClF = 0;
0680 for (const auto cluster : lfhcalClustersF) {
0681 if (iClF == iClFHigh) {
0682 hRecFClusterEcalib_Ehigh_eta->Fill(mcenergy, cluster.getEnergy() / mcenergy, mceta);
0683 hRecFClusterNCells_Ehigh_eta->Fill(mcenergy, cluster.getNhits(), mceta);
0684 }
0685 iClF++;
0686 }
0687
0688
0689
0690
0691 int iECl = 0;
0692 float highestEEmCl = 0;
0693 int iEClHigh = 0;
0694
0695 if (enableECalCluster) {
0696 try {
0697 const auto& fEMCClustersF = *(event->GetCollection<edm4eic::Cluster>("EcalEndcapPClusters"));
0698 m_log->info("-----> found fEMCClustersF:", fEMCClustersF.size());
0699 for (const auto cluster : fEMCClustersF) {
0700 if (iECl < maxNCluster && enableTreeCluster) {
0701 t_fEMC_cluster_E[iECl] = cluster.getEnergy();
0702 t_fEMC_cluster_NCells[iECl] = (int)cluster.getNhits();
0703 t_fEMC_cluster_Eta[iECl] = (-1.) * std::log(std::tan(cluster.getIntrinsicTheta() / 2.));
0704 t_fEMC_cluster_Phi[iECl] = cluster.getIntrinsicPhi();
0705 }
0706
0707 if (cluster.getEnergy() > highestEEmCl) {
0708 iEClHigh = iECl;
0709 highestEEmCl = cluster.getEnergy();
0710 }
0711 hRecFEmClusterEcalib_E_eta->Fill(mcenergy, cluster.getEnergy() / mcenergy, mceta);
0712 iECl++;
0713 }
0714 t_fEMC_clusters_N = iECl;
0715 hRecFEmNClusters_E_eta->Fill(mcenergy, iECl, mceta);
0716
0717
0718 iECl = 0;
0719 for (const auto cluster : fEMCClustersF) {
0720 if (iECl == iEClHigh) {
0721 hRecFEmClusterEcalib_Ehigh_eta->Fill(mcenergy, cluster.getEnergy() / mcenergy, mceta);
0722 }
0723 iECl++;
0724 }
0725 } catch (...) {
0726 std::cout << "ECal clusters not in output" << std::endl;
0727 enableECalCluster = false;
0728 }
0729 }
0730
0731
0732
0733 if (enableTree) {
0734 t_lFHCal_towers_N = (int)input_tower_recSav.size();
0735 for (int iCell = 0; iCell < (int)input_tower_recSav.size(); iCell++) {
0736 m_log->trace("{} \t {} \t {} \t {} \t {} \t {}", input_tower_recSav.at(iCell).cellIDx,
0737 input_tower_recSav.at(iCell).cellIDy, input_tower_recSav.at(iCell).cellIDz,
0738 input_tower_recSav.at(iCell).energy,
0739 input_tower_recSav.at(iCell).tower_clusterIDA,
0740 input_tower_recSav.at(iCell).tower_clusterIDB);
0741
0742 t_lFHCal_towers_cellE[iCell] = input_tower_recSav.at(iCell).energy;
0743 t_lFHCal_towers_cellT[iCell] = input_tower_recSav.at(iCell).time;
0744 t_lFHCal_towers_cellIDx[iCell] = (short)input_tower_recSav.at(iCell).cellIDx;
0745 t_lFHCal_towers_cellIDy[iCell] = (short)input_tower_recSav.at(iCell).cellIDy;
0746 t_lFHCal_towers_cellIDz[iCell] = (short)input_tower_recSav.at(iCell).cellIDz;
0747 t_lFHCal_towers_clusterIDA[iCell] = (short)input_tower_recSav.at(iCell).tower_clusterIDA;
0748 t_lFHCal_towers_clusterIDB[iCell] = (short)input_tower_recSav.at(iCell).tower_clusterIDB;
0749 t_lFHCal_towers_cellTrueID[iCell] = input_tower_recSav.at(iCell).tower_trueID;
0750 }
0751
0752 event_tree->Fill();
0753
0754 t_lFHCal_towers_N = 0;
0755 for (Int_t itow = 0; itow < maxNTowers; itow++) {
0756 t_lFHCal_towers_cellE[itow] = 0;
0757 t_lFHCal_towers_cellT[itow] = 0;
0758 t_lFHCal_towers_cellIDx[itow] = 0;
0759 t_lFHCal_towers_cellIDy[itow] = 0;
0760 t_lFHCal_towers_cellIDz[itow] = 0;
0761 t_lFHCal_towers_clusterIDA[itow] = 0;
0762 t_lFHCal_towers_clusterIDB[itow] = 0;
0763 t_lFHCal_towers_cellTrueID[itow] = 0;
0764 }
0765 }
0766
0767
0768
0769
0770 if (enableTreeCluster) {
0771 cluster_tree->Fill();
0772
0773 t_mc_N = 0;
0774 t_lFHCal_clusters_N = 0;
0775 t_fEMC_clusters_N = 0;
0776 for (Int_t iMC = 0; iMC < maxNMC; iMC++) {
0777 t_mc_E[iMC] = 0;
0778 t_mc_Phi[iMC] = 0;
0779 t_mc_Eta[iMC] = 0;
0780 }
0781 for (Int_t iCl = 0; iCl < maxNCluster; iCl++) {
0782 t_lFHCal_cluster_E[iCl] = 0;
0783 t_lFHCal_cluster_NCells[iCl] = 0;
0784 t_lFHCal_cluster_Eta[iCl] = 0;
0785 t_lFHCal_cluster_Phi[iCl] = 0;
0786 t_fEMC_cluster_E[iCl] = 0;
0787 t_fEMC_cluster_NCells[iCl] = 0;
0788 t_fEMC_cluster_Eta[iCl] = 0;
0789 t_fEMC_cluster_Phi[iCl] = 0;
0790 }
0791 }
0792 }
0793
0794
0795
0796
0797 void lfhcal_studiesProcessor::Finish() {
0798 std::cout << "------> LFHCal " << nEventsWithCaloHits << " with calo info present" << std::endl;
0799
0800
0801 if (enableTree) {
0802 delete[] t_lFHCal_towers_cellE;
0803 delete[] t_lFHCal_towers_cellT;
0804 delete[] t_lFHCal_towers_cellIDx;
0805 delete[] t_lFHCal_towers_cellIDy;
0806 delete[] t_lFHCal_towers_cellIDz;
0807 delete[] t_lFHCal_towers_clusterIDA;
0808 delete[] t_lFHCal_towers_clusterIDB;
0809 delete[] t_lFHCal_towers_cellTrueID;
0810 }
0811
0812 if (enableTreeCluster) {
0813 delete[] t_mc_E;
0814 delete[] t_mc_Phi;
0815 delete[] t_mc_Eta;
0816 delete[] t_lFHCal_cluster_E;
0817 delete[] t_lFHCal_cluster_NCells;
0818 delete[] t_lFHCal_cluster_Phi;
0819 delete[] t_lFHCal_cluster_Eta;
0820 delete[] t_fEMC_cluster_E;
0821 delete[] t_fEMC_cluster_NCells;
0822 delete[] t_fEMC_cluster_Eta;
0823 delete[] t_fEMC_cluster_Phi;
0824 }
0825 }
