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