Back to home page

EIC code displayed by LXR
 
 

    

Warning, file /EICrecon/src/benchmarks/reconstruction/lfhcal_studies/lfhcal_studiesProcessor.cc was not indexed or was modified since last indexation (in which case cross-reference links may be missing, inaccurate or erroneous).

0001 // Copyright 2023, Friederike Bock
0002 // Subject to the terms in the LICENSE file found in the top-level directory.
0003 //
0004 //  Sections Copyright (C) 2023 Friederike Bock
0005 //  under SPDX-License-Identifier: LGPL-3.0-or-later
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 <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 "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 // InitWithGlobalRootLock
0043 //******************************************************************************************//
0044 void lfhcal_studiesProcessor::Init() {
0045   std::string plugin_name = ("lfhcal_studies");
0046 
0047   // ===============================================================================================
0048   // Get JANA application and seup general variables
0049   // ===============================================================================================
0050   auto app = GetApplication();
0051 
0052   m_log = app->GetService<Log_service>()->logger(plugin_name);
0053 
0054   // Ask service locator for the DD4hep geometry
0055   auto dd4hep_service = app->GetService<DD4hep_service>();
0056 
0057   // Ask service locator a file to write histograms to
0058   auto root_file_service = app->GetService<RootFile_service>();
0059 
0060   // Get TDirectory for histograms root file
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   // Create a directory for this plugin. And subdirectories for series of histograms
0068   // ===============================================================================================
0069   m_dir_main = file->mkdir(plugin_name.c_str());
0070 
0071   // ===============================================================================================
0072   // Simulations hists
0073   // ===============================================================================================
0074   hMCEnergyVsEta        = new TH2D("hMCEnergyVsEta", "; E (GeV); #eta", 1500, 0., 150., 500, 0, 5);
0075   hMCEnergyVsEta->SetDirectory(m_dir_main);
0076 
0077   // ===============================================================================================
0078   // Sum cell clusters rec histos
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   hPosCaloHitsZX        = new TH2D("hPosCaloHitsZX", "; Z (cm); X (cm)", 200, 300., 500., 400, -400., 400.);
0088   hPosCaloHitsZY        = new TH2D("hPosCaloHitsZY", "; Z (cm); Y (cm)", 200, 300., 500., 400, -400., 400.);
0089   hClusterEcalib_E_eta->SetDirectory(m_dir_main);
0090   hClusterNCells_E_eta->SetDirectory(m_dir_main);
0091   hClusterEcalib_E_phi->SetDirectory(m_dir_main);
0092   hPosCaloHitsXY->SetDirectory(m_dir_main);
0093   hPosCaloHitsZX->SetDirectory(m_dir_main);
0094   hPosCaloHitsZY->SetDirectory(m_dir_main);
0095 
0096   // ===============================================================================================
0097   // Sum cell clusters sim histos
0098   // ===============================================================================================
0099   hClusterESimcalib_E_eta = new TH3D("hClusterESimcalib_E_eta", "; E_{MC} (GeV); E_{rec,sim hit}/E_{MC}; #eta" ,
0100                                       1500, 0., 150.0, 200, 0., 2.0, 50, 0, 5);
0101   hClusterSimNCells_E_eta = new TH3D("hClusterSimNCells_E_eta", "; E_{MC} (GeV); N_{cells, sim}; #eta",
0102                                       1500, 0., 150.0, 500, -0.5, 499.5, 50, 0, 5);
0103   hClusterESimcalib_E_phi = new TH3D("hClusterESimcalib_E_phi", "; E_{MC} (GeV); E_{rec,sim hit}/E_{MC}; #varphi (rad)" ,
0104                                       1500, 0., 150.0, 200, 0., 2.0, 360 , -TMath::Pi(), TMath::Pi());
0105   hCellESim_layerX        = new TH2D("hCellESim_layerX", "; #cell ID X; E_{rec,sim hit} (GeV)" , 240, -0.5, 239.5, 5000, 0, 1);
0106   hCellESim_layerY        = new TH2D("hCellESim_layerY", "; #cell ID Y; E_{rec,sim hit} (GeV)" , 240, -0.5, 239.5, 5000, 0, 1);
0107   hCellESim_layerZ        = new TH2D("hCellESim_layerZ", "; #cell ID Z; E_{rec,sim hit} (GeV)" , 70, -0.5, 69.5, 5000, 0, 1);
0108   hCellTSim_layerZ        = new TH2D("hCellTSim_layerZ", "; #cell ID Z; t_{rec,sim hit} (GeV)" , 70, -0.5, 69.5, 5000, 0, 1000);
0109   hPosCaloSimHitsXY       = new TH2D("hPosCaloSimHitsXY", "; X (cm); Y (cm)", 400, -400., 400., 400, -400., 400.);
0110   hPosCaloSimHitsZX       = new TH2D("hPosCaloSimHitsZX", "; Z (cm); X (cm)", 200, 300., 500., 400, -400., 400.);
0111   hPosCaloSimHitsZY       = new TH2D("hPosCaloSimHitsZY", "; Z (cm); Y (cm)", 200, 300., 500., 400, -400., 400.);
0112   hClusterESimcalib_E_eta->SetDirectory(m_dir_main);
0113   hClusterSimNCells_E_eta->SetDirectory(m_dir_main);
0114   hClusterESimcalib_E_phi->SetDirectory(m_dir_main);
0115   hCellESim_layerX->SetDirectory(m_dir_main);
0116   hCellESim_layerY->SetDirectory(m_dir_main);
0117   hCellESim_layerZ->SetDirectory(m_dir_main);
0118   hCellTSim_layerZ->SetDirectory(m_dir_main);
0119   hPosCaloSimHitsXY->SetDirectory(m_dir_main);
0120   hPosCaloSimHitsZX->SetDirectory(m_dir_main);
0121   hPosCaloSimHitsZY->SetDirectory(m_dir_main);
0122 
0123   // ===============================================================================================
0124   // rec cluster MA clusters histos
0125   // ===============================================================================================
0126   hRecClusterEcalib_E_eta     = new TH3D("hRecClusterEcalib_E_eta", "; E_{MC} (GeV); E_{rec,rec clus}/E_{MC}; #eta",
0127                                           1500, 0., 150.0, 200, 0., 2.0, 50, 0, 5);
0128   hRecNClusters_E_eta         = new TH3D("hRecNClusters_E_eta", "; E_{MC} (GeV); N_{rec cl.}; #eta",
0129                                           1500, 0., 150.0, 10, -0.5, 9.5, 50, 0, 5);
0130   // rec cluster highest
0131   hRecClusterEcalib_Ehigh_eta = new TH3D("hRecClusterEcalib_Ehigh_eta", "; E_{MC} (GeV); E_{rec,rec clus high.}/E_{MC}; #eta",
0132                                           1500, 0., 150.0, 200, 0., 2.0, 50, 0, 5);
0133   hRecClusterNCells_Ehigh_eta = new TH3D("hRecClusterNCells_Ehigh_eta", "; E_{MC} (GeV); N_{cells, rec cl., high.}; #eta",
0134                                           1500, 0., 150.0, 500, -0.5, 499.5, 50, 0, 5);
0135   hRecClusterEcalib_E_eta->SetDirectory(m_dir_main);
0136   hRecNClusters_E_eta->SetDirectory(m_dir_main);
0137   hRecClusterEcalib_Ehigh_eta->SetDirectory(m_dir_main);
0138   hRecClusterNCells_Ehigh_eta->SetDirectory(m_dir_main);
0139 
0140   // ===============================================================================================
0141   // rec cluster framework Island clusters histos
0142   // ===============================================================================================
0143   hRecFClusterEcalib_E_eta      = new TH3D("hRecFClusterEcalib_E_eta", "; E_{MC} (GeV); E_{rec,island clus}/E_{MC}; #eta",
0144                                             1500, 0., 150.0, 200, 0., 2.0, 50, 0, 5);
0145   hRecFNClusters_E_eta          = new TH3D("hRecFNClusters_E_eta", "; E_{MC} (GeV); N_{rec f. cl.}; #eta",
0146                                             1500, 0., 150.0, 10, -0.5, 9.5, 50, 0, 5);
0147   // rec cluster framework highest
0148   hRecFClusterEcalib_Ehigh_eta  = new TH3D("hRecFClusterEcalib_Ehigh_eta", "; E_{MC} (GeV); E_{rec,island clus high.}/E_{MC}; #eta",
0149                                             1500, 0., 150.0, 200, 0., 2.0, 50, 0, 5);
0150   hRecFClusterNCells_Ehigh_eta  = new TH3D("hRecFClusterNCells_Ehigh_eta", "; E_{MC} (GeV); N_{cells, rec f. cl., high.}; #eta",
0151                                             1500, 0., 150.0, 500, -0.5, 499.5, 50, 0, 5);
0152   hRecFClusterEcalib_E_eta->SetDirectory(m_dir_main);
0153   hRecFNClusters_E_eta->SetDirectory(m_dir_main);
0154   hRecFClusterEcalib_Ehigh_eta->SetDirectory(m_dir_main);
0155   hRecFClusterNCells_Ehigh_eta->SetDirectory(m_dir_main);
0156 
0157   // ===============================================================================================
0158   // FEcal rec cluster framework Island clusters histos
0159   // ===============================================================================================
0160   hRecFEmClusterEcalib_E_eta      = new TH3D("hRecFEmClusterEcalib_E_eta", "; E_{MC} (GeV); E_{Ecal, rec,island clus}/E_{MC}; #eta",
0161                                               1500, 0., 150.0, 200, 0., 2.0, 50, 0, 5);
0162   hRecFEmNClusters_E_eta          = new TH3D("hRecFEmNClusters_E_eta", "; E_{MC} (GeV); N_{Ecal, rec f. cl.}; #eta",
0163                                               1500, 0., 150.0, 10, -0.5, 9.5, 50, 0, 5);
0164   // rec cluster framework highest
0165   hRecFEmClusterEcalib_Ehigh_eta  = new TH3D("hRecFEmClusterEcalib_Ehigh_eta", "; E_{MC} (GeV); E_{Ecal, rec,island clus high.}/E_{MC}; #eta",
0166                                               1500, 0., 150.0, 200, 0., 2.0, 50, 0, 5);
0167   hRecFEmClusterEcalib_E_eta->SetDirectory(m_dir_main);
0168   hRecFEmNClusters_E_eta->SetDirectory(m_dir_main);
0169   hRecFEmClusterEcalib_Ehigh_eta->SetDirectory(m_dir_main);
0170 
0171   // ===============================================================================================
0172   // Sampling fraction
0173   // ===============================================================================================
0174   hSamplingFractionEta    = new TH2D("hSamplingFractionEta", "; #eta; f", 400, 1., 5., 500, 0., 0.2);
0175   hSamplingFractionEta->SetDirectory(m_dir_main);
0176 
0177   // ===============================================================================================
0178   // Tree for clusterizer studies
0179   // ===============================================================================================
0180   if (enableTree){
0181     event_tree                = new TTree("event_tree", "event_tree");
0182     event_tree->SetDirectory(m_dir_main);
0183 
0184     t_lFHCal_towers_cellE       = new float[maxNTowers];
0185     t_lFHCal_towers_cellT       = new float[maxNTowers];
0186     t_lFHCal_towers_cellIDx     = new short[maxNTowers];
0187     t_lFHCal_towers_cellIDy     = new short[maxNTowers];
0188     t_lFHCal_towers_cellIDz     = new short[maxNTowers];
0189     t_lFHCal_towers_clusterIDA  = new short[maxNTowers];
0190     t_lFHCal_towers_clusterIDB  = new short[maxNTowers];
0191     t_lFHCal_towers_cellTrueID  = new int[maxNTowers];
0192 
0193     // towers LFHCAL
0194     event_tree->Branch("tower_LFHCAL_N", &t_lFHCal_towers_N, "tower_LFHCAL_N/I");
0195     event_tree->Branch("tower_LFHCAL_E", t_lFHCal_towers_cellE, "tower_LFHCAL_E[tower_LFHCAL_N]/F");
0196     event_tree->Branch("tower_LFHCAL_T", t_lFHCal_towers_cellT, "tower_LFHCAL_T[tower_LFHCAL_N]/F");
0197     event_tree->Branch("tower_LFHCAL_ix", t_lFHCal_towers_cellIDx, "tower_LFHCAL_ix[tower_LFHCAL_N]/S");
0198     event_tree->Branch("tower_LFHCAL_iy", t_lFHCal_towers_cellIDy, "tower_LFHCAL_iy[tower_LFHCAL_N]/S");
0199     event_tree->Branch("tower_LFHCAL_iz", t_lFHCal_towers_cellIDz, "tower_LFHCAL_iz[tower_LFHCAL_N]/S");
0200     event_tree->Branch("tower_LFHCAL_clusIDA", t_lFHCal_towers_clusterIDA, "tower_LFHCAL_clusIDA[tower_LFHCAL_N]/S");
0201     event_tree->Branch("tower_LFHCAL_clusIDB", t_lFHCal_towers_clusterIDB, "tower_LFHCAL_clusIDB[tower_LFHCAL_N]/S");
0202     event_tree->Branch("tower_LFHCAL_trueID", t_lFHCal_towers_cellTrueID, "tower_LFHCAL_trueID[tower_LFHCAL_N]/I");
0203   }
0204 
0205   // ===============================================================================================
0206   // Tree for cluster studies
0207   // ===============================================================================================
0208   if (enableTreeCluster){
0209     cluster_tree          = new TTree("cluster_tree", "cluster_tree");
0210     cluster_tree->SetDirectory(m_dir_main);
0211 
0212     t_mc_E                  = new float[maxNMC];
0213     t_mc_Phi                = new float[maxNMC];
0214     t_mc_Eta                = new float[maxNMC];
0215     t_lFHCal_cluster_E      = new float[maxNCluster];
0216     t_lFHCal_cluster_NCells = new int[maxNCluster];
0217     t_lFHCal_cluster_Phi    = new float[maxNCluster];
0218     t_lFHCal_cluster_Eta    = new float[maxNCluster];
0219     t_fEMC_cluster_E       = new float[maxNCluster];
0220     t_fEMC_cluster_NCells  = new int[maxNCluster];
0221     t_fEMC_cluster_Eta     = new float[maxNCluster];
0222     t_fEMC_cluster_Phi     = new float[maxNCluster];
0223 
0224     // MC particles
0225     cluster_tree->Branch("mc_N", &t_mc_N, "mc_N/I");
0226     cluster_tree->Branch("mc_E", t_mc_E, "mc_E[mc_N]/F");
0227     cluster_tree->Branch("mc_Phi", t_mc_Phi, "mc_Phi[mc_N]/F");
0228     cluster_tree->Branch("mc_Eta", t_mc_Eta, "mc_Eta[mc_N]/F");
0229     // clusters LFHCAL
0230     cluster_tree->Branch("cluster_LFHCAL_N", &t_lFHCal_clusters_N, "cluster_LFHCAL_N/I");
0231     cluster_tree->Branch("cluster_LFHCAL_E", t_lFHCal_cluster_E, "cluster_LFHCAL_E[cluster_LFHCAL_N]/F");
0232     cluster_tree->Branch("cluster_LFHCAL_Ncells", t_lFHCal_cluster_NCells, "cluster_LFHCAL_Ncells[cluster_LFHCAL_N]/I");
0233     cluster_tree->Branch("cluster_LFHCAL_Eta", t_lFHCal_cluster_Eta, "cluster_LFHCAL_Eta[cluster_LFHCAL_N]/F");
0234     cluster_tree->Branch("cluster_LFHCAL_Phi", t_lFHCal_cluster_Phi, "cluster_LFHCAL_Phi[cluster_LFHCAL_N]/F");
0235     // clusters FECal
0236     cluster_tree->Branch("cluster_FECAL_N", &t_fEMC_clusters_N, "cluster_FECAL_N/I");
0237     cluster_tree->Branch("cluster_FECAL_E", t_fEMC_cluster_E, "cluster_FECAL_E[cluster_FECAL_N]/F");
0238     cluster_tree->Branch("cluster_FECAL_Ncells", t_fEMC_cluster_NCells, "cluster_FECAL_Ncells[cluster_FECAL_N]/I");
0239     cluster_tree->Branch("cluster_FECAL_Eta", t_fEMC_cluster_Eta, "cluster_FECAL_Eta[cluster_FECAL_N]/F");
0240     cluster_tree->Branch("cluster_FECAL_Phi", t_fEMC_cluster_Phi, "cluster_FECAL_Phi[cluster_FECAL_N]/F");
0241   }
0242 
0243   std::cout << __PRETTY_FUNCTION__ << " " << __LINE__ << std::endl;
0244   auto detector = dd4hep_service->detector();
0245   std::cout << "--------------------------\nID specification:\n";
0246   try {
0247     m_decoder = detector->readout("LFHCALHits").idSpec().decoder();
0248     std::cout <<"1st: "<< m_decoder << std::endl;
0249     iLx = m_decoder->index("towerx");
0250     iLy = m_decoder->index("towery");
0251     iLz = m_decoder->index("layerz");
0252     iPassive = m_decoder->index("passive");
0253     std::cout << "full list: " << " " << m_decoder->fieldDescription() << std::endl;
0254   } catch (...) {
0255       std::cout <<"2nd: "  << m_decoder << std::endl;
0256       m_log->error("readoutClass not in the output");
0257       throw std::runtime_error("readoutClass not in the output.");
0258   }
0259 
0260 }
0261 
0262 //******************************************************************************************//
0263 // ProcessSequential
0264 //******************************************************************************************//
0265 void lfhcal_studiesProcessor::Process(const std::shared_ptr<const JEvent>& event) {
0266 // void lfhcal_studiesProcessor::ProcessSequential(const std::shared_ptr<const JEvent>& event) {
0267 
0268   // ===============================================================================================
0269   // process MC particles
0270   // ===============================================================================================
0271   const auto &mcParticles = *(event->GetCollection<edm4hep::MCParticle>("MCParticles"));
0272   double mceta    = 0;
0273   double mcphi    = 0;
0274   double mcp      = 0;
0275   double mcenergy = 0;
0276   int iMC         = 0;
0277   for (auto mcparticle : mcParticles) {
0278     if (mcparticle.getGeneratorStatus() != 1)
0279       continue;
0280     const auto& mom = mcparticle.getMomentum();
0281     // get particle energy
0282     mcenergy = mcparticle.getEnergy();
0283     //determine mceta from momentum
0284     mceta = -log(tan(atan2(sqrt(mom.x * mom.x + mom.y * mom.y), mom.z) / 2.));
0285     // determine mcphi from momentum
0286     mcphi = atan2(mom.y, mom.x);
0287     // determine mc momentum
0288     mcp = sqrt(mom.x * mom.x + mom.y * mom.y + mom.z * mom.z);
0289     m_log->trace("MC particle:{} \t {} \t {} \t totmom: {} phi {} eta {}", mom.x, mom.y, mom.z, mcp, mcphi, mceta);
0290     hMCEnergyVsEta->Fill(mcp,mceta);
0291 
0292     if (enableTreeCluster){
0293       if (iMC < maxNMC){
0294         t_mc_E[iMC]   = mcenergy;
0295         t_mc_Phi[iMC] = mcphi;
0296         t_mc_Eta[iMC] = mceta;
0297       }
0298     }
0299     iMC++;
0300   }
0301   if (enableTreeCluster) t_mc_N = iMC;
0302   // ===============================================================================================
0303   // process sim hits
0304   // ===============================================================================================
0305   std::vector<towersStrct> input_tower_sim;
0306   int nCaloHitsSim = 0;
0307   float sumActiveCaloEnergy = 0;
0308   float sumPassiveCaloEnergy = 0;
0309   const auto &simHits = *(event->GetCollection<edm4hep::SimCalorimeterHit>(nameSimHits));
0310   for (const auto caloHit : simHits) {
0311     float x         = caloHit.getPosition().x / 10.;
0312     float y         = caloHit.getPosition().y / 10.;
0313     float z         = caloHit.getPosition().z / 10.;
0314     uint64_t cellID = caloHit.getCellID();
0315     float energy    = caloHit.getEnergy();
0316     double time = std::numeric_limits<double>::max();
0317     for (const auto& c : caloHit.getContributions()) {
0318         if (c.getTime() <= time) {
0319             time = c.getTime();
0320         }
0321     }
0322 
0323     auto detector_module_x  = m_decoder->get(cellID, 1);
0324     auto detector_module_y  = m_decoder->get(cellID, 2);
0325     auto detector_passive = m_decoder->get(cellID, iPassive);
0326     auto detector_layer_x = m_decoder->get(cellID, iLx);
0327     auto detector_layer_y = m_decoder->get(cellID, iLy);
0328     long detector_layer_rz = -1;
0329     if (isLFHCal)
0330       detector_layer_rz = m_decoder->get(cellID, 7);
0331     auto detector_layer_z = m_decoder->get(cellID, iLz);
0332     if(detector_passive == 0) {
0333       sumActiveCaloEnergy += energy;
0334     } else {
0335       sumPassiveCaloEnergy += energy;
0336     }
0337 
0338     if (detector_passive > 0) continue;
0339     // calc cell IDs
0340     long cellIDx = -1;
0341     long cellIDy = -1;
0342     long cellIDz = -1;
0343     if (isLFHCal){
0344       cellIDx = 54ll*2 - detector_module_x * 2 + detector_layer_x;
0345       cellIDy = 54ll*2 - detector_module_y * 2 + detector_layer_y;
0346       cellIDz = detector_layer_rz*10+detector_layer_z;
0347     }
0348     nCaloHitsSim++;
0349 
0350     hPosCaloSimHitsXY->Fill(x, y);
0351     hPosCaloSimHitsZX->Fill(z, x);
0352     hPosCaloSimHitsZY->Fill(z, y);
0353 
0354     hCellESim_layerZ->Fill(cellIDz, energy);
0355     hCellESim_layerX->Fill(cellIDx, energy);
0356     hCellESim_layerY->Fill(cellIDy, energy);
0357     hCellTSim_layerZ->Fill(cellIDz, time);
0358 
0359     //loop over input_tower_sim and find if there is already a tower with the same cellID
0360     bool found = false;
0361     for (auto& tower : input_tower_sim) {
0362       if (tower.cellID == cellID) {
0363         tower.energy += energy;
0364         found = true;
0365         break;
0366       }
0367     }
0368     if (!found) {
0369       towersStrct tempstructT;
0370       tempstructT.energy        = energy;
0371       tempstructT.time          = time;
0372       tempstructT.posx          = x;
0373       tempstructT.posy          = y;
0374       tempstructT.posz          = z;
0375       tempstructT.cellID        = cellID;
0376       tempstructT.cellIDx       = cellIDx;
0377       tempstructT.cellIDy       = cellIDy;
0378       tempstructT.cellIDz       = cellIDz;
0379       tempstructT.tower_trueID  = 0; //TODO how to get trueID?
0380       input_tower_sim.push_back(tempstructT);
0381     }
0382   }
0383 
0384 
0385   // ===============================================================================================
0386   // read rec hits & fill structs
0387   // ===============================================================================================
0388   const auto &recHits = *(event->GetCollection<edm4eic::CalorimeterHit>(nameRecHits));
0389   int nCaloHitsRec = 0;
0390   std::vector<towersStrct> input_tower_rec;
0391   std::vector<towersStrct> input_tower_recSav;
0392   // process rec hits
0393   for (const auto caloHit : recHits) {
0394     float x         = caloHit.getPosition().x / 10.;
0395     float y         = caloHit.getPosition().y / 10.;
0396     float z         = caloHit.getPosition().z / 10.;
0397     uint64_t cellID = caloHit.getCellID();
0398     float energy    = caloHit.getEnergy();
0399     float time      = caloHit.getTime();
0400 
0401     auto detector_module_x  = m_decoder->get(cellID, 1);
0402     auto detector_module_y  = m_decoder->get(cellID, 2);
0403     auto detector_passive = m_decoder->get(cellID, iPassive);
0404     auto detector_layer_x = m_decoder->get(cellID, iLx);
0405     auto detector_layer_y = m_decoder->get(cellID, iLy);
0406     int detector_layer_rz = -1;
0407     if (isLFHCal){
0408       detector_layer_rz   = m_decoder->get(cellID, 7);
0409     }
0410 
0411     if (detector_passive > 0) continue;
0412 
0413     // calc cell IDs
0414     long cellIDx = -1;
0415     long cellIDy = -1;
0416     if (isLFHCal){
0417       cellIDx = 54ll*2 - detector_module_x * 2 + detector_layer_x;
0418       cellIDy = 54ll*2 - detector_module_y * 2 + detector_layer_y;
0419     }
0420 
0421     hPosCaloHitsXY->Fill(x, y);
0422     hPosCaloHitsZX->Fill(z, x);
0423     hPosCaloHitsZY->Fill(z, y);
0424 
0425     nCaloHitsRec++;
0426 
0427     //loop over input_tower_rec and find if there is already a tower with the same cellID
0428     bool found = false;
0429     for (auto& tower : input_tower_rec) {
0430       if (tower.cellID == cellID) {
0431         tower.energy += energy;
0432         found = true;
0433         break;
0434       }
0435     }
0436     if (!found) {
0437       towersStrct tempstructT;
0438       tempstructT.energy        = energy;
0439       tempstructT.time          = time;
0440       tempstructT.posx          = x;
0441       tempstructT.posy          = y;
0442       tempstructT.posz          = z;
0443       tempstructT.cellID        = cellID;
0444       tempstructT.cellIDx       = cellIDx;
0445       tempstructT.cellIDy       = cellIDy;
0446       if (isLFHCal)
0447         tempstructT.cellIDz       = detector_layer_rz;
0448       tempstructT.tower_trueID  = 0; //TODO how to get trueID?
0449       input_tower_rec.push_back(tempstructT);
0450       input_tower_recSav.push_back(tempstructT);
0451     }
0452   }
0453   m_log->trace("LFHCal mod: nCaloHits sim  {}\t rec {}", nCaloHitsSim, nCaloHitsRec);
0454 
0455   if (nCaloHitsRec > 0) nEventsWithCaloHits++;
0456 
0457   // ===============================================================================================
0458   // sort tower arrays
0459   // ===============================================================================================
0460   hSamplingFractionEta->Fill(mceta, sumActiveCaloEnergy / (sumActiveCaloEnergy+sumPassiveCaloEnergy));
0461   std::sort(input_tower_rec.begin(), input_tower_rec.end(), &acompare);
0462   std::sort(input_tower_recSav.begin(), input_tower_recSav.end(), &acompare);
0463   std::sort(input_tower_sim.begin(), input_tower_sim.end(), &acompare);
0464 
0465   // ===============================================================================================
0466   // calculated summed hit energy for rec and sim hits
0467   // ===============================================================================================
0468 
0469   // rec hits
0470   double tot_energyRecHit = 0;
0471   for (auto& tower : input_tower_rec) {
0472     tower.energy = tower.energy; // calibrate
0473     tot_energyRecHit += tower.energy;
0474   }
0475 
0476   double samplingFractionFe = 0.037;
0477   double samplingFractionW  = 0.019;
0478   int minCellIDzDiffSamp    = 5;
0479   // sim hits
0480   double tot_energySimHit = 0;
0481   for (auto& tower : input_tower_sim) {
0482     if (tower.cellIDz < minCellIDzDiffSamp)
0483       tower.energy = tower.energy/samplingFractionW; // calibrate
0484     else
0485       tower.energy = tower.energy/samplingFractionFe; // calibrate
0486     tot_energySimHit += tower.energy;
0487   }
0488   m_log->trace("Mc E: {} \t eta: {} \t sim E rec: {}\t rec E rec: {}", mcenergy, mceta, tot_energySimHit, tot_energyRecHit);
0489 
0490   // ===============================================================================================
0491   // Fill summed hits histos
0492   // ===============================================================================================
0493   // rec hits
0494   hClusterNCells_E_eta->Fill(mcenergy, nCaloHitsRec, mceta);
0495   hClusterEcalib_E_eta->Fill(mcenergy, tot_energyRecHit/mcenergy, mceta);
0496   hClusterEcalib_E_phi->Fill(mcenergy, tot_energyRecHit/mcenergy, mcphi);
0497   // sim hits
0498   hClusterSimNCells_E_eta->Fill(mcenergy, nCaloHitsSim, mceta);
0499   hClusterESimcalib_E_eta->Fill(mcenergy, tot_energySimHit/mcenergy, mceta);
0500   hClusterESimcalib_E_phi->Fill(mcenergy, tot_energySimHit/mcenergy, mcphi);
0501 
0502   // ===============================================================================================
0503   // MA clusterization
0504   // ===============================================================================================
0505   int removedCells  = 0;
0506   float minAggE     = 0.001;
0507   float seedE       = 0.100;
0508 
0509   if (!input_tower_rec.empty()){
0510 
0511     // clean up rec array for clusterization
0512     while (input_tower_rec.at(input_tower_rec.size()-1).energy < minAggE ){
0513       input_tower_rec.pop_back();
0514       removedCells++;
0515     }
0516     m_log->trace("removed {} with E < {} GeV", removedCells, minAggE);
0517 
0518     int nclusters = 0;
0519     // vector of clusters
0520     std::vector<clustersStrct> clusters_calo;
0521     // vector of towers within the currently found cluster
0522     std::vector<towersStrct> cluster_towers;
0523     while (!input_tower_rec.empty() ) {
0524       cluster_towers.clear();
0525       clustersStrct tempstructC;
0526       // always start with highest energetic tower
0527       if(input_tower_rec.at(0).energy > seedE){
0528         m_log->trace("seed: {}\t {} \t {}", input_tower_rec.at(0).energy, input_tower_rec.at(0).cellIDx, input_tower_rec.at(0).cellIDy,  input_tower_rec.at(0).cellIDz);
0529         tempstructC = findMACluster(seedE, minAggE, input_tower_rec, cluster_towers);
0530 
0531         // determine remaining cluster properties from its towers
0532         float* showershape_eta_phi = CalculateM02andWeightedPosition(cluster_towers, tempstructC.cluster_E, 4.5);
0533         tempstructC.cluster_M02 = showershape_eta_phi[0];
0534         tempstructC.cluster_M20 = showershape_eta_phi[1];
0535         tempstructC.cluster_Eta = showershape_eta_phi[2];
0536         tempstructC.cluster_Phi = showershape_eta_phi[3];
0537         tempstructC.cluster_X = showershape_eta_phi[4];
0538         tempstructC.cluster_Y = showershape_eta_phi[5];
0539         tempstructC.cluster_Z = showershape_eta_phi[6];
0540         tempstructC.cluster_towers = cluster_towers;
0541         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 );
0542         clusters_calo.push_back(tempstructC);
0543 
0544         clusters_calo.push_back(tempstructC);
0545 
0546         nclusters++;
0547       } else {
0548         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);
0549         input_tower_rec.clear();
0550       }
0551     }
0552 
0553     // -----------------------------------------------------------------------------------------------
0554     // --------------------------- Fill LFHCal MA clusters in tree and hists -------------------------
0555     // -----------------------------------------------------------------------------------------------
0556     std::sort(clusters_calo.begin(), clusters_calo.end(), &acompareCl);
0557     m_log->info("-----> found {} clusters" , clusters_calo.size());
0558     hRecNClusters_E_eta->Fill(mcenergy, clusters_calo.size(), mceta);
0559     int iCl = 0;
0560     for (const auto cluster : clusters_calo) {
0561       if (iCl < maxNCluster && enableTreeCluster){
0562         t_lFHCal_cluster_E[iCl]       = (float)cluster.cluster_E;
0563         t_lFHCal_cluster_NCells[iCl]  = (int)cluster.cluster_NTowers;
0564         t_lFHCal_cluster_Eta[iCl]     = (float)cluster.cluster_Eta;
0565         t_lFHCal_cluster_Phi[iCl]     = (float)cluster.cluster_Phi;
0566       }
0567       hRecClusterEcalib_E_eta->Fill(mcenergy, cluster.cluster_E/mcenergy, mceta);
0568       for (int iCell = 0; iCell < (int)cluster.cluster_towers.size(); iCell++){
0569         int pSav = 0;
0570         while(cluster.cluster_towers.at(iCell).cellID !=  input_tower_recSav.at(pSav).cellID && pSav < (int)input_tower_recSav.size() ) pSav++;
0571         if (cluster.cluster_towers.at(iCell).cellID == input_tower_recSav.at(pSav).cellID)
0572           input_tower_recSav.at(pSav).tower_clusterIDA = iCl;
0573       }
0574 
0575       if (iCl == 0){
0576         hRecClusterEcalib_Ehigh_eta->Fill(mcenergy, cluster.cluster_E/mcenergy, mceta);
0577         hRecClusterNCells_Ehigh_eta->Fill(mcenergy, cluster.cluster_NTowers, mceta);
0578       }
0579       iCl++;
0580       m_log->trace("MA cluster {}:\t {} \t {}", iCl, cluster.cluster_E, cluster.cluster_NTowers);
0581     }
0582     if (iCl < maxNCluster && enableTreeCluster) t_lFHCal_clusters_N = (int)iCl;
0583 
0584     clusters_calo.clear();
0585   } else {
0586     hRecNClusters_E_eta->Fill(mcenergy, 0., mceta);
0587     if (enableTreeCluster) t_lFHCal_clusters_N = 0;
0588   }
0589 
0590   // ===============================================================================================
0591   // ------------------------------- Fill LFHCAl Island clusters in hists --------------------------
0592   // ===============================================================================================
0593   int iClF          = 0;
0594   float highestEFr  = 0;
0595   int iClFHigh      = 0;
0596 
0597   const auto &lfhcalClustersF = *(event->GetCollection<edm4eic::Cluster>(nameClusters));
0598   for (const auto cluster : lfhcalClustersF) {
0599     if (cluster.getEnergy() > highestEFr){
0600       iClFHigh    = iClF;
0601       highestEFr  = cluster.getEnergy();
0602     }
0603     hRecFClusterEcalib_E_eta->Fill(mcenergy, cluster.getEnergy()/mcenergy, mceta);
0604     m_log->trace("Island cluster {}:\t {} \t {}", iClF, cluster.getEnergy(), cluster.getNhits());
0605     iClF++;
0606   }
0607   hRecFNClusters_E_eta->Fill(mcenergy, iClF, mceta);
0608   // fill hists for highest Island cluster
0609   iClF          = 0;
0610   for (const auto cluster : lfhcalClustersF) {
0611     if (iClF == iClFHigh){
0612       hRecFClusterEcalib_Ehigh_eta->Fill(mcenergy, cluster.getEnergy()/mcenergy, mceta);
0613       hRecFClusterNCells_Ehigh_eta->Fill(mcenergy, cluster.getNhits(), mceta);
0614     }
0615     iClF++;
0616   }
0617 
0618 
0619   // ===============================================================================================
0620   // ------------------------------- Fill FECal Island clusters in hists --------------------------
0621   // ===============================================================================================
0622   int iECl            = 0;
0623   float highestEEmCl  = 0;
0624   int iEClHigh        = 0;
0625 
0626   if (enableECalCluster){
0627     try {
0628       const auto &fEMCClustersF = *(event->GetCollection<edm4eic::Cluster>("EcalEndcapPClusters"));
0629       m_log->info("-----> found fEMCClustersF:" , fEMCClustersF.size());
0630       for (const auto cluster : fEMCClustersF) {
0631         if (iECl < maxNCluster && enableTreeCluster){
0632             t_fEMC_cluster_E[iECl]       = (float)cluster.getEnergy();
0633             t_fEMC_cluster_NCells[iECl]  = (int)cluster.getNhits();
0634             t_fEMC_cluster_Eta[iECl]     = (-1.) * std::log(std::tan((float)cluster.getIntrinsicTheta() / 2.));
0635             t_fEMC_cluster_Phi[iECl]     = (float)cluster.getIntrinsicPhi();
0636         }
0637 
0638         if (cluster.getEnergy() > highestEEmCl){
0639           iEClHigh      = iECl;
0640           highestEEmCl  = cluster.getEnergy();
0641         }
0642         hRecFEmClusterEcalib_E_eta->Fill(mcenergy, cluster.getEnergy()/mcenergy, mceta);
0643         iECl++;
0644       }
0645       t_fEMC_clusters_N  = iECl;
0646       hRecFEmNClusters_E_eta->Fill(mcenergy, iECl, mceta);
0647 
0648       // fill hists for highest Island cluster
0649       iECl          = 0;
0650       for (const auto cluster : fEMCClustersF) {
0651         if (iECl == iEClHigh){
0652           hRecFEmClusterEcalib_Ehigh_eta->Fill(mcenergy, cluster.getEnergy()/mcenergy, mceta);
0653         }
0654         iECl++;
0655       }
0656     } catch (...){
0657       std::cout << "ECal clusters not in output" << std::endl;
0658       enableECalCluster = false;
0659     }
0660   }
0661   // ===============================================================================================
0662   // Write clusterizer tree & clean-up variables
0663   // ===============================================================================================
0664   if (enableTree){
0665     t_lFHCal_towers_N = (int)input_tower_recSav.size();
0666     for (int iCell = 0; iCell < (int)input_tower_recSav.size(); iCell++){
0667       m_log->trace("{} \t {} \t {} \t {} \t {} \t {}", input_tower_recSav.at(iCell).cellIDx, input_tower_recSav.at(iCell).cellIDy, input_tower_recSav.at(iCell).cellIDz , input_tower_recSav.at(iCell).energy, input_tower_recSav.at(iCell).tower_clusterIDA, input_tower_recSav.at(iCell).tower_clusterIDB  );
0668 
0669       t_lFHCal_towers_cellE[iCell]      = (float)input_tower_recSav.at(iCell).energy;
0670       t_lFHCal_towers_cellT[iCell]      = (float)input_tower_recSav.at(iCell).time;
0671       t_lFHCal_towers_cellIDx[iCell]    = (short)input_tower_recSav.at(iCell).cellIDx;
0672       t_lFHCal_towers_cellIDy[iCell]    = (short)input_tower_recSav.at(iCell).cellIDy;
0673       t_lFHCal_towers_cellIDz[iCell]    = (short)input_tower_recSav.at(iCell).cellIDz;
0674       t_lFHCal_towers_clusterIDA[iCell] = (short)input_tower_recSav.at(iCell).tower_clusterIDA;
0675       t_lFHCal_towers_clusterIDB[iCell] = (short)input_tower_recSav.at(iCell).tower_clusterIDB;
0676       t_lFHCal_towers_cellTrueID[iCell] = (int)input_tower_recSav.at(iCell).tower_trueID;
0677     }
0678 
0679     event_tree->Fill();
0680 
0681     t_lFHCal_towers_N = 0;
0682     for (Int_t itow = 0; itow < maxNTowers; itow++){
0683       t_lFHCal_towers_cellE[itow]       = 0;
0684       t_lFHCal_towers_cellT[itow]       = 0;
0685       t_lFHCal_towers_cellIDx[itow]     = 0;
0686       t_lFHCal_towers_cellIDy[itow]     = 0;
0687       t_lFHCal_towers_cellIDz[itow]     = 0;
0688       t_lFHCal_towers_clusterIDA[itow]  = 0;
0689       t_lFHCal_towers_clusterIDB[itow]  = 0;
0690       t_lFHCal_towers_cellTrueID[itow]  = 0;
0691     }
0692   }
0693 
0694   // ===============================================================================================
0695   // Write cluster tree & clean-up variables
0696   // ===============================================================================================
0697   if (enableTreeCluster){
0698     cluster_tree->Fill();
0699 
0700     t_mc_N                = 0;
0701     t_lFHCal_clusters_N   = 0;
0702     t_fEMC_clusters_N    = 0;
0703     for (Int_t iMC = 0; iMC < maxNMC; iMC++){
0704       t_mc_E[iMC]                   = 0;
0705       t_mc_Phi[iMC]                 = 0;
0706       t_mc_Eta[iMC]                 = 0;
0707     }
0708     for (Int_t iCl = 0; iCl < maxNCluster; iCl++){
0709       t_lFHCal_cluster_E[iCl]       = 0;
0710       t_lFHCal_cluster_NCells[iCl]  = 0;
0711       t_lFHCal_cluster_Eta[iCl]     = 0;
0712       t_lFHCal_cluster_Phi[iCl]     = 0;
0713       t_fEMC_cluster_E[iCl]        = 0;
0714       t_fEMC_cluster_NCells[iCl]   = 0;
0715       t_fEMC_cluster_Eta[iCl]      = 0;
0716       t_fEMC_cluster_Phi[iCl]      = 0;
0717     }
0718   }
0719 
0720 }
0721 
0722 //******************************************************************************************//
0723 // Finish
0724 //******************************************************************************************//
0725 void lfhcal_studiesProcessor::Finish() {
0726   std::cout << "------> LFHCal " << nEventsWithCaloHits << " with calo info present"<< std::endl;
0727   // Do any final calculations here.
0728 
0729   if (enableTree) {
0730     delete[] t_lFHCal_towers_cellE;
0731     delete[] t_lFHCal_towers_cellT;
0732     delete[] t_lFHCal_towers_cellIDx;
0733     delete[] t_lFHCal_towers_cellIDy;
0734     delete[] t_lFHCal_towers_cellIDz;
0735     delete[] t_lFHCal_towers_clusterIDA;
0736     delete[] t_lFHCal_towers_clusterIDB;
0737     delete[] t_lFHCal_towers_cellTrueID;
0738   }
0739 
0740   if (enableTreeCluster) {
0741     delete[] t_mc_E;
0742     delete[] t_mc_Phi;
0743     delete[] t_mc_Eta;
0744     delete[] t_lFHCal_cluster_E;
0745     delete[] t_lFHCal_cluster_NCells;
0746     delete[] t_lFHCal_cluster_Phi;
0747     delete[] t_lFHCal_cluster_Eta;
0748     delete[] t_fEMC_cluster_E;
0749     delete[] t_fEMC_cluster_NCells;
0750     delete[] t_fEMC_cluster_Eta;
0751     delete[] t_fEMC_cluster_Phi;
0752   }
0753 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!