File indexing completed on 2026-10-15 08:21:04
0001 #include "ParticleFlowEval.h"
0002 #include "CaloRawClusterEval.h"
0003 #include "CaloTruthEval.h"
0004 #include "SvtxTrackEval.h"
0005
0006 #include <calobase/RawClusterv1.h>
0007 #include <calobase/RawClusterContainer.h>
0008 #include <calobase/RawTowerContainer.h>
0009 #include <calobase/TowerInfo.h>
0010 #include <calobase/TowerInfoContainer.h>
0011 #include <calobase/TowerInfoDefs.h>
0012
0013 #include <g4main/PHG4Particlev2.h>
0014 #include <g4main/PHG4Shower.h>
0015
0016 #include <particleflowreco/ParticleFlowElement.h>
0017
0018 #include <phool/getClass.h>
0019
0020 #include <tr1/functional>
0021
0022 #include <cassert>
0023 #include <cfloat>
0024 #include <climits>
0025 #include <cmath>
0026 #include <iostream>
0027 #include <map>
0028 #include <set>
0029 #include <string>
0030
0031 class RawTower;
0032
0033 ParticleFlowEval::ParticleFlowEval(PHCompositeNode* topNode)
0034 : _clustereval_emcal(topNode, "CEMC")
0035 , _clustereval_ihcal(topNode, "HCALIN")
0036 , _clustereval_ohcal(topNode, "HCALOUT")
0037 , _trackeval(topNode)
0038 , _track(nullptr)
0039 {
0040 _topo_ihcal = new RawClusterv1();
0041 _topo_ohcal = new RawClusterv1();
0042 _topo_emcal.clear();
0043
0044
0045
0046
0047
0048 }
0049
0050 ParticleFlowEval::~ParticleFlowEval()
0051 {
0052 delete _topo_ihcal;
0053 delete _topo_ohcal;
0054 }
0055
0056 void ParticleFlowEval::next_event(PHCompositeNode* topNode)
0057 {
0058 _clustereval_emcal.next_event(topNode);
0059 _clustereval_ihcal.next_event(topNode);
0060 _clustereval_ohcal.next_event(topNode);
0061 _trackeval.next_event(topNode);
0062 }
0063
0064 void ParticleFlowEval::SeparateHCals(RawCluster *hcluster)
0065 {
0066 RawCluster::TowerConstRange begin_end = hcluster->get_towers();
0067 for(RawCluster::TowerConstIterator towerIt = begin_end.first; towerIt != begin_end.second; ++towerIt)
0068 {
0069 if(RawTowerDefs::decode_caloid(towerIt->first) == RawTowerDefs::CalorimeterId::HCALOUT)
0070 {
0071 _topo_ohcal->addTower(towerIt->first, towerIt->second);
0072 }
0073 else
0074 {
0075 _topo_ihcal->addTower(towerIt->first, towerIt->second);
0076 }
0077 }
0078 }
0079
0080 PHG4Particle* ParticleFlowEval::get_truth_particle_flow(ParticleFlowElement* pflow)
0081 {
0082 _pflow = pflow;
0083
0084 PHG4Particle *g4particle = nullptr;
0085
0086 _track = _pflow->get_track();
0087 _topo_emcal = _pflow->get_eclusters();
0088
0089 RawCluster *hcal_cluster = _pflow->get_hcluster();
0090 if(hcal_cluster)
0091 {
0092 SeparateHCals(hcal_cluster);
0093 }
0094
0095 if(_pflow->get_type() == ParticleFlowElement::PFLOWTYPE::UNMATCHED_CHARGED_HADRON)
0096 {
0097 g4particle = _trackeval.max_truth_particle_by_nclusters(_track);
0098 }
0099
0100
0101 if(_pflow->get_type() == ParticleFlowElement::PFLOWTYPE::MATCHED_CHARGED_HADRON)
0102 {
0103 g4particle = _trackeval.max_truth_particle_by_nclusters(_track);
0104 }
0105
0106 if(_pflow->get_type() == ParticleFlowElement::PFLOWTYPE::UNMATCHED_EM_PARTICLE)
0107 {
0108 g4particle = max_truth_primary_particle_by_energy_emcal();
0109 }
0110
0111 if(_pflow->get_type() == ParticleFlowElement::PFLOWTYPE::UNMATCHED_NEUTRAL_HADRON)
0112 {
0113 g4particle = max_truth_primary_particle_by_energy_hcal();
0114 }
0115
0116 if(_pflow->get_type() == ParticleFlowElement::PFLOWTYPE::LEFTOVER_EM_PARTICLE)
0117 {
0118 std::set<PHG4Particle*> primaries = all_truth_primary_particles_emcal();
0119
0120 primaries.merge(all_truth_primary_particles_hcal());
0121
0122 std::vector<SvtxTrack*> assoc_tracks = _pflow->get_assoc_tracks();
0123
0124 for(std::vector<SvtxTrack*>::iterator itr = assoc_tracks.begin(); itr != assoc_tracks.end(); itr++)
0125 {
0126 PHG4Particle *primary_track = _trackeval.max_truth_particle_by_nclusters((*itr));
0127 primaries.erase(primary_track);
0128 }
0129
0130 std::map<int, std::pair<float, PHG4Particle*>> primary_energy;
0131
0132 for(std::set<PHG4Particle*>::iterator itr = primaries.begin(); itr != primaries.end(); itr++)
0133 {
0134 if(primary_energy.find((*itr)->get_track_id()) != primary_energy.end())
0135 {
0136 for(unsigned int i = 0; i < _topo_emcal.size(); i++)
0137 {
0138 primary_energy[(*itr)->get_track_id()].first += _clustereval_emcal.get_energy_contribution(_topo_emcal.at(i), (*itr));
0139 }
0140 primary_energy[(*itr)->get_track_id()].first += _clustereval_ihcal.get_energy_contribution(_topo_ihcal, (*itr));
0141 primary_energy[(*itr)->get_track_id()].first += _clustereval_ohcal.get_energy_contribution(_topo_ohcal, (*itr));
0142 }
0143 else
0144 {
0145 primary_energy[(*itr)->get_track_id()].first = 0.;
0146 for(unsigned int i = 0; i < _topo_emcal.size(); i++)
0147 {
0148 primary_energy[(*itr)->get_track_id()].first += _clustereval_emcal.get_energy_contribution(_topo_emcal.at(i), (*itr));
0149 }
0150 primary_energy[(*itr)->get_track_id()].first += _clustereval_ihcal.get_energy_contribution(_topo_ihcal, (*itr));
0151 primary_energy[(*itr)->get_track_id()].first += _clustereval_ohcal.get_energy_contribution(_topo_ohcal, (*itr));
0152 primary_energy[(*itr)->get_track_id()].second = (*itr);
0153 }
0154 }
0155
0156 int max_index = -1;
0157 float max_energy = 0.;
0158
0159 for(std::map<int, std::pair<float, PHG4Particle*>>::iterator itr = primary_energy.begin(); itr != primary_energy.end(); ++itr)
0160 {
0161 if(itr->second.first > max_energy)
0162 {
0163 max_energy = itr->second.first;
0164 max_index = itr->first;
0165 }
0166 }
0167
0168 g4particle = primary_energy[max_index].second;
0169 }
0170
0171 return g4particle;
0172
0173 }
0174
0175 std::set<PHG4Shower*> ParticleFlowEval::all_truth_primary_showers_emcal()
0176 {
0177 std::set<PHG4Shower*> all_showers;
0178 for(unsigned int i = 0; i < _topo_emcal.size(); i++)
0179 {
0180 std::set<PHG4Shower*> showers = _clustereval_emcal.all_truth_primary_showers(_topo_emcal.at(i));
0181
0182 all_showers.merge(showers);
0183 }
0184 return all_showers;
0185 }
0186
0187 std::set<PHG4Shower*> ParticleFlowEval::all_truth_primary_showers_emcal(RawCluster *cluster)
0188 {
0189 if(std::find(_topo_emcal.begin(), _topo_emcal.end(), cluster) != _topo_emcal.end())
0190 {
0191 return _clustereval_emcal.all_truth_primary_showers(cluster);
0192 }
0193 else
0194 {
0195 return std::set<PHG4Shower*>();
0196 }
0197
0198 }
0199
0200 std::set<PHG4Shower*> ParticleFlowEval::all_truth_primary_showers_hcal()
0201 {
0202 std::set<PHG4Shower*> all_showers = all_truth_primary_showers_ihcal();
0203
0204 all_showers.merge(all_truth_primary_showers_ohcal());
0205
0206 return all_showers;
0207 }
0208
0209 std::set<PHG4Shower*> ParticleFlowEval::all_truth_primary_showers_ihcal()
0210 {
0211 return _clustereval_ihcal.all_truth_primary_showers(_topo_ihcal);
0212 }
0213
0214 std::set<PHG4Shower*> ParticleFlowEval::all_truth_primary_showers_ohcal()
0215 {
0216 return _clustereval_ohcal.all_truth_primary_showers(_topo_ohcal);
0217 }
0218
0219 PHG4Shower* ParticleFlowEval::max_truth_primary_shower_by_energy_emcal()
0220 {
0221 if(_topo_emcal.size() == 0) return nullptr;
0222
0223 std::map<int, std::pair<float, PHG4Shower*>> shower_energy;
0224
0225 for(unsigned int i = 0; i < _topo_emcal.size(); i++)
0226 {
0227 std::set<PHG4Shower*> showers = _clustereval_emcal.all_truth_primary_showers(_topo_emcal.at(i));
0228
0229 for(std::set<PHG4Shower*>::iterator itr = showers.begin(); itr != showers.end(); itr++)
0230 {
0231 if(shower_energy.find((*itr)->get_id()) != shower_energy.end())
0232 {
0233 shower_energy[(*itr)->get_id()].first += _clustereval_emcal.get_energy_contribution(_topo_emcal.at(i), (*itr));
0234 }
0235 else
0236 {
0237 shower_energy[(*itr)->get_id()].first = _clustereval_emcal.get_energy_contribution(_topo_emcal.at(i), (*itr));
0238 shower_energy[(*itr)->get_id()].second = (*itr);
0239 }
0240 }
0241 }
0242
0243 int max_index = -1;
0244 int max_energy = -1;
0245
0246 for(std::map<int, std::pair<float, PHG4Shower*>>::iterator itr = shower_energy.begin(); itr != shower_energy.end(); ++itr)
0247 {
0248 if(itr->second.first > max_energy)
0249 {
0250 max_energy = itr->second.first;
0251 max_index = itr->first;
0252 }
0253 }
0254
0255 return shower_energy[max_index].second;
0256 }
0257
0258 PHG4Shower* ParticleFlowEval::max_truth_primary_shower_by_energy_emcal(RawCluster *cluster)
0259 {
0260 return _clustereval_emcal.max_truth_primary_shower_by_energy(cluster);
0261 }
0262
0263 PHG4Shower* ParticleFlowEval::max_truth_primary_shower_by_energy_hcal()
0264 {
0265 PHG4Shower *max_shower = nullptr;
0266
0267 std::set<PHG4Shower*> showers = all_truth_primary_showers_hcal();
0268
0269 float max_shower_energy = 0.;
0270
0271 for(std::set<PHG4Shower*>::iterator itr = showers.begin(); itr != showers.end(); itr++)
0272 {
0273 float shower_energy = _clustereval_ihcal.get_energy_contribution(_topo_ihcal, (*itr));
0274 shower_energy += _clustereval_ohcal.get_energy_contribution(_topo_ohcal, (*itr));
0275
0276 if(shower_energy > max_shower_energy)
0277 {
0278 max_shower_energy = shower_energy;
0279 max_shower = (*itr);
0280 }
0281 }
0282
0283 return max_shower;
0284 }
0285
0286 PHG4Shower* ParticleFlowEval::max_truth_primary_shower_by_energy_ihcal()
0287 {
0288 return _clustereval_ihcal.max_truth_primary_shower_by_energy(_topo_ihcal);
0289 }
0290
0291 PHG4Shower* ParticleFlowEval::max_truth_primary_shower_by_energy_ohcal()
0292 {
0293 return _clustereval_ihcal.max_truth_primary_shower_by_energy(_topo_ohcal);
0294 }
0295
0296 std::set<PHG4Particle*> ParticleFlowEval::all_truth_primary_particles_emcal()
0297 {
0298 std::set<PHG4Particle*> all_primaries;
0299
0300 for(unsigned int i = 0; i < _topo_emcal.size(); i++)
0301 {
0302 std::set<PHG4Particle*> primaries = all_truth_primary_particles_emcal(_topo_emcal.at(i));
0303
0304 for(std::set<PHG4Particle*>::iterator itr = primaries.begin(); itr != primaries.end(); itr++)
0305 {
0306 if(all_primaries.find((*itr)) != all_primaries.end())
0307 {
0308 continue;
0309 }
0310 else
0311 {
0312 all_primaries.insert(*itr);
0313 }
0314 }
0315 }
0316
0317 return all_primaries;
0318 }
0319
0320 std::set<PHG4Particle*> ParticleFlowEval::all_truth_primary_particles_emcal(RawCluster *cluster)
0321 {
0322 return _clustereval_emcal.all_truth_primary_particles(cluster);
0323 }
0324
0325 std::set<PHG4Particle*> ParticleFlowEval::all_truth_primary_particles_hcal()
0326 {
0327 std::set<PHG4Particle*> all_primaries = all_truth_primary_particles_ihcal();
0328
0329 all_primaries.merge(all_truth_primary_particles_ohcal());
0330
0331 return all_primaries;
0332 }
0333
0334 std::set<PHG4Particle*> ParticleFlowEval::all_truth_primary_particles_ihcal()
0335 {
0336 return _clustereval_ihcal.all_truth_primary_particles(_topo_ihcal);
0337 }
0338
0339 std::set<PHG4Particle*> ParticleFlowEval::all_truth_primary_particles_ohcal()
0340 {
0341 return _clustereval_ohcal.all_truth_primary_particles(_topo_ohcal);
0342 }
0343
0344 PHG4Particle* ParticleFlowEval::max_truth_primary_particle_by_energy_emcal()
0345 {
0346 if(_topo_emcal.size() == 0)
0347 {
0348 return nullptr;
0349 }
0350
0351 std::map<int, std::pair<float, PHG4Particle*>> primary_energy;
0352
0353 for(unsigned int i = 0; i < _topo_emcal.size(); i++)
0354 {
0355 std::set<PHG4Particle*> primaries = _clustereval_emcal.all_truth_primary_particles(_topo_emcal.at(i));
0356
0357 for(std::set<PHG4Particle*>::iterator itr = primaries.begin(); itr != primaries.end(); itr++)
0358 {
0359 if(primary_energy.find((*itr)->get_track_id()) != primary_energy.end())
0360 {
0361 primary_energy[(*itr)->get_track_id()].first += _clustereval_emcal.get_energy_contribution(_topo_emcal.at(i), (*itr));
0362 }
0363 else
0364 {
0365 primary_energy[(*itr)->get_track_id()].first = _clustereval_emcal.get_energy_contribution(_topo_emcal.at(i), (*itr));
0366 primary_energy[(*itr)->get_track_id()].second = (*itr);
0367 }
0368 }
0369 }
0370
0371 int max_index = -1;
0372 int max_energy = -1;
0373
0374 for(std::map<int, std::pair<float, PHG4Particle*>>::iterator itr = primary_energy.begin(); itr != primary_energy.end(); ++itr)
0375 {
0376 if(itr->second.first > max_energy)
0377 {
0378 max_energy = itr->second.first;
0379 max_index = itr->first;
0380 }
0381 }
0382
0383 return primary_energy[max_index].second;
0384 }
0385
0386 PHG4Particle* ParticleFlowEval::max_truth_primary_particle_by_energy_emcal(RawCluster* cluster)
0387 {
0388 return _clustereval_emcal.max_truth_primary_particle_by_energy(cluster);
0389 }
0390
0391 PHG4Particle* ParticleFlowEval::max_truth_primary_particle_by_energy_hcal()
0392 {
0393 PHG4Particle *max_primary = nullptr;
0394
0395 std::set<PHG4Particle*> primaries = all_truth_primary_particles_hcal();
0396
0397 float max_primary_energy = 0.;
0398
0399 for(std::set<PHG4Particle*>::iterator itr = primaries.begin(); itr != primaries.end(); itr++)
0400 {
0401 float primary_energy = _clustereval_ihcal.get_energy_contribution(_topo_ihcal, (*itr));
0402 primary_energy += _clustereval_ohcal.get_energy_contribution(_topo_ohcal, (*itr));
0403
0404 if(primary_energy > max_primary_energy)
0405 {
0406 max_primary_energy = primary_energy;
0407 max_primary = (*itr);
0408 }
0409 }
0410
0411 return max_primary;
0412 }
0413
0414 PHG4Particle* ParticleFlowEval::max_truth_primary_particle_by_energy_ihcal()
0415 {
0416 return _clustereval_ihcal.max_truth_primary_particle_by_energy(_topo_ihcal);
0417 }
0418
0419 PHG4Particle* ParticleFlowEval::max_truth_primary_particle_by_energy_ohcal()
0420 {
0421 return _clustereval_ohcal.max_truth_primary_particle_by_energy(_topo_ohcal);
0422 }
0423
0424 PHG4Particle* ParticleFlowEval::max_truth_particle_by_nclusters_track()
0425 {
0426 return _trackeval.max_truth_particle_by_nclusters(_track);
0427 }
