Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-01-30 10:30:54 

0001 // ----------------------------------------------------------------------------
0002 // 'CheckClusterSplittingProcessor.cc'
0003 // Derek Anderson
0004 // 04.11.2024
0005 //
0006 // EICrecon plugin to collect some metrics on cluster splitting
0007 // across several calorimeters.
0008 // ----------------------------------------------------------------------------
0009 
0010 #define CHECKCLUSTERSPLITTING_CC
0011 
0012 // plugin definition
0013 #include "CheckClusterSplittingProcessor.h"
0014 
0015 // the following just makes this a JANA plugin
0016 extern "C" {
0017   void InitPlugin(JApplication *app) {
0018     InitJANAPlugin(app);
0019     app -> Add( new CheckClusterSplittingProcessor );
0020   }
0021 }
0022 
0023 // make common namespaces implicit
0024 using namespace std;
0025 
0026 
0027 
0028 // public methods -------------------------------------------------------------
0029 
0030 void CheckClusterSplittingProcessor::InitWithGlobalRootLock() {
0031 
0032   // create output directory
0033   auto rootfile_svc = GetApplication() -> GetService<RootFile_service>();
0034   auto rootfile     = rootfile_svc -> GetHistFile();
0035   rootfile -> mkdir("CheckClusterSplitting") -> cd();
0036 
0037   // turn on errors
0038   TH1::SetDefaultSumw2(true);
0039   TH2::SetDefaultSumw2(true);
0040 
0041   // create output histograms
0042   BuildEvtHistograms();
0043   BuildTrkHistograms();
0044   BuildCalHistograms();
0045   return;
0046 
0047 }  // end 'InitWithGlobalRootLock()'
0048 
0049 
0050 
0051 void CheckClusterSplittingProcessor::ProcessSequential(const shared_ptr<const JEvent>& event) {
0052 
0053   // grab particles and track projections
0054   const auto& particles   = *(event -> GetCollection<edm4eic::ReconstructedParticle>(m_config.parCollect.data()));
0055   const auto& projections = *(event -> GetCollection<edm4eic::TrackSegment>(m_config.projCollect.data()));
0056 
0057   // get particle energy
0058   double ePar = 0.;
0059   for (auto particle : particles) {
0060     if (particle.getType() == 1) {
0061       ePar = particle.getEnergy();
0062       break;
0063     }
0064   }  // end particle loop
0065 
0066   // loop over calorimeters
0067   for (int iCalo = 0; iCalo < m_config.vecCalos.size(); iCalo++) {
0068 
0069     // grab collection
0070     const auto& clusters = *(event -> GetCollection<edm4eic::Cluster>(m_config.vecCalos[iCalo].first.data()));
0071     const int   idCalo   = GetCaloID(iCalo);
0072 
0073     // grab relevant track projections
0074     GetProjections(projections, idCalo);
0075     if (m_vecTrkProj.size() <= 0) continue;
0076 
0077     // determine lead cluster energy and cluster sum
0078     double eSum  = 0.;
0079     double eLead = 0.;
0080     for (auto cluster : clusters) {
0081       if (cluster.getEnergy() > eLead) {
0082         eLead = cluster.getEnergy();
0083       }
0084       eSum += cluster.getEnergy();
0085     }  // end cluster loop
0086     if (eSum <= 0.) continue;
0087 
0088     // run calculations
0089     int32_t nClusters  = 0;
0090     int32_t nClustAb90 = 0;
0091     int32_t nClustBe90 = 0;
0092     int32_t nClustRec  = 0;
0093     for (auto cluster : clusters) {
0094 
0095       // create struct to hold histogram content
0096       Hist::CalContent cContent;
0097       cContent.GetInfo(cluster);
0098 
0099       // skip hit-less clusters
0100       if (cluster.getHits().size() <= 0) continue;
0101 
0102       // calculate spatial extent of cluster
0103       FPair widths   = GetClusterWidths(cluster);
0104       FPair etaRange = GetEtaRange(cluster);
0105       FPair phiRange = GetPhiRange(cluster);
0106 
0107       // find closest projection
0108       bool    foundMatch = false;
0109       double  drMin      = 999.;
0110       double  drSigMin   = 999.;
0111       int32_t nProj      = 0;
0112       int32_t iMatch     = -1;
0113       int32_t iProject   = 0;
0114       for (auto project : m_vecTrkProj) {
0115 
0116         // get projection eta/phi
0117         const double projPhi = atan2(project.position.y, project.position.x);
0118         const double projEta = edm4hep::utils::eta(project.position);
0119 
0120         // check if in cluster eta/phi extent
0121         const bool isInEtaRange = ((projEta >= etaRange.first) && (projEta < etaRange.second));
0122         const bool isInPhiRange = ((projPhi >= phiRange.first) && (projPhi < phiRange.second));
0123         if (!isInEtaRange || !isInPhiRange) {
0124           continue;
0125         } else {
0126           ++nProj;
0127         }
0128 
0129         // grab projection info
0130         Hist::TrkContent tContent;
0131         tContent.GetInfo(project);
0132 
0133         // get distances
0134         const double dEta = (projEta - cContent.eta);
0135         const double dPhi = (projPhi - cContent.phi);
0136 
0137         // grab remaining projection info
0138         tContent.dr     = hypot(dEta, dPhi);
0139         tContent.drSig  = hypot(dEta/widths.first, dPhi/widths.second);
0140         tContent.eOverP = cluster.getEnergy() / tContent.p;
0141         tContent.deposit = tContent.p * m_config.vecAvgEP[iCalo];
0142 
0143         // fill all projection histograms
0144         FillTrkHistograms(iCalo, Hist::TType::TrkAll, tContent);
0145 
0146         // if projection is closer to cluster, set current match
0147         if (tContent.dr <= drMin) {
0148           drMin      = tContent.dr;
0149           drSigMin   = tContent.drSig;
0150           iMatch     = iProject;
0151           foundMatch = true;
0152         }
0153         ++iProject;
0154 
0155       }  // end projection loop
0156 
0157       // grab info of matched projection
0158       Hist::TrkContent mContent;
0159       if (foundMatch) {
0160         mContent.GetInfo( m_vecTrkProj.at(iMatch) );
0161         mContent.dr      = drMin;
0162         mContent.drSig   = drSigMin;
0163         mContent.eOverP  = cluster.getEnergy() / mContent.p;
0164         mContent.deposit = mContent.p * m_config.vecAvgEP[iCalo];
0165         FillTrkHistograms(iCalo, Hist::TType::TrkMatch, mContent);
0166       }
0167 
0168       // calculate significance of cluster energy
0169       const double sigma = (cluster.getEnergy() - mContent.deposit) / m_config.vecSigEP[iCalo];
0170 
0171       // grab remaining cluster information
0172       cContent.nProj  = nProj;
0173       cContent.frac   = cluster.getEnergy() / eSum;
0174       cContent.purity = cluster.getEnergy() / ePar;
0175       cContent.eOverP = foundMatch ? mContent.eOverP : numeric_limits<double>::max();
0176       cContent.sigma  = foundMatch ? sigma : numeric_limits<double>::max();
0177 
0178       // check purity
0179       const bool isBelow90 = (cContent.purity < 0.9);
0180 
0181       // if purity < 90%, see if we can recover 90%
0182       bool    didRecover = false;
0183       float   drRecover  = numeric_limits<double>::max();
0184       int32_t nRecover   = numeric_limits<int32_t>::max();
0185       if (isBelow90) {
0186 
0187         // build lists of clusters vs. energy, delta-r, and if already checked
0188         GetClustersAndDr(clusters, cContent, cluster.getObjectID().index);
0189 
0190         // now try to recover if there are other clusters to check
0191         const uint32_t nToCheck = m_mapClustToDr.size();
0192         if (nToCheck > 0) {
0193           TryToRecover(cContent.ene, ePar, didRecover, drRecover, nRecover);
0194         }
0195       }  // end if purity below 90%
0196 
0197       // grab recover info
0198       cContent.nAdd90  = nRecover;
0199       cContent.drAdd90 = drRecover;
0200       if (didRecover) {
0201         cContent.perAdd90 = (double) nRecover / (double) clusters.size();
0202         ++nClustRec;
0203       } else {
0204         cContent.perAdd90 = numeric_limits<double>::max();
0205       }
0206 
0207       // fill cluster histograms
0208       if (isBelow90) {
0209         FillCalHistograms(iCalo, Hist::CType::CalBelow90, cContent);
0210         ++nClustBe90;
0211       } else {
0212         FillCalHistograms(iCalo, Hist::CType::CalAbove90, cContent);
0213         ++nClustAb90;
0214       }
0215       FillCalHistograms(iCalo, Hist::CType::CalAll, cContent);
0216       ++nClusters;
0217 
0218     }  // end cluster loop
0219 
0220     // set event level-information
0221     Hist::EvtContent eContent;
0222     eContent.nTrk   = m_vecTrkProj.size();
0223     eContent.nClust = nClusters;
0224     eContent.nCl90  = nClustAb90;
0225     eContent.nClB90 = nClustBe90;
0226     eContent.nClRec = nClustRec;
0227     eContent.ePar   = ePar;
0228     eContent.eLead  = eLead;
0229     eContent.eTot   = eSum;
0230     eContent.lOverT = eLead / eSum;
0231     eContent.lOverP = eLead / ePar;
0232     eContent.tOverP = eSum / ePar;
0233 
0234     // fill event-level histograms
0235     FillEvtHistograms(iCalo, Hist::EType::EvtAll, eContent);
0236 
0237   }  // end calorimeter loop
0238   return;
0239 
0240 }  // end 'ProcessSequential(shared_ptr<JEvent>&)'
0241 
0242 
0243 
0244 void CheckClusterSplittingProcessor::FinishWithGlobalRootLock() {
0245 
0246   /* nothing to do */
0247   return;
0248 
0249 }  // end 'FinishWithGlobalRootLock()'
0250 
0251 
0252 
0253 // private methods ------------------------------------------------------------
0254 
0255 void CheckClusterSplittingProcessor::BuildEvtHistograms() {
0256 
0257   // histogram binning/labales
0258   VecAxisDef vecEvtAxes = {
0259     make_pair("N_{trk}",                  make_tuple(100, -0.5, 99.5)),
0260     make_pair("N_{clust}",                make_tuple(100, -0.5, 99.5)),
0261     make_pair("E_{par} [GeV]",            make_tuple(100, 0., 50.)),
0262     make_pair("E_{lead} [GeV]",           make_tuple(100, 0., 50.)),
0263     make_pair("#SigmaE_{clust} [GeV]",    make_tuple(100, 0. ,50.)),
0264     make_pair("E_{lead}/#SigmaE_{clust}", make_tuple(250, 0., 5.)),
0265     make_pair("E_{lead}/E_{par}",         make_tuple(250, 0., 5.)),
0266     make_pair("#SigmaE_{clust}/E_{par}",  make_tuple(250, 0., 5.)),
0267     make_pair("N_{clust}^{>90%}",         make_tuple(100, -0.5, 99.5)),
0268     make_pair("N_{clust}^{<90%}",         make_tuple(100, -0.5, 99.5)),
0269     make_pair("N_{clust}^{recovered}",    make_tuple(100, -0.5, 99.5))
0270   };
0271 
0272   // type labels
0273   vector<string> vecEvtTypes = {
0274     "All"
0275   };
0276 
0277   // histogram definitions
0278   VecHistDef1D vecEvtDef1D = {
0279     make_tuple("EvtNTrk",            vecEvtAxes[Hist::EVar::EvtNTr]),
0280     make_tuple("EvtNClust",          vecEvtAxes[Hist::EVar::EvtNCl]),
0281     make_tuple("EvtParEne",          vecEvtAxes[Hist::EVar::EvtPE]),
0282     make_tuple("EvtLeadEne",         vecEvtAxes[Hist::EVar::EvtLE]),
0283     make_tuple("EvtTotalEne",        vecEvtAxes[Hist::EVar::EvtTE]),
0284     make_tuple("EvtLeadTotEneFrac",  vecEvtAxes[Hist::EVar::EvtLT]),
0285     make_tuple("EvtLeadParEneFrac",  vecEvtAxes[Hist::EVar::EvtLP]),
0286     make_tuple("EvtTotParEneFrac",   vecEvtAxes[Hist::EVar::EvtTP]),
0287     make_tuple("EvtNClustAbove90",   vecEvtAxes[Hist::EVar::EvtNA90]),
0288     make_tuple("EvtNClustBelow90",   vecEvtAxes[Hist::EVar::EvtNB90]),
0289     make_tuple("EvtNClustRecovered", vecEvtAxes[Hist::EVar::EvtNR])
0290   };
0291 
0292   // make histograms
0293   m_vecEvtH1D.resize( m_config.vecCalos.size() );
0294   for (size_t iCalo = 0; iCalo < m_config.vecCalos.size(); iCalo++) {
0295 
0296     m_vecEvtH1D[iCalo].resize( vecEvtTypes.size() );
0297     for (size_t iEType = 0; iEType < vecEvtTypes.size(); iEType++) {
0298 
0299       // make 1d histograms
0300       for (auto evtDef1D : vecEvtDef1D) {
0301 
0302         // make name
0303         string name("h");
0304         name += m_config.vecCalos[iCalo].second;
0305         name += vecEvtTypes[iEType];
0306         name += get<0>(evtDef1D);
0307 
0308         // make title
0309         string title(";");
0310         title += get<1>(evtDef1D).first;
0311         title += ";counts";
0312 
0313         // create histogram
0314         m_vecEvtH1D[iCalo][iEType].push_back(
0315           new TH1D(
0316             name.data(),
0317             title.data(),
0318             get<0>( get<1>(evtDef1D).second ),
0319             get<1>( get<1>(evtDef1D).second ),
0320             get<2>( get<1>(evtDef1D).second )
0321           )
0322         );
0323       }  // end 1d hist def loop
0324     }  // end type loop
0325   }  // end calo loop
0326   return;
0327 
0328 }  // end 'BuildEvtHistograms()'
0329 
0330 
0331 
0332 void CheckClusterSplittingProcessor::BuildTrkHistograms() {
0333 
0334   // histogram binning/labales
0335   VecAxisDef vecTrkAxes = {
0336     make_pair("r_{x} [mm]",      make_tuple(3000, -3000., 3000.)),
0337     make_pair("r_{y} [mm]",      make_tuple(3000, -3000., 3000.)),
0338     make_pair("r_{z} [mm]",      make_tuple(3000, -3000., 3000.)),
0339     make_pair("#eta",            make_tuple(100, -5., 5.)),
0340     make_pair("#varphi",         make_tuple(180, -3.15, 3.15)),
0341     make_pair("p [GeV/c]",       make_tuple(100, 0., 50.)),
0342     make_pair("#Deltar",         make_tuple(500, 0., 5.)),
0343     make_pair("#Delta#tilde{r}", make_tuple(500, 0., 5.)),
0344     make_pair("E/p",             make_tuple(250, 0., 5.)),
0345     make_pair("p #times <E/p>",  make_tuple(100, 0., 50.))
0346   };
0347 
0348   // type labels
0349   vector<string> vecTrkTypes = {
0350     "All",
0351     "Match"
0352   };
0353 
0354   // histogram definitions
0355   VecHistDef1D vecTrkDef1D = {
0356     make_tuple("TrkPosX",      vecTrkAxes[Hist::TVar::TrkRX]),
0357     make_tuple("TrkPosY",      vecTrkAxes[Hist::TVar::TrkRY]),
0358     make_tuple("TrkPosZ",      vecTrkAxes[Hist::TVar::TrkRZ]),
0359     make_tuple("TrkEta",       vecTrkAxes[Hist::TVar::TrkH]),
0360     make_tuple("TrkPhi",       vecTrkAxes[Hist::TVar::TrkF]),
0361     make_tuple("TrkMom",       vecTrkAxes[Hist::TVar::TrkP]),
0362     make_tuple("TrkDeltaR",    vecTrkAxes[Hist::TVar::TrkDR]),
0363     make_tuple("TrkDeltaRSig", vecTrkAxes[Hist::TVar::TrkDRS]),
0364     make_tuple("TrkEoverP",    vecTrkAxes[Hist::TVar::TrkEP]),
0365     make_tuple("TrkDeposit",   vecTrkAxes[Hist::TVar::TrkDep])
0366   };
0367   VecHistDef2D vecTrkDef2D = {
0368     make_tuple("TrkPosYvsX",   vecTrkAxes[Hist::TVar::TrkRX],  vecTrkAxes[Hist::TVar::TrkRY]),
0369     make_tuple("TrkMomVsEta",  vecTrkAxes[Hist::TVar::TrkH],   vecTrkAxes[Hist::TVar::TrkP]),
0370     make_tuple("TrkEtaVsPhi",  vecTrkAxes[Hist::TVar::TrkF],   vecTrkAxes[Hist::TVar::TrkH]),
0371     make_tuple("TrkEPvsDR",    vecTrkAxes[Hist::TVar::TrkDR],  vecTrkAxes[Hist::TVar::TrkEP]),
0372     make_tuple("TrkEPvsDRSig", vecTrkAxes[Hist::TVar::TrkDRS], vecTrkAxes[Hist::TVar::TrkEP])
0373   };
0374 
0375   // make histograms
0376   m_vecTrkH1D.resize( m_config.vecCalos.size() );
0377   m_vecTrkH2D.resize( m_config.vecCalos.size() );
0378   for (size_t iCalo = 0; iCalo < m_config.vecCalos.size(); iCalo++) {
0379 
0380     m_vecTrkH1D[iCalo].resize( vecTrkTypes.size() );
0381     m_vecTrkH2D[iCalo].resize( vecTrkTypes.size() );
0382     for (size_t iTType = 0; iTType < vecTrkTypes.size(); iTType++) {
0383 
0384       // make 1d histograms
0385       for (auto trkDef1D : vecTrkDef1D) {
0386 
0387         // make name
0388         string name("h");
0389         name += m_config.vecCalos[iCalo].second;
0390         name += vecTrkTypes[iTType];
0391         name += get<0>(trkDef1D);
0392 
0393         // make title
0394         string title(";");
0395         title += get<1>(trkDef1D).first;
0396         title += ";counts";
0397 
0398         // create histogram
0399         m_vecTrkH1D[iCalo][iTType].push_back(
0400           new TH1D(
0401             name.data(),
0402             title.data(),
0403             get<0>( get<1>(trkDef1D).second ),
0404             get<1>( get<1>(trkDef1D).second ),
0405             get<2>( get<1>(trkDef1D).second )
0406           )
0407         );
0408       }  // end 1d hist def loop
0409 
0410       // make 2d histograms
0411       for (auto trkDef2D : vecTrkDef2D) {
0412 
0413         // make name
0414         string name("h");
0415         name += m_config.vecCalos[iCalo].second;
0416         name += vecTrkTypes[iTType];
0417         name += get<0>(trkDef2D);
0418 
0419         // make title
0420         string title(";");
0421         title += get<1>(trkDef2D).first;
0422         title += ";";
0423         title += get<2>(trkDef2D).first;
0424         title += ";counts";
0425 
0426         // create histogram
0427         m_vecTrkH2D[iCalo][iTType].push_back(
0428           new TH2D(
0429             name.data(),
0430             title.data(),
0431             get<0>( get<1>(trkDef2D).second ),
0432             get<1>( get<1>(trkDef2D).second ),
0433             get<2>( get<1>(trkDef2D).second ),
0434             get<0>( get<2>(trkDef2D).second ),
0435             get<1>( get<2>(trkDef2D).second ),
0436             get<2>( get<2>(trkDef2D).second )
0437           )
0438         );
0439       }  // end 2d hist def loop
0440     }  // end type loop
0441   }  // end calo loop
0442   return;
0443 
0444 }  // end 'BuildTrkHistograms()'
0445 
0446 
0447 
0448 void CheckClusterSplittingProcessor::BuildCalHistograms() {
0449 
0450   // histogram binning/labales
0451   VecAxisDef vecCalAxes = {
0452     make_pair("N_{hit}",                   make_tuple(50, -0.5, 49.5)),
0453     make_pair("N_{proj}",                  make_tuple(20, -0.5, 19.5)),
0454     make_pair("r_{x} [mm]",                make_tuple(3000, -3000., 3000.)),
0455     make_pair("r_{y} [mm]",                make_tuple(3000, -3000., 3000.)),
0456     make_pair("r_{z} [mm]",                make_tuple(3000, -3000., 3000.)),
0457     make_pair("#eta",                      make_tuple(100, -5., 5.)),
0458     make_pair("#varphi",                   make_tuple(180, -3.15, 3.15)),
0459     make_pair("E [GeV]",                   make_tuple(100, 0., 50.)),
0460     make_pair("E_{clust}/#SigmaE_{clust}", make_tuple(250, 0., 5.)),
0461     make_pair("E_{clust}/E_{par}",         make_tuple(250, 0., 5.)),
0462     make_pair("E/p",                       make_tuple(250, 0., 5.)),
0463     make_pair("S(E_{clust})",              make_tuple(20, -10., 10.)),
0464     make_pair("N_{90}",                    make_tuple(100, -0.5, 99.5)),
0465     make_pair("N_{90}/N_{clust}",          make_tuple(20, 0., 2.)),
0466     make_pair("#Deltar_{90}",              make_tuple(60, 0., 3.))
0467   };
0468 
0469   // type labels
0470   vector<string> vecCalTypes = {
0471     "All",
0472     "AboveNinety",
0473     "BelowNinety"
0474   };
0475 
0476   // histogram definitions
0477   VecHistDef1D vecCalDef1D = {
0478     make_tuple("ClustNHit",    vecCalAxes[Hist::CVar::CalNH]),
0479     make_tuple("ClustNProj",   vecCalAxes[Hist::CVar::CalNP]),
0480     make_tuple("ClustPosX",    vecCalAxes[Hist::CVar::CalRX]),
0481     make_tuple("ClustPosY",    vecCalAxes[Hist::CVar::CalRY]),
0482     make_tuple("ClustPosZ",    vecCalAxes[Hist::CVar::CalRZ]),
0483     make_tuple("ClustEta",     vecCalAxes[Hist::CVar::CalH]),
0484     make_tuple("ClustPhi",     vecCalAxes[Hist::CVar::CalF]),
0485     make_tuple("ClustEne",     vecCalAxes[Hist::CVar::CalE]),
0486     make_tuple("ClustEpsilon", vecCalAxes[Hist::CVar::CalEps]),
0487     make_tuple("ClustPurity",  vecCalAxes[Hist::CVar::CalRho]),
0488     make_tuple("ClustEoverP",  vecCalAxes[Hist::CVar::CalEP]),
0489     make_tuple("ClustCutSig",  vecCalAxes[Hist::CVar::CalSig]),
0490     make_tuple("ClustNAdd90",  vecCalAxes[Hist::CVar::CalN90]),
0491     make_tuple("ClustPAdd90",  vecCalAxes[Hist::CVar::CalP90]),
0492     make_tuple("ClustDrAdd90", vecCalAxes[Hist::CVar::CalDR90])
0493   };
0494   VecHistDef2D vecCalDef2D = {
0495     make_tuple("ClustPosYvsX",  vecCalAxes[Hist::CVar::CalRX], vecCalAxes[Hist::CVar::CalRY]),
0496     make_tuple("ClustEneVsEta", vecCalAxes[Hist::CVar::CalH],  vecCalAxes[Hist::CVar::CalE]),
0497     make_tuple("ClustEneVsPhi", vecCalAxes[Hist::CVar::CalF],  vecCalAxes[Hist::CVar::CalF]),
0498     make_tuple("ClustEtaVsPhi", vecCalAxes[Hist::CVar::CalF],  vecCalAxes[Hist::CVar::CalH]),
0499     make_tuple("ClustPurVsEne", vecCalAxes[Hist::CVar::CalE],  vecCalAxes[Hist::CVar::CalRho]),
0500     make_tuple("ClustN90VsEp",  vecCalAxes[Hist::CVar::CalEP], vecCalAxes[Hist::CVar::CalN90]),
0501     make_tuple("ClustP90VsEp",  vecCalAxes[Hist::CVar::CalEP], vecCalAxes[Hist::CVar::CalP90]),
0502     make_tuple("ClustDR90VsEp", vecCalAxes[Hist::CVar::CalEP], vecCalAxes[Hist::CVar::CalDR90])
0503   };
0504 
0505   // make histograms
0506   m_vecCalH1D.resize( m_config.vecCalos.size() );
0507   m_vecCalH2D.resize( m_config.vecCalos.size() );
0508   for (size_t iCalo = 0; iCalo < m_config.vecCalos.size(); iCalo++) {
0509 
0510     m_vecCalH1D[iCalo].resize( vecCalTypes.size() );
0511     m_vecCalH2D[iCalo].resize( vecCalTypes.size() );
0512     for (size_t iCType = 0; iCType < vecCalTypes.size(); iCType++) {
0513 
0514       // make 1d histograms
0515       for (auto calDef1D : vecCalDef1D) {
0516 
0517         // make name
0518         string name("h");
0519         name += m_config.vecCalos[iCalo].second;
0520         name += vecCalTypes[iCType];
0521         name += get<0>(calDef1D);
0522 
0523         // make title
0524         string title(";");
0525         title += get<1>(calDef1D).first;
0526         title += ";counts";
0527 
0528         // create histogram
0529         m_vecCalH1D[iCalo][iCType].push_back(
0530           new TH1D(
0531             name.data(),
0532             title.data(),
0533             get<0>( get<1>(calDef1D).second ),
0534             get<1>( get<1>(calDef1D).second ),
0535             get<2>( get<1>(calDef1D).second )
0536           )
0537         );
0538       }  // end 1d hist def loop
0539 
0540       // make 2d histograms
0541       for (auto calDef2D : vecCalDef2D) {
0542 
0543         // make name
0544         string name("h");
0545         name += m_config.vecCalos[iCalo].second;
0546         name += vecCalTypes[iCType];
0547         name += get<0>(calDef2D);
0548 
0549         // make title
0550         string title(";");
0551         title += get<1>(calDef2D).first;
0552         title += ";";
0553         title += get<2>(calDef2D).first;
0554         title += ";counts";
0555 
0556         // create histogram
0557         m_vecCalH2D[iCalo][iCType].push_back(
0558           new TH2D(
0559             name.data(),
0560             title.data(),
0561             get<0>( get<1>(calDef2D).second ),
0562             get<1>( get<1>(calDef2D).second ),
0563             get<2>( get<1>(calDef2D).second ),
0564             get<0>( get<2>(calDef2D).second ),
0565             get<1>( get<2>(calDef2D).second ),
0566             get<2>( get<2>(calDef2D).second )
0567           )
0568         );
0569       }  // end 2d hist def loop
0570     }  // end type label loop
0571   }  // end calo label loop
0572   return;
0573 
0574 }  // end 'BuildHistograms()'
0575 
0576 
0577 
0578 void CheckClusterSplittingProcessor::FillEvtHistograms(const int calo, const int type, const Hist::EvtContent& content) {
0579 
0580   // fill 1d hisgorams
0581   m_vecEvtH1D.at(calo).at(type).at(Hist::E1D::EvtNTrk)       -> Fill(content.nTrk);
0582   m_vecEvtH1D.at(calo).at(type).at(Hist::E1D::EvtNClust)     -> Fill(content.nClust);
0583   m_vecEvtH1D.at(calo).at(type).at(Hist::E1D::EvtParEne)     -> Fill(content.ePar);
0584   m_vecEvtH1D.at(calo).at(type).at(Hist::E1D::EvtLeadEne)    -> Fill(content.eLead);
0585   m_vecEvtH1D.at(calo).at(type).at(Hist::E1D::EvtTotEne)     -> Fill(content.eTot);
0586   m_vecEvtH1D.at(calo).at(type).at(Hist::E1D::EvtLeadDivTot) -> Fill(content.lOverT);
0587   m_vecEvtH1D.at(calo).at(type).at(Hist::E1D::EvtLeadDivPar) -> Fill(content.lOverP);
0588   m_vecEvtH1D.at(calo).at(type).at(Hist::E1D::EvtTotDivPar)  -> Fill(content.tOverP);
0589   m_vecEvtH1D.at(calo).at(type).at(Hist::E1D::EvtNClAb90)    -> Fill(content.nCl90);
0590   m_vecEvtH1D.at(calo).at(type).at(Hist::E1D::EvtNClBe90)    -> Fill(content.nClB90);
0591   m_vecEvtH1D.at(calo).at(type).at(Hist::E1D::EvtNClRec)     -> Fill(content.nClRec);
0592   return;
0593 
0594 }  // end 'FillEvtHistograms(int, int, Hist::EvtContent&)'
0595 
0596 
0597 
0598 void CheckClusterSplittingProcessor::FillTrkHistograms(const int calo, const int type, const Hist::TrkContent& content) {
0599 
0600   // fill 1d histograms
0601   m_vecTrkH1D.at(calo).at(type).at(Hist::T1D::TrkPosX)        -> Fill(content.rx);
0602   m_vecTrkH1D.at(calo).at(type).at(Hist::T1D::TrkPosY)        -> Fill(content.ry);
0603   m_vecTrkH1D.at(calo).at(type).at(Hist::T1D::TrkPosZ)        -> Fill(content.rz);
0604   m_vecTrkH1D.at(calo).at(type).at(Hist::T1D::TrkEta)         -> Fill(content.eta);
0605   m_vecTrkH1D.at(calo).at(type).at(Hist::T1D::TrkPhi)         -> Fill(content.phi);
0606   m_vecTrkH1D.at(calo).at(type).at(Hist::T1D::TrkMom)         -> Fill(content.p);
0607   m_vecTrkH1D.at(calo).at(type).at(Hist::T1D::TrkDeltaR)      -> Fill(content.dr);
0608   m_vecTrkH1D.at(calo).at(type).at(Hist::T1D::TrkDeltaRScale) -> Fill(content.drSig);
0609   m_vecTrkH1D.at(calo).at(type).at(Hist::T1D::TrkEOverP)      -> Fill(content.eOverP);
0610   m_vecTrkH1D.at(calo).at(type).at(Hist::T1D::TrkDeposit)     -> Fill(content.deposit);
0611 
0612   // fill 2d histograms
0613   m_vecTrkH2D.at(calo).at(type).at(Hist::T2D::TrkPosYvsX)   -> Fill(content.rx,    content.rx);
0614   m_vecTrkH2D.at(calo).at(type).at(Hist::T2D::TrkMomVsEta)  -> Fill(content.eta,   content.p);
0615   m_vecTrkH2D.at(calo).at(type).at(Hist::T2D::TrkEtaVsPhi)  -> Fill(content.phi,   content.eta);
0616   m_vecTrkH2D.at(calo).at(type).at(Hist::T2D::TrkEPvsDR)    -> Fill(content.dr,    content.eOverP);
0617   m_vecTrkH2D.at(calo).at(type).at(Hist::T2D::TrkEPvsDRSig) -> Fill(content.drSig, content.eOverP);
0618   return;
0619 
0620 }  // end 'FillTrkHisotgrams(int, int, Hist::TrkContent&)'
0621 
0622 
0623 
0624 void CheckClusterSplittingProcessor::FillCalHistograms(const int calo, const int type, const Hist::CalContent& content) {
0625 
0626   // fill 1d histograms
0627   m_vecCalH1D.at(calo).at(type).at(Hist::C1D::CalNHit)     -> Fill(content.nHit);
0628   m_vecCalH1D.at(calo).at(type).at(Hist::C1D::CalNProj)    -> Fill(content.nProj);
0629   m_vecCalH1D.at(calo).at(type).at(Hist::C1D::CalPosX)     -> Fill(content.rx);
0630   m_vecCalH1D.at(calo).at(type).at(Hist::C1D::CalPosY)     -> Fill(content.ry);
0631   m_vecCalH1D.at(calo).at(type).at(Hist::C1D::CalPosZ)     -> Fill(content.rz);
0632   m_vecCalH1D.at(calo).at(type).at(Hist::C1D::CalEta)      -> Fill(content.eta);
0633   m_vecCalH1D.at(calo).at(type).at(Hist::C1D::CalPhi)      -> Fill(content.phi);
0634   m_vecCalH1D.at(calo).at(type).at(Hist::C1D::CalEne)      -> Fill(content.ene);
0635   m_vecCalH1D.at(calo).at(type).at(Hist::C1D::CalFraction) -> Fill(content.frac);
0636   m_vecCalH1D.at(calo).at(type).at(Hist::C1D::CalPurity)   -> Fill(content.purity);
0637   m_vecCalH1D.at(calo).at(type).at(Hist::C1D::CalEOverP)   -> Fill(content.eOverP);
0638   m_vecCalH1D.at(calo).at(type).at(Hist::C1D::CalSigma)    -> Fill(content.sigma);
0639   m_vecCalH1D.at(calo).at(type).at(Hist::C1D::CalNAdd90)   -> Fill(content.nAdd90);
0640   m_vecCalH1D.at(calo).at(type).at(Hist::C1D::CalPAdd90)   -> Fill(content.perAdd90);
0641   m_vecCalH1D.at(calo).at(type).at(Hist::C1D::CalDRAdd90)  -> Fill(content.drAdd90);
0642 
0643   // fill 2d histograms
0644   m_vecCalH2D.at(calo).at(type).at(Hist::C2D::CalPosYvsX)  -> Fill(content.rx,     content.ry);
0645   m_vecCalH2D.at(calo).at(type).at(Hist::C2D::CalEneVsEta) -> Fill(content.eta,    content.ene);
0646   m_vecCalH2D.at(calo).at(type).at(Hist::C2D::CalEneVsPhi) -> Fill(content.phi,    content.ene);
0647   m_vecCalH2D.at(calo).at(type).at(Hist::C2D::CalEtaVsPhi) -> Fill(content.phi,    content.eta);
0648   m_vecCalH2D.at(calo).at(type).at(Hist::C2D::CalPurVsEne) -> Fill(content.ene,    content.purity);
0649   m_vecCalH2D.at(calo).at(type).at(Hist::C2D::CalN90VsEP)  -> Fill(content.eOverP, content.nAdd90);
0650   m_vecCalH2D.at(calo).at(type).at(Hist::C2D::CalP90VsEP)  -> Fill(content.eOverP, content.perAdd90);
0651   m_vecCalH2D.at(calo).at(type).at(Hist::C2D::CalDR90VsEP) -> Fill(content.eOverP, content.drAdd90);
0652   return;
0653 
0654 }  // end 'FillCalHistograms(int, int, Hist::CalContent&)'
0655 
0656 
0657 
0658 void CheckClusterSplittingProcessor::SaveHistograms() {
0659 
0660   // make directory and save histograms
0661   for (size_t iCalo = 0; iCalo < m_config.vecCalos.size(); iCalo++) {
0662 
0663     // save event histograms
0664     for (size_t iEType = 0; iEType < m_vecEvtH1D[iCalo].size(); iEType++) {
0665       for (auto evt1D : m_vecEvtH1D[iCalo][iEType]) {
0666         evt1D -> Write();
0667       }
0668     }
0669 
0670     // save track histograms
0671     for (size_t iTType = 0; iTType < m_vecTrkH1D[iCalo].size(); iTType++) {
0672       for (auto trk1D : m_vecTrkH1D[iCalo][iTType]) {
0673         trk1D -> Write();
0674       }
0675       for (auto trk2D : m_vecTrkH2D[iCalo][iTType]) {
0676         trk2D -> Write();
0677       }
0678     }
0679 
0680     // save cluster histograms
0681     for (size_t iCType = 0; iCType < m_vecTrkH1D[iCalo].size(); iCType++) {
0682       for (auto cal1D : m_vecCalH1D[iCalo][iCType]) {
0683         cal1D -> Write();
0684       }
0685       for (auto cal2D : m_vecCalH2D[iCalo][iCType]) {
0686         cal2D -> Write();
0687       }
0688     }
0689   }  // end calo loop
0690   return;
0691 
0692 }  // end 'SaveHistograms()'
0693 
0694 
0695 
0696 void CheckClusterSplittingProcessor::TryToRecover(const double eStart, const double ePar, bool& didRecover, float& drRecover, int32_t& nRecover) {
0697 
0698   // get no. of clusters to check
0699   const uint32_t nToCheck = m_mapClustToDr.size();
0700 
0701   // iterate through clusters to check until recovered or nothing left to check
0702   bool     checkedAll = false;
0703   bool     recovered  = false;
0704   double   drChecked  = 0.;
0705   double   eChecked   = eStart;
0706   uint32_t nChecked   = 1;
0707   uint32_t nStep      = 1;
0708   while (!recovered && !checkedAll) {
0709 
0710     // get new radius to check in,
0711     //   and exit if larger than max
0712     const double drToCheck = nStep * m_config.drStep;
0713     if (drToCheck > m_config.drMax) {
0714       break;
0715     }
0716 
0717     // loop through list
0718     for (auto idAndDr : m_mapClustToDr) {
0719 
0720       // skip if already checked
0721       if (m_mapClustToChecked[idAndDr.first]) {
0722         continue;
0723       }
0724 
0725       // skip if not in radius to check
0726       if (idAndDr.second > drToCheck) {
0727         continue;
0728       }
0729 
0730      // if in raidus add to running energy sum
0731      eChecked += m_mapClustToEne[idAndDr.first];
0732 
0733      // and flag as used, increment counters
0734      m_mapClustToChecked[idAndDr.first] = true;
0735      ++nChecked;
0736     }  // end to-check loop
0737 
0738     // now check if recovered 90%
0739     const double newPurity = eChecked / ePar;
0740     if (newPurity >= 0.9) {
0741       recovered = true;
0742     }
0743 
0744     // and check if all clusters have been checked
0745     if (nChecked == nToCheck) {
0746       checkedAll = true;
0747     }
0748 
0749     // update coounters
0750     drChecked = drToCheck;
0751     ++nStep;
0752 
0753   }  // end while not recovered and not checked all clusters
0754 
0755   // if recover successful, update otuputs
0756   if (recovered) {
0757     didRecover = recovered;
0758     drRecover  = drChecked;
0759     nRecover   = nChecked;
0760   }
0761   return;
0762 
0763 }  // end 'TryToRecover(double, double, bool&, float&, int32_t&)'
0764 
0765 
0766 
0767 void CheckClusterSplittingProcessor::GetClustersAndDr(const edm4eic::ClusterCollection& others, const Hist::CalContent& reference, const int idRef) {
0768 
0769   // loop through clusters
0770   m_mapClustToChecked.clear();
0771   m_mapClustToDr.clear();
0772   m_mapClustToEne.clear();
0773   for (auto other : others) {
0774 
0775     // skip this cluster
0776     const int id = other.getObjectID().index;
0777     if (id == idRef) {
0778       continue;
0779     }
0780 
0781     // grab cluster info
0782     Hist::CalContent cOther;
0783     cOther.GetInfo(other);
0784 
0785     // get distances
0786     const double dEta = (cOther.eta - reference.eta);
0787     const double dPhi = (cOther.phi - reference.phi);
0788     const double dR   = hypot(dEta, dPhi);
0789 
0790     // add to lists
0791     m_mapClustToChecked.insert( {id, false} );
0792     m_mapClustToEne.insert( {id, cOther.ene} );
0793     m_mapClustToDr.insert( {id, dR} );
0794 
0795   }  // end cluster loop
0796   return;
0797 
0798 }  // end 'GetClustersAndDr(edm4eic::ClusterCollection&, Hist::CalContent&, int)'
0799 
0800 
0801 
0802 void CheckClusterSplittingProcessor::GetProjections(const edm4eic::TrackSegmentCollection& projections, const int calo) {
0803 
0804   // make sure vector is empty
0805   m_vecTrkProj.clear();
0806 
0807   // collect projections
0808   for (auto project : projections) {
0809     for (auto point : project.getPoints()) {
0810       const bool isInSystem  = (point.system  == calo);
0811       const bool isAtFace    = (point.surface == 1);
0812       if (isInSystem && isAtFace) {
0813         m_vecTrkProj.push_back(point);
0814       }
0815     }
0816   }
0817   return;
0818 
0819 }  // end 'GetProjections(edm4eic::TrackSegmentCollection&)'
0820 
0821 
0822 
0823 int CheckClusterSplittingProcessor::GetCaloID(const int iCalo) {
0824 
0825   // get detector element
0826   auto detector = GetApplication() -> GetService<DD4hep_service>() -> detector();
0827   auto element  = detector -> detector(m_config.vecCalos.at(iCalo).second);
0828   return element.id();
0829 
0830 }  // end 'GetCaloID(int)'
0831 
0832 
0833 
0834 FPair CheckClusterSplittingProcessor::GetClusterWidths(const edm4eic::Cluster& cluster) {
0835 
0836   // get eta, phi for cluster
0837   const float hClust = edm4hep::utils::eta( cluster.getPosition() );
0838   const float fClust = atan2( cluster.getPosition().y, cluster.getPosition().x );
0839 
0840   float hWidth2 = 0.;
0841   float fWidth2 = 0.;
0842   for (auto hit : cluster.getHits()) {
0843 
0844     // get eta, phi for hit
0845     const float hHit = edm4hep::utils::eta( hit.getPosition() );
0846     const float fHit = atan2( hit.getPosition().y, hit.getPosition().x );
0847 
0848     // increment sum of squares
0849     hWidth2 += pow(hHit - hClust, 2.);
0850     fWidth2 += pow(fHit - fClust, 2.);
0851   }
0852 
0853   // calculate widths
0854   const float hWidth = sqrt(hWidth2 / cluster.getHits().size());
0855   const float fWidth = sqrt(fWidth2 / cluster.getHits().size());
0856 
0857   // if nHit = 1, return some min size
0858   //   FIXME this should be the cell dimensions
0859   FPair widths  = make_pair(hWidth, fWidth);
0860   FPair minimum = make_pair(0.05,   0.05);
0861   return (cluster.getHits().size() > 1) ? widths : minimum;
0862 
0863 }  // end 'GetClusterWidths(edm4eic::Cluster&)'
0864 
0865 
0866 
0867 FPair CheckClusterSplittingProcessor::GetEtaRange(const edm4eic::Cluster& cluster) {
0868 
0869   // get hits
0870   auto hits = cluster.getHits();
0871 
0872   // get hit with min eta
0873   auto minEtaHit = min_element(
0874     hits.begin(),
0875     hits.end(),
0876     [](const auto minEtaHit1, const auto minEtaHit2) {
0877       return ( edm4hep::utils::eta(minEtaHit1.getPosition()) < edm4hep::utils::eta(minEtaHit2.getPosition()) );
0878     }
0879   );
0880 
0881   // get hit with max eta
0882   auto maxEtaHit = max_element(
0883     hits.begin(),
0884     hits.end(),
0885     [](const auto maxEtaHit1, const auto maxEtaHit2) {
0886       return ( edm4hep::utils::eta(maxEtaHit1.getPosition()) < edm4hep::utils::eta(maxEtaHit2.getPosition()) );
0887     }
0888   );
0889   return make_pair(
0890     edm4hep::utils::eta(minEtaHit -> getPosition()),
0891     edm4hep::utils::eta(maxEtaHit -> getPosition())
0892   );
0893 
0894 }  // 'GetEtaRange(edm4eic::Cluster& cluster)'
0895 
0896 
0897 
0898 FPair CheckClusterSplittingProcessor::GetPhiRange(const edm4eic::Cluster& cluster) {
0899 
0900   // get hits
0901   auto hits = cluster.getHits();
0902 
0903   // get hit with min phi
0904   auto minPhiHit = min_element(
0905     hits.begin(),
0906     hits.end(),
0907     [](const auto minPhiHit1, const auto minPhiHit2) {
0908       return ( atan2(minPhiHit1.getPosition().y, minPhiHit1.getPosition().x) < atan2(minPhiHit2.getPosition().y, minPhiHit2.getPosition().x) );
0909     }
0910   );
0911 
0912   // get hit with max phi
0913   auto maxPhiHit = max_element(
0914     hits.begin(),
0915     hits.end(),
0916     [](const auto maxPhiHit1, const auto maxPhiHit2) {
0917       return ( atan2(maxPhiHit1.getPosition().y, maxPhiHit1.getPosition().x) < atan2(maxPhiHit2.getPosition().y, maxPhiHit2.getPosition().x) );
0918     }
0919   );
0920   return make_pair(
0921     atan2(minPhiHit -> getPosition().y, minPhiHit -> getPosition().x),
0922     atan2(maxPhiHit -> getPosition().y, maxPhiHit -> getPosition().x)
0923   );
0924 
0925 }  // 'GetPhiRange(edm4eic::Cluster& cluster)'
0926 
0927 // end ------------------------------------------------------------------------
[ 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!