Warning, file /EICrecon/src/benchmarks/reconstruction/femc_studies/femc_studiesProcessor.cc was not indexed
or was modified since last indexation (in which case cross-reference links may be missing, inaccurate or erroneous).
0001
0002
0003
0004
0005
0006
0007 #include "femc_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 "benchmarks/reconstruction/lfhcal_studies/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 femc_studiesProcessor::Init() {
0047 std::string plugin_name = ("femc_studies");
0048
0049
0050
0051
0052 auto* app = GetApplication();
0053
0054 m_log = app->GetService<Log_service>()->logger(plugin_name);
0055
0056
0057 auto root_file_service = app->GetService<RootFile_service>();
0058
0059
0060 auto dd4hep_service = app->GetService<DD4hep_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 hClusterEcalib_E_eta->SetDirectory(m_dir_main);
0093 hClusterNCells_E_eta->SetDirectory(m_dir_main);
0094 hClusterEcalib_E_phi->SetDirectory(m_dir_main);
0095 hPosCaloHitsXY->SetDirectory(m_dir_main);
0096
0097
0098
0099
0100 hClusterESimcalib_E_eta =
0101 new TH3D("hClusterESimcalib_E_eta", "; E_{MC} (GeV); E_{rec,sim hit}/E_{MC}; #eta", 1500, 0.,
0102 150.0, 200, 0., 2.0, 50, 0, 5);
0103 hClusterSimNCells_E_eta =
0104 new TH3D("hClusterSimNCells_E_eta", "; E_{MC} (GeV); N_{cells, sim}; #eta", 1500, 0., 150.0,
0105 500, -0.5, 499.5, 50, 0, 5);
0106 hClusterESimcalib_E_phi =
0107 new TH3D("hClusterESimcalib_E_phi", "; E_{MC} (GeV); E_{rec,sim hit}/E_{MC}; #varphi (rad)",
0108 1500, 0., 150.0, 200, 0., 2.0, 360, -TMath::Pi(), TMath::Pi());
0109 hCellESim_layerX = new TH2D("hCellESim_layerX", "; #cell ID X; E_{rec,sim hit} (GeV)", 500, -0.5,
0110 499.5, 5000, 0, 1);
0111 hCellESim_layerY = new TH2D("hCellESim_layerY", "; #cell ID Y; E_{rec,sim hit} (GeV)", 500, -0.5,
0112 499.5, 5000, 0, 1);
0113 hCellTSim_layerX = new TH2D("hCellTSim_layerX", "; #cell ID X; t_{rec,sim hit} (GeV)", 500, -0.5,
0114 499.5, 5000, 0, 1000);
0115 hPosCaloSimHitsXY =
0116 new TH2D("hPosCaloSimHitsXY", "; X (cm); Y (cm)", 400, -400., 400., 400, -400., 400.);
0117 hClusterESimcalib_E_eta->SetDirectory(m_dir_main);
0118 hClusterSimNCells_E_eta->SetDirectory(m_dir_main);
0119 hClusterESimcalib_E_phi->SetDirectory(m_dir_main);
0120 hCellESim_layerX->SetDirectory(m_dir_main);
0121 hCellESim_layerY->SetDirectory(m_dir_main);
0122 hCellTSim_layerX->SetDirectory(m_dir_main);
0123 hPosCaloSimHitsXY->SetDirectory(m_dir_main);
0124
0125
0126
0127
0128 hRecClusterEcalib_E_eta =
0129 new TH3D("hRecClusterEcalib_E_eta", "; E_{MC} (GeV); E_{rec,rec clus}/E_{MC}; #eta", 1500, 0.,
0130 150.0, 200, 0., 2.0, 50, 0, 5);
0131 hRecNClusters_E_eta = new TH3D("hRecNClusters_E_eta", "; E_{MC} (GeV); N_{rec cl.}; #eta", 1500,
0132 0., 150.0, 10, -0.5, 9.5, 50, 0, 5);
0133
0134 hRecClusterEcalib_Ehigh_eta =
0135 new TH3D("hRecClusterEcalib_Ehigh_eta", "; E_{MC} (GeV); E_{rec,rec clus high.}/E_{MC}; #eta",
0136 1500, 0., 150.0, 200, 0., 2.0, 50, 0, 5);
0137 hRecClusterNCells_Ehigh_eta =
0138 new TH3D("hRecClusterNCells_Ehigh_eta", "; E_{MC} (GeV); N_{cells, rec cl., high.}; #eta",
0139 1500, 0., 150.0, 500, -0.5, 499.5, 50, 0, 5);
0140 hRecClusterEcalib_E_eta->SetDirectory(m_dir_main);
0141 hRecNClusters_E_eta->SetDirectory(m_dir_main);
0142 hRecClusterEcalib_Ehigh_eta->SetDirectory(m_dir_main);
0143 hRecClusterNCells_Ehigh_eta->SetDirectory(m_dir_main);
0144
0145
0146
0147
0148 hRecFClusterEcalib_E_eta =
0149 new TH3D("hRecFClusterEcalib_E_eta", "; E_{MC} (GeV); E_{rec,island clus}/E_{MC}; #eta", 1500,
0150 0., 150.0, 200, 0., 2.0, 50, 0, 5);
0151 hRecFNClusters_E_eta = new TH3D("hRecFNClusters_E_eta", "; E_{MC} (GeV); N_{rec f. cl.}; #eta",
0152 1500, 0., 150.0, 10, -0.5, 9.5, 50, 0, 5);
0153
0154 hRecFClusterEcalib_Ehigh_eta = new TH3D("hRecFClusterEcalib_Ehigh_eta",
0155 "; E_{MC} (GeV); E_{rec,island clus high.}/E_{MC}; #eta",
0156 1500, 0., 150.0, 200, 0., 2.0, 50, 0, 5);
0157 hRecFClusterNCells_Ehigh_eta =
0158 new TH3D("hRecFClusterNCells_Ehigh_eta", "; E_{MC} (GeV); N_{cells, rec f. cl., high.}; #eta",
0159 1500, 0., 150.0, 500, -0.5, 499.5, 50, 0, 5);
0160 hRecFClusterEcalib_E_eta->SetDirectory(m_dir_main);
0161 hRecFNClusters_E_eta->SetDirectory(m_dir_main);
0162 hRecFClusterEcalib_Ehigh_eta->SetDirectory(m_dir_main);
0163 hRecFClusterNCells_Ehigh_eta->SetDirectory(m_dir_main);
0164
0165
0166
0167
0168 hSamplingFractionEta = new TH2D("hSamplingFractionEta", "; #eta; f", 400, 1., 5., 500, 0., 0.2);
0169 hSamplingFractionEta->SetDirectory(m_dir_main);
0170
0171
0172
0173
0174 if (enableTree) {
0175 event_tree = new TTree("event_tree", "event_tree");
0176 event_tree->SetDirectory(m_dir_main);
0177
0178 t_fEMC_towers_cellE = new float[maxNTowers];
0179 t_fEMC_towers_cellT = new float[maxNTowers];
0180 t_fEMC_towers_cellIDx = new short[maxNTowers];
0181 t_fEMC_towers_cellIDy = new short[maxNTowers];
0182 t_fEMC_towers_clusterIDA = new short[maxNTowers];
0183 t_fEMC_towers_clusterIDB = new short[maxNTowers];
0184 t_fEMC_towers_cellTrueID = new int[maxNTowers];
0185
0186
0187 event_tree->Branch("tower_FEMC_N", &t_fEMC_towers_N, "tower_FEMC_N/I");
0188 event_tree->Branch("tower_FEMC_E", t_fEMC_towers_cellE, "tower_FEMC_E[tower_FEMC_N]/F");
0189 event_tree->Branch("tower_FEMC_T", t_fEMC_towers_cellT, "tower_FEMC_T[tower_FEMC_N]/F");
0190 event_tree->Branch("tower_FEMC_ix", t_fEMC_towers_cellIDx, "tower_FEMC_ix[tower_FEMC_N]/S");
0191 event_tree->Branch("tower_FEMC_iy", t_fEMC_towers_cellIDy, "tower_FEMC_iy[tower_FEMC_N]/S");
0192 event_tree->Branch("tower_FEMC_clusIDA", t_fEMC_towers_clusterIDA,
0193 "tower_FEMC_clusIDA[tower_FEMC_N]/S");
0194 event_tree->Branch("tower_FEMC_clusIDB", t_fEMC_towers_clusterIDB,
0195 "tower_FEMC_clusIDB[tower_FEMC_N]/S");
0196 event_tree->Branch("tower_FEMC_trueID", t_fEMC_towers_cellTrueID,
0197 "tower_FEMC_trueID[tower_FEMC_N]/I");
0198 }
0199
0200
0201
0202
0203 if (enableTreeCluster) {
0204 cluster_tree = new TTree("cluster_tree", "cluster_tree");
0205 cluster_tree->SetDirectory(m_dir_main);
0206
0207 t_mc_E = new float[maxNMC];
0208 t_mc_Phi = new float[maxNMC];
0209 t_mc_Eta = new float[maxNMC];
0210 t_fEMC_cluster_E = new float[maxNCluster];
0211 t_fEMC_cluster_NCells = new int[maxNCluster];
0212 t_fEMC_cluster_Phi = new float[maxNCluster];
0213 t_fEMC_cluster_Eta = new float[maxNCluster];
0214
0215
0216 cluster_tree->Branch("mc_N", &t_mc_N, "mc_N/I");
0217 cluster_tree->Branch("mc_E", t_mc_E, "mc_E[mc_N]/F");
0218 cluster_tree->Branch("mc_Phi", t_mc_Phi, "mc_Phi[mc_N]/F");
0219 cluster_tree->Branch("mc_Eta", t_mc_Eta, "mc_Eta[mc_N]/F");
0220
0221 cluster_tree->Branch("cluster_FEMC_N", &t_fEMC_clusters_N, "cluster_FEMC_N/I");
0222 cluster_tree->Branch("cluster_FEMC_E", t_fEMC_cluster_E, "cluster_FEMC_E[cluster_FEMC_N]/F");
0223 cluster_tree->Branch("cluster_FEMC_Ncells", t_fEMC_cluster_NCells,
0224 "cluster_FEMC_Ncells[cluster_FEMC_N]/I");
0225 cluster_tree->Branch("cluster_FEMC_Eta", t_fEMC_cluster_Eta,
0226 "cluster_FEMC_Eta[cluster_FEMC_N]/F");
0227 cluster_tree->Branch("cluster_FEMC_Phi", t_fEMC_cluster_Phi,
0228 "cluster_FEMC_Phi[cluster_FEMC_N]/F");
0229 }
0230
0231 std::cout << __PRETTY_FUNCTION__ << " " << __LINE__ << std::endl;
0232 auto detector = dd4hep_service->detector();
0233 std::cout << "--------------------------\nID specification:\n";
0234 try {
0235 m_decoder = detector->readout("EcalEndcapPHits").idSpec().decoder();
0236 std::cout << "1st: " << m_decoder << std::endl;
0237 std::cout << "full list: "
0238 << " " << m_decoder->fieldDescription() << std::endl;
0239 } catch (...) {
0240 std::cout << "2nd: " << m_decoder << std::endl;
0241 m_log->error("readoutClass not in the output");
0242 throw std::runtime_error("readoutClass not in the output.");
0243 }
0244 }
0245
0246
0247
0248
0249 void femc_studiesProcessor::Process(const std::shared_ptr<const JEvent>& event) {
0250
0251
0252
0253
0254
0255 const auto& mcParticles = *(event->GetCollection<edm4hep::MCParticle>("MCParticles"));
0256 double mceta = 0;
0257 double mcphi = 0;
0258 double mcp = 0;
0259 double mcenergy = 0;
0260 int iMC = 0;
0261 for (const auto mcparticle : mcParticles) {
0262 if (mcparticle.getGeneratorStatus() != 1) {
0263 continue;
0264 }
0265 const auto& mom = mcparticle.getMomentum();
0266
0267 mcenergy = mcparticle.getEnergy();
0268
0269 mceta = -log(tan(atan2(sqrt(mom.x * mom.x + mom.y * mom.y), mom.z) / 2.));
0270
0271 mcphi = atan2(mom.y, mom.x);
0272
0273 mcp = sqrt(mom.x * mom.x + mom.y * mom.y + mom.z * mom.z);
0274 m_log->trace("MC particle:{} \t {} \t {} \t totmom: {} phi {} eta {}", mom.x, mom.y, mom.z, mcp,
0275 mcphi, mceta);
0276 hMCEnergyVsEta->Fill(mcp, mceta);
0277
0278 if (enableTreeCluster) {
0279 if (iMC < maxNMC) {
0280 t_mc_E[iMC] = mcenergy;
0281 t_mc_Phi[iMC] = mcphi;
0282 t_mc_Eta[iMC] = mceta;
0283 }
0284 }
0285 iMC++;
0286 }
0287 if (enableTreeCluster) {
0288 t_mc_N = iMC;
0289 }
0290
0291
0292
0293 std::vector<towersStrct> input_tower_sim;
0294 int nCaloHitsSim = 0;
0295 float sumActiveCaloEnergy = 0;
0296 float sumPassiveCaloEnergy = 0;
0297 const auto& simHits = *(event->GetCollection<edm4hep::SimCalorimeterHit>(nameSimHits));
0298 for (const auto caloHit : simHits) {
0299 float x = caloHit.getPosition().x / 10.;
0300 float y = caloHit.getPosition().y / 10.;
0301 float z = caloHit.getPosition().z / 10.;
0302 uint64_t cellID = caloHit.getCellID();
0303 float energy = caloHit.getEnergy();
0304 double time = std::numeric_limits<double>::max();
0305 for (const auto& c : caloHit.getContributions()) {
0306 time = std::min<double>(c.getTime(), time);
0307 }
0308
0309 auto detector_layer_x = floor((x + 246) / 2.5);
0310 auto detector_layer_y = floor((y + 246) / 2.5);
0311 auto detector_passive = 0;
0312 if (detector_passive == 0) {
0313 sumActiveCaloEnergy += energy;
0314 } else {
0315 sumPassiveCaloEnergy += energy;
0316 }
0317
0318 if (detector_passive > 0) {
0319 continue;
0320 }
0321
0322 int cellIDx = detector_layer_x;
0323 int cellIDy = detector_layer_y;
0324 int cellIDz = 0;
0325 nCaloHitsSim++;
0326
0327 hPosCaloSimHitsXY->Fill(x, y);
0328 hCellESim_layerX->Fill(cellIDx, energy);
0329 hCellESim_layerY->Fill(cellIDy, energy);
0330 hCellTSim_layerX->Fill(cellIDx, time);
0331
0332
0333 bool found = false;
0334 for (auto& tower : input_tower_sim) {
0335 if (tower.cellID == static_cast<decltype(tower.cellID)>(cellID)) {
0336 tower.energy += energy;
0337 found = true;
0338 break;
0339 }
0340 }
0341 if (!found) {
0342 towersStrct tempstructT;
0343 tempstructT.energy = energy;
0344 tempstructT.time = time;
0345 tempstructT.posx = x;
0346 tempstructT.posy = y;
0347 tempstructT.posz = z;
0348 tempstructT.cellID = cellID;
0349 tempstructT.cellIDx = cellIDx;
0350 tempstructT.cellIDy = cellIDy;
0351 tempstructT.cellIDz = cellIDz;
0352 tempstructT.tower_trueID = 0;
0353 input_tower_sim.push_back(tempstructT);
0354 }
0355 }
0356
0357
0358
0359
0360 const auto& recHits = *(event->GetCollection<edm4eic::CalorimeterHit>(nameRecHits));
0361 int nCaloHitsRec = 0;
0362 std::vector<towersStrct> input_tower_rec;
0363 std::vector<towersStrct> input_tower_recSav;
0364
0365 for (const auto caloHit : recHits) {
0366 float x = caloHit.getPosition().x / 10.;
0367 float y = caloHit.getPosition().y / 10.;
0368 float z = caloHit.getPosition().z / 10.;
0369 uint64_t cellID = caloHit.getCellID();
0370 float energy = caloHit.getEnergy();
0371 float time = caloHit.getTime();
0372
0373 auto detector_passive = 0;
0374 auto detector_layer_x = floor((x + 246) / 2.5);
0375 auto detector_layer_y = floor((y + 246) / 2.5);
0376 if (detector_passive > 0) {
0377 continue;
0378 }
0379
0380
0381 int cellIDx = detector_layer_x;
0382 int cellIDy = detector_layer_y;
0383
0384 hPosCaloHitsXY->Fill(x, y);
0385 nCaloHitsRec++;
0386
0387
0388 bool found = false;
0389 for (auto& tower : input_tower_rec) {
0390 if (tower.cellID == static_cast<decltype(tower.cellID)>(cellID)) {
0391 tower.energy += energy;
0392 found = true;
0393 break;
0394 }
0395 }
0396 if (!found) {
0397 towersStrct tempstructT;
0398 tempstructT.energy = energy;
0399 tempstructT.time = time;
0400 tempstructT.posx = x;
0401 tempstructT.posy = y;
0402 tempstructT.posz = z;
0403 tempstructT.cellID = cellID;
0404 tempstructT.cellIDx = cellIDx;
0405 tempstructT.cellIDy = cellIDy;
0406 tempstructT.tower_trueID = 0;
0407 input_tower_rec.push_back(tempstructT);
0408 input_tower_recSav.push_back(tempstructT);
0409 }
0410 }
0411 m_log->trace("FEMC mod: nCaloHits sim {}\t rec {}", nCaloHitsSim, nCaloHitsRec);
0412 if (nCaloHitsRec > 0) {
0413 nEventsWithCaloHits++;
0414 }
0415
0416
0417
0418
0419 hSamplingFractionEta->Fill(mceta,
0420 sumActiveCaloEnergy / (sumActiveCaloEnergy + sumPassiveCaloEnergy));
0421 std::ranges::sort(input_tower_rec, &acompare);
0422 std::ranges::sort(input_tower_recSav, &acompare);
0423 std::ranges::sort(input_tower_sim, &acompare);
0424
0425
0426
0427
0428
0429
0430 double tot_energyRecHit = 0;
0431 for (auto& tower : input_tower_rec) {
0432 tower.energy = tower.energy;
0433 tot_energyRecHit += tower.energy;
0434 }
0435
0436 double samplingFraction = 1.;
0437
0438 double tot_energySimHit = 0;
0439 for (auto& tower : input_tower_sim) {
0440 tower.energy = tower.energy / samplingFraction;
0441 tot_energySimHit += tower.energy;
0442 }
0443 m_log->trace("Mc E: {} \t eta: {} \t sim E rec: {}\t rec E rec: {}", mcenergy, mceta,
0444 tot_energySimHit, tot_energyRecHit);
0445
0446
0447
0448
0449
0450 hClusterNCells_E_eta->Fill(mcenergy, nCaloHitsRec, mceta);
0451 hClusterEcalib_E_eta->Fill(mcenergy, tot_energyRecHit / mcenergy, mceta);
0452 hClusterEcalib_E_phi->Fill(mcenergy, tot_energyRecHit / mcenergy, mcphi);
0453
0454 hClusterSimNCells_E_eta->Fill(mcenergy, nCaloHitsSim, mceta);
0455 hClusterESimcalib_E_eta->Fill(mcenergy, tot_energySimHit / mcenergy, mceta);
0456 hClusterESimcalib_E_phi->Fill(mcenergy, tot_energySimHit / mcenergy, mcphi);
0457
0458
0459
0460
0461 int removedCells = 0;
0462 float minAggE = 0.001;
0463 float seedE = 0.20;
0464
0465 if (!input_tower_rec.empty()) {
0466
0467
0468 while (input_tower_rec.at(input_tower_rec.size() - 1).energy < minAggE) {
0469 input_tower_rec.pop_back();
0470 removedCells++;
0471 }
0472 m_log->trace("removed {} with E < {} GeV", removedCells, minAggE);
0473
0474 int nclusters = 0;
0475
0476 std::vector<clustersStrct> clusters_calo;
0477
0478 std::vector<towersStrct> cluster_towers;
0479 while (!input_tower_rec.empty()) {
0480 cluster_towers.clear();
0481 clustersStrct tempstructC;
0482
0483 if (input_tower_rec.at(0).energy > seedE) {
0484 m_log->trace("seed: {}\t {} \t {}", input_tower_rec.at(0).energy,
0485 input_tower_rec.at(0).cellIDx, input_tower_rec.at(0).cellIDy);
0486 tempstructC = findMACluster(seedE, minAggE, input_tower_rec, cluster_towers, 0.1);
0487
0488
0489 float* showershape_eta_phi =
0490 CalculateM02andWeightedPosition(cluster_towers, tempstructC.cluster_E, 4.5);
0491
0492 tempstructC.cluster_M02 = showershape_eta_phi[0];
0493 tempstructC.cluster_M20 = showershape_eta_phi[1];
0494 tempstructC.cluster_Eta = showershape_eta_phi[2];
0495 tempstructC.cluster_Phi = showershape_eta_phi[3];
0496 tempstructC.cluster_X = showershape_eta_phi[4];
0497 tempstructC.cluster_Y = showershape_eta_phi[5];
0498 tempstructC.cluster_Z = showershape_eta_phi[6];
0499
0500 tempstructC.cluster_towers = cluster_towers;
0501 m_log->trace("---------> \t {} \tcluster with E = {} \tEta: {} \tPhi: {} \tX: {} \tY: {} "
0502 "\tZ: {} \tntowers: {} \ttrueID: {}",
0503 nclusters, tempstructC.cluster_E, tempstructC.cluster_Eta,
0504 tempstructC.cluster_Phi, tempstructC.cluster_X, tempstructC.cluster_Y,
0505 tempstructC.cluster_Z, tempstructC.cluster_NTowers,
0506 tempstructC.cluster_trueID);
0507 clusters_calo.push_back(tempstructC);
0508
0509 nclusters++;
0510 } else {
0511 m_log->trace("remaining: {} largest: {} \t {} \t {} \t {}", (int)input_tower_rec.size(),
0512 input_tower_rec.at(0).energy, input_tower_rec.at(0).cellIDx,
0513 input_tower_rec.at(0).cellIDy, input_tower_rec.at(0).cellIDz);
0514 input_tower_rec.clear();
0515 }
0516 }
0517
0518
0519
0520
0521 std::ranges::sort(clusters_calo, &acompareCl);
0522 m_log->info("-----> found {} clusters", clusters_calo.size());
0523 hRecNClusters_E_eta->Fill(mcenergy, clusters_calo.size(), mceta);
0524 int iCl = 0;
0525 for (const auto& cluster : clusters_calo) {
0526 if (iCl < maxNCluster && enableTreeCluster) {
0527 t_fEMC_cluster_E[iCl] = (float)cluster.cluster_E;
0528 t_fEMC_cluster_NCells[iCl] = (int)cluster.cluster_NTowers;
0529 t_fEMC_cluster_Eta[iCl] = (float)cluster.cluster_Eta;
0530 t_fEMC_cluster_Phi[iCl] = (float)cluster.cluster_Phi;
0531 }
0532 hRecClusterEcalib_E_eta->Fill(mcenergy, cluster.cluster_E / mcenergy, mceta);
0533 for (const auto cluster_tower : cluster.cluster_towers) {
0534 int pSav = 0;
0535 while (cluster_tower.cellID != input_tower_recSav.at(pSav).cellID &&
0536 pSav < (int)input_tower_recSav.size()) {
0537 pSav++;
0538 }
0539 if (cluster_tower.cellID == input_tower_recSav.at(pSav).cellID) {
0540 input_tower_recSav.at(pSav).tower_clusterIDA = iCl;
0541 }
0542 }
0543
0544 if (iCl == 0) {
0545 hRecClusterEcalib_Ehigh_eta->Fill(mcenergy, cluster.cluster_E / mcenergy, mceta);
0546 hRecClusterNCells_Ehigh_eta->Fill(mcenergy, cluster.cluster_NTowers, mceta);
0547 }
0548 iCl++;
0549 m_log->trace("MA cluster {}:\t {} \t {}", iCl, cluster.cluster_E, cluster.cluster_NTowers);
0550 }
0551 if (iCl < maxNCluster && enableTreeCluster) {
0552 t_fEMC_clusters_N = iCl;
0553 }
0554
0555 clusters_calo.clear();
0556 } else {
0557 hRecNClusters_E_eta->Fill(mcenergy, 0., mceta);
0558 if (enableTreeCluster) {
0559 t_fEMC_clusters_N = 0;
0560 }
0561 }
0562
0563
0564
0565
0566 int iClF = 0;
0567 float highestEFr = 0;
0568 int iClFHigh = 0;
0569
0570 const auto& fecalClustersF = *(event->GetCollection<edm4eic::Cluster>(nameClusters));
0571 for (const auto cluster : fecalClustersF) {
0572 if (cluster.getEnergy() > highestEFr) {
0573 iClFHigh = iClF;
0574 highestEFr = cluster.getEnergy();
0575 }
0576 hRecFClusterEcalib_E_eta->Fill(mcenergy, cluster.getEnergy() / mcenergy, mceta);
0577 m_log->trace("Island cluster {}:\t {} \t {}", iClF, cluster.getEnergy(), cluster.getNhits());
0578
0579 iClF++;
0580 }
0581 hRecFNClusters_E_eta->Fill(mcenergy, iClF, mceta);
0582
0583 iClF = 0;
0584 for (const auto cluster : fecalClustersF) {
0585 if (iClF == iClFHigh) {
0586 hRecFClusterEcalib_Ehigh_eta->Fill(mcenergy, cluster.getEnergy() / mcenergy, mceta);
0587 hRecFClusterNCells_Ehigh_eta->Fill(mcenergy, cluster.getNhits(), mceta);
0588 }
0589 iClF++;
0590 }
0591
0592
0593
0594
0595 if (enableTree) {
0596 t_fEMC_towers_N = (int)input_tower_recSav.size();
0597 for (int iCell = 0; iCell < (int)input_tower_recSav.size(); iCell++) {
0598 m_log->trace("{} \t {} \t {} \t {} \t {}", input_tower_recSav.at(iCell).cellIDx,
0599 input_tower_recSav.at(iCell).cellIDy, input_tower_recSav.at(iCell).energy,
0600 input_tower_recSav.at(iCell).tower_clusterIDA,
0601 input_tower_recSav.at(iCell).tower_clusterIDB);
0602
0603 t_fEMC_towers_cellE[iCell] = input_tower_recSav.at(iCell).energy;
0604 t_fEMC_towers_cellT[iCell] = input_tower_recSav.at(iCell).time;
0605 t_fEMC_towers_cellIDx[iCell] = (short)input_tower_recSav.at(iCell).cellIDx;
0606 t_fEMC_towers_cellIDy[iCell] = (short)input_tower_recSav.at(iCell).cellIDy;
0607 t_fEMC_towers_clusterIDA[iCell] = (short)input_tower_recSav.at(iCell).tower_clusterIDA;
0608 t_fEMC_towers_clusterIDB[iCell] = (short)input_tower_recSav.at(iCell).tower_clusterIDB;
0609 t_fEMC_towers_cellTrueID[iCell] = input_tower_recSav.at(iCell).tower_trueID;
0610 }
0611
0612 event_tree->Fill();
0613
0614 t_fEMC_towers_N = 0;
0615 for (Int_t itow = 0; itow < maxNTowers; itow++) {
0616 t_fEMC_towers_cellE[itow] = 0;
0617 t_fEMC_towers_cellT[itow] = 0;
0618 t_fEMC_towers_cellIDx[itow] = 0;
0619 t_fEMC_towers_cellIDy[itow] = 0;
0620 t_fEMC_towers_clusterIDA[itow] = 0;
0621 t_fEMC_towers_clusterIDB[itow] = 0;
0622 t_fEMC_towers_cellTrueID[itow] = 0;
0623 }
0624 }
0625
0626
0627
0628
0629 if (enableTreeCluster) {
0630 cluster_tree->Fill();
0631
0632 t_mc_N = 0;
0633 t_fEMC_clusters_N = 0;
0634 t_fEMC_clusters_N = 0;
0635 for (Int_t iMC = 0; iMC < maxNMC; iMC++) {
0636 t_mc_E[iMC] = 0;
0637 t_mc_Phi[iMC] = 0;
0638 t_mc_Eta[iMC] = 0;
0639 }
0640 for (Int_t iCl = 0; iCl < maxNCluster; iCl++) {
0641 t_fEMC_cluster_E[iCl] = 0;
0642 t_fEMC_cluster_NCells[iCl] = 0;
0643 t_fEMC_cluster_Eta[iCl] = 0;
0644 t_fEMC_cluster_Phi[iCl] = 0;
0645 }
0646 }
0647 }
0648
0649
0650
0651
0652 void femc_studiesProcessor::Finish() {
0653 std::cout << "------> FEMC " << nEventsWithCaloHits << " with calo info present" << std::endl;
0654 if (enableTreeCluster) {
0655 cluster_tree->Write();
0656 }
0657
0658
0659 if (enableTree) {
0660 delete[] t_fEMC_towers_cellE;
0661 delete[] t_fEMC_towers_cellT;
0662 delete[] t_fEMC_towers_cellIDx;
0663 delete[] t_fEMC_towers_cellIDy;
0664 delete[] t_fEMC_towers_clusterIDA;
0665 delete[] t_fEMC_towers_clusterIDB;
0666 delete[] t_fEMC_towers_cellTrueID;
0667 }
0668
0669 if (enableTreeCluster) {
0670 delete[] t_mc_E;
0671 delete[] t_mc_Phi;
0672 delete[] t_mc_Eta;
0673 delete[] t_fEMC_cluster_E;
0674 delete[] t_fEMC_cluster_NCells;
0675 delete[] t_fEMC_cluster_Phi;
0676 delete[] t_fEMC_cluster_Eta;
0677 }
0678 }
