Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-12-16 09:20:18 

0001 #include "ParticleFlowReco.h"
0002 
0003 #include "ParticleFlowElementContainer.h"
0004 #include "ParticleFlowElementv1.h"
0005 
0006 #include <globalvertex/GlobalVertex.h>
0007 #include <globalvertex/GlobalVertexMap.h>
0008 
0009 #include <calobase/RawCluster.h>
0010 #include <calobase/RawClusterContainer.h>
0011 #include <calobase/RawClusterUtility.h>
0012 #include <calobase/RawTowerGeom.h>
0013 #include <calobase/RawTowerGeomContainer.h>
0014 
0015 #include <trackbase_historic/SvtxTrack.h>
0016 #include <trackbase_historic/SvtxTrackMap.h>
0017 #include <trackbase_historic/SvtxTrackState.h>
0018 
0019 #include <fun4all/Fun4AllReturnCodes.h>
0020 
0021 #include <phool/PHCompositeNode.h>
0022 #include <phool/PHRandomSeed.h>
0023 #include <phool/getClass.h>
0024 
0025 #include <TLorentzVector.h>
0026 
0027 #include <gsl/gsl_randist.h>
0028 #include <gsl/gsl_rng.h>  // for gsl_rng_uniform_pos
0029 
0030 #include <cmath>
0031 #include <iostream>
0032 
0033 // examine second value of std::pair, sort by smallest
0034 bool sort_by_pair_second_lowest(const std::pair<int, float> &a, const std::pair<int, float> &b)
0035 {
0036   return (a.second < b.second);
0037 }
0038 
0039 float ParticleFlowReco::calculate_dR(float eta1, float eta2, float phi1, float phi2)
0040 {
0041   float deta = eta1 - eta2;
0042   float dphi = phi1 - phi2;
0043   while (dphi > M_PI)
0044   {
0045     dphi -= 2 * M_PI;
0046   }
0047   while (dphi < -M_PI)
0048   {
0049     dphi += 2 * M_PI;
0050   }
0051   return sqrt(pow(deta, 2) + pow(dphi, 2));
0052 }
0053 
0054 std::pair<float, float> ParticleFlowReco::get_expected_signature(int trk)
0055 {
0056   float response = (0.553437 + 0.0572246 * log(_pflow_TRK_p[trk])) * _pflow_TRK_p[trk];
0057   float resolution = sqrt(pow(0.119123, 2) + (pow(0.312361, 2) / _pflow_TRK_p[trk])) * _pflow_TRK_p[trk];
0058 
0059   std::pair<float, float> expected_signature(response, resolution);
0060 
0061   return expected_signature;
0062 }
0063 
0064 //____________________________________________________________________________..
0065 ParticleFlowReco::ParticleFlowReco(const std::string &name)
0066   : SubsysReco(name)
0067 {
0068 }
0069 
0070 //____________________________________________________________________________..
0071 int ParticleFlowReco::InitRun(PHCompositeNode *topNode)
0072 {
0073   return CreateNode(topNode);
0074 }
0075 
0076 //____________________________________________________________________________..
0077 int ParticleFlowReco::process_event(PHCompositeNode *topNode)
0078 {
0079   if (Verbosity() > 0)
0080   {
0081     std::cout << "ParticleFlowReco::process_event with Nsigma = " << _energy_match_Nsigma << std::endl;
0082   }
0083   // get handle to pflow node
0084   ParticleFlowElementContainer *pflowContainer = findNode::getClass<ParticleFlowElementContainer>(topNode, "ParticleFlowElements");
0085   if (!pflowContainer)
0086   {
0087     std::cout << " ERROR -- can't find ParticleFlowElements node after it should have been created" << std::endl;
0088     return Fun4AllReturnCodes::ABORTEVENT;
0089   }
0090 
0091   // used for indexing PFlow elements in container
0092   int global_pflow_index = 0;
0093 
0094   // read in tower geometries
0095   RawTowerGeomContainer *geomEM = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC");
0096   RawTowerGeomContainer *geomIH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0097   RawTowerGeomContainer *geomOH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0098 
0099   if (!geomEM || !geomIH || !geomOH)
0100   {
0101     std::cout << "ParticleFlowReco::process_event : FATAL ERROR, cannot find tower geometry containers" << std::endl;
0102     return Fun4AllReturnCodes::ABORTEVENT;
0103   }
0104 
0105   // read in clusters
0106   RawClusterContainer *clustersEM = findNode::getClass<RawClusterContainer>(topNode, "TOPOCLUSTER_EMCAL");
0107   RawClusterContainer *clustersHAD = findNode::getClass<RawClusterContainer>(topNode, "TOPOCLUSTER_HCAL");
0108 
0109   if (!clustersEM)
0110   {
0111     std::cout << "ParticleFlowReco::process_event : FATAL ERROR, cannot find cluster container TOPOCLUSTER_EMCAL" << std::endl;
0112     return Fun4AllReturnCodes::ABORTEVENT;
0113   }
0114   if (!clustersHAD)
0115   {
0116     std::cout << "ParticleFlowReco::process_event : FATAL ERROR, cannot find cluster container TOPOCLUSTER_HCAL" << std::endl;
0117     return Fun4AllReturnCodes::ABORTEVENT;
0118   }
0119 
0120   // reset internal particle-flow representation
0121   _pflow_TRK_p.clear();
0122   _pflow_TRK_eta.clear();
0123   _pflow_TRK_phi.clear();
0124   _pflow_TRK_match_EM.clear();
0125   _pflow_TRK_match_HAD.clear();
0126   _pflow_TRK_addtl_match_EM.clear();
0127   _pflow_TRK_trk.clear();
0128   _pflow_TRK_EMproj_phi.clear();
0129   _pflow_TRK_EMproj_eta.clear();
0130   _pflow_TRK_HADproj_eta.clear();
0131   _pflow_TRK_HADproj_phi.clear();
0132 
0133   _pflow_EM_E.clear();
0134   _pflow_EM_eta.clear();
0135   _pflow_EM_phi.clear();
0136   _pflow_EM_tower_eta.clear();
0137   _pflow_EM_tower_phi.clear();
0138   _pflow_EM_match_HAD.clear();
0139   _pflow_EM_match_TRK.clear();
0140   _pflow_EM_cluster.clear();
0141 
0142   _pflow_HAD_E.clear();
0143   _pflow_HAD_eta.clear();
0144   _pflow_HAD_phi.clear();
0145   _pflow_HAD_tower_eta.clear();
0146   _pflow_HAD_tower_phi.clear();
0147   _pflow_HAD_match_EM.clear();
0148   _pflow_HAD_match_TRK.clear();
0149   _pflow_HAD_cluster.clear();
0150 
0151   GlobalVertexMap *vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0152   GlobalVertex *vertex = nullptr;
0153 
0154   if (vertexmap)
0155   {
0156     if (!vertexmap->empty())
0157     {
0158       vertex = (vertexmap->begin()->second);
0159     }
0160   }
0161 
0162   if (Verbosity() > 2)
0163   {
0164     std::cout << "ParticleFlowReco::process_event : initial population of TRK, EM, HAD objects " << std::endl;
0165   }
0166 
0167   // read in tracks with > 0.5 GeV
0168   {
0169     SvtxTrackMap *trackmap = findNode::getClass<SvtxTrackMap>(topNode, _track_map_name);
0170 
0171     float cemcradius = geomEM->get_radius();
0172     float ohcalradius = geomOH->get_radius();
0173 
0174     for (auto &iter : *trackmap)
0175     {
0176       SvtxTrack *track = iter.second;
0177 
0178       if(_only_crossing_zero && (track->get_crossing() != 0))
0179       {
0180         continue;
0181       }
0182 
0183       if (track->get_pt() < 0.5)
0184       {
0185         continue;
0186       }
0187 
0188       if (std::fabs(track->get_eta()) > 1.1)
0189       {
0190         continue;
0191       }
0192 
0193       if (Verbosity() > 2)
0194       {
0195         std::cout << "Track with p= " << track->get_p() << ", eta / phi = "
0196                   << track->get_eta() << " / " << track->get_phi()
0197                   << std::endl;
0198       }
0199 
0200       _pflow_TRK_trk.push_back(track);
0201       _pflow_TRK_p.push_back(track->get_p());
0202       _pflow_TRK_eta.push_back(track->get_eta());
0203       _pflow_TRK_phi.push_back(track->get_phi());
0204       _pflow_TRK_match_EM.emplace_back();
0205       _pflow_TRK_match_HAD.emplace_back();
0206       _pflow_TRK_addtl_match_EM.emplace_back();
0207 
0208       SvtxTrackState *cemcstate = track->get_state(cemcradius);
0209       SvtxTrackState *ohstate = track->get_state(ohcalradius);
0210       /// Get the track projections. If they failed for some reason, just use the track
0211       /// phi and eta values at the point of closest approach
0212       if (cemcstate)
0213       {
0214         _pflow_TRK_EMproj_phi.push_back(std::atan2(cemcstate->get_y(), cemcstate->get_x()));
0215         _pflow_TRK_EMproj_eta.push_back(std::asinh(cemcstate->get_z() / std::sqrt((cemcstate->get_x() * cemcstate->get_x()) + (cemcstate->get_y() * cemcstate->get_y()))));
0216       }
0217       else
0218       {
0219         _pflow_TRK_EMproj_phi.push_back(track->get_phi());
0220         _pflow_TRK_EMproj_eta.push_back(track->get_eta());
0221       }
0222       if (ohstate)
0223       {
0224         _pflow_TRK_HADproj_phi.push_back(std::atan2(ohstate->get_y(), ohstate->get_x()));
0225         _pflow_TRK_HADproj_eta.push_back(std::asinh(ohstate->get_z() / std::sqrt((ohstate->get_x() * ohstate->get_x()) + (ohstate->get_y() * ohstate->get_y()))));
0226       }
0227       else
0228       {
0229         _pflow_TRK_HADproj_phi.push_back(track->get_phi());
0230         _pflow_TRK_HADproj_eta.push_back(track->get_eta());
0231       }
0232     }
0233 
0234   }  //
0235 
0236   // read in EMCal topoClusters with E > 0.2 GeV
0237   {
0238     RawClusterContainer::ConstRange begin_end = clustersEM->getClusters();
0239     for (RawClusterContainer::ConstIterator hiter = begin_end.first; hiter != begin_end.second; ++hiter)
0240     {
0241       float cluster_E = hiter->second->get_energy();
0242       if (cluster_E < 0.2)
0243       {
0244         continue;
0245       }
0246 
0247       float cluster_phi = hiter->second->get_phi();
0248       /// default assume at vx_z = 0
0249       float cluster_theta = M_PI / 2.0 - std::atan2(hiter->second->get_z(), hiter->second->get_r());
0250       float cluster_eta = -1 * log(tan(cluster_theta / 2.0));
0251 
0252       if (vertex)
0253       {
0254         cluster_eta = RawClusterUtility::GetPseudorapidity(*(hiter->second), CLHEP::Hep3Vector(vertex->get_x(), vertex->get_y(), vertex->get_z()));
0255       }
0256 
0257       _pflow_EM_E.push_back(cluster_E);
0258       _pflow_EM_eta.push_back(cluster_eta);
0259       _pflow_EM_phi.push_back(cluster_phi);
0260       _pflow_EM_cluster.push_back(hiter->second);
0261       _pflow_EM_match_HAD.emplace_back();
0262       _pflow_EM_match_TRK.emplace_back();
0263 
0264       if (Verbosity() > 5 && cluster_E > 0.2)
0265       {
0266         std::cout << " EM topoCluster with E = " << cluster_E << ", eta / phi = " << cluster_eta << " / " << cluster_phi << " , nTow = " << hiter->second->getNTowers() << std::endl;
0267       }
0268 
0269       std::vector<float> this_cluster_tower_eta;
0270       std::vector<float> this_cluster_tower_phi;
0271 
0272       // read in towers
0273       RawCluster::TowerConstRange begin_end_towers = hiter->second->get_towers();
0274       for (RawCluster::TowerConstIterator iter = begin_end_towers.first; iter != begin_end_towers.second; ++iter)
0275       {
0276         if (RawTowerDefs::decode_caloid(iter->first) == RawTowerDefs::CalorimeterId::CEMC)
0277         {
0278           RawTowerGeom *tower_geom = geomEM->get_tower_geometry(iter->first);
0279 
0280           this_cluster_tower_phi.push_back(tower_geom->get_phi());
0281           this_cluster_tower_eta.push_back(tower_geom->get_eta());
0282         }
0283         else
0284         {
0285           std::cout << "ParticleFlowReco::process_event : FATAL ERROR , EM topoClusters seem to contain HCal towers" << std::endl;
0286           return Fun4AllReturnCodes::ABORTEVENT;
0287         }
0288       }  // close tower loop
0289 
0290       _pflow_EM_tower_eta.push_back(this_cluster_tower_eta);
0291       _pflow_EM_tower_phi.push_back(this_cluster_tower_phi);
0292 
0293     }  // close cluster loop
0294 
0295   }  // close
0296 
0297   // read in HCal topoClusters with E > 0.2 GeV
0298   {
0299     RawClusterContainer::ConstRange begin_end = clustersHAD->getClusters();
0300     for (RawClusterContainer::ConstIterator hiter = begin_end.first; hiter != begin_end.second; ++hiter)
0301     {
0302       float cluster_E = hiter->second->get_energy();
0303       if (cluster_E < 0.2)
0304       {
0305         continue;
0306       }
0307 
0308       float cluster_phi = hiter->second->get_phi();
0309       // for now, assume event at vx_z = 0
0310       float cluster_theta = M_PI / 2.0 - std::atan2(hiter->second->get_z(), hiter->second->get_r());
0311       float cluster_eta = -1 * log(tan(cluster_theta / 2.0));
0312       if (vertex)
0313       {
0314         cluster_eta = RawClusterUtility::GetPseudorapidity(*(hiter->second), CLHEP::Hep3Vector(vertex->get_x(), vertex->get_y(), vertex->get_z()));
0315       }
0316 
0317       _pflow_HAD_E.push_back(cluster_E);
0318       _pflow_HAD_eta.push_back(cluster_eta);
0319       _pflow_HAD_phi.push_back(cluster_phi);
0320       _pflow_HAD_cluster.push_back(hiter->second);
0321 
0322       _pflow_HAD_match_EM.emplace_back();
0323       _pflow_HAD_match_TRK.emplace_back();
0324 
0325       if (Verbosity() > 5 && cluster_E > 0.2)
0326       {
0327         std::cout << " HAD topoCluster with E = " << cluster_E << ", eta / phi = " << cluster_eta << " / " << cluster_phi << " , nTow = " << hiter->second->getNTowers() << std::endl;
0328       }
0329 
0330       std::vector<float> this_cluster_tower_eta;
0331       std::vector<float> this_cluster_tower_phi;
0332 
0333       // read in towers
0334       RawCluster::TowerConstRange begin_end_towers = hiter->second->get_towers();
0335       for (RawCluster::TowerConstIterator iter = begin_end_towers.first; iter != begin_end_towers.second; ++iter)
0336       {
0337         if (RawTowerDefs::decode_caloid(iter->first) == RawTowerDefs::CalorimeterId::HCALIN)
0338         {
0339           RawTowerGeom *tower_geom = geomIH->get_tower_geometry(iter->first);
0340 
0341           this_cluster_tower_phi.push_back(tower_geom->get_phi());
0342           this_cluster_tower_eta.push_back(tower_geom->get_eta());
0343         }
0344 
0345         else if (RawTowerDefs::decode_caloid(iter->first) == RawTowerDefs::CalorimeterId::HCALOUT)
0346         {
0347           RawTowerGeom *tower_geom = geomOH->get_tower_geometry(iter->first);
0348 
0349           this_cluster_tower_phi.push_back(tower_geom->get_phi());
0350           this_cluster_tower_eta.push_back(tower_geom->get_eta());
0351         }
0352         else
0353         {
0354           std::cout << "ParticleFlowReco::process_event : FATAL ERROR , HCal topoClusters seem to contain EM towers" << std::endl;
0355           return Fun4AllReturnCodes::ABORTEVENT;
0356         }
0357 
0358       }  // close tower loop
0359 
0360       _pflow_HAD_tower_eta.push_back(this_cluster_tower_eta);
0361       _pflow_HAD_tower_phi.push_back(this_cluster_tower_phi);
0362 
0363     }  // close cluster loop
0364 
0365   }  // close
0366 
0367   // BEGIN LINKING STEP
0368 
0369   // Link TRK -> EM (best match, but keep reserve of others), and TRK -> HAD (best match)
0370   if (Verbosity() > 2)
0371   {
0372     std::cout << "ParticleFlowReco::process_event : TRK -> EM and TRK -> HAD linking " << std::endl;
0373   }
0374 
0375   for (unsigned int trk = 0; trk < _pflow_TRK_p.size(); trk++)
0376   {
0377     if (Verbosity() > 10)
0378     {
0379       std::cout << " TRK " << trk << " with p / eta / phi = " << _pflow_TRK_p[trk] << " / " << _pflow_TRK_eta[trk] << " / " << _pflow_TRK_phi[trk] << std::endl;
0380     }
0381 
0382     // TRK -> EM link
0383     float min_em_dR = 0.2;
0384     int min_em_index = -1;
0385 
0386     for (unsigned int em = 0; em < _pflow_EM_E.size(); em++)
0387     {
0388       float dR = calculate_dR(_pflow_TRK_EMproj_eta[trk], _pflow_EM_eta[em], _pflow_TRK_EMproj_phi[trk], _pflow_EM_phi[em]);
0389 
0390       if (dR > 0.2)
0391       {
0392         continue;
0393       }
0394 
0395       bool has_overlap = false;
0396 
0397       for (unsigned int tow = 0; tow < _pflow_EM_tower_eta.at(em).size(); tow++)
0398       {
0399         float tower_eta = _pflow_EM_tower_eta.at(em).at(tow);
0400         float tower_phi = _pflow_EM_tower_phi.at(em).at(tow);
0401 
0402         float deta = tower_eta - _pflow_TRK_EMproj_eta[trk];
0403         float dphi = tower_phi - _pflow_TRK_EMproj_phi[trk];
0404         if (dphi > M_PI)
0405         {
0406           dphi -= 2 * M_PI;
0407         }
0408         if (dphi < -M_PI)
0409         {
0410           dphi += 2 * M_PI;
0411         }
0412 
0413         if (std::fabs(deta) < 0.025 * 2.5 && std::fabs(dphi) < 0.025 * 2.5)
0414         {
0415           has_overlap = true;
0416           break;
0417         }
0418       }
0419 
0420       if (has_overlap)
0421       {
0422         if (Verbosity() > 5)
0423         {
0424           std::cout << " -> possible match to EM " << em << " with dR = " << dR << std::endl;
0425         }
0426 
0427         _pflow_TRK_addtl_match_EM.at(trk).emplace_back(em, dR);
0428       }
0429       else
0430       {
0431         if (Verbosity() > 5)
0432         {
0433           std::cout << " -> no match to EM " << em << " (even though dR = " << dR << " )" << std::endl;
0434         }
0435       }
0436     }
0437 
0438     // sort possible matches
0439 
0440     std::sort(_pflow_TRK_addtl_match_EM.at(trk).begin(), _pflow_TRK_addtl_match_EM.at(trk).end(), sort_by_pair_second_lowest);
0441     if (Verbosity() > 10)
0442     {
0443       for (auto &n : _pflow_TRK_addtl_match_EM.at(trk))
0444       {
0445         std::cout << " -> sorted list of matches, EM / dR = " << n.first << " / " << n.second << std::endl;
0446       }
0447     }
0448 
0449     if (!_pflow_TRK_addtl_match_EM.at(trk).empty())
0450     {
0451       min_em_index = _pflow_TRK_addtl_match_EM.at(trk).at(0).first;
0452       min_em_dR = _pflow_TRK_addtl_match_EM.at(trk).at(0).second;
0453       // delete best matched element
0454       _pflow_TRK_addtl_match_EM.at(trk).erase(_pflow_TRK_addtl_match_EM.at(trk).begin());
0455     }
0456 
0457     if (min_em_index > -1)
0458     {
0459       _pflow_EM_match_TRK.at(min_em_index).push_back(trk);
0460       _pflow_TRK_match_EM.at(trk).push_back(min_em_index);
0461 
0462       if (Verbosity() > 5)
0463       {
0464         std::cout << " -> matched EM " << min_em_index << " with pt / eta / phi = " << _pflow_EM_E.at(min_em_index) << " / " << _pflow_EM_eta.at(min_em_index) << " / " << _pflow_EM_phi.at(min_em_index) << ", dR = " << min_em_dR;
0465         std::cout << " ( " << _pflow_TRK_addtl_match_EM.at(trk).size() << " other possible matches ) " << std::endl;
0466       }
0467     }
0468     else
0469     {
0470       if (Verbosity() > 5)
0471       {
0472         std::cout << " -> no EM match! ( best dR = " << min_em_dR << " ) " << std::endl;
0473       }
0474     }
0475 
0476     // TRK -> HAD link
0477     float min_had_dR = 0.2;
0478     int min_had_index = -1;
0479     float max_had_pt = 0;
0480 
0481     // TODO: sequential linking should better happen here -- i.e. allow EM-matched HAD's into the possible pool
0482     for (unsigned int had = 0; had < _pflow_HAD_E.size(); had++)
0483     {
0484       float dR = calculate_dR(_pflow_TRK_HADproj_eta[trk], _pflow_HAD_eta[had], _pflow_TRK_HADproj_phi[trk], _pflow_HAD_phi[had]);
0485 
0486       if (dR > 0.5)
0487       {
0488         continue;
0489       }
0490 
0491       bool has_overlap = false;
0492 
0493       for (unsigned int tow = 0; tow < _pflow_HAD_tower_eta.at(had).size(); tow++)
0494       {
0495         float tower_eta = _pflow_HAD_tower_eta.at(had).at(tow);
0496         float tower_phi = _pflow_HAD_tower_phi.at(had).at(tow);
0497 
0498         float deta = tower_eta - _pflow_TRK_HADproj_eta[trk];
0499         float dphi = tower_phi - _pflow_TRK_HADproj_phi[trk];
0500         if (dphi > M_PI)
0501         {
0502           dphi -= 2 * M_PI;
0503         }
0504         if (dphi < -M_PI)
0505         {
0506           dphi += 2 * M_PI;
0507         }
0508 
0509         if (std::fabs(deta) < 0.1 * 1.5 && std::fabs(dphi) < 0.1 * 1.5)
0510         {
0511           has_overlap = true;
0512           break;
0513         }
0514       }
0515 
0516       if (has_overlap)
0517       {
0518         if (Verbosity() > 5)
0519         {
0520           std::cout << " -> possible match to HAD " << had << " with dR = " << dR << std::endl;
0521         }
0522 
0523         if (_pflow_HAD_E.at(had) > max_had_pt)
0524         {
0525           max_had_pt = _pflow_HAD_E.at(had);
0526           min_had_index = had;
0527           min_had_dR = dR;
0528         }
0529       }
0530       else
0531       {
0532         if (Verbosity() > 5)
0533         {
0534           std::cout << " -> no match to HAD " << had << " (even though dR = " << dR << " )" << std::endl;
0535         }
0536       }
0537     }
0538 
0539     if (min_had_index > -1)
0540     {
0541       _pflow_HAD_match_TRK.at(min_had_index).push_back(trk);
0542       _pflow_TRK_match_HAD.at(trk).push_back(min_had_index);
0543 
0544       if (Verbosity() > 5)
0545       {
0546         std::cout << " -> matched HAD " << min_had_index << " with pt / eta / phi = " << _pflow_HAD_E.at(min_had_index) << " / " << _pflow_HAD_eta.at(min_had_index) << " / " << _pflow_HAD_phi.at(min_had_index) << ", dR = " << min_had_dR << std::endl;
0547       }
0548     }
0549     else
0550     {
0551       if (Verbosity() > 5)
0552       {
0553         std::cout << " -> no HAD match! ( best dR = " << min_had_dR << " ) " << std::endl;
0554       }
0555     }
0556   }
0557 
0558   // EM->HAD linking
0559   if (Verbosity() > 2)
0560   {
0561     std::cout << "ParticleFlowReco::process_event : EM -> HAD linking " << std::endl;
0562   }
0563 
0564   for (unsigned int em = 0; em < _pflow_EM_E.size(); em++)
0565   {
0566     if (Verbosity() > 10)
0567     {
0568       std::cout << " EM with E / eta / phi = " << _pflow_EM_E[em] << " / " << _pflow_EM_eta[em] << " / " << _pflow_EM_phi[em] << std::endl;
0569     }
0570 
0571     // TRK -> HAD link
0572     float min_had_dR = 0.2;
0573     int min_had_index = -1;
0574     float max_had_pt = 0;
0575 
0576     for (unsigned int had = 0; had < _pflow_HAD_E.size(); had++)
0577     {
0578       float dR = calculate_dR(_pflow_EM_eta[em], _pflow_HAD_eta[had], _pflow_EM_phi[em], _pflow_HAD_phi[had]);
0579       if (dR > 0.5)
0580       {
0581         continue;
0582       }
0583 
0584       bool has_overlap = false;
0585 
0586       for (unsigned int tow = 0; tow < _pflow_HAD_tower_eta.at(had).size(); tow++)
0587       {
0588         float tower_eta = _pflow_HAD_tower_eta.at(had).at(tow);
0589         float tower_phi = _pflow_HAD_tower_phi.at(had).at(tow);
0590 
0591         float deta = tower_eta - _pflow_EM_eta[em];
0592         float dphi = tower_phi - _pflow_EM_phi[em];
0593         if (dphi > M_PI)
0594         {
0595           dphi -= 2 * M_PI;
0596         }
0597         if (dphi < -M_PI)
0598         {
0599           dphi += 2 * M_PI;
0600         }
0601 
0602         if (std::fabs(deta) < 0.1 * 1.5 && std::fabs(dphi) < 0.1 * 1.5)
0603         {
0604           has_overlap = true;
0605           break;
0606         }
0607       }
0608 
0609       if (has_overlap)
0610       {
0611         if (Verbosity() > 5)
0612         {
0613           std::cout << " -> possible match to HAD " << had << " with dR = " << dR << std::endl;
0614         }
0615 
0616         if (_pflow_HAD_E.at(had) > max_had_pt)
0617         {
0618           max_had_pt = _pflow_HAD_E.at(had);
0619           min_had_index = had;
0620           min_had_dR = dR;
0621         }
0622       }
0623       else
0624       {
0625         if (Verbosity() > 5)
0626         {
0627           std::cout << " -> no match to HAD " << had << " (even though dR = " << dR << " )" << std::endl;
0628         }
0629       }
0630     }
0631 
0632     if (min_had_index > -1)
0633     {
0634       _pflow_HAD_match_EM.at(min_had_index).push_back(em);
0635       _pflow_EM_match_HAD.at(em).push_back(min_had_index);
0636 
0637       if (Verbosity() > 5)
0638       {
0639         std::cout << " -> matched HAD with E / eta / phi = " << _pflow_HAD_E.at(min_had_index) << " / " << _pflow_HAD_eta.at(min_had_index) << " / " << _pflow_HAD_phi.at(min_had_index) << ", dR = " << min_had_dR << std::endl;
0640       }
0641     }
0642     else
0643     {
0644       if (Verbosity() > 5)
0645       {
0646         std::cout << " -> no HAD match! ( best dR = " << min_had_dR << " ) " << std::endl;
0647       }
0648     }
0649   }
0650 
0651   // SEQUENTIAL MATCHING: if TRK -> EM and EM -> HAD, ensure that TRK -> HAD
0652   if (Verbosity() > 2)
0653   {
0654     std::cout << "ParticleFlowReco::process_event : sequential TRK -> EM && EM -> HAD ==> TRK -> HAD matching " << std::endl;
0655   }
0656 
0657   for (unsigned int trk = 0; trk < _pflow_TRK_p.size(); trk++)
0658   {
0659     // go through all matched EMs
0660     for (unsigned int i = 0; i < _pflow_TRK_match_EM.at(trk).size(); i++)
0661     {
0662       int em = _pflow_TRK_match_EM.at(trk).at(i);
0663 
0664       // if this EM has a matched HAD...
0665       for (unsigned int j = 0; j < _pflow_EM_match_HAD.at(em).size(); j++)
0666       {
0667         int had = _pflow_EM_match_HAD.at(em).at(j);
0668 
0669         // and the TRK is NOT matched to this HAD...
0670         bool is_trk_matched_to_HAD = false;
0671         for (int existing_had : _pflow_TRK_match_HAD.at(trk))
0672         {
0673           if (had == existing_had)
0674           {
0675             is_trk_matched_to_HAD = true;
0676           }
0677         }
0678 
0679         // if this is the case, create TRK->HAD link
0680         if (!is_trk_matched_to_HAD)
0681         {
0682           _pflow_TRK_match_HAD.at(trk).push_back(had);
0683           _pflow_HAD_match_TRK.at(had).push_back(trk);
0684 
0685           if (Verbosity() > 5)
0686           {
0687             std::cout << " TRK " << trk << " with pt / eta / phi = " << _pflow_TRK_p.at(trk) << " / " << _pflow_TRK_eta.at(trk) << " / " << _pflow_TRK_phi.at(trk) << std::endl;
0688             std::cout << " -> sequential match to HAD " << had << " through EM " << j << std::endl;
0689           }
0690         }
0691 
0692       }  // close the HAD loop
0693 
0694     }  // close the EM loop
0695 
0696   }  // close the TRK loop
0697 
0698   // TRK->EM->HAD removal
0699   if (Verbosity() > 2)
0700   {
0701     std::cout << "ParticleFlowReco::process_event : resolve TRK(s) + EM(s) -> HAD systems " << std::endl;
0702   }
0703 
0704   for (unsigned int had = 0; had < _pflow_HAD_E.size(); had++)
0705   {
0706     // only consider HAD with matched tracks ... others we will deal with later
0707     if (_pflow_HAD_match_TRK.at(had).empty())
0708     {
0709       continue;
0710     }
0711 
0712     if (Verbosity() > 5)
0713     {
0714       std::cout << " HAD " << had << " with E / eta / phi = " << _pflow_HAD_E.at(had) << " / " << _pflow_HAD_eta.at(had) << " / " << _pflow_HAD_phi.at(had) << std::endl;
0715     }
0716 
0717     // setup for Sum-pT^trk -> calo prediction
0718     float total_TRK_p = 0;
0719     float total_expected_E = 0;
0720     float total_expected_E_var = 0;
0721 
0722     // begin with this HAD calo energy
0723     float total_EMHAD_E = _pflow_HAD_E.at(had);
0724 
0725     std::vector<RawCluster *> matchedEClusters;
0726 
0727     // iterate over the EMs matched to this HAD
0728     for (int em : _pflow_HAD_match_EM.at(had))
0729     {
0730       // ensure there is at least one track matched to this EM
0731       if (_pflow_EM_match_TRK.at(em).empty())
0732       {
0733         continue;
0734       }
0735 
0736       // add it to the total calo E
0737       total_EMHAD_E += _pflow_EM_E.at(em);
0738       matchedEClusters.push_back(_pflow_EM_cluster.at(em));
0739       if (Verbosity() > 5)
0740       {
0741         std::cout << " -> -> LINKED EM " << em << " with E / eta / phi = " << _pflow_EM_E.at(em) << " / " << _pflow_EM_eta.at(em) << " / " << _pflow_EM_phi.at(em) << std::endl;
0742       }
0743     }
0744 
0745     // iterate over the TRKs matched to this HAD
0746     for (unsigned int j = 0; j < _pflow_HAD_match_TRK.at(had).size(); j++)
0747     {
0748       int trk = _pflow_HAD_match_TRK.at(had).at(j);
0749 
0750       if (Verbosity() > 5)
0751       {
0752         std::cout << " -> -> LINKED TRK " << trk << " with p / eta / phi = " << _pflow_TRK_p.at(trk) << " / " << _pflow_TRK_eta.at(trk) << " / " << _pflow_TRK_phi.at(trk) << std::endl;
0753       }
0754 
0755       total_TRK_p += _pflow_TRK_p.at(trk);
0756 
0757       std::pair<float, float> expected_signature = get_expected_signature(trk);
0758 
0759       float expected_E_mean = expected_signature.first;
0760       float expected_E_sigma = expected_signature.second;
0761 
0762       if (Verbosity() > 5)
0763       {
0764         std::cout << " -> -> -> expected calo signature is " << expected_E_mean << " +/- " << expected_E_sigma << std::endl;
0765       }
0766 
0767       total_expected_E += expected_E_mean;
0768       total_expected_E_var += pow(expected_E_sigma, 2);
0769 
0770       // add PFlow element for each track
0771       ParticleFlowElement *pflow = new ParticleFlowElementv1();
0772 
0773       // assume pion mass
0774       TLorentzVector tlv;
0775       tlv.SetPtEtaPhiM(_pflow_TRK_p[trk] / cosh(_pflow_TRK_eta[trk]), _pflow_TRK_eta[trk], _pflow_TRK_phi[trk], 0.135);
0776 
0777       pflow->set_px(tlv.Px());
0778       pflow->set_py(tlv.Py());
0779       pflow->set_pz(tlv.Pz());
0780       pflow->set_e(tlv.E());
0781       pflow->set_track(_pflow_TRK_trk[trk]);
0782       pflow->set_eclusters(matchedEClusters);
0783       pflow->set_hcluster(_pflow_HAD_cluster.at(had));
0784       pflow->set_id(global_pflow_index);
0785       pflow->set_type(ParticleFlowElement::PFLOWTYPE::MATCHED_CHARGED_HADRON);
0786 
0787       pflowContainer->AddParticleFlowElement(global_pflow_index, pflow);
0788       global_pflow_index++;
0789     }
0790     // Track + E+HCal PF elements are created
0791 
0792     // process compatibility of fit
0793     float total_expected_E_err = std::sqrt(total_expected_E_var);
0794 
0795     if (Verbosity() > 5)
0796     {
0797       std::cout << " -> Total track Sum p = " << total_TRK_p << " , expected calo Sum E = " << total_expected_E << " +/- " << total_expected_E_err << " , observed EM+HAD Sum E = " << total_EMHAD_E << std::endl;
0798     }
0799 
0800     // if Sum pT > calo, add in additional possible matched EMs associated with tracks until that is no longer the case
0801 
0802     if (total_expected_E > total_EMHAD_E)
0803     {
0804       if (Verbosity() > 5)
0805       {
0806         std::cout << " -> Expected E > Observed E, looking for additional potential TRK->EM matches" << std::endl;
0807       }
0808 
0809       std::map<int, float> additional_EMs;
0810 
0811       for (int trk : _pflow_HAD_match_TRK.at(had))
0812       {
0813         int addtl_matches = _pflow_TRK_addtl_match_EM.at(trk).size();
0814 
0815         if (Verbosity() > 10)
0816         {
0817           std::cout << " -> -> TRK " << trk << " has " << addtl_matches << " additional matches! " << std::endl;
0818         }
0819 
0820         for (auto &n : _pflow_TRK_addtl_match_EM.at(trk))
0821         {
0822           if (Verbosity() > 10)
0823           {
0824             std::cout << " -> -> -> additional match to EM = " << n.first << " with dR = " << n.second << std::endl;
0825           }
0826 
0827           float existing_dR = 0.21;
0828           int counts = additional_EMs.count(n.first);
0829           if (counts > 0)
0830           {
0831             existing_dR = additional_EMs[n.first];
0832           }
0833           if (n.second < existing_dR)
0834           {
0835             additional_EMs[n.first] = n.second;
0836           }
0837         }
0838       }
0839 
0840       // map now assured to have only minimal dR values for each possible additional EM
0841       // translate the map to a vector of pairs, then sort by smallest dR
0842 
0843       std::vector<std::pair<int, float> > additional_EMs_vec;
0844 
0845       additional_EMs_vec.reserve(additional_EMs.size());
0846       for (auto &x : additional_EMs)
0847       {
0848         additional_EMs_vec.emplace_back(x.first, x.second);
0849       }
0850 
0851       std::sort(additional_EMs_vec.begin(), additional_EMs_vec.end(), sort_by_pair_second_lowest);
0852 
0853       if (Verbosity() > 5)
0854       {
0855         std::cout << " -> Sorting the set of potential additional EMs " << std::endl;
0856       }
0857 
0858       // now add in additional EMs until there are none left or it is no longer the case that Sum pT > calo
0859 
0860       int n_EM_added = 0;
0861       while (!additional_EMs_vec.empty() && total_expected_E > total_EMHAD_E)
0862       {
0863         int new_EM = additional_EMs_vec.at(0).first;
0864 
0865         if (Verbosity() > 5)
0866         {
0867           std::cout << " -> adding EM " << new_EM << " ( dR = " << additional_EMs_vec.at(0).second << " to the system (should not see it as orphan below)" << std::endl;
0868         }
0869 
0870         // for now, just make the first HAD-linked track point to this new EM, and vice versa
0871         _pflow_EM_match_TRK.at(new_EM).push_back(_pflow_HAD_match_TRK.at(had).at(0));
0872         _pflow_TRK_match_EM.at(_pflow_HAD_match_TRK.at(had).at(0)).push_back(new_EM);
0873 
0874         // add to expected calo
0875         total_EMHAD_E += _pflow_EM_E.at(new_EM);
0876 
0877         // erase lowest-dR EM
0878         additional_EMs_vec.erase(additional_EMs_vec.begin());
0879 
0880         n_EM_added++;
0881       }
0882 
0883       if (Verbosity() > 5)
0884       {
0885         if (n_EM_added > 0)
0886         {
0887           std::cout << "After adding N = " << n_EM_added << " any additional EMs : " << std::endl;
0888           std::cout << "-> Total track Sum p = " << total_TRK_p << " , expected calo Sum E = " << total_expected_E << " +/- " << total_expected_E_err << " , observed EM+HAD Sum E = " << total_EMHAD_E << std::endl;
0889         }
0890         else
0891         {
0892           std::cout << "No additional EMs found, continuing hypothesis check" << std::endl;
0893         }
0894       }
0895     }
0896 
0897     if (total_expected_E + _energy_match_Nsigma * total_expected_E_err > total_EMHAD_E)
0898     {
0899       if (Verbosity() > 5)
0900       {
0901         std::cout << " -> -> calo compatible within Nsigma = " << _energy_match_Nsigma << " , remove and keep tracks " << std::endl;
0902       }
0903 
0904       // PFlow elements already created from tracks above, no more needs to be done
0905     }
0906     else
0907     {
0908       float residual_energy = total_EMHAD_E - total_expected_E;
0909 
0910       if (Verbosity() > 5)
0911       {
0912         std::cout << " -> -> calo not compatible, create leftover cluster with " << residual_energy << std::endl;
0913       }
0914 
0915       // create additional PFlow element (tracks already created above)
0916       ParticleFlowElement *pflow = new ParticleFlowElementv1();
0917 
0918       // assume no mass, but could update to use K0L mass(?)
0919       TLorentzVector tlv;
0920       tlv.SetPtEtaPhiM(residual_energy / cosh(_pflow_HAD_eta[had]), _pflow_HAD_eta[had], _pflow_HAD_phi[had], 0);
0921 
0922       pflow->set_px(tlv.Px());
0923       pflow->set_py(tlv.Py());
0924       pflow->set_pz(tlv.Pz());
0925       pflow->set_e(tlv.E());
0926       pflow->set_track(nullptr);
0927       pflow->set_eclusters(matchedEClusters);
0928       pflow->set_hcluster(_pflow_HAD_cluster.at(had));
0929       pflow->set_id(global_pflow_index);
0930       pflow->set_type(ParticleFlowElement::PFLOWTYPE::LEFTOVER_EM_PARTICLE);
0931 
0932       pflowContainer->AddParticleFlowElement(global_pflow_index, pflow);
0933       global_pflow_index++;
0934     }
0935 
0936   }  // close HAD loop
0937 
0938   // TRK->EM removal
0939 
0940   if (Verbosity() > 2)
0941   {
0942     std::cout << "ParticleFlowReco::process_event : resolve TRK(s) -> EM(s) ( + no HAD) systems " << std::endl;
0943   }
0944 
0945   for (unsigned int em = 0; em < _pflow_EM_E.size(); em++)
0946   {
0947     // only consider EM with matched tracks, but no matched HADs
0948     if (!_pflow_EM_match_HAD.at(em).empty())
0949     {
0950       continue;
0951     }
0952     if (_pflow_EM_match_TRK.at(em).empty())
0953     {
0954       continue;
0955     }
0956 
0957     if (Verbosity() > 5)
0958     {
0959       std::cout << " EM " << em << " with E / eta / phi = " << _pflow_EM_E.at(em) << " / " << _pflow_EM_eta.at(em) << " / " << _pflow_EM_phi.at(em) << std::endl;
0960     }
0961 
0962     // setup for Sum-pT^trk -> calo prediction
0963     float total_TRK_p = 0;
0964     float total_expected_E = 0;
0965     float total_expected_E_var = 0;
0966 
0967     // begin with this EM calo energy
0968     float total_EM_E = _pflow_EM_E.at(em);
0969 
0970     // iterate over the TRKs matched to this EM
0971     for (unsigned int j = 0; j < _pflow_EM_match_TRK.at(em).size(); j++)
0972     {
0973       int trk = _pflow_EM_match_TRK.at(em).at(j);
0974 
0975       if (Verbosity() > 5)
0976       {
0977         std::cout << " -> -> LINKED TRK with p / eta / phi = " << _pflow_TRK_p.at(trk) << " / " << _pflow_TRK_eta.at(trk) << " / " << _pflow_TRK_phi.at(trk) << std::endl;
0978       }
0979 
0980       total_TRK_p += _pflow_TRK_p.at(trk);
0981 
0982       std::pair<float, float> expected_signature = get_expected_signature(trk);
0983 
0984       float expected_E_mean = expected_signature.first;
0985       float expected_E_sigma = expected_signature.second;
0986 
0987       if (Verbosity() > 5)
0988       {
0989         std::cout << " -> -> -> expected calo signature is " << expected_E_mean << " +/- " << expected_E_sigma << std::endl;
0990       }
0991 
0992       total_expected_E += expected_E_mean;
0993       total_expected_E_var += pow(expected_E_sigma, 2);
0994 
0995       // add PFlow element for each track
0996       ParticleFlowElement *pflow = new ParticleFlowElementv1();
0997 
0998       // assume pion mass
0999       TLorentzVector tlv;
1000       tlv.SetPtEtaPhiM(_pflow_TRK_p[trk] / cosh(_pflow_TRK_eta[trk]), _pflow_TRK_eta[trk], _pflow_TRK_phi[trk], 0.135);
1001 
1002       std::vector<RawCluster *> eclus;
1003       eclus.push_back(_pflow_EM_cluster.at(em));
1004 
1005       pflow->set_px(tlv.Px());
1006       pflow->set_py(tlv.Py());
1007       pflow->set_pz(tlv.Pz());
1008       pflow->set_e(tlv.E());
1009       pflow->set_track(_pflow_TRK_trk.at(trk));
1010       pflow->set_eclusters(eclus);
1011       pflow->set_hcluster(nullptr);
1012       pflow->set_id(global_pflow_index);
1013       pflow->set_type(ParticleFlowElement::PFLOWTYPE::MATCHED_CHARGED_HADRON);
1014 
1015       pflowContainer->AddParticleFlowElement(global_pflow_index, pflow);
1016       global_pflow_index++;
1017     }
1018 
1019     // process compatibility of fit
1020     float total_expected_E_err = std::sqrt(total_expected_E_var);
1021 
1022     if (Verbosity() > 5)
1023     {
1024       std::cout << " -> Total track Sum p = " << total_TRK_p << " , expected calo Sum E = " << total_expected_E << " +/- " << total_expected_E_err << " , observed EM Sum E = " << total_EM_E << std::endl;
1025     }
1026 
1027     if (total_expected_E + _energy_match_Nsigma * total_expected_E_err > total_EM_E)
1028     {
1029       if (Verbosity() > 5)
1030       {
1031         std::cout << " -> -> calo compatible within Nsigma = " << _energy_match_Nsigma << "  , remove and keep tracks " << std::endl;
1032       }
1033 
1034       // PFlow elements already created from tracks above, no more needs to be done
1035     }
1036     else
1037     {
1038       float residual_energy = total_EM_E - total_expected_E;
1039 
1040       if (Verbosity() > 5)
1041       {
1042         std::cout << " -> -> calo not compatible, create leftover cluster with " << residual_energy << std::endl;
1043       }
1044 
1045       // create additional PFlow element (tracks already created above)
1046       ParticleFlowElement *pflow = new ParticleFlowElementv1();
1047 
1048       // assume no mass, but could update to use K0L mass(?)
1049       TLorentzVector tlv;
1050       tlv.SetPtEtaPhiM(residual_energy / cosh(_pflow_EM_eta[em]), _pflow_EM_eta[em], _pflow_EM_phi[em], 0);
1051 
1052       std::vector<RawCluster *> eclus;
1053       eclus.push_back(_pflow_EM_cluster.at(em));
1054 
1055       pflow->set_px(tlv.Px());
1056       pflow->set_py(tlv.Py());
1057       pflow->set_pz(tlv.Pz());
1058       pflow->set_e(tlv.E());
1059       pflow->set_eclusters(eclus);
1060       pflow->set_hcluster(nullptr);
1061       pflow->set_track(nullptr);
1062       pflow->set_id(global_pflow_index);
1063       pflow->set_type(ParticleFlowElement::PFLOWTYPE::LEFTOVER_EM_PARTICLE);
1064 
1065       pflowContainer->AddParticleFlowElement(global_pflow_index, pflow);
1066       global_pflow_index++;
1067     }
1068 
1069   }  // close EM loop
1070 
1071   // now remove unmatched elements
1072 
1073   if (Verbosity() > 2)
1074   {
1075     std::cout << "ParticleFlowReco::process_event : remove TRK-unlinked EMs and HADs " << std::endl;
1076   }
1077 
1078   for (unsigned int em = 0; em < _pflow_EM_E.size(); em++)
1079   {
1080     // only consider EMs withOUT matched tracks ... we have dealt with the matched cases above
1081     if (!_pflow_EM_match_TRK.at(em).empty())
1082     {
1083       continue;
1084     }
1085 
1086     if (Verbosity() > 5)
1087     {
1088       std::cout << " unmatched EM " << em << " with E / eta / phi = " << _pflow_EM_E.at(em) << " / " << _pflow_EM_eta.at(em) << " / " << _pflow_EM_phi.at(em) << std::endl;
1089     }
1090 
1091     // add PFlow element for this EM
1092     ParticleFlowElement *pflow = new ParticleFlowElementv1();
1093 
1094     // assume massless, could be updated to use K0L
1095     TLorentzVector tlv;
1096     tlv.SetPtEtaPhiM(_pflow_EM_E[em] / cosh(_pflow_EM_eta[em]), _pflow_EM_eta[em], _pflow_EM_phi[em], 0);
1097 
1098     std::vector<RawCluster *> eclus;
1099     eclus.push_back(_pflow_EM_cluster.at(em));
1100 
1101     pflow->set_px(tlv.Px());
1102     pflow->set_py(tlv.Py());
1103     pflow->set_pz(tlv.Pz());
1104     pflow->set_e(tlv.E());
1105     pflow->set_eclusters(eclus);
1106     pflow->set_hcluster(nullptr);
1107     pflow->set_track(nullptr);
1108     pflow->set_id(global_pflow_index);
1109     pflow->set_type(ParticleFlowElement::PFLOWTYPE::UNMATCHED_EM_PARTICLE);
1110 
1111     pflowContainer->AddParticleFlowElement(global_pflow_index, pflow);
1112     global_pflow_index++;
1113 
1114   }  // close EM loop
1115 
1116   for (unsigned int had = 0; had < _pflow_HAD_E.size(); had++)
1117   {
1118     // only consider HADs withOUT matched tracks ... we have dealt with the matched cases above
1119     if (!_pflow_HAD_match_TRK.at(had).empty())
1120     {
1121       continue;
1122     }
1123 
1124     if (Verbosity() > 5)
1125     {
1126       std::cout << " unmatched HAD " << had << " with E / eta / phi = " << _pflow_HAD_E.at(had) << " / " << _pflow_HAD_eta.at(had) << " / " << _pflow_HAD_phi.at(had) << std::endl;
1127     }
1128 
1129     // add PFlow element for this HAD
1130     ParticleFlowElement *pflow = new ParticleFlowElementv1();
1131 
1132     // assume massless, could be updated to use K0L
1133     TLorentzVector tlv;
1134     tlv.SetPtEtaPhiM(_pflow_HAD_E[had] / cosh(_pflow_HAD_eta[had]), _pflow_HAD_eta[had], _pflow_HAD_phi[had], 0);
1135 
1136     pflow->set_px(tlv.Px());
1137     pflow->set_py(tlv.Py());
1138     pflow->set_pz(tlv.Pz());
1139     pflow->set_e(tlv.E());
1140     pflow->set_track(nullptr);
1141     pflow->set_eclusters(std::vector<RawCluster *>());
1142     pflow->set_hcluster(_pflow_HAD_cluster.at(had));
1143     pflow->set_id(global_pflow_index);
1144     pflow->set_type(ParticleFlowElement::PFLOWTYPE::UNMATCHED_NEUTRAL_HADRON);
1145 
1146     pflowContainer->AddParticleFlowElement(global_pflow_index, pflow);
1147     global_pflow_index++;
1148 
1149   }  // close HAD loop
1150 
1151   for (unsigned int trk = 0; trk < _pflow_TRK_p.size(); trk++)
1152   {
1153     // only consider TRKs withOUT matched EM or HAD
1154     if (!_pflow_TRK_match_EM.at(trk).empty() || !_pflow_TRK_match_HAD.at(trk).empty())
1155     {
1156       continue;
1157     }
1158 
1159     if (Verbosity() > 5)
1160     {
1161       std::cout << " unmatched TRK " << trk << " with p / eta / phi = " << _pflow_TRK_p.at(trk) << " / " << _pflow_TRK_eta.at(trk) << " / " << _pflow_TRK_phi.at(trk) << std::endl;
1162     }
1163 
1164     // add PFlow element for this TRK
1165     ParticleFlowElement *pflow = new ParticleFlowElementv1();
1166 
1167     // assume massless, could be updated to use K0L
1168     TLorentzVector tlv;
1169     tlv.SetPtEtaPhiM(_pflow_TRK_p[trk] / cosh(_pflow_TRK_eta[trk]), _pflow_TRK_eta[trk], _pflow_TRK_phi[trk], 0.135);
1170 
1171     pflow->set_px(tlv.Px());
1172     pflow->set_py(tlv.Py());
1173     pflow->set_pz(tlv.Pz());
1174     pflow->set_e(tlv.E());
1175     pflow->set_track(_pflow_TRK_trk.at(trk));
1176     pflow->set_eclusters(std::vector<RawCluster *>());
1177     pflow->set_hcluster(nullptr);
1178     pflow->set_id(global_pflow_index);
1179     pflow->set_type(ParticleFlowElement::PFLOWTYPE::UNMATCHED_CHARGED_HADRON);
1180 
1181     pflowContainer->AddParticleFlowElement(global_pflow_index, pflow);
1182     global_pflow_index++;
1183 
1184   }  // close TRK loop
1185 
1186   // DEBUG: print out all PFLow elements
1187   if (Verbosity() > 5)
1188   {
1189     std::cout << "ParticleFlowReco::process_event: summary of PFlow objects " << std::endl;
1190 
1191     ParticleFlowElementContainer::ConstRange begin_end = pflowContainer->getParticleFlowElements();
1192     for (ParticleFlowElementContainer::ConstIterator hiter = begin_end.first; hiter != begin_end.second; ++hiter)
1193     {
1194       hiter->second->identify();
1195     }
1196   }
1197 
1198   return Fun4AllReturnCodes::EVENT_OK;
1199 }
1200 
1201 int ParticleFlowReco::CreateNode(PHCompositeNode *topNode)
1202 {
1203   PHNodeIterator iter(topNode);
1204 
1205   // Looking for the DST node
1206   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
1207   if (!dstNode)
1208   {
1209     std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
1210     return Fun4AllReturnCodes::ABORTRUN;
1211   }
1212 
1213   // store the PFlow elements under a sub-node directory
1214   PHCompositeNode *pflowNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "PARTICLEFLOW"));
1215   if (!pflowNode)
1216   {
1217     pflowNode = new PHCompositeNode("PARTICLEFLOW");
1218     dstNode->addNode(pflowNode);
1219   }
1220 
1221   // create the ParticleFlowElementContainer node...
1222   ParticleFlowElementContainer *pflowElementContainer = findNode::getClass<ParticleFlowElementContainer>(topNode, "ParticleFlowElements");
1223   if (!pflowElementContainer)
1224   {
1225     pflowElementContainer = new ParticleFlowElementContainer();
1226     PHIODataNode<PHObject> *pflowElementNode = new PHIODataNode<PHObject>(pflowElementContainer, "ParticleFlowElements", "PHObject");
1227     pflowNode->addNode(pflowElementNode);
1228   }
1229   else
1230   {
1231     std::cout << PHWHERE << "::ERROR - ParticleFlowElements node alerady exists, but should not" << std::endl;
1232     exit(-1);
1233   }
1234 
1235   return Fun4AllReturnCodes::EVENT_OK;
1236 }
[ 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!