File indexing completed on 2025-09-17 08:07:26
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 if (c.getTime() <= time) {
0363 time = c.getTime();
0364 }
0365 }
0366
0367 auto detector_module_x = m_decoder->get(cellID, 1);
0368 auto detector_module_y = m_decoder->get(cellID, 2);
0369 auto detector_passive = m_decoder->get(cellID, iPassive);
0370 auto detector_layer_x = m_decoder->get(cellID, iLx);
0371 auto detector_layer_y = m_decoder->get(cellID, iLy);
0372 long detector_layer_rz = -1;
0373 if (isLFHCal) {
0374 detector_layer_rz = m_decoder->get(cellID, 7);
0375 }
0376 auto detector_layer_z = m_decoder->get(cellID, iLz);
0377 if (detector_passive == 0) {
0378 sumActiveCaloEnergy += energy;
0379 } else {
0380 sumPassiveCaloEnergy += energy;
0381 }
0382
0383 if (detector_passive > 0) {
0384 continue;
0385 }
0386
0387 long cellIDx = -1;
0388 long cellIDy = -1;
0389 long cellIDz = -1;
0390 if (isLFHCal) {
0391 cellIDx = 54LL * 2 - detector_module_x * 2 + detector_layer_x;
0392 cellIDy = 54LL * 2 - detector_module_y * 2 + detector_layer_y;
0393 cellIDz = detector_layer_rz * 10 + detector_layer_z;
0394 }
0395 nCaloHitsSim++;
0396
0397 hPosCaloSimHitsXY->Fill(x, y);
0398 hPosCaloSimHitsZX->Fill(z, x);
0399 hPosCaloSimHitsZY->Fill(z, y);
0400
0401 hCellESim_layerZ->Fill(cellIDz, energy);
0402 hCellESim_layerX->Fill(cellIDx, energy);
0403 hCellESim_layerY->Fill(cellIDy, energy);
0404 hCellTSim_layerZ->Fill(cellIDz, time);
0405
0406
0407 bool found = false;
0408 for (auto& tower : input_tower_sim) {
0409 if (tower.cellID == static_cast<decltype(tower.cellID)>(cellID)) {
0410 tower.energy += energy;
0411 found = true;
0412 break;
0413 }
0414 }
0415 if (!found) {
0416 towersStrct tempstructT;
0417 tempstructT.energy = energy;
0418 tempstructT.time = time;
0419 tempstructT.posx = x;
0420 tempstructT.posy = y;
0421 tempstructT.posz = z;
0422 tempstructT.cellID = cellID;
0423 tempstructT.cellIDx = cellIDx;
0424 tempstructT.cellIDy = cellIDy;
0425 tempstructT.cellIDz = cellIDz;
0426 tempstructT.tower_trueID = 0;
0427 input_tower_sim.push_back(tempstructT);
0428 }
0429 }
0430
0431
0432
0433
0434 const auto& recHits = *(event->GetCollection<edm4eic::CalorimeterHit>(nameRecHits));
0435 int nCaloHitsRec = 0;
0436 std::vector<towersStrct> input_tower_rec;
0437 std::vector<towersStrct> input_tower_recSav;
0438
0439 for (const auto caloHit : recHits) {
0440 float x = caloHit.getPosition().x / 10.;
0441 float y = caloHit.getPosition().y / 10.;
0442 float z = caloHit.getPosition().z / 10.;
0443 uint64_t cellID = caloHit.getCellID();
0444 float energy = caloHit.getEnergy();
0445 float time = caloHit.getTime();
0446
0447 auto detector_module_x = m_decoder->get(cellID, 1);
0448 auto detector_module_y = m_decoder->get(cellID, 2);
0449 auto detector_passive = m_decoder->get(cellID, iPassive);
0450 auto detector_layer_x = m_decoder->get(cellID, iLx);
0451 auto detector_layer_y = m_decoder->get(cellID, iLy);
0452 int detector_layer_rz = -1;
0453 if (isLFHCal) {
0454 detector_layer_rz = m_decoder->get(cellID, 7);
0455 }
0456
0457 if (detector_passive > 0) {
0458 continue;
0459 }
0460
0461
0462 long cellIDx = -1;
0463 long cellIDy = -1;
0464 if (isLFHCal) {
0465 cellIDx = 54LL * 2 - detector_module_x * 2 + detector_layer_x;
0466 cellIDy = 54LL * 2 - detector_module_y * 2 + detector_layer_y;
0467 }
0468
0469 hPosCaloHitsXY->Fill(x, y);
0470 hPosCaloHitsZX->Fill(z, x);
0471 hPosCaloHitsZY->Fill(z, y);
0472
0473 nCaloHitsRec++;
0474
0475
0476 bool found = false;
0477 for (auto& tower : input_tower_rec) {
0478 if (tower.cellID == static_cast<decltype(tower.cellID)>(cellID)) {
0479 tower.energy += energy;
0480 found = true;
0481 break;
0482 }
0483 }
0484 if (!found) {
0485 towersStrct tempstructT;
0486 tempstructT.energy = energy;
0487 tempstructT.time = time;
0488 tempstructT.posx = x;
0489 tempstructT.posy = y;
0490 tempstructT.posz = z;
0491 tempstructT.cellID = cellID;
0492 tempstructT.cellIDx = cellIDx;
0493 tempstructT.cellIDy = cellIDy;
0494 if (isLFHCal) {
0495 tempstructT.cellIDz = detector_layer_rz;
0496 }
0497 tempstructT.tower_trueID = 0;
0498 input_tower_rec.push_back(tempstructT);
0499 input_tower_recSav.push_back(tempstructT);
0500 }
0501 }
0502 m_log->trace("LFHCal mod: nCaloHits sim {}\t rec {}", nCaloHitsSim, nCaloHitsRec);
0503
0504 if (nCaloHitsRec > 0) {
0505 nEventsWithCaloHits++;
0506 }
0507
0508
0509
0510
0511 hSamplingFractionEta->Fill(mceta,
0512 sumActiveCaloEnergy / (sumActiveCaloEnergy + sumPassiveCaloEnergy));
0513 std::sort(input_tower_rec.begin(), input_tower_rec.end(), &acompare);
0514 std::sort(input_tower_recSav.begin(), input_tower_recSav.end(), &acompare);
0515 std::sort(input_tower_sim.begin(), input_tower_sim.end(), &acompare);
0516
0517
0518
0519
0520
0521
0522 double tot_energyRecHit = 0;
0523 for (auto& tower : input_tower_rec) {
0524 tower.energy = tower.energy;
0525 tot_energyRecHit += tower.energy;
0526 }
0527
0528 double samplingFractionFe = 0.037;
0529 double samplingFractionW = 0.019;
0530 int minCellIDzDiffSamp = 5;
0531
0532 double tot_energySimHit = 0;
0533 for (auto& tower : input_tower_sim) {
0534 if (tower.cellIDz < minCellIDzDiffSamp) {
0535 tower.energy = tower.energy / samplingFractionW;
0536 } else {
0537 tower.energy = tower.energy / samplingFractionFe;
0538 }
0539 tot_energySimHit += tower.energy;
0540 }
0541 m_log->trace("Mc E: {} \t eta: {} \t sim E rec: {}\t rec E rec: {}", mcenergy, mceta,
0542 tot_energySimHit, tot_energyRecHit);
0543
0544
0545
0546
0547
0548 hClusterNCells_E_eta->Fill(mcenergy, nCaloHitsRec, mceta);
0549 hClusterEcalib_E_eta->Fill(mcenergy, tot_energyRecHit / mcenergy, mceta);
0550 hClusterEcalib_E_phi->Fill(mcenergy, tot_energyRecHit / mcenergy, mcphi);
0551
0552 hClusterSimNCells_E_eta->Fill(mcenergy, nCaloHitsSim, mceta);
0553 hClusterESimcalib_E_eta->Fill(mcenergy, tot_energySimHit / mcenergy, mceta);
0554 hClusterESimcalib_E_phi->Fill(mcenergy, tot_energySimHit / mcenergy, mcphi);
0555
0556
0557
0558
0559 int removedCells = 0;
0560 float minAggE = 0.001;
0561 float seedE = 0.100;
0562
0563 if (!input_tower_rec.empty()) {
0564
0565
0566 while (input_tower_rec.at(input_tower_rec.size() - 1).energy < minAggE) {
0567 input_tower_rec.pop_back();
0568 removedCells++;
0569 }
0570 m_log->trace("removed {} with E < {} GeV", removedCells, minAggE);
0571
0572 int nclusters = 0;
0573
0574 std::vector<clustersStrct> clusters_calo;
0575
0576 std::vector<towersStrct> cluster_towers;
0577 while (!input_tower_rec.empty()) {
0578 cluster_towers.clear();
0579 clustersStrct tempstructC;
0580
0581 if (input_tower_rec.at(0).energy > seedE) {
0582 m_log->trace("seed: {}\t {} \t {}", input_tower_rec.at(0).energy,
0583 input_tower_rec.at(0).cellIDx, input_tower_rec.at(0).cellIDy,
0584 input_tower_rec.at(0).cellIDz);
0585 tempstructC = findMACluster(seedE, minAggE, input_tower_rec, cluster_towers);
0586
0587
0588 float* showershape_eta_phi =
0589 CalculateM02andWeightedPosition(cluster_towers, tempstructC.cluster_E, 4.5);
0590
0591 tempstructC.cluster_M02 = showershape_eta_phi[0];
0592 tempstructC.cluster_M20 = showershape_eta_phi[1];
0593 tempstructC.cluster_Eta = showershape_eta_phi[2];
0594 tempstructC.cluster_Phi = showershape_eta_phi[3];
0595 tempstructC.cluster_X = showershape_eta_phi[4];
0596 tempstructC.cluster_Y = showershape_eta_phi[5];
0597 tempstructC.cluster_Z = showershape_eta_phi[6];
0598
0599 tempstructC.cluster_towers = cluster_towers;
0600 m_log->trace("---------> \t {} \tcluster with E = {} \tEta: {} \tPhi: {} \tX: {} \tY: {} "
0601 "\tZ: {} \tntowers: {} \ttrueID: {}",
0602 nclusters, tempstructC.cluster_E, tempstructC.cluster_Eta,
0603 tempstructC.cluster_Phi, tempstructC.cluster_X, tempstructC.cluster_Y,
0604 tempstructC.cluster_Z, tempstructC.cluster_NTowers,
0605 tempstructC.cluster_trueID);
0606 clusters_calo.push_back(tempstructC);
0607
0608 clusters_calo.push_back(tempstructC);
0609
0610 nclusters++;
0611 } else {
0612 m_log->trace("remaining: {} largest: {} \t {} \t {} \t {}", (int)input_tower_rec.size(),
0613 input_tower_rec.at(0).energy, input_tower_rec.at(0).cellIDx,
0614 input_tower_rec.at(0).cellIDy, input_tower_rec.at(0).cellIDz);
0615 input_tower_rec.clear();
0616 }
0617 }
0618
0619
0620
0621
0622 std::sort(clusters_calo.begin(), clusters_calo.end(), &acompareCl);
0623 m_log->info("-----> found {} clusters", clusters_calo.size());
0624 hRecNClusters_E_eta->Fill(mcenergy, clusters_calo.size(), mceta);
0625 int iCl = 0;
0626 for (const auto& cluster : clusters_calo) {
0627 if (iCl < maxNCluster && enableTreeCluster) {
0628 t_lFHCal_cluster_E[iCl] = (float)cluster.cluster_E;
0629 t_lFHCal_cluster_NCells[iCl] = (int)cluster.cluster_NTowers;
0630 t_lFHCal_cluster_Eta[iCl] = (float)cluster.cluster_Eta;
0631 t_lFHCal_cluster_Phi[iCl] = (float)cluster.cluster_Phi;
0632 }
0633 hRecClusterEcalib_E_eta->Fill(mcenergy, cluster.cluster_E / mcenergy, mceta);
0634 for (const auto& cluster_tower : cluster.cluster_towers) {
0635 int pSav = 0;
0636 while (cluster_tower.cellID != input_tower_recSav.at(pSav).cellID &&
0637 pSav < (int)input_tower_recSav.size()) {
0638 pSav++;
0639 }
0640 if (cluster_tower.cellID == input_tower_recSav.at(pSav).cellID) {
0641 input_tower_recSav.at(pSav).tower_clusterIDA = iCl;
0642 }
0643 }
0644
0645 if (iCl == 0) {
0646 hRecClusterEcalib_Ehigh_eta->Fill(mcenergy, cluster.cluster_E / mcenergy, mceta);
0647 hRecClusterNCells_Ehigh_eta->Fill(mcenergy, cluster.cluster_NTowers, mceta);
0648 }
0649 iCl++;
0650 m_log->trace("MA cluster {}:\t {} \t {}", iCl, cluster.cluster_E, cluster.cluster_NTowers);
0651 }
0652 if (iCl < maxNCluster && enableTreeCluster) {
0653 t_lFHCal_clusters_N = (int)iCl;
0654 }
0655
0656 clusters_calo.clear();
0657 } else {
0658 hRecNClusters_E_eta->Fill(mcenergy, 0., mceta);
0659 if (enableTreeCluster) {
0660 t_lFHCal_clusters_N = 0;
0661 }
0662 }
0663
0664
0665
0666
0667 int iClF = 0;
0668 float highestEFr = 0;
0669 int iClFHigh = 0;
0670
0671 const auto& lfhcalClustersF = *(event->GetCollection<edm4eic::Cluster>(nameClusters));
0672 for (const auto cluster : lfhcalClustersF) {
0673 if (cluster.getEnergy() > highestEFr) {
0674 iClFHigh = iClF;
0675 highestEFr = cluster.getEnergy();
0676 }
0677 hRecFClusterEcalib_E_eta->Fill(mcenergy, cluster.getEnergy() / mcenergy, mceta);
0678 m_log->trace("Island cluster {}:\t {} \t {}", iClF, cluster.getEnergy(), cluster.getNhits());
0679 iClF++;
0680 }
0681 hRecFNClusters_E_eta->Fill(mcenergy, iClF, mceta);
0682
0683 iClF = 0;
0684 for (const auto cluster : lfhcalClustersF) {
0685 if (iClF == iClFHigh) {
0686 hRecFClusterEcalib_Ehigh_eta->Fill(mcenergy, cluster.getEnergy() / mcenergy, mceta);
0687 hRecFClusterNCells_Ehigh_eta->Fill(mcenergy, cluster.getNhits(), mceta);
0688 }
0689 iClF++;
0690 }
0691
0692
0693
0694
0695 int iECl = 0;
0696 float highestEEmCl = 0;
0697 int iEClHigh = 0;
0698
0699 if (enableECalCluster) {
0700 try {
0701 const auto& fEMCClustersF = *(event->GetCollection<edm4eic::Cluster>("EcalEndcapPClusters"));
0702 m_log->info("-----> found fEMCClustersF:", fEMCClustersF.size());
0703 for (const auto cluster : fEMCClustersF) {
0704 if (iECl < maxNCluster && enableTreeCluster) {
0705 t_fEMC_cluster_E[iECl] = (float)cluster.getEnergy();
0706 t_fEMC_cluster_NCells[iECl] = (int)cluster.getNhits();
0707 t_fEMC_cluster_Eta[iECl] =
0708 (-1.) * std::log(std::tan((float)cluster.getIntrinsicTheta() / 2.));
0709 t_fEMC_cluster_Phi[iECl] = (float)cluster.getIntrinsicPhi();
0710 }
0711
0712 if (cluster.getEnergy() > highestEEmCl) {
0713 iEClHigh = iECl;
0714 highestEEmCl = cluster.getEnergy();
0715 }
0716 hRecFEmClusterEcalib_E_eta->Fill(mcenergy, cluster.getEnergy() / mcenergy, mceta);
0717 iECl++;
0718 }
0719 t_fEMC_clusters_N = iECl;
0720 hRecFEmNClusters_E_eta->Fill(mcenergy, iECl, mceta);
0721
0722
0723 iECl = 0;
0724 for (const auto cluster : fEMCClustersF) {
0725 if (iECl == iEClHigh) {
0726 hRecFEmClusterEcalib_Ehigh_eta->Fill(mcenergy, cluster.getEnergy() / mcenergy, mceta);
0727 }
0728 iECl++;
0729 }
0730 } catch (...) {
0731 std::cout << "ECal clusters not in output" << std::endl;
0732 enableECalCluster = false;
0733 }
0734 }
0735
0736
0737
0738 if (enableTree) {
0739 t_lFHCal_towers_N = (int)input_tower_recSav.size();
0740 for (int iCell = 0; iCell < (int)input_tower_recSav.size(); iCell++) {
0741 m_log->trace("{} \t {} \t {} \t {} \t {} \t {}", input_tower_recSav.at(iCell).cellIDx,
0742 input_tower_recSav.at(iCell).cellIDy, input_tower_recSav.at(iCell).cellIDz,
0743 input_tower_recSav.at(iCell).energy,
0744 input_tower_recSav.at(iCell).tower_clusterIDA,
0745 input_tower_recSav.at(iCell).tower_clusterIDB);
0746
0747 t_lFHCal_towers_cellE[iCell] = (float)input_tower_recSav.at(iCell).energy;
0748 t_lFHCal_towers_cellT[iCell] = (float)input_tower_recSav.at(iCell).time;
0749 t_lFHCal_towers_cellIDx[iCell] = (short)input_tower_recSav.at(iCell).cellIDx;
0750 t_lFHCal_towers_cellIDy[iCell] = (short)input_tower_recSav.at(iCell).cellIDy;
0751 t_lFHCal_towers_cellIDz[iCell] = (short)input_tower_recSav.at(iCell).cellIDz;
0752 t_lFHCal_towers_clusterIDA[iCell] = (short)input_tower_recSav.at(iCell).tower_clusterIDA;
0753 t_lFHCal_towers_clusterIDB[iCell] = (short)input_tower_recSav.at(iCell).tower_clusterIDB;
0754 t_lFHCal_towers_cellTrueID[iCell] = (int)input_tower_recSav.at(iCell).tower_trueID;
0755 }
0756
0757 event_tree->Fill();
0758
0759 t_lFHCal_towers_N = 0;
0760 for (Int_t itow = 0; itow < maxNTowers; itow++) {
0761 t_lFHCal_towers_cellE[itow] = 0;
0762 t_lFHCal_towers_cellT[itow] = 0;
0763 t_lFHCal_towers_cellIDx[itow] = 0;
0764 t_lFHCal_towers_cellIDy[itow] = 0;
0765 t_lFHCal_towers_cellIDz[itow] = 0;
0766 t_lFHCal_towers_clusterIDA[itow] = 0;
0767 t_lFHCal_towers_clusterIDB[itow] = 0;
0768 t_lFHCal_towers_cellTrueID[itow] = 0;
0769 }
0770 }
0771
0772
0773
0774
0775 if (enableTreeCluster) {
0776 cluster_tree->Fill();
0777
0778 t_mc_N = 0;
0779 t_lFHCal_clusters_N = 0;
0780 t_fEMC_clusters_N = 0;
0781 for (Int_t iMC = 0; iMC < maxNMC; iMC++) {
0782 t_mc_E[iMC] = 0;
0783 t_mc_Phi[iMC] = 0;
0784 t_mc_Eta[iMC] = 0;
0785 }
0786 for (Int_t iCl = 0; iCl < maxNCluster; iCl++) {
0787 t_lFHCal_cluster_E[iCl] = 0;
0788 t_lFHCal_cluster_NCells[iCl] = 0;
0789 t_lFHCal_cluster_Eta[iCl] = 0;
0790 t_lFHCal_cluster_Phi[iCl] = 0;
0791 t_fEMC_cluster_E[iCl] = 0;
0792 t_fEMC_cluster_NCells[iCl] = 0;
0793 t_fEMC_cluster_Eta[iCl] = 0;
0794 t_fEMC_cluster_Phi[iCl] = 0;
0795 }
0796 }
0797 }
0798
0799
0800
0801
0802 void lfhcal_studiesProcessor::Finish() {
0803 std::cout << "------> LFHCal " << nEventsWithCaloHits << " with calo info present" << std::endl;
0804
0805
0806 if (enableTree) {
0807 delete[] t_lFHCal_towers_cellE;
0808 delete[] t_lFHCal_towers_cellT;
0809 delete[] t_lFHCal_towers_cellIDx;
0810 delete[] t_lFHCal_towers_cellIDy;
0811 delete[] t_lFHCal_towers_cellIDz;
0812 delete[] t_lFHCal_towers_clusterIDA;
0813 delete[] t_lFHCal_towers_clusterIDB;
0814 delete[] t_lFHCal_towers_cellTrueID;
0815 }
0816
0817 if (enableTreeCluster) {
0818 delete[] t_mc_E;
0819 delete[] t_mc_Phi;
0820 delete[] t_mc_Eta;
0821 delete[] t_lFHCal_cluster_E;
0822 delete[] t_lFHCal_cluster_NCells;
0823 delete[] t_lFHCal_cluster_Phi;
0824 delete[] t_lFHCal_cluster_Eta;
0825 delete[] t_fEMC_cluster_E;
0826 delete[] t_fEMC_cluster_NCells;
0827 delete[] t_fEMC_cluster_Eta;
0828 delete[] t_fEMC_cluster_Phi;
0829 }
0830 }
