g4simulation/g4eval/CaloTruthEval.cc

0001 #include "CaloTruthEval.h"
0002
0003 #include <g4main/PHG4HitContainer.h>
0004 #include <g4main/PHG4HitDefs.h>
0005 #include <g4main/PHG4Shower.h>
0006 #include <g4main/PHG4TruthInfoContainer.h>
0007
0008 #include <phool/getClass.h>
0009
0010 #include <cassert>
0011 #include <cmath>
0012 #include <iostream>
0013 #include <map>
0014 #include <set>
0015 #include <utility>
0016
0017 CaloTruthEval::CaloTruthEval(PHCompositeNode* topNode, const std::string& caloname)
0018   : _basetrutheval(topNode)
0019   , _caloname(caloname)
0020   , _caloid(PHG4HitDefs::get_volume_id(caloname))
0021 {
0022   get_node_pointers(topNode);
0023 }
0024
0025 CaloTruthEval::~CaloTruthEval()
0026 {
0027   if (_verbosity > 0)
0028   {
0029     if ((_errors > 0) || (_verbosity > 1))
0030     {
0031       std::cout << "CaloTruthEval::~CaloTruthEval() - Error Count: " << _errors << std::endl;
0032     }
0033   }
0034 }
0035
0036 void CaloTruthEval::next_event(PHCompositeNode* topNode)
0037 {
0038   _cache_get_shower_energy_deposit.clear();
0039   _cache_all_truth_hits_g4shower.clear();
0040   _cache_all_truth_hits_g4particle.clear();
0041   _cache_get_primary_particle_g4hit.clear();
0042   _cache_get_shower_hits_from_primary.clear();
0043
0044   _basetrutheval.next_event(topNode);
0045
0046   get_node_pointers(topNode);
0047 }
0048
0049 bool CaloTruthEval::has_reduced_node_pointers()
0050 {
0051   if (!_basetrutheval.has_reduced_node_pointers())
0052   {
0053     return false;
0054   }
0055
0056   if (_strict)
0057   {
0058     assert(_truthinfo);
0059   }
0060   else if (!_truthinfo)
0061   {
0062     return false;
0063   }
0064
0065   return true;
0066 }
0067
0068 PHG4Shower* CaloTruthEval::get_primary_shower(PHG4Shower* shower)
0069 {
0070   return _basetrutheval.get_primary_shower(shower);
0071 }
0072
0073 PHG4Shower* CaloTruthEval::get_primary_shower(PHG4Particle* particle)
0074 {
0075   return _basetrutheval.get_primary_shower(particle);
0076 }
0077
0078 PHG4Particle* CaloTruthEval::get_primary_particle(PHG4Shower* shower)
0079 {
0080   return _basetrutheval.get_primary_particle(shower);
0081 }
0082
0083 PHG4Particle* CaloTruthEval::get_primary_particle(PHG4Particle* particle)
0084 {
0085   return _basetrutheval.get_primary_particle(particle);
0086 }
0087
0088 int CaloTruthEval::get_embed(PHG4Particle* particle)
0089 {
0090   return _basetrutheval.get_embed(particle);
0091 }
0092
0093 PHG4VtxPoint* CaloTruthEval::get_vertex(PHG4Particle* particle)
0094 {
0095   return _basetrutheval.get_vertex(particle);
0096 }
0097
0098 bool CaloTruthEval::are_same_shower(PHG4Shower* s1, PHG4Shower* s2)
0099 {
0100   return _basetrutheval.are_same_shower(s1, s2);
0101 }
0102
0103 bool CaloTruthEval::are_same_particle(PHG4Particle* p1, PHG4Particle* p2)
0104 {
0105   return _basetrutheval.are_same_particle(p1, p2);
0106 }
0107
0108 bool CaloTruthEval::are_same_vertex(PHG4VtxPoint* vtx1, PHG4VtxPoint* vtx2)
0109 {
0110   return _basetrutheval.are_same_vertex(vtx1, vtx2);
0111 }
0112
0113 bool CaloTruthEval::is_primary(PHG4Shower* shower)
0114 {
0115   return _basetrutheval.is_primary(shower);
0116 }
0117
0118 bool CaloTruthEval::is_primary(PHG4Particle* particle)
0119 {
0120   return _basetrutheval.is_primary(particle);
0121 }
0122
0123 float CaloTruthEval::get_shower_energy_deposit(PHG4Particle* primary)
0124 {
0125   if (!has_reduced_node_pointers())
0126   {
0127     ++_errors;
0128     return NAN;
0129   }
0130
0131   if (_strict)
0132   {
0133     assert(primary);
0134   }
0135   else if (!primary)
0136   {
0137     ++_errors;
0138     return NAN;
0139   }
0140
0141   if (!is_primary(primary))
0142   {
0143     return NAN;
0144   }
0145
0146   primary = get_primary_particle(primary);
0147
0148   if (_strict)
0149   {
0150     assert(primary);
0151   }
0152   else if (!primary)
0153   {
0154     ++_errors;
0155     return NAN;
0156   }
0157
0158   if (_do_cache)
0159   {
0160     std::map<PHG4Particle*, float>::iterator iter =
0161         _cache_get_shower_energy_deposit.find(primary);
0162     if (iter != _cache_get_shower_energy_deposit.end())
0163     {
0164       return iter->second;
0165     }
0166   }
0167
0168   float shower_e = 0.0;
0169   PHG4Shower* shower = get_primary_shower(primary);
0170   if (shower)
0171   {
0172     shower_e = shower->get_edep(_g4hit_container_id);
0173   }
0174
0175   if (_do_cache)
0176   {
0177     _cache_get_shower_energy_deposit.insert(std::make_pair(primary, shower_e));
0178   }
0179
0180   return shower_e;
0181 }
0182
0183 float CaloTruthEval::get_shower_eh_ratio(PHG4Particle* primary)
0184 {
0185   if (!has_full_node_pointers())
0186   {
0187     ++_errors;
0188     return NAN;
0189   }
0190
0191   if (_strict)
0192   {
0193     assert(primary);
0194   }
0195   else if (!primary)
0196   {
0197     ++_errors;
0198     return NAN;
0199   }
0200
0201   if (!is_primary(primary))
0202   {
0203     return NAN;
0204   }
0205
0206   PHG4Shower* shower = get_primary_shower(primary);
0207   if (!shower)
0208   {
0209     return 0.0;
0210   }
0211
0212   float ratio = shower->get_eh_ratio(get_caloid());
0213
0214   return ratio;
0215 }
0216
0217 bool CaloTruthEval::has_full_node_pointers()
0218 {
0219   if (!_basetrutheval.has_full_node_pointers())
0220   {
0221     return false;
0222   }
0223
0224   if (_strict)
0225   {
0226     assert(_truthinfo);
0227   }
0228   else if (!_truthinfo)
0229   {
0230     return false;
0231   }
0232
0233   if (_strict)
0234   {
0235     assert(_g4hits);
0236   }
0237   if (!_g4hits)
0238   {
0239     return false;
0240   }
0241
0242   return true;
0243 }
0244
0245 std::set<PHG4Hit*> CaloTruthEval::all_truth_hits(PHG4Shower* shower)
0246 {
0247   if (!has_full_node_pointers())
0248   {
0249     ++_errors;
0250     return std::set<PHG4Hit*>();
0251   }
0252
0253   if (_strict)
0254   {
0255     assert(shower);
0256   }
0257   else if (!shower || !_g4hits)
0258   {
0259     ++_errors;
0260     return std::set<PHG4Hit*>();
0261   }
0262
0263   if (_do_cache)
0264   {
0265     std::map<PHG4Shower*, std::set<PHG4Hit*> >::iterator iter =
0266         _cache_all_truth_hits_g4shower.find(shower);
0267     if (iter != _cache_all_truth_hits_g4shower.end())
0268     {
0269       return iter->second;
0270     }
0271   }
0272
0273   std::set<PHG4Hit*> truth_hits;
0274
0275   // loop over all g4hits on the shower
0276   PHG4Shower::HitIdIter iter = shower->find_g4hit_id(_g4hit_container_id);
0277   if (iter != shower->end_g4hit_id())
0278   {
0279     return truth_hits;
0280   }
0281
0282   for (unsigned long long jter : iter->second)
0283   {
0284     PHG4Hit* g4hit = _g4hits->findHit(jter);
0285
0286     if (_strict)
0287     {
0288       assert(g4hit);
0289     }
0290     else if (!g4hit)
0291     {
0292       ++_errors;
0293     }
0294
0295     if (g4hit)
0296     {
0297       truth_hits.insert(g4hit);
0298     }
0299   }
0300
0301   if (_do_cache)
0302   {
0303     _cache_all_truth_hits_g4shower.insert(std::make_pair(shower, truth_hits));
0304   }
0305
0306   return truth_hits;
0307 }
0308
0309 std::set<PHG4Hit*> CaloTruthEval::all_truth_hits(PHG4Particle* particle)
0310 {
0311   if (!has_full_node_pointers())
0312   {
0313     ++_errors;
0314     return std::set<PHG4Hit*>();
0315   }
0316
0317   if (_strict)
0318   {
0319     assert(particle);
0320   }
0321   else if (!particle)
0322   {
0323     ++_errors;
0324     return std::set<PHG4Hit*>();
0325   }
0326   if (!_g4hits)
0327   {
0328     ++_errors;
0329     return std::set<PHG4Hit*>();
0330   }
0331
0332   if (_do_cache)
0333   {
0334     std::map<PHG4Particle*, std::set<PHG4Hit*> >::iterator iter =
0335         _cache_all_truth_hits_g4particle.find(particle);
0336     if (iter != _cache_all_truth_hits_g4particle.end())
0337     {
0338       return iter->second;
0339     }
0340   }
0341
0342   std::set<PHG4Hit*> truth_hits;
0343
0344   // loop over all the g4hits
0345   for (PHG4HitContainer::ConstIterator g4iter = _g4hits->getHits().first;
0346        g4iter != _g4hits->getHits().second;
0347        ++g4iter)
0348   {
0349     PHG4Hit* g4hit = g4iter->second;
0350     if (is_g4hit_from_particle(g4hit, particle))
0351     {
0352       continue;
0353     }
0354     truth_hits.insert(g4hit);
0355   }
0356
0357   if (_do_cache)
0358   {
0359     _cache_all_truth_hits_g4particle.insert(std::make_pair(particle, truth_hits));
0360   }
0361
0362   return truth_hits;
0363 }
0364
0365 PHG4Particle* CaloTruthEval::get_parent_particle(PHG4Hit* g4hit)
0366 {
0367   return _basetrutheval.get_particle(g4hit);
0368 }
0369
0370 PHG4Particle* CaloTruthEval::get_primary_particle(PHG4Hit* g4hit)
0371 {
0372   if (!has_full_node_pointers())
0373   {
0374     ++_errors;
0375     return nullptr;
0376   }
0377
0378   if (_strict)
0379   {
0380     assert(g4hit);
0381   }
0382   else if (!g4hit)
0383   {
0384     ++_errors;
0385     return nullptr;
0386   }
0387
0388   if (_do_cache)
0389   {
0390     std::map<PHG4Hit*, PHG4Particle*>::iterator iter =
0391         _cache_get_primary_particle_g4hit.find(g4hit);
0392     if (iter != _cache_get_primary_particle_g4hit.end())
0393     {
0394       return iter->second;
0395     }
0396   }
0397
0398   PHG4Particle* primary = _basetrutheval.get_primary_particle(g4hit);
0399
0400   if (_do_cache)
0401   {
0402     _cache_get_primary_particle_g4hit.insert(std::make_pair(g4hit, primary));
0403   }
0404
0405   if (_strict)
0406   {
0407     assert(primary);
0408   }
0409   else if (!primary)
0410   {
0411     ++_errors;
0412   }
0413
0414   return primary;
0415 }
0416
0417 bool CaloTruthEval::is_g4hit_from_particle(PHG4Hit* g4hit, PHG4Particle* particle)
0418 {
0419   return _basetrutheval.is_g4hit_from_particle(g4hit, particle);
0420 }
0421
0422 std::set<PHG4Hit*> CaloTruthEval::get_shower_hits_from_primary(PHG4Particle* primary)
0423 {
0424   if (!has_full_node_pointers())
0425   {
0426     ++_errors;
0427     return std::set<PHG4Hit*>();
0428   }
0429
0430   if (_strict)
0431   {
0432     assert(primary);
0433   }
0434   else if (!primary)
0435   {
0436     ++_errors;
0437     return std::set<PHG4Hit*>();
0438   }
0439
0440   if (!is_primary(primary))
0441   {
0442     return std::set<PHG4Hit*>();
0443   }
0444
0445   primary = get_primary_particle(primary);
0446
0447   if (_strict)
0448   {
0449     assert(primary);
0450   }
0451   else if (!primary)
0452   {
0453     ++_errors;
0454     return std::set<PHG4Hit*>();
0455   }
0456
0457   if (_do_cache)
0458   {
0459     std::map<PHG4Particle*, std::set<PHG4Hit*> >::iterator iter =
0460         _cache_get_shower_hits_from_primary.find(primary);
0461     if (iter != _cache_get_shower_hits_from_primary.end())
0462     {
0463       return iter->second;
0464     }
0465   }
0466
0467   std::set<PHG4Hit*> truth_hits;
0468
0469   PHG4Shower* shower = get_primary_shower(primary);
0470   if (shower)
0471   {
0472     truth_hits = all_truth_hits(shower);
0473   }
0474
0475   if (_do_cache)
0476   {
0477     _cache_get_shower_hits_from_primary.insert(std::make_pair(primary, truth_hits));
0478   }
0479
0480   return truth_hits;
0481 }
0482
0483 void CaloTruthEval::get_node_pointers(PHCompositeNode* topNode)
0484 {
0485   // need things off of the DST...
0486   _truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0487
0488   std::string name = "G4HIT_" + _caloname;
0489   _g4hits = findNode::getClass<PHG4HitContainer>(topNode, name.c_str());
0490   _g4hit_container_id = PHG4HitDefs::get_volume_id(name);
0491
0492   return;
0493 }