packages/HFTrackEfficiency/HFTrackEfficiency.cc

0001 #include "HFTrackEfficiency.h"
0002
0003 #include <decayfinder/DecayFinderContainerBase.h>
0004
0005 #include <trackbase_historic/PHG4ParticleSvtxMap.h>
0006 #include <trackbase_historic/SvtxTrack.h>
0007 #include <trackbase_historic/SvtxTrackMap.h>
0008 #include <trackbase_historic/SvtxTrackMap_v2.h>
0009
0010 #include <g4main/PHG4Particle.h>
0011 #include <g4main/PHG4TruthInfoContainer.h>
0012 #include <g4main/PHG4VtxPoint.h>
0013
0014 #include <phhepmc/PHHepMCGenEvent.h>
0015 #include <phhepmc/PHHepMCGenEventMap.h>
0016
0017 #include <fun4all/Fun4AllReturnCodes.h>
0018
0019 #include <phool/PHCompositeNode.h>
0020 #include <phool/getClass.h>
0021
0022 #include <TDatabasePDG.h>
0023 #include <TFile.h>
0024 #include <TTree.h>
0025
0026 #include <CLHEP/Vector/LorentzVector.h>
0027 #include <CLHEP/Vector/ThreeVector.h>
0028
0029 #include <HepMC/GenEvent.h>
0030 #include <HepMC/GenParticle.h>
0031 #include <HepMC/GenVertex.h>  // for GenVertex::particle_iterator
0032 #include <HepMC/SimpleVector.h>
0033
0034 #include <algorithm>
0035 #include <cmath>
0036 #include <iomanip>
0037 #include <iostream>
0038 #include <map>
0039 #include <set>
0040
0041 //____________________________________________________________________________..
0042 HFTrackEfficiency::HFTrackEfficiency(const std::string &name)
0043   : SubsysReco(name)
0044   , m_input_track_map_node_name("SvtxTrackMap")
0045   , m_output_track_map_node_name("HFSelected")
0046   , m_outfile_name("outputHFTrackEff.root")
0047   , m_truthRecoMatchPercent(5.)
0048   , m_nDaughters(2)
0049 {
0050 }
0051
0052 //____________________________________________________________________________..
0053 int HFTrackEfficiency::Init(PHCompositeNode *topNode)
0054 {
0055   m_truthRecoMatchPercent /= 100.;
0056
0057   if (m_write_track_map)
0058   {
0059     PHNodeIterator nodeIter(topNode);
0060
0061     PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(nodeIter.findFirst("PHCompositeNode", "DST"));
0062     if (!dstNode)
0063     {
0064       dstNode = new PHCompositeNode("DST");
0065       topNode->addNode(dstNode);
0066       std::cout << "DST node added" << std::endl;
0067     }
0068
0069     outputNodeName = m_output_track_map_node_name + "_SvtxTrackMap";
0070     m_output_trackMap = new SvtxTrackMap_v2();
0071     PHIODataNode<PHObject> *outputTrackNode = new PHIODataNode<PHObject>(m_output_trackMap, outputNodeName, "PHObject");
0072     dstNode->addNode(outputTrackNode);
0073     std::cout << outputNodeName << " node added" << std::endl;
0074   }
0075
0076   return Fun4AllReturnCodes::EVENT_OK;
0077 }
0078
0079 //____________________________________________________________________________..
0080 int HFTrackEfficiency::process_event(PHCompositeNode *topNode)
0081 {
0082   PHNodeIterator nodeIter(topNode);
0083   PHCompositeNode *dstNode = static_cast<PHCompositeNode *>(nodeIter.findFirst("PHCompositeNode", "DST"));
0084   assert(dstNode);
0085
0086   std::string df_node_name = m_df_module_name + "_DecayMap";
0087   m_decayMap = findNode::getClass<DecayFinderContainerBase>(topNode, df_node_name);
0088   if (!m_decayMap)
0089   {
0090     std::cout << __FILE__ << ": Missing node " << df_node_name << std::endl;
0091     return Fun4AllReturnCodes::ABORTRUN;
0092   }
0093
0094   m_input_trackMap = findNode::getClass<SvtxTrackMap>(topNode, m_input_track_map_node_name);
0095   if (!m_input_trackMap)
0096   {
0097     std::cout << __FILE__ << ": Missing node " << m_input_track_map_node_name << std::endl;
0098     return Fun4AllReturnCodes::ABORTRUN;
0099   }
0100
0101   m_geneventmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
0102   if (!m_geneventmap)
0103   {
0104     std::cout << __FILE__ << ": Missing node PHHepMCGenEventMap" << std::endl;
0105   }
0106
0107   m_truthInfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0108   if (!m_truthInfo)
0109   {
0110     if (Verbosity() >= VERBOSITY_MORE)
0111     {
0112       std::cout << __FILE__ << ": Missing node G4TruthInfo" << std::endl;
0113     }
0114   }
0115
0116   m_dst_truth_reco_map = findNode::getClass<PHG4ParticleSvtxMap>(topNode, "PHG4ParticleSvtxMap");
0117   if (m_dst_truth_reco_map)
0118   {
0119     if (Verbosity() >= VERBOSITY_MORE)
0120     {
0121       std::cout << __FILE__ << ": PHG4ParticleSvtxMap found, truth matching will be more accurate" << std::endl;
0122     }
0123   }
0124   else
0125   {
0126     if (Verbosity() >= VERBOSITY_MORE)
0127     {
0128       std::cout << __FILE__ << ": PHG4ParticleSvtxMap not found, reverting to true matching by momentum relations. Truth matching will be less accurate" << std::endl;
0129     }
0130   }
0131
0132   if (m_decay_descriptor.empty() && !m_decayMap->empty())
0133   {
0134     getDecayDescriptor();
0135     getNDaughters();
0136   }
0137
0138   bool reconstructableDecay = false;
0139   for (auto &iter : *m_decayMap)
0140   {
0141     Decay decay = iter.second;
0142     ++m_counter_allDecays;
0143
0144     m_all_tracks_reconstructed = findTracks(topNode, decay);
0145     if (m_all_tracks_reconstructed)
0146     {
0147       ++m_counter_acceptedDecays;
0148       reconstructableDecay = true;
0149     }
0150
0151     if (m_write_nTuple)
0152     {
0153       if (m_counter_allDecays == 1)
0154       {
0155         initializeBranches();
0156       }
0157       m_tree->Fill();
0158       resetBranches();
0159     }
0160   }
0161
0162   if (m_triggerOnDecay && !reconstructableDecay)
0163   {
0164     if (Verbosity() >= VERBOSITY_MORE)
0165     {
0166       std::cout << "No decays were reconstructable in this event, skipping" << std::endl;
0167     }
0168     return Fun4AllReturnCodes::ABORTEVENT;
0169   }
0170
0171   return Fun4AllReturnCodes::EVENT_OK;
0172 }
0173
0174 int HFTrackEfficiency::End(PHCompositeNode * /*topNode*/)
0175 {
0176   PrintEff();
0177
0178   if (m_write_nTuple && m_counter_allDecays != 0)
0179   {
0180     m_outfile->Write();
0181     m_outfile->Close();
0182     delete m_outfile;
0183   }
0184
0185   return Fun4AllReturnCodes::EVENT_OK;
0186 }
0187
0188 bool HFTrackEfficiency::findTracks(PHCompositeNode *topNode, Decay decay)
0189 {
0190   int trackableParticles[] = {11, 13, 211, 321, 2212};
0191   bool recoTrackFound = false;
0192   std::set<SvtxTrack *> selectedTracks;
0193
0194   CLHEP::HepLorentzVector motherRecoLV;
0195   CLHEP::HepLorentzVector daughterSumTrueLV;
0196   CLHEP::HepLorentzVector *motherTrueLV = new CLHEP::HepLorentzVector();
0197   CLHEP::HepLorentzVector *daughterTrueLV = new CLHEP::HepLorentzVector();
0198
0199   HepMC::GenEvent *theEvent = nullptr;
0200   if (m_geneventmap)
0201   {
0202     m_genevt = m_geneventmap->get(decay[0].first.first);
0203     assert(m_genevt);
0204
0205     theEvent = m_genevt->getEvent();
0206     HepMC::GenParticle *mother = theEvent->barcode_to_particle(decay[0].first.second);
0207     assert(mother);
0208     if (Verbosity() >= VERBOSITY_MORE)
0209     {
0210       mother->print();
0211     }
0212
0213     m_true_mother_pT = mother->momentum().perp();
0214     m_true_mother_p = std::sqrt(std::pow(mother->momentum().px(), 2) + std::pow(mother->momentum().py(), 2) + std::pow(mother->momentum().pz(), 2));  // Must have an old HepMC build, no mag function
0215     m_true_mother_eta = mother->momentum().eta();
0216
0217     HepMC::GenVertex *thisVtx = mother->production_vertex();
0218     m_primary_vtx_x = thisVtx->point3d().x();
0219     m_primary_vtx_y = thisVtx->point3d().y();
0220     m_primary_vtx_z = thisVtx->point3d().z();
0221   }
0222
0223   for (unsigned int i = 1; i < decay.size(); ++i)
0224   {
0225     m_dst_track = nullptr;
0226     int truth_ID = -1;
0227     if (std::find(std::begin(trackableParticles), std::end(trackableParticles),
0228                   std::abs(decay[i].second)) != std::end(trackableParticles))
0229     {
0230       if (theEvent && decay[i].first.second > -1)
0231       {
0232         HepMC::GenParticle *daughterHepMC = theEvent->barcode_to_particle(decay[i].first.second);
0233         if (Verbosity() >= VERBOSITY_MORE)
0234         {
0235           daughterHepMC->print();
0236         }
0237
0238         daughterTrueLV->setVectM(CLHEP::Hep3Vector(daughterHepMC->momentum().px(), daughterHepMC->momentum().py(), daughterHepMC->momentum().pz()), getParticleMass(decay[i].second));
0239         daughterSumTrueLV += *daughterTrueLV;
0240
0241         m_true_track_PID[i - 1] = daughterHepMC->pdg_id();
0242
0243         // Now get the decay vertex position
0244         HepMC::GenVertex *thisVtx = daughterHepMC->production_vertex();
0245         m_secondary_vtx_x = thisVtx->point3d().x();
0246         m_secondary_vtx_y = thisVtx->point3d().y();
0247         m_secondary_vtx_z = thisVtx->point3d().z();
0248
0249         // We need the G4 ID, not the HepMC ID to use the truth/reco map
0250         if (m_dst_truth_reco_map)
0251         {
0252           PHG4TruthInfoContainer::ConstRange range = m_truthInfo->GetParticleRange();
0253
0254           for (PHG4TruthInfoContainer::ConstIterator iter = range.first; iter != range.second; ++iter)
0255           {
0256             PHG4Particle *daughterG4 = iter->second;
0257
0258             if (std::abs(daughterG4->get_px() - daughterTrueLV->x()) <= 5e-3 &&
0259                 std::abs(daughterG4->get_py() - daughterTrueLV->y()) <= 5e-3 &&
0260                 std::abs(daughterG4->get_pz() - daughterTrueLV->z()) <= 5e-3 && daughterG4->get_pid() == decay[i].second)
0261             {
0262               truth_ID = daughterG4->get_track_id();
0263               break;
0264             }
0265           }
0266         }
0267       }
0268       else
0269       {
0270         PHG4TruthInfoContainer::ConstRange range = m_truthInfo->GetParticleRange();
0271
0272         for (PHG4TruthInfoContainer::ConstIterator iter = range.first; iter != range.second; ++iter)
0273         {
0274           PHG4Particle *daughterG4 = iter->second;
0275
0276           PHG4Particle *motherG4 = nullptr;
0277           if (daughterG4->get_parent_id() != 0)
0278           {
0279             motherG4 = m_truthInfo->GetParticle(daughterG4->get_parent_id());
0280           }
0281           else
0282           {
0283             continue;
0284           }
0285
0286           if (motherG4->get_pid() == decay[0].second && motherG4->get_barcode() == decay[0].first.second && daughterG4->get_pid() == decay[i].second && daughterG4->get_barcode() == decay[i].first.second)
0287           {
0288             if (Verbosity() >= VERBOSITY_MORE)
0289             {
0290               daughterG4->identify();
0291             }
0292
0293             CLHEP::Hep3Vector *mother3Vector = new CLHEP::Hep3Vector(motherG4->get_px(), motherG4->get_py(), motherG4->get_pz());
0294             motherTrueLV->setVectM((*mother3Vector), getParticleMass(decay[0].second));
0295             m_true_mother_pT = motherTrueLV->perp();
0296             m_true_mother_p = mother3Vector->mag();
0297             m_true_mother_eta = motherTrueLV->pseudoRapidity();
0298
0299             PHG4VtxPoint *thisVtx = m_truthInfo->GetVtx(motherG4->get_vtx_id());
0300             m_primary_vtx_x = thisVtx->get_x();
0301             m_primary_vtx_y = thisVtx->get_y();
0302             m_primary_vtx_z = thisVtx->get_z();
0303
0304             daughterTrueLV->setVectM(CLHEP::Hep3Vector(daughterG4->get_px(), daughterG4->get_py(), daughterG4->get_pz()), getParticleMass(decay[i].second));
0305             daughterSumTrueLV += *daughterTrueLV;
0306
0307             // Now get the decay vertex position
0308             thisVtx = m_truthInfo->GetVtx(daughterG4->get_vtx_id());
0309             m_secondary_vtx_x = thisVtx->get_x();
0310             m_secondary_vtx_y = thisVtx->get_y();
0311             m_secondary_vtx_z = thisVtx->get_z();
0312
0313             m_true_track_PID[i - 1] = daughterG4->get_pid();
0314             truth_ID = daughterG4->get_track_id();
0315
0316             delete mother3Vector;
0317           }
0318         }
0319       }
0320
0321       m_true_track_pT[i - 1] = (float) daughterTrueLV->perp();
0322       m_true_track_eta[i - 1] = (float) daughterTrueLV->pseudoRapidity();
0323       m_min_true_track_pT = std::min(m_true_track_pT[i - 1], m_min_true_track_pT);
0324       m_max_true_track_pT = std::max(m_true_track_pT[i - 1], m_max_true_track_pT);
0325
0326       if (m_dst_truth_reco_map && truth_ID >= 0)
0327       {
0328         std::map<float, std::set<unsigned int>> reco_set = m_dst_truth_reco_map->get(truth_ID);
0329         if (reco_set.empty())
0330         {
0331           continue;
0332         }
0333         const auto &best_weight = reco_set.rbegin();
0334         if (best_weight->second.empty())
0335         {
0336           continue;
0337         }
0338         unsigned int best_reco_id = *best_weight->second.rbegin();
0339         m_dst_track = m_input_trackMap->get(best_reco_id);
0340         if (m_dst_track)
0341         {
0342           m_used_truth_reco_map[i - 1] = true;
0343           recoTrackFound = true;
0344         }
0345       }
0346       else
0347       {
0348         for (auto &iter : *m_input_trackMap)
0349         {
0350           m_dst_track = iter.second;
0351           float delta_px = (m_dst_track->get_px() - daughterTrueLV->px()) / daughterTrueLV->px();
0352           float delta_py = (m_dst_track->get_py() - daughterTrueLV->py()) / daughterTrueLV->py();
0353           float delta_pz = (m_dst_track->get_pz() - daughterTrueLV->pz()) / daughterTrueLV->pz();
0354
0355           if (std::abs(delta_px) <= m_truthRecoMatchPercent && std::abs(delta_py) <= m_truthRecoMatchPercent && std::abs(delta_pz) <= m_truthRecoMatchPercent)
0356           {
0357             recoTrackFound = true;
0358             break;
0359           }
0360         }
0361       }
0362
0363       if (recoTrackFound)
0364       {
0365         selectedTracks.insert(m_dst_track);
0366         if (Verbosity() >= VERBOSITY_MORE)
0367         {
0368           m_dst_track->identify();
0369         }
0370         m_reco_track_exists[i - 1] = true;
0371         m_reco_track_pT[i - 1] = m_dst_track->get_pt();
0372         m_reco_track_eta[i - 1] = m_dst_track->get_eta();
0373         m_reco_track_chi2nDoF[i - 1] = m_dst_track->get_chisq() / m_dst_track->get_ndf();
0374         if (m_dst_track->get_silicon_seed())
0375         {
0376           m_reco_track_silicon_seeds[i - 1] = static_cast<int>(m_dst_track->get_silicon_seed()->size_cluster_keys());
0377         }
0378         else
0379         {
0380           m_reco_track_silicon_seeds[i - 1] = 0;
0381         }
0382         m_reco_track_tpc_seeds[i - 1] = static_cast<int>(m_dst_track->get_tpc_seed()->size_cluster_keys());
0383         m_min_reco_track_pT = std::min(m_reco_track_pT[i - 1], m_min_reco_track_pT);
0384         m_max_reco_track_pT = std::max(m_reco_track_pT[i - 1], m_max_reco_track_pT);
0385
0386         CLHEP::HepLorentzVector *daughterRecoLV = new CLHEP::HepLorentzVector();
0387         daughterRecoLV->setVectM(CLHEP::Hep3Vector(m_dst_track->get_px(), m_dst_track->get_py(), m_dst_track->get_pz()), getParticleMass(m_true_track_PID[i - 1]));
0388
0389         motherRecoLV += *daughterRecoLV;
0390         delete daughterRecoLV;
0391       }
0392
0393       recoTrackFound = false;
0394     }
0395   }
0396
0397   m_true_mother_mass = daughterSumTrueLV.m();
0398   bool foundDecay = true;
0399   if (selectedTracks.size() == m_nDaughters)
0400   {
0401     m_reco_mother_mass = motherRecoLV.m();
0402     if (m_write_track_map)
0403     {
0404       m_output_trackMap = findNode::getClass<SvtxTrackMap>(topNode, outputNodeName);
0405       for (const auto &track : selectedTracks)
0406       {
0407         m_output_trackMap->insertWithKey(track, track->get_id());
0408       }
0409     }
0410   }
0411   else
0412   {
0413     foundDecay = false;
0414   }
0415
0416   selectedTracks.clear();
0417
0418   delete motherTrueLV;
0419   delete daughterTrueLV;
0420
0421   return foundDecay;
0422 }
0423
0424 void HFTrackEfficiency::initializeBranches()
0425 {
0426   m_outfile = new TFile(m_outfile_name.c_str(), "RECREATE");
0427   m_tree = new TTree("HFTrackEfficiency", "HFTrackEfficiency");
0428   m_tree->OptimizeBaskets();
0429   m_tree->SetAutoSave(-5e6);  // Save the output file every 5MB
0430
0431   m_tree->Branch("all_tracks_reconstructed", &m_all_tracks_reconstructed, "all_tracks_reconstructed/O");
0432   m_tree->Branch("true_mother_mass", &m_true_mother_mass, "true_mother_mass/F");
0433   m_tree->Branch("reco_mother_mass", &m_reco_mother_mass, "reco_mother_mass/F");
0434   m_tree->Branch("true_mother_pT", &m_true_mother_pT, "true_mother_pT/F");
0435   m_tree->Branch("true_mother_p", &m_true_mother_p, "true_mother_p/F");
0436   m_tree->Branch("true_mother_eta", &m_true_mother_eta, "true_mother_eta/F");
0437   m_tree->Branch("min_true_track_pT", &m_min_true_track_pT, "min_true_track_pT/F");
0438   m_tree->Branch("min_reco_track_pT", &m_min_reco_track_pT, "min_reco_track_pT/F");
0439   m_tree->Branch("max_true_track_pT", &m_max_true_track_pT, "max_true_track_pT/F");
0440   m_tree->Branch("max_reco_track_pT", &m_max_reco_track_pT, "max_reco_track_pT/F");
0441
0442   for (unsigned int iTrack = 0; iTrack < m_nDaughters; ++iTrack)
0443   {
0444     std::string daughter_number = "track_" + std::to_string(iTrack + 1);
0445     m_tree->Branch("reco_" + TString(daughter_number) + "_exists", &m_reco_track_exists[iTrack], "reco_" + TString(daughter_number) + "_exists/O");
0446     m_tree->Branch("reco_" + TString(daughter_number) + "_used_truth_reco_map", &m_used_truth_reco_map[iTrack], "reco_" + TString(daughter_number) + "_used_truth_reco_map/O");
0447     m_tree->Branch("true_" + TString(daughter_number) + "_pT", &m_true_track_pT[iTrack], "true_" + TString(daughter_number) + "_pT/F");
0448     m_tree->Branch("reco_" + TString(daughter_number) + "_pT", &m_reco_track_pT[iTrack], "reco_" + TString(daughter_number) + "_pT/F");
0449     m_tree->Branch("true_" + TString(daughter_number) + "_eta", &m_true_track_eta[iTrack], "true_" + TString(daughter_number) + "_eta/F");
0450     m_tree->Branch("reco_" + TString(daughter_number) + "_eta", &m_reco_track_eta[iTrack], "reco_" + TString(daughter_number) + "_eta/F");
0451     m_tree->Branch("true_" + TString(daughter_number) + "_PID", &m_true_track_PID[iTrack], "true_" + TString(daughter_number) + "_PID/F");
0452     m_tree->Branch("reco_" + TString(daughter_number) + "_chi2nDoF", &m_reco_track_chi2nDoF[iTrack], "reco_" + TString(daughter_number) + "_chi2nDoF/F");
0453     m_tree->Branch("reco_" + TString(daughter_number) + "_silicon_seeds", &m_reco_track_silicon_seeds[iTrack], "reco_" + TString(daughter_number) + "_silicon_seeds/I");
0454     m_tree->Branch("reco_" + TString(daughter_number) + "_tpc_seeds", &m_reco_track_tpc_seeds[iTrack], "reco_" + TString(daughter_number) + "_tpc_seeds/I");
0455   }
0456
0457   m_tree->Branch("true_primary_vertex_x", &m_primary_vtx_x, "true_primary_vertex_x/F");
0458   m_tree->Branch("true_primary_vertex_y", &m_primary_vtx_y, "true_primary_vertex_y/F");
0459   m_tree->Branch("true_primary_vertex_z", &m_primary_vtx_z, "true_primary_vertex_z/F");
0460   m_tree->Branch("true_secondary_vertex_x", &m_secondary_vtx_x, "true_secondary_vertex_x/F");
0461   m_tree->Branch("true_secondary_vertex_y", &m_secondary_vtx_y, "true_secondary_vertex_y/F");
0462   m_tree->Branch("true_secondary_vertex_z", &m_secondary_vtx_z, "true_secondary_vertex_z/F");
0463 }
0464
0465 void HFTrackEfficiency::resetBranches()
0466 {
0467   m_all_tracks_reconstructed = false;
0468   m_true_mother_mass = std::numeric_limits<float>::quiet_NaN();
0469   m_reco_mother_mass = std::numeric_limits<float>::quiet_NaN();
0470   m_true_mother_pT = std::numeric_limits<float>::quiet_NaN();
0471   m_true_mother_p = std::numeric_limits<float>::quiet_NaN();
0472   m_true_mother_eta = std::numeric_limits<float>::quiet_NaN();
0473   m_min_true_track_pT = std::numeric_limits<float>::max();
0474   m_min_reco_track_pT = std::numeric_limits<float>::max();
0475   m_max_true_track_pT = -1 * std::numeric_limits<float>::max();
0476   m_max_reco_track_pT = -1 * std::numeric_limits<float>::max();
0477   for (unsigned int iTrack = 0; iTrack < m_nDaughters; ++iTrack)
0478   {
0479     m_reco_track_exists[iTrack] = false;
0480     m_used_truth_reco_map[iTrack] = false;
0481     m_true_track_pT[iTrack] = std::numeric_limits<float>::quiet_NaN();
0482     m_reco_track_pT[iTrack] = std::numeric_limits<float>::quiet_NaN();
0483     m_true_track_eta[iTrack] = std::numeric_limits<float>::quiet_NaN();
0484     m_reco_track_eta[iTrack] = std::numeric_limits<float>::quiet_NaN();
0485     m_true_track_PID[iTrack] = std::numeric_limits<float>::quiet_NaN();
0486     m_reco_track_chi2nDoF[iTrack] = std::numeric_limits<float>::quiet_NaN();
0487     m_reco_track_silicon_seeds[iTrack] = 0;
0488     m_reco_track_tpc_seeds[iTrack] = 0;
0489   }
0490
0491   m_primary_vtx_x = std::numeric_limits<float>::quiet_NaN();
0492   m_primary_vtx_y = std::numeric_limits<float>::quiet_NaN();
0493   m_primary_vtx_z = std::numeric_limits<float>::quiet_NaN();
0494   m_secondary_vtx_x = std::numeric_limits<float>::quiet_NaN();
0495   m_secondary_vtx_y = std::numeric_limits<float>::quiet_NaN();
0496   m_secondary_vtx_z = std::numeric_limits<float>::quiet_NaN();
0497 }
0498
0499 void HFTrackEfficiency::getDecayDescriptor()
0500 {
0501   Decay decay = m_decayMap->begin()->second;
0502   m_decay_descriptor = getParticleName(decay[0].second);
0503   m_decay_descriptor += " ->";
0504   for (unsigned int i = 1; i < decay.size(); ++i)
0505   {
0506     m_decay_descriptor += " ";
0507     m_decay_descriptor += getParticleName(decay[i].second);
0508   }
0509 }
0510
0511 void HFTrackEfficiency::getNDaughters()
0512 {
0513   int trackableParticles[] = {11, 13, 211, 321, 2212};
0514   m_nDaughters = 0;
0515   Decay decay = m_decayMap->begin()->second;
0516   for (unsigned int i = 1; i < decay.size(); ++i)
0517   {
0518     if (std::find(std::begin(trackableParticles), std::end(trackableParticles),
0519                   std::abs(decay[i].second)) != std::end(trackableParticles))
0520     {
0521       ++m_nDaughters;
0522     }
0523   }
0524 }
0525
0526 std::string HFTrackEfficiency::getParticleName(const int PDGID)
0527 {
0528   return TDatabasePDG::Instance()->GetParticle(PDGID)->GetName();
0529 }
0530
0531 float HFTrackEfficiency::getParticleMass(const int PDGID)
0532 {
0533   return TDatabasePDG::Instance()->GetParticle(PDGID)->Mass();
0534 }
0535
0536 //____________________________________________________________________________..
0537 void HFTrackEfficiency::PrintEff()
0538 {
0539   std::cout << "\n--------------- Heavy Flavor Tracking Efficiency ---------------" << std::endl;
0540   std::cout << "Tracking Efficiency for " << m_decay_descriptor << std::endl;
0541   std::cout << "Number of decays fully in acceptance: " << m_counter_allDecays << std::endl;
0542   std::cout << "Number of decays with all tracks reconstructed: " << m_counter_acceptedDecays << std::endl;
0543   std::cout << "Efficiency: " << std::setw(4) << (float) m_counter_acceptedDecays * 100 / float(m_counter_allDecays) << "%" << std::endl;
0544   std::cout << "----------------------------------------------------------------\n"
0545             << std::endl;
0546 }