Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ElectronID.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file ElectronID.cc
1 #include "ElectronID.h"
2 
3 #include <cstdlib>
4 #include <cstdio>
5 #include <iomanip>
6 #include <fstream>
7 
8 #include <phool/PHCompositeNode.h>
9 #include <phool/recoConsts.h>
10 #include <phool/PHIODataNode.h> // for PHIODataNode
11 #include <phool/PHNode.h> // for PHNode
12 #include <phool/PHNodeIterator.h>
13 #include <phool/PHObject.h> // for PHObject
14 #include <phool/getClass.h>
15 #include <phool/phool.h> // for PHWHERE
16 #include <phool/PHRandomSeed.h>
17 #include <phool/getClass.h>
18 
19 #include <fun4all/Fun4AllReturnCodes.h>
20 
21 #include <trackbase_historic/SvtxTrackMap.h>
22 #include <trackbase_historic/SvtxTrack.h>
23 #include <trackbase_historic/SvtxVertex.h>
24 #include <trackbase_historic/SvtxVertexMap.h>
25 #include <trackbase/TrkrDefs.h>
26 
27 #include <calobase/RawClusterContainer.h>
28 #include <calobase/RawCluster.h>
29 #include <calobase/RawClusterv1.h>
30 
31 #include <globalvertex/GlobalVertexMap.h>
32 #include <globalvertex/GlobalVertex.h>
33 
34 #include <g4main/PHG4TruthInfoContainer.h>
35 #include <g4main/PHG4Particle.h>
36 #include <g4main/PHG4VtxPoint.h>
37 
38 #include "trackpidassoc/TrackPidAssoc.h"
39 
40 //TMVA class
41 #include <vector>
42 #include <iostream>
43 #include <map>
44 #include <string>
45 #include "TMVA/Tools.h"
46 #include "TMVA/Reader.h"
47 #include "TMVA/MethodCuts.h"
48 #include "TMVA/Factory.h"
49 #include "TMVA/DataLoader.h"
50 
51 
52 // gsl
53 #include <gsl/gsl_randist.h>
54 #include <gsl/gsl_rng.h>
55 
56 #include <TVector3.h>
57 #include <TMatrixFfwd.h> // for TMatrixF
58 
59 #include <iostream> // for operator<<, basic_ostream
60 #include <set> // for _Rb_tree_iterator, set
61 #include <utility> // for pair
62 
63 class PHCompositeNode;
64 
65 using namespace std;
66 
67 ElectronID::ElectronID(const std::string& name, const std::string &filename) : SubsysReco(name)
68 {
69  OutputNtupleFile=nullptr;
70  OutputFileName=filename;
71  EventNumber=0;
72  output_ntuple = true;
73 
75  EMOP_lowerlimit = 0.7;
76  EMOP_higherlimit = 100.0;
77  HOP_lowerlimit = 0.0;
78  HinOEM_higherlimit = 100.0;
79  Pt_lowerlimit = 0.0;
80  Pt_higherlimit = 100.0;
81  Nmvtx_lowerlimit = 2;
82  Nintt_lowerlimit = 0;
83  Ntpc_lowerlimit = 0;
84  Nquality_higherlimit = 100;
85  Ntpc_lowerlimit = 20;
86  Nquality_higherlimit = 5.;
87  PROB_cut = 0.;
88 
90  BDT_cut_p = 0.0;
91  BDT_cut_n = 0.0;
92  ISUSE_BDT_p =0;
93  ISUSE_BDT_n =0;
94 
95  // unsigned int _nlayers_maps = 3;
96  // unsigned int _nlayers_intt = 4;
97  // unsigned int _nlayers_tpc = 48;
98 }
99 
100 ElectronID::~ElectronID()
101 {
102 }
103 
104 int ElectronID::Init(PHCompositeNode *topNode)
105 {
106 
107  if(output_ntuple) {
108 
109  OutputNtupleFile = new TFile(OutputFileName.c_str(),"RECREATE");
110  std::cout << "PairMaker::Init: output file " << OutputFileName.c_str() << " opened." << endl;
111  ntpBDTresponse = new TNtuple("ntpbeforecut","","select_p:select_n");
112 
113  ntpbeforecut = new TNtuple("ntpbeforecut","","p:pt:cemce3x3:hcaline3x3:hcaloute3x3:cemce3x3overp:hcaline3x3overcemce3x3:hcale3x3overp:charge:pid:quality:e_cluster:EventNumber:z:vtxid:nmvtx:nintt:ntpc:cemc_prob:cemc_ecore:cemc_chi2");
114  ntpcutEMOP = new TNtuple("ntpcutEMOP","","p:pt:cemce3x3:hcaline3x3:hcaloute3x3:charge:pid:quality:e_cluster:EventNumber:z:vtxid:cemc_prob:cemc_ecore");
115  ntpcutHOP = new TNtuple("ntpcutHOP","","p:pt:cemce3x3:hcaline3x3:hcaloute3x3:charge:pid:quality:e_cluster:EventNumber:z:vtxid:cemc_prob:cemc_ecore");
116  ntpcutEMOP_HinOEM = new TNtuple("ntpcutEMOP_HinOEM","","p:pt:cemce3x3:hcaline3x3:hcaloute3x3:charge:pid:quality:e_cluster:EventNumber:z:vtxid:cemc_prob:cemc_ecore");
117  ntpcutEMOP_HinOEM_Pt = new TNtuple("ntpcutEMOP_HinOEM_Pt","","p:pt:cemce3x3:hcaline3x3:hcaloute3x3:charge:pid:quality:e_cluster:EventNumber:z:vtxid:cemc_prob:cemc_ecore");
118  ntpcutEMOP_HinOEM_Pt_read = new TNtuple("ntpcutEMOP_HinOEM_Pt_read","","p:pt:cemce3x3:hcaline3x3:hcaloute3x3:charge:pid:quality:e_cluster:EventNumber:z:vtxid:trackid:cemc_prob:cemc_ecore");
119  ntpcutBDT_read = new TNtuple("ntpcutBDT_read","","p:pt:cemce3x3:hcaline3x3:hcaloute3x3:charge:pid:quality:e_cluster:EventNumber:z:vtxid:trackid:cemc_prob:cemc_ecore:EOP:HOM:cemc_chi2");
120  }
121  else {
122  PHNodeIterator iter(topNode);
123  PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
124  if (!dstNode)
125  {
126  cerr << PHWHERE << " ERROR: Can not find DST node." << endl;
127  return Fun4AllReturnCodes::ABORTEVENT;
128  }
129  }
130 
131  return Fun4AllReturnCodes::EVENT_OK;
132 
133 }
134 
135 int ElectronID::InitRun(PHCompositeNode* topNode)
136 {
137  int ret = GetNodes(topNode);
138  if (ret != Fun4AllReturnCodes::EVENT_OK) return ret;
139 
140  return Fun4AllReturnCodes::EVENT_OK;
141 }
142 
143 int ElectronID::process_event(PHCompositeNode* topNode)
144 {
145  EventNumber++;
146  float ntp[30];
147 
148  SvtxTrackMap *trackmap = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
149  if(!trackmap) {
150  cerr << PHWHERE << " ERROR: Can not find SvtxTrackMap node." << endl;
151  return Fun4AllReturnCodes::ABORTEVENT;
152  }
153 
154  cout<<"EventNumber ===================== " << EventNumber-1 << endl;
155  if(EventNumber==1) topNode->print();
156 
157  SvtxVertexMap *vtxmap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
158  if(!vtxmap) {
159  cout << "SvtxVertexMap node not found!" << endl;
160  return Fun4AllReturnCodes::ABORTEVENT;
161  }
162  //cout << "Number of SvtxVertexMap entries = " << vtxmap->size() << endl;
163 
164  RawClusterContainer* cemc_cluster_container = findNode::getClass<RawClusterContainer>(topNode, "CLUSTER_CEMC");
165  if(!cemc_cluster_container) {
166  cerr << PHWHERE << " ERROR: Can not find CLUSTER_CEMC node." << endl;
167  return Fun4AllReturnCodes::ABORTEVENT;
168  }
169 
170 
171 // Truth info
172  PHG4TruthInfoContainer* truth_container = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
173  if(!truth_container) {
174  cerr << PHWHERE << " ERROR: Can not find G4TruthInfo node." << endl;
175  return Fun4AllReturnCodes::ABORTEVENT;
176  }
177  PHG4TruthInfoContainer::ConstRange range = truth_container->GetPrimaryParticleRange();
178  cout << "number of MC particles = " << truth_container->size() << " " << truth_container->GetNumPrimaryVertexParticles() << endl;
179 
180  int mycount = 0;
181  for (PHG4TruthInfoContainer::ConstIterator iter = range.first; iter != range.second; ++iter)
182  {
183  PHG4Particle* g4particle = iter->second;
184  mycount++;
185  int gflavor = g4particle->get_pid();
186  double gpx = g4particle->get_px();
187  double gpy = g4particle->get_py();
188  double gpz= g4particle->get_pz();
189  double gpt = sqrt(gpx*gpx+gpy*gpy);
190  double phi = atan2(gpy,gpx);
191  double eta = asinh(gpz/gpt);
192  int primid = g4particle->get_primary_id();
193  int parentid = g4particle->get_parent_id();
194  int trackid = g4particle->get_track_id();
195  if(trackid>truth_container->GetNumPrimaryVertexParticles()-50) cout << trackid << " " << parentid << " " << primid << " " << gflavor << " " << gpt << " " << phi << " " << eta << endl;
196  }
197  cout << "mycount = " << mycount << endl;
198 // end Truth
199 
200 
201 //MVA method setup**********************
202 
203  // This loads the library
204  TMVA::Tools::Instance();
205 
206  // Default MVA methods with trained + tested
207  std::map<std::string,int> Use_p;
208  std::map<std::string,int> Use_n;
209  // Boosted Decision Trees
210  Use_p["BDT"] = ISUSE_BDT_p; // uses Adaptive Boost; for positive particles
211  Use_n["BDT"] = ISUSE_BDT_n; // uses Adaptive Boost; for negative particles
212 
213  // Create the Reader object
214 
215  TMVA::Reader *reader_positive = new TMVA::Reader( "!Color:!Silent" );
216  TMVA::Reader *reader_negative = new TMVA::Reader( "!Color:!Silent" );
217 
218  // Create a set of variables and declare them to the reader
219  // - the variable names MUST corresponds in name and type to those given in the weight file(s) used
220  Float_t var1_p, var2_p, var3_p;
221  Float_t var1_n, var2_n, var3_n;
222  reader_positive->AddVariable( "var1_p", &var1_p );
223  reader_positive->AddVariable( "var2_p", &var2_p );
224  reader_positive->AddVariable( "var3_p", &var3_p );
225 
226  reader_negative->AddVariable( "var1_n", &var1_n );
227  reader_negative->AddVariable( "var2_n", &var2_n );
228  reader_negative->AddVariable( "var3_n", &var3_n );
229 
230  // Book the MVA methods
231  TString dir_p, dir_n;
232  dir_p = "dataset/Weights_positive/"; // weights for positive particles
233  dir_n = "dataset/Weights_negative/"; // weights for negative particles
234  TString prefix = "TMVAClassification";
235 
236  // Book method(s)
237  for (std::map<std::string,int>::iterator it = Use_p.begin(); it != Use_p.end(); it++) {
238  if (it->second) {
239  TString methodName_p = TString(it->first) + TString(" method");
240  TString weightfile_p = dir_p + prefix + TString("_") + TString(it->first) + TString(".weights.xml");
241  reader_positive->BookMVA( methodName_p, weightfile_p );
242  }
243  }
244 
245  for (std::map<std::string,int>::iterator it = Use_n.begin(); it != Use_n.end(); it++) {
246  if (it->second) {
247  TString methodName_n = TString(it->first) + TString(" method");
248  TString weightfile_n = dir_n + prefix + TString("_") + TString(it->first) + TString(".weights.xml");
249  reader_negative->BookMVA( methodName_n, weightfile_n );
250  }
251  }
252 
253 //end of MVA method setup*******************
254 
255  int nmvtx = 0;
256  int nintt = 0;
257  int ntpc = 0;
258  cout << _nlayers_maps << " " << _nlayers_intt << " " << _nlayers_tpc<<endl;
259  // get the tracks
260  for(SvtxTrackMap::Iter it = _track_map->begin(); it != _track_map->end(); ++it)
261  {
262  SvtxTrack *track = it->second;
263 
264  PHG4Particle* g4particle = findMCmatch(track, truth_container);
265  int gflavor = g4particle->get_pid();
266  if(gflavor==0) continue;
267 
268  nmvtx = 0;
269  nintt = 0;
270  ntpc = 0;
271 
272  for (SvtxTrack::ConstClusterKeyIter iter = track->begin_cluster_keys(); iter != track->end_cluster_keys(); ++iter)
273  {
274  TrkrDefs::cluskey cluser_key = *iter;
275  int trackerid = TrkrDefs::getTrkrId(cluser_key);
276  // cout << "trackerid= " << trackerid << endl;
277  if(trackerid==0) nmvtx++;
278  if(trackerid==1) nintt++;
279  if(trackerid==2) ntpc++;
280  //unsigned int layer = TrkrDefs::getLayer(cluser_key);
281  //if (_nlayers_maps > 0 && layer < _nlayers_maps) nmvtx++;
282  //if (_nlayers_intt > 0 && layer >= _nlayers_maps && layer < _nlayers_maps + _nlayers_intt) nintt++;
283  //if (_nlayers_tpc > 0 && layer >= (_nlayers_maps + _nlayers_intt) && layer < (_nlayers_maps + _nlayers_intt + _nlayers_tpc)) ntpc++;
284  }
285 
286 
287  double px = track->get_px();
288  double py = track->get_py();
289 
290  double z = track->get_z();
291 
292  unsigned int vtxid = track->get_vertex_id();
293 
294  double mom = track->get_p();
295  double pt = sqrt(px*px + py*py);
296  int charge = track->get_charge();
297  int pid = it->first;
298  float quality = track->get_quality();
299 
300  double e_cluster = track->get_cal_cluster_e(SvtxTrack::CAL_LAYER::CEMC);
301 
302  double e_cemc_3x3 = track->get_cal_energy_3x3(SvtxTrack::CAL_LAYER::CEMC);
303  double e_hcal_in_3x3 = track->get_cal_energy_3x3(SvtxTrack::CAL_LAYER::HCALIN);
304  double e_hcal_out_3x3 = track->get_cal_energy_3x3(SvtxTrack::CAL_LAYER::HCALOUT);
305 
306  unsigned int cemc_clusid = track->get_cal_cluster_id(SvtxTrack::CAL_LAYER::CEMC);
307  double cemc_prob = 0.;
308  double cemc_ecore = 0.;
309  double cemc_chi2 = 99999.;
310  if(cemc_clusid<99999) {
311  RawCluster* cemc_cluster = cemc_cluster_container->getCluster(cemc_clusid);
312  cemc_prob = cemc_cluster->get_prob();
313  cemc_ecore = cemc_cluster->get_ecore();
314  cemc_chi2 = cemc_cluster->get_chi2();
315  }
316 
317  // CEMC E/p cut
318  double cemceoverp = e_cemc_3x3 / mom;
319  //double cemceoverp = cemc_ecore / mom;
320  // HCaline/CEMCe cut
321  double hcalineovercemce = e_hcal_in_3x3 / e_cemc_3x3;
322  // HCal E/p cut
323  double hcaleoverp = (e_hcal_in_3x3 + e_hcal_out_3x3) / mom;
324 
325  //MVA valuables
326  if(cemceoverp>0.0 && cemceoverp<10.0 && hcalineovercemce>0.0 && hcalineovercemce<10.0){
327  if(charge>0){
328  var1_p = cemceoverp;
329  var2_p = hcalineovercemce;
330  var3_p = cemc_chi2;
331  }
332  if(charge<0){
333  var1_n = cemceoverp;
334  var2_n = hcalineovercemce;
335  var3_n = cemc_chi2;
336  }
337  }
338 
339  if (Use_p["BDT"] && Use_n["BDT"] && quality < Nquality_higherlimit && nmvtx >= Nmvtx_lowerlimit && nintt >= Nintt_lowerlimit && ntpc >= Ntpc_lowerlimit && pt > Pt_lowerlimit && pt < Pt_higherlimit &&cemceoverp>0.0 && cemceoverp<10.0 && hcalineovercemce>0.0 && hcalineovercemce<10.0) {
340  float select_p=reader_positive->EvaluateMVA("BDT method");
341  float select_n=reader_negative->EvaluateMVA("BDT method");
342  ntp[0] = select_p;
343  ntp[1] = select_n;
344  if(output_ntuple) { ntpBDTresponse -> Fill(ntp); }
345  if(select_p>BDT_cut_p && select_n>BDT_cut_n){
346  // add to the association map
347  _track_pid_assoc->addAssoc(TrackPidAssoc::electron, it->second->get_id());
348  }
349 
350  }// end of BDT cut
351  else{ // for traditional cuts
352 
353  ntp[0] = mom;
354  ntp[1] = pt;
355  ntp[2] = e_cemc_3x3;
356  ntp[3] = e_hcal_in_3x3;
357  ntp[4] = e_hcal_out_3x3;
358  ntp[5] = cemceoverp;
359  ntp[6] = hcalineovercemce;
360  ntp[7] = hcaleoverp;
361  ntp[8] = charge;
362  ntp[9] = pid;
363  ntp[10] = quality;
364  ntp[11] = e_cluster;
365  ntp[12] = EventNumber;
366  ntp[13] = z;
367  ntp[14] = vtxid;
368  ntp[15] = nmvtx;
369  ntp[16] = nintt;
370  ntp[17] = ntpc;
371  ntp[18] = cemc_prob;
372  ntp[19] = cemc_ecore;
373  ntp[20] = cemc_chi2;
374  if(output_ntuple) { ntpbeforecut -> Fill(ntp); }
375 
376  //std::cout << " Pt_lowerlimit " << Pt_lowerlimit << " Pt_higherlimit " << Pt_higherlimit << " HOP_lowerlimit " << HOP_lowerlimit <<std::endl;
377  //std::cout << " EMOP_lowerlimit " << EMOP_lowerlimit << " EMOP_higherlimit " << EMOP_higherlimit << " HinOEM_higherlimit " << HinOEM_higherlimit <<std::endl;
378 
380  if(cemceoverp > EMOP_lowerlimit && cemceoverp < EMOP_higherlimit && quality < Nquality_higherlimit && nmvtx >= Nmvtx_lowerlimit && nintt >= Nintt_lowerlimit && ntpc >= Ntpc_lowerlimit)
381  {
382 
383  ntp[0] = mom;
384  ntp[1] = pt;
385  ntp[2] = e_cemc_3x3;
386  ntp[3] = e_hcal_in_3x3;
387  ntp[4] = e_hcal_out_3x3;
388  ntp[5] = charge;
389  ntp[6] = pid;
390  ntp[7] = quality;
391  ntp[8] = e_cluster;
392  ntp[9] = EventNumber;
393  ntp[10] = z;
394  ntp[11] = vtxid;
395  ntp[12] = cemc_prob;
396  ntp[13] = cemc_ecore;
397  if(output_ntuple) { ntpcutEMOP -> Fill(ntp); }
398 
399  if(hcalineovercemce < HinOEM_higherlimit)
400  {
401 
402  ntp[0] = mom;
403  ntp[1] = pt;
404  ntp[2] = e_cemc_3x3;
405  ntp[3] = e_hcal_in_3x3;
406  ntp[4] = e_hcal_out_3x3;
407  ntp[5] = charge;
408  ntp[6] = pid;
409  ntp[7] = quality;
410  ntp[8] = e_cluster;
411  ntp[9] = EventNumber;
412  ntp[10] = z;
413  ntp[11] = vtxid;
414  ntp[12] = cemc_prob;
415  ntp[13] = cemc_ecore;
416  if(output_ntuple) { ntpcutEMOP_HinOEM -> Fill(ntp); }
417 
418  if( pt > Pt_lowerlimit && pt < Pt_higherlimit)
419  {
420 
421  ntp[0] = mom;
422  ntp[1] = pt;
423  ntp[2] = e_cemc_3x3;
424  ntp[3] = e_hcal_in_3x3;
425  ntp[4] = e_hcal_out_3x3;
426  ntp[5] = charge;
427  ntp[6] = pid;
428  ntp[7] = quality;
429  ntp[8] = e_cluster;
430  ntp[9] = EventNumber;
431  ntp[10] = z;
432  ntp[11] = vtxid;
433  ntp[12] = cemc_prob;
434  ntp[13] = cemc_ecore;
435  if(output_ntuple) { ntpcutEMOP_HinOEM_Pt -> Fill(ntp); }
436 
437  if(Verbosity() > 0) {
438  std::cout << " Track " << it->first << " identified as electron " << " mom " << mom << " e_cemc_3x3 " << e_cemc_3x3 << " cemceoverp " << cemceoverp << " e_hcal_in_3x3 " << e_hcal_in_3x3 << " e_hcal_out_3x3 " << e_hcal_out_3x3 << std::endl;
439  }
440 
441  // add to the association map
442  _track_pid_assoc->addAssoc(TrackPidAssoc::electron, it->second->get_id());
443  }
444  }
445  }
446  }//end of BDT else
447 
449 
450  // if(hcaleoverp > HOP_lowerlimit && quality < 10)
451  if(hcaleoverp > HOP_lowerlimit && quality < 5)
452  {
453 
454  ntp[0] = mom;
455  ntp[1] = pt;
456  ntp[2] = e_cemc_3x3;
457  ntp[3] = e_hcal_in_3x3;
458  ntp[4] = e_hcal_out_3x3;
459  ntp[5] = charge;
460  ntp[6] = pid;
461  ntp[7] = quality;
462  ntp[8] = e_cluster;
463  ntp[9] = EventNumber;
464  ntp[10] = z;
465  ntp[11] = vtxid;
466  ntp[12] = cemc_prob;
467  ntp[13] = cemc_ecore;
468  if(output_ntuple) { ntpcutHOP -> Fill(ntp); }
469 
470  if(Verbosity() > 0) {
471  std::cout << " Track " << it->first << " identified as hadron " << " mom " << mom << " e_cemc_3x3 " << e_cemc_3x3 << " hcaleoverp " << hcaleoverp << " e_hcal_in_3x3 " << e_hcal_in_3x3 << " e_hcal_out_3x3 " << e_hcal_out_3x3 << std::endl;
472  }
473 
474  // add to the association map
475  _track_pid_assoc->addAssoc(TrackPidAssoc::hadron, it->second->get_id());
476  }
477 
478 
479  }//end of evet loop.
480 
481  delete reader_positive;
482  delete reader_negative;
483 
484  // Read back the association map
485  if(Verbosity() > 1)
486  std::cout << "Read back the association map electron entries" << std::endl;
487  auto electrons = _track_pid_assoc->getTracks(TrackPidAssoc::electron);
488  for(auto it = electrons.first; it != electrons.second; ++it)
489  {
490  SvtxTrack *tr = _track_map->get(it->second);
491 
492  double z = tr->get_z();
493  unsigned int vtxid = tr->get_vertex_id();
494  double mom = tr->get_p();
495  double px = tr->get_px();
496  double py = tr->get_py();
497  double pt = sqrt(px*px + py*py);
498  int charge = tr->get_charge();
499  int pid = it->first;
500  int tr_id = it->second;
501  int quality = tr->get_quality();
502 
503  double e_cluster = tr->get_cal_cluster_e(SvtxTrack::CAL_LAYER::CEMC);
504 
505  double e_cemc_3x3 = tr->get_cal_energy_3x3(SvtxTrack::CAL_LAYER::CEMC);
506  double e_hcal_in_3x3 = tr->get_cal_energy_3x3(SvtxTrack::CAL_LAYER::HCALIN);
507  double e_hcal_out_3x3 = tr->get_cal_energy_3x3(SvtxTrack::CAL_LAYER::HCALOUT);
508 
509  double EOP = e_cemc_3x3/mom;
510  double HOM = e_hcal_in_3x3/e_cemc_3x3;
511 
512 
513  unsigned int cemc_clusid = tr->get_cal_cluster_id(SvtxTrack::CAL_LAYER::CEMC);
514  double cemc_prob = 0.;
515  double cemc_ecore = 0.;
516  double cemc_chi2 = 99999.;
517  if(cemc_clusid<99999) {
518  RawCluster* cemc_cluster = cemc_cluster_container->getCluster(cemc_clusid);
519  cemc_prob = cemc_cluster->get_prob();
520  cemc_ecore = cemc_cluster->get_ecore();
521  cemc_chi2 = cemc_cluster->get_chi2();
522  }
523 
524  ntp[0] = mom;
525  ntp[1] = pt;
526  ntp[2] = e_cemc_3x3;
527  ntp[3] = e_hcal_in_3x3;
528  ntp[4] = e_hcal_out_3x3;
529  ntp[5] = charge;
530  ntp[6] = pid;
531  ntp[7] = quality;
532  ntp[8] = e_cluster;
533  ntp[9] = EventNumber;
534  ntp[10] = z;
535  ntp[11] = vtxid;
536  ntp[12] = tr_id;
537  ntp[13] = cemc_prob;
538  ntp[14] = cemc_ecore;
539  ntp[15] = EOP;
540  ntp[16] = HOM;
541  ntp[17] = cemc_chi2;
542  //if(output_ntuple) { ntpcutEMOP_HinOEM_Pt_read -> Fill(ntp); }
543  if(output_ntuple) { ntpcutBDT_read -> Fill(ntp); }
544 
545  if(Verbosity() > 1)
546  std::cout << " pid " << it->first << " track ID " << it->second << " mom " << mom << std::endl;
547  }
548 
549  if(Verbosity() > 1)
550  std::cout << "Read back the association map hadron entries" << std::endl;
551  auto hadrons = _track_pid_assoc->getTracks(TrackPidAssoc::hadron);
552  for(auto it = hadrons.first; it != hadrons.second; ++it)
553  {
554  SvtxTrack *tr = _track_map->get(it->second);
555  double p = tr->get_p();
556 
557  if(Verbosity() > 1)
558  std::cout << " pid " << it->first << " track ID " << it->second << " mom " << p << std::endl;
559  }
560 
561 
562  return Fun4AllReturnCodes::EVENT_OK;
563 }
564 
565 
566 int ElectronID::GetNodes(PHCompositeNode* topNode)
567 {
568  _track_map = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
569 
570  // if the TrackPidAssoc node is not present, create it
571  _track_pid_assoc = findNode::getClass<TrackPidAssoc>(topNode, "TrackPidAssoc");
572  if(!_track_pid_assoc)
573  {
574  PHNodeIterator iter(topNode);
575  PHCompositeNode* dstNode = static_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
576  if (!dstNode)
577  {
578  cerr << PHWHERE << "DST Node missing, quit!" << endl;
579  return Fun4AllReturnCodes::ABORTEVENT;
580  }
581 
582  // Get the SVTX node
583  PHNodeIterator iter_dst(dstNode);
584  PHCompositeNode* tb_node =
585  dynamic_cast<PHCompositeNode*>(iter_dst.findFirst("PHCompositeNode", "SVTX"));
586  if (!tb_node)
587  {
588  cout << PHWHERE << "SVTX node missing, quit!" << endl;
589  return Fun4AllReturnCodes::ABORTEVENT;
590  }
591 
592  // now add the new node
593  _track_pid_assoc = new TrackPidAssoc;
594  PHIODataNode<PHObject>* assoc_node = new PHIODataNode<PHObject>(_track_pid_assoc, "TrackPidAssoc", "PHObject");
595  tb_node->addNode(assoc_node);
596  if (Verbosity() > 0)
597  cout << PHWHERE << "Svtx/TrackPidAssoc node added" << endl;
598  }
599 
600  return Fun4AllReturnCodes::EVENT_OK;
601 }
602 
603 
604 PHG4Particle* ElectronID::findMCmatch(SvtxTrack* track, PHG4TruthInfoContainer* truth_container)
605 {
606  double px = track->get_px();
607  double py = track->get_py();
608  double pz = track->get_pz();
609  double pt = sqrt(px*px+py*py);
610  double phi = atan2(py,px);
611  double eta = asinh(pz/pt);
612  PHG4TruthInfoContainer::ConstRange range = truth_container->GetPrimaryParticleRange();
613 // cout << "number of MC particles = " << truth_container->size() << " " << truth_container->GetNumPrimaryVertexParticles() << endl;
614 
615  double thedistance = 9999.;
616  PHG4Particle* matchedMC = NULL;
617 
618  for (PHG4TruthInfoContainer::ConstIterator iter = range.first; iter != range.second; ++iter)
619  {
620  PHG4Particle* g4particle = iter->second;
621  int trackid = g4particle->get_track_id();
622  if(trackid<=truth_container->GetNumPrimaryVertexParticles()-50) continue; // only embedded particles
623  //int gflavor = g4particle->get_pid();
624  double gpx = g4particle->get_px();
625  double gpy = g4particle->get_py();
626  double gpz= g4particle->get_pz();
627  double gpt = sqrt(gpx*gpx+gpy*gpy);
628  double gphi = atan2(gpy,gpx);
629  double geta = asinh(gpz/gpt);
630  //int primid = g4particle->get_primary_id();
631  //int parentid = g4particle->get_parent_id();
632  if(sqrt(pow(gphi-phi,2)+pow(geta-eta,2))<thedistance) {
633  thedistance = sqrt(pow(gphi-phi,2)+pow(geta-eta,2));
634  matchedMC = g4particle;
635  }
636  }
637 
638  return matchedMC;
639 }
640 
641 
642 int ElectronID::End(PHCompositeNode * /*topNode*/)
643 {
644 if(output_ntuple) {
645  OutputNtupleFile -> cd();
646  OutputNtupleFile -> Write();
647  OutputNtupleFile -> Close();
648 }
649 
650  cout << "************END************" << endl;
651  return Fun4AllReturnCodes::EVENT_OK;
652 }
653