Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
SoftLeptonTaggingTruth.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file SoftLeptonTaggingTruth.C
1 #include "SoftLeptonTaggingTruth.h"
2 
3 #include <fun4all/SubsysReco.h>
4 #include <fun4all/Fun4AllServer.h>
5 #include <fun4all/PHTFileServer.h>
6 #include <phool/PHCompositeNode.h>
7 #include <fun4all/Fun4AllReturnCodes.h>
8 #include <phool/getClass.h>
9 
10 #include <phool/PHCompositeNode.h>
11 
12 #include <g4main/PHG4TruthInfoContainer.h>
13 #include <g4main/PHG4VtxPoint.h>
14 #include <g4main/PHG4Particle.h>
15 #include <g4main/PHG4HitDefs.h>
16 
17 #include <g4eval/JetEvalStack.h>
18 #include <g4eval/JetTruthEval.h>
19 
20 #include <TString.h>
21 #include <TFile.h>
22 #include <TMath.h>
23 #include <TH1F.h>
24 #include <TH2F.h>
25 #include <TVector3.h>
26 #include <TLorentzVector.h>
27 #include <TAxis.h>
28 #include <TLine.h>
29 #include <TDatabasePDG.h>
30 
31 #include <exception>
32 #include <stdexcept>
33 #include <iostream>
34 #include <vector>
35 #include <set>
36 #include <algorithm>
37 #include <cassert>
38 #include <cmath>
39 
40 using namespace std;
41 
42 SoftLeptonTaggingTruth::SoftLeptonTaggingTruth(const std::string & truth_jet,
43  enu_flags flags) :
44  SubsysReco("SoftLeptonTaggingTruth_" + truth_jet), //
45  _jetevalstacks(), //
46  _truth_jet(truth_jet), _reco_jets(), _truth_container(NULL), _flags(flags), //
47  eta_range(-1, 1), //
48  _jet_match_dR(.7), _jet_match_dca(.2), _jet_match_E_Ratio(.0)
49 {
50 
51 }
52 
53 SoftLeptonTaggingTruth::~SoftLeptonTaggingTruth()
54 {
55 }
56 
57 int
58 SoftLeptonTaggingTruth::InitRun(PHCompositeNode *topNode)
59 {
60  _truth_container = findNode::getClass<PHG4TruthInfoContainer>(topNode,
61  "G4TruthInfo");
62  if (!_truth_container)
63  {
64  cout << "SoftLeptonTaggingTruth::InitRun - Fatal Error - "
65  << "unable to find DST node " << "G4TruthInfo" << endl;
66  assert(_truth_container);
67  }
68 
69  if (flag(kProcessTruthSpectrum))
70  {
71  if (not _jettrutheval)
72  _jettrutheval = shared_ptr < JetTruthEval
73  > (new JetTruthEval(topNode, _truth_jet));
74 
75  assert(_jettrutheval);
76  _jettrutheval->set_strict(true);
77  _jettrutheval->set_verbosity(verbosity + 1);
78  }
79 
80  return Fun4AllReturnCodes::EVENT_OK;
81 }
82 
83 int
84 SoftLeptonTaggingTruth::End(PHCompositeNode *topNode)
85 {
86 
87  return Fun4AllReturnCodes::EVENT_OK;
88 }
89 
90 int
91 SoftLeptonTaggingTruth::Init(PHCompositeNode *topNode)
92 {
93 
94  Fun4AllHistoManager *hm = getHistoManager();
95  assert(hm);
96 
97  if (flag(kProcessTruthSpectrum))
98  {
99  if (verbosity >= 1)
100  cout << "SoftLeptonTaggingTruth::Init - Process TruthSpectrum "
101  << _truth_jet << endl;
102  Init_Spectrum(topNode, _truth_jet);
103  }
104 
105  return Fun4AllReturnCodes::EVENT_OK;
106 }
107 
108 int
109 SoftLeptonTaggingTruth::process_event(PHCompositeNode *topNode)
110 {
111 
112  if (verbosity > 2)
113  cout << "SoftLeptonTaggingTruth::process_event() entered" << endl;
114 
115  for (jetevalstacks_map::iterator it_jetevalstack = _jetevalstacks.begin();
116  it_jetevalstack != _jetevalstacks.end(); ++it_jetevalstack)
117  {
118  assert(it_jetevalstack->second);
119  it_jetevalstack->second->next_event(topNode);
120  }
121  if (_jettrutheval)
122  _jettrutheval->next_event(topNode);
123 
124  if (flag(kProcessTruthSpectrum))
125  {
126  if (verbosity >= 1)
127  cout << "SoftLeptonTaggingTruth::process_event - Process TruthSpectrum "
128  << _truth_jet << endl;
129  process_Spectrum(topNode, _truth_jet, false);
130  }
131 
132  return Fun4AllReturnCodes::EVENT_OK;
133 }
134 
136 void
137 SoftLeptonTaggingTruth::set_eta_range(double low, double high)
138 {
139  if (low > high)
140  swap(low, high);
141  assert(low < high);
142  // eliminate zero range
143 
144  eta_range.first = low;
145  eta_range.second = high;
146 }
147 
150 TString
151 SoftLeptonTaggingTruth::get_eta_range_str(const char * eta_name) const
152 {
153  assert(eta_name);
154  return TString(
155  Form("%.1f < %s < %.1f", eta_range.first, eta_name, eta_range.second));
156 }
157 
159 bool
160 SoftLeptonTaggingTruth::jet_acceptance_cut(const Jet * jet) const
161 {
162  assert(jet);
163  bool eta_cut = (jet->get_eta() >= eta_range.first)
164  and (jet->get_eta() <= eta_range.second);
165  return eta_cut;
166 }
167 
168 string
169 SoftLeptonTaggingTruth::get_histo_prefix(const std::string & src_jet_name,
170  const std::string & reco_jet_name)
171 {
172  std::string histo_prefix = "h_QAG4SimJet_";
173 
174  if (src_jet_name.length() > 0)
175  {
176  histo_prefix += src_jet_name;
177  histo_prefix += "_";
178  }
179  if (reco_jet_name.length() > 0)
180  {
181  histo_prefix += reco_jet_name;
182  histo_prefix += "_";
183  }
184 
185  return histo_prefix;
186 }
187 
188 int
189 SoftLeptonTaggingTruth::Init_Spectrum(PHCompositeNode *topNode,
190  const std::string & jet_name)
191 {
192 
193  Fun4AllHistoManager *hm = getHistoManager();
194  assert(hm);
195 
196  // normalization plot with counts
197  const int norm_size = 3;
198 
199  TH1D * h_norm = new TH1D(
200  TString(get_histo_prefix(jet_name)) + "Normalization",
201  " Normalization;Item;Count", norm_size, .5, norm_size + .5);
202  int i = 1;
203 
204  h_norm->GetXaxis()->SetBinLabel(i++, "Event");
205  h_norm->GetXaxis()->SetBinLabel(i++, "Inclusive Jets");
206  h_norm->GetXaxis()->SetBinLabel(i++, "Leading Jets");
207  assert(norm_size >= i - 1);
208  h_norm->GetXaxis()->LabelsOption("v");
209  hm->registerHisto(h_norm);
210 
211  hm->registerHisto(
212  new TH1F(
213  //
214  TString(get_histo_prefix(jet_name)) + "Leading_Et", //
215  TString(jet_name) + " leading jet Et, " + get_eta_range_str()
216  + ";E_{T} (GeV)", 100, 0, 100));
217  hm->registerHisto(
218  new TH1F(
219  //
220  TString(get_histo_prefix(jet_name)) + "Leading_eta", //
221  TString(jet_name) + " leading jet #eta, " + get_eta_range_str()
222  + ";#eta", 50, -1, 1));
223  hm->registerHisto(
224  new TH1F(
225  //
226  TString(get_histo_prefix(jet_name)) + "Leading_phi", //
227  TString(jet_name) + " leading jet #phi, " + get_eta_range_str()
228  + ";#phi", 50, -M_PI, M_PI));
229 
230  TH1F * lcomp = new TH1F(
231  //
232  TString(get_histo_prefix(jet_name)) + "Leading_CompSize", //
233  TString(jet_name) + " leading jet # of component, " + get_eta_range_str()
234  + ";Number of component;", 100, 1, 3000);
235  useLogBins(lcomp->GetXaxis());
236  hm->registerHisto(lcomp);
237 
238  hm->registerHisto(
239  new TH1F(
240  //
241  TString(get_histo_prefix(jet_name)) + "Leading_Mass", //
242  TString(jet_name) + " leading jet mass, " + get_eta_range_str()
243  + ";Jet Mass (GeV);", 100, 0, 20));
244  hm->registerHisto(
245  new TH1F(
246  //
247  TString(get_histo_prefix(jet_name)) + "Leading_CEMC_Ratio", //
248  TString(jet_name) + " leading jet EMCal ratio, " + get_eta_range_str()
249  + ";Energy ratio CEMC/Total;", 100, 0, 1.01));
250  hm->registerHisto(
251  new TH1F(
252  //
253  TString(get_histo_prefix(jet_name)) + "Leading_CEMC_HCalIN_Ratio", //
254  TString(jet_name) + " leading jet EMCal+HCal_{IN} ratio, "
255  + get_eta_range_str() + ";Energy ratio (CEMC + HCALIN)/Total;",
256  100, 0, 1.01));
257 
258  // reco jet has no definition for leakages, since leakage is not reconstructed as part of jet energy.
259  // It is only available for truth jets
260  hm->registerHisto(
261  new TH1F(
262  //
263  TString(get_histo_prefix(jet_name)) + "Leading_Leakage_Ratio", //
264  TString(jet_name) + " leading jet leakage ratio, "
265  + get_eta_range_str() + ";Energy ratio, Back leakage/Total;", 100,
266  0, 1.01));
267 
268  TH1F * h = new TH1F(
269  //
270  TString(get_histo_prefix(jet_name)) + "Inclusive_E", //
271  TString(jet_name) + " inclusive jet E, " + get_eta_range_str()
272  + ";Total jet energy (GeV)", 100, 1e-3, 100);
273  useLogBins(h->GetXaxis());
274  hm->registerHisto(h);
275  hm->registerHisto(
276  new TH1F(
277  //
278  TString(get_histo_prefix(jet_name)) + "Inclusive_eta", //
279  TString(jet_name) + " inclusive jet #eta, " + get_eta_range_str()
280  + ";#eta;Jet energy density", 50, -1, 1));
281  hm->registerHisto(
282  new TH1F(
283  //
284  TString(get_histo_prefix(jet_name)) + "Inclusive_phi", //
285  TString(jet_name) + " inclusive jet #phi, " + get_eta_range_str()
286  + ";#phi;Jet energy density", 50, -M_PI, M_PI));
287 
288  TH2F * h2 = NULL;
289 
290  h2 = new TH2F(
291  //
292  TString(get_histo_prefix(jet_name)) + "Leading_Trk_PtRel", //
293  TString(jet_name) + " leading jet, " + get_eta_range_str()
294  + ";j_{T} (GeV/c);Particle Selection", 100, 0, 10, 4, 0.5, 4.5);
295  i = 1;
296  h2->GetYaxis()->SetBinLabel(i++, "Electron");
297  h2->GetYaxis()->SetBinLabel(i++, "Muon");
298  h2->GetYaxis()->SetBinLabel(i++, "Other Charged");
299  h2->GetYaxis()->SetBinLabel(i++, "Other Neutral");
300  hm->registerHisto(h2);
301 
302  h2 = new TH2F(
303  //
304  TString(get_histo_prefix(jet_name)) + "Leading_Trk_DCA2D", //
305  TString(jet_name) + " leading jet, " + get_eta_range_str()
306  + ";DCA_{2D} (cm);Particle Selection", 200, -.2, .2, 4, 0.5, 4.5);
307  i = 1;
308  h2->GetYaxis()->SetBinLabel(i++, "Electron");
309  h2->GetYaxis()->SetBinLabel(i++, "Muon");
310  h2->GetYaxis()->SetBinLabel(i++, "Other Charged");
311  h2->GetYaxis()->SetBinLabel(i++, "Other Neutral");
312  hm->registerHisto(h2);
313 
314  h2 = new TH2F(
315  //
316  TString(get_histo_prefix(jet_name)) + "Leading_Trk_PoverETotal", //
317  TString(jet_name) + " leading jet, " + get_eta_range_str()
318  + ";|P_{part.}|c/E_{Jet};Particle Selection", 120, 0, 1.2, 4, 0.5,
319  4.5);
320  i = 1;
321  h2->GetYaxis()->SetBinLabel(i++, "Electron");
322  h2->GetYaxis()->SetBinLabel(i++, "Muon");
323  h2->GetYaxis()->SetBinLabel(i++, "Other Charged");
324  h2->GetYaxis()->SetBinLabel(i++, "Other Neutral");
325  hm->registerHisto(h2);
326 
327  h2 = new TH2F(
328  //
329  TString(get_histo_prefix(jet_name)) + "Leading_Trk_dR", //
330  TString(jet_name) + " leading jet, " + get_eta_range_str()
331  + ";#DeltaR;Particle Selection", 100, 0, 1, 4, 0.5, 4.5);
332  i = 1;
333  h2->GetYaxis()->SetBinLabel(i++, "Electron");
334  h2->GetYaxis()->SetBinLabel(i++, "Muon");
335  h2->GetYaxis()->SetBinLabel(i++, "Other Charged");
336  h2->GetYaxis()->SetBinLabel(i++, "Other Neutral");
337  hm->registerHisto(h2);
338 
339  return Fun4AllReturnCodes::EVENT_OK;
340 }
341 
342 int
343 SoftLeptonTaggingTruth::process_Spectrum(PHCompositeNode *topNode,
344  const std::string & jet_name, const bool is_reco_jet)
345 {
346  JetMap* jets = findNode::getClass<JetMap>(topNode, jet_name.c_str());
347  if (!jets)
348  {
349  cout
350  << "SoftLeptonTaggingTruth::process_Spectrum - Error can not find DST JetMap node "
351  << jet_name << endl;
352  exit(-1);
353  }
354 
355  Fun4AllHistoManager *hm = getHistoManager();
356  assert(hm);
357 
358  TH1D* h_norm = dynamic_cast<TH1D*>(hm->getHisto(
359  get_histo_prefix(jet_name) + "Normalization"));
360  assert(h_norm);
361  h_norm->Fill("Event", 1);
362  h_norm->Fill("Inclusive Jets", jets->size());
363 
364  TH1F * ie = dynamic_cast<TH1F*>(hm->getHisto(
365  (get_histo_prefix(jet_name)) + "Inclusive_E" //
366  ));
367  assert(ie);
368  TH1F * ieta = dynamic_cast<TH1F*>(hm->getHisto(
369  (get_histo_prefix(jet_name)) + "Inclusive_eta" //
370  ));
371  assert(ieta);
372  TH1F * iphi = dynamic_cast<TH1F*>(hm->getHisto(
373  (get_histo_prefix(jet_name)) + "Inclusive_phi" //
374  ));
375  assert(iphi);
376 
377  Jet* leading_jet = NULL;
378  double max_et = 0;
379  for (JetMap::Iter iter = jets->begin(); iter != jets->end(); ++iter)
380  {
381  Jet* jet = iter->second;
382  assert(jet);
383 
384  if (not jet_acceptance_cut(jet))
385  continue;
386 
387  if (jet->get_et() > max_et)
388  {
389  leading_jet = jet;
390  max_et = jet->get_et();
391  }
392 
393  ie->Fill(jet->get_e());
394  ieta->Fill(jet->get_eta());
395  iphi->Fill(jet->get_phi());
396  }
397 
398  if (leading_jet)
399  {
400  if (verbosity)
401  {
402  cout
403  << "SoftLeptonTaggingTruth::process_Spectrum - processing leading jet with # comp = "
404  << leading_jet->size_comp() << endl;
405  leading_jet->identify();
406  }
407 
408  h_norm->Fill("Leading Jets", 1);
409 
410  TH1F * let = dynamic_cast<TH1F*>(hm->getHisto(
411  (get_histo_prefix(jet_name)) + "Leading_Et" //
412  ));
413  assert(let);
414  TH1F * leta = dynamic_cast<TH1F*>(hm->getHisto(
415  (get_histo_prefix(jet_name)) + "Leading_eta" //
416  ));
417  assert(leta);
418  TH1F * lphi = dynamic_cast<TH1F*>(hm->getHisto(
419  (get_histo_prefix(jet_name)) + "Leading_phi" //
420  ));
421  assert(lphi);
422 
423  TH1F * lcomp = dynamic_cast<TH1F*>(hm->getHisto(
424  (get_histo_prefix(jet_name)) + "Leading_CompSize" //
425  ));
426  assert(lcomp);
427  TH1F * lmass = dynamic_cast<TH1F*>(hm->getHisto(
428  (get_histo_prefix(jet_name)) + "Leading_Mass" //
429  ));
430  assert(lmass);
431  TH1F * lcemcr = dynamic_cast<TH1F*>(hm->getHisto(
432  (get_histo_prefix(jet_name)) + "Leading_CEMC_Ratio" //
433  ));
434  assert(lcemcr);
435  TH1F * lemchcalr = dynamic_cast<TH1F*>(hm->getHisto(
436  (get_histo_prefix(jet_name)) + "Leading_CEMC_HCalIN_Ratio" //
437  ));
438  assert(lemchcalr);
439  TH1F * lleak = dynamic_cast<TH1F*>(hm->getHisto(
440  (get_histo_prefix(jet_name)) + "Leading_Leakage_Ratio" //
441  ));
442  assert(lleak);
443 
444  let->Fill(leading_jet->get_et());
445  leta->Fill(leading_jet->get_eta());
446  lphi->Fill(leading_jet->get_phi());
447  lcomp->Fill(leading_jet->size_comp());
448  lmass->Fill(leading_jet->get_mass());
449 
450  if (is_reco_jet)
451  { // this is a reco jet
452  jetevalstacks_map::iterator it_stack = _jetevalstacks.find(jet_name);
453  assert(it_stack != _jetevalstacks.end());
454  shared_ptr<JetEvalStack> eval_stack = it_stack->second;
455  assert(eval_stack);
456  JetRecoEval* recoeval = eval_stack->get_reco_eval();
457  assert(recoeval);
458 
459  lcemcr->Fill( //
460  (recoeval->get_energy_contribution(leading_jet, Jet::CEMC_TOWER) //
461  +//
462  recoeval->get_energy_contribution(leading_jet,
463  Jet::CEMC_CLUSTER) //
464  ) / leading_jet->get_e());
465  lemchcalr->Fill( //
466  (recoeval->get_energy_contribution(leading_jet, Jet::CEMC_TOWER) //
467  +//
468  recoeval->get_energy_contribution(leading_jet,
469  Jet::CEMC_CLUSTER) //
470  +//
471  recoeval->get_energy_contribution(leading_jet,
472  Jet::HCALIN_TOWER) //
473  +//
474  recoeval->get_energy_contribution(leading_jet,
475  Jet::HCALIN_CLUSTER) //
476  ) / leading_jet->get_e());
477  // reco jet has no definition for leakages, since leakage is not reconstructed as part of jet energy.
478  // It is only available for truth jets
479  }
480  else
481  { // this is a truth jet
482  assert(_jettrutheval);
483 
484  double cemc_e = 0;
485  double hcalin_e = 0;
486  double bh_e = 0;
487 
488  set<PHG4Shower*> showers = _jettrutheval->all_truth_showers(
489  leading_jet);
490 
491  for (set<PHG4Shower*>::const_iterator it = showers.begin();
492  it != showers.end(); ++it)
493  {
494  cemc_e += (*it)->get_edep(
495  PHG4HitDefs::get_volume_id("G4HIT_CEMC"));
496  cemc_e += (*it)->get_edep(
497  PHG4HitDefs::get_volume_id("G4HIT_CEMC_ELECTRONICS"));
498  cemc_e += (*it)->get_edep(
499  PHG4HitDefs::get_volume_id("G4HIT_ABSORBER_CEMC"));
500 
501  hcalin_e += (*it)->get_edep(
502  PHG4HitDefs::get_volume_id("G4HIT_HCALIN"));
503  hcalin_e += (*it)->get_edep(
504  PHG4HitDefs::get_volume_id("G4HIT_ABSORBER_HCALIN"));
505 
506  bh_e += (*it)->get_edep(PHG4HitDefs::get_volume_id("G4HIT_BH_1"));
507  bh_e += (*it)->get_edep(
508  PHG4HitDefs::get_volume_id("G4HIT_BH_FORWARD_PLUS"));
509  bh_e += (*it)->get_edep(
510  PHG4HitDefs::get_volume_id("G4HIT_BH_FORWARD_NEG"));
511  }
512 
513  lcemcr->Fill(cemc_e / leading_jet->get_e());
514  lemchcalr->Fill((cemc_e + hcalin_e) / leading_jet->get_e());
515  lleak->Fill(bh_e / leading_jet->get_e());
516 
517  }
518 
519  TH2F * lptrel = dynamic_cast<TH2F*>(hm->getHisto(
520  (get_histo_prefix(jet_name)) + "Leading_Trk_PtRel" //
521  ));
522  assert(lptrel);
523  TH2F * ldca2d = dynamic_cast<TH2F*>(hm->getHisto(
524  (get_histo_prefix(jet_name)) + "Leading_Trk_DCA2D" //
525  ));
526  assert(ldca2d);
527  TH2F * lpoe = dynamic_cast<TH2F*>(hm->getHisto(
528  (get_histo_prefix(jet_name)) + "Leading_Trk_PoverETotal" //
529  ));
530  assert(lpoe);
531  TH2F * ldr = dynamic_cast<TH2F*>(hm->getHisto(
532  (get_histo_prefix(jet_name)) + "Leading_Trk_dR" //
533  ));
534  assert(ldr);
535 
536  // jet vector
537  TVector3 j_v(leading_jet->get_px(), leading_jet->get_py(),
538  leading_jet->get_pz());
539 
540  PHG4TruthInfoContainer::ConstRange primary_range =
541  _truth_container->GetPrimaryParticleRange();
542  for (PHG4TruthInfoContainer::ConstIterator particle_iter =
543  primary_range.first; particle_iter != primary_range.second;
544  ++particle_iter)
545  {
546  PHG4Particle *particle = particle_iter->second;
547  assert(particle);
548  const PHG4VtxPoint* primary_vtx = //
549  _truth_container->GetPrimaryVtx(particle->get_vtx_id());
550  assert(primary_vtx);
551 
552  TVector3 p_v(particle->get_px(), particle->get_py(),
553  particle->get_pz());
554 
555  TVector3 v_v(primary_vtx->get_x(), primary_vtx->get_y(),
556  primary_vtx->get_z());
557 
558  const double dR = sqrt(
559  (p_v.Eta() - j_v.Eta()) * (p_v.Eta() - j_v.Eta())
560  + (p_v.Phi() - j_v.Phi()) * (p_v.Phi() - j_v.Phi()));
561  const double E_Ratio = p_v.Mag() / leading_jet->get_e();
562 
563  const double charge3 = TDatabasePDG::Instance()->GetParticle(
564  particle->get_pid())->Charge();
565  TVector3 dca_v(cos(p_v.Phi() + TMath::Pi() / 2),
566  sin(p_v.Phi() + TMath::Pi() / 2), 0);
567  const double dca2d = v_v.Dot(dca_v) * copysign(1, charge3);
568 
569  if (dR > _jet_match_dR)
570  continue;
571  if (abs(dca2d) > _jet_match_dca)
572  continue;
573  if (E_Ratio < _jet_match_E_Ratio)
574  continue;
575 
576 
577  int histo_pid = 0;
578  if (abs(particle->get_pid())
579  == TDatabasePDG::Instance()->GetParticle("e-")->PdgCode())
580  histo_pid = ldr->GetYaxis()->FindBin("Electron");
581  else if (abs(particle->get_pid())
582  == TDatabasePDG::Instance()->GetParticle("mu-")->PdgCode())
583  histo_pid = ldr->GetYaxis()->FindBin("Muon");
584  else if (charge3 != 0)
585  histo_pid = ldr->GetYaxis()->FindBin("Other Charged");
586  else
587  histo_pid = ldr->GetYaxis()->FindBin("Other Neutral");
588  assert(histo_pid);
589 
590  if (verbosity)
591  {
592  cout
593  << "SoftLeptonTaggingTruth::process_Spectrum - particle hist ID "
594  << histo_pid << ", charge x 3 = " << charge3 << ", dR = "
595  << dR << ", E_Ratio = " << E_Ratio << "/"
596  << _jet_match_E_Ratio << " for ";
597  particle->identify();
598  }
599 
600  const double pT_rel = p_v.Dot(j_v) / j_v.Mag();
601 
602  ldr->Fill(dR, histo_pid);
603  lpoe->Fill(E_Ratio, histo_pid);
604  lptrel->Fill(pT_rel, histo_pid);
605  ldca2d->Fill(dca2d, histo_pid);
606  } // for (PHG4TruthInfoContainer::ConstIterator particle_iter =
607 
608  } // if (leading_jet)
609 
610  return Fun4AllReturnCodes::EVENT_OK;
611 }
612 
614 Fun4AllHistoManager *
615 SoftLeptonTaggingTruth::getHistoManager()
616 {
617 
618  Fun4AllServer *se = Fun4AllServer::instance();
619  Fun4AllHistoManager *hm = se->getHistoManager("histSoftLeptonTaggingTruth");
620 
621  if (not hm)
622  {
623 // cout
624 // << "QAHistManagerDef::get_HistoManager - Making Fun4AllHistoManager EMCalAna_HISTOS"
625 // << endl;
626  hm = new Fun4AllHistoManager("histSoftLeptonTaggingTruth");
627  se->registerHistoManager(hm);
628  }
629 
630  assert(hm);
631 
632  return hm;
633 }
634 
636 void
637 SoftLeptonTaggingTruth::useLogBins(TAxis * axis)
638 {
639  assert(axis);
640  assert(axis->GetXmin()>0);
641  assert( axis->GetXmax()>0);
642 
643  const int bins = axis->GetNbins();
644 
645  Axis_t from = log10(axis->GetXmin());
646  Axis_t to = log10(axis->GetXmax());
647  Axis_t width = (to - from) / bins;
648  vector<Axis_t> new_bins(bins + 1);
649 
650  for (int i = 0; i <= bins; i++)
651  {
652  new_bins[i] = TMath::Power(10, from + i * width);
653  }
654  axis->Set(bins, new_bins.data());
655 
656 }