Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
BJetModule.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file BJetModule.cc
1 
8 #include "BJetModule.h"
9 
10 #include <hfjettruthgeneration/HFJetDefs.h>
11 
12 #include <Acts/Definitions/Algebra.hpp>
13 
14 #include <fun4all/Fun4AllReturnCodes.h>
15 #include <fun4all/Fun4AllServer.h>
16 
17 #include <phool/PHCompositeNode.h>
18 #include <phool/getClass.h>
19 
20 #include <g4main/PHG4Hit.h>
21 #include <g4main/PHG4Particle.h>
22 #include <g4main/PHG4TruthInfoContainer.h>
23 #include <g4main/PHG4VtxPoint.h>
24 
25 #include <g4jets/Jet.h>
26 #include <g4jets/JetMap.h>
27 
28 #include <trackbase_historic/SvtxTrack.h>
29 #include <trackbase_historic/SvtxTrackMap.h>
30 #include <trackbase_historic/SvtxVertex.h>
31 #include <trackbase_historic/SvtxVertexMap.h>
32 
33 #include <g4eval/JetEvalStack.h>
34 #include <g4eval/SvtxEvalStack.h>
35 
36 #include <HepMC/GenEvent.h>
37 #include <HepMC/GenVertex.h>
38 #include <phhepmc/PHHepMCGenEvent.h>
39 #include <phhepmc/PHHepMCGenEventMap.h>
40 
41 #include "TDatabasePDG.h"
42 #include "TLorentzVector.h"
43 
44 #include <iostream>
45 #include <memory>
46 
47 //#define _DEBUG_
48 
49 #define LogDebug(exp) std::cout << "DEBUG: " << __FILE__ << ": " << __LINE__ << ": " << exp << "\n"
50 #define LogError(exp) std::cout << "ERROR: " << __FILE__ << ": " << __LINE__ << ": " << exp << "\n"
51 #define LogWarning(exp) std::cout << "WARNING: " << __FILE__ << ": " << __LINE__ << ": " << exp << "\n"
52 
53 using namespace std;
54 
55 BJetModule::BJetModule(const string &name, const string &out)
56  : SubsysReco(name)
57  , _foutname(out)
58  , _truthjetmap_name("AntiKt_Truth_r04")
59  , _recojetmap_name("AntiKt_Tower_r04")
60  , _trackmap_name("SvtxTrackMap")
61  , _vertexmap_name("SvtxVertexMap")
62  , _embedding_id(1)
63 {
64 }
65 
66 int BJetModule::Init(PHCompositeNode *topNode)
67 {
68  _ievent = 0;
69 
70  _b_event = -1;
71 
72  reset_tree_vars();
73 
74  _f = new TFile(_foutname.c_str(), "RECREATE");
75 
76  _tree = new TTree("T", "a withered deciduous tree");
77 
78  _tree->Branch("event", &_b_event, "event/I");
79  setBranches();
80 
81  return 0;
82 }
83 
84 double calc_dca(const TVector3 &track_point, const TVector3 &track_direction, const TVector3 &vertex)
85 {
86  TVector3 VP = vertex - track_point;
87  double d = -99;
88 
89  if (track_direction.Mag() > 0)
90  {
91  d = (track_direction.Cross(VP)).Mag() / track_direction.Mag();
92  }
93 
94  return d;
95 }
96 
97 bool calc_dca3d_line(
98  double &dca_xy, double &dca_z, /*cm*/
99  const TVector3 &track_point, const TVector3 &track_direction, const TVector3 &vertex)
100 {
101  dca_xy = NAN;
102  dca_z = NAN;
103 
104  if (track_direction.Mag() == 0) return false;
105 
106  TVector3 PV = track_point - vertex;
107  if (PV.Mag() < 0.000001)
108  {
109  dca_xy = 0;
110  dca_z = 0;
111  return true;
112  }
113 
114  TVector3 PVxMom = track_direction.Cross(PV);
115 
116  TVector3 PCA = track_direction;
117 
118  if (PVxMom.Mag() < 0.000001)
119  {
120  std::cout << __FILE__ << ": " << __LINE__ << ": PVxMom.Mag2() < 0.000001" << std::endl;
121  return false;
122  }
123 
124  PCA.Rotate(TMath::PiOver2(), PVxMom); // direction
125 
126  PCA.SetMag(1.);
127 
128  PCA.SetMag(PV.Dot(PCA));
129 
130  TVector3 r = track_direction.Cross(TVector3(0, 0, 1));
131  r.SetMag(1.);
132 
133  dca_xy = PCA.Dot(r);
134 
135  dca_z = PCA.Z();
136 
137  return true;
138 }
139 
140 bool calc_dca_circle_line(
141  double &dca_xy, double &dca_z, /*cm*/
142  const TVector3 &track_point /*cm*/, const TVector3 &track_mom /*GeV/c*/, const TVector3 &vertex, /*cm*/
143  const double &q = 1. /*e*/, const double &B = 1.4 /*Tesla*/)
144 {
145  double z0 = track_point.Z() - vertex.Z();
146  double r0 = (track_point - vertex).Perp();
147  double pz = track_mom.Z();
148  double pt = track_mom.Perp();
149 
150  if (pt == 0 || B == 0 || q == 0)
151  {
152  LogWarning("pt == 0 || B ==0 || q==0");
153  dca_z = NAN;
154  dca_xy = NAN;
155  return false;
156  }
157 
158  dca_z = z0 - pz / pt * r0;
159 
160  const double ten_over_c = 333.564095198; // GeV/c in cm*e*T
161  double R = pt / B / q * ten_over_c; // radius in cm
162 
163  // make 2d versions
164  TVector3 track_point_2d(track_point.X(), track_point.Y(), 0);
165  TVector3 track_mom_2d(track_mom.X(), track_mom.Y(), 0);
166  TVector3 vertex_2d(vertex.X(), vertex.Y(), 0);
167 
168  // calc track circle center
169  TVector3 track_2d_circle_center = track_mom_2d; //copy mom_2d
170  track_2d_circle_center.RotateZ(TMath::PiOver2()); //rotate Pi/2
171  track_2d_circle_center.SetMag(R); // scale to R, R could be negtive
172  track_2d_circle_center += track_point_2d; //shift base to point_2d
173 
174  dca_xy = TMath::Abs(R) - (track_2d_circle_center - vertex_2d).Mag();
175 
176  return true;
177 }
178 
179 int BJetModule::process_event(PHCompositeNode *topNode)
180 {
181  _b_event = _ievent;
182 
183  reset_tree_vars();
184 
185  PHHepMCGenEventMap *geneventmap = findNode::getClass<PHHepMCGenEventMap>(
186  topNode, "PHHepMCGenEventMap");
187  if (!geneventmap)
188  {
189  std::cout << PHWHERE
190  << " - Fatal error - missing node PHHepMCGenEventMap"
191  << std::endl;
192  return Fun4AllReturnCodes::ABORTRUN;
193  }
194 
195  PHHepMCGenEvent *genevt = geneventmap->get(_embedding_id);
196  if (!genevt)
197  {
198  std::cout << PHWHERE
199  << " - Fatal error - node PHHepMCGenEventMap missing subevent with embedding ID "
200  << _embedding_id;
201  std::cout << ". Print PHHepMCGenEventMap:";
202  geneventmap->identify();
203  return Fun4AllReturnCodes::ABORTRUN;
204  }
205 
206  //PHHepMCGenEvent *genevt = findNode::getClass<PHHepMCGenEvent>(topNode,"PHHepMCGenEvent");
207  HepMC::GenEvent *theEvent = genevt->getEvent();
208  //theEvent->print();
209 
210  PHG4TruthInfoContainer *truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
211 
212  PHG4VtxPoint *first_point = truthinfo->GetPrimaryVtx(truthinfo->GetPrimaryVertexIndex());
213  _b_truth_vertex_n = 0;
214  _b_truth_vertex_x[_b_truth_vertex_n] = first_point->get_x();
215  _b_truth_vertex_y[_b_truth_vertex_n] = first_point->get_y();
216  _b_truth_vertex_z[_b_truth_vertex_n] = first_point->get_z();
217  TVector3 truth_primary_vertex(first_point->get_x(), first_point->get_y(), first_point->get_z());
218  ++_b_truth_vertex_n;
219 
220  auto jet_eval_stack = unique_ptr<JetEvalStack>(new JetEvalStack(topNode, _recojetmap_name, _truthjetmap_name));
221  if (!jet_eval_stack)
222  {
223  LogError("!jet_eval_stack");
224  return Fun4AllReturnCodes::ABORTRUN;
225  }
226 
227  JetRecoEval *jet_reco_eval = jet_eval_stack->get_reco_eval();
228  if (!jet_reco_eval)
229  {
230  LogError("!jet_reco_eval");
231  return Fun4AllReturnCodes::ABORTRUN;
232  }
233 
234  JetMap *truth_jets = findNode::getClass<JetMap>(topNode, _truthjetmap_name);
235  JetMap *reco_jets = findNode::getClass<JetMap>(topNode, _recojetmap_name);
236 
237  _b_truthjet_n = 0;
238  for (JetMap::Iter iter = truth_jets->begin(); iter != truth_jets->end();
239  ++iter)
240  {
241  Jet *truth_jet = iter->second;
242  if (truth_jet->get_pt() < 10 || fabs(truth_jet->get_eta()) > 2)
243  continue;
244 
245  _b_truthjet_pt[_b_truthjet_n] = truth_jet->get_pt();
246  _b_truthjet_phi[_b_truthjet_n] = truth_jet->get_phi();
247  _b_truthjet_eta[_b_truthjet_n] = truth_jet->get_eta();
248 
249  int jet_flavor = -999;
250 
251  jet_flavor = truth_jet->get_property(static_cast<Jet::PROPERTY>(prop_JetPartonFlavor));
252  if (abs(jet_flavor) < 100)
253  _b_truthjet_parton_flavor[_b_truthjet_n] = jet_flavor;
254 
255  jet_flavor = truth_jet->get_property(static_cast<Jet::PROPERTY>(prop_JetHadronFlavor));
256  if (abs(jet_flavor) < 100)
257  _b_truthjet_hadron_flavor[_b_truthjet_n] = jet_flavor;
258  //cout << "DEBUG: " << __LINE__ << endl;
259  //auto reco_jet = unique_ptr<Jet>(jet_reco_eval->best_jet_from(truth_jet));
260  Jet *reco_jet = jet_reco_eval->best_jet_from(truth_jet);
261 
262  //cout << "DEBUG: " << __LINE__ << endl;
263  if (!reco_jet)
264  {
266  Jet* matchedjet = nullptr;
267  float mindr = std::numeric_limits<float>::max();
268  for (JetMap::Iter riter = reco_jets->begin(); riter != reco_jets->end();
269  ++riter)
270  {
271  Jet* mjet = riter->second;
272  if(mjet->get_pt() < 5)
273  { continue; }
274  float dr = dR(mjet->get_eta(), truth_jet->get_eta(),
275  mjet->get_phi(), truth_jet->get_phi());
276  if(dr < mindr)
277  {
278  mindr = dr;
279  matchedjet = mjet;
280  }
281  }
282 
283  if(matchedjet)
284  {
285  _b_recojet_valid[_b_truthjet_n] = 1;
286  _b_recojet_pt[_b_truthjet_n] = matchedjet->get_pt();
287  _b_recojet_phi[_b_truthjet_n] = matchedjet->get_phi();
288  _b_recojet_eta[_b_truthjet_n] = matchedjet->get_eta();
289  }
290  else
291  {
292  _b_recojet_valid[_b_truthjet_n] = 0;
293  }
294  }
295  else
296  {
297  _b_recojet_valid[_b_truthjet_n] = 1;
298  _b_recojet_pt[_b_truthjet_n] = reco_jet->get_pt();
299  _b_recojet_phi[_b_truthjet_n] = reco_jet->get_phi();
300  _b_recojet_eta[_b_truthjet_n] = reco_jet->get_eta();
301  }
302 
303  _b_truthjet_n++;
304  //cout << "DEBUG: " << __LINE__ << endl;
305  }
306 
307  PHG4TruthInfoContainer::Range range = truthinfo->GetPrimaryParticleRange();
308  //PHG4TruthInfoContainer::Range range = truthinfo->GetParticleRange();
309 
310  _b_particle_n = 0;
311  for (auto iter = range.first; iter != range.second; ++iter)
312  {
313  PHG4Particle *g4particle = iter->second; // You may ask yourself, why second?
314 
315  unsigned int embed_id = truthinfo->isEmbeded(g4particle->get_track_id());
316 
317  TLorentzVector t;
318  t.SetPxPyPzE(g4particle->get_px(), g4particle->get_py(), g4particle->get_pz(), g4particle->get_e());
319 
320  float truth_pt = t.Pt();
321  if (truth_pt < 0.5) continue;
322  float truth_eta = t.Eta();
323  if (fabs(truth_eta) > 1.1) continue;
324  float truth_phi = t.Phi();
325  int truth_pid = g4particle->get_pid();
326 
327  //original cuts
328  /*
329  if (truth_pid == 22 || truth_pid == 2112 || truth_pid == -2112
330  || truth_pid == 130)
331  continue;
332  if (truth_pid == 2112 || truth_pid == -2112 || truth_pid == 130)
333  continue;
334  // save high-pT photons
335  //if (truth_pid == 22 && truth_pt < 20) continue;
336  if (truth_pid == 12 || truth_pid == -12 || truth_pid == 16
337  || truth_pid == -16 || truth_pid == 14 || truth_pid == -14)
338  continue;
339  */
340 
341  //only keeps stable charged particles
342  if (!(abs(truth_pid) == 211 || abs(truth_pid) == 321 || abs(truth_pid) == 2212 || abs(truth_pid) == 11 || abs(truth_pid) == 13))
343  continue;
344 
345  _b_particle_pid[_b_particle_n] = truth_pid;
346  _b_particle_pt[_b_particle_n] = truth_pt;
347  _b_particle_eta[_b_particle_n] = truth_eta;
348  _b_particle_phi[_b_particle_n] = truth_phi;
349  _b_particle_embed[_b_particle_n] = embed_id;
350 
351  PHG4VtxPoint *point = truthinfo->GetVtx(g4particle->get_vtx_id());
352  _b_particle_vertex_x[_b_particle_n] = point->get_x();
353  _b_particle_vertex_y[_b_particle_n] = point->get_y();
354  _b_particle_vertex_z[_b_particle_n] = point->get_z();
355 
356  TVector3 track_point(point->get_x(), point->get_y(), point->get_z());
357  TVector3 track_mom(g4particle->get_px(), g4particle->get_py(), g4particle->get_pz());
358 
359  double truth_dca_xy = NAN;
360  double truth_dca_z = NAN;
361 
362  //double charge = TDatabasePDG::Instance()->GetParticle(truth_pid)->Charge()/3.;
363  //calc_dca_circle_line(truth_dca_xy, truth_dca_z, track_point, track_mom, truth_primary_vertex, charge, 1.4);
364 
365  calc_dca3d_line(truth_dca_xy, truth_dca_z, track_point, track_mom, truth_primary_vertex);
366 
367 #ifdef _DEBUG_
368  {
369  cout << "track_point: --------------" << endl;
370  track_point.Print();
371 
372  cout << "track_mom: --------------" << endl;
373  track_mom.Print();
374 
375  cout << "truth_primary_vertex: --------------" << endl;
376  truth_primary_vertex.Print();
377 
378  cout
379  << __LINE__
380  << ": " << TDatabasePDG::Instance()->GetParticle(truth_pid)->GetName()
381  << ": truth_dca_xy: " << truth_dca_xy
382  << ": truth_dca_z: " << truth_dca_z
383  << endl
384  << endl;
385  }
386 #endif
387 
388  _b_particle_dca_xy[_b_particle_n] = truth_dca_xy;
389  _b_particle_dca_z[_b_particle_n] = truth_dca_z;
390 
391  // for (HepMC::GenEvent::particle_const_iterator p =
392  // theEvent->particles_begin(); p != theEvent->particles_end();
393  // ++p) {
394  //
395  // if ((*p)->pdg_id() != truth_pid)
396  // continue;
397  //
398  // if (fabs(truth_pt - (*p)->momentum().perp()) > 0.001)
399  // continue;
400  // if (fabs(truth_eta - (*p)->momentum().pseudoRapidity()) > 0.001)
401  // continue;
402  // if (fabs(truth_phi - (*p)->momentum().phi()) > 0.001)
403  // continue;
404  //
405  // HepMC::GenVertex *production_vertex = (*p)->production_vertex();
406  // _b_particle_dca_xy[_b_particle_n] = production_vertex->point3d().perp();
407  //
408  // {
409  // cout
410  // << __LINE__
411  // <<": From HepMC: {"
412  // << production_vertex->point3d().x() <<", "
413  // << production_vertex->point3d().y() <<", "
414  // << production_vertex->point3d().z() <<"}"
415  // <<endl<<endl;
416  // }
417  // }
418 
419  _b_particle_n++;
420  }
421 
422  SvtxTrackMap *trackmap = findNode::getClass<SvtxTrackMap>(topNode, _trackmap_name.c_str());
423 
424  SvtxVertexMap *vertexmap = findNode::getClass<SvtxVertexMap>(topNode, _vertexmap_name.c_str());
425 
426  auto svtxevalstack = unique_ptr<SvtxEvalStack>(new SvtxEvalStack(topNode));
427  if (!svtxevalstack)
428  {
429  LogError("!svtxevalstack");
430  return Fun4AllReturnCodes::ABORTRUN;
431  }
432 
433  svtxevalstack->next_event(topNode);
434  SvtxTrackEval *trackeval = svtxevalstack->get_track_eval();
435 
436  float vertex_x = -99;
437  float vertex_y = -99;
438  float vertex_z = -99;
439  // float vertex_err_x = -99;
440  // float vertex_err_y = -99;
441  //float vertex_err_z = -99;
442 
443  _b_track_n = 0;
444 
445  for (SvtxTrackMap::Iter iter = trackmap->begin(); iter != trackmap->end();
446  ++iter)
447  {
448  SvtxTrack *track = iter->second;
449 
450  float track_pt = track->get_pt();
451  float track_eta = track->get_eta();
452  float track_phi = track->get_phi();
453 
454  if (track_pt < 0.5)
455  continue;
456  if (fabs(track_eta) > 1.1)
457  continue;
458 
459  //std::set<PHG4Hit*> assoc_hits = trackeval->all_truth_hits(track);//TODO
460  int nmaps = 0;
461 
462  unsigned int nclusters = track->size_cluster_keys();
463  unsigned int nclusters_by_layer = 0;
464 
465  TrackSeed *siliconSeed = track->get_silicon_seed();
466  for (auto iter = siliconSeed->begin_cluster_keys();
467  iter != siliconSeed->end_cluster_keys(); ++iter)
468  {
469  TrkrDefs::cluskey key = *iter;
470  unsigned int cluster_layer = TrkrDefs::getLayer(key);
471 
472  if (TrkrDefs::getTrkrId(key) == TrkrDefs::TrkrId::mvtxId)
473  {
474  ++nmaps;
475  }
476  nclusters_by_layer |= (0x3FFFFFFF & (0x1 << cluster_layer));
477  }
478 
479  TrackSeed *tpcSeed = track->get_silicon_seed();
480  for (auto iter = tpcSeed->begin_cluster_keys();
481  iter != tpcSeed->end_cluster_keys(); ++iter)
482  {
483  unsigned int cluster_layer = TrkrDefs::getLayer(*iter);
484  nclusters_by_layer |= (0x3FFFFFFF & (0x1 << cluster_layer));
485  }
486 
487  PHG4Particle *g4particle = trackeval->max_truth_particle_by_nclusters(
488  track);
489 
490  unsigned int truth_nclusters = trackeval->get_nclusters_contribution(
491  track, g4particle);
492  unsigned int truth_nclusters_by_layer = trackeval->get_nclusters_contribution_by_layer(
493  track, g4particle);
494 
495  unsigned int truth_embed_id = truthinfo->isEmbeded(
496  g4particle->get_track_id());
497  bool truth_is_primary = truthinfo->is_primary(g4particle);
498  int truth_pid = g4particle->get_pid();
499 
500  //TVector3 g4particle_mom(g4particle->get_px(),g4particle->get_py(),g4particle->get_pz());
501  //_b_track_best_pt[_b_track_n] = g4particle_mom.Pt();
502 
503  //int truth_parent_id = g4particle->get_parent_id();
504  //int truth_primary_id = g4particle->get_primary_id();
505 
506  TVector3 track_point(track->get_x(), track->get_y(), track->get_z());
507  TVector3 track_direction(track->get_px(), track->get_py(), track->get_pz());
508  SvtxVertex* svertex = vertexmap->get(track->get_vertex_id());
509  vertex_x = svertex->get_x();
510  vertex_y = svertex->get_y();
511  vertex_z = svertex->get_z();
512  TVector3 vertex(vertex_x, vertex_y, vertex_z);
513 
514  TVector3 track_point_2d(track->get_x(), track->get_y(), 0);
515  TVector3 track_direction_2d(track->get_px(), track->get_py(), 0);
516  TVector3 vertex_2d(vertex_x, vertex_y, 0);
517  float dca3d_xy = NAN;
518  float dca3d_xy_error = NAN;
519  float dca3d_z = NAN;
520  float dca3d_z_error = NAN;
521  calc3DDCA(track, vertexmap, dca3d_xy, dca3d_xy_error, dca3d_z, dca3d_z_error);
522 
523  float dca2d_calc = calc_dca(track_point_2d, track_direction_2d, vertex_2d);
524  float dca2d_calc_truth = calc_dca(track_point_2d, track_direction_2d, TVector3(0, 0, 0));
525 
526  float dca3d_calc = calc_dca(track_point, track_direction, vertex);
527  float dca3d_calc_truth = calc_dca(track_point, track_direction, TVector3(0, 0, 0));
528 
529  // float dca2d_calc_error = sqrt(
530  // dca2d_error*dca2d_error +
531  // vertex_err_x*vertex_err_x +
532  // vertex_err_y*vertex_err_y);
533 
534  _b_track_pca_x[_b_track_n] = track->get_x() - vertex_x;
535  _b_track_pca_y[_b_track_n] = track->get_y() - vertex_y;
536  _b_track_pca_z[_b_track_n] = track->get_z() - vertex_z;
537 
538  TVector3 dca_vector(
539  _b_track_pca_x[_b_track_n],
540  _b_track_pca_y[_b_track_n],
541  _b_track_pca_z[_b_track_n]);
542 
543  _b_track_pca_phi[_b_track_n] = dca_vector.Phi();
544 
545  TLorentzVector t;
546  t.SetPxPyPzE(g4particle->get_px(), g4particle->get_py(),
547  g4particle->get_pz(), g4particle->get_e());
548 
549  float truth_pt = t.Pt();
550  float truth_eta = t.Eta();
551  float truth_phi = t.Phi();
552 
553  PHG4VtxPoint *point = truthinfo->GetVtx(g4particle->get_vtx_id());
554  _b_track_best_vertex_x[_b_track_n] = point->get_x();
555  _b_track_best_vertex_y[_b_track_n] = point->get_y();
556  _b_track_best_vertex_z[_b_track_n] = point->get_z();
557 
558  TVector3 track_best_point(point->get_x(), point->get_y(), point->get_z());
559  TVector3 track_best_mom(g4particle->get_px(), g4particle->get_py(), g4particle->get_pz());
560 
561  double truth_dca_xy = NAN;
562  double truth_dca_z = NAN;
563 
564  //double charge = TDatabasePDG::Instance()->GetParticle(truth_pid)->Charge()/3.;
565  //calc_dca_circle_line(truth_dca_xy, truth_dca_z, track_best_point, track_best_mom, truth_primary_vertex, charge, 1.4);
566 
567  calc_dca3d_line(truth_dca_xy, truth_dca_z, track_best_point, track_best_mom, truth_primary_vertex);
568 
569  _b_track_best_dca_xy[_b_track_n] = truth_dca_xy;
570  _b_track_best_dca_z[_b_track_n] = truth_dca_z;
571 
572  //_b_track_best_parent_pid[_b_track_n] = truthinfo->GetParticle(g4particle->get_parent_id())->get_pid();
573  int truth_in = -1;
574  int truth_out = -1;
575  int nhepmc = 0;
576  for (HepMC::GenEvent::particle_const_iterator p =
577  theEvent->particles_begin();
578  p != theEvent->particles_end();
579  ++p)
580  {
581  if ((*p)->pdg_id() != truth_pid)
582  continue;
583 
584  if (fabs(truth_pt - (*p)->momentum().perp()) > 0.001)
585  continue;
586  if (fabs(truth_eta - (*p)->momentum().pseudoRapidity()) > 0.001)
587  continue;
588  if (fabs(truth_phi - (*p)->momentum().phi()) > 0.001)
589  continue;
590 
591  HepMC::GenVertex *production_vertex = (*p)->production_vertex();
592 
593  truth_in = production_vertex->particles_in_size();
594  truth_out = production_vertex->particles_out_size();
595 
596  HepMC::GenVertex::particles_in_const_iterator first_parent =
597  production_vertex->particles_in_const_begin();
598  _b_track_best_parent_pid[_b_track_n] = (*first_parent)->pdg_id();
599 
600  nhepmc++;
601  }
602 
603  _b_track_pt[_b_track_n] = track_pt;
604  _b_track_eta[_b_track_n] = track_eta;
605  _b_track_phi[_b_track_n] = track_phi;
606 
607  _b_track_dca3d_xy[_b_track_n] = dca3d_xy;
608  _b_track_dca3d_xy_error[_b_track_n] = dca3d_xy_error;
609 
610  _b_track_dca3d_z[_b_track_n] = dca3d_z;
611  _b_track_dca3d_z_error[_b_track_n] = dca3d_z_error;
612 
613  _b_track_dca2d_calc[_b_track_n] = dca2d_calc;
614  _b_track_dca2d_calc_truth[_b_track_n] = dca2d_calc_truth;
615 
616  _b_track_dca3d_calc[_b_track_n] = dca3d_calc;
617  _b_track_dca3d_calc_truth[_b_track_n] = dca3d_calc_truth;
618 
619  // _b_track_pca_phi[_b_track_n] = dca_phi;
620  // _b_track_pca_x[_b_track_n] = dca_x;
621  // _b_track_pca_y[_b_track_n] = dca_y;
622  // _b_track_pca_z[_b_track_n] = dca_z;
623 
624  _b_track_quality[_b_track_n] = track->get_quality();
625  _b_track_chisq[_b_track_n] = track->get_chisq();
626  _b_track_ndf[_b_track_n] = track->get_ndf();
627 
628  _b_track_nmaps[_b_track_n] = nmaps;
629 
630  _b_track_nclusters[_b_track_n] = nclusters;
631  _b_track_nclusters_by_layer[_b_track_n] = nclusters_by_layer;
632 
633  _b_track_best_nclusters[_b_track_n] = truth_nclusters;
634  _b_track_best_nclusters_by_layer[_b_track_n] = truth_nclusters_by_layer;
635  _b_track_best_embed[_b_track_n] = truth_embed_id;
636  _b_track_best_primary[_b_track_n] = truth_is_primary;
637  _b_track_best_pid[_b_track_n] = truth_pid;
638  _b_track_best_pt[_b_track_n] = truth_pt;
639 
640  _b_track_best_in[_b_track_n] = truth_in;
641  _b_track_best_out[_b_track_n] = truth_out;
642 
643  _b_track_n++;
644  }
645 
646  _tree->Fill();
647 
648  _ievent++;
649 
650  //delete jet_eval_stack;
651  //delete svtxevalstack;
652 
653  return 0;
654 }
655 
656 void BJetModule::calc3DDCA(SvtxTrack* track, SvtxVertexMap* vertexmap,
657  float& dca3d_xy, float& dca3d_xy_error,
658  float& dca3d_z, float& dca3d_z_error)
659 {
660  Acts::Vector3 pos(track->get_x(),
661  track->get_y(),
662  track->get_z());
663  Acts::Vector3 mom(track->get_px(),
664  track->get_py(),
665  track->get_pz());
666 
667  auto vtxid = track->get_vertex_id();
668  auto svtxVertex = vertexmap->get(vtxid);
669  Acts::Vector3 vertex(svtxVertex->get_x(),
670  svtxVertex->get_y(),
671  svtxVertex->get_z());
672 
673  pos -= vertex;
674 
675  Acts::ActsSymMatrix<3> posCov;
676  for(int i = 0; i < 3; ++i)
677  {
678  for(int j = 0; j < 3; ++j)
679  {
680  posCov(i, j) = track->get_error(i, j);
681  }
682  }
683 
684  Acts::Vector3 r = mom.cross(Acts::Vector3(0.,0.,1.));
685  float phi = atan2(r(1), r(0));
686 
687  Acts::RotationMatrix3 rot;
688  Acts::RotationMatrix3 rot_T;
689  rot(0,0) = cos(phi);
690  rot(0,1) = -sin(phi);
691  rot(0,2) = 0;
692  rot(1,0) = sin(phi);
693  rot(1,1) = cos(phi);
694  rot(1,2) = 0;
695  rot(2,0) = 0;
696  rot(2,1) = 0;
697  rot(2,2) = 1;
698 
699  rot_T = rot.transpose();
700 
701  Acts::Vector3 pos_R = rot * pos;
702  Acts::ActsSymMatrix<3> rotCov = rot * posCov * rot_T;
703 
704  dca3d_xy= pos_R(0);
705  dca3d_z = pos_R(2);
706  dca3d_xy_error = sqrt(rotCov(0,0));
707  dca3d_z_error = sqrt(rotCov(2,2));
708 }
709 
710 int BJetModule::reset_tree_vars()
711 {
712  _b_truthjet_n = 0;
713 
714  for (int n = 0; n < 10; n++)
715  {
716  _b_truthjet_parton_flavor[n] = -9999;
717  _b_truthjet_hadron_flavor[n] = -9999;
718 
719  _b_truthjet_pt[n] = -99;
720  _b_truthjet_eta[n] = -99;
721  _b_truthjet_phi[n] = -99;
722 
723  _b_recojet_valid[n] = 0;
724  _b_truthjet_pt[n] = -99;
725  _b_truthjet_eta[n] = -99;
726  _b_truthjet_phi[n] = -99;
727  }
728 
729  _b_particle_n = 0;
730  _b_track_n = 0;
731 
732  for (int n = 0; n < 1000; n++)
733  {
734  _b_particle_pt[n] = -1;
735  _b_particle_eta[n] = -1;
736  _b_particle_phi[n] = -1;
737  _b_particle_pid[n] = 0;
738  _b_particle_embed[n] = 0;
739 
740  _b_particle_vertex_x[n] = -1;
741  _b_particle_vertex_y[n] = -1;
742  _b_particle_vertex_z[n] = -1;
743 
744  _b_particle_dca_xy[n] = -1;
745  _b_particle_dca_z[n] = -1;
746 
747  _b_track_pt[n] = -1;
748  _b_track_eta[n] = -1;
749  _b_track_phi[n] = -1;
750 
751  _b_track_nmaps[n] = 0;
752 
753  _b_track_nclusters[n] = 0;
754  _b_track_nclusters_by_layer[n] = 0;
755 
756  _b_track_dca2d[n] = -1;
757  _b_track_dca2d_error[n] = -1;
758  _b_track_dca3d_xy[n] = -1;
759  _b_track_dca3d_xy_error[n] = -1;
760  _b_track_dca3d_z[n] = -1;
761  _b_track_dca3d_z_error[n] = -1;
762  _b_track_dca2d_calc[n] = -1;
763  _b_track_dca2d_calc_truth[n] = -1;
764  _b_track_dca3d_calc[n] = -1;
765  _b_track_dca3d_calc_truth[n] = -1;
766 
767  _b_track_pca_phi[n] = -99;
768  _b_track_pca_x[n] = -99;
769  _b_track_pca_y[n] = -99;
770  _b_track_pca_z[n] = -99;
771 
772  _b_track_quality[n] = -99;
773  _b_track_chisq[n] = -99;
774  _b_track_ndf[n] = -99;
775 
776  _b_track_best_nclusters[n] = 0;
777  _b_track_best_nclusters_by_layer[n] = 0;
778  _b_track_best_embed[n] = 0;
779  _b_track_best_primary[n] = false;
780  _b_track_best_pid[n] = 0;
781 
782  _b_track_best_in[n] = 0;
783  _b_track_best_out[n] = 0;
784  _b_track_best_parent_pid[n] = 0;
785 
786  _b_track_best_vertex_x[n] = -99;
787  _b_track_best_vertex_y[n] = -99;
788  _b_track_best_vertex_z[n] = -99;
789  _b_track_best_dca_xy[n] = -99;
790  _b_track_best_dca_z[n] = -99;
791  }
792 
793  return 0;
794 }
795 
796 
797 void BJetModule::setBranches()
798 {
799 
800  _tree->Branch("truth_vertex_n", &_b_truth_vertex_n, "truth_vertex_n/I");
801  _tree->Branch("truth_vertex_x", _b_truth_vertex_x, "truth_vertex_x[truth_vertex_n]/F");
802  _tree->Branch("truth_vertex_y", _b_truth_vertex_y, "truth_vertex_y[truth_vertex_n]/F");
803  _tree->Branch("truth_vertex_z", _b_truth_vertex_z, "truth_vertex_z[truth_vertex_n]/F");
804 
805  _tree->Branch("truthjet_n", &_b_truthjet_n, "truthjet_n/I");
806  _tree->Branch("truthjet_parton_flavor", _b_truthjet_parton_flavor, "truthjet_parton_flavor[truthjet_n]/I");
807  _tree->Branch("truthjet_hadron_flavor", _b_truthjet_hadron_flavor, "truthjet_hadron_flavor[truthjet_n]/I");
808  _tree->Branch("truthjet_pt", _b_truthjet_pt, "truthjet_pt[truthjet_n]/F");
809  _tree->Branch("truthjet_eta", _b_truthjet_eta, "truthjet_eta[truthjet_n]/F");
810  _tree->Branch("truthjet_phi", _b_truthjet_phi, "truthjet_phi[truthjet_n]/F");
811 
812  _tree->Branch("recojet_valid", _b_recojet_valid, "recojet_valid[truthjet_n]/I");
813  _tree->Branch("recojet_pt", _b_recojet_pt, "recojet_pt[truthjet_n]/F");
814  _tree->Branch("recojet_eta", _b_recojet_eta, "recojet_eta[truthjet_n]/F");
815  _tree->Branch("recojet_phi", _b_recojet_phi, "recojet_phi[truthjet_n]/F");
816 
817  _tree->Branch("particle_n", &_b_particle_n, "particle_n/I");
818  _tree->Branch("particle_pt", _b_particle_pt, "particle_pt[particle_n]/F");
819  _tree->Branch("particle_eta", _b_particle_eta, "particle_eta[particle_n]/F");
820  _tree->Branch("particle_phi", _b_particle_phi, "particle_phi[particle_n]/F");
821  _tree->Branch("particle_pid", _b_particle_pid, "particle_pid[particle_n]/I");
822  _tree->Branch("particle_embed", _b_particle_embed, "particle_embed[particle_n]/i");
823 
824  _tree->Branch("particle_vertex_x", _b_particle_vertex_x, "particle_vertex_x[particle_n]/F");
825  _tree->Branch("particle_vertex_y", _b_particle_vertex_y, "particle_vertex_y[particle_n]/F");
826  _tree->Branch("particle_vertex_z", _b_particle_vertex_z, "particle_vertex_z[particle_n]/F");
827  _tree->Branch("particle_dca_xy", _b_particle_dca_xy, "particle_dca_xy[particle_n]/F");
828  _tree->Branch("particle_dca_z", _b_particle_dca_z, "particle_dca_z[particle_n]/F");
829 
830  _tree->Branch("track_n", &_b_track_n, "track_n/I");
831  _tree->Branch("track_pt", _b_track_pt, "track_pt[track_n]/F");
832  _tree->Branch("track_eta", _b_track_eta, "track_eta[track_n]/F");
833  _tree->Branch("track_phi", _b_track_phi, "track_phi[track_n]/F");
834 
835  _tree->Branch("track_dca2d", _b_track_dca2d, "track_dca2d[track_n]/F");
836  _tree->Branch("track_dca2d_error", _b_track_dca2d_error, "track_dca2d_error[track_n]/F");
837 
838  _tree->Branch("track_dca3d_xy", _b_track_dca3d_xy, "track_dca3d_xy[track_n]/F");
839  _tree->Branch("track_dca3d_xy_error", _b_track_dca3d_xy_error, "track_dca3d_xy_error[track_n]/F");
840 
841  _tree->Branch("track_dca3d_z", _b_track_dca3d_z, "track_dca3d_z[track_n]/F");
842  _tree->Branch("track_dca3d_z_error", _b_track_dca3d_z_error, "track_dca3d_z_error[track_n]/F");
843 
844  _tree->Branch("track_dca2d_calc", _b_track_dca2d_calc, "track_dca2d_calc[track_n]/F");
845  _tree->Branch("track_dca2d_calc_truth", _b_track_dca2d_calc_truth, "track_dca2d_calc_truth[track_n]/F");
846 
847  _tree->Branch("track_dca3d_calc", _b_track_dca3d_calc, "track_dca3d_calc[track_n]/F");
848  _tree->Branch("track_dca3d_calc_truth", _b_track_dca3d_calc_truth, "track_dca3d_calc_truth[track_n]/F");
849 
850  _tree->Branch("track_pca_phi", _b_track_pca_phi, "track_pca_phi[track_n]/F");
851  _tree->Branch("track_pca_x", _b_track_pca_x, "track_pca_x[track_n]/F");
852  _tree->Branch("track_pca_y", _b_track_pca_y, "track_pca_y[track_n]/F");
853  _tree->Branch("track_pca_z", _b_track_pca_z, "track_pca_z[track_n]/F");
854 
855  _tree->Branch("track_quality", _b_track_quality, "track_quality[track_n]/F");
856  _tree->Branch("track_chisq", _b_track_chisq, "track_chisq[track_n]/F");
857  _tree->Branch("track_ndf", _b_track_ndf, "track_ndf[track_n]/I");
858 
859  _tree->Branch("track_nmaps", _b_track_nmaps, "track_nmaps[track_n]/I");
860 
861  _tree->Branch("track_nclusters", _b_track_nclusters, "track_nclusters[track_n]/i");
862  _tree->Branch("track_nclusters_by_layer", _b_track_nclusters_by_layer, "track_nclusters_by_layer[track_n]/i");
863 
864  _tree->Branch("track_best_nclusters", _b_track_best_nclusters, "track_best_nclusters[track_n]/i");
865  _tree->Branch("track_best_nclusters_by_layer", _b_track_best_nclusters_by_layer, "track_best_nclusters_by_layer[track_n]/i");
866 
867  _tree->Branch("track_best_embed", _b_track_best_embed, "track_best_embed[track_n]/i");
868  _tree->Branch("track_best_primary", _b_track_best_primary, "track_best_primary[track_n]/O");
869  _tree->Branch("track_best_pid", _b_track_best_pid, "track_best_pid[track_n]/I");
870  _tree->Branch("track_best_pt", _b_track_best_pt, "track_best_pt[track_n]/F");
871  _tree->Branch("track_best_in", _b_track_best_in, "track_best_in[track_n]/I");
872  _tree->Branch("track_best_out", _b_track_best_out, "track_best_out[track_n]/I");
873  _tree->Branch("track_best_parent_pid", _b_track_best_parent_pid, "track_best_parent_pid[track_n]/I");
874 
875  _tree->Branch("track_best_vertex_x", _b_track_best_vertex_x, "track_best_vertex_x[track_n]/F");
876  _tree->Branch("track_best_vertex_y", _b_track_best_vertex_y, "track_best_vertex_y[track_n]/F");
877  _tree->Branch("track_best_vertex_z", _b_track_best_vertex_z, "track_best_vertex_z[track_n]/F");
878 
879  _tree->Branch("track_best_dca_xy", _b_track_best_dca_xy, "track_best_dca_xy[track_n]/F");
880  _tree->Branch("track_best_dca_z", _b_track_best_dca_z, "track_best_dca_z[track_n]/F");
881 
882 }
883 int BJetModule::End(PHCompositeNode *topNode)
884 {
885  //_f->ls();
886  //_tree->Write();
887  _f->Write();
888  _f->Close();
889 
890  return 0;
891 }