Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2026-10-15 08:21:05 

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 #include <globalvertex/SvtxVertex.h>
0009 #include <globalvertex/SvtxVertexMap.h>
0010 
0011 #include <calobase/RawCluster.h>
0012 #include <calobase/RawClusterContainer.h>
0013 #include <calobase/RawClusterUtility.h>
0014 #include <calobase/RawTowerGeom.h>
0015 #include <calobase/RawTowerGeomContainer.h>
0016 
0017 #include <trackbase_historic/SvtxTrack.h>
0018 #include <trackbase_historic/SvtxTrackMap.h>
0019 #include <trackbase_historic/SvtxTrackState.h>
0020 
0021 #include <fun4all/Fun4AllReturnCodes.h>
0022 
0023 #include <phool/PHCompositeNode.h>
0024 #include <phool/PHRandomSeed.h>
0025 #include <phool/getClass.h>
0026 
0027 #include <CLHEP/Vector/ThreeVector.h>
0028 
0029 #include <TLorentzVector.h>
0030 #include <TVectorD.h>
0031 
0032 #include <gsl/gsl_randist.h>
0033 #include <gsl/gsl_rng.h>  // for gsl_rng_uniform_pos
0034 
0035 #include <cmath>
0036 #include <iostream>
0037 
0038 // examine second value of std::pair, sort by smallest
0039 bool sort_by_pair_second_lowest(const std::pair<int, float> &a, const std::pair<int, float> &b)
0040 {
0041   return (a.second < b.second);
0042 }
0043 
0044 float ParticleFlowReco::calculate_dR(float eta1, float eta2, float phi1, float phi2)
0045 {
0046   float deta = eta1 - eta2;
0047   float dphi = phi1 - phi2;
0048   while (dphi > M_PI)
0049   {
0050     dphi -= 2 * M_PI;
0051   }
0052   while (dphi < -M_PI)
0053   {
0054     dphi += 2 * M_PI;
0055   }
0056   return sqrt(pow(deta, 2) + pow(dphi, 2));
0057 }
0058 
0059 std::pair<float, float> ParticleFlowReco::get_expected_signature(int trk)
0060 {
0061   float response = (0.553437 + 0.0572246 * log(_pflow_TRK_p[trk])) * _pflow_TRK_p[trk];
0062   float resolution = sqrt(pow(0.119123, 2) + pow(0.312361, 2) / _pflow_TRK_p[trk]) * _pflow_TRK_p[trk];
0063 
0064   std::pair<float, float> expected_signature(response, resolution);
0065 
0066   return expected_signature;
0067 }
0068 
0069 //____________________________________________________________________________..
0070 ParticleFlowReco::ParticleFlowReco(const std::string &name)
0071   : SubsysReco(name)
0072   , _only_crossing_zero(true)
0073   , _dont_use_global_vertex(false)
0074   , _energy_match_Nsigma(1.5)
0075 {
0076 }
0077 
0078 //____________________________________________________________________________..
0079 int ParticleFlowReco::InitRun(PHCompositeNode *topNode)
0080 {
0081   return CreateNode(topNode);
0082 }
0083 
0084 //____________________________________________________________________________..
0085 int ParticleFlowReco::process_event(PHCompositeNode *topNode)
0086 {
0087   if(Verbosity() > 0)
0088   {
0089     std::cout << "ParticleFlowReco::process_event with Nsigma = " << _energy_match_Nsigma << std::endl;
0090   }
0091   // get handle to pflow node
0092   _pflowContainer = findNode::getClass<ParticleFlowElementContainer>(topNode, "ParticleFlowElements");
0093   if(!_pflowContainer)
0094   {
0095     std::cout << " ERROR -- can't find ParticleFlowElements node after it should have been created" << std::endl;
0096     return Fun4AllReturnCodes::ABORTEVENT;
0097   }
0098 
0099   // used for indexing PFlow elements in container
0100   _global_pflow_index = 0;
0101 
0102   // read in tower geometries
0103   _geomEM = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC");
0104   _geomIH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0105   _geomOH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0106 
0107   if(_isEMCalDetailedGeo)
0108   {
0109     _geomEM = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC_DETAILED");
0110 
0111     if(!_geomEM)
0112     {
0113       std::cout << "ParticleFlowReco::process_event : FATAL ERROR, cannot find TOWERGEOM_CEMC_DETAILED" << std::endl;
0114       return Fun4AllReturnCodes::ABORTEVENT;
0115     }
0116   }
0117   else
0118   {
0119     if(!_geomEM)
0120     {
0121       std::cout << "ParticleFlowReco::process_event : FATAL ERROR, cannot find TOWERGEOM_CEMC" << std::endl;
0122       return Fun4AllReturnCodes::ABORTEVENT;
0123     }
0124   }
0125 
0126   if(!_geomIH)
0127   {
0128     std::cout << "ParticleFlowReco::process_event : FATAL ERROR, cannot find TOWERGEOM_HCALIN" << std::endl;
0129     return Fun4AllReturnCodes::ABORTEVENT;
0130   }
0131 
0132   if(!_geomOH)
0133   {
0134     std::cout << "ParticleFlowReco::process_event : FATAL ERROR, cannot find TOWERGEOM_HCALOUT" << std::endl;
0135     return Fun4AllReturnCodes::ABORTEVENT;
0136   }
0137 
0138   _trackmap = findNode::getClass<SvtxTrackMap>(topNode, _track_map_name);
0139   if(!_trackmap)
0140   {
0141     std::cout << "ParticleFlowReco::process_event : FATAL ERROR, cannot find SvtxTrackMap" << std::endl;
0142     return Fun4AllReturnCodes::ABORTEVENT;
0143   }
0144 
0145   _clustersEM = findNode::getClass<RawClusterContainer>(topNode, "TOPOCLUSTER_EMCAL");
0146 
0147   if(!_clustersEM)
0148   {
0149     std::cout << "ParticleFlowReco::process_event : FATAL ERROR, cannot find cluster container TOPOCLUSTER_EMCAL" << std::endl;
0150     return Fun4AllReturnCodes::ABORTEVENT;
0151   }
0152 
0153   _clustersHAD = findNode::getClass<RawClusterContainer>(topNode, "TOPOCLUSTER_HCAL");
0154 
0155   if(!_clustersHAD)
0156   {
0157     std::cout << "ParticleFlowReco::process_event : FATAL ERROR, cannot find cluster container TOPOCLUSTER_HCAL" << std::endl;
0158     return Fun4AllReturnCodes::ABORTEVENT;
0159   }
0160 
0161   _global_vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0162   if((!_global_vertexmap) && (_dont_use_global_vertex == false))
0163   {
0164     std::cout << "ParticleFlowReco::process_event : FATAL ERROR, cannot find GlobalVertexMap" << std::endl;
0165     return Fun4AllReturnCodes::ABORTEVENT;
0166   }
0167 
0168   _svtx_vertexmap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
0169   if(!_svtx_vertexmap)
0170   {
0171     std::cout << "ParticleFlowReco::process_event : FATAL ERROR, cannot find SvtxVertexMap" << std::endl;
0172     return Fun4AllReturnCodes::ABORTEVENT;
0173   }
0174 
0175 
0176 
0177   ResetRepresentation();
0178 
0179   if(Verbosity() > 2)
0180   {
0181     std::cout << "ParticleFlowReco::process_event : initial population of TRK, EM, HAD objects " << std::endl;
0182   }
0183 
0184   PrepareTracks();
0185 
0186   PrepareEMCalClusters();
0187 
0188   PrepareHCalClusters();
0189 
0190   // BEGIN LINKING STEP
0191 
0192   // Link TRK -> EM (best match, but keep reserve of others), and TRK -> HAD (best match)
0193   if(Verbosity() > 2)
0194   {
0195     std::cout << "ParticleFlowReco::process_event : TRK -> EM and TRK -> HAD linking " << std::endl;
0196   }
0197 
0198   LinkTrackToCaloClusters();
0199 
0200   // EM->HAD linking
0201   if(Verbosity() > 2)
0202   {
0203     std::cout << "ParticleFlowReco::process_event : EM -> HAD linking " << std::endl;
0204   }
0205 
0206   LinkEMCalToHCalClusters();
0207 
0208   // SEQUENTIAL MATCHING: if TRK -> EM and EM -> HAD, ensure that TRK -> HAD
0209   if(Verbosity() > 2)
0210   {
0211     std::cout << "ParticleFlowReco::process_event : sequential TRK -> EM && EM -> HAD ==> TRK -> HAD matching " << std::endl;
0212   }
0213 
0214   //LinkTrackToHCalViaEMCal();
0215 
0216   // TRK->EM->HAD removal
0217   if(Verbosity() > 2)
0218   {
0219     std::cout << "ParticleFlowReco::process_event : resolve TRK(s) + EM(s) -> HAD systems " << std::endl;
0220   }
0221 
0222   ProcessMatchedHCalClusters();
0223 
0224   // TRK->EM removal
0225 
0226   if(Verbosity() > 2)
0227   {
0228     std::cout << "ParticleFlowReco::process_event : resolve TRK(s) -> EM(s) ( + no HAD) systems " << std::endl;
0229   }
0230 
0231   ProcessMatchedEMCalClusters();
0232 
0233   // now remove unmatched elements
0234 
0235   if(Verbosity() > 2)
0236   {
0237     std::cout << "ParticleFlowReco::process_event : remove TRK-unlinked EMs and HADs " << std::endl;
0238   }
0239 
0240   ProcessClusterOnlyMatch();
0241 
0242   ProcessUnmatchedEMCalClusters();
0243 
0244   ProcessUnmatchedHCalClusters();
0245 
0246   ProcessUnmatchedTracks();
0247 
0248   // DEBUG: print out all PFLow elements
0249   if(Verbosity() > 5)
0250   {
0251     std::cout << "ParticleFlowReco::process_event: summary of PFlow objects " << std::endl;
0252 
0253     ParticleFlowElementContainer::ConstRange begin_end = _pflowContainer->getParticleFlowElements();
0254     for(ParticleFlowElementContainer::ConstIterator hiter = begin_end.first; hiter != begin_end.second; ++hiter)
0255     {
0256       hiter->second->identify();
0257     }
0258   }
0259 
0260   return Fun4AllReturnCodes::EVENT_OK;
0261 }
0262 
0263 void ParticleFlowReco::ResetRepresentation()
0264 {
0265   // reset internal particle-flow representation
0266   _pflow_TRK_p.clear();
0267   _pflow_TRK_vertex.clear();
0268   _pflow_TRK_eta.clear();
0269   _pflow_TRK_phi.clear();
0270   _pflow_TRK_match_EM.clear();
0271   _pflow_TRK_match_HAD.clear();
0272   _pflow_TRK_addtl_match_EM.clear();
0273   _pflow_TRK_trk.clear();
0274   _pflow_TRK_projections.clear();
0275   /*
0276   _pflow_TRK_EMproj_x.clear();
0277   _pflow_TRK_EMproj_y.clear();
0278   _pflow_TRK_EMproj_z.clear();
0279   _pflow_TRK_HADproj_x.clear();
0280   _pflow_TRK_HADproj_y.clear();
0281   _pflow_TRK_HADproj_z.clear();
0282   */
0283 
0284   _pflow_EM_E.clear();
0285   _pflow_EM_x.clear();
0286   _pflow_EM_y.clear();
0287   _pflow_EM_z.clear();
0288   _pflow_EM_tower_x.clear();
0289   _pflow_EM_tower_y.clear();
0290   _pflow_EM_tower_z.clear();
0291   _pflow_EM_match_HAD.clear();
0292   _pflow_EM_match_TRK.clear();
0293   _pflow_EM_cluster.clear();
0294 
0295   _pflow_HAD_E.clear();
0296   _pflow_HAD_x.clear();
0297   _pflow_HAD_y.clear();
0298   _pflow_HAD_z.clear();
0299   _pflow_HAD_tower_x.clear();
0300   _pflow_HAD_tower_y.clear();
0301   _pflow_HAD_tower_z.clear();
0302   _pflow_HAD_match_EM.clear();
0303   _pflow_HAD_match_TRK.clear();
0304   _pflow_HAD_cluster.clear();
0305 }
0306 
0307 void ParticleFlowReco::PrepareTracks()
0308 {
0309   for(auto &iter : *_trackmap)
0310   {
0311     SvtxTrack *track = iter.second;
0312 
0313     if(!isAcceptableTrack(track))
0314     {
0315       continue;
0316     }
0317 
0318     if(Verbosity() > 2)
0319     {
0320       std::cout << "Track with p= " << track->get_p() << ", eta / phi = "
0321                 << track->get_eta() << " / " << track->get_phi()
0322                 << std::endl;
0323     }
0324 
0325     _pflow_TRK_trk.push_back(track);
0326     _pflow_TRK_vertex.push_back(track->get_vertex_id());
0327     _pflow_TRK_p.push_back(track->get_p());
0328     _pflow_TRK_eta.push_back(track->get_eta());
0329     _pflow_TRK_phi.push_back(track->get_phi());
0330     _pflow_TRK_match_EM.emplace_back();
0331     _pflow_TRK_match_HAD.emplace_back();
0332     _pflow_TRK_addtl_match_EM.emplace_back();
0333 
0334     std::vector<SvtxTrackState*> track_states(6);
0335 
0336     track_states[0] = track->get_state(_geomEM->get_radius());
0337     track_states[1] = track->get_state(_geomEM->get_radius()+_geomEM->get_thickness());
0338     track_states[2] = track->get_state(_geomIH->get_radius());
0339     track_states[3] = track->get_state(_geomIH->get_radius()+_geomIH->get_thickness());
0340     track_states[4] = track->get_state(_geomOH->get_radius());
0341     track_states[5] = track->get_state(_geomOH->get_radius()+_geomOH->get_thickness());
0342 
0343     /// Get the track projections. If they failed for some reason, just use the track
0344     /// phi and eta values at the point of closest approach
0345     std::vector<CLHEP::Hep3Vector> this_projection(6);
0346 
0347     for(unsigned int s = 0; s < this_projection.size(); s++)
0348     {
0349       if(track_states[s])
0350       {
0351         this_projection.at(s) = CLHEP::Hep3Vector(track_states[s]->get_x(), track_states[s]->get_y(), track_states[s]->get_z());
0352       }
0353       else
0354       {
0355         this_projection.at(s) = CLHEP::Hep3Vector(NAN, NAN, NAN);
0356       }
0357     }
0358 
0359     _pflow_TRK_projections.push_back(this_projection);
0360   }
0361 }
0362 
0363 void ParticleFlowReco::PrepareEMCalClusters()
0364 {
0365   RawClusterContainer::ConstRange begin_end = _clustersEM->getClusters();
0366   for(RawClusterContainer::ConstIterator hiter = begin_end.first; hiter != begin_end.second; ++hiter)
0367   {
0368     float cluster_E = hiter->second->get_energy();
0369     if(cluster_E < 0.2)
0370     {
0371       continue;
0372     }
0373 
0374     _pflow_EM_E.push_back(cluster_E);
0375     _pflow_EM_x.push_back(hiter->second->get_x());
0376     _pflow_EM_y.push_back(hiter->second->get_y());
0377     _pflow_EM_z.push_back(hiter->second->get_z());
0378     _pflow_EM_cluster.push_back(hiter->second);
0379     _pflow_EM_match_HAD.emplace_back();
0380     _pflow_EM_match_TRK.emplace_back();
0381 
0382     std::vector<float> this_cluster_tower_x;
0383     std::vector<float> this_cluster_tower_y;
0384     std::vector<float> this_cluster_tower_z;
0385 
0386     // read in towers
0387     RawCluster::TowerConstRange begin_end_towers = hiter->second->get_towers();
0388     for(RawCluster::TowerConstIterator iter = begin_end_towers.first; iter != begin_end_towers.second; ++iter)
0389     {
0390       if(RawTowerDefs::decode_caloid(iter->first) == RawTowerDefs::CalorimeterId::CEMC)
0391       {
0392         RawTowerGeom *tower_geom = _geomEM->get_tower_geometry(iter->first);
0393 
0394         this_cluster_tower_x.push_back(tower_geom->get_center_x());
0395         this_cluster_tower_y.push_back(tower_geom->get_center_y());
0396         this_cluster_tower_z.push_back(tower_geom->get_center_z());
0397       }
0398       else
0399       {
0400         std::cout << "ParticleFlowReco::process_event : FATAL ERROR , EM topoClusters seem to contain HCal towers" << std::endl;
0401       }
0402     }  // close tower loop
0403 
0404     _pflow_EM_tower_x.push_back(this_cluster_tower_x);
0405     _pflow_EM_tower_y.push_back(this_cluster_tower_y);
0406     _pflow_EM_tower_z.push_back(this_cluster_tower_z);
0407 
0408   }  // close cluster loop
0409 }
0410 
0411 void ParticleFlowReco::PrepareHCalClusters()
0412 {
0413   RawClusterContainer::ConstRange begin_end = _clustersHAD->getClusters();
0414   for(RawClusterContainer::ConstIterator hiter = begin_end.first; hiter != begin_end.second; ++hiter)
0415   {
0416     float cluster_E = hiter->second->get_energy();
0417     if(cluster_E < 0.2)
0418     {
0419       continue;
0420     }
0421 
0422     _pflow_HAD_E.push_back(cluster_E);
0423     _pflow_HAD_x.push_back(hiter->second->get_x());
0424     _pflow_HAD_y.push_back(hiter->second->get_y());
0425     _pflow_HAD_z.push_back(hiter->second->get_z());
0426     _pflow_HAD_cluster.push_back(hiter->second);
0427 
0428     _pflow_HAD_match_EM.emplace_back();
0429     _pflow_HAD_match_TRK.emplace_back();
0430 
0431     std::vector<float> this_cluster_tower_x;
0432     std::vector<float> this_cluster_tower_y;
0433     std::vector<float> this_cluster_tower_z;
0434 
0435     // read in towers
0436     RawCluster::TowerConstRange begin_end_towers = hiter->second->get_towers();
0437     for(RawCluster::TowerConstIterator iter = begin_end_towers.first; iter != begin_end_towers.second; ++iter)
0438     {
0439       if(RawTowerDefs::decode_caloid(iter->first) == RawTowerDefs::CalorimeterId::HCALIN)
0440       {
0441         RawTowerGeom *tower_geom = _geomIH->get_tower_geometry(iter->first);
0442 
0443         this_cluster_tower_x.push_back(tower_geom->get_center_x());
0444         this_cluster_tower_y.push_back(tower_geom->get_center_y());
0445         this_cluster_tower_z.push_back(tower_geom->get_center_z());
0446       }
0447 
0448       else if(RawTowerDefs::decode_caloid(iter->first) == RawTowerDefs::CalorimeterId::HCALOUT)
0449       {
0450         RawTowerGeom *tower_geom = _geomOH->get_tower_geometry(iter->first);
0451 
0452         this_cluster_tower_x.push_back(tower_geom->get_center_x());
0453         this_cluster_tower_y.push_back(tower_geom->get_center_y());
0454         this_cluster_tower_z.push_back(tower_geom->get_center_z());
0455       }
0456       else
0457       {
0458         std::cout << "ParticleFlowReco::process_event : FATAL ERROR , HCal topoClusters seem to contain EM towers" << std::endl;
0459       }
0460 
0461     }  // close tower loop
0462 
0463     _pflow_HAD_tower_x.push_back(this_cluster_tower_x);
0464     _pflow_HAD_tower_y.push_back(this_cluster_tower_y);
0465     _pflow_HAD_tower_z.push_back(this_cluster_tower_z);
0466 
0467   }  // close cluster loop
0468 }
0469 
0470 void ParticleFlowReco::LinkTrackToCaloClusters()
0471 {
0472   for(unsigned int trk = 0; trk < _pflow_TRK_p.size(); trk++)
0473   {
0474     if(Verbosity() > 10)
0475     {
0476       std::cout << " TRK " << trk << " with p / eta / phi = " << _pflow_TRK_p[trk] << " / " << _pflow_TRK_eta[trk] << " / " << _pflow_TRK_phi[trk] << std::endl;
0477     }
0478 
0479     //CLHEP::Hep3Vector trackproj_emcal(_pflow_TRK_EMproj_x[trk], _pflow_TRK_EMproj_y[trk], _pflow_TRK_EMproj_z[trk]);
0480 
0481     //CLHEP::Hep3Vector trackproj_hcal(_pflow_TRK_HADproj_x[trk], _pflow_TRK_HADproj_y[trk], _pflow_TRK_HADproj_z[trk]);
0482 
0483     CLHEP::Hep3Vector vertex = GetVertexForCalorimeter();
0484 
0485     // TRK -> EM link
0486     float min_em_dR = 0.2;
0487     int min_em_index = -1;
0488 
0489     for(unsigned int em = 0; em < _pflow_EM_E.size(); em++)
0490     {
0491       CLHEP::Hep3Vector emcal_cluster(_pflow_EM_x[em], _pflow_EM_y[em], _pflow_EM_z[em]);
0492       emcal_cluster = emcal_cluster - vertex;
0493 
0494       CLHEP::Hep3Vector track_proj = Minimize_PhiEta(emcal_cluster, trk, vertex, false);
0495 
0496       //If the projection doesn't reach the calorimeter, we should not try to match
0497       if(std::isnan(track_proj.x()))
0498       {
0499         continue;
0500       }
0501 
0502       float dR = fabs(track_proj.deltaR(emcal_cluster));
0503 
0504       if(dR > 0.2)
0505       {
0506         continue;
0507       }
0508 
0509       bool has_overlap = false;
0510 
0511       if(_use_overlap_match)
0512       {
0513         for(unsigned int tow = 0; tow < _pflow_EM_tower_x.at(em).size(); tow++)
0514         {
0515           CLHEP::Hep3Vector tower(_pflow_EM_tower_x.at(em).at(tow), _pflow_EM_tower_y.at(em).at(tow), _pflow_EM_tower_z.at(em).at(tow));
0516           tower = tower - vertex;
0517 
0518           float deta = track_proj.eta() - tower.eta();
0519           float dphi = track_proj.deltaPhi(tower);
0520 
0521           if(fabs(deta) < 0.025 * _trk_emc_overlap_eta && fabs(dphi) < 0.025 * _trk_emc_overlap_phi)
0522           {
0523             has_overlap = true;
0524             break;
0525           }
0526         }
0527       }
0528       else
0529       {
0530         has_overlap = true;
0531       }
0532 
0533       if(has_overlap)
0534       {
0535         if(Verbosity() > 5)
0536         {
0537           std::cout << " -> possible match to EM " << em << " with dR = " << dR << std::endl;
0538         }
0539 
0540         _pflow_TRK_addtl_match_EM.at(trk).push_back(std::pair<int, float>(em, dR));
0541       }
0542       else
0543       {
0544         if(Verbosity() > 5)
0545         {
0546           std::cout << " -> no match to EM " << em << " (even though dR = " << dR << " )" << std::endl;
0547         }
0548       }
0549     }
0550 
0551     // sort possible matches
0552 
0553     std::sort(_pflow_TRK_addtl_match_EM.at(trk).begin(), _pflow_TRK_addtl_match_EM.at(trk).end(), sort_by_pair_second_lowest);
0554     if(Verbosity() > 10)
0555     {
0556       for(auto &n : _pflow_TRK_addtl_match_EM.at(trk))
0557       {
0558         std::cout << " -> sorted list of matches, EM / dR = " << n.first << " / " << n.second << std::endl;
0559       }
0560     }
0561 
0562     if(_pflow_TRK_addtl_match_EM.at(trk).size() > 0)
0563     {
0564       min_em_index = _pflow_TRK_addtl_match_EM.at(trk).at(0).first;
0565       min_em_dR = _pflow_TRK_addtl_match_EM.at(trk).at(0).second;
0566       // delete best matched element
0567       _pflow_TRK_addtl_match_EM.at(trk).erase(_pflow_TRK_addtl_match_EM.at(trk).begin());
0568     }
0569 
0570     if(min_em_index > -1)
0571     {
0572       _pflow_EM_match_TRK.at(min_em_index).push_back(trk);
0573       _pflow_TRK_match_EM.at(trk).push_back(min_em_index);
0574 
0575       if(Verbosity() > 5)
0576       {
0577         std::cout << " -> matched EM " << min_em_index << " with pt / x / y / z = " << _pflow_EM_E.at(min_em_index) << " / " << _pflow_EM_x.at(min_em_index) << " / " << _pflow_EM_y.at(min_em_index) << " / " << _pflow_EM_z.at(min_em_index) << ", dR = " << min_em_dR;
0578         std::cout << " ( " << _pflow_TRK_addtl_match_EM.at(trk).size() << " other possible matches ) " << std::endl;
0579       }
0580     }
0581     else
0582     {
0583       if(Verbosity() > 5)
0584       {
0585         std::cout << " -> no EM match! ( best dR = " << min_em_dR << " ) " << std::endl;
0586       }
0587     }
0588 
0589     // TRK -> HAD link
0590     float min_had_dR = 0.2;
0591     int min_had_index = -1;
0592     float max_had_pt = 0;
0593 
0594     // TODO: sequential linking should better happen here -- i.e. allow EM-matched HAD's into the possible pool
0595     for(unsigned int had = 0; had < _pflow_HAD_E.size(); had++)
0596     {
0597       CLHEP::Hep3Vector hcal_cluster(_pflow_HAD_x[had], _pflow_HAD_y[had], _pflow_HAD_z[had]);
0598       hcal_cluster = hcal_cluster - vertex;
0599 
0600       CLHEP::Hep3Vector track_proj = Minimize_PhiEta(hcal_cluster, trk, vertex, true);
0601 
0602       //If the projection doesn't reach the calorimeter, we should not try to match
0603       if(std::isnan(track_proj.x()))
0604       {
0605         continue;
0606       }
0607 
0608       float dR = fabs(track_proj.deltaR(hcal_cluster));
0609 
0610       if(dR > 0.5)
0611       {
0612         continue;
0613       }
0614 
0615       bool has_overlap = false;
0616 
0617       if(_use_overlap_match)
0618       {
0619         for(unsigned int tow = 0; tow < _pflow_HAD_tower_x.at(had).size(); tow++)
0620         {
0621           CLHEP::Hep3Vector tower(_pflow_HAD_tower_x.at(had).at(tow), _pflow_HAD_tower_y.at(had).at(tow), _pflow_HAD_tower_z.at(had).at(tow));
0622           tower = tower - vertex;
0623 
0624           float deta = track_proj.eta() - tower.eta();
0625           float dphi = track_proj.deltaPhi(tower);
0626 
0627           if(fabs(deta) < 0.1 * _trk_had_overlap_eta && fabs(dphi) < 0.1 * _trk_had_overlap_phi)
0628           {
0629             has_overlap = true;
0630             break;
0631           }
0632         }
0633       }
0634       else
0635       {
0636         has_overlap = true;
0637       }
0638 
0639       if(has_overlap)
0640       {
0641         if(Verbosity() > 5)
0642         {
0643           std::cout << " -> possible match to HAD " << had << " with dR = " << dR << std::endl;
0644         }
0645 
0646         if(_pflow_HAD_E.at(had) > max_had_pt)
0647         {
0648           max_had_pt = _pflow_HAD_E.at(had);
0649           min_had_index = had;
0650           min_had_dR = dR;
0651         }
0652       }
0653       else
0654       {
0655         if(Verbosity() > 5)
0656         {
0657           std::cout << " -> no match to HAD " << had << " (even though dR = " << dR << " )" << std::endl;
0658         }
0659       }
0660     }
0661 
0662     if(min_had_index > -1)
0663     {
0664       _pflow_HAD_match_TRK.at(min_had_index).push_back(trk);
0665       _pflow_TRK_match_HAD.at(trk).push_back(min_had_index);
0666 
0667       if(Verbosity() > 5)
0668       {
0669         std::cout << " -> matched HAD " << min_had_index << " with E / x / y / z = " << _pflow_HAD_E.at(min_had_index) << " / " << _pflow_HAD_x.at(min_had_index) << " / " << _pflow_HAD_y.at(min_had_index) << " / " << _pflow_HAD_z.at(min_had_index) << ", dR = " << min_had_dR << std::endl;
0670       }
0671     }
0672     else
0673     {
0674       if(Verbosity() > 5)
0675       {
0676         std::cout << " -> no HAD match! ( best dR = " << min_had_dR << " ) " << std::endl;
0677       }
0678     }
0679   }
0680 }
0681 
0682 void ParticleFlowReco::LinkEMCalToHCalClusters()
0683 {
0684 
0685   for(unsigned int had = 0; had < _pflow_HAD_E.size(); had++)
0686   {
0687     CLHEP::Hep3Vector hcal_cluster(_pflow_HAD_x[had], _pflow_HAD_y[had], _pflow_HAD_z[had]);
0688 
0689     for(unsigned int em = 0; em < _pflow_EM_E.size(); em++)
0690     {
0691       CLHEP::Hep3Vector emcal_cluster(_pflow_EM_x[em], _pflow_EM_y[em], _pflow_EM_z[em]);
0692 
0693       if(Verbosity() > 10)
0694       {
0695         std::cout << " EM with E / x / y / z = " << _pflow_EM_E[em] << " / " << _pflow_EM_x[em] << " / " << _pflow_EM_y[em] << " / " << _pflow_EM_z[em] << std::endl;
0696       }
0697 
0698       float dR = fabs(emcal_cluster.deltaR(hcal_cluster));
0699 
0700       if(dR > 0.5)
0701       {
0702         continue;
0703       }
0704 
0705       bool has_overlap = false;
0706 
0707       if(_use_overlap_match)
0708       {
0709         for(unsigned int tow = 0; tow < _pflow_EM_tower_x.at(em).size(); tow++)
0710         {
0711           CLHEP::Hep3Vector tower(_pflow_EM_tower_x.at(em).at(tow), _pflow_EM_tower_y.at(em).at(tow), _pflow_EM_tower_z.at(em).at(tow));
0712 
0713           float deta = hcal_cluster.eta() - tower.eta();
0714           float dphi = hcal_cluster.deltaPhi(tower);
0715 
0716           if(fabs(deta) < 0.1 * _emc_had_overlap_eta && fabs(dphi) < 0.1 * _emc_had_overlap_phi)
0717           {
0718             has_overlap = true;
0719             break;
0720           }
0721         }
0722       }
0723       else
0724       {
0725         has_overlap = true;
0726       }
0727 
0728       if(has_overlap)
0729       {
0730         if(Verbosity() > 5)
0731         {
0732           std::cout << " -> possible EM match to HAD " << em << " with dR = " << dR << std::endl;
0733         }
0734 /*
0735         if(dR < min_dR) //We are going to look for neutral hadrons, match the closest cluster
0736         {
0737           min_em_index = em;
0738           min_dR = dR;
0739         }
0740         */
0741         //We associate all EM matched to HAD, later, after track matching, we selected the closest to HAD
0742         _pflow_HAD_match_EM.at(had).push_back(em);
0743         _pflow_EM_match_HAD.at(em).push_back(had);
0744       }
0745       else
0746       {
0747         if(Verbosity() > 5)
0748         {
0749           std::cout << " -> no match to HAD " << had << " (even though dR = " << dR << " )" << std::endl;
0750         }
0751       }
0752     }
0753 /*
0754     if(min_em_index > -1)
0755     {
0756       _pflow_HAD_match_EM.at(had).push_back(min_em_index);
0757       _pflow_EM_match_HAD.at(min_em_index).push_back(had);
0758 
0759       if(Verbosity() > 5)
0760       {
0761         std::cout << " -> matched EM to HAD with E / x / y / z = " << _pflow_EM_E.at(min_em_index) << " / " << _pflow_EM_x.at(min_em_index) << " / " << _pflow_EM_y.at(min_em_index) << " / " << _pflow_EM_z.at(min_em_index) << ", dR = " << min_dR << std::endl;
0762       }
0763     }
0764     else
0765     {
0766       if(Verbosity() > 5)
0767       {
0768         std::cout << " -> no EM match! ( best dR = " << min_dR << " ) " << std::endl;
0769       }
0770     }
0771     */
0772   }
0773 }
0774 
0775 void ParticleFlowReco::LinkTrackToHCalViaEMCal()
0776 {
0777   for(unsigned int trk = 0; trk < _pflow_TRK_p.size(); trk++)
0778   {
0779     // go through all matched EMs
0780     for(unsigned int i = 0; i < _pflow_TRK_match_EM.at(trk).size(); i++)
0781     {
0782       int em = _pflow_TRK_match_EM.at(trk).at(i);
0783 
0784       // if this EM has a matched HAD...
0785       for(unsigned int j = 0; j < _pflow_EM_match_HAD.at(em).size(); j++)
0786       {
0787         int had = _pflow_EM_match_HAD.at(em).at(j);
0788 
0789         // and the TRK is NOT matched to this HAD...
0790         bool is_trk_matched_to_HAD = false;
0791         for(int existing_had : _pflow_TRK_match_HAD.at(trk))
0792         {
0793           if(had == existing_had)
0794           {
0795             is_trk_matched_to_HAD = true;
0796           }
0797         }
0798 
0799         // if this is the case, create TRK->HAD link
0800         if(!is_trk_matched_to_HAD)
0801         {
0802           _pflow_TRK_match_HAD.at(trk).push_back(had);
0803           _pflow_HAD_match_TRK.at(had).push_back(trk);
0804 
0805           if(Verbosity() > 5)
0806           {
0807             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;
0808             std::cout << " -> sequential match to HAD " << had << " through EM " << j << std::endl;
0809           }
0810         }
0811 
0812       }  // close the HAD loop
0813 
0814     }  // close the EM loop
0815 
0816   }  // close the TRK loop
0817 }
0818 
0819 void ParticleFlowReco::ProcessMatchedHCalClusters()
0820 {
0821   for(unsigned int had = 0; had < _pflow_HAD_E.size(); had++)
0822   {
0823     // only consider HAD with matched tracks ... others we will deal with later
0824     if(_pflow_HAD_match_TRK.at(had).size() == 0)
0825     {
0826       continue;
0827     }
0828 
0829     if(Verbosity() > 5)
0830     {
0831       std::cout << " HAD " << had << " with E / x / y / z = " << _pflow_HAD_E.at(had) << " / " << _pflow_HAD_x.at(had) << " / " << _pflow_HAD_y.at(had) << " / " << _pflow_HAD_z.at(had) << std::endl;
0832     }
0833 
0834     // setup for Sum-pT^trk -> calo prediction
0835     float total_TRK_p = 0;
0836     float total_expected_E = 0;
0837     float total_expected_E_var = 0;
0838 
0839     // begin with this HAD calo energy
0840     float total_EMHAD_E = _pflow_HAD_E.at(had);
0841 
0842     std::vector<RawCluster*> matchedEClusters_this_had_via_trk;
0843 
0844     //Remove later
0845     //std::vector<int> matchedEClustersIDs;
0846 /*
0847     // iterate over the EMs matched to this HAD
0848     for(int em : _pflow_HAD_match_EM.at(had))
0849     {
0850       // ensure there is at least one track matched to this EM
0851       if(_pflow_EM_match_TRK.at(em).size() == 0)
0852       {
0853         continue;
0854       }
0855 
0856       // add it to the total calo E
0857       total_EMHAD_E += _pflow_EM_E.at(em);
0858       matchedEClusters.push_back(_pflow_EM_cluster.at(em));
0859       //matchedEClustersIDs.push_back(em);
0860       if(Verbosity() > 5)
0861       {
0862         std::cout << " -> -> LINKED EM " << em << " with E / x / y / z = " << _pflow_EM_E.at(em) << " / " << _pflow_EM_x.at(em) << " / " << _pflow_EM_y.at(em) << " / " << _pflow_EM_z.at(em) << std::endl;
0863       }
0864     }
0865 */
0866     //In case of EM leftover
0867     std::vector<SvtxTrack*> assoc_tracks;
0868 
0869     // iterate over the TRKs matched to this HAD
0870     for(unsigned int j = 0; j < _pflow_HAD_match_TRK.at(had).size(); j++)
0871     {
0872       int trk = _pflow_HAD_match_TRK.at(had).at(j);
0873 
0874       if(Verbosity() > 5)
0875       {
0876         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;
0877       }
0878 
0879       total_TRK_p += _pflow_TRK_p.at(trk);
0880 
0881       std::pair<float, float> expected_signature = get_expected_signature(trk);
0882 
0883       float expected_E_mean = expected_signature.first;
0884       float expected_E_sigma = expected_signature.second;
0885 
0886       if(Verbosity() > 5)
0887       {
0888         std::cout << " -> -> -> expected calo signature is " << expected_E_mean << " +/- " << expected_E_sigma << std::endl;
0889       }
0890 
0891       total_expected_E += expected_E_mean;
0892       total_expected_E_var += pow(expected_E_sigma, 2);
0893 
0894       std::vector<RawCluster *> matchedEClusters_this_track;
0895 
0896       for(unsigned int em = 0; em < _pflow_TRK_match_EM.at(trk).size(); em++)
0897       {
0898         total_EMHAD_E += _pflow_EM_E.at(_pflow_TRK_match_EM.at(trk).at(em));
0899         matchedEClusters_this_track.push_back(_pflow_EM_cluster.at(_pflow_TRK_match_EM.at(trk).at(em)));
0900         if(std::find(matchedEClusters_this_had_via_trk.begin(), matchedEClusters_this_had_via_trk.end(), _pflow_EM_cluster.at(_pflow_TRK_match_EM.at(trk).at(em))) == matchedEClusters_this_had_via_trk.end())
0901         {
0902           matchedEClusters_this_had_via_trk.push_back(_pflow_EM_cluster.at(_pflow_TRK_match_EM.at(trk).at(em)));
0903         }
0904       }
0905 
0906       // add PFlow element for each track
0907       ParticleFlowElement *pflow = new ParticleFlowElementv1();
0908 
0909       // assume pion mass
0910       float pfE = std::sqrt((_pflow_TRK_p[trk]*_pflow_TRK_p[trk]) + (0.135*0.135));
0911 
0912       CLHEP::Hep3Vector pflow_vector;
0913       pflow_vector.setREtaPhi(_pflow_TRK_p[trk], _pflow_TRK_eta[trk], _pflow_TRK_phi[trk]);
0914 
0915       pflow->set_px(pflow_vector.getX());
0916       pflow->set_py(pflow_vector.getY());
0917       pflow->set_pz(pflow_vector.getZ());
0918       pflow->set_e(pfE);
0919       pflow->set_track(_pflow_TRK_trk[trk]);
0920       pflow->set_eclusters(matchedEClusters_this_track);
0921       pflow->set_hcluster(_pflow_HAD_cluster.at(had));
0922       pflow->set_id(_global_pflow_index);
0923       pflow->set_type(ParticleFlowElement::PFLOWTYPE::MATCHED_CHARGED_HADRON);
0924 /*
0925       std::cout << "Creating particle flow element" << std::endl;
0926       std::cout << "PF ID: " << _global_pflow_index << " as MATCHED_CHARGED_HADRON" << std::endl;
0927       std::cout << "Track ID: " << trk << std::endl;
0928       for(unsigned int iemc = 0; iemc < matchedEClustersIDs.size(); iemc++)
0929       {
0930         std::cout << "EMCal ID: " << matchedEClustersIDs.at(iemc) << std::endl;
0931       }
0932       std::cout << "HCal ID: " << had << std::endl;
0933 */
0934 
0935       _pflowContainer->AddParticleFlowElement(_global_pflow_index, pflow);
0936       _global_pflow_index++;
0937       assoc_tracks.push_back(_pflow_TRK_trk[trk]);
0938     }
0939     // Track + E+HCal PF elements are created
0940 
0941     // process compatibility of fit
0942     float total_expected_E_err = sqrt(total_expected_E_var);
0943 
0944     if(Verbosity() > 5)
0945     {
0946       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;
0947     }
0948 
0949     // if Sum pT > calo, add in additional possible matched EMs associated with tracks until that is no longer the case
0950     //Remove later
0951     //std::vector<int> additional_EMs_IDs;
0952 /*
0953     if(total_expected_E > total_EMHAD_E)
0954     {
0955       if(Verbosity() > 5)
0956       {
0957         std::cout << " -> Expected E > Observed E, looking for additional potential TRK->EM matches" << std::endl;
0958       }
0959 
0960       std::map<int, float> additional_EMs;
0961 
0962       for(int trk : _pflow_HAD_match_TRK.at(had))
0963       {
0964         int addtl_matches = _pflow_TRK_addtl_match_EM.at(trk).size();
0965 
0966         if(Verbosity() > 10)
0967         {
0968           std::cout << " -> -> TRK " << trk << " has " << addtl_matches << " additional matches! " << std::endl;
0969         }
0970 
0971         for(auto &n : _pflow_TRK_addtl_match_EM.at(trk))
0972         {
0973           if(Verbosity() > 10)
0974           {
0975             std::cout << " -> -> -> additional match to EM = " << n.first << " with dR = " << n.second << std::endl;
0976           }
0977 
0978           float existing_dR = 0.21;
0979           int counts = additional_EMs.count(n.first);
0980           if(counts > 0)
0981           {
0982             existing_dR = additional_EMs[n.first];
0983           }
0984           if(n.second < existing_dR)
0985           {
0986             additional_EMs[n.first] = n.second;
0987           }
0988         }
0989       }
0990 
0991       // map now assured to have only minimal dR values for each possible additional EM
0992       // translate the map to a vector of pairs, then sort by smallest dR
0993 
0994       std::vector<std::pair<int, float> > additional_EMs_vec;
0995 
0996       additional_EMs_vec.reserve(additional_EMs.size());
0997       for(auto &x : additional_EMs)
0998       {
0999         additional_EMs_vec.emplace_back(x.first, x.second);
1000       }
1001 
1002       std::sort(additional_EMs_vec.begin(), additional_EMs_vec.end(), sort_by_pair_second_lowest);
1003 
1004       if(Verbosity() > 5)
1005       {
1006         std::cout << " -> Sorting the set of potential additional EMs " << std::endl;
1007       }
1008 
1009       // now add in additional EMs until there are none left or it is no longer the case that Sum pT > calo
1010 
1011       int n_EM_added = 0;
1012       while (additional_EMs_vec.size() != 0 && total_expected_E > total_EMHAD_E)
1013       {
1014         int new_EM = additional_EMs_vec.at(0).first;
1015 
1016         if(Verbosity() > 5)
1017         {
1018           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;
1019         }
1020 
1021         // for now, just make the first HAD-linked track point to this new EM, and vice versa
1022         _pflow_EM_match_TRK.at(new_EM).push_back(_pflow_HAD_match_TRK.at(had).at(0));
1023         _pflow_TRK_match_EM.at(_pflow_HAD_match_TRK.at(had).at(0)).push_back(new_EM);
1024 
1025         // add to expected calo
1026         total_EMHAD_E += _pflow_EM_E.at(new_EM);
1027 
1028         // erase lowest-dR EM
1029         additional_EMs_vec.erase(additional_EMs_vec.begin());
1030 
1031         //additional_EMs_IDs.push_back(new_EM);
1032 
1033         n_EM_added++;
1034       }
1035 
1036       if(Verbosity() > 5)
1037       {
1038         if(n_EM_added > 0)
1039         {
1040           std::cout << "After adding N = " << n_EM_added << " any additional EMs : " << std::endl;
1041           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;
1042         }
1043         else
1044         {
1045           std::cout << "No additional EMs found, continuing hypothesis check" << std::endl;
1046         }
1047       }
1048     }
1049 */
1050     if(total_expected_E + _energy_match_Nsigma * total_expected_E_err > total_EMHAD_E)
1051     {
1052       if(Verbosity() > 5)
1053       {
1054         std::cout << " -> -> calo compatible within Nsigma = " << _energy_match_Nsigma << " , remove and keep tracks " << std::endl;
1055       }
1056 
1057       // PFlow elements already created from tracks above, no more needs to be done
1058     }
1059     else
1060     {
1061       float residual_energy = total_EMHAD_E - total_expected_E;
1062 
1063       if(Verbosity() > 5)
1064       {
1065         std::cout << " -> -> calo not compatible, create leftover cluster with " << residual_energy << std::endl;
1066       }
1067 
1068       // create additional PFlow element (tracks already created above)
1069       ParticleFlowElement *pflow = new ParticleFlowElementv1();
1070 
1071       // assume no mass, but could update to use K0L mass(?)
1072       CLHEP::Hep3Vector pflow_eta_phi(_pflow_HAD_x[had], _pflow_HAD_y[had], _pflow_HAD_z[had]);
1073 
1074       CLHEP::Hep3Vector pflow_vector;
1075       pflow_vector.setREtaPhi(residual_energy, pflow_eta_phi.eta(), pflow_eta_phi.phi());
1076 
1077       pflow->set_px(pflow_vector.getX());
1078       pflow->set_py(pflow_vector.getY());
1079       pflow->set_pz(pflow_vector.getZ());
1080       pflow->set_e(residual_energy);
1081       pflow->set_track(nullptr);
1082       pflow->set_eclusters(matchedEClusters_this_had_via_trk);
1083       pflow->set_hcluster(_pflow_HAD_cluster.at(had));
1084       pflow->set_assoc_tracks(assoc_tracks);
1085       pflow->set_id(_global_pflow_index);
1086       pflow->set_type(ParticleFlowElement::PFLOWTYPE::LEFTOVER_EM_PARTICLE);
1087 /*
1088       std::cout << "Creating particle flow element" << std::endl;
1089       std::cout << "PF ID: " << _global_pflow_index << " as LEFTOVER_EM_PARTICLE" << std::endl;
1090       std::cout << "Track ID: " << "XXX" << std::endl;
1091       for(unsigned int iemc = 0; iemc < matchedEClustersIDs.size(); iemc++)
1092       {
1093         std::cout << "EMCal ID: " << matchedEClustersIDs.at(iemc) << std::endl;
1094       }
1095       for(unsigned int iemc = 0; iemc < additional_EMs_IDs.size(); iemc++)
1096       {
1097         std::cout << "Additional EMCal ID: " << additional_EMs_IDs.at(iemc) << std::endl;
1098       }
1099       std::cout << "HCal ID: " << had << std::endl;
1100 */
1101       _pflowContainer->AddParticleFlowElement(_global_pflow_index, pflow);
1102       _global_pflow_index++;
1103     }
1104 
1105   }  // close HAD loop
1106 }
1107 
1108 void ParticleFlowReco::ProcessMatchedEMCalClusters()
1109 {
1110   for(unsigned int em = 0; em < _pflow_EM_E.size(); em++)
1111   {
1112     // only consider EM with matched tracks, but no matched HADs
1113     if(_pflow_EM_match_HAD.at(em).size() != 0)
1114     {
1115       continue;
1116     }
1117     if(_pflow_EM_match_TRK.at(em).size() == 0)
1118     {
1119       continue;
1120     }
1121 
1122     if(Verbosity() > 5)
1123     {
1124       std::cout << " EM " << em << " with E / x / y / z = " << _pflow_EM_E.at(em) << " / " << _pflow_EM_x.at(em) << " / " << _pflow_EM_y.at(em) << " / " << _pflow_EM_z.at(em) << std::endl;
1125     }
1126 
1127     // setup for Sum-pT^trk -> calo prediction
1128     float total_TRK_p = 0;
1129     float total_expected_E = 0;
1130     float total_expected_E_var = 0;
1131 
1132     // begin with this EM calo energy
1133     float total_EM_E = _pflow_EM_E.at(em);
1134 
1135     //In case of EM leftover
1136     std::vector<SvtxTrack*> assoc_tracks;
1137 
1138     // iterate over the TRKs matched to this EM
1139     for(unsigned int j = 0; j < _pflow_EM_match_TRK.at(em).size(); j++)
1140     {
1141       int trk = _pflow_EM_match_TRK.at(em).at(j);
1142 
1143       //Exclude tracks that were already matched to HADs
1144       if(_pflow_TRK_match_HAD.at(trk).size() != 0)
1145       {
1146         continue;
1147       }
1148 
1149       if(Verbosity() > 5)
1150       {
1151         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;
1152       }
1153 
1154       total_TRK_p += _pflow_TRK_p.at(trk);
1155 
1156       std::pair<float, float> expected_signature = get_expected_signature(trk);
1157 
1158       float expected_E_mean = expected_signature.first;
1159       float expected_E_sigma = expected_signature.second;
1160 
1161       if(Verbosity() > 5)
1162       {
1163         std::cout << " -> -> -> expected calo signature is " << expected_E_mean << " +/- " << expected_E_sigma << std::endl;
1164       }
1165 
1166       total_expected_E += expected_E_mean;
1167       total_expected_E_var += pow(expected_E_sigma, 2);
1168 
1169       // add PFlow element for each track
1170       ParticleFlowElement *pflow = new ParticleFlowElementv1();
1171 
1172       // assume pion mass
1173       float pfE = std::sqrt((_pflow_TRK_p[trk]*_pflow_TRK_p[trk]) + (0.135*0.135));
1174 
1175       CLHEP::Hep3Vector pflow_vector;
1176       pflow_vector.setREtaPhi(_pflow_TRK_p[trk], _pflow_TRK_eta[trk], _pflow_TRK_phi[trk]);
1177 
1178       std::vector<RawCluster *> eclus;
1179       eclus.push_back(_pflow_EM_cluster.at(em));
1180 
1181       pflow->set_px(pflow_vector.getX());
1182       pflow->set_py(pflow_vector.getY());
1183       pflow->set_pz(pflow_vector.getZ());
1184       pflow->set_e(pfE);
1185       pflow->set_track(_pflow_TRK_trk.at(trk));
1186       pflow->set_eclusters(eclus);
1187       pflow->set_hcluster(nullptr);
1188       pflow->set_id(_global_pflow_index);
1189       pflow->set_type(ParticleFlowElement::PFLOWTYPE::MATCHED_CHARGED_HADRON);
1190 /*
1191       std::cout << "Creating particle flow element" << std::endl;
1192       std::cout << "PF ID: " << _global_pflow_index << " as MATCHED_CHARGED_HADRON (EMCal+Track)" << std::endl;
1193       std::cout << "Track ID: " << trk << std::endl;
1194       std::cout << "EMCal ID: " << em << std::endl;
1195       std::cout << "HCal ID: " << "XXX" << std::endl;
1196 */
1197       _pflowContainer->AddParticleFlowElement(_global_pflow_index, pflow);
1198       _global_pflow_index++;
1199       assoc_tracks.push_back(_pflow_TRK_trk.at(trk));
1200     }
1201 
1202     // process compatibility of fit
1203     float total_expected_E_err = sqrt(total_expected_E_var);
1204 
1205     if(Verbosity() > 5)
1206     {
1207       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;
1208     }
1209 
1210     if(total_expected_E + _energy_match_Nsigma * total_expected_E_err > total_EM_E)
1211     {
1212       if(Verbosity() > 5)
1213       {
1214         std::cout << " -> -> calo compatible within Nsigma = " << _energy_match_Nsigma << "  , remove and keep tracks " << std::endl;
1215       }
1216 
1217       // PFlow elements already created from tracks above, no more needs to be done
1218     }
1219     else
1220     {
1221       float residual_energy = total_EM_E - total_expected_E;
1222 
1223       if(Verbosity() > 5)
1224       {
1225         std::cout << " -> -> calo not compatible, create leftover cluster with " << residual_energy << std::endl;
1226       }
1227 
1228       // create additional PFlow element (tracks already created above)
1229       ParticleFlowElement *pflow = new ParticleFlowElementv1();
1230 
1231       // assume no mass, but could update to use K0L mass(?)
1232       CLHEP::Hep3Vector pflow_eta_phi(_pflow_EM_x[em], _pflow_EM_y[em], _pflow_EM_z[em]);
1233 
1234       CLHEP::Hep3Vector pflow_vector;
1235       pflow_vector.setREtaPhi(residual_energy, pflow_eta_phi.eta(), pflow_eta_phi.phi());
1236 
1237       std::vector<RawCluster *> eclus;
1238       eclus.push_back(_pflow_EM_cluster.at(em));
1239 
1240       pflow->set_px(pflow_vector.getX());
1241       pflow->set_py(pflow_vector.getY());
1242       pflow->set_pz(pflow_vector.getZ());
1243       pflow->set_e(residual_energy);
1244       pflow->set_eclusters(eclus);
1245       pflow->set_hcluster(nullptr);
1246       pflow->set_track(nullptr);
1247       pflow->set_assoc_tracks(assoc_tracks);
1248       pflow->set_id(_global_pflow_index);
1249       pflow->set_type(ParticleFlowElement::PFLOWTYPE::LEFTOVER_EM_PARTICLE);
1250 /*
1251       std::cout << "Creating particle flow element" << std::endl;
1252       std::cout << "PF ID: " << _global_pflow_index << " as LEFTOVER_EM_PARTICLE" << std::endl;
1253       std::cout << "Track ID: " << "XXX" << std::endl;
1254       std::cout << "EMCal ID: " << em << std::endl;
1255       std::cout << "HCal ID: " << "XXX" << std::endl;
1256 */
1257       _pflowContainer->AddParticleFlowElement(_global_pflow_index, pflow);
1258       _global_pflow_index++;
1259     }
1260 
1261   }  // close EM loop
1262 }
1263 
1264 void ParticleFlowReco::ProcessClusterOnlyMatch()
1265 {
1266   CLHEP::Hep3Vector vertex = GetVertexForCalorimeter();
1267 
1268 
1269   for(unsigned int had = 0; had < _pflow_HAD_E.size(); had++)
1270   {
1271     // only consider HADs withOUT matched tracks ... we have dealt with the matched cases above
1272     if(_pflow_HAD_match_TRK.at(had).size() != 0)
1273     {
1274       continue;
1275     }
1276 
1277     //Don't look at clusters that were not matched with larger DeltaR
1278     if(_pflow_HAD_match_EM.at(had).size() == 0)
1279     {
1280       continue;
1281     }
1282 
1283     CLHEP::Hep3Vector hcal_cluster(_pflow_HAD_x[had], _pflow_HAD_y[had], _pflow_HAD_z[had]);
1284     hcal_cluster = hcal_cluster - vertex;
1285 
1286     //float matched_dR = 0.;
1287     int matched_index = -1;
1288     float min_dR = 999.;
1289 
1290     for(unsigned int em = 0; em < _pflow_EM_E.size(); em++)
1291     {
1292       // only consider EMs withOUT matched tracks ... we have dealt with the matched cases above
1293       if(_pflow_EM_match_TRK.at(em).size() != 0)
1294       {
1295         continue;
1296       }
1297 
1298       //Don't look at clusters that were not matched with larger DeltaR
1299       if(_pflow_EM_match_HAD.at(em).size() == 0)
1300       {
1301         continue;
1302       }
1303 
1304       CLHEP::Hep3Vector emcal_cluster(_pflow_EM_x[em], _pflow_EM_y[em], _pflow_EM_z[em]);
1305       emcal_cluster = emcal_cluster - vertex;
1306 
1307       float dR = fabs(emcal_cluster.deltaR(hcal_cluster));
1308 
1309       if(dR > 0.2)
1310       {
1311         continue;
1312       }
1313 
1314       if(dR < min_dR)
1315       {
1316         min_dR = dR;
1317         matched_index = em;
1318         //matched_dR = emcal_cluster.deltaR(hcal_cluster);
1319       }
1320 
1321     }
1322 
1323     if(matched_index < 0) //Previous HAD-EM match had also track matching or no match found
1324     {
1325       continue;
1326     }
1327 
1328     // add PFlow element for this EM+HAD
1329     ParticleFlowElement *pflow = new ParticleFlowElementv1();
1330 
1331     // assume massless, could be updated to use K0L
1332     // assuming EM kinematics due to higher granularity of the EMCal, this has to be further studied!
1333     CLHEP::Hep3Vector pflow_eta_phi(_pflow_EM_x[matched_index], _pflow_EM_y[matched_index], _pflow_EM_z[matched_index]);
1334     pflow_eta_phi = pflow_eta_phi - vertex;
1335 
1336     CLHEP::Hep3Vector pflow_vector;
1337     pflow_vector.setREtaPhi(_pflow_EM_E[matched_index]+_pflow_HAD_E[had], pflow_eta_phi.eta(), pflow_eta_phi.phi());
1338 
1339     std::vector<RawCluster *> eclus;
1340     eclus.push_back(_pflow_EM_cluster.at(matched_index));
1341 
1342     pflow->set_px(pflow_vector.getX());
1343     pflow->set_py(pflow_vector.getY());
1344     pflow->set_pz(pflow_vector.getZ());
1345     pflow->set_e(_pflow_EM_E[matched_index]);
1346     pflow->set_eclusters(eclus);
1347     pflow->set_hcluster(_pflow_HAD_cluster.at(had));
1348     pflow->set_track(nullptr);
1349     pflow->set_id(_global_pflow_index);
1350     pflow->set_type(ParticleFlowElement::PFLOWTYPE::UNMATCHED_NEUTRAL_HADRON);
1351 /*
1352     std::cout << "Creating particle flow element" << std::endl;
1353     std::cout << "PF ID: " << _global_pflow_index << " as UNMATCHED_NEUTRAL_HADRON" << std::endl;
1354     std::cout << "Track ID: " << "XXX" << std::endl;
1355     std::cout << "EMCal ID: " << matched_index << std::endl;
1356     std::cout << "HCal ID: " << had << std::endl;
1357 */
1358     _pflowContainer->AddParticleFlowElement(_global_pflow_index, pflow);
1359     _global_pflow_index++;
1360   }
1361 
1362   /*
1363 
1364         if(fabs(emcal_cluster.deltaR(hcal_cluster)) > 0.5)
1365         {
1366           std::cout << "No match" << std::endl;
1367           continue;
1368         }
1369 
1370         std::cout << "EMCal matched to HCal with dPhi: " << emcal_cluster.deltaPhi(hcal_cluster) << " dEta: " << emcal_cluster.eta()-hcal_cluster.eta() << std::endl;
1371         std::cout << "EMCal ID: " << em << " phi: " << emcal_cluster.phi() << " eta: " << emcal_cluster.eta() << std::endl;
1372         std::cout << "HCal ID: " << had << " phi: " << hcal_cluster.phi() << " eta: " << hcal_cluster.eta() << std::endl;
1373         for(unsigned int im = 0; im < _pflow_EM_match_HAD.at(em).size(); im++)
1374         {
1375           std::cout << "Previously Linked HAD ID = " << _pflow_EM_match_HAD.at(em).at(im) << std::endl;
1376         }
1377         */
1378 }
1379 
1380 void ParticleFlowReco::ProcessUnmatchedEMCalClusters()
1381 {
1382   for(unsigned int em = 0; em < _pflow_EM_E.size(); em++)
1383   {
1384     // only consider EMs withOUT matched tracks ... we have dealt with the matched cases above
1385     if(_pflow_EM_match_TRK.at(em).size() != 0)
1386     {
1387       continue;
1388     }
1389 
1390     // only consider EMs withOUT matched HADs
1391     if(_pflow_EM_match_HAD.at(em).size() != 0)
1392     {
1393       continue;
1394     }
1395 
1396     if(Verbosity() > 5)
1397     {
1398       std::cout << " unmatched EM " << em << " with E / x / y / z = " << _pflow_EM_E.at(em) << " / " << _pflow_EM_x.at(em) << " / " << _pflow_EM_y.at(em) << " / " << _pflow_EM_z.at(em) << std::endl;
1399     }
1400 
1401     // add PFlow element for this EM
1402     ParticleFlowElement *pflow = new ParticleFlowElementv1();
1403 
1404     // assume massless, could be updated to use K0L
1405     CLHEP::Hep3Vector vertex = GetVertexForCalorimeter();
1406     CLHEP::Hep3Vector pflow_eta_phi(_pflow_EM_x[em], _pflow_EM_y[em], _pflow_EM_z[em]);
1407     pflow_eta_phi = pflow_eta_phi - vertex;
1408 
1409     CLHEP::Hep3Vector pflow_vector;
1410     pflow_vector.setREtaPhi(_pflow_EM_E[em], pflow_eta_phi.eta(), pflow_eta_phi.phi());
1411 
1412     std::vector<RawCluster *> eclus;
1413     eclus.push_back(_pflow_EM_cluster.at(em));
1414 
1415     pflow->set_px(pflow_vector.getX());
1416     pflow->set_py(pflow_vector.getY());
1417     pflow->set_pz(pflow_vector.getZ());
1418     pflow->set_e(_pflow_EM_E[em]);
1419     pflow->set_eclusters(eclus);
1420     pflow->set_hcluster(nullptr);
1421     pflow->set_track(nullptr);
1422     pflow->set_id(_global_pflow_index);
1423     pflow->set_type(ParticleFlowElement::PFLOWTYPE::UNMATCHED_EM_PARTICLE);
1424 /*
1425     std::cout << "Creating particle flow element" << std::endl;
1426     std::cout << "PF ID: " << _global_pflow_index << " as UNMATCHED_EM_PARTICLE" << std::endl;
1427     std::cout << "Track ID: " << "XXX" << std::endl;
1428     std::cout << "EMCal ID: " << em << std::endl;
1429     std::cout << "HCal ID: " << "XXX" << std::endl;
1430 */
1431     _pflowContainer->AddParticleFlowElement(_global_pflow_index, pflow);
1432     _global_pflow_index++;
1433 
1434   }  // close EM loop
1435 }
1436 
1437 void ParticleFlowReco::ProcessUnmatchedHCalClusters()
1438 {
1439   for(unsigned int had = 0; had < _pflow_HAD_E.size(); had++)
1440   {
1441     // only consider HADs withOUT matched tracks ... we have dealt with the matched cases above
1442     if(_pflow_HAD_match_TRK.at(had).size() != 0)
1443     {
1444       continue;
1445     }
1446 
1447     // only consider HADs withOUT matched EMs
1448     if(_pflow_HAD_match_EM.at(had).size() != 0)
1449     {
1450       continue;
1451     }
1452 
1453     if(Verbosity() > 5)
1454     {
1455       std::cout << " unmatched HAD " << had << " with E / x / y / z = " << _pflow_HAD_E.at(had) << " / " << _pflow_HAD_x.at(had) << " / " << _pflow_HAD_y.at(had) << " / " << _pflow_HAD_z.at(had) << std::endl;
1456     }
1457 
1458     // add PFlow element for this HAD
1459     ParticleFlowElement *pflow = new ParticleFlowElementv1();
1460 
1461     // assume massless, could be updated to use K0L
1462     CLHEP::Hep3Vector vertex = GetVertexForCalorimeter();
1463     CLHEP::Hep3Vector pflow_eta_phi(_pflow_HAD_x[had], _pflow_HAD_y[had], _pflow_HAD_z[had]);
1464     pflow_eta_phi = pflow_eta_phi - vertex;
1465 
1466     CLHEP::Hep3Vector pflow_vector;
1467     pflow_vector.setREtaPhi(_pflow_HAD_E[had], pflow_eta_phi.eta(), pflow_eta_phi.phi());
1468 
1469     pflow->set_px(pflow_vector.getX());
1470     pflow->set_py(pflow_vector.getY());
1471     pflow->set_pz(pflow_vector.getZ());
1472     pflow->set_e(_pflow_HAD_E[had]);
1473     pflow->set_track(nullptr);
1474     pflow->set_eclusters(std::vector<RawCluster *>());
1475     pflow->set_hcluster(_pflow_HAD_cluster.at(had));
1476     pflow->set_id(_global_pflow_index);
1477     pflow->set_type(ParticleFlowElement::PFLOWTYPE::UNMATCHED_NEUTRAL_HADRON);
1478 /*
1479     std::cout << "Creating particle flow element" << std::endl;
1480     std::cout << "PF ID: " << _global_pflow_index << " as UNMATCHED_NEUTRAL_HADRON" << std::endl;
1481     std::cout << "Track ID: " << "XXX" << std::endl;
1482     std::cout << "EMCal ID: " << "XXX" << std::endl;
1483     std::cout << "HCal ID: " << had << std::endl;
1484 */
1485     _pflowContainer->AddParticleFlowElement(_global_pflow_index, pflow);
1486     _global_pflow_index++;
1487 
1488   }  // close HAD loop
1489 }
1490 
1491 void ParticleFlowReco::ProcessUnmatchedTracks()
1492 {
1493   for(unsigned int trk = 0; trk < _pflow_TRK_p.size(); trk++)
1494   {
1495     // only consider TRKs withOUT matched EM or HAD
1496     if(_pflow_TRK_match_EM.at(trk).size() != 0 || _pflow_TRK_match_HAD.at(trk).size() != 0)
1497     {
1498       continue;
1499     }
1500 
1501     if(Verbosity() > 5)
1502     {
1503       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;
1504     }
1505 
1506     // add PFlow element for this TRK
1507     ParticleFlowElement *pflow = new ParticleFlowElementv1();
1508 
1509     // assume massless, could be updated to use K0L
1510     float pfE = std::sqrt((_pflow_TRK_p[trk]*_pflow_TRK_p[trk]) + (0.135*0.135));
1511 
1512     CLHEP::Hep3Vector pflow_vector;
1513     pflow_vector.setREtaPhi(_pflow_TRK_p[trk], _pflow_TRK_eta[trk], _pflow_TRK_phi[trk]);
1514 
1515     pflow->set_px(pflow_vector.getX());
1516     pflow->set_py(pflow_vector.getY());
1517     pflow->set_pz(pflow_vector.getZ());
1518     pflow->set_e(pfE);
1519     pflow->set_track(_pflow_TRK_trk.at(trk));
1520     pflow->set_eclusters(std::vector<RawCluster *>());
1521     pflow->set_hcluster(nullptr);
1522     pflow->set_id(_global_pflow_index);
1523     pflow->set_type(ParticleFlowElement::PFLOWTYPE::UNMATCHED_CHARGED_HADRON);
1524 /*
1525     std::cout << "Creating particle flow element" << std::endl;
1526     std::cout << "PF ID: " << _global_pflow_index << " as UNMATCHED_CHARGED_HADRON" << std::endl;
1527     std::cout << "Track ID: " << trk << std::endl;
1528     std::cout << "EMCal ID: " << "XXX" << std::endl;
1529     std::cout << "HCal ID: " << "XXX" << std::endl;
1530 */
1531     _pflowContainer->AddParticleFlowElement(_global_pflow_index, pflow);
1532     _global_pflow_index++;
1533 
1534   }  // close TRK loop
1535 }
1536 
1537 CLHEP::Hep3Vector ParticleFlowReco::GetVertexForCalorimeter()
1538 {
1539   CLHEP::Hep3Vector calo_vertex(0., 0., 0.);
1540 
1541   if(_dont_use_global_vertex)
1542   {
1543     size_t max_ntracks_vertex_id = 0;
1544 
1545     for(unsigned int ivertex = 0; ivertex < _svtx_vertexmap->size(); ivertex++)
1546     {
1547       SvtxVertex *vertex = _svtx_vertexmap->get(ivertex);
1548       if(vertex->get_beam_crossing() != 0)
1549       {
1550         continue;
1551       }
1552 
1553       if(vertex->size_tracks() > max_ntracks_vertex_id)
1554       {
1555         max_ntracks_vertex_id = vertex->size_tracks();
1556         calo_vertex.setX(vertex->get_x());
1557         calo_vertex.setY(vertex->get_y());
1558         calo_vertex.setZ(vertex->get_z());
1559       }
1560     }
1561   }
1562   else
1563   {
1564     int max_ntracks_vertex_id = 0;
1565 
1566     for(auto &iter : *_global_vertexmap)
1567     {
1568       GlobalVertex *global_vtx = iter.second;
1569 
1570       if(global_vtx->get_beam_crossing() != 0)
1571       {
1572         continue;
1573       }
1574 
1575       for(GlobalVertex::VertexIter v_iter = global_vtx->begin_vertexes(); v_iter != global_vtx->end_vertexes(); ++v_iter)
1576       {
1577         auto vertex = *v_iter;
1578 
1579         //In case there is only MBD vertex, use it
1580         if(vertex.first == GlobalVertex::VTXTYPE::MBD)
1581         {
1582           calo_vertex.setX(global_vtx->get_x());
1583           calo_vertex.setY(global_vtx->get_y());
1584           calo_vertex.setZ(global_vtx->get_z());
1585         }
1586 
1587         for(auto v : vertex.second)
1588         {
1589           int n_tracks = 0;
1590           for(Vertex::TrackIter trk_iter = v->begin_tracks(); trk_iter != v->end_tracks(); ++trk_iter)
1591           {
1592             n_tracks++;
1593           }
1594           if(n_tracks > max_ntracks_vertex_id)
1595           {
1596             max_ntracks_vertex_id = n_tracks;
1597             calo_vertex.setX(global_vtx->get_x());
1598             calo_vertex.setY(global_vtx->get_y());
1599             calo_vertex.setZ(global_vtx->get_z());
1600           }
1601         }
1602 
1603       }
1604 
1605     }
1606 
1607   }
1608 
1609   return calo_vertex;
1610 }
1611 
1612 std::pair<int, int> ParticleFlowReco::GetCaloLimits(CLHEP::Hep3Vector &cluster, int trackID)
1613 {
1614   int inner_layer = -1;
1615   int outer_layer = -1;
1616 
1617   for(unsigned int icalo = 0; icalo < _pflow_TRK_projections.at(trackID).size(); icalo++)
1618   {
1619     if(std::isnan(_pflow_TRK_projections.at(trackID).at(icalo).x()))
1620     {
1621       break;
1622     }
1623 
1624     if(_pflow_TRK_projections.at(trackID).at(icalo).perp() < cluster.perp())
1625     {
1626       inner_layer = icalo;
1627     }
1628     else
1629     {
1630       outer_layer = icalo;
1631       break;
1632     }
1633   }
1634 
1635   return std::make_pair(inner_layer, outer_layer);
1636 }
1637 
1638 CLHEP::Hep3Vector ParticleFlowReco::Minimize_PhiEta(CLHEP::Hep3Vector &cluster, int trackID, CLHEP::Hep3Vector &vertex, bool isHCal)
1639 {
1640   CLHEP::Hep3Vector projection(NAN, NAN, NAN);
1641 
1642   std::pair<int, int> calo_limits = std::make_pair(0, 1);
1643 
1644   if(isHCal)
1645   {
1646     calo_limits = GetCaloLimits(cluster, trackID);
1647     if(calo_limits.first < 2)
1648     {
1649       return projection;
1650     }
1651   }
1652 
1653   if((calo_limits.first >= 0) && (calo_limits.second >= 0)) //has both states
1654   {
1655     TVectorD p1(2), p2(2), q(2);
1656 
1657     p1[0] = (_pflow_TRK_projections.at(trackID).at(calo_limits.first) - vertex).phi();
1658     p1[1] = (_pflow_TRK_projections.at(trackID).at(calo_limits.first) - vertex).eta();
1659 
1660     p2[0] = (_pflow_TRK_projections.at(trackID).at(calo_limits.second) - vertex).phi();
1661     p2[1] = (_pflow_TRK_projections.at(trackID).at(calo_limits.second) - vertex).eta();
1662 
1663     q[0] = cluster.phi();
1664     q[1] = cluster.eta();
1665 
1666     TVectorD d = p2; d -= p1;        // direction vector (p2 - p1)
1667     TVectorD r = p1; r -= q;         // difference (p1 - q)
1668 
1669     double a = d * d;                // dot(d, d)
1670     double b = 2.0 * (d * r);        // 2 * dot(d, p1 - q)
1671 
1672     double lambda = -b / (2.0 * a);  // unconstrained min
1673     lambda = std::max(0.0, std::min(1.0, lambda)); // clamp to [0, 1]
1674 
1675     projection = (1.-lambda)*_pflow_TRK_projections.at(trackID).at(calo_limits.first) + lambda*_pflow_TRK_projections.at(trackID).at(calo_limits.second);
1676   }
1677   else if((calo_limits.first >= 0) && (calo_limits.second < 0))
1678   {
1679     projection = _pflow_TRK_projections.at(trackID).at(calo_limits.first);
1680   }
1681 
1682   return projection;
1683 }
1684 
1685 bool ParticleFlowReco::isAcceptableTrack(SvtxTrack *track)
1686 {
1687   if(_only_crossing_zero && (track->get_crossing() != 0))
1688   {
1689     return false;
1690   }
1691 
1692   if(track->get_pt() < _min_track_pt)
1693   {
1694     return false;
1695   }
1696 
1697   if(fabs(track->get_eta()) > _acc_track_eta)
1698   {
1699     return false;
1700   }
1701 
1702   return true;
1703 }
1704 
1705 int ParticleFlowReco::CreateNode(PHCompositeNode *topNode)
1706 {
1707   PHNodeIterator iter(topNode);
1708 
1709   // Looking for the DST node
1710   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
1711   if(!dstNode)
1712   {
1713     std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
1714     return Fun4AllReturnCodes::ABORTRUN;
1715   }
1716 
1717   // store the PFlow elements under a sub-node directory
1718   PHCompositeNode *pflowNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "PARTICLEFLOW"));
1719   if(!pflowNode)
1720   {
1721     pflowNode = new PHCompositeNode("PARTICLEFLOW");
1722     dstNode->addNode(pflowNode);
1723   }
1724 
1725   // create the ParticleFlowElementContainer node...
1726   ParticleFlowElementContainer *pflowElementContainer = findNode::getClass<ParticleFlowElementContainer>(topNode, "ParticleFlowElements");
1727   if(!pflowElementContainer)
1728   {
1729     pflowElementContainer = new ParticleFlowElementContainer();
1730     PHIODataNode<PHObject> *pflowElementNode = new PHIODataNode<PHObject>(pflowElementContainer, "ParticleFlowElements", "PHObject");
1731     pflowNode->addNode(pflowElementNode);
1732   }
1733   else
1734   {
1735     std::cout << PHWHERE << "::ERROR - ParticleFlowElements node alerady exists, but should not" << std::endl;
1736     exit(-1);
1737   }
1738 
1739   return Fun4AllReturnCodes::EVENT_OK;
1740 }
[ 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!