Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
BDiJetModule.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file BDiJetModule.C
1 
8 #include "BDiJetModule.h"
9 
10 #include <phool/getClass.h>
11 #include <phool/phool.h>
12 #include <phool/PHCompositeNode.h>
13 #include <phool/PHIODataNode.h>
14 #include <phool/PHNodeIterator.h>
15 
16 #include <phgeom/PHGeomUtility.h>
17 #include <phfield/PHFieldUtility.h>
18 
19 #include <fun4all/Fun4AllServer.h>
20 #include <fun4all/Fun4AllReturnCodes.h>
21 #include <fun4all/PHTFileServer.h>
22 
23 #include <g4main/PHG4TruthInfoContainer.h>
24 #include <g4main/PHG4Particle.h>
25 #include <g4main/PHG4Hit.h>
26 #include <g4main/PHG4HitContainer.h>
27 
28 #include <g4hough/SvtxCluster.h>
29 #include <g4hough/SvtxClusterMap.h>
30 #include <g4hough/SvtxHitMap.h>
31 #include <g4hough/SvtxTrack.h>
32 #include <g4hough/SvtxVertex.h>
33 #include <g4hough/SvtxTrackMap.h>
34 #include <g4hough/SvtxVertexMap.h>
35 
36 #include <g4jets/JetMap.h>
37 #include <g4jets/Jet.h>
38 
39 #include <g4detectors/PHG4CellContainer.h>
40 #include <g4detectors/PHG4CylinderGeomContainer.h>
41 #include <g4detectors/PHG4Cell.h>
42 #include <g4detectors/PHG4CylinderGeom_MAPS.h>
43 #include <g4detectors/PHG4CylinderGeom_Siladders.h>
44 
45 #include <phgenfit/Fitter.h>
46 #include <phgenfit/PlanarMeasurement.h>
47 #include <phgenfit/Track.h>
48 #include <phgenfit/SpacepointMeasurement.h>
49 
50 //GenFit
51 #include <GenFit/FieldManager.h>
52 #include <GenFit/GFRaveConverters.h>
53 #include <GenFit/GFRaveVertex.h>
54 #include <GenFit/GFRaveVertexFactory.h>
55 #include <GenFit/MeasuredStateOnPlane.h>
56 #include <GenFit/RKTrackRep.h>
57 #include <GenFit/StateOnPlane.h>
58 #include <GenFit/Track.h>
59 
60 //Rave
61 #include <rave/Version.h>
62 #include <rave/Track.h>
63 #include <rave/VertexFactory.h>
64 #include <rave/ConstantMagneticField.h>
65 
66 #include <phhepmc/PHHepMCGenEvent.h>
67 #include <HepMC/GenEvent.h>
68 #include <HepMC/GenVertex.h>
69 
70 #include <HFJetTruthGeneration/HFJetDefs.h>
71 
72 #include "TTree.h"
73 #include "TFile.h"
74 #include "TLorentzVector.h"
75 
76 
77 #include <iostream>
78 #include <map>
79 #include <utility>
80 #include <vector>
81 #include <memory>
82 
83 
84 
85 #define LogDebug(exp) std::cout<<"DEBUG: " <<__FILE__<<": "<<__LINE__<<": "<< exp <<"\n"
86 #define LogError(exp) std::cout<<"ERROR: " <<__FILE__<<": "<<__LINE__<<": "<< exp <<"\n"
87 #define LogWarning(exp) std::cout<<"WARNING: " <<__FILE__<<": "<<__LINE__<<": "<< exp <<"\n"
88 
89 
90 BDiJetModule::BDiJetModule(const std::string &name,
91  const std::string &ofName) :
92  SubsysReco(name),
93  _jetmap_truth(NULL),
94  _clustermap(NULL),
95  _trackmap(NULL),
96  _vertexmap(NULL),
97  _verbose(false),
98  _write_tree(true),
99  _ana_truth(true),
100  _ana_reco(true),
101  _do_evt_display(false),
102  _use_ladder_geom(false),
103  _cut_jet(true),
104  _cut_Ncluster(false),
105  _foutname(ofName),
106  _f(NULL),
107  _tree(NULL),
108  _primary_pid_guess(211),
109  _cut_min_pT(0.1),
110  _cut_chi2_ndf(5.0),
111  _cut_jet_pT(10.0),
112  _cut_jet_eta(0.7),
113  _cut_jet_R(0.4),
114  _track_fitting_alg_name("DafRef"),
115  _fitter(NULL),
116  _vertexing_method("avf-smoothing:1"),
117  _vertex_finder(NULL)
118 {
119 
120 
121 }
122 
123 int BDiJetModule::Init(PHCompositeNode *topNode)
124 {
125 
126  //-- set up the tree
127  if ( _write_tree )
128  {
129  ResetVariables();
130  InitTree();
131  }
132 
133  return 0;
134 
135 }
136 
137 int BDiJetModule::InitRun(PHCompositeNode *topNode)
138 {
139 
140  TGeoManager* tgeo_manager = PHGeomUtility::GetTGeoManager(topNode);
141  PHField * field = PHFieldUtility::GetFieldMapNode(nullptr, topNode);
142 
143  _fitter = PHGenFit::Fitter::getInstance(tgeo_manager, field,
144  _track_fitting_alg_name, "RKTrackRep", _do_evt_display);
145 
146  if (!_fitter) {
147  std::cerr << PHWHERE << std::endl;
148  return Fun4AllReturnCodes::ABORTRUN;
149  }
150 
151  _vertex_finder = new genfit::GFRaveVertexFactory(verbosity);
152  _vertex_finder->setMethod(_vertexing_method.data());
153 
154  if (!_vertex_finder) {
155  std::cerr << PHWHERE << std::endl;
156  return Fun4AllReturnCodes::ABORTRUN;
157  }
158 
159 
160  //-- set counters
161  _ievent = -1; // since we count at the beginning of the event, start at -1
162  _b_event = -1;
163 
164 
165  return 0;
166 
167 }
168 
169 int BDiJetModule::process_event(PHCompositeNode *topNode)
170 {
171 
172  // count event
173  _ievent++;
174  _b_event = _ievent;
175 
176  //-- reset tree variables
177  if ( _write_tree )
178  ResetVariables();
179 
180  //------
181  //-- Get the nodes
182  //------
183  int retval = GetNodes(topNode);
184  if ( retval != Fun4AllReturnCodes::EVENT_OK )
185  return retval;
186 
187  //------
188  //-- Analyze truth jets
189  //------
190  if ( _ana_truth )
191  {
192  _truthjet_n = 0;
193 
194  int ijet_t = 0;
195  for (JetMap::Iter iter = _jetmap_truth->begin(); iter != _jetmap_truth->end(); ++iter)
196  {
197  Jet* this_jet = iter->second;
198 
199  float this_pt = this_jet->get_pt();
200  float this_phi = this_jet->get_phi();
201  float this_eta = this_jet->get_eta();
202 
203  // this_jet->get_property(static_cast<Jet::PROPERTY>(prop_JetPartonFlavor)) != 5)
204  if (this_jet->get_pt() < _cut_jet_pT || fabs(this_eta) > _cut_jet_eta)
205  continue;
206 
207  if ( _write_tree )
208  {
209  _truthjet_parton_flavor[_truthjet_n] = this_jet->get_property(static_cast<Jet::PROPERTY>(prop_JetPartonFlavor));
210  _truthjet_hadron_flavor[_truthjet_n] = this_jet->get_property(static_cast<Jet::PROPERTY>(prop_JetHadronFlavor));
211  _truthjet_pt[_truthjet_n] = this_pt;
212  _truthjet_phi[_truthjet_n] = this_phi;
213  _truthjet_eta[_truthjet_n] = this_eta;
214  }
215 
216  if (_verbose)
217  std::cout << " truth jet #" << ijet_t << ", pt / eta / phi = "
218  << this_pt << " / " << this_eta << " / " << this_phi
219  << std::endl;
220 
221  _truthjet_n++;
222  }
223 
224  // only care if the event has a jet of interest
225  if ( _truthjet_n < 1 )
226  return Fun4AllReturnCodes::ABORTEVENT;
227 
228  } // _ana_truth
229 
230  //------
231  //-- Analyze reconstructed information
232  //------
233  if ( _ana_reco )
234  {
235 
237  std::vector<genfit::Track*> rf_gf_tracks;
238  rf_gf_tracks.clear();
239  std::vector<PHGenFit::Track*> rf_phgf_tracks;
240  rf_phgf_tracks.clear();
241 
242  std::map<unsigned int, unsigned int> svtxtrk_gftrk_map;
243  std::map<unsigned int, unsigned int> svtxtrk_nmvtx_map;
244  std::map<unsigned int, unsigned int> svtxtrk_nintt_map;
245 
246  SvtxVertex *vertex = _vertexmap->get(0);
247 
248  for (SvtxTrackMap::Iter iter = _trackmap->begin();
249  iter != _trackmap->end();
250  ++iter)
251  {
252  SvtxTrack* svtx_track = iter->second;
253  if (!svtx_track)
254  continue;
255  if (!(svtx_track->get_pt() > _cut_min_pT))
256  continue;
257  if ((svtx_track->get_chisq() / svtx_track->get_ndf()) > _cut_chi2_ndf)
258  continue;
259 
260  int n_MVTX = 0, n_INTT = 0;
261  for (SvtxTrack::ConstClusterIter iter2 = svtx_track->begin_clusters();
262  iter2 != svtx_track->end_clusters();
263  iter2++)
264  {
265 
266  SvtxCluster* cluster = _clustermap->get(*iter2);
267  unsigned int layer = cluster->get_layer();
268 
269  if (layer < 3) n_MVTX++;
270  else if (layer < 7) n_INTT++;
271  } // iter2 over SvtxTrack::ConstClusterIter
272 
273  //cout << "n MVTX: " << n_MVTX << ", n INTT: " << n_INTT << std::endl;
274 
275  /*
276  //if (n_MVTX<2 || n_INTT<2)
277  if (n_MVTX<2)
278  continue;
279  */
280 
281  PHGenFit::Track* rf_phgf_track = MakeGenFitTrack(topNode, svtx_track, vertex);
282 
283  if (rf_phgf_track)
284  {
285  svtxtrk_gftrk_map.insert(std::pair<int, int>(svtx_track->get_id(), rf_phgf_tracks.size()));
286  svtxtrk_nmvtx_map.insert(std::pair<int, int>(svtx_track->get_id(), n_MVTX));
287  svtxtrk_nintt_map.insert(std::pair<int, int>(svtx_track->get_id(), n_INTT));
288  rf_phgf_tracks.push_back(rf_phgf_track);
289 
290  if ( _cut_Ncluster )
291  {
292  if ( n_MVTX >= 1 && n_INTT >= 2 )
293  rf_gf_tracks.push_back(rf_phgf_track->getGenFitTrack());
294  }
295  else
296  {
297  rf_gf_tracks.push_back(rf_phgf_track->getGenFitTrack());
298  }
299 
300  }
301  } // iter over SvtxTrackMap
302 
304  std::vector<genfit::GFRaveVertex*> rave_vertices;
305  rave_vertices.clear();
307 
308  if ( _vertexing_method.compare("avr-smoothing:1") != 0 )
309  _vertex_finder->setMethod(_vertexing_method.data());
310  if (rf_gf_tracks.size() >= 2)
311  {
312  try {
313  _vertex_finder->findVertices(&rave_vertices, rf_gf_tracks);
314  } catch (...) {
315  std::cout << PHWHERE << "GFRaveVertexFactory::findVertices failed!";
316  }
317  }
318 
319  FillVertexMap(rave_vertices, rf_gf_tracks);
320 
322  if ( _vertexing_method.compare("avr-smoothing:1") != 0 )
323  _vertex_finder->setMethod("avr-smoothing:1");
324  std::vector<genfit::GFRaveVertex*> rave_vertices_jet;
325  rave_vertices_jet.clear();
326 
328  for (JetMap::ConstIter iter = _jetmap_truth->begin(); iter != _jetmap_truth->end(); iter++)
329  {
330  Jet *jet = iter->second;
331 
332  float jet_pT = jet->get_pt();
333  float jet_eta = jet->get_eta();
334  float jet_phi = jet->get_phi();
335 
336  if ( jet_pT < _cut_jet_pT ) continue;
337  if ( fabs(jet_eta) > _cut_jet_eta ) continue;
338 
339  int counter_r10 = 0, counter_miss = 0;
340 
341  std::vector<genfit::Track*> rf_gf_tracks_jet;
342  rf_gf_tracks_jet.clear();
343  // vector<genfit::Track*> rf_gf_tracks_jet_pT05;
344  // rf_gf_tracks_jet_pT05.clear();
345  // vector<genfit::Track*> rf_gf_tracks_jet_pT10;
346  // rf_gf_tracks_jet_pT10.clear();
347  // vector<genfit::Track*> rf_gf_tracks_jet_pT15;
348  // rf_gf_tracks_jet_pT15.clear();
349  // vector<genfit::Track*> rf_gf_tracks_jet_pT20;
350  // rf_gf_tracks_jet_pT20.clear();
351 
352  //cout << "JET ETA: " << jet_eta << ", JET PHI: " << jet_phi << std::endl;
353 
354  for (SvtxTrackMap::ConstIter iter2 = _trackmap->begin(); iter2 != _trackmap->end(); iter2++)
355  {
356  SvtxTrack* svtx_track = iter2->second;
357 
358  float trk_phi = svtx_track->get_phi();
359  float trk_eta = svtx_track->get_eta();
360  // float trk_pT = svtx_track->get_pt();
361 
362  float sin_phi = sin(jet_phi - trk_phi);
363  float cos_phi = cos(jet_phi - trk_phi);
364  float dphi = atan2(sin_phi, cos_phi);
365 
366  if (sqrt((jet_eta - trk_eta) * (jet_eta - trk_eta) + dphi * dphi) < 1.0)
367  {
368  counter_r10++;
369  }
370 
371  if (sqrt((jet_eta - trk_eta) * (jet_eta - trk_eta) + dphi * dphi) > _cut_jet_R) continue;
372 
373  if (svtxtrk_gftrk_map.find(svtx_track->get_id()) != svtxtrk_gftrk_map.end()) {
374  unsigned int trk_index = svtxtrk_gftrk_map[svtx_track->get_id()];
375  unsigned int nclus_mvtx = svtxtrk_nmvtx_map[svtx_track->get_id()];
376  unsigned int nclus_intt = svtxtrk_nintt_map[svtx_track->get_id()];
377 
378  //cout << "NCLUS MVTX: " << nclus_mvtx << std::endl;
379 
380  PHGenFit::Track* rf_phgf_track = rf_phgf_tracks.at(trk_index);
381 
382  if ( _cut_Ncluster ) {
383  if ( nclus_mvtx >= 1 && nclus_intt >= 2 ) {
384  rf_gf_tracks_jet.push_back(rf_phgf_track->getGenFitTrack());
385  // if (trk_pT > 0.5) rf_gf_tracks_jet_pT05.push_back(rf_phgf_track->getGenFitTrack());
386  // if (trk_pT > 1.0) rf_gf_tracks_jet_pT10.push_back(rf_phgf_track->getGenFitTrack());
387  // if (trk_pT > 1.5) rf_gf_tracks_jet_pT15.push_back(rf_phgf_track->getGenFitTrack());
388  // if (trk_pT > 2.0) rf_gf_tracks_jet_pT20.push_back(rf_phgf_track->getGenFitTrack());
389  }
390  } else {
391  rf_gf_tracks_jet.push_back(rf_phgf_track->getGenFitTrack());
392  // if (trk_pT > 0.5) rf_gf_tracks_jet_pT05.push_back(rf_phgf_track->getGenFitTrack());
393  // if (trk_pT > 1.0) rf_gf_tracks_jet_pT10.push_back(rf_phgf_track->getGenFitTrack());
394  // if (trk_pT > 1.5) rf_gf_tracks_jet_pT15.push_back(rf_phgf_track->getGenFitTrack());
395  // if (trk_pT > 2.0) rf_gf_tracks_jet_pT20.push_back(rf_phgf_track->getGenFitTrack());
396  }
397  } else {
398  counter_miss++;
399  }
400  }//trackmap
401 
404  if (rf_gf_tracks_jet.size() > 1) {
405  try {
406  _vertex_finder->findVertices(&rave_vertices_jet, rf_gf_tracks_jet);
407  } catch (...) {
408  std::cout << PHWHERE << "GFRaveVertexFactory::findVertices failed!";
409  }
410  }
411 
412  //cout << "N MISS: " << counter_miss << std::endl;
413  if ( verbosity > 0 )
414  {
415  std::cout << "JET PT: " << jet_pT
416  << ", N TRK: " << int(jet->size_comp())
417  << ", CUT10: " << counter_r10
418  << ", CUT04: " << rf_gf_tracks_jet.size()
419  << ", N VTX: " << rave_vertices_jet.size()
420  << std::endl;
421  }
422 
423  if ( _write_tree )
424  {
425  rv_sv_pT00_nvtx[rv_sv_njets] = rave_vertices_jet.size();
426  for (int ivtx = 0; ivtx < int(rave_vertices_jet.size()); ivtx++) {
427  genfit::GFRaveVertex* rave_vtx = rave_vertices_jet[ivtx];
428  rv_sv_pT00_vtx_x[rv_sv_njets][ivtx] = rave_vtx->getPos().X();
429  rv_sv_pT00_vtx_y[rv_sv_njets][ivtx] = rave_vtx->getPos().Y();
430  rv_sv_pT00_vtx_z[rv_sv_njets][ivtx] = rave_vtx->getPos().Z();
431 
432  rv_sv_pT00_vtx_ex[rv_sv_njets][ivtx] = sqrt(rave_vtx->getCov()[0][0]);
433  rv_sv_pT00_vtx_ey[rv_sv_njets][ivtx] = sqrt(rave_vtx->getCov()[1][1]);
434  rv_sv_pT00_vtx_ez[rv_sv_njets][ivtx] = sqrt(rave_vtx->getCov()[2][2]);
435 
436  rv_sv_pT00_vtx_ntrk[rv_sv_njets][ivtx] = (int)rave_vtx->getNTracks();
437 
438  float vtx_mass, vtx_px, vtx_py, vtx_pz;
439  rv_sv_pT00_vtx_ntrk_good[rv_sv_njets][ivtx] = GetSVMass_mom(rave_vtx, vtx_mass, vtx_px, vtx_py, vtx_pz);
440 
441  //cout << "N TRK: " << rv_sv_pT00_vtx_ntrk[rv_sv_njets][ivtx] << ", GOOD: " << rv_sv_pT00_vtx_ntrk_good[rv_sv_njets][ivtx] << std::endl;
442 
443  TVector3 vec1(vtx_px, vtx_py, vtx_pz);
444  TVector3 vec2(rv_sv_pT00_vtx_x[rv_sv_njets][ivtx] - rv_prim_vtx[0], rv_sv_pT00_vtx_y[rv_sv_njets][ivtx] - rv_prim_vtx[1], rv_sv_pT00_vtx_z[rv_sv_njets][ivtx] - rv_prim_vtx[2]);
445  float theta = vec1.Angle(vec2);
446  float A = vec1.Mag() * sin(theta);
447  float vtx_mass_corr = sqrt(vtx_mass * vtx_mass + A * A) + A;
448 
449  rv_sv_pT00_vtx_mass[rv_sv_njets][ivtx] = vtx_mass;
450  rv_sv_pT00_vtx_mass_corr[rv_sv_njets][ivtx] = vtx_mass_corr;
451  rv_sv_pT00_vtx_pT[rv_sv_njets][ivtx] = sqrt(vtx_px * vtx_px + vtx_py * vtx_py);
452  }
453  }
454 
455  for (genfit::GFRaveVertex* vertex : rave_vertices_jet) {
456  delete vertex;
457  }
458  rave_vertices_jet.clear();
459 
460 
462  rv_sv_jet_id[rv_sv_njets] = jet->get_id();
463  rv_sv_jet_pT[rv_sv_njets] = jet_pT;
464  rv_sv_jet_phi[rv_sv_njets] = jet_phi;
465  if (jet->has_property(static_cast<Jet::PROPERTY>(prop_JetPartonFlavor))) {
466  rv_sv_jet_prop[rv_sv_njets][0] = int(jet->get_property(static_cast<Jet::PROPERTY>(prop_JetPartonFlavor)));
467  rv_sv_jet_prop[rv_sv_njets][1] = int(jet->get_property(static_cast<Jet::PROPERTY>(prop_JetHadronFlavor)));
468  }
469 
470  rv_sv_njets++;
471 
472  }
473 
474 
475 
476 
477  } // _ana_reco
478 
479 
480  //-- Cleanup
481 
482  // fill tree
483  if ( _write_tree )
484  _tree->Fill();
485 
486  //-- End
487  return 0;
488 
489 }
490 
491 int BDiJetModule::End(PHCompositeNode * topNode)
492 {
493  if ( _write_tree )
494  {
495  _f->Write();
496  _f->Close();
497  }
498 
499  return 0;
500 }
501 
502 /*
503  * Initialize the TTree
504  */
505 void BDiJetModule::InitTree()
506 {
507  _f = new TFile(_foutname.c_str(), "RECREATE");
508 
509  _tree = new TTree("ttree", "a withered deciduous tree");
510 
511  _tree->Branch("event", &_b_event, "event/I");
512 
513  if ( _ana_truth )
514  {
515  _tree->Branch("truthjet_n", &_truthjet_n, "truthjet_n/I");
516  _tree->Branch("truthjet_parton_flavor", _truthjet_parton_flavor, "truthjet_parton_flavor[truthjet_n]/I");
517  _tree->Branch("truthjet_hadron_flavor", _truthjet_hadron_flavor, "truthjet_hadron_flavor[truthjet_n]/I");
518  _tree->Branch("truthjet_pt", _truthjet_pt, "truthjet_pt[truthjet_n]/F");
519  _tree->Branch("truthjet_eta", _truthjet_eta,
520  "truthjet_eta[truthjet_n]/F");
521  _tree->Branch("truthjet_phi", _truthjet_phi,
522  "truthjet_phi[truthjet_n]/F");
523  }
524 
525  if ( _ana_reco )
526  {
527  _tree->Branch("gf_prim_vtx", gf_prim_vtx, "gf_prim_vtx[3]/F");
528  _tree->Branch("gf_prim_vtx_err", gf_prim_vtx_err, "gf_prim_vtx_err[3]/F");
529  _tree->Branch("gf_prim_vtx_ntrk", &gf_prim_vtx_ntrk, "gf_prim_vtx_ntrk/I");
530  _tree->Branch("rv_prim_vtx", rv_prim_vtx, "rv_prim_vtx[3]/F");
531  _tree->Branch("rv_prim_vtx_err", rv_prim_vtx_err, "rv_prim_vtx_err[3]/F");
532  _tree->Branch("rv_prim_vtx_ntrk", &rv_prim_vtx_ntrk, "rv_prim_vtx_ntrk/I");
533 
534  _tree->Branch("rv_sv_njets", &rv_sv_njets, "rv_sv_njets/I");
535  _tree->Branch("rv_sv_jet_id", rv_sv_jet_id, "rv_sv_jet_id[rv_sv_njets]/I");
536  _tree->Branch("rv_sv_jet_prop", rv_sv_jet_prop, "rv_sv_jet_prop[rv_sv_njets][2]/I");
537  _tree->Branch("rv_sv_jet_pT", rv_sv_jet_pT, "rv_sv_jet_pT[rv_sv_njets]/F");
538  _tree->Branch("rv_sv_jet_phi", rv_sv_jet_phi, "rv_sv_jet_phi[rv_sv_njets]/F");
539 
540  _tree->Branch("rv_sv_pT00_nvtx", rv_sv_pT00_nvtx, "rv_sv_pT00_nvtx[rv_sv_njets]/I");
541  _tree->Branch("rv_sv_pT00_vtx_x", rv_sv_pT00_vtx_x, "rv_sv_pT00_vtx_x[rv_sv_njets][30]/F");
542  _tree->Branch("rv_sv_pT00_vtx_y", rv_sv_pT00_vtx_y, "rv_sv_pT00_vtx_y[rv_sv_njets][30]/F");
543  _tree->Branch("rv_sv_pT00_vtx_z", rv_sv_pT00_vtx_z, "rv_sv_pT00_vtx_z[rv_sv_njets][30]/F");
544  _tree->Branch("rv_sv_pT00_vtx_ex", rv_sv_pT00_vtx_ex, "rv_sv_pT00_vtx_ex[rv_sv_njets][30]/F");
545  _tree->Branch("rv_sv_pT00_vtx_ey", rv_sv_pT00_vtx_ey, "rv_sv_pT00_vtx_ey[rv_sv_njets][30]/F");
546  _tree->Branch("rv_sv_pT00_vtx_ez", rv_sv_pT00_vtx_ez, "rv_sv_pT00_vtx_ez[rv_sv_njets][30]/F");
547  _tree->Branch("rv_sv_pT00_vtx_ntrk", rv_sv_pT00_vtx_ntrk, "rv_sv_pT00_vtx_ntrk[rv_sv_njets][30]/I");
548  _tree->Branch("rv_sv_pT00_vtx_ntrk_good", rv_sv_pT00_vtx_ntrk_good, "rv_sv_pT00_vtx_ntrk_good[rv_sv_njets][30]/I");
549  _tree->Branch("rv_sv_pT00_vtx_mass", rv_sv_pT00_vtx_mass, "rv_sv_pT00_vtx_mass[rv_sv_njets][30]/F");
550  _tree->Branch("rv_sv_pT00_vtx_mass_corr", rv_sv_pT00_vtx_mass_corr, "rv_sv_pT00_vtx_mass_corr[rv_sv_njets][30]/F");
551  _tree->Branch("rv_sv_pT00_vtx_pT", rv_sv_pT00_vtx_pT, "rv_sv_pT00_vtx_pT[rv_sv_njets][30]/F");
552  }
553 
554 }
555 
556 
557 /*
558  * Reset all the TTree variables
559  * Should be called for each event
560  */
561 void BDiJetModule::ResetVariables()
562 {
563 
564  //-- Truth jet info
565  _truthjet_n = 0;
566  for (int n = 0; n < 10; n++)
567  {
568  _truthjet_parton_flavor[n] = -9999;
569  _truthjet_hadron_flavor[n] = -9999;
570 
571  _truthjet_pt[n] = -1;
572  _truthjet_eta[n] = -1;
573  _truthjet_phi[n] = -1;
574  }
575 
576  //-- Reco info
577  gf_prim_vtx_ntrk = rv_prim_vtx_ntrk = 0;
578  for (int i = 0; i < 3; i++)
579  {
580  gf_prim_vtx[i] = gf_prim_vtx_err[i] = -999;
581  rv_prim_vtx[i] = rv_prim_vtx_err[i] = -999;
582  }//i
583 
584  rv_sv_njets = 0;
585 
586  for (int ijet = 0; ijet < 10; ijet++)
587  {
588  rv_sv_jet_id[ijet] = -999;
589  rv_sv_jet_prop[ijet][0] = -999;
590  rv_sv_jet_prop[ijet][1] = -999;
591  rv_sv_jet_pT[ijet] = -999;
592  rv_sv_jet_phi[ijet] = -999;
593 
594  rv_sv_pT00_nvtx[ijet] = 0;
595 
596  for (int ivtx = 0; ivtx < 30; ivtx++) {
597  rv_sv_pT00_vtx_ntrk[ijet][ivtx] = 0;
598  rv_sv_pT00_vtx_ntrk_good[ijet][ivtx] = 0;
599 
600  rv_sv_pT00_vtx_mass[ijet][ivtx] = -999;
601  rv_sv_pT00_vtx_mass_corr[ijet][ivtx] = -999;
602  rv_sv_pT00_vtx_pT[ijet][ivtx] = -999;
603 
604  rv_sv_pT00_vtx_x[ijet][ivtx] = -999;
605  rv_sv_pT00_vtx_y[ijet][ivtx] = -999;
606  rv_sv_pT00_vtx_z[ijet][ivtx] = -999;
607  rv_sv_pT00_vtx_ex[ijet][ivtx] = -999;
608  rv_sv_pT00_vtx_ey[ijet][ivtx] = -999;
609  rv_sv_pT00_vtx_ez[ijet][ivtx] = -999;
610 
611  }//ivtx
612  }//ijet
613 
614 }
615 
616 /*
617  * Get all the necessary nodes from the node tree
618  */
619 int BDiJetModule::GetNodes(PHCompositeNode * topNode)
620 {
621 
622  // Truth jet node
623  _jetmap_truth = findNode::getClass<JetMap>(topNode,
624  "AntiKt_Truth_r04");
625  if ( !_jetmap_truth && _ana_truth && _ievent < 2)
626  {
627  std::cout << PHWHERE << " AntiKt_Truth_r04 node not found ... aborting" << std::endl;
628  return Fun4AllReturnCodes::ABORTEVENT;
629  }
630 
631  // Input Svtx Clusters
632  _clustermap = findNode::getClass<SvtxClusterMap>(topNode, "SvtxClusterMap");
633  if (!_clustermap && _ana_reco && _ievent < 2) {
634  std::cout << PHWHERE << " SvtxClusterMap node not found on node tree"
635  << std::endl;
636  return Fun4AllReturnCodes::ABORTEVENT;
637  }
638 
639  // Input Svtx Tracks
640  _trackmap = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
641  if (!_trackmap && _ana_reco && _ievent < 2) {
642  std::cout << PHWHERE << " SvtxClusterMap node not found on node tree"
643  << std::endl;
644  return Fun4AllReturnCodes::ABORTEVENT;
645  }
646 
647  // Input Svtx Vertices
648  _vertexmap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
649  if (!_vertexmap && _ana_reco && _ievent < 2) {
650  std::cout << PHWHERE << " SvtxVertexrMap node not found on node tree"
651  << std::endl;
652  return Fun4AllReturnCodes::ABORTEVENT;
653  }
654 
655 
656 
657  return Fun4AllReturnCodes::EVENT_OK;
658 
659 }
660 
661 
662 /*
663  * Load reconstructed track into GenFit track for input to Rave
664  */
665 PHGenFit::Track* BDiJetModule::MakeGenFitTrack(PHCompositeNode * topNode,
666  const SvtxTrack * intrack,
667  const SvtxVertex * vertex)
668 {
669 
670 
671  if (!intrack) {
672  std::cerr << PHWHERE << " Input SvtxTrack is NULL!" << std::endl;
673  return NULL;
674  }
675 
676  SvtxHitMap* hitsmap = NULL;
677  hitsmap = findNode::getClass<SvtxHitMap>(topNode, "SvtxHitMap");
678  if (!hitsmap) {
679  std::cout << PHWHERE << "ERROR: Can't find node SvtxHitMap" << std::endl;
680  return NULL;
681  }
682 
683  PHG4CellContainer* cells_svtx = findNode::getClass<PHG4CellContainer>(topNode, "G4CELL_SVTX");
684  PHG4CellContainer* cells_intt = findNode::getClass<PHG4CellContainer>(topNode, "G4CELL_SILICON_TRACKER");
685  PHG4CellContainer* cells_maps = findNode::getClass<PHG4CellContainer>(topNode, "G4CELL_MAPS");
686 
687  if (_use_ladder_geom and !cells_svtx and !cells_intt and !cells_maps) {
688  std::cout << PHWHERE << "No PHG4CellContainer found!" << std::endl;
689  return NULL;
690  }
691 
692  PHG4CylinderGeomContainer* geom_container_intt = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_SILICON_TRACKER");
693  PHG4CylinderGeomContainer* geom_container_maps = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_MAPS");
694 
695  if (_use_ladder_geom and !geom_container_intt and !geom_container_maps) {
696  std::cout << PHWHERE << "No PHG4CylinderGeomContainer found!" << std::endl;
697  return NULL;
698  }
699 
700  // prepare seed
701  TVector3 seed_mom(100, 0, 0);
702  TVector3 seed_pos(0, 0, 0);
703  TMatrixDSym seed_cov(6);
704  for (int i = 0; i < 6; i++) {
705  for (int j = 0; j < 6; j++) {
706  seed_cov[i][j] = 100.;
707  }
708  }
709 
710  /*
711  TVector3 seed_pos(intrack->get_x(), intrack->get_y(), intrack->get_z());
712  TVector3 seed_mom(intrack->get_px(), intrack->get_py(), intrack->get_pz());
713  TMatrixDSym seed_cov(6);
714  for (int i=0; i<6; i++){
715  for (int j=0; j<6; j++){
716  seed_cov[i][j] = intrack->get_error(i,j);
717  }
718  }
719  */
720 
721  // Create measurements
722  std::vector<PHGenFit::Measurement*> measurements;
723 
724 
725  std::map<float, unsigned int> m_r_cluster_id;
726 
727  for (auto iter = intrack->begin_clusters(); iter != intrack->end_clusters(); ++iter) {
728  unsigned int cluster_id = *iter;
729  SvtxCluster* cluster = _clustermap->get(cluster_id);
730  float x = cluster->get_x();
731  float y = cluster->get_y();
732  float r = sqrt(x * x + y * y);
733  m_r_cluster_id.insert(std::pair<float, unsigned int>(r, cluster_id));
734  }
735 
736  for (auto iter = m_r_cluster_id.begin(); iter != m_r_cluster_id.end(); ++iter) {
737  //for (SvtxTrack::ConstClusterIter iter = intrack->begin_clusters(); iter != intrack->end_clusters(); ++iter){
738  unsigned int cluster_id = iter->second;
739  //unsigned int cluster_id = *iter;
740  SvtxCluster* cluster = _clustermap->get(cluster_id);
741  if (!cluster) {
742  LogError("No cluster Found!");
743  continue;
744  }
745 
746  TVector3 pos(cluster->get_x(), cluster->get_y(), cluster->get_z());
747 
748  seed_mom.SetPhi(pos.Phi());
749  seed_mom.SetTheta(pos.Theta());
750 
751  //TODO use u, v explicitly?
752  TVector3 n(cluster->get_x(), cluster->get_y(), 0);
753 
754  unsigned int begin_hit_id = *(cluster->begin_hits());
755  SvtxHit* svtxhit = hitsmap->find(begin_hit_id)->second;
756 
757  PHG4Cell* cell_svtx = nullptr;
758  PHG4Cell* cell_intt = nullptr;
759  PHG4Cell* cell_maps = nullptr;
760 
761  if (_use_ladder_geom and cells_svtx) cell_svtx = cells_svtx->findCell(svtxhit->get_cellid());
762  if (_use_ladder_geom and cells_intt) cell_intt = cells_intt->findCell(svtxhit->get_cellid());
763  if (_use_ladder_geom and cells_maps) cell_maps = cells_maps->findCell(svtxhit->get_cellid());
764  if (_use_ladder_geom and !(cell_svtx or cell_intt or cell_maps)) {
765  if (verbosity >= 0)
766  LogError("!(cell_svtx or cell_intt or cell_maps)");
767  continue;
768  }
769 
770  //float phi_tilt[7] = {0.304, 0.304, 0.304, 0.244, 0.244, 0.209, 0.201};
771  unsigned int layer = cluster->get_layer();
772 
773  //NEW
774  if (_use_ladder_geom) {
775  if (cell_maps) {
776  PHG4Cell* cell = cell_maps;
777 
778  int stave_index = cell->get_stave_index();
779  int half_stave_index = cell->get_half_stave_index();
780  int module_index = cell->get_module_index();
781  int chip_index = cell->get_chip_index();
782 
783  double ladder_location[3] = { 0.0, 0.0, 0.0 };
784  PHG4CylinderGeom_MAPS *geom = (PHG4CylinderGeom_MAPS*) geom_container_maps->GetLayerGeom(layer);
785  // returns the center of the sensor in world coordinates - used to get the ladder phi location
786  geom->find_sensor_center(stave_index, half_stave_index, module_index, chip_index, ladder_location);
787  n.SetXYZ(ladder_location[0], ladder_location[1], 0);
788  n.RotateZ(geom->get_stave_phi_tilt());
789  } else if (cell_intt) {
790  PHG4Cell* cell = cell_intt;
791  PHG4CylinderGeom_Siladders* geom = (PHG4CylinderGeom_Siladders*) geom_container_intt->GetLayerGeom(layer);
792  double hit_location[3] = { 0.0, 0.0, 0.0 };
793  geom->find_segment_center(cell->get_ladder_z_index(), cell->get_ladder_phi_index(), hit_location);
794 
795  n.SetXYZ(hit_location[0], hit_location[1], 0);
796  n.RotateZ(geom->get_strip_phi_tilt());
797  }
798  }
799 
800  /*
801  //OLD
802  if ((cells_maps and geom_container_maps) and layer < 3) {
803  unsigned int begin_hit_id = *(cluster->begin_hits());
804  SvtxHit* hit = hitsmap->find(begin_hit_id)->second;
805  PHG4Cell* cell = cells_maps->findCell(hit->get_cellid());
806  int stave_index = cell->get_stave_index();
807  int half_stave_index = cell->get_half_stave_index();
808  int module_index = cell->get_module_index();
809  int chip_index = cell->get_chip_index();
810 
811  double ladder_location[3] = {0.0, 0.0, 0.0};
812  PHG4CylinderGeom_MAPS *geom = (PHG4CylinderGeom_MAPS*) geom_container_maps->GetLayerGeom(layer);
813  geom->find_sensor_center(stave_index, half_stave_index, module_index, chip_index, ladder_location);
814  n.SetXYZ(ladder_location[0], ladder_location[1], 0);
815  n.RotateZ(phi_tilt[layer]);
816  } else if ((cells_intt and geom_container_intt) and pos.Perp() < 30.) {
817  unsigned int begin_hit_id = *(cluster->begin_hits());
818  SvtxHit* hit = hitsmap->find(begin_hit_id)->second;
819  PHG4Cell* cell = cells_intt->findCell(hit->get_cellid());
820  PHG4CylinderGeom_Siladders* geom = (PHG4CylinderGeom_Siladders*) geom_container_intt->GetLayerGeom(layer);
821  double hit_location[3] = { 0.0, 0.0, 0.0 };
822  geom->find_segment_center(cell->get_ladder_z_index(),cell->get_ladder_phi_index(), hit_location);
823 
824  n.SetXYZ(hit_location[0], hit_location[1], 0);
825  n.RotateZ(phi_tilt[layer]);
826  }
827  */
828 
829  PHGenFit::Measurement* meas = new PHGenFit::PlanarMeasurement(pos, n,
830  cluster->get_rphi_error(), cluster->get_z_error());
831 
832  measurements.push_back(meas);
833  }
834 
835 
836  //TODO Add multiple TrackRep choices.
837  genfit::AbsTrackRep* rep = new genfit::RKTrackRep(_primary_pid_guess);
838  PHGenFit::Track* track(new PHGenFit::Track(rep, seed_pos, seed_mom, seed_cov));
839  track->addMeasurements(measurements);
840 
841  if (_fitter->processTrack(track, false) != 0) {
842  if (verbosity >= 1)
843  LogWarning("Track fitting failed");
844  return NULL;
845  }
846 
847  TVector3 vertex_position(0, 0, 0);
848  if (vertex) {
849  vertex_position.SetXYZ(vertex->get_x(), vertex->get_y(), vertex->get_z());
850  }
851 
852  std::shared_ptr<genfit::MeasuredStateOnPlane> gf_state_vertex_ca = NULL;
853  try {
854  gf_state_vertex_ca = std::shared_ptr < genfit::MeasuredStateOnPlane> (track->extrapolateToPoint(vertex_position));
855  } catch (...) {
856  if (verbosity >= 2)
857  LogWarning("extrapolateToPoint failed!");
858  return NULL;
859  }
860 
861  TVector3 mom = gf_state_vertex_ca->getMom();
862  TMatrixDSym cov = gf_state_vertex_ca->get6DCov();
863 
864  //cout << "OUT Ex: " << sqrt(cov[0][0]) << ", Ey: " << sqrt(cov[1][1]) << ", Ez: " << sqrt(cov[2][2]) << std::endl;
865 
866  //cout << "IN Px: " << intrack->get_px() << ", Py: " << intrack->get_py() << ", Pz: " << intrack->get_pz() << std::endl;
867  //cout << "OUT Px: " << mom.X() << ", Py: " << mom.Y() << ", Pz: " << mom.Z() << std::endl;
868 
869  return track;
870 }
871 
872 /*
873  * Fill SvtxVertexMap from GFRaveVertexes and Tracks
874  */
875 void BDiJetModule::FillVertexMap(
876  const std::vector<genfit::GFRaveVertex*>& rave_vertices,
877  const std::vector<genfit::Track*>& gf_tracks) {
878 
879  for (unsigned int ivtx = 0; ivtx < rave_vertices.size(); ++ivtx)
880  {
881  genfit::GFRaveVertex* rave_vtx = rave_vertices[ivtx];
882 
883  //cout << "V0 x: " << rave_vtx->getPos().X() << ", y: " << rave_vtx->getPos().Y() << ", z: " << rave_vtx->getPos().Z() << std::endl;
884  rv_prim_vtx[0] = rave_vtx->getPos().X();
885  rv_prim_vtx[1] = rave_vtx->getPos().Y();
886  rv_prim_vtx[2] = rave_vtx->getPos().Z();
887 
888  rv_prim_vtx_err[0] = sqrt(rave_vtx->getCov()[0][0]);
889  rv_prim_vtx_err[1] = sqrt(rave_vtx->getCov()[1][1]);
890  rv_prim_vtx_err[2] = sqrt(rave_vtx->getCov()[2][2]);
891 
892  rv_prim_vtx_ntrk = rave_vtx->getNTracks();
893  } // ivtx
894 
895  for (SvtxVertexMap::Iter iter = _vertexmap->begin();
896  iter != _vertexmap->end();
897  ++iter)
898  {
899  SvtxVertex *svtx_vertex = iter->second;
900 
901  //cout << "V1 x: " << svtx_vertex->get_x() << ", y: " << svtx_vertex->get_y() << ", z: " << svtx_vertex->get_z() << std::endl;
902  gf_prim_vtx[0] = svtx_vertex->get_x();
903  gf_prim_vtx[1] = svtx_vertex->get_y();
904  gf_prim_vtx[2] = svtx_vertex->get_z();
905 
906  gf_prim_vtx_err[0] = sqrt(svtx_vertex->get_error(0, 0));
907  gf_prim_vtx_err[1] = sqrt(svtx_vertex->get_error(1, 1));
908  gf_prim_vtx_err[2] = sqrt(svtx_vertex->get_error(2, 2));
909 
910  gf_prim_vtx_ntrk = int(svtx_vertex->size_tracks());
911  } // iter on SvtxVertexMap
912 
913  return;
914 }
915 
916 int BDiJetModule::GetSVMass_mom(
917  const genfit::GFRaveVertex* rave_vtx,
918  float & vtx_mass,
919  float & vtx_px,
920  float & vtx_py,
921  float & vtx_pz
922 ) {
923 
924  float sum_E = 0, sum_px = 0, sum_py = 0, sum_pz = 0;
925 
926  int N_good = 0;
927 
928  for (unsigned int itrk = 0; itrk < rave_vtx->getNTracks(); itrk++) {
929  TVector3 mom_trk = rave_vtx->getParameters(itrk)->getMom();
930 
931  double w_trk = rave_vtx->getParameters(itrk)->getWeight();
932 
933  sum_px += mom_trk.X();
934  sum_py += mom_trk.Y();
935  sum_pz += mom_trk.Z();
936  sum_E += sqrt(mom_trk.Mag2() + 0.140 * 0.140);
937 
938  //cout << "W: " << w_trk << std::endl;
939  if ( w_trk > 0.6 ) {
940  N_good++;
941  }
942 
943  }//itrk
944 
945  vtx_mass = sqrt(sum_E * sum_E - sum_px * sum_px - sum_py * sum_py - sum_pz * sum_pz);
946  vtx_px = sum_px;
947  vtx_py = sum_py;
948  vtx_pz = sum_pz;
949 
950  //cout << "Mass: " << vtx_mass << ", pT: " << vtx_pT << std::endl;
951  return N_good;
952 }
953 
954 
955