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