ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ParticleFlowReco.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file ParticleFlowReco.cc
1 #include "ParticleFlowReco.h"
2 
5 
6 
7 #include <calobase/RawCluster.h>
8 #include <calobase/RawClusterContainer.h>
9 #include <calobase/RawTower.h>
10 #include <calobase/RawTowerContainer.h>
11 #include <calobase/RawTowerGeom.h>
12 #include <calobase/RawTowerGeomContainer.h>
13 
15 #include <g4main/PHG4Particle.h>
16 
18 
19 #include <phool/PHCompositeNode.h>
20 #include <phool/PHRandomSeed.h>
21 #include <phool/getClass.h>
22 
23 
24 #include <TLorentzVector.h>
25 
26 #include <gsl/gsl_randist.h>
27 #include <gsl/gsl_rng.h> // for gsl_rng_uniform_pos
28 
29 #include <iostream>
30 
31 // examine second value of std::pair, sort by smallest
32 bool sort_by_pair_second_lowest( const std::pair<int,float> &a, const std::pair<int,float> &b)
33 {
34  return (a.second < b.second);
35 }
36 
37 float ParticleFlowReco::calculate_dR( float eta1, float eta2, float phi1, float phi2 ) {
38 
39  float deta = eta1 - eta2;
40  float dphi = phi1 - phi2;
41  while ( dphi > 3.14159 ) dphi -= 2 * 3.14159;
42  while ( dphi < -3.14159 ) dphi += 2 * 3.14159;
43  return sqrt( pow( deta, 2 ) + pow( dphi ,2 ) );
44 
45 }
46 
47 std::pair<float, float> ParticleFlowReco::get_expected_signature( int trk ) {
48 
49  float response = ( 0.553437 + 0.0572246 * log( _pflow_TRK_p[ trk ] ) ) * _pflow_TRK_p[ trk ];
50  float resolution = sqrt( pow( 0.119123 , 2 ) + pow( 0.312361 , 2 ) / _pflow_TRK_p[ trk ] ) * _pflow_TRK_p[ trk ] ;
51 
52  std::pair<float, float> expected_signature( response , resolution );
53 
54  return expected_signature;
55 
56 }
57 
58 //____________________________________________________________________________..
60  SubsysReco(name),
61  _energy_match_Nsigma( 1.5 ) ,
62  _emulate_efficiency( 1.1 )
63 {
64 
65  _tr_eff = gsl_rng_alloc(gsl_rng_mt19937);
66  gsl_rng_set(_tr_eff,PHRandomSeed());
67 }
68 
69 //____________________________________________________________________________..
71 {
72  gsl_rng_free(_tr_eff);
73 }
74 
75 //____________________________________________________________________________..
77 {
79 }
80 
81 //____________________________________________________________________________..
83 {
84  return CreateNode(topNode);
85 
86 }
87 
88 //____________________________________________________________________________..
90 {
91  if (Verbosity() > 0)
92  {
93  std::cout << "ParticleFlowReco::process_event with Nsigma = " << _energy_match_Nsigma << " , emulate efficiency = " << _emulate_efficiency << std::endl;
94  }
95  // get handle to pflow node
96  ParticleFlowElementContainer *pflowContainer = findNode::getClass<ParticleFlowElementContainer>(topNode, "ParticleFlowElements");
97  if (!pflowContainer) {
98  std::cout << " ERROR -- can't find ParticleFlowElements node after it should have been created" << std::endl;
100  }
101 
102  // used for indexing PFlow elements in container
103  int global_pflow_index = 0;
104 
105  // read in towers
106  RawTowerContainer *towersEM = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_CEMC");
107  RawTowerContainer *towersIH = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALIN");
108  RawTowerContainer *towersOH = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALOUT");
109 
110  if ( !towersEM || !towersIH || !towersOH ) {
111  std::cout << "ParticleFlowReco::process_event : FATAL ERROR, cannot find tower containers" << std::endl;
113  }
114 
115  // read in tower geometries
116  RawTowerGeomContainer *geomEM = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC");
117  RawTowerGeomContainer *geomIH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
118  RawTowerGeomContainer *geomOH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
119 
120  if ( !geomEM || !geomIH || !geomOH ) {
121  std::cout << "ParticleFlowReco::process_event : FATAL ERROR, cannot find tower geometry containers" << std::endl;
123  }
124 
125  // read in clusters
126  RawClusterContainer *clustersEM = findNode::getClass<RawClusterContainer>(topNode, "TOPOCLUSTER_EMCAL");
127  RawClusterContainer *clustersHAD = findNode::getClass<RawClusterContainer>(topNode, "TOPOCLUSTER_HCAL");
128 
129  if ( !clustersEM ) {
130  std::cout << "ParticleFlowReco::process_event : FATAL ERROR, cannot find cluster container TOPOCLUSTER_EMCAL" << std::endl;
132  }
133  if ( !clustersHAD ) {
134  std::cout << "ParticleFlowReco::process_event : FATAL ERROR, cannot find cluster container TOPOCLUSTER_HCAL" << std::endl;
136  }
137 
138  // reset internal particle-flow representation
139  _pflow_TRK_p.clear();
140  _pflow_TRK_eta.clear();
141  _pflow_TRK_phi.clear();
142  _pflow_TRK_match_EM.clear();
143  _pflow_TRK_match_HAD.clear();
145 
146  _pflow_EM_E.clear();
147  _pflow_EM_eta.clear();
148  _pflow_EM_phi.clear();
149  _pflow_EM_tower_eta.clear();
150  _pflow_EM_tower_phi.clear();
151  _pflow_EM_match_HAD.clear();
152  _pflow_EM_match_TRK.clear();
153 
154  _pflow_HAD_E.clear();
155  _pflow_HAD_eta.clear();
156  _pflow_HAD_phi.clear();
157  _pflow_HAD_tower_eta.clear();
158  _pflow_HAD_tower_phi.clear();
159  _pflow_HAD_match_EM.clear();
160  _pflow_HAD_match_TRK.clear();
161 
162 
163  if ( Verbosity() > 2 )
164  std::cout << "ParticleFlowReco::process_event : initial population of TRK, EM, HAD objects " << std::endl;
165 
166  // read in tracks with > 0.5 GeV
167  // currently just do this from truth
168  {
169 
170  PHG4TruthInfoContainer* truthinfo = findNode::getClass <PHG4TruthInfoContainer> (topNode, "G4TruthInfo");
172 
173  for (PHG4TruthInfoContainer::ConstIterator iter = range.first; iter != range.second; ++iter) {
174  PHG4Particle* g4particle = iter->second;
175 
176  TLorentzVector t;
177  t.SetPxPyPzE( g4particle->get_px(), g4particle->get_py(), g4particle->get_pz(), g4particle->get_e() );
178 
179  float truth_pt = t.Pt();
180  float truth_p = t.P();
181  if (truth_pt < 0.5)
182  continue; // only keep pt > 0.5 GeV
183 
184  float truth_eta = t.Eta();
185  if (fabs (truth_eta) > 1.1)
186  continue; // only keep |eta| < 1.1
187 
188  float truth_phi = t.Phi();
189 
190  int truth_pid = abs( g4particle->get_pid() ); // particle species
191 
192  // only keep charged truth particles
193  if ( truth_pid != 211 && truth_pid != 321 && truth_pid != 2212 ) continue;
194 
195  // if emulating finite efficiency, roll the dice here
196  if ( _emulate_efficiency < 1.0 ) {
197 
198  float this_eff = gsl_rng_uniform_pos(_tr_eff);
199 
200  if ( this_eff > _emulate_efficiency ) {
201  // lost this track due to inefficiency
202  continue;
203  }
204 
205  }
206 
207  _pflow_TRK_p.push_back( truth_p );
208  _pflow_TRK_eta.push_back( truth_eta );
209  _pflow_TRK_phi.push_back( truth_phi );
210 
211  _pflow_TRK_match_EM.push_back( std::vector<int>() );
212  _pflow_TRK_match_HAD.push_back( std::vector<int>() );
213 
214  _pflow_TRK_addtl_match_EM.push_back( std::vector< std::pair<int,float> >() );
215 
216  if ( Verbosity() > 5 && truth_pt > 0.5 )
217  std::cout << " TRK with p / pT = " << truth_p << " / " << truth_pt << " , eta / phi = " << truth_eta << " / " << truth_phi << std::endl;
218 
219  } // close truth paticle loop
220 
221  } //
222 
223 
224 
225  // read in EMCal topoClusters with E > 0.2 GeV
226  {
227  RawClusterContainer::ConstRange begin_end = clustersEM->getClusters();
228  for ( RawClusterContainer::ConstIterator hiter = begin_end.first; hiter != begin_end.second; ++hiter)
229  {
230  float cluster_E = hiter->second->get_energy();
231  if ( cluster_E < 0.2 ) continue;
232 
233  float cluster_phi = hiter->second->get_phi();
234  // for now, assume event at vx_z = 0
235  float cluster_theta = 3.14159 / 2.0 - atan2( hiter->second->get_z() , hiter->second->get_r() );
236  float cluster_eta = -1 * log( tan( cluster_theta / 2.0 ) );
237 
238  _pflow_EM_E.push_back( cluster_E );
239  _pflow_EM_eta.push_back( cluster_eta );
240  _pflow_EM_phi.push_back( cluster_phi );
241 
242  _pflow_EM_match_HAD.push_back( std::vector<int>() );
243  _pflow_EM_match_TRK.push_back( std::vector<int>() );
244 
245  if ( Verbosity() > 5 && cluster_E > 0.2 )
246  std::cout << " EM topoCluster with E = " << cluster_E << ", eta / phi = " << cluster_eta << " / " << cluster_phi << " , nTow = " << hiter->second->getNTowers() << std::endl;
247 
248  std::vector<float> this_cluster_tower_eta;
249  std::vector<float> this_cluster_tower_phi;
250 
251  // read in towers
252  RawCluster::TowerConstRange begin_end_towers = hiter->second->get_towers();
253  for (RawCluster::TowerConstIterator iter = begin_end_towers.first; iter != begin_end_towers.second; ++iter) {
254 
256  RawTower* tower = towersEM->getTower(iter->first);
257  RawTowerGeom *tower_geom = geomEM->get_tower_geometry(tower->get_key());
258 
259  this_cluster_tower_phi.push_back( tower_geom->get_phi() );
260  this_cluster_tower_eta.push_back( tower_geom->get_eta() );
261  }
262  else {
263  std::cout << "ParticleFlowReco::process_event : FATAL ERROR , EM topoClusters seem to contain HCal towers" << std::endl;
265  }
266  } // close tower loop
267 
268  _pflow_EM_tower_eta.push_back( this_cluster_tower_eta );
269  _pflow_EM_tower_phi.push_back( this_cluster_tower_phi );
270 
271  } // close cluster loop
272 
273  } // close
274 
275  // read in HCal topoClusters with E > 0.2 GeV
276  {
277  RawClusterContainer::ConstRange begin_end = clustersHAD->getClusters();
278  for ( RawClusterContainer::ConstIterator hiter = begin_end.first; hiter != begin_end.second; ++hiter)
279  {
280  float cluster_E = hiter->second->get_energy();
281  if ( cluster_E < 0.2 ) continue;
282 
283  float cluster_phi = hiter->second->get_phi();
284  // for now, assume event at vx_z = 0
285  float cluster_theta = 3.14159 / 2.0 - atan2( hiter->second->get_z() , hiter->second->get_r() );
286  float cluster_eta = -1 * log( tan( cluster_theta / 2.0 ) );
287 
288  _pflow_HAD_E.push_back( cluster_E );
289  _pflow_HAD_eta.push_back( cluster_eta );
290  _pflow_HAD_phi.push_back( cluster_phi );
291 
292  _pflow_HAD_match_EM.push_back( std::vector<int>() );
293  _pflow_HAD_match_TRK.push_back( std::vector<int>() );
294 
295  if ( Verbosity() > 5 && cluster_E > 0.2 )
296  std::cout << " HAD topoCluster with E = " << cluster_E << ", eta / phi = " << cluster_eta << " / " << cluster_phi << " , nTow = " << hiter->second->getNTowers() << std::endl;
297 
298  std::vector<float> this_cluster_tower_eta;
299  std::vector<float> this_cluster_tower_phi;
300 
301  // read in towers
302  RawCluster::TowerConstRange begin_end_towers = hiter->second->get_towers();
303  for (RawCluster::TowerConstIterator iter = begin_end_towers.first; iter != begin_end_towers.second; ++iter) {
304 
306 
307  RawTower* tower = towersIH->getTower(iter->first);
308  RawTowerGeom *tower_geom = geomIH->get_tower_geometry(tower->get_key());
309 
310  this_cluster_tower_phi.push_back( tower_geom->get_phi() );
311  this_cluster_tower_eta.push_back( tower_geom->get_eta() );
312  }
313 
315 
316  RawTower* tower = towersOH->getTower(iter->first);
317  RawTowerGeom *tower_geom = geomOH->get_tower_geometry(tower->get_key());
318 
319  this_cluster_tower_phi.push_back( tower_geom->get_phi() );
320  this_cluster_tower_eta.push_back( tower_geom->get_eta() );
321  } else {
322  std::cout << "ParticleFlowReco::process_event : FATAL ERROR , HCal topoClusters seem to contain EM towers" << std::endl;
324  }
325 
326 
327  } // close tower loop
328 
329  _pflow_HAD_tower_eta.push_back( this_cluster_tower_eta );
330  _pflow_HAD_tower_phi.push_back( this_cluster_tower_phi );
331 
332  } // close cluster loop
333 
334  } // close
335 
336  // BEGIN LINKING STEP
337 
338  // Link TRK -> EM (best match, but keep reserve of others), and TRK -> HAD (best match)
339  if ( Verbosity() > 2 )
340  std::cout << "ParticleFlowReco::process_event : TRK -> EM and TRK -> HAD linking " << std::endl;
341 
342  for (unsigned int trk = 0; trk < _pflow_TRK_p.size() ; trk++ ) {
343 
344  if ( Verbosity() > 10 )
345  std::cout << " TRK " << trk << " with p / eta / phi = " << _pflow_TRK_p[ trk ] << " / " << _pflow_TRK_eta[ trk ] << " / " << _pflow_TRK_phi[ trk ] << std::endl;
346 
347  // TRK -> EM link
348  float min_em_dR = 0.2;
349  int min_em_index = -1;
350 
351  for (unsigned int em = 0 ; em < _pflow_EM_E.size() ; em++) {
352 
353  float dR = calculate_dR( _pflow_TRK_eta[ trk ] , _pflow_EM_eta[ em ] , _pflow_TRK_phi[ trk ] , _pflow_EM_phi[ em ] );
354  if ( dR > 0.2 ) continue;
355 
356  bool has_overlap = false;
357 
358  for (unsigned int tow = 0; tow < _pflow_EM_tower_eta.at( em ).size() ; tow++) {
359 
360  float tower_eta = _pflow_EM_tower_eta.at( em ).at( tow );
361  float tower_phi = _pflow_EM_tower_phi.at( em ).at( tow );
362 
363  float deta = tower_eta - _pflow_TRK_eta[ trk ];
364  float dphi = tower_phi - _pflow_TRK_phi[ trk ];
365  if ( dphi > 3.14159 ) dphi -= 2 * 3.14159;
366  if ( dphi < -3.14159 ) dphi += 2 * 3.14159;
367 
368  if ( fabs( deta ) < 0.025 * 2.5 && fabs( dphi ) < 0.025 * 2.5 ) {
369  has_overlap = true;
370  break;
371  }
372 
373  }
374 
375  if ( has_overlap ) {
376 
377  if ( Verbosity() > 5 )
378  std::cout << " -> possible match to EM " << em << " with dR = " << dR << std::endl;
379 
380  _pflow_TRK_addtl_match_EM.at( trk ).push_back( std::pair<int,float>( em, dR ) );
381 
382  } else {
383 
384  if ( Verbosity() > 5 )
385  std::cout << " -> no match to EM " << em << " (even though dR = " << dR << " )" << std::endl;
386 
387  }
388 
389  }
390 
391  // sort possible matches
392 
393  std::sort( _pflow_TRK_addtl_match_EM.at( trk ).begin(), _pflow_TRK_addtl_match_EM.at( trk ).end(), sort_by_pair_second_lowest );
394  if ( Verbosity() > 10 ) {
395  for (unsigned int n = 0; n < _pflow_TRK_addtl_match_EM.at( trk ).size(); n++) {
396  std::cout << " -> sorted list of matches, EM / dR = " << _pflow_TRK_addtl_match_EM.at( trk ).at( n ).first << " / " << _pflow_TRK_addtl_match_EM.at( trk ).at( n ).second << std::endl;
397  }
398  }
399 
400  if ( _pflow_TRK_addtl_match_EM.at( trk ).size() > 0 ) {
401  min_em_index = _pflow_TRK_addtl_match_EM.at( trk ).at( 0 ).first;
402  min_em_dR = _pflow_TRK_addtl_match_EM.at( trk ).at( 0 ).second;
403  // delete best matched element
404  _pflow_TRK_addtl_match_EM.at( trk ).erase( _pflow_TRK_addtl_match_EM.at( trk ).begin() );
405  }
406 
407 
408  if ( min_em_index > -1 ) {
409  _pflow_EM_match_TRK.at( min_em_index ).push_back( trk );
410  _pflow_TRK_match_EM.at( trk ).push_back( min_em_index );
411 
412  if ( Verbosity() > 5 ) {
413  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;
414  std::cout << " ( " << _pflow_TRK_addtl_match_EM.at( trk ).size() << " other possible matches ) " << std::endl;
415  }
416 
417  } else {
418 
419  if ( Verbosity() > 5 )
420  std::cout << " -> no EM match! ( best dR = " << min_em_dR << " ) " << std::endl;
421  }
422 
423  // TRK -> HAD link
424  float min_had_dR = 0.2;
425  int min_had_index = -1;
426  float max_had_pt = 0;
427 
428  // TODO: sequential linking should better happen here -- i.e. allow EM-matched HAD's into the possible pool
429  for (unsigned int had = 0 ; had < _pflow_HAD_E.size() ; had++) {
430 
431  float dR = calculate_dR( _pflow_TRK_eta[ trk ] , _pflow_HAD_eta[ had ] , _pflow_TRK_phi[ trk ] , _pflow_HAD_phi[ had ] );
432  if ( dR > 0.5 ) continue;
433 
434  bool has_overlap = false;
435 
436  for (unsigned int tow = 0; tow < _pflow_HAD_tower_eta.at( had ).size() ; tow++) {
437 
438  float tower_eta = _pflow_HAD_tower_eta.at( had ).at( tow );
439  float tower_phi = _pflow_HAD_tower_phi.at( had ).at( tow );
440 
441  float deta = tower_eta - _pflow_TRK_eta[ trk ];
442  float dphi = tower_phi - _pflow_TRK_phi[ trk ];
443  if ( dphi > 3.14159 ) dphi -= 2 * 3.14159;
444  if ( dphi < -3.14159 ) dphi += 2 * 3.14159;
445 
446  if ( fabs( deta ) < 0.1 * 1.5 && fabs( dphi ) < 0.1 * 1.5 ) {
447  has_overlap = true;
448  break;
449  }
450 
451  }
452 
453  if ( has_overlap ) {
454 
455  if ( Verbosity() > 5 )
456  std::cout << " -> possible match to HAD " << had << " with dR = " << dR << std::endl;
457 
458  if ( _pflow_HAD_E.at( had ) > max_had_pt ) {
459  max_had_pt = _pflow_HAD_E.at( had );
460  min_had_index = had;
461  min_had_dR = dR;
462  }
463 
464  } else {
465 
466  if ( Verbosity() > 5 )
467  std::cout << " -> no match to HAD " << had << " (even though dR = " << dR << " )" << std::endl;
468 
469  }
470 
471  }
472 
473  if ( min_had_index > -1 ) {
474  _pflow_HAD_match_TRK.at( min_had_index ).push_back( trk );
475  _pflow_TRK_match_HAD.at( trk ).push_back( min_had_index );
476 
477  if ( Verbosity() > 5 )
478  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;
479 
480  } else {
481  if ( Verbosity() > 5 )
482  std::cout << " -> no HAD match! ( best dR = " << min_had_dR << " ) " << std::endl;
483  }
484 
485 
486  }
487 
488 
489  // EM->HAD linking
490  if ( Verbosity() > 2 )
491  std::cout << "ParticleFlowReco::process_event : EM -> HAD linking " << std::endl;
492 
493  for (unsigned int em = 0; em < _pflow_EM_E.size() ; em++ ) {
494 
495  if ( Verbosity() > 10 )
496  std::cout << " EM with E / eta / phi = " << _pflow_EM_E[ em ] << " / " << _pflow_EM_eta[ em ] << " / " << _pflow_EM_phi[ em ] << std::endl;
497 
498  // TRK -> HAD link
499  float min_had_dR = 0.2;
500  int min_had_index = -1;
501  float max_had_pt = 0;
502 
503  for (unsigned int had = 0 ; had < _pflow_HAD_E.size() ; had++) {
504 
505  float dR = calculate_dR( _pflow_EM_eta[ em ] , _pflow_HAD_eta[ had ] , _pflow_EM_phi[ em ] , _pflow_HAD_phi[ had ] );
506  if ( dR > 0.5 ) continue;
507 
508  bool has_overlap = false;
509 
510  for (unsigned int tow = 0; tow < _pflow_HAD_tower_eta.at( had ).size() ; tow++) {
511 
512  float tower_eta = _pflow_HAD_tower_eta.at( had ).at( tow );
513  float tower_phi = _pflow_HAD_tower_phi.at( had ).at( tow );
514 
515  float deta = tower_eta - _pflow_EM_eta[ em ];
516  float dphi = tower_phi - _pflow_EM_phi[ em ];
517  if ( dphi > 3.14159 ) dphi -= 2 * 3.14159;
518  if ( dphi < -3.14159 ) dphi += 2 * 3.14159;
519 
520  if ( fabs( deta ) < 0.1 * 1.5 && fabs( dphi ) < 0.1 * 1.5 ) {
521  has_overlap = true;
522  break;
523  }
524 
525  }
526 
527  if ( has_overlap ) {
528 
529  if ( Verbosity() > 5 )
530  std::cout << " -> possible match to HAD " << had << " with dR = " << dR << std::endl;
531 
532  if ( _pflow_HAD_E.at( had ) > max_had_pt ) {
533  max_had_pt = _pflow_HAD_E.at( had );
534  min_had_index = had;
535  min_had_dR = dR;
536  }
537 
538  } else {
539 
540  if ( Verbosity() > 5 )
541  std::cout << " -> no match to HAD " << had << " (even though dR = " << dR << " )" << std::endl;
542 
543  }
544 
545  }
546 
547  if ( min_had_index > -1 ) {
548  _pflow_HAD_match_EM.at( min_had_index ).push_back( em );
549  _pflow_EM_match_HAD.at( em ).push_back( min_had_index );
550 
551  if ( Verbosity() > 5 )
552  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;
553 
554  } else {
555  if ( Verbosity() > 5 )
556  std::cout << " -> no HAD match! ( best dR = " << min_had_dR << " ) " << std::endl;
557  }
558 
559  }
560 
561 
562  // SEQUENTIAL MATCHING: if TRK -> EM and EM -> HAD, ensure that TRK -> HAD
563  if ( Verbosity() > 2 )
564  std::cout << "ParticleFlowReco::process_event : sequential TRK -> EM && EM -> HAD ==> TRK -> HAD matching " << std::endl;
565 
566  for (unsigned int trk = 0; trk < _pflow_TRK_p.size() ; trk++ ) {
567 
568  // go through all matched EMs
569  for (unsigned int i = 0; i < _pflow_TRK_match_EM.at( trk ).size(); i++) {
570 
571  int em = _pflow_TRK_match_EM.at( trk ).at( i );
572 
573  // if this EM has a matched HAD...
574  for (unsigned int j = 0; j < _pflow_EM_match_HAD.at( em ).size() ; j++) {
575 
576  int had = _pflow_EM_match_HAD.at( em ).at( j );
577 
578  // and the TRK is NOT matched to this HAD...
579  bool is_trk_matched_to_HAD = false;
580  for (unsigned int k = 0; k < _pflow_TRK_match_HAD.at( trk ).size(); k++) {
581  int existing_had = _pflow_TRK_match_HAD.at( trk ).at( k );
582  if ( had == existing_had ) is_trk_matched_to_HAD = true;
583  }
584 
585  // if this is the case, create TRK->HAD link
586  if ( ! is_trk_matched_to_HAD ) {
587  _pflow_TRK_match_HAD.at( trk ).push_back( had );
588  _pflow_HAD_match_TRK.at( had ).push_back( trk );
589 
590  if ( Verbosity() > 5 ) {
591  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;
592  std::cout << " -> sequential match to HAD " << had << " through EM " << j << std::endl;
593  }
594 
595  }
596 
597  } // close the HAD loop
598 
599  } // close the EM loop
600 
601  } // close the TRK loop
602 
603  // TRK->EM->HAD removal
604  if ( Verbosity() > 2 )
605  std::cout << "ParticleFlowReco::process_event : resolve TRK(s) + EM(s) -> HAD systems " << std::endl;
606 
607  for (unsigned int had = 0; had < _pflow_HAD_E.size() ; had++ ) {
608 
609  // only consider HAD with matched tracks ... others we will deal with later
610  if ( _pflow_HAD_match_TRK.at( had ).size() == 0 ) continue;
611 
612  if ( Verbosity() > 5 ) {
613  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;
614  }
615 
616  // setup for Sum-pT^trk -> calo prediction
617  float total_TRK_p = 0;
618  float total_expected_E = 0;
619  float total_expected_E_var = 0;
620 
621  // begin with this HAD calo energy
622  float total_EMHAD_E = _pflow_HAD_E.at( had );
623 
624  // iterate over the EMs matched to this HAD
625  for (unsigned int j = 0; j < _pflow_HAD_match_EM.at( had ).size() ; j++ ) {
626 
627  int em = _pflow_HAD_match_EM.at( had ).at( j );
628 
629  // ensure there is at least one track matched to this EM
630  if ( _pflow_EM_match_TRK.at( em ).size() == 0 ) continue;
631 
632  // add it to the total calo E
633  total_EMHAD_E += _pflow_EM_E.at( em );
634 
635  if ( Verbosity() > 5 ) {
636  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;
637  }
638 
639  }
640 
641  // iterate over the TRKs matched to this HAD
642  for (unsigned int j = 0 ; j < _pflow_HAD_match_TRK.at( had ).size() ; j++) {
643 
644  int trk = _pflow_HAD_match_TRK.at( had ).at( j );
645 
646  if ( Verbosity() > 5 ) {
647  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;
648  }
649 
650  total_TRK_p += _pflow_TRK_p.at( trk );
651 
652  std::pair<float, float> expected_signature = get_expected_signature( trk );
653 
654  float expected_E_mean = expected_signature.first;
655  float expected_E_sigma = expected_signature.second;
656 
657  if ( Verbosity() > 5 ) {
658  std::cout << " -> -> -> expected calo signature is " << expected_E_mean << " +/- " << expected_E_sigma << std::endl;
659  }
660 
661  total_expected_E += expected_E_mean;
662  total_expected_E_var += pow( expected_E_sigma , 2 );
663 
664  // add PFlow element for each track
666 
667  // assume pion mass
668  TLorentzVector tlv; tlv.SetPtEtaPhiM( _pflow_TRK_p[ trk ] / cosh( _pflow_TRK_eta[ trk ] ) , _pflow_TRK_eta[ trk ] , _pflow_TRK_phi[ trk ] , 0.135 );
669 
670  pflow->set_px( tlv.Px() );
671  pflow->set_py( tlv.Py() );
672  pflow->set_pz( tlv.Pz() );
673  pflow->set_e( tlv.E() );
674  pflow->set_id( global_pflow_index );
675  pflow->set_type( ParticleFlowElement::PFLOWTYPE::MATCHED_CHARGED_HADRON );
676 
677  pflowContainer->AddParticleFlowElement( global_pflow_index, pflow );
678  global_pflow_index++;
679 
680  }
681 
682 
683  // process compatibility of fit
684  float total_expected_E_err = sqrt( total_expected_E_var );
685 
686  if ( Verbosity() > 5 ) {
687  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;
688  }
689 
690  // if Sum pT > calo, add in additional possible matched EMs associated with tracks until that is no longer the case
691 
692  if ( total_expected_E > total_EMHAD_E ) {
693 
694  if ( Verbosity() > 5 )
695  std::cout << " -> Expected E > Observed E, looking for additional potential TRK->EM matches" << std::endl;
696 
697  std::map<int, float> additional_EMs;
698 
699  for (unsigned int j = 0 ; j < _pflow_HAD_match_TRK.at( had ).size() ; j++) {
700 
701  int trk = _pflow_HAD_match_TRK.at( had ).at( j );
702 
703  int addtl_matches = _pflow_TRK_addtl_match_EM.at( trk ).size();
704 
705  if ( Verbosity() > 10 )
706  std::cout << " -> -> TRK " << trk << " has " << addtl_matches << " additional matches! " << std::endl;
707 
708  for (unsigned int n = 0 ; n < _pflow_TRK_addtl_match_EM.at( trk ).size() ; n++ ) {
709  if ( Verbosity() > 10 )
710  std::cout << " -> -> -> additional match to EM = " << _pflow_TRK_addtl_match_EM.at( trk ).at( n ).first << " with dR = " << _pflow_TRK_addtl_match_EM.at( trk ).at( n ).second << std::endl;
711 
712  float existing_dR = 0.21;
713  int counts = additional_EMs.count( _pflow_TRK_addtl_match_EM.at( trk ).at( n ).first );
714  if ( counts > 0 ) {
715  existing_dR = additional_EMs[ _pflow_TRK_addtl_match_EM.at( trk ).at( n ).first ];
716  }
717  if ( _pflow_TRK_addtl_match_EM.at( trk ).at( n ).second < existing_dR )
718  additional_EMs[ _pflow_TRK_addtl_match_EM.at( trk ).at( n ).first ] = _pflow_TRK_addtl_match_EM.at( trk ).at( n ).second;
719  }
720 
721  }
722 
723  // map now assured to have only minimal dR values for each possible additional EM
724  // translate the map to a vector of pairs, then sort by smallest dR
725 
726  std::vector< std::pair<int,float> > additional_EMs_vec;
727 
728  for (auto& x : additional_EMs) {
729  additional_EMs_vec.push_back( std::pair<int,float>( x.first , x.second ) );
730  }
731 
732  std::sort( additional_EMs_vec.begin(), additional_EMs_vec.end(), sort_by_pair_second_lowest );
733 
734  if ( Verbosity() > 5 )
735  std::cout << " -> Sorting the set of potential additional EMs " << std::endl;
736 
737  // now add in additional EMs until there are none left or it is no longer the case that Sum pT > calo
738 
739  int n_EM_added = 0;
740  while ( additional_EMs_vec.size() != 0 && total_expected_E > total_EMHAD_E ) {
741 
742  int new_EM = additional_EMs_vec.at( 0 ).first;
743 
744  if ( Verbosity() > 5 )
745  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;
746 
747  // for now, just make the first HAD-linked track point to this new EM, and vice versa
748  _pflow_EM_match_TRK.at( new_EM ).push_back( _pflow_HAD_match_TRK.at( had ).at( 0 ) );
749  _pflow_TRK_match_EM.at( _pflow_HAD_match_TRK.at( had ).at( 0 ) ).push_back( new_EM );
750 
751  // add to expected calo
752  total_EMHAD_E += _pflow_EM_E.at( new_EM );
753 
754  // erase lowest-dR EM
755  additional_EMs_vec.erase( additional_EMs_vec.begin() );
756 
757  n_EM_added++;
758  }
759 
760  if ( Verbosity() > 5) {
761  if ( n_EM_added > 0 ) {
762  std::cout << "After adding N = " << n_EM_added << " any additional EMs : " << std::endl;
763  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;
764  }
765  else {
766  std::cout << "No additional EMs found, continuing hypothesis check" << std::endl;
767  }
768  }
769  }
770 
771 
772 
773  if ( total_expected_E + _energy_match_Nsigma * total_expected_E_err > total_EMHAD_E ) {
774 
775  if ( Verbosity() > 5 ) {
776  std::cout << " -> -> calo compatible within Nsigma = " << _energy_match_Nsigma << " , remove and keep tracks " << std::endl;
777  }
778 
779  // PFlow elements already created from tracks above, no more needs to be done
780 
781  } else {
782 
783  float residual_energy = total_EMHAD_E - total_expected_E;
784 
785  if ( Verbosity() > 5 ) {
786  std::cout << " -> -> calo not compatible, create leftover cluster with " << residual_energy << std::endl;
787  }
788 
789  // create additional PFlow element (tracks already created above)
791 
792  // assume no mass, but could update to use K0L mass(?)
793  TLorentzVector tlv; tlv.SetPtEtaPhiM( residual_energy / cosh( _pflow_HAD_eta[ had ] ) , _pflow_HAD_eta[ had ] , _pflow_HAD_phi[ had ] , 0 );
794 
795  pflow->set_px( tlv.Px() );
796  pflow->set_py( tlv.Py() );
797  pflow->set_pz( tlv.Pz() );
798  pflow->set_e( tlv.E() );
799  pflow->set_id( global_pflow_index );
800  pflow->set_type( ParticleFlowElement::PFLOWTYPE::LEFTOVER_EM_PARTICLE );
801 
802  pflowContainer->AddParticleFlowElement( global_pflow_index, pflow );
803  global_pflow_index++;
804 
805  }
806 
807  } // close HAD loop
808 
809  // TRK->EM removal
810 
811  if ( Verbosity() > 2 )
812  std::cout << "ParticleFlowReco::process_event : resolve TRK(s) -> EM(s) ( + no HAD) systems " << std::endl;
813 
814  for (unsigned int em = 0; em < _pflow_EM_E.size() ; em++ ) {
815 
816  // only consider EM with matched tracks, but no matched HADs
817  if ( _pflow_EM_match_HAD.at( em ).size() != 0 ) continue;
818  if ( _pflow_EM_match_TRK.at( em ).size() == 0 ) continue;
819 
820  if ( Verbosity() > 5 ) {
821  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;
822  }
823 
824  // setup for Sum-pT^trk -> calo prediction
825  float total_TRK_p = 0;
826  float total_expected_E = 0;
827  float total_expected_E_var = 0;
828 
829  // begin with this EM calo energy
830  float total_EM_E = _pflow_EM_E.at( em );
831 
832  // iterate over the TRKs matched to this EM
833  for (unsigned int j = 0 ; j < _pflow_EM_match_TRK.at( em ).size() ; j++) {
834 
835  int trk = _pflow_EM_match_TRK.at( em ).at( j );
836 
837  if ( Verbosity() > 5 ) {
838  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;
839  }
840 
841  total_TRK_p += _pflow_TRK_p.at( trk );
842 
843  std::pair<float, float> expected_signature = get_expected_signature( trk );
844 
845  float expected_E_mean = expected_signature.first;
846  float expected_E_sigma = expected_signature.second;
847 
848  if ( Verbosity() > 5 ) {
849  std::cout << " -> -> -> expected calo signature is " << expected_E_mean << " +/- " << expected_E_sigma << std::endl;
850  }
851 
852  total_expected_E += expected_E_mean;
853  total_expected_E_var += pow( expected_E_sigma , 2 );
854 
855  // add PFlow element for each track
857 
858  // assume pion mass
859  TLorentzVector tlv; tlv.SetPtEtaPhiM( _pflow_TRK_p[ trk ] / cosh( _pflow_TRK_eta[ trk ] ) , _pflow_TRK_eta[ trk ] , _pflow_TRK_phi[ trk ] , 0.135 );
860 
861  pflow->set_px( tlv.Px() );
862  pflow->set_py( tlv.Py() );
863  pflow->set_pz( tlv.Pz() );
864  pflow->set_e( tlv.E() );
865  pflow->set_id( global_pflow_index );
866  pflow->set_type( ParticleFlowElement::PFLOWTYPE::MATCHED_CHARGED_HADRON );
867 
868  pflowContainer->AddParticleFlowElement( global_pflow_index, pflow );
869  global_pflow_index++;
870 
871  }
872 
873  // process compatibility of fit
874  float total_expected_E_err = sqrt( total_expected_E_var );
875 
876  if ( Verbosity() > 5 ) {
877  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;
878  }
879 
880  if ( total_expected_E + _energy_match_Nsigma * total_expected_E_err > total_EM_E ) {
881 
882  if ( Verbosity() > 5 ) {
883  std::cout << " -> -> calo compatible within Nsigma = " << _energy_match_Nsigma << " , remove and keep tracks " << std::endl;
884  }
885 
886  // PFlow elements already created from tracks above, no more needs to be done
887 
888  } else {
889 
890  float residual_energy = total_EM_E - total_expected_E;
891 
892  if ( Verbosity() > 5 ) {
893  std::cout << " -> -> calo not compatible, create leftover cluster with " << residual_energy << std::endl;
894  }
895 
896  // create additional PFlow element (tracks already created above)
898 
899  // assume no mass, but could update to use K0L mass(?)
900  TLorentzVector tlv; tlv.SetPtEtaPhiM( residual_energy / cosh( _pflow_EM_eta[ em ] ) , _pflow_EM_eta[ em ] , _pflow_EM_phi[ em ] , 0 );
901 
902  pflow->set_px( tlv.Px() );
903  pflow->set_py( tlv.Py() );
904  pflow->set_pz( tlv.Pz() );
905  pflow->set_e( tlv.E() );
906  pflow->set_id( global_pflow_index );
907  pflow->set_type( ParticleFlowElement::PFLOWTYPE::LEFTOVER_EM_PARTICLE );
908 
909  pflowContainer->AddParticleFlowElement( global_pflow_index, pflow );
910  global_pflow_index++;
911 
912  }
913 
914  } // close EM loop
915 
916  // now remove unmatched elements
917 
918  if ( Verbosity() > 2 )
919  std::cout << "ParticleFlowReco::process_event : remove TRK-unlinked EMs and HADs " << std::endl;
920 
921  for (unsigned int em = 0; em < _pflow_EM_E.size() ; em++ ) {
922 
923  // only consider EMs withOUT matched tracks ... we have dealt with the matched cases above
924  if ( _pflow_EM_match_TRK.at( em ).size() != 0 ) continue;
925 
926  if ( Verbosity() > 5 ) {
927  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;
928  }
929 
930  // add PFlow element for this EM
932 
933  // assume massless, could be updated to use K0L
934  TLorentzVector tlv; tlv.SetPtEtaPhiM( _pflow_EM_E[ em ] / cosh( _pflow_EM_eta[ em ] ) , _pflow_EM_eta[ em ] , _pflow_EM_phi[ em ] , 0 );
935 
936  pflow->set_px( tlv.Px() );
937  pflow->set_py( tlv.Py() );
938  pflow->set_pz( tlv.Pz() );
939  pflow->set_e( tlv.E() );
940  pflow->set_id( global_pflow_index );
941  pflow->set_type( ParticleFlowElement::PFLOWTYPE::UNMATCHED_EM_PARTICLE );
942 
943  pflowContainer->AddParticleFlowElement( global_pflow_index, pflow );
944  global_pflow_index++;
945 
946 
947  } // close EM loop
948 
949  for (unsigned int had = 0; had < _pflow_HAD_E.size() ; had++ ) {
950 
951  // only consider HADs withOUT matched tracks ... we have dealt with the matched cases above
952  if ( _pflow_HAD_match_TRK.at( had ).size() != 0 ) continue;
953 
954  if ( Verbosity() > 5 ) {
955  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;
956  }
957 
958  // add PFlow element for this HAD
960 
961  // assume massless, could be updated to use K0L
962  TLorentzVector tlv; tlv.SetPtEtaPhiM( _pflow_HAD_E[ had ] / cosh( _pflow_HAD_eta[ had ] ) , _pflow_HAD_eta[ had ] , _pflow_HAD_phi[ had ] , 0 );
963 
964  pflow->set_px( tlv.Px() );
965  pflow->set_py( tlv.Py() );
966  pflow->set_pz( tlv.Pz() );
967  pflow->set_e( tlv.E() );
968  pflow->set_id( global_pflow_index );
969  pflow->set_type( ParticleFlowElement::PFLOWTYPE::UNMATCHED_NEUTRAL_HADRON );
970 
971  pflowContainer->AddParticleFlowElement( global_pflow_index, pflow );
972  global_pflow_index++;
973 
974 
975  } // close HAD loop
976 
977  for (unsigned int trk = 0; trk < _pflow_TRK_p.size() ; trk++ ) {
978 
979  // only consider TRKs withOUT matched EM or HAD
980  if ( _pflow_TRK_match_EM.at( trk ).size() != 0 || _pflow_TRK_match_HAD.at( trk ).size() != 0 ) continue;
981 
982  if ( Verbosity() > 5 ) {
983  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;
984  }
985 
986  // add PFlow element for this TRK
988 
989  // assume massless, could be updated to use K0L
990  TLorentzVector tlv; tlv.SetPtEtaPhiM( _pflow_TRK_p[ trk ] / cosh( _pflow_TRK_eta[ trk ] ) , _pflow_TRK_eta[ trk ] , _pflow_TRK_phi[ trk ] , 0.135 );
991 
992  pflow->set_px( tlv.Px() );
993  pflow->set_py( tlv.Py() );
994  pflow->set_pz( tlv.Pz() );
995  pflow->set_e( tlv.E() );
996  pflow->set_id( global_pflow_index );
997  pflow->set_type( ParticleFlowElement::PFLOWTYPE::UNMATCHED_CHARGED_HADRON );
998 
999  pflowContainer->AddParticleFlowElement( global_pflow_index, pflow );
1000  global_pflow_index++;
1001 
1002  } // close TRK loop
1003 
1004  // DEBUG: print out all PFLow elements
1005  if ( Verbosity() > 5 ) {
1006  std::cout << "ParticleFlowReco::process_event: summary of PFlow objects " << std::endl;
1007 
1009  for ( ParticleFlowElementContainer::ConstIterator hiter = begin_end.first; hiter != begin_end.second; ++hiter) {
1010  hiter->second->identify();
1011  }
1012 
1013  }
1014 
1016 }
1017 
1019 {
1020  PHNodeIterator iter(topNode);
1021 
1022  // Looking for the DST node
1023  PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
1024  if (!dstNode)
1025  {
1026  std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
1028  }
1029 
1030  // store the PFlow elements under a sub-node directory
1031  PHCompositeNode *pflowNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "PARTICLEFLOW"));
1032  if (!pflowNode)
1033  {
1034  pflowNode = new PHCompositeNode("PARTICLEFLOW");
1035  dstNode->addNode( pflowNode );
1036  }
1037 
1038  // create the ParticleFlowElementContainer node...
1039  ParticleFlowElementContainer *pflowElementContainer = findNode::getClass<ParticleFlowElementContainer>(topNode, "ParticleFlowElements");
1040  if (!pflowElementContainer)
1041  {
1042  pflowElementContainer = new ParticleFlowElementContainer();
1043  PHIODataNode<PHObject> *pflowElementNode = new PHIODataNode<PHObject>(pflowElementContainer, "ParticleFlowElements", "PHObject");
1044  pflowNode->addNode( pflowElementNode );
1045  }
1046  else
1047  {
1048  std::cout << PHWHERE << "::ERROR - ParticleFlowElements node alerady exists, but should not" << std::endl;
1049  exit(-1);
1050  }
1051 
1053 }
1054 
1055 //____________________________________________________________________________..
1057 {
1058  if (Verbosity() > 0)
1059  {
1060  std::cout << "ParticleFlowReco::ResetEvent(PHCompositeNode *topNode) Resetting internal structures, prepare for next event" << std::endl;
1061  }
1063 }
1064 
1065 //____________________________________________________________________________..
1066 int ParticleFlowReco::EndRun(const int runnumber)
1067 {
1068  if (Verbosity() > 0)
1069  {
1070  std::cout << "ParticleFlowReco::EndRun(const int runnumber) Ending Run for Run " << runnumber << std::endl;
1071  }
1073 }
1074 
1075 //____________________________________________________________________________..
1077 {
1078  if (Verbosity() > 0)
1079  {
1080  std::cout << "ParticleFlowReco::End(PHCompositeNode *topNode) This is the End..." << std::endl;
1081  }
1083 }
1084 
1085 //____________________________________________________________________________..
1087 {
1088  if (Verbosity() > 0)
1089  {
1090  std::cout << "ParticleFlowReco::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
1091  }
1093 }
1094 
1095 //____________________________________________________________________________..
1096 void ParticleFlowReco::Print(const std::string &what) const
1097 {
1098  std::cout << "ParticleFlowReco::Print(const std::string &what) const Printing info for " << what << std::endl;
1099 }