g4simulation/g4eval/JetTruthEval.cc

0001 #include "JetTruthEval.h"
0002
0003 #include "CaloTruthEval.h"
0004 #include "SvtxTruthEval.h"
0005
0006 #include <jetbase/Jet.h>
0007 #include <jetbase/JetContainer.h>
0008
0009 #include <g4main/PHG4TruthInfoContainer.h>
0010
0011 #include <phool/getClass.h>
0012 #include <phool/phool.h>
0013
0014 #include <cassert>
0015 #include <cstdlib>
0016 #include <iostream>
0017 #include <map>
0018 #include <utility>
0019
0020 JetTruthEval::JetTruthEval(PHCompositeNode* topNode,
0021                            const std::string& truthjetname)
0022   : _truthjetname(truthjetname)
0023   , _svtxevalstack(topNode)
0024   , _cemcevalstack(topNode, "CEMC")
0025   , _hcalinevalstack(topNode, "HCALIN")
0026   , _hcaloutevalstack(topNode, "HCALOUT")
0027   , _femcevalstack(topNode, "FEMC")
0028   , _fhcalevalstack(topNode, "FHCAL")
0029   , _eemcevalstack(topNode, "EEMC")
0030 {
0031   get_node_pointers(topNode);
0032 }
0033
0034 JetTruthEval::~JetTruthEval()
0035 {
0036   if (_verbosity > 0)
0037   {
0038     if ((_errors > 0) || (_verbosity > 1))
0039     {
0040       std::cout << "JetTruthEval::~JetTruthEval() - Error Count: " << _errors << std::endl;
0041     }
0042   }
0043 }
0044
0045 void JetTruthEval::next_event(PHCompositeNode* topNode)
0046 {
0047   _svtxevalstack.next_event(topNode);
0048   _cemcevalstack.next_event(topNode);
0049   _hcalinevalstack.next_event(topNode);
0050   _hcaloutevalstack.next_event(topNode);
0051   _femcevalstack.next_event(topNode);
0052   _fhcalevalstack.next_event(topNode);
0053
0054   _cache_all_truth_particles.clear();
0055   _cache_all_truth_showers.clear();
0056   _cache_all_truth_hits.clear();
0057   _cache_get_truth_jet.clear();
0058
0059   get_node_pointers(topNode);
0060 }
0061
0062 std::set<PHG4Particle*> JetTruthEval::all_truth_particles(Jet* truthjet)
0063 {
0064   if (_strict)
0065   {
0066     assert(truthjet);
0067   }
0068   else if (!truthjet)
0069   {
0070     ++_errors;
0071     return std::set<PHG4Particle*>();
0072   }
0073
0074   if (_do_cache)
0075   {
0076     std::map<Jet*, std::set<PHG4Particle*> >::iterator iter =
0077         _cache_all_truth_particles.find(truthjet);
0078     if (iter != _cache_all_truth_particles.end())
0079     {
0080       return iter->second;
0081     }
0082   }
0083
0084   std::set<PHG4Particle*> truth_particles;
0085
0086   // loop over all the entries in the truthjet
0087   /* for (Jet::ConstIter iter = truthjet->comp_begin(); */
0088   /*      iter != truthjet->comp_end(); */
0089   /*      ++iter) */
0090   for (const auto& iter : truthjet->get_comp_vec())  // vector of pair<Jet::SRC, unsigned int>
0091   {
0092     Jet::SRC source = iter.first;
0093     unsigned int index = iter.second;
0094     if (source != Jet::PARTICLE)
0095     {
0096       std::cout << PHWHERE << " truth jet contains something other than particles!" << std::endl;
0097       exit(-1);
0098     }
0099
0100     PHG4Particle* truth_particle = _truthinfo->GetParticle(index);
0101
0102     if (_strict)
0103     {
0104       assert(truth_particle);
0105     }
0106     else if (!truth_particle)
0107     {
0108       ++_errors;
0109       continue;
0110     }
0111
0112     truth_particles.insert(truth_particle);
0113   }
0114
0115   if (_do_cache)
0116   {
0117     _cache_all_truth_particles.insert(std::make_pair(truthjet, truth_particles));
0118   }
0119
0120   return truth_particles;
0121 }
0122
0123 std::set<PHG4Shower*> JetTruthEval::all_truth_showers(Jet* truthjet)
0124 {
0125   if (_strict)
0126   {
0127     assert(truthjet);
0128   }
0129   else if (!truthjet)
0130   {
0131     ++_errors;
0132     return std::set<PHG4Shower*>();
0133   }
0134
0135   if (_do_cache)
0136   {
0137     std::map<Jet*, std::set<PHG4Shower*> >::iterator iter =
0138         _cache_all_truth_showers.find(truthjet);
0139     if (iter != _cache_all_truth_showers.end())
0140     {
0141       return iter->second;
0142     }
0143   }
0144
0145   std::set<PHG4Shower*> truth_showers;
0146
0147   std::set<PHG4Particle*> truth_particles = all_truth_particles(truthjet);
0148
0149   // loop over all the entries in the truthjet
0150   for (auto particle : truth_particles)
0151   {
0152     if (_strict)
0153     {
0154       assert(particle);
0155     }
0156     else if (!particle)
0157     {
0158       ++_errors;
0159       continue;
0160     }
0161
0162     // any calo truth eval module would work here...
0163     CaloTruthEval* cemc_truth_eval = _cemcevalstack.get_truth_eval();
0164     PHG4Shower* shower = cemc_truth_eval->get_primary_shower(particle);
0165     if (shower)
0166     {
0167       truth_showers.insert(shower);
0168     }
0169   }
0170
0171   if (_do_cache)
0172   {
0173     _cache_all_truth_showers.insert(std::make_pair(truthjet, truth_showers));
0174   }
0175
0176   return truth_showers;
0177 }
0178
0179 std::set<PHG4Hit*> JetTruthEval::all_truth_hits(Jet* truthjet)
0180 {
0181   if (_strict)
0182   {
0183     assert(truthjet);
0184   }
0185   else if (!truthjet)
0186   {
0187     ++_errors;
0188     return std::set<PHG4Hit*>();
0189   }
0190
0191   if (_do_cache)
0192   {
0193     std::map<Jet*, std::set<PHG4Hit*> >::iterator iter =
0194         _cache_all_truth_hits.find(truthjet);
0195     if (iter != _cache_all_truth_hits.end())
0196     {
0197       return iter->second;
0198     }
0199   }
0200
0201   std::set<PHG4Hit*> truth_hits;
0202
0203   std::set<PHG4Particle*> truth_particles = all_truth_particles(truthjet);
0204
0205   // loop over all the entries in the truthjet
0206   for (auto particle : truth_particles)
0207   {
0208     if (_strict)
0209     {
0210       assert(particle);
0211     }
0212     else if (!particle)
0213     {
0214       ++_errors;
0215       continue;
0216     }
0217
0218     // ask the svtx truth eval to backtrack the particles to g4hits
0219     SvtxTruthEval* svtx_truth_eval = _svtxevalstack.get_truth_eval();
0220     std::set<PHG4Hit*> svtx_g4hits = svtx_truth_eval->all_truth_hits(particle);
0221
0222     for (auto g4hit : svtx_g4hits)
0223     {
0224       if (_strict)
0225       {
0226         assert(g4hit);
0227       }
0228       else if (!g4hit)
0229       {
0230         ++_errors;
0231         continue;
0232       }
0233
0234       truth_hits.insert(g4hit);
0235     }
0236   }
0237
0238   if (_do_cache)
0239   {
0240     _cache_all_truth_hits.insert(std::make_pair(truthjet, truth_hits));
0241   }
0242
0243   return truth_hits;
0244 }
0245
0246 Jet* JetTruthEval::get_truth_jet(PHG4Particle* particle)
0247 {
0248   if (_strict)
0249   {
0250     assert(particle);
0251   }
0252   else if (!particle)
0253   {
0254     ++_errors;
0255     return nullptr;
0256   }
0257
0258   if (_do_cache)
0259   {
0260     std::map<PHG4Particle*, Jet*>::iterator iter =
0261         _cache_get_truth_jet.find(particle);
0262     if (iter != _cache_get_truth_jet.end())
0263     {
0264       return iter->second;
0265     }
0266   }
0267
0268   Jet* truth_jet = nullptr;
0269
0270   // loop over all jets and look for this particle...
0271   for (auto candidate : *_truthjets)
0272   {
0273     /* Jet* candidate = _truthjet.second; */
0274
0275     // loop over all consituents and look for this particle
0276     for (const std::pair<Jet::SRC, unsigned int>& jter
0277         : candidate->get_comp_vec())
0278     {
0279       unsigned int index = jter.second;
0280
0281       PHG4Particle* constituent = _truthinfo->GetParticle(index);
0282       if (_strict)
0283       {
0284         assert(constituent);
0285       }
0286       else if (!constituent)
0287       {
0288         ++_errors;
0289         continue;
0290       }
0291
0292       if (get_svtx_eval_stack()->get_truth_eval()->are_same_particle(constituent, particle))
0293       {
0294         truth_jet = candidate;
0295         break;
0296       }
0297     }
0298
0299     if (truth_jet)
0300     {
0301       break;
0302     }
0303   }
0304
0305   if (_do_cache)
0306   {
0307     _cache_get_truth_jet.insert(std::make_pair(particle, truth_jet));
0308   }
0309
0310   return truth_jet;
0311 }
0312
0313 void JetTruthEval::get_node_pointers(PHCompositeNode* topNode)
0314 {
0315   _truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0316   if (!_truthinfo)
0317   {
0318     std::cout << PHWHERE << " ERROR: Can't find G4TruthInfo" << std::endl;
0319     exit(-1);
0320   }
0321
0322   _truthjets = findNode::getClass<JetContainer>(topNode, _truthjetname);
0323   if (!_truthjets)
0324   {
0325     std::cout << PHWHERE << " ERROR: Can't find " << _truthjetname << std::endl;
0326     exit(-1);
0327   }
0328
0329   return;
0330 }