Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
FastTrackingEval.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file FastTrackingEval.cc
1 
8 #include "FastTrackingEval.h"
9 
10 #include <phool/phool.h>
11 #include <phool/getClass.h>
12 #include <fun4all/Fun4AllReturnCodes.h>
13 #include <g4main/PHG4HitContainer.h>
14 #include <g4main/PHG4TruthInfoContainer.h>
15 #include <g4main/PHG4Particle.h>
16 #include <g4main/PHG4Hit.h>
17 #include <g4main/PHG4VtxPoint.h>
18 #include <fun4all/PHTFileServer.h>
19 #include <fun4all/Fun4AllServer.h>
20 
21 /*
22 #include <g4hough/SvtxVertexMap.h>
23 #include <g4hough/SvtxVertex.h>
24 #include <g4hough/SvtxTrackMap.h>
25 #include <g4hough/SvtxTrack.h>
26 #include <g4hough/SvtxTrack_FastSim.h>
27 #include <g4hough/SvtxClusterMap.h>
28 #include <g4hough/SvtxCluster.h>
29 #include <g4hough/SvtxHitMap.h>
30 #include <g4hough/SvtxHit.h>
31 */
32 
33 #include <trackbase_historic/SvtxVertexMap.h>
34 #include <trackbase_historic/SvtxVertex.h>
35 #include <trackbase_historic/SvtxTrackMap.h>
36 #include <trackbase_historic/SvtxTrack.h>
37 #include <trackbase_historic/SvtxTrack_FastSim.h>
38 //#include <trackbase_historic/SvtxClusterMap.h>
39 //#include <trackbase_historic/SvtxCluster.h>
40 //#include <trackbase_historic/SvtxHitMap.h>
41 //#include <trackbase_historic/SvtxHit.h>
42 
43 #include <g4eval/SvtxEvalStack.h>
44 #include <g4eval/SvtxTrackEval.h>
45 #include <g4eval/SvtxClusterEval.h>
46 #include <g4eval/SvtxTruthEval.h>
47 #include <g4eval/SvtxVertexEval.h>
48 #include <g4eval/SvtxHitEval.h>
49 
50 #include <TTree.h>
51 #include <TH2D.h>
52 #include <TVector3.h>
53 
54 #include <math.h>
55 
56 #include <iostream>
57 
58 #define LogDebug(exp) std::cout<<"DEBUG: " <<__FILE__<<": "<<__LINE__<<": "<< exp <<"\n"
59 #define LogError(exp) std::cout<<"ERROR: " <<__FILE__<<": "<<__LINE__<<": "<< exp <<"\n"
60 #define LogWarning(exp) std::cout<<"WARNING: " <<__FILE__<<": "<<__LINE__<<": "<< exp <<"\n"
61 
62 using namespace std;
63 
64 //----------------------------------------------------------------------------//
65 //-- Constructor:
66 //-- simple initialization
67 //----------------------------------------------------------------------------//
68 FastTrackingEval::FastTrackingEval(const string &name, const string &filename, const string &trackmapname) :
69  SubsysReco(name),
70  _outfile_name(filename),
71  _trackmapname(trackmapname),
72  _event(0),
73  _flags(NONE),
74  _eval_tree_tracks( NULL)
75 {
76 }
77 
78 //----------------------------------------------------------------------------//
79 //-- Init():
80 //-- Intialize all histograms, trees, and ntuples
81 //----------------------------------------------------------------------------//
82 int FastTrackingEval::Init(PHCompositeNode *topNode) {
83  cout << PHWHERE << " Openning file " << _outfile_name << endl;
84  PHTFileServer::get().open(_outfile_name, "RECREATE");
85 
86  // create TTree
87  _eval_tree_tracks = new TTree("tracks", "FastSim Eval => tracks");
88  _eval_tree_tracks->Branch("event", &event, "event/I");
89  _eval_tree_tracks->Branch("gtrackID", &gtrackID, "gtrackID/I");
90  _eval_tree_tracks->Branch("gflavor", &gflavor, "gflavor/I");
91  _eval_tree_tracks->Branch("gpx", &gpx, "gpx/F");
92  _eval_tree_tracks->Branch("gpy", &gpy, "gpy/F");
93  _eval_tree_tracks->Branch("gpz", &gpz, "gpz/F");
94  _eval_tree_tracks->Branch("gpt", &gpt, "gpt/F");
95  _eval_tree_tracks->Branch("gp", &gp, "gp/F");
96  _eval_tree_tracks->Branch("gtheta", &gtheta, "gtheta/F");
97  _eval_tree_tracks->Branch("geta", &geta, "geta/F");
98  _eval_tree_tracks->Branch("gphi", &gphi, "gphi/F");
99  _eval_tree_tracks->Branch("gvx", &gvx, "gvx/F");
100  _eval_tree_tracks->Branch("gvy", &gvy, "gvy/F");
101  _eval_tree_tracks->Branch("gvz", &gvz, "gvz/F");
102  _eval_tree_tracks->Branch("trackID", &trackID, "trackID/I");
103  _eval_tree_tracks->Branch("charge", &charge, "charge/I");
104  _eval_tree_tracks->Branch("nhits", &nhits, "nhits/I");
105  _eval_tree_tracks->Branch("px", &px, "px/F");
106  _eval_tree_tracks->Branch("py", &py, "py/F");
107  _eval_tree_tracks->Branch("pz", &pz, "pz/F");
108  _eval_tree_tracks->Branch("pt", &pt, "pt/F");
109  _eval_tree_tracks->Branch("p", &p, "p/F");
110  _eval_tree_tracks->Branch("theta", &theta, "theta/F");
111  _eval_tree_tracks->Branch("eta", &eta, "eta/F");
112  _eval_tree_tracks->Branch("phi", &phi, "phi/F");
113  _eval_tree_tracks->Branch("dca2d", &dca2d, "dca2d/F");
114 
115  _h2d_Delta_mom_vs_truth_eta = new TH2D("_h2d_Delta_mom_vs_truth_eta",
116  "#frac{#Delta p}{truth p} vs. truth #eta", 54, -4.5, +4.5, 1000, -1,
117  1);
118 
119  _h2d_Delta_mom_vs_truth_mom = new TH2D("_h2d_Delta_mom_vs_truth_mom",
120  "#frac{#Delta p}{truth p} vs. truth p", 41, -0.5, 40.5, 1000, -1,
121  1);
122 
123  return Fun4AllReturnCodes::EVENT_OK;
124 }
125 
126 //----------------------------------------------------------------------------//
127 //-- process_event():
128 //-- Call user instructions for every event.
129 //-- This function contains the analysis structure.
130 //----------------------------------------------------------------------------//
131 int FastTrackingEval::process_event(PHCompositeNode *topNode) {
132  _event++;
133  if (Verbosity() >= 2 and _event % 1000 == 0)
134  cout << PHWHERE << "Events processed: " << _event << endl;
135 
136  //std::cout << "Opening nodes" << std::endl;
137  GetNodes(topNode);
138 
139  //std::cout << "Filling trees" << std::endl;
140  fill_tree(topNode);
141  //std::cout << "DONE" << std::endl;
142 
143  return Fun4AllReturnCodes::EVENT_OK;
144 }
145 
146 //----------------------------------------------------------------------------//
147 //-- End():
148 //-- End method, wrap everything up
149 //----------------------------------------------------------------------------//
150 int FastTrackingEval::End(PHCompositeNode *topNode) {
151 
152  PHTFileServer::get().cd(_outfile_name);
153 
154  _eval_tree_tracks->Write();
155 
156  _h2d_Delta_mom_vs_truth_eta->Write();
157  _h2d_Delta_mom_vs_truth_mom->Write();
158 
159  //PHTFileServer::get().close();
160 
161  return Fun4AllReturnCodes::EVENT_OK;
162 }
163 
164 //----------------------------------------------------------------------------//
165 //-- fill_tree():
166 //-- Fill the trees with truth, track fit, and cluster information
167 //----------------------------------------------------------------------------//
168 void FastTrackingEval::fill_tree(PHCompositeNode *topNode) {
169 
170  // Make sure to reset all the TTree variables before trying to set them.
171  reset_variables();
172  //std::cout << "A1" << std::endl;
173  event = _event;
174 
175  if (!_truth_container) {
176  LogError("_truth_container not found!");
177  return;
178  }
179 
180  if (!_trackmap) {
181  LogError("_trackmap not found!");
182  return;
183  }
184 
185  PHG4TruthInfoContainer::ConstRange range =
186  _truth_container->GetPrimaryParticleRange();
187  //std::cout << "A2" << std::endl;
188  for (PHG4TruthInfoContainer::ConstIterator truth_itr = range.first;
189  truth_itr != range.second; ++truth_itr) {
190 
191  PHG4Particle* g4particle = truth_itr->second;
192  if(!g4particle) {
193  LogDebug("");
194  continue;
195  }
196  //std::cout << "B1" << std::endl;
197 
198  SvtxTrack_FastSim* track = NULL;
199 
200  //std::cout << "TRACKmap size " << _trackmap->size() << std::endl;
201  for (SvtxTrackMap::ConstIter track_itr = _trackmap->begin();
202  track_itr != _trackmap->end();
203  track_itr++) {
204  //std::cout << "TRACK * " << track_itr->first << std::endl;
205  SvtxTrack_FastSim* temp = dynamic_cast<SvtxTrack_FastSim*>(track_itr->second);
206  if(!temp) {
207  std::cout << "ERROR CASTING PARTICLE!" << std::endl;
208  continue;
209  }
210  //std::cout << " PARTICLE!" << std::endl;
211 
212  if ((temp->get_truth_track_id() - g4particle->get_track_id()) == 0) {
213  track = temp;
214  }
215  }
216 
217  //std::cout << "B2" << std::endl;
218  gtrackID = g4particle->get_track_id();
219  gflavor = g4particle->get_pid();
220 
221  gpx = g4particle->get_px();
222  gpy = g4particle->get_py();
223  gpz = g4particle->get_pz();
224  gpt = sqrt(gpx*gpx+gpy*gpy);
225  gp = sqrt(gpx*gpx+gpy*gpy+gpz*gpz);
226  gtheta = atan2(gpt,gpz);
227  geta=-1.*log(tan(gtheta/2.));
228  gphi = atan(gpy/gpx);
229 
230  if (track) {
231  //std::cout << "C1" << std::endl;
232  trackID = track->get_id();
233  charge = track->get_charge();
234  nhits = track->size_clusters();
235 
236  px = track->get_px();
237  py = track->get_py();
238  pz = track->get_pz();
239  pt = sqrt(px*px+py*py);
240  p = sqrt(px*px+py*py+pz*pz);
241  theta = atan2(pt,pz);
242  eta=-1.*log(tan(theta/2.));
243  phi = atan(py/px);
244  dca2d = track->get_dca2d();
245 
246  TVector3 truth_mom(gpx,gpy,gpz);
247  TVector3 reco_mom(px, py, pz);
248  //std::cout << "C2" << std::endl;
249 
250  _h2d_Delta_mom_vs_truth_mom->Fill(truth_mom.Mag(), (reco_mom.Mag()-truth_mom.Mag())/truth_mom.Mag());
251  _h2d_Delta_mom_vs_truth_eta->Fill(truth_mom.Eta(), (reco_mom.Mag()-truth_mom.Mag())/truth_mom.Mag());
252  }
253  //std::cout << "B3" << std::endl;
254 
255  _eval_tree_tracks->Fill();
256 
257  }
258  //std::cout << "A3" << std::endl;
259 
260  return;
261 
262 }
263 
264 //----------------------------------------------------------------------------//
265 //-- reset_variables():
266 //-- Reset all the tree variables to their default values.
267 //-- Needs to be called at the start of every event
268 //----------------------------------------------------------------------------//
269 void FastTrackingEval::reset_variables() {
270  event = -9999;
271 
272  //-- truth
273  gtrackID = -9999;
274  gflavor = -9999;
275  gpx = -9999;
276  gpy = -9999;
277  gpz = -9999;
278  gpt = -9999;
279  gp = -9999;
280  gtheta = -9999;
281  geta = -9999;
282  gphi = -9999;
283  gvx = -9999;
284  gvy = -9999;
285  gvz = -9999;
286 
287  //-- reco
288  trackID = -9999;
289  charge = -9999;
290  nhits = -9999;
291  px = -9999;
292  py = -9999;
293  pz = -9999;
294  pt = -9999;
295  p = -9999;
296  theta = -9999;
297  eta = -9999;
298  phi = -9999;
299  dca2d = -9999;
300 }
301 
302 //----------------------------------------------------------------------------//
303 //-- GetNodes():
304 //-- Get all the all the required nodes off the node tree
305 //----------------------------------------------------------------------------//
306 int FastTrackingEval::GetNodes(PHCompositeNode * topNode) {
307  //DST objects
308  //Truth container
309  _truth_container = findNode::getClass<PHG4TruthInfoContainer>(topNode,
310  "G4TruthInfo");
311  if (!_truth_container && _event < 2) {
312  cout << PHWHERE << " PHG4TruthInfoContainer node not found on node tree"
313  << endl;
314  return Fun4AllReturnCodes::ABORTEVENT;
315  }
316 
317  _trackmap = findNode::getClass<SvtxTrackMap>(topNode,
318  _trackmapname);
319  //std::cout << _trackmapname.c_str() << std::endl;
320  if (!_trackmap && _event < 2) {
321  cout << PHWHERE << "SvtxTrackMap node with name "
322  << _trackmapname
323  <<" not found on node tree"
324  << endl;
325  return Fun4AllReturnCodes::ABORTEVENT;
326  }
327 
328  return Fun4AllReturnCodes::EVENT_OK;
329 }
330