Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-12-15 09:25:36 

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008 
0009 #pragma once
0010 
0011 #include <array>
0012 #include <string>
0013 #include <vector>
0014 
0015 #include <TCanvas.h>
0016 #include <TEfficiency.h>
0017 #include <TFile.h>
0018 #include <TTree.h>
0019 
0020 struct RecoTrackInfo {
0021   double eta;
0022   double pt;
0023   unsigned int nMajorityHits;
0024   unsigned int nMeasurements;
0025 };
0026 
0027 struct ParticleInfo {
0028   ULong64_t particleId = 0;
0029   double eta = 0;
0030   double p = 0;
0031   double pt = 0;
0032   UShort_t nHits = 0;
0033 };
0034 
0035 namespace {
0036 
0037 inline std::uint64_t hashBarcodeComponents(std::uint32_t vertexPrimary,
0038                                            std::uint32_t vertexSecondary,
0039                                            std::uint32_t particle,
0040                                            std::uint32_t generation,
0041                                            std::uint32_t subParticle) {
0042   auto hashCombine = [](std::uint64_t seed, std::uint64_t value) {
0043     constexpr std::uint64_t kMagic = 0x9e3779b97f4a7c15ull;
0044     seed ^= value + kMagic + (seed << 6) + (seed >> 2);
0045     return seed;
0046   };
0047 
0048   std::uint64_t hash = 0xcbf29ce484222325ull;
0049   hash = hashCombine(hash, vertexPrimary);
0050   hash = hashCombine(hash, vertexSecondary);
0051   hash = hashCombine(hash, particle);
0052   hash = hashCombine(hash, generation);
0053   hash = hashCombine(hash, subParticle);
0054   return hash;
0055 }
0056 
0057 }  // namespace
0058 
0059 /// Helper for reading tree
0060 ///
0061 struct TreeReader {
0062   // The constructor
0063   explicit TreeReader(TTree* tree_) : tree(tree_) {}
0064 
0065   // Get entry
0066   void getEntry(unsigned int i) const {
0067     if (entryNumbers.size() > i) {
0068       tree->GetEntry(entryNumbers[i]);
0069     } else {
0070       tree->GetEntry(i);
0071     }
0072   };
0073 
0074   // The tree being read
0075   TTree* tree = nullptr;
0076 
0077  protected:
0078   /// The entry numbers for accessing events in increased order (there could be
0079   /// multiple entries corresponding to one event number)
0080   std::vector<long long> entryNumbers = {};
0081 };
0082 
0083 /// Struct used for reading track states written out by the
0084 /// RootTrackStatesWriter
0085 ///
0086 struct TrackStatesReader : public TreeReader {
0087   // Delete the default constructor
0088   TrackStatesReader() = delete;
0089 
0090   // The constructor
0091   TrackStatesReader(TTree* tree_, bool sortEvents) : TreeReader(tree_) {
0092     tree->SetBranchAddress("event_nr", &eventId);
0093     tree->SetBranchAddress("eLOC0_prt", &LOC0_prt);
0094     tree->SetBranchAddress("eLOC1_prt", &LOC1_prt);
0095     tree->SetBranchAddress("ePHI_prt", &PHI_prt);
0096     tree->SetBranchAddress("eTHETA_prt", &THETA_prt);
0097     tree->SetBranchAddress("eQOP_prt", &QOP_prt);
0098     tree->SetBranchAddress("eT_prt", &T_prt);
0099     tree->SetBranchAddress("eLOC0_flt", &LOC0_flt);
0100     tree->SetBranchAddress("eLOC1_flt", &LOC1_flt);
0101     tree->SetBranchAddress("ePHI_flt", &PHI_flt);
0102     tree->SetBranchAddress("eTHETA_flt", &THETA_flt);
0103     tree->SetBranchAddress("eQOP_flt", &QOP_flt);
0104     tree->SetBranchAddress("eT_flt", &T_flt);
0105     tree->SetBranchAddress("eLOC0_smt", &LOC0_smt);
0106     tree->SetBranchAddress("eLOC1_smt", &LOC1_smt);
0107     tree->SetBranchAddress("ePHI_smt", &PHI_smt);
0108     tree->SetBranchAddress("eTHETA_smt", &THETA_smt);
0109     tree->SetBranchAddress("eQOP_smt", &QOP_smt);
0110     tree->SetBranchAddress("eT_smt", &T_smt);
0111 
0112     tree->SetBranchAddress("res_eLOC0_prt", &res_LOC0_prt);
0113     tree->SetBranchAddress("res_eLOC1_prt", &res_LOC1_prt);
0114     tree->SetBranchAddress("res_ePHI_prt", &res_PHI_prt);
0115     tree->SetBranchAddress("res_eTHETA_prt", &res_THETA_prt);
0116     tree->SetBranchAddress("res_eQOP_prt", &res_QOP_prt);
0117     tree->SetBranchAddress("res_eT_prt", &res_T_prt);
0118     tree->SetBranchAddress("res_eLOC0_flt", &res_LOC0_flt);
0119     tree->SetBranchAddress("res_eLOC1_flt", &res_LOC1_flt);
0120     tree->SetBranchAddress("res_ePHI_flt", &res_PHI_flt);
0121     tree->SetBranchAddress("res_eTHETA_flt", &res_THETA_flt);
0122     tree->SetBranchAddress("res_eQOP_flt", &res_QOP_flt);
0123     tree->SetBranchAddress("res_eT_flt", &res_T_flt);
0124     tree->SetBranchAddress("res_eLOC0_smt", &res_LOC0_smt);
0125     tree->SetBranchAddress("res_eLOC1_smt", &res_LOC1_smt);
0126     tree->SetBranchAddress("res_ePHI_smt", &res_PHI_smt);
0127     tree->SetBranchAddress("res_eTHETA_smt", &res_THETA_smt);
0128     tree->SetBranchAddress("res_eQOP_smt", &res_QOP_smt);
0129     tree->SetBranchAddress("res_eT_smt", &res_T_smt);
0130 
0131     tree->SetBranchAddress("pull_eLOC0_prt", &pull_LOC0_prt);
0132     tree->SetBranchAddress("pull_eLOC1_prt", &pull_LOC1_prt);
0133     tree->SetBranchAddress("pull_ePHI_prt", &pull_PHI_prt);
0134     tree->SetBranchAddress("pull_eTHETA_prt", &pull_THETA_prt);
0135     tree->SetBranchAddress("pull_eQOP_prt", &pull_QOP_prt);
0136     tree->SetBranchAddress("pull_eT_prt", &pull_T_prt);
0137     tree->SetBranchAddress("pull_eLOC0_flt", &pull_LOC0_flt);
0138     tree->SetBranchAddress("pull_eLOC1_flt", &pull_LOC1_flt);
0139     tree->SetBranchAddress("pull_ePHI_flt", &pull_PHI_flt);
0140     tree->SetBranchAddress("pull_eTHETA_flt", &pull_THETA_flt);
0141     tree->SetBranchAddress("pull_eQOP_flt", &pull_QOP_flt);
0142     tree->SetBranchAddress("pull_eT_flt", &pull_T_flt);
0143     tree->SetBranchAddress("pull_eLOC0_smt", &pull_LOC0_smt);
0144     tree->SetBranchAddress("pull_eLOC1_smt", &pull_LOC1_smt);
0145     tree->SetBranchAddress("pull_ePHI_smt", &pull_PHI_smt);
0146     tree->SetBranchAddress("pull_eTHETA_smt", &pull_THETA_smt);
0147     tree->SetBranchAddress("pull_eQOP_smt", &pull_QOP_smt);
0148     tree->SetBranchAddress("pull_eT_smt", &pull_T_smt);
0149 
0150     tree->SetBranchAddress("g_x_prt", &g_x_prt);
0151     tree->SetBranchAddress("g_y_prt", &g_y_prt);
0152     tree->SetBranchAddress("g_z_prt", &g_z_prt);
0153     tree->SetBranchAddress("g_x_flt", &g_x_flt);
0154     tree->SetBranchAddress("g_y_flt", &g_y_flt);
0155     tree->SetBranchAddress("g_z_flt", &g_z_flt);
0156     tree->SetBranchAddress("g_x_smt", &g_x_smt);
0157     tree->SetBranchAddress("g_y_smt", &g_y_smt);
0158     tree->SetBranchAddress("g_z_smt", &g_z_smt);
0159 
0160     tree->SetBranchAddress("nStates", &nStates);
0161     tree->SetBranchAddress("nMeasurements", &nMeasurements);
0162     tree->SetBranchAddress("volume_id", &volume_id);
0163     tree->SetBranchAddress("layer_id", &layer_id);
0164     tree->SetBranchAddress("module_id", &module_id);
0165     tree->SetBranchAddress("predicted", &predicted);
0166     tree->SetBranchAddress("filtered", &filtered);
0167     tree->SetBranchAddress("smoothed", &smoothed);
0168 
0169     // It's not necessary if you just need to read one file, but please do it to
0170     // synchronize events if multiple root files are read
0171     if (sortEvents) {
0172       tree->SetEstimate(tree->GetEntries() + 1);
0173       entryNumbers.resize(tree->GetEntries());
0174       tree->Draw("event_nr", "", "goff");
0175       // Sort to get the entry numbers of the ordered events
0176       TMath::Sort(tree->GetEntries(), tree->GetV1(), entryNumbers.data(),
0177                   false);
0178     }
0179   }
0180 
0181   // The variables
0182   std::uint32_t eventId = 0;
0183   std::vector<float>* LOC0_prt =
0184       new std::vector<float>;  ///< predicted parameter local x
0185   std::vector<float>* LOC1_prt =
0186       new std::vector<float>;  ///< predicted parameter local y
0187   std::vector<float>* PHI_prt =
0188       new std::vector<float>;  ///< predicted parameter phi
0189   std::vector<float>* THETA_prt =
0190       new std::vector<float>;  ///< predicted parameter theta
0191   std::vector<float>* QOP_prt =
0192       new std::vector<float>;  ///< predicted parameter q/p
0193   std::vector<float>* T_prt =
0194       new std::vector<float>;  ///< predicted parameter t
0195   std::vector<float>* LOC0_flt =
0196       new std::vector<float>;  ///< filtered parameter local x
0197   std::vector<float>* LOC1_flt =
0198       new std::vector<float>;  ///< filtered parameter local y
0199   std::vector<float>* PHI_flt =
0200       new std::vector<float>;  ///< filtered parameter phi
0201   std::vector<float>* THETA_flt =
0202       new std::vector<float>;  ///< filtered parameter theta
0203   std::vector<float>* QOP_flt =
0204       new std::vector<float>;  ///< filtered parameter q/p
0205   std::vector<float>* T_flt = new std::vector<float>;  ///< filtered parameter t
0206   std::vector<float>* LOC0_smt =
0207       new std::vector<float>;  ///< smoothed parameter local x
0208   std::vector<float>* LOC1_smt =
0209       new std::vector<float>;  ///< smoothed parameter local y
0210   std::vector<float>* PHI_smt =
0211       new std::vector<float>;  ///< smoothed parameter phi
0212   std::vector<float>* THETA_smt =
0213       new std::vector<float>;  ///< smoothed parameter theta
0214   std::vector<float>* QOP_smt =
0215       new std::vector<float>;  ///< smoothed parameter q/p
0216   std::vector<float>* T_smt = new std::vector<float>;  ///< smoothed parameter t
0217 
0218   std::vector<float>* res_LOC0_prt =
0219       new std::vector<float>;  ///< residual of predicted parameter local x
0220   std::vector<float>* res_LOC1_prt =
0221       new std::vector<float>;  ///< residual of predicted parameter local y
0222   std::vector<float>* res_PHI_prt =
0223       new std::vector<float>;  ///< residual of predicted parameter phi
0224   std::vector<float>* res_THETA_prt =
0225       new std::vector<float>;  ///< residual of predicted parameter theta
0226   std::vector<float>* res_QOP_prt =
0227       new std::vector<float>;  ///< residual of predicted parameter q/p
0228   std::vector<float>* res_T_prt =
0229       new std::vector<float>;  ///< residual of predicted parameter t
0230   std::vector<float>* res_LOC0_flt =
0231       new std::vector<float>;  ///< residual of filtered parameter local x
0232   std::vector<float>* res_LOC1_flt =
0233       new std::vector<float>;  ///< residual of filtered parameter local y
0234   std::vector<float>* res_PHI_flt =
0235       new std::vector<float>;  ///< residual of filtered parameter phi
0236   std::vector<float>* res_THETA_flt =
0237       new std::vector<float>;  ///< residual of filtered parameter theta
0238   std::vector<float>* res_QOP_flt =
0239       new std::vector<float>;  ///< residual of filtered parameter q/p
0240   std::vector<float>* res_T_flt =
0241       new std::vector<float>;  ///< residual of filtered parameter t
0242   std::vector<float>* res_LOC0_smt =
0243       new std::vector<float>;  ///< residual of smoothed parameter local x
0244   std::vector<float>* res_LOC1_smt =
0245       new std::vector<float>;  ///< residual of smoothed parameter local y
0246   std::vector<float>* res_PHI_smt =
0247       new std::vector<float>;  ///< residual of smoothed parameter phi
0248   std::vector<float>* res_THETA_smt =
0249       new std::vector<float>;  ///< residual of smoothed parameter theta
0250   std::vector<float>* res_QOP_smt =
0251       new std::vector<float>;  ///< residual of smoothed parameter q/p
0252   std::vector<float>* res_T_smt =
0253       new std::vector<float>;  ///< residual of smoothed parameter t
0254 
0255   std::vector<float>* pull_LOC0_prt =
0256       new std::vector<float>;  ///< pull of predicted parameter local x
0257   std::vector<float>* pull_LOC1_prt =
0258       new std::vector<float>;  ///< pull of predicted parameter local y
0259   std::vector<float>* pull_PHI_prt =
0260       new std::vector<float>;  ///< pull of predicted parameter phi
0261   std::vector<float>* pull_THETA_prt =
0262       new std::vector<float>;  ///< pull of predicted parameter theta
0263   std::vector<float>* pull_QOP_prt =
0264       new std::vector<float>;  ///< pull of predicted parameter q/p
0265   std::vector<float>* pull_T_prt =
0266       new std::vector<float>;  ///< pull of predicted parameter t
0267   std::vector<float>* pull_LOC0_flt =
0268       new std::vector<float>;  ///< pull of filtered parameter local x
0269   std::vector<float>* pull_LOC1_flt =
0270       new std::vector<float>;  ///< pull of filtered parameter local y
0271   std::vector<float>* pull_PHI_flt =
0272       new std::vector<float>;  ///< pull of filtered parameter phi
0273   std::vector<float>* pull_THETA_flt =
0274       new std::vector<float>;  ///< pull of filtered parameter theta
0275   std::vector<float>* pull_QOP_flt =
0276       new std::vector<float>;  ///< pull of filtered parameter q/p
0277   std::vector<float>* pull_T_flt =
0278       new std::vector<float>;  ///< pull of filtered parameter t
0279   std::vector<float>* pull_LOC0_smt =
0280       new std::vector<float>;  ///< pull of smoothed parameter local x
0281   std::vector<float>* pull_LOC1_smt =
0282       new std::vector<float>;  ///< pull of smoothed parameter local y
0283   std::vector<float>* pull_PHI_smt =
0284       new std::vector<float>;  ///< pull of smoothed parameter phi
0285   std::vector<float>* pull_THETA_smt =
0286       new std::vector<float>;  ///< pull of smoothed parameter theta
0287   std::vector<float>* pull_QOP_smt =
0288       new std::vector<float>;  ///< pull of smoothed parameter q/p
0289   std::vector<float>* pull_T_smt =
0290       new std::vector<float>;  ///< pull of smoothed parameter t
0291 
0292   std::vector<float>* g_x_prt = new std::vector<float>;
0293   std::vector<float>* g_y_prt = new std::vector<float>;
0294   std::vector<float>* g_z_prt = new std::vector<float>;
0295   std::vector<float>* g_x_flt = new std::vector<float>;
0296   std::vector<float>* g_y_flt = new std::vector<float>;
0297   std::vector<float>* g_z_flt = new std::vector<float>;
0298   std::vector<float>* g_x_smt = new std::vector<float>;
0299   std::vector<float>* g_y_smt = new std::vector<float>;
0300   std::vector<float>* g_z_smt = new std::vector<float>;
0301 
0302   std::vector<int>* volume_id = new std::vector<int>;  ///< volume_id
0303   std::vector<int>* layer_id = new std::vector<int>;   ///< layer_id
0304   std::vector<int>* module_id = new std::vector<int>;  ///< module_id
0305 
0306   std::vector<bool>* predicted = new std::vector<bool>;  ///< prediction status
0307   std::vector<bool>* filtered = new std::vector<bool>;   ///< filtering status
0308   std::vector<bool>* smoothed = new std::vector<bool>;   ///< smoothing status
0309 
0310   unsigned int nStates = 0, nMeasurements = 0;
0311 };
0312 
0313 /// Struct used for reading track summary info written out by the
0314 /// RootTrackSummaryWriter
0315 ///
0316 struct TrackSummaryReader : public TreeReader {
0317   // Delete the default constructor
0318   TrackSummaryReader() = delete;
0319 
0320   // The constructor
0321   TrackSummaryReader(TTree* tree_, bool sortEvents) : TreeReader(tree_) {
0322     tree->SetBranchAddress("event_nr", &eventId);
0323     tree->SetBranchAddress("nStates", &nStates);
0324     tree->SetBranchAddress("nMeasurements", &nMeasurements);
0325     tree->SetBranchAddress("nOutliers", &nOutliers);
0326     tree->SetBranchAddress("nHoles", &nHoles);
0327     tree->SetBranchAddress("chi2Sum", &chi2Sum);
0328     tree->SetBranchAddress("measurementChi2", &measurementChi2);
0329     tree->SetBranchAddress("NDF", &NDF);
0330     tree->SetBranchAddress("measurementVolume", &measurementVolume);
0331     tree->SetBranchAddress("measurementLayer", &measurementLayer);
0332     tree->SetBranchAddress("outlierVolume", &outlierVolume);
0333     tree->SetBranchAddress("outlierLayer", &outlierLayer);
0334     tree->SetBranchAddress("nMajorityHits", &nMajorityHits);
0335     tree->SetBranchAddress("nSharedHits", &nSharedHits);
0336     if (tree->GetBranch("majorityParticleId") != nullptr) {
0337       majorityParticleId =
0338           new std::vector<std::vector<std::uint32_t>>;
0339       tree->SetBranchAddress("majorityParticleId", &majorityParticleId);
0340       hasCombinedMajorityParticleId = true;
0341     } else {
0342       majorityParticleVertexPrimary = new std::vector<std::uint32_t>;
0343       majorityParticleVertexSecondary = new std::vector<std::uint32_t>;
0344       majorityParticleParticle = new std::vector<std::uint32_t>;
0345       majorityParticleGeneration = new std::vector<std::uint32_t>;
0346       majorityParticleSubParticle = new std::vector<std::uint32_t>;
0347       tree->SetBranchAddress("majorityParticleId_vertex_primary",
0348                              &majorityParticleVertexPrimary);
0349       tree->SetBranchAddress("majorityParticleId_vertex_secondary",
0350                              &majorityParticleVertexSecondary);
0351       tree->SetBranchAddress("majorityParticleId_particle",
0352                              &majorityParticleParticle);
0353       tree->SetBranchAddress("majorityParticleId_generation",
0354                              &majorityParticleGeneration);
0355       tree->SetBranchAddress("majorityParticleId_sub_particle",
0356                              &majorityParticleSubParticle);
0357     }
0358 
0359     tree->SetBranchAddress("hasFittedParams", &hasFittedParams);
0360 
0361     tree->SetBranchAddress("t_theta", &t_theta);
0362     tree->SetBranchAddress("t_phi", &t_phi);
0363     tree->SetBranchAddress("t_eta", &t_eta);
0364     tree->SetBranchAddress("t_p", &t_p);
0365     tree->SetBranchAddress("t_pT", &t_pT);
0366     tree->SetBranchAddress("t_d0", &t_d0);
0367     tree->SetBranchAddress("t_z0", &t_z0);
0368     tree->SetBranchAddress("t_charge", &t_charge);
0369     tree->SetBranchAddress("t_time", &t_time);
0370 
0371     tree->SetBranchAddress("eLOC0_fit", &eLOC0_fit);
0372     tree->SetBranchAddress("eLOC1_fit", &eLOC1_fit);
0373     tree->SetBranchAddress("ePHI_fit", &ePHI_fit);
0374     tree->SetBranchAddress("eTHETA_fit", &eTHETA_fit);
0375     tree->SetBranchAddress("eQOP_fit", &eQOP_fit);
0376     tree->SetBranchAddress("eT_fit", &eT_fit);
0377     tree->SetBranchAddress("err_eLOC0_fit", &err_eLOC0_fit);
0378     tree->SetBranchAddress("err_eLOC1_fit", &err_eLOC1_fit);
0379     tree->SetBranchAddress("err_ePHI_fit", &err_ePHI_fit);
0380     tree->SetBranchAddress("err_eTHETA_fit", &err_eTHETA_fit);
0381     tree->SetBranchAddress("err_eQOP_fit", &err_eQOP_fit);
0382     tree->SetBranchAddress("err_eT_fit", &err_eT_fit);
0383 
0384     // It's not necessary if you just need to read one file, but please do it to
0385     // synchronize events if multiple root files are read
0386     if (sortEvents) {
0387       tree->SetEstimate(tree->GetEntries() + 1);
0388       entryNumbers.resize(tree->GetEntries());
0389       tree->Draw("event_nr", "", "goff");
0390       // Sort to get the entry numbers of the ordered events
0391       TMath::Sort(tree->GetEntries(), tree->GetV1(), entryNumbers.data(),
0392                   false);
0393     }
0394   }
0395 
0396   // The variables
0397   std::uint32_t eventId = 0;
0398   std::vector<unsigned int>* nStates = new std::vector<unsigned int>;
0399   std::vector<unsigned int>* nMeasurements = new std::vector<unsigned int>;
0400   std::vector<unsigned int>* nOutliers = new std::vector<unsigned int>;
0401   std::vector<unsigned int>* nHoles = new std::vector<unsigned int>;
0402   std::vector<unsigned int>* nSharedHits = new std::vector<unsigned int>;
0403   std::vector<float>* chi2Sum = new std::vector<float>;
0404   std::vector<unsigned int>* NDF = new std::vector<unsigned int>;
0405   std::vector<std::vector<double>>* measurementChi2 =
0406       new std::vector<std::vector<double>>;
0407   std::vector<std::vector<double>>* outlierChi2 =
0408       new std::vector<std::vector<double>>;
0409   std::vector<std::vector<double>>* measurementVolume =
0410       new std::vector<std::vector<double>>;
0411   std::vector<std::vector<double>>* measurementLayer =
0412       new std::vector<std::vector<double>>;
0413   std::vector<std::vector<double>>* outlierVolume =
0414       new std::vector<std::vector<double>>;
0415   std::vector<std::vector<double>>* outlierLayer =
0416       new std::vector<std::vector<double>>;
0417   std::vector<unsigned int>* nMajorityHits = new std::vector<unsigned int>;
0418   std::vector<std::vector<std::uint32_t>>* majorityParticleId = nullptr;
0419   std::vector<std::uint32_t>* majorityParticleVertexPrimary = nullptr;
0420   std::vector<std::uint32_t>* majorityParticleVertexSecondary = nullptr;
0421   std::vector<std::uint32_t>* majorityParticleParticle = nullptr;
0422   std::vector<std::uint32_t>* majorityParticleGeneration = nullptr;
0423   std::vector<std::uint32_t>* majorityParticleSubParticle = nullptr;
0424   bool hasCombinedMajorityParticleId = false;
0425 
0426   std::uint64_t majorityParticleHash(std::size_t idx) const {
0427     if (hasCombinedMajorityParticleId && majorityParticleId != nullptr) {
0428       const auto& components = majorityParticleId->at(idx);
0429       auto comp = [&](std::size_t i) -> std::uint32_t {
0430         return (components.size() > i) ? components[i] : 0u;
0431       };
0432       return hashBarcodeComponents(comp(0), comp(1), comp(2), comp(3), comp(4));
0433     }
0434     auto safeAt = [](const std::vector<std::uint32_t>* vec, std::size_t i) {
0435       return (vec != nullptr && vec->size() > i) ? vec->at(i) : 0u;
0436     };
0437     return hashBarcodeComponents(safeAt(majorityParticleVertexPrimary, idx),
0438                                  safeAt(majorityParticleVertexSecondary, idx),
0439                                  safeAt(majorityParticleParticle, idx),
0440                                  safeAt(majorityParticleGeneration, idx),
0441                                  safeAt(majorityParticleSubParticle, idx));
0442   }
0443 
0444   std::vector<bool>* hasFittedParams = new std::vector<bool>;
0445 
0446   // True parameters
0447   std::vector<float>* t_d0 = new std::vector<float>;
0448   std::vector<float>* t_z0 = new std::vector<float>;
0449   std::vector<float>* t_phi = new std::vector<float>;
0450   std::vector<float>* t_theta = new std::vector<float>;
0451   std::vector<float>* t_eta = new std::vector<float>;
0452   std::vector<float>* t_p = new std::vector<float>;
0453   std::vector<float>* t_pT = new std::vector<float>;
0454   std::vector<float>* t_time = new std::vector<float>;
0455   std::vector<int>* t_charge = new std::vector<int>;
0456 
0457   // Estimated parameters
0458   std::vector<float>* eLOC0_fit = new std::vector<float>;
0459   std::vector<float>* eLOC1_fit = new std::vector<float>;
0460   std::vector<float>* ePHI_fit = new std::vector<float>;
0461   std::vector<float>* eTHETA_fit = new std::vector<float>;
0462   std::vector<float>* eQOP_fit = new std::vector<float>;
0463   std::vector<float>* eT_fit = new std::vector<float>;
0464 
0465   std::vector<float>* err_eLOC0_fit = new std::vector<float>;
0466   std::vector<float>* err_eLOC1_fit = new std::vector<float>;
0467   std::vector<float>* err_ePHI_fit = new std::vector<float>;
0468   std::vector<float>* err_eTHETA_fit = new std::vector<float>;
0469   std::vector<float>* err_eQOP_fit = new std::vector<float>;
0470   std::vector<float>* err_eT_fit = new std::vector<float>;
0471 };
0472 
0473 /// Struct used for reading particles written out by the
0474 /// RootTrackFinderNTupleWriter
0475 ///
0476 struct ParticleReader : public TreeReader {
0477   // Delete the default constructor
0478   ParticleReader() = delete;
0479 
0480   // The constructor
0481   ParticleReader(TTree* tree_, bool sortEvents) : TreeReader(tree_) {
0482     tree->SetBranchAddress("event_id", &eventId);
0483     if (tree->GetBranch("particle_id") != nullptr) {
0484       combinedParticleId = new std::vector<std::uint32_t>;
0485       tree->SetBranchAddress("particle_id", &combinedParticleId);
0486       hasCombinedParticleId = true;
0487     } else {
0488       tree->SetBranchAddress("particle_id_vertex_primary",
0489                              &particleVertexPrimary);
0490       tree->SetBranchAddress("particle_id_vertex_secondary",
0491                              &particleVertexSecondary);
0492       tree->SetBranchAddress("particle_id_particle", &particleParticle);
0493       tree->SetBranchAddress("particle_id_generation", &particleGeneration);
0494       tree->SetBranchAddress("particle_id_sub_particle",
0495                              &particleSubParticle);
0496     }
0497     tree->SetBranchAddress("particle_type", &particleType);
0498     tree->SetBranchAddress("vx", &vx);
0499     tree->SetBranchAddress("vy", &vy);
0500     tree->SetBranchAddress("vz", &vz);
0501     tree->SetBranchAddress("vt", &vt);
0502     tree->SetBranchAddress("px", &px);
0503     tree->SetBranchAddress("py", &py);
0504     tree->SetBranchAddress("pz", &pz);
0505     tree->SetBranchAddress("m", &m);
0506     tree->SetBranchAddress("q", &q);
0507     tree->SetBranchAddress("nhits", &nHits);
0508     tree->SetBranchAddress("ntracks", &nTracks);
0509     tree->SetBranchAddress("ntracks_majority", &nTracksMajority);
0510 
0511     // It's not necessary if you just need to read one file, but please do it to
0512     // synchronize events if multiple root files are read
0513     if (sortEvents) {
0514       tree->SetEstimate(tree->GetEntries() + 1);
0515       entryNumbers.resize(tree->GetEntries());
0516       tree->Draw("event_id", "", "goff");
0517       // Sort to get the entry numbers of the ordered events
0518       TMath::Sort(tree->GetEntries(), tree->GetV1(), entryNumbers.data(),
0519                   false);
0520     }
0521   }
0522 
0523   // Get all the particles with requested event id
0524   std::vector<ParticleInfo> getParticles(
0525       const std::uint32_t& eventNumber) const {
0526     // Find the start entry and the batch size for this event
0527     std::string eventNumberStr = std::to_string(eventNumber);
0528     std::string findStartEntry = "event_id<" + eventNumberStr;
0529     std::string findParticlesSize = "event_id==" + eventNumberStr;
0530     std::size_t startEntry = tree->GetEntries(findStartEntry.c_str());
0531     std::size_t nParticles = tree->GetEntries(findParticlesSize.c_str());
0532     if (nParticles == 0) {
0533       throw std::invalid_argument(
0534           "No particles found. Please check the input file.");
0535     }
0536     std::vector<ParticleInfo> particles;
0537     particles.reserve(nParticles);
0538     for (unsigned int i = 0; i < nParticles; ++i) {
0539       getEntry(startEntry + i);
0540       auto pt = std::hypot(px, py);
0541       auto p = std::hypot(pt, pz);
0542       auto eta = std::atanh(pz / p * 1.);
0543 
0544       std::uint32_t vp = particleVertexPrimary;
0545       std::uint32_t vs = particleVertexSecondary;
0546       std::uint32_t particle = particleParticle;
0547       std::uint32_t generation = particleGeneration;
0548       std::uint32_t subParticle = particleSubParticle;
0549       if (hasCombinedParticleId && combinedParticleId != nullptr) {
0550         auto comp = [&](std::size_t idx) -> std::uint32_t {
0551           return (combinedParticleId->size() > idx) ? combinedParticleId->at(idx)
0552                                                     : 0u;
0553         };
0554         vp = comp(0);
0555         vs = comp(1);
0556         particle = comp(2);
0557         generation = comp(3);
0558         subParticle = comp(4);
0559       }
0560 
0561       const auto barcodeHash =
0562           hashBarcodeComponents(vp, vs, particle, generation, subParticle);
0563       particles.push_back({barcodeHash, eta, p, pt, nHits});
0564     }
0565     return particles;
0566   }
0567 
0568   // The variables
0569   ULong64_t eventId = 0;
0570   std::vector<std::uint32_t>* combinedParticleId = nullptr;
0571   bool hasCombinedParticleId = false;
0572   std::uint32_t particleVertexPrimary = 0;
0573   std::uint32_t particleVertexSecondary = 0;
0574   std::uint32_t particleParticle = 0;
0575   std::uint32_t particleGeneration = 0;
0576   std::uint32_t particleSubParticle = 0;
0577   Int_t particleType = 0;
0578   float vx = 0, vy = 0, vz = 0;
0579   float vt = 0;
0580   float px = 0, py = 0, pz = 0;
0581   float m = 0;
0582   float q = 0;
0583   UShort_t nHits = 0;
0584   UShort_t nTracks = 0;
0585   UShort_t nTracksMajority = 0;
0586 };
[ 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!