Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-05 08:17:57 

0001 #include "JetRecoEval.h"
0002 
0003 #include "CaloEvalStack.h"
0004 #include "CaloRawClusterEval.h"
0005 #include "CaloRawTowerEval.h"  // for CaloRawTowerEval
0006 #include "JetTruthEval.h"
0007 #include "SvtxEvalStack.h"
0008 #include "SvtxTrackEval.h"
0009 
0010 #include <calobase/RawCluster.h>
0011 #include <calobase/RawClusterContainer.h>
0012 #include <calobase/RawTower.h>
0013 #include <calobase/RawTowerContainer.h>
0014 
0015 #include <jetbase/Jet.h>
0016 #include <jetbase/JetContainer.h>
0017 
0018 #include <g4main/PHG4Particle.h>
0019 
0020 #include <trackbase_historic/SvtxTrack.h>
0021 #include <trackbase_historic/SvtxTrackMap.h>
0022 
0023 #include <phool/getClass.h>
0024 #include <phool/phool.h>
0025 
0026 #include <cassert>
0027 #include <cfloat>
0028 #include <cmath>
0029 #include <cstdlib>
0030 #include <iostream>
0031 #include <map>
0032 #include <set>
0033 #include <string>
0034 
0035 JetRecoEval::JetRecoEval(PHCompositeNode* topNode,
0036                          const std::string& recojetname,
0037                          const std::string& truthjetname)
0038   : _jettrutheval(topNode, truthjetname)
0039   , _recojetname(recojetname)
0040   , _truthjetname(truthjetname)
0041 {
0042   get_node_pointers(topNode);
0043 }
0044 
0045 JetRecoEval::~JetRecoEval()
0046 {
0047   if (_verbosity > 0)
0048   {
0049     if ((_errors > 0) || (_verbosity > 1))
0050     {
0051       std::cout << "JetRecoEval::~JetRecoEval() - Error Count: " << _errors << std::endl;
0052     }
0053   }
0054 }
0055 
0056 void JetRecoEval::next_event(PHCompositeNode* topNode)
0057 {
0058   _cache_all_truth_showers.clear();
0059   _cache_all_truth_particles.clear();
0060   _cache_all_truth_jets.clear();
0061   _cache_max_truth_jet_by_energy.clear();
0062   _cache_all_jets_from.clear();
0063   _cache_best_jet_from.clear();
0064   _cache_get_energy_contribution.clear();
0065   _cache_get_energy_contribution_src.clear();
0066   _cache_all_truth_hits.clear();
0067 
0068   _jettrutheval.next_event(topNode);
0069 
0070   get_node_pointers(topNode);
0071 }
0072 
0073 void JetRecoEval::set_track_nodename(const std::string& name)
0074 {
0075   m_TrackNodeName = name;
0076   _jettrutheval.set_track_nodename(name);
0077 }
0078 
0079 std::set<PHG4Shower*> JetRecoEval::all_truth_showers(Jet* recojet)
0080 {
0081   if (_strict)
0082   {
0083     assert(recojet);
0084   }
0085   else if (!recojet)
0086   {
0087     ++_errors;
0088     return std::set<PHG4Shower*>();
0089   }
0090 
0091   if (_do_cache)
0092   {
0093     std::map<Jet*, std::set<PHG4Shower*> >::iterator iter =
0094         _cache_all_truth_showers.find(recojet);
0095     if (iter != _cache_all_truth_showers.end())
0096     {
0097       return iter->second;
0098     }
0099   }
0100 
0101   std::set<PHG4Shower*> truth_showers;
0102 
0103   // loop over all the jet constituents, backtrack each reco object to the
0104   // truth hits and combine with other consituents
0105   for (auto jter : recojet->get_comp_vec())
0106   {
0107       Jet::SRC source = jter.first;
0108       unsigned int index = jter.second;
0109 
0110     std::set<PHG4Shower*> new_showers;
0111 
0112     if (source == Jet::TRACK)
0113     {
0114       if (!_trackmap)
0115       {
0116         std::cout << PHWHERE << "ERROR: can't find SvtxTrackMap" << std::endl;
0117         exit(-1);
0118       }
0119 
0120       // SvtxTrack* track = _trackmap->get(index);
0121 
0122       // if (_strict) {assert(track);}
0123       // else if (!track) {++_errors; continue;}
0124 
0125       // new_showers = get_svtx_eval_stack()->get_track_eval()->all_truth_showers(track);
0126     }
0127 
0128     else if (source == Jet::CEMC_TOWER)
0129     {
0130       if (!_cemctowers)
0131       {
0132         std::cout << PHWHERE << "ERROR: can't find TOWER_CEMC" << std::endl;
0133         exit(-1);
0134       }
0135 
0136       RawTower* tower = _cemctowers->getTower(index);
0137 
0138       if (_strict)
0139       {
0140         assert(tower);
0141       }
0142       else if (!tower)
0143       {
0144         ++_errors;
0145         continue;
0146       }
0147 
0148       new_showers = get_cemc_eval_stack()->get_rawtower_eval()->all_truth_primary_showers(tower);
0149     }
0150     else if (source == Jet::CEMC_CLUSTER)
0151     {
0152       if (!_cemcclusters)
0153       {
0154         std::cout << PHWHERE << "ERROR: can't find CLUSTER_CEMC" << std::endl;
0155         exit(-1);
0156       }
0157 
0158       RawCluster* cluster = _cemcclusters->getCluster(index);
0159 
0160       if (_strict)
0161       {
0162         assert(cluster);
0163       }
0164       else if (!cluster)
0165       {
0166         ++_errors;
0167         continue;
0168       }
0169 
0170       new_showers = get_cemc_eval_stack()->get_rawcluster_eval()->all_truth_primary_showers(cluster);
0171     }
0172     else if (source == Jet::EEMC_TOWER)
0173     {
0174       if (!_eemctowers)
0175       {
0176         std::cout << PHWHERE << "ERROR: can't find TOWER_EEMC" << std::endl;
0177         exit(-1);
0178       }
0179 
0180       RawTower* tower = _eemctowers->getTower(index);
0181 
0182       if (_strict)
0183       {
0184         assert(tower);
0185       }
0186       else if (!tower)
0187       {
0188         ++_errors;
0189         continue;
0190       }
0191 
0192       new_showers = get_eemc_eval_stack()->get_rawtower_eval()->all_truth_primary_showers(tower);
0193     }
0194     else if (source == Jet::EEMC_CLUSTER)
0195     {
0196       if (!_eemcclusters)
0197       {
0198         std::cout << PHWHERE << "ERROR: can't find CLUSTER_EEMC" << std::endl;
0199         exit(-1);
0200       }
0201 
0202       RawCluster* cluster = _eemcclusters->getCluster(index);
0203 
0204       if (_strict)
0205       {
0206         assert(cluster);
0207       }
0208       else if (!cluster)
0209       {
0210         ++_errors;
0211         continue;
0212       }
0213 
0214       new_showers = get_eemc_eval_stack()->get_rawcluster_eval()->all_truth_primary_showers(cluster);
0215     }
0216     else if (source == Jet::HCALIN_TOWER)
0217     {
0218       if (!_hcalintowers)
0219       {
0220         std::cout << PHWHERE << "ERROR: can't find TOWER_HCALIN" << std::endl;
0221         exit(-1);
0222       }
0223 
0224       RawTower* tower = _hcalintowers->getTower(index);
0225 
0226       if (_strict)
0227       {
0228         assert(tower);
0229       }
0230       else if (!tower)
0231       {
0232         ++_errors;
0233         continue;
0234       }
0235 
0236       new_showers = get_hcalin_eval_stack()->get_rawtower_eval()->all_truth_primary_showers(tower);
0237     }
0238     else if (source == Jet::HCALIN_CLUSTER)
0239     {
0240       if (!_hcalinclusters)
0241       {
0242         std::cout << PHWHERE << "ERROR: can't find CLUSTER_HCALIN" << std::endl;
0243         exit(-1);
0244       }
0245 
0246       RawCluster* cluster = _hcalinclusters->getCluster(index);
0247 
0248       if (_strict)
0249       {
0250         assert(cluster);
0251       }
0252       else if (!cluster)
0253       {
0254         ++_errors;
0255         continue;
0256       }
0257 
0258       new_showers = get_hcalin_eval_stack()->get_rawcluster_eval()->all_truth_primary_showers(cluster);
0259     }
0260     else if (source == Jet::HCALOUT_TOWER)
0261     {
0262       if (!_hcalouttowers)
0263       {
0264         std::cout << PHWHERE << "ERROR: can't find TOWER_HCALOUT" << std::endl;
0265         exit(-1);
0266       }
0267 
0268       RawTower* tower = _hcalouttowers->getTower(index);
0269 
0270       if (_strict)
0271       {
0272         assert(tower);
0273       }
0274       else if (!tower)
0275       {
0276         ++_errors;
0277         continue;
0278       }
0279 
0280       new_showers = get_hcalout_eval_stack()->get_rawtower_eval()->all_truth_primary_showers(tower);
0281     }
0282     else if (source == Jet::HCALOUT_CLUSTER)
0283     {
0284       if (!_hcaloutclusters)
0285       {
0286         std::cout << PHWHERE << "ERROR: can't find CLUSTER_HCALOUT" << std::endl;
0287         exit(-1);
0288       }
0289 
0290       RawCluster* cluster = _hcaloutclusters->getCluster(index);
0291 
0292       if (_strict)
0293       {
0294         assert(cluster);
0295       }
0296       else if (!cluster)
0297       {
0298         ++_errors;
0299         continue;
0300       }
0301 
0302       new_showers = get_hcalout_eval_stack()->get_rawcluster_eval()->all_truth_primary_showers(cluster);
0303     }
0304     else if (source == Jet::FEMC_TOWER)
0305     {
0306       if (!_femctowers)
0307       {
0308         std::cout << PHWHERE << "ERROR: can't find TOWER_FEMC" << std::endl;
0309         exit(-1);
0310       }
0311 
0312       RawTower* tower = _femctowers->getTower(index);
0313 
0314       if (_strict)
0315       {
0316         assert(tower);
0317       }
0318       else if (!tower)
0319       {
0320         ++_errors;
0321         continue;
0322       }
0323 
0324       new_showers = get_femc_eval_stack()->get_rawtower_eval()->all_truth_primary_showers(tower);
0325     }
0326     else if (source == Jet::FEMC_CLUSTER)
0327     {
0328       if (!_femcclusters)
0329       {
0330         std::cout << PHWHERE << "ERROR: can't find CLUSTER_FEMC" << std::endl;
0331         exit(-1);
0332       }
0333 
0334       RawCluster* cluster = _femcclusters->getCluster(index);
0335 
0336       if (_strict)
0337       {
0338         assert(cluster);
0339       }
0340       else if (!cluster)
0341       {
0342         ++_errors;
0343         continue;
0344       }
0345 
0346       new_showers = get_femc_eval_stack()->get_rawcluster_eval()->all_truth_primary_showers(cluster);
0347     }
0348     else if (source == Jet::FHCAL_TOWER)
0349     {
0350       if (!_fhcaltowers)
0351       {
0352         std::cout << PHWHERE << "ERROR: can't find TOWER_FHCAL" << std::endl;
0353         exit(-1);
0354       }
0355 
0356       RawTower* tower = _fhcaltowers->getTower(index);
0357 
0358       if (_strict)
0359       {
0360         assert(tower);
0361       }
0362       else if (!tower)
0363       {
0364         ++_errors;
0365         continue;
0366       }
0367 
0368       new_showers = get_fhcal_eval_stack()->get_rawtower_eval()->all_truth_primary_showers(tower);
0369     }
0370     else if (source == Jet::FHCAL_CLUSTER)
0371     {
0372       if (!_fhcalclusters)
0373       {
0374         std::cout << PHWHERE << "ERROR: can't find CLUSTER_FHCAL" << std::endl;
0375         exit(-1);
0376       }
0377 
0378       RawCluster* cluster = _fhcalclusters->getCluster(index);
0379 
0380       if (_strict)
0381       {
0382         assert(cluster);
0383       }
0384       else if (!cluster)
0385       {
0386         ++_errors;
0387         continue;
0388       }
0389 
0390       new_showers = get_fhcal_eval_stack()->get_rawcluster_eval()->all_truth_primary_showers(cluster);
0391     }
0392 
0393     for (auto new_shower : new_showers)
0394     {
0395       truth_showers.insert(new_shower);
0396     }
0397   }
0398 
0399   if (_do_cache)
0400   {
0401     _cache_all_truth_showers.insert(std::make_pair(recojet, truth_showers));
0402   }
0403 
0404   return truth_showers;
0405 }
0406 
0407 std::set<PHG4Particle*> JetRecoEval::all_truth_particles(Jet* recojet)
0408 {
0409   if (_strict)
0410   {
0411     assert(recojet);
0412   }
0413   else if (!recojet)
0414   {
0415     ++_errors;
0416     return std::set<PHG4Particle*>();
0417   }
0418 
0419   if (_do_cache)
0420   {
0421     std::map<Jet*, std::set<PHG4Particle*> >::iterator iter =
0422         _cache_all_truth_particles.find(recojet);
0423     if (iter != _cache_all_truth_particles.end())
0424     {
0425       return iter->second;
0426     }
0427   }
0428 
0429   std::set<PHG4Particle*> truth_particles;
0430 
0431   // loop over all the jet constituents, backtrack each reco object to the
0432   // truth hits and combine with other consituents
0433   for (auto jter : recojet->get_comp_vec())
0434   {
0435       Jet::SRC source = jter.first;
0436       unsigned int index = jter.second;
0437 
0438     std::set<PHG4Particle*> new_particles;
0439 
0440     if (source == Jet::TRACK)
0441     {
0442       if (!_trackmap)
0443       {
0444         std::cout << PHWHERE << "ERROR: can't find TrackMap" << std::endl;
0445         exit(-1);
0446       }
0447 
0448       SvtxTrack* track = _trackmap->get(index);
0449 
0450       if (_strict)
0451       {
0452         assert(track);
0453       }
0454       else if (!track)
0455       {
0456         ++_errors;
0457         continue;
0458       }
0459 
0460       new_particles = get_svtx_eval_stack()->get_track_eval()->all_truth_particles(track);
0461     }
0462     else if (source == Jet::CEMC_TOWER)
0463     {
0464       if (!_cemctowers)
0465       {
0466         std::cout << PHWHERE << "ERROR: can't find TOWER_CEMC" << std::endl;
0467         exit(-1);
0468       }
0469 
0470       RawTower* tower = _cemctowers->getTower(index);
0471 
0472       if (_strict)
0473       {
0474         assert(tower);
0475       }
0476       else if (!tower)
0477       {
0478         ++_errors;
0479         continue;
0480       }
0481 
0482       new_particles = get_cemc_eval_stack()->get_rawtower_eval()->all_truth_primary_particles(tower);
0483     }
0484     else if (source == Jet::CEMC_CLUSTER)
0485     {
0486       if (!_cemcclusters)
0487       {
0488         std::cout << PHWHERE << "ERROR: can't find CLUSTER_CEMC" << std::endl;
0489         exit(-1);
0490       }
0491 
0492       RawCluster* cluster = _cemcclusters->getCluster(index);
0493 
0494       if (_strict)
0495       {
0496         assert(cluster);
0497       }
0498       else if (!cluster)
0499       {
0500         ++_errors;
0501         continue;
0502       }
0503 
0504       new_particles = get_cemc_eval_stack()->get_rawcluster_eval()->all_truth_primary_particles(cluster);
0505     }
0506     else if (source == Jet::EEMC_TOWER)
0507     {
0508       if (!_eemctowers)
0509       {
0510         std::cout << PHWHERE << "ERROR: can't find TOWER_EEMC" << std::endl;
0511         exit(-1);
0512       }
0513 
0514       RawTower* tower = _eemctowers->getTower(index);
0515 
0516       if (_strict)
0517       {
0518         assert(tower);
0519       }
0520       else if (!tower)
0521       {
0522         ++_errors;
0523         continue;
0524       }
0525 
0526       new_particles = get_eemc_eval_stack()->get_rawtower_eval()->all_truth_primary_particles(tower);
0527     }
0528     else if (source == Jet::EEMC_CLUSTER)
0529     {
0530       if (!_eemcclusters)
0531       {
0532         std::cout << PHWHERE << "ERROR: can't find CLUSTER_EEMC" << std::endl;
0533         exit(-1);
0534       }
0535 
0536       RawCluster* cluster = _eemcclusters->getCluster(index);
0537 
0538       if (_strict)
0539       {
0540         assert(cluster);
0541       }
0542       else if (!cluster)
0543       {
0544         ++_errors;
0545         continue;
0546       }
0547 
0548       new_particles = get_eemc_eval_stack()->get_rawcluster_eval()->all_truth_primary_particles(cluster);
0549     }
0550     else if (source == Jet::HCALIN_TOWER)
0551     {
0552       if (!_hcalintowers)
0553       {
0554         std::cout << PHWHERE << "ERROR: can't find TOWER_HCALIN" << std::endl;
0555         exit(-1);
0556       }
0557 
0558       RawTower* tower = _hcalintowers->getTower(index);
0559 
0560       if (_strict)
0561       {
0562         assert(tower);
0563       }
0564       else if (!tower)
0565       {
0566         ++_errors;
0567         continue;
0568       }
0569 
0570       new_particles = get_hcalin_eval_stack()->get_rawtower_eval()->all_truth_primary_particles(tower);
0571     }
0572     else if (source == Jet::HCALIN_CLUSTER)
0573     {
0574       if (!_hcalinclusters)
0575       {
0576         std::cout << PHWHERE << "ERROR: can't find CLUSTER_HCALIN" << std::endl;
0577         exit(-1);
0578       }
0579 
0580       RawCluster* cluster = _hcalinclusters->getCluster(index);
0581 
0582       if (_strict)
0583       {
0584         assert(cluster);
0585       }
0586       else if (!cluster)
0587       {
0588         ++_errors;
0589         continue;
0590       }
0591 
0592       new_particles = get_hcalin_eval_stack()->get_rawcluster_eval()->all_truth_primary_particles(cluster);
0593     }
0594     else if (source == Jet::HCALOUT_TOWER)
0595     {
0596       if (!_hcalouttowers)
0597       {
0598         std::cout << PHWHERE << "ERROR: can't find TOWER_HCALOUT" << std::endl;
0599         exit(-1);
0600       }
0601 
0602       RawTower* tower = _hcalouttowers->getTower(index);
0603 
0604       if (_strict)
0605       {
0606         assert(tower);
0607       }
0608       else if (!tower)
0609       {
0610         ++_errors;
0611         continue;
0612       }
0613 
0614       new_particles = get_hcalout_eval_stack()->get_rawtower_eval()->all_truth_primary_particles(tower);
0615     }
0616     else if (source == Jet::HCALOUT_CLUSTER)
0617     {
0618       if (!_hcaloutclusters)
0619       {
0620         std::cout << PHWHERE << "ERROR: can't find CLUSTER_HCALOUT" << std::endl;
0621         exit(-1);
0622       }
0623 
0624       RawCluster* cluster = _hcaloutclusters->getCluster(index);
0625 
0626       if (_strict)
0627       {
0628         assert(cluster);
0629       }
0630       else if (!cluster)
0631       {
0632         ++_errors;
0633         continue;
0634       }
0635 
0636       new_particles = get_hcalout_eval_stack()->get_rawcluster_eval()->all_truth_primary_particles(cluster);
0637     }
0638     else if (source == Jet::FEMC_TOWER)
0639     {
0640       if (!_femctowers)
0641       {
0642         std::cout << PHWHERE << "ERROR: can't find TOWER_FEMC" << std::endl;
0643         exit(-1);
0644       }
0645 
0646       RawTower* tower = _femctowers->getTower(index);
0647 
0648       if (_strict)
0649       {
0650         assert(tower);
0651       }
0652       else if (!tower)
0653       {
0654         ++_errors;
0655         continue;
0656       }
0657 
0658       new_particles = get_femc_eval_stack()->get_rawtower_eval()->all_truth_primary_particles(tower);
0659     }
0660     else if (source == Jet::FEMC_CLUSTER)
0661     {
0662       if (!_femcclusters)
0663       {
0664         std::cout << PHWHERE << "ERROR: can't find CLUSTER_FEMC" << std::endl;
0665         exit(-1);
0666       }
0667 
0668       RawCluster* cluster = _femcclusters->getCluster(index);
0669 
0670       if (_strict)
0671       {
0672         assert(cluster);
0673       }
0674       else if (!cluster)
0675       {
0676         ++_errors;
0677         continue;
0678       }
0679 
0680       new_particles = get_femc_eval_stack()->get_rawcluster_eval()->all_truth_primary_particles(cluster);
0681     }
0682     else if (source == Jet::FHCAL_TOWER)
0683     {
0684       if (!_fhcaltowers)
0685       {
0686         std::cout << PHWHERE << "ERROR: can't find TOWER_FHCAL" << std::endl;
0687         exit(-1);
0688       }
0689 
0690       RawTower* tower = _fhcaltowers->getTower(index);
0691 
0692       if (_strict)
0693       {
0694         assert(tower);
0695       }
0696       else if (!tower)
0697       {
0698         ++_errors;
0699         continue;
0700       }
0701 
0702       new_particles = get_fhcal_eval_stack()->get_rawtower_eval()->all_truth_primary_particles(tower);
0703     }
0704     else if (source == Jet::FHCAL_CLUSTER)
0705     {
0706       if (!_fhcalclusters)
0707       {
0708         std::cout << PHWHERE << "ERROR: can't find CLUSTER_FHCAL" << std::endl;
0709         exit(-1);
0710       }
0711 
0712       RawCluster* cluster = _fhcalclusters->getCluster(index);
0713 
0714       if (_strict)
0715       {
0716         assert(cluster);
0717       }
0718       else if (!cluster)
0719       {
0720         ++_errors;
0721         continue;
0722       }
0723 
0724       new_particles = get_fhcal_eval_stack()->get_rawcluster_eval()->all_truth_primary_particles(cluster);
0725     }
0726 
0727     for (auto new_particle : new_particles)
0728     {
0729       truth_particles.insert(new_particle);
0730     }
0731   }
0732 
0733   if (_do_cache)
0734   {
0735     _cache_all_truth_particles.insert(std::make_pair(recojet, truth_particles));
0736   }
0737 
0738   return truth_particles;
0739 }
0740 
0741 std::set<Jet*> JetRecoEval::all_truth_jets(Jet* recojet)
0742 {
0743   if (_strict)
0744   {
0745     assert(recojet);
0746   }
0747   else if (!recojet)
0748   {
0749     ++_errors;
0750     return std::set<Jet*>();
0751   }
0752 
0753   if (_do_cache)
0754   {
0755     std::map<Jet*, std::set<Jet*> >::iterator iter =
0756         _cache_all_truth_jets.find(recojet);
0757     if (iter != _cache_all_truth_jets.end())
0758     {
0759       return iter->second;
0760     }
0761   }
0762 
0763   std::set<Jet*> truth_jets;
0764 
0765   // get all truth particles (this can include muons and other truth excludes)...
0766   std::set<PHG4Particle*> particles = all_truth_particles(recojet);
0767 
0768   // backtrack from the truth particles to the truth jets...
0769   for (auto particle : particles)
0770   {
0771     if (_strict)
0772     {
0773       assert(particle);
0774     }
0775     else if (!particle)
0776     {
0777       ++_errors;
0778       continue;
0779     }
0780 
0781     Jet* truth_jet = _jettrutheval.get_truth_jet(particle);
0782     if (!truth_jet)
0783     {
0784       continue;
0785     }
0786 
0787     truth_jets.insert(truth_jet);
0788   }
0789 
0790   if (_do_cache)
0791   {
0792     _cache_all_truth_jets.insert(std::make_pair(recojet, truth_jets));
0793   }
0794 
0795   return truth_jets;
0796 }
0797 
0798 Jet* JetRecoEval::max_truth_jet_by_energy(Jet* recojet)
0799 {
0800   if (_strict)
0801   {
0802     assert(recojet);
0803   }
0804   else if (!recojet)
0805   {
0806     ++_errors;
0807     return nullptr;
0808   }
0809 
0810   if (_do_cache)
0811   {
0812     std::map<Jet*, Jet*>::iterator iter =
0813         _cache_max_truth_jet_by_energy.find(recojet);
0814     if (iter != _cache_max_truth_jet_by_energy.end())
0815     {
0816       return iter->second;
0817     }
0818   }
0819 
0820   Jet* truthjet = nullptr;
0821   float max_energy = FLT_MAX * -1.0;
0822 
0823   std::set<Jet*> truthjets = all_truth_jets(recojet);
0824   for (auto candidate : truthjets)
0825   {
0826     if (_strict)
0827     {
0828       assert(candidate);
0829     }
0830     else if (!candidate)
0831     {
0832       ++_errors;
0833       continue;
0834     }
0835 
0836     float energy = get_energy_contribution(recojet, candidate);
0837     if (energy > max_energy)
0838     {
0839       truthjet = candidate;
0840       max_energy = energy;
0841     }
0842   }
0843 
0844   if (_do_cache)
0845   {
0846     _cache_max_truth_jet_by_energy.insert(std::make_pair(recojet, truthjet));
0847   }
0848 
0849   return truthjet;
0850 }
0851 
0852 std::set<Jet*> JetRecoEval::all_jets_from(Jet* truthjet)
0853 {
0854   if (_strict)
0855   {
0856     assert(truthjet);
0857   }
0858   else if (!truthjet)
0859   {
0860     ++_errors;
0861     return std::set<Jet*>();
0862   }
0863 
0864   if (_do_cache)
0865   {
0866     std::map<Jet*, std::set<Jet*> >::iterator iter =
0867         _cache_all_jets_from.find(truthjet);
0868     if (iter != _cache_all_jets_from.end())
0869     {
0870       return iter->second;
0871     }
0872   }
0873 
0874   std::set<Jet*> recojets;
0875 
0876   // loop over all reco jets
0877   for (auto recojet : *_recojets)
0878   {
0879     /* Jet* recojet = _recojet.second; */
0880 
0881     // if this jet back tracks to the truth jet
0882     std::set<Jet*> truthcandidates = all_truth_jets(recojet);
0883     for (auto truthcandidate : truthcandidates)
0884     {
0885       if (_strict)
0886       {
0887         assert(truthcandidate);
0888       }
0889       else if (!truthcandidate)
0890       {
0891         ++_errors;
0892         continue;
0893       }
0894 
0895       if (truthcandidate->get_id() == truthjet->get_id())
0896       {
0897         recojets.insert(recojet);
0898       }
0899     }
0900   }
0901 
0902   if (_do_cache)
0903   {
0904     _cache_all_jets_from.insert(std::make_pair(truthjet, recojets));
0905   }
0906 
0907   return recojets;
0908 }
0909 
0910 Jet* JetRecoEval::best_jet_from(Jet* truthjet)
0911 {
0912 
0913   if (_strict)
0914   {
0915     assert(truthjet);
0916   }
0917   else if (!truthjet)
0918   {
0919     ++_errors;
0920     return nullptr;
0921   }
0922 
0923   if (_do_cache)
0924   {
0925     std::map<Jet*, Jet*>::iterator iter =
0926         _cache_best_jet_from.find(truthjet);
0927     if (iter != _cache_best_jet_from.end())
0928     {
0929       return iter->second;
0930     }
0931   }
0932 
0933   Jet* bestrecojet = nullptr;
0934   float max_energy = FLT_MAX * -1.0;
0935 
0936   std::set<Jet*> recojets = all_jets_from(truthjet);
0937   for (auto recojet : recojets)
0938   {
0939     if (_strict)
0940     {
0941       assert(recojet);
0942     }
0943     else if (!recojet)
0944     {
0945       ++_errors;
0946       continue;
0947     }
0948 
0949     float energy = get_energy_contribution(recojet, truthjet);
0950     if (energy > max_energy)
0951     {
0952       bestrecojet = recojet;
0953       max_energy = energy;
0954     }
0955   }
0956 
0957   if (_do_cache)
0958   {
0959     _cache_best_jet_from.insert(std::make_pair(truthjet, bestrecojet));
0960   }
0961   return bestrecojet;
0962 }
0963 
0964 Jet* JetRecoEval::unique_reco_jet_from_truth(Jet* truthjet)
0965 {
0966   if (_strict)
0967   {
0968     assert(truthjet);
0969   }
0970   else if (!truthjet)
0971   {
0972     ++_errors;
0973     return nullptr;
0974   }
0975 
0976   Jet* recojet = best_jet_from(truthjet);
0977 
0978   if (recojet)
0979   {
0980     Jet* back_matching = max_truth_jet_by_energy(recojet);
0981 
0982     if (back_matching->get_id() == truthjet->get_id())
0983     {
0984       return recojet;  // uniquely matched
0985     }
0986     else
0987     {
0988       return nullptr;
0989     }
0990   }
0991   else
0992   {
0993     return nullptr;
0994   }
0995 }
0996 
0997 Jet* JetRecoEval::unique_truth_jet_from_reco(Jet* recojet)
0998 {
0999   if (_strict)
1000   {
1001     assert(recojet);
1002   }
1003   else if (!recojet)
1004   {
1005     ++_errors;
1006     return nullptr;
1007   }
1008 
1009   Jet* truthjet = max_truth_jet_by_energy(recojet);
1010 
1011   if (truthjet)
1012   {
1013     Jet* back_matching = best_jet_from(truthjet);
1014 
1015     if (back_matching->get_id() == recojet->get_id())
1016     {
1017       return truthjet;  // uniquely matched
1018     }
1019     else
1020     {
1021       return nullptr;
1022     }
1023   }
1024   else
1025   {
1026     return nullptr;
1027   }
1028 }
1029 
1030 // overlap calculations
1031 float JetRecoEval::get_energy_contribution(Jet* recojet, Jet* truthjet)
1032 {
1033   if (_strict)
1034   {
1035     assert(recojet);
1036     assert(truthjet);
1037   }
1038   else if (!recojet || !truthjet)
1039   {
1040     ++_errors;
1041     return NAN;
1042   }
1043 
1044   if (_do_cache)
1045   {
1046     std::map<std::pair<Jet*, Jet*>, float>::iterator iter =
1047         _cache_get_energy_contribution.find(std::make_pair(recojet, truthjet));
1048     if (iter != _cache_get_energy_contribution.end())
1049     {
1050       return iter->second;
1051     }
1052   }
1053 
1054   float energy_contribution = 0.0;
1055 
1056   std::set<PHG4Particle*> truthjetcomp = get_truth_eval()->all_truth_particles(truthjet);
1057   // loop over all truthjet constituents
1058   for (auto truthparticle : truthjetcomp)
1059   {
1060     if (_strict)
1061     {
1062       assert(truthparticle);
1063     }
1064     else if (!truthparticle)
1065     {
1066       ++_errors;
1067       continue;
1068     }
1069 
1070     for (auto jter : recojet->get_comp_vec())
1071     {
1072       Jet::SRC source = jter.first;
1073       unsigned int index = jter.second;
1074 
1075       float energy = 0.0;
1076 
1077       if (source == Jet::TRACK)
1078       {
1079         SvtxTrack* track = _trackmap->get(index);
1080 
1081         if (_strict)
1082         {
1083           assert(track);
1084         }
1085         else if (!track)
1086         {
1087           ++_errors;
1088           continue;
1089         }
1090 
1091         PHG4Particle* maxtruthparticle = get_svtx_eval_stack()->get_track_eval()->max_truth_particle_by_nclusters(track);
1092 
1093         if (maxtruthparticle == nullptr)
1094         {
1095           // in extreme rare cases, noise hits can make a track with no maxtruthparticle matched
1096           energy = 0;
1097         }
1098         else if (maxtruthparticle->get_track_id() == truthparticle->get_track_id())
1099         {
1100           energy = track->get_p();
1101         }
1102       }
1103       else if (source == Jet::CEMC_TOWER)
1104       {
1105         RawTower* tower = _cemctowers->getTower(index);
1106 
1107         if (_strict)
1108         {
1109           assert(tower);
1110         }
1111         else if (!tower)
1112         {
1113           ++_errors;
1114           continue;
1115         }
1116 
1117         energy = get_cemc_eval_stack()->get_rawtower_eval()->get_energy_contribution(tower, truthparticle);
1118       }
1119       else if (source == Jet::CEMC_CLUSTER)
1120       {
1121         RawCluster* cluster = _cemcclusters->getCluster(index);
1122 
1123         if (_strict)
1124         {
1125           assert(cluster);
1126         }
1127         else if (!cluster)
1128         {
1129           ++_errors;
1130           continue;
1131         }
1132 
1133         energy = get_cemc_eval_stack()->get_rawcluster_eval()->get_energy_contribution(cluster, truthparticle);
1134       }
1135       else if (source == Jet::EEMC_TOWER)
1136       {
1137         RawTower* tower = _eemctowers->getTower(index);
1138 
1139         if (_strict)
1140         {
1141           assert(tower);
1142         }
1143         else if (!tower)
1144         {
1145           ++_errors;
1146           continue;
1147         }
1148 
1149         energy = get_eemc_eval_stack()->get_rawtower_eval()->get_energy_contribution(tower, truthparticle);
1150       }
1151       else if (source == Jet::EEMC_CLUSTER)
1152       {
1153         RawCluster* cluster = _eemcclusters->getCluster(index);
1154 
1155         if (_strict)
1156         {
1157           assert(cluster);
1158         }
1159         else if (!cluster)
1160         {
1161           ++_errors;
1162           continue;
1163         }
1164 
1165         energy = get_eemc_eval_stack()->get_rawcluster_eval()->get_energy_contribution(cluster, truthparticle);
1166       }
1167       else if (source == Jet::HCALIN_TOWER)
1168       {
1169         RawTower* tower = _hcalintowers->getTower(index);
1170 
1171         if (_strict)
1172         {
1173           assert(tower);
1174         }
1175         else if (!tower)
1176         {
1177           ++_errors;
1178           continue;
1179         }
1180 
1181         energy = get_hcalin_eval_stack()->get_rawtower_eval()->get_energy_contribution(tower, truthparticle);
1182       }
1183       else if (source == Jet::HCALIN_CLUSTER)
1184       {
1185         RawCluster* cluster = _hcalinclusters->getCluster(index);
1186 
1187         if (_strict)
1188         {
1189           assert(cluster);
1190         }
1191         else if (!cluster)
1192         {
1193           ++_errors;
1194           continue;
1195         }
1196 
1197         energy = get_hcalin_eval_stack()->get_rawcluster_eval()->get_energy_contribution(cluster, truthparticle);
1198       }
1199       else if (source == Jet::HCALOUT_TOWER)
1200       {
1201         RawTower* tower = _hcalouttowers->getTower(index);
1202 
1203         if (_strict)
1204         {
1205           assert(tower);
1206         }
1207         else if (!tower)
1208         {
1209           ++_errors;
1210           continue;
1211         }
1212 
1213         energy = get_hcalout_eval_stack()->get_rawtower_eval()->get_energy_contribution(tower, truthparticle);
1214       }
1215       else if (source == Jet::HCALOUT_CLUSTER)
1216       {
1217         RawCluster* cluster = _hcaloutclusters->getCluster(index);
1218 
1219         if (_strict)
1220         {
1221           assert(cluster);
1222         }
1223         else if (!cluster)
1224         {
1225           ++_errors;
1226           continue;
1227         }
1228 
1229         energy = get_hcalout_eval_stack()->get_rawcluster_eval()->get_energy_contribution(cluster, truthparticle);
1230       }
1231       else if (source == Jet::FEMC_TOWER)
1232       {
1233         RawTower* tower = _femctowers->getTower(index);
1234 
1235         if (_strict)
1236         {
1237           assert(tower);
1238         }
1239         else if (!tower)
1240         {
1241           ++_errors;
1242           continue;
1243         }
1244 
1245         energy = get_femc_eval_stack()->get_rawtower_eval()->get_energy_contribution(tower, truthparticle);
1246       }
1247       else if (source == Jet::FEMC_CLUSTER)
1248       {
1249         RawCluster* cluster = _femcclusters->getCluster(index);
1250 
1251         if (_strict)
1252         {
1253           assert(cluster);
1254         }
1255         else if (!cluster)
1256         {
1257           ++_errors;
1258           continue;
1259         }
1260 
1261         energy = get_femc_eval_stack()->get_rawcluster_eval()->get_energy_contribution(cluster, truthparticle);
1262       }
1263       else if (source == Jet::FHCAL_TOWER)
1264       {
1265         RawTower* tower = _fhcaltowers->getTower(index);
1266 
1267         if (_strict)
1268         {
1269           assert(tower);
1270         }
1271         else if (!tower)
1272         {
1273           ++_errors;
1274           continue;
1275         }
1276 
1277         energy = get_fhcal_eval_stack()->get_rawtower_eval()->get_energy_contribution(tower, truthparticle);
1278       }
1279       else if (source == Jet::FHCAL_CLUSTER)
1280       {
1281         RawCluster* cluster = _fhcalclusters->getCluster(index);
1282 
1283         if (_strict)
1284         {
1285           assert(cluster);
1286         }
1287         else if (!cluster)
1288         {
1289           ++_errors;
1290           continue;
1291         }
1292 
1293         energy = get_fhcal_eval_stack()->get_rawcluster_eval()->get_energy_contribution(cluster, truthparticle);
1294       }
1295 
1296       energy_contribution += energy;
1297     }
1298   }
1299 
1300   if (_do_cache)
1301   {
1302     _cache_get_energy_contribution.insert(std::make_pair(std::make_pair(recojet, truthjet), energy_contribution));
1303   }
1304 
1305   return energy_contribution;
1306 }
1307 
1308 // overlap calculations
1309 float JetRecoEval::get_energy_contribution(Jet* recojet, Jet::SRC src)
1310 {
1311   if (_strict)
1312   {
1313     assert(recojet);
1314   }
1315   else if (!recojet)
1316   {
1317     ++_errors;
1318     return NAN;
1319   }
1320 
1321   if (_do_cache)
1322   {
1323     std::map<std::pair<Jet*, Jet::SRC>, float>::iterator iter =
1324         _cache_get_energy_contribution_src.find(std::make_pair(recojet, src));
1325     if (iter != _cache_get_energy_contribution_src.end())
1326     {
1327       return iter->second;
1328     }
1329   }
1330   
1331   float energy = 0.0;
1332 
1333   // loop over all recojet constituents
1334   for (Jet::ITER_comp_vec jter = recojet->comp_begin(src);
1335        jter != recojet->comp_end(src);
1336        ++jter)
1337   {
1338     Jet::SRC source = jter->first;
1339     assert(source == src);  // jet container consistency check
1340     unsigned int index = jter->second;
1341 
1342     if (source == Jet::TRACK)
1343     {
1344       SvtxTrack* track = _trackmap->get(index);
1345       energy += track->get_p();
1346     }
1347     else if (source == Jet::CEMC_TOWER)
1348     {
1349       RawTower* tower = _cemctowers->getTower(index);
1350 
1351       if (_strict)
1352       {
1353         assert(tower);
1354       }
1355       else if (!tower)
1356       {
1357         ++_errors;
1358         continue;
1359       }
1360 
1361       energy += tower->get_energy();
1362     }
1363     else if (source == Jet::CEMC_CLUSTER)
1364     {
1365       RawCluster* cluster = _cemcclusters->getCluster(index);
1366 
1367       if (_strict)
1368       {
1369         assert(cluster);
1370       }
1371       else if (!cluster)
1372       {
1373         ++_errors;
1374         continue;
1375       }
1376 
1377       energy += cluster->get_energy();
1378     }
1379     else if (source == Jet::EEMC_TOWER)
1380     {
1381       RawTower* tower = _eemctowers->getTower(index);
1382 
1383       if (_strict)
1384       {
1385         assert(tower);
1386       }
1387       else if (!tower)
1388       {
1389         ++_errors;
1390         continue;
1391       }
1392 
1393       energy += tower->get_energy();
1394     }
1395     else if (source == Jet::EEMC_CLUSTER)
1396     {
1397       RawCluster* cluster = _eemcclusters->getCluster(index);
1398 
1399       if (_strict)
1400       {
1401         assert(cluster);
1402       }
1403       else if (!cluster)
1404       {
1405         ++_errors;
1406         continue;
1407       }
1408 
1409       energy += cluster->get_energy();
1410     }
1411     else if (source == Jet::HCALIN_TOWER)
1412     {
1413       RawTower* tower = _hcalintowers->getTower(index);
1414 
1415       if (_strict)
1416       {
1417         assert(tower);
1418       }
1419       else if (!tower)
1420       {
1421         ++_errors;
1422         continue;
1423       }
1424 
1425       energy += tower->get_energy();
1426     }
1427     else if (source == Jet::HCALIN_CLUSTER)
1428     {
1429       RawCluster* cluster = _hcalinclusters->getCluster(index);
1430 
1431       if (_strict)
1432       {
1433         assert(cluster);
1434       }
1435       else if (!cluster)
1436       {
1437         ++_errors;
1438         continue;
1439       }
1440 
1441       energy += cluster->get_energy();
1442     }
1443     else if (source == Jet::HCALOUT_TOWER)
1444     {
1445       RawTower* tower = _hcalouttowers->getTower(index);
1446 
1447       if (_strict)
1448       {
1449         assert(tower);
1450       }
1451       else if (!tower)
1452       {
1453         ++_errors;
1454         continue;
1455       }
1456 
1457       energy += tower->get_energy();
1458     }
1459     else if (source == Jet::HCALOUT_CLUSTER)
1460     {
1461       RawCluster* cluster = _hcaloutclusters->getCluster(index);
1462 
1463       if (_strict)
1464       {
1465         assert(cluster);
1466       }
1467       else if (!cluster)
1468       {
1469         ++_errors;
1470         continue;
1471       }
1472 
1473       energy += cluster->get_energy();
1474     }
1475     else if (source == Jet::FEMC_TOWER)
1476     {
1477       RawTower* tower = _femctowers->getTower(index);
1478 
1479       if (_strict)
1480       {
1481         assert(tower);
1482       }
1483       else if (!tower)
1484       {
1485         ++_errors;
1486         continue;
1487       }
1488 
1489       energy += tower->get_energy();
1490     }
1491     else if (source == Jet::FEMC_CLUSTER)
1492     {
1493       RawCluster* cluster = _femcclusters->getCluster(index);
1494 
1495       if (_strict)
1496       {
1497         assert(cluster);
1498       }
1499       else if (!cluster)
1500       {
1501         ++_errors;
1502         continue;
1503       }
1504 
1505       energy += cluster->get_energy();
1506     }
1507     else if (source == Jet::FHCAL_TOWER)
1508     {
1509       RawTower* tower = _fhcaltowers->getTower(index);
1510 
1511       if (_strict)
1512       {
1513         assert(tower);
1514       }
1515       else if (!tower)
1516       {
1517         ++_errors;
1518         continue;
1519       }
1520 
1521       energy += tower->get_energy();
1522     }
1523     else if (source == Jet::FHCAL_CLUSTER)
1524     {
1525       RawCluster* cluster = _fhcalclusters->getCluster(index);
1526 
1527       if (_strict)
1528       {
1529         assert(cluster);
1530       }
1531       else if (!cluster)
1532       {
1533         ++_errors;
1534         continue;
1535       }
1536 
1537       energy += cluster->get_energy();
1538     }  // else if (source == Jet::FHCAL_CLUSTER)
1539 
1540   }  // for (Jet::ConstIter jter = recojet->lower_bound_comp(src);
1541 
1542   if (_do_cache)
1543   {
1544     _cache_get_energy_contribution_src.insert(std::make_pair(std::make_pair(recojet, src), energy));
1545   }
1546 
1547   return energy;
1548 }
1549 
1550 std::set<PHG4Hit*> JetRecoEval::all_truth_hits(Jet* recojet)
1551 {
1552   if (_strict)
1553   {
1554     assert(recojet);
1555   }
1556   else if (!recojet)
1557   {
1558     ++_errors;
1559     return std::set<PHG4Hit*>();
1560   }
1561 
1562   if (_do_cache)
1563   {
1564     std::map<Jet*, std::set<PHG4Hit*> >::iterator iter =
1565         _cache_all_truth_hits.find(recojet);
1566     if (iter != _cache_all_truth_hits.end())
1567     {
1568       return iter->second;
1569     }
1570   }
1571 
1572   std::set<PHG4Hit*> truth_hits;
1573 
1574   // loop over all the jet constituents, backtrack each reco object to the
1575   // truth hits and combine with other consituents
1576 
1577   for (auto jter : recojet->get_comp_vec())
1578   {
1579       Jet::SRC source = jter.first;
1580       unsigned int index = jter.second;
1581 
1582     std::set<PHG4Hit*> new_hits;
1583 
1584     if (source == Jet::TRACK)
1585     {
1586       if (!_trackmap)
1587       {
1588         std::cout << PHWHERE << "ERROR: can't find SvtxTrackMap" << std::endl;
1589         exit(-1);
1590       }
1591 
1592       SvtxTrack* track = _trackmap->get(index);
1593 
1594       if (_strict)
1595       {
1596         assert(track);
1597       }
1598       else if (!track)
1599       {
1600         ++_errors;
1601         continue;
1602       }
1603 
1604       new_hits = get_svtx_eval_stack()->get_track_eval()->all_truth_hits(track);
1605     }
1606     else if (source == Jet::CEMC_TOWER)
1607     {
1608       if (!_cemctowers)
1609       {
1610         std::cout << PHWHERE << "ERROR: can't find TOWER_CEMC" << std::endl;
1611         exit(-1);
1612       }
1613 
1614       RawTower* tower = _cemctowers->getTower(index);
1615 
1616       if (_strict)
1617       {
1618         assert(tower);
1619       }
1620       else if (!tower)
1621       {
1622         ++_errors;
1623         continue;
1624       }
1625 
1626       new_hits = get_cemc_eval_stack()->get_rawtower_eval()->all_truth_hits(tower);
1627     }
1628     else if (source == Jet::CEMC_CLUSTER)
1629     {
1630       if (!_cemcclusters)
1631       {
1632         std::cout << PHWHERE << "ERROR: can't find CLUSTER_CEMC" << std::endl;
1633         exit(-1);
1634       }
1635 
1636       RawCluster* cluster = _cemcclusters->getCluster(index);
1637 
1638       if (_strict)
1639       {
1640         assert(cluster);
1641       }
1642       else if (!cluster)
1643       {
1644         ++_errors;
1645         continue;
1646       }
1647 
1648       new_hits = get_cemc_eval_stack()->get_rawcluster_eval()->all_truth_hits(cluster);
1649     }
1650     else if (source == Jet::EEMC_TOWER)
1651     {
1652       if (!_eemctowers)
1653       {
1654         std::cout << PHWHERE << "ERROR: can't find TOWER_EEMC" << std::endl;
1655         exit(-1);
1656       }
1657 
1658       RawTower* tower = _eemctowers->getTower(index);
1659 
1660       if (_strict)
1661       {
1662         assert(tower);
1663       }
1664       else if (!tower)
1665       {
1666         ++_errors;
1667         continue;
1668       }
1669 
1670       new_hits = get_eemc_eval_stack()->get_rawtower_eval()->all_truth_hits(tower);
1671     }
1672     else if (source == Jet::EEMC_CLUSTER)
1673     {
1674       if (!_eemcclusters)
1675       {
1676         std::cout << PHWHERE << "ERROR: can't find CLUSTER_EEMC" << std::endl;
1677         exit(-1);
1678       }
1679 
1680       RawCluster* cluster = _eemcclusters->getCluster(index);
1681 
1682       if (_strict)
1683       {
1684         assert(cluster);
1685       }
1686       else if (!cluster)
1687       {
1688         ++_errors;
1689         continue;
1690       }
1691 
1692       new_hits = get_eemc_eval_stack()->get_rawcluster_eval()->all_truth_hits(cluster);
1693     }
1694     else if (source == Jet::HCALIN_TOWER)
1695     {
1696       if (!_hcalintowers)
1697       {
1698         std::cout << PHWHERE << "ERROR: can't find TOWER_HCALIN" << std::endl;
1699         exit(-1);
1700       }
1701 
1702       RawTower* tower = _hcalintowers->getTower(index);
1703 
1704       if (_strict)
1705       {
1706         assert(tower);
1707       }
1708       else if (!tower)
1709       {
1710         ++_errors;
1711         continue;
1712       }
1713 
1714       new_hits = get_hcalin_eval_stack()->get_rawtower_eval()->all_truth_hits(tower);
1715     }
1716     else if (source == Jet::HCALIN_CLUSTER)
1717     {
1718       if (!_hcalinclusters)
1719       {
1720         std::cout << PHWHERE << "ERROR: can't find CLUSTER_HCALIN" << std::endl;
1721         exit(-1);
1722       }
1723 
1724       RawCluster* cluster = _hcalinclusters->getCluster(index);
1725 
1726       if (_strict)
1727       {
1728         assert(cluster);
1729       }
1730       else if (!cluster)
1731       {
1732         ++_errors;
1733         continue;
1734       }
1735 
1736       new_hits = get_hcalin_eval_stack()->get_rawcluster_eval()->all_truth_hits(cluster);
1737     }
1738     else if (source == Jet::HCALOUT_TOWER)
1739     {
1740       if (!_hcalouttowers)
1741       {
1742         std::cout << PHWHERE << "ERROR: can't find TOWER_HCALOUT" << std::endl;
1743         exit(-1);
1744       }
1745 
1746       RawTower* tower = _hcalouttowers->getTower(index);
1747 
1748       if (_strict)
1749       {
1750         assert(tower);
1751       }
1752       else if (!tower)
1753       {
1754         ++_errors;
1755         continue;
1756       }
1757 
1758       new_hits = get_hcalout_eval_stack()->get_rawtower_eval()->all_truth_hits(tower);
1759     }
1760     else if (source == Jet::HCALOUT_CLUSTER)
1761     {
1762       if (!_hcaloutclusters)
1763       {
1764         std::cout << PHWHERE << "ERROR: can't find CLUSTER_HCALOUT" << std::endl;
1765         exit(-1);
1766       }
1767 
1768       RawCluster* cluster = _hcaloutclusters->getCluster(index);
1769 
1770       if (_strict)
1771       {
1772         assert(cluster);
1773       }
1774       else if (!cluster)
1775       {
1776         ++_errors;
1777         continue;
1778       }
1779 
1780       new_hits = get_hcalout_eval_stack()->get_rawcluster_eval()->all_truth_hits(cluster);
1781     }
1782     else if (source == Jet::FEMC_TOWER)
1783     {
1784       if (!_femctowers)
1785       {
1786         std::cout << PHWHERE << "ERROR: can't find TOWER_FEMC" << std::endl;
1787         exit(-1);
1788       }
1789 
1790       RawTower* tower = _femctowers->getTower(index);
1791 
1792       if (_strict)
1793       {
1794         assert(tower);
1795       }
1796       else if (!tower)
1797       {
1798         ++_errors;
1799         continue;
1800       }
1801 
1802       new_hits = get_femc_eval_stack()->get_rawtower_eval()->all_truth_hits(tower);
1803     }
1804     else if (source == Jet::FEMC_CLUSTER)
1805     {
1806       if (!_femcclusters)
1807       {
1808         std::cout << PHWHERE << "ERROR: can't find CLUSTER_FEMC" << std::endl;
1809         exit(-1);
1810       }
1811 
1812       RawCluster* cluster = _femcclusters->getCluster(index);
1813 
1814       if (_strict)
1815       {
1816         assert(cluster);
1817       }
1818       else if (!cluster)
1819       {
1820         ++_errors;
1821         continue;
1822       }
1823 
1824       new_hits = get_femc_eval_stack()->get_rawcluster_eval()->all_truth_hits(cluster);
1825     }
1826     else if (source == Jet::FHCAL_TOWER)
1827     {
1828       if (!_fhcaltowers)
1829       {
1830         std::cout << PHWHERE << "ERROR: can't find TOWER_FHCAL" << std::endl;
1831         exit(-1);
1832       }
1833 
1834       RawTower* tower = _fhcaltowers->getTower(index);
1835 
1836       if (_strict)
1837       {
1838         assert(tower);
1839       }
1840       else if (!tower)
1841       {
1842         ++_errors;
1843         continue;
1844       }
1845 
1846       new_hits = get_fhcal_eval_stack()->get_rawtower_eval()->all_truth_hits(tower);
1847     }
1848     else if (source == Jet::FHCAL_CLUSTER)
1849     {
1850       if (!_fhcalclusters)
1851       {
1852         std::cout << PHWHERE << "ERROR: can't find CLUSTER_FHCAL" << std::endl;
1853         exit(-1);
1854       }
1855 
1856       RawCluster* cluster = _fhcalclusters->getCluster(index);
1857 
1858       if (_strict)
1859       {
1860         assert(cluster);
1861       }
1862       else if (!cluster)
1863       {
1864         ++_errors;
1865         continue;
1866       }
1867 
1868       new_hits = get_fhcal_eval_stack()->get_rawcluster_eval()->all_truth_hits(cluster);
1869     }
1870 
1871     for (auto new_hit : new_hits)
1872     {
1873       truth_hits.insert(new_hit);
1874     }
1875   }
1876 
1877   if (_do_cache)
1878   {
1879     _cache_all_truth_hits.insert(std::make_pair(recojet, truth_hits));
1880   }
1881 
1882   return truth_hits;
1883 }
1884 
1885 void JetRecoEval::get_node_pointers(PHCompositeNode* topNode)
1886 {
1887   // need things off of the DST...
1888   _recojets = findNode::getClass<JetContainer>(topNode, _recojetname.c_str());
1889   if (!_recojets)
1890   {
1891     std::cout << PHWHERE << " ERROR: Can't find " << _recojetname << std::endl;
1892     exit(-1);
1893   }
1894 
1895   _truthjets = findNode::getClass<JetContainer>(topNode, _truthjetname.c_str());
1896   if (!_truthjets)
1897   {
1898     std::cout << PHWHERE << " ERROR: Can't find " << _truthjetname << std::endl;
1899     exit(-1);
1900   }
1901 
1902   _trackmap = findNode::getClass<SvtxTrackMap>(topNode, m_TrackNodeName);
1903   if (!_trackmap)
1904   {
1905     _trackmap = findNode::getClass<SvtxTrackMap>(topNode, "TrackMap");
1906   }
1907   _cemctowers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_CEMC");
1908   _hcalintowers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALIN");
1909   _hcalouttowers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALOUT");
1910   _femctowers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_FEMC");
1911   _fhcaltowers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_FHCAL");
1912   _eemctowers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_EEMC");
1913   _cemcclusters = findNode::getClass<RawClusterContainer>(topNode, "CLUSTER_CEMC");
1914   _hcalinclusters = findNode::getClass<RawClusterContainer>(topNode, "CLUSTER_HCALIN");
1915   _hcaloutclusters = findNode::getClass<RawClusterContainer>(topNode, "CLUSTER_HCALOUT");
1916   _femcclusters = findNode::getClass<RawClusterContainer>(topNode, "CLUSTER_FEMC");
1917   _fhcalclusters = findNode::getClass<RawClusterContainer>(topNode, "CLUSTER_FHCAL");
1918   _eemcclusters = findNode::getClass<RawClusterContainer>(topNode, "CLUSTER_EEMC");
1919 
1920   return;
1921 }
[ 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!