File indexing completed on 2025-07-11 07:53:39
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/JEvent.h>
0014 #include <JANA/Services/JGlobalRootLock.h>
0015 #include <RtypesCore.h>
0016 #include <TMath.h>
0017 #include <algorithm>
0018 #include <cmath>
0019 #include <cstdint>
0020 #include <edm4eic/CalorimeterHitCollection.h>
0021 #include <edm4eic/ClusterCollection.h>
0022 #include <edm4hep/CaloHitContributionCollection.h>
0023 #include <edm4hep/MCParticleCollection.h>
0024 #include <edm4hep/SimCalorimeterHitCollection.h>
0025 #include <edm4hep/Vector3d.h>
0026 #include <edm4hep/Vector3f.h>
0027 #include <fmt/core.h>
0028 #include <gsl/pointers>
0029 #include <iostream>
0030 #include <limits>
0031 #include <map>
0032 #include <podio/RelationRange.h>
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 if (c.getTime() <= time) {
0361 time = c.getTime();
0362 }
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::sort(input_tower_rec.begin(), input_tower_rec.end(), &acompare);
0512 std::sort(input_tower_recSav.begin(), input_tower_recSav.end(), &acompare);
0513 std::sort(input_tower_sim.begin(), input_tower_sim.end(), &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::sort(clusters_calo.begin(), clusters_calo.end(), &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 = (int)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] = (float)cluster.getEnergy();
0704 t_fEMC_cluster_NCells[iECl] = (int)cluster.getNhits();
0705 t_fEMC_cluster_Eta[iECl] =
0706 (-1.) * std::log(std::tan((float)cluster.getIntrinsicTheta() / 2.));
0707 t_fEMC_cluster_Phi[iECl] = (float)cluster.getIntrinsicPhi();
0708 }
0709
0710 if (cluster.getEnergy() > highestEEmCl) {
0711 iEClHigh = iECl;
0712 highestEEmCl = cluster.getEnergy();
0713 }
0714 hRecFEmClusterEcalib_E_eta->Fill(mcenergy, cluster.getEnergy() / mcenergy, mceta);
0715 iECl++;
0716 }
0717 t_fEMC_clusters_N = iECl;
0718 hRecFEmNClusters_E_eta->Fill(mcenergy, iECl, mceta);
0719
0720
0721 iECl = 0;
0722 for (const auto cluster : fEMCClustersF) {
0723 if (iECl == iEClHigh) {
0724 hRecFEmClusterEcalib_Ehigh_eta->Fill(mcenergy, cluster.getEnergy() / mcenergy, mceta);
0725 }
0726 iECl++;
0727 }
0728 } catch (...) {
0729 std::cout << "ECal clusters not in output" << std::endl;
0730 enableECalCluster = false;
0731 }
0732 }
0733
0734
0735
0736 if (enableTree) {
0737 t_lFHCal_towers_N = (int)input_tower_recSav.size();
0738 for (int iCell = 0; iCell < (int)input_tower_recSav.size(); iCell++) {
0739 m_log->trace("{} \t {} \t {} \t {} \t {} \t {}", input_tower_recSav.at(iCell).cellIDx,
0740 input_tower_recSav.at(iCell).cellIDy, input_tower_recSav.at(iCell).cellIDz,
0741 input_tower_recSav.at(iCell).energy,
0742 input_tower_recSav.at(iCell).tower_clusterIDA,
0743 input_tower_recSav.at(iCell).tower_clusterIDB);
0744
0745 t_lFHCal_towers_cellE[iCell] = (float)input_tower_recSav.at(iCell).energy;
0746 t_lFHCal_towers_cellT[iCell] = (float)input_tower_recSav.at(iCell).time;
0747 t_lFHCal_towers_cellIDx[iCell] = (short)input_tower_recSav.at(iCell).cellIDx;
0748 t_lFHCal_towers_cellIDy[iCell] = (short)input_tower_recSav.at(iCell).cellIDy;
0749 t_lFHCal_towers_cellIDz[iCell] = (short)input_tower_recSav.at(iCell).cellIDz;
0750 t_lFHCal_towers_clusterIDA[iCell] = (short)input_tower_recSav.at(iCell).tower_clusterIDA;
0751 t_lFHCal_towers_clusterIDB[iCell] = (short)input_tower_recSav.at(iCell).tower_clusterIDB;
0752 t_lFHCal_towers_cellTrueID[iCell] = (int)input_tower_recSav.at(iCell).tower_trueID;
0753 }
0754
0755 event_tree->Fill();
0756
0757 t_lFHCal_towers_N = 0;
0758 for (Int_t itow = 0; itow < maxNTowers; itow++) {
0759 t_lFHCal_towers_cellE[itow] = 0;
0760 t_lFHCal_towers_cellT[itow] = 0;
0761 t_lFHCal_towers_cellIDx[itow] = 0;
0762 t_lFHCal_towers_cellIDy[itow] = 0;
0763 t_lFHCal_towers_cellIDz[itow] = 0;
0764 t_lFHCal_towers_clusterIDA[itow] = 0;
0765 t_lFHCal_towers_clusterIDB[itow] = 0;
0766 t_lFHCal_towers_cellTrueID[itow] = 0;
0767 }
0768 }
0769
0770
0771
0772
0773 if (enableTreeCluster) {
0774 cluster_tree->Fill();
0775
0776 t_mc_N = 0;
0777 t_lFHCal_clusters_N = 0;
0778 t_fEMC_clusters_N = 0;
0779 for (Int_t iMC = 0; iMC < maxNMC; iMC++) {
0780 t_mc_E[iMC] = 0;
0781 t_mc_Phi[iMC] = 0;
0782 t_mc_Eta[iMC] = 0;
0783 }
0784 for (Int_t iCl = 0; iCl < maxNCluster; iCl++) {
0785 t_lFHCal_cluster_E[iCl] = 0;
0786 t_lFHCal_cluster_NCells[iCl] = 0;
0787 t_lFHCal_cluster_Eta[iCl] = 0;
0788 t_lFHCal_cluster_Phi[iCl] = 0;
0789 t_fEMC_cluster_E[iCl] = 0;
0790 t_fEMC_cluster_NCells[iCl] = 0;
0791 t_fEMC_cluster_Eta[iCl] = 0;
0792 t_fEMC_cluster_Phi[iCl] = 0;
0793 }
0794 }
0795 }
0796
0797
0798
0799
0800 void lfhcal_studiesProcessor::Finish() {
0801 std::cout << "------> LFHCal " << nEventsWithCaloHits << " with calo info present" << std::endl;
0802
0803
0804 if (enableTree) {
0805 delete[] t_lFHCal_towers_cellE;
0806 delete[] t_lFHCal_towers_cellT;
0807 delete[] t_lFHCal_towers_cellIDx;
0808 delete[] t_lFHCal_towers_cellIDy;
0809 delete[] t_lFHCal_towers_cellIDz;
0810 delete[] t_lFHCal_towers_clusterIDA;
0811 delete[] t_lFHCal_towers_clusterIDB;
0812 delete[] t_lFHCal_towers_cellTrueID;
0813 }
0814
0815 if (enableTreeCluster) {
0816 delete[] t_mc_E;
0817 delete[] t_mc_Phi;
0818 delete[] t_mc_Eta;
0819 delete[] t_lFHCal_cluster_E;
0820 delete[] t_lFHCal_cluster_NCells;
0821 delete[] t_lFHCal_cluster_Phi;
0822 delete[] t_lFHCal_cluster_Eta;
0823 delete[] t_fEMC_cluster_E;
0824 delete[] t_fEMC_cluster_NCells;
0825 delete[] t_fEMC_cluster_Eta;
0826 delete[] t_fEMC_cluster_Phi;
0827 }
0828 }
