FairRoot/PandaRoot
runOnlineDisplayMCCheckFaster2.C
Go to the documentation of this file.
1 #include "drawhelpers.hxx"
2 
3 void SaveAndUpdateHisto(TH1* currenthisto, TFile& storagefile)
4 {
5  if (storagefile.Get(currenthisto->GetName()) != 0) {
6  currenthisto->Add((TH1*)storagefile.Get(currenthisto->GetName()));
7  }
8  cout << currenthisto->GetName() << ": " << currenthisto->GetEntries() << endl;
9  currenthisto->Write(currenthisto->GetName(), TObject::kWriteDelete);
10 }
11 
12 int runOnlineDisplayMCCheckFaster2(Int_t maximumTime = 2050) {
13  gROOT->Reset();
14  // Verbosity level (0=quiet, 1=event level, 2=track level, 3=debug)
15  Int_t iVerbose = 0;
16 
17  enum DetectorId {
18  kDCH,kDRC,kDSK,kEMC,kGEM,kLUMI,kMDT,kMVD,kRPC,kSTT,kTPC,kTOF,kFTS,kHYPG,kHYP};
19 
20  TStopwatch timer;
21  timer.Start();
22  gROOT->LoadMacro("$VMCWORKDIR/gconfig/rootlogon.C");
23  rootlogon();
24 
25  TDatabasePDG *pdg = new TDatabasePDG();
26 
27  // Reco Parameters
28  Double_t stthitlifetime = 245;
29 
30  TString simFileName = "Sim_Dpm_500.root";
31  TString parFileName = "Sim_Dpm_500_params.root";
32  TString digiFileName = "Sim_Dpm_500_digi.root";
33  TString recoFileName = "Sim_Dpm_500_reco.root";
34  TString outFileName = "Sim_Dpm_500_streamdisplay.root";
35  TString histofilename = "onlinetrkmccheck.root";
36 
37  TFile filedigi(digiFileName.Data());
38  TFile filereco(recoFileName.Data());
39  TFile filerecopixel(recoFileName.Data());
40 
41  TTree *treedigi = (TTree*) filedigi.Get("pndsim");
42  TTree *recotree = (TTree*) filereco.Get("pndsim");
43 
44  TFile simfile(simFileName.Data());
45  TTree* mctree = (TTree*) simfile.Get("pndsim");
46  TClonesArray* mcarray = new TClonesArray("PndMCTrack");
47  mctree->SetBranchAddress("MCTrack", &mcarray);
48 
49  FairRunAna *fRun = new FairRunAna();
50  fRun->SetInputFile(simFileName.Data());
51  fRun->AddFriend(recoFileName.Data());
52  fRun->AddFriend(digiFileName.Data());
53  fRun->SetOutputFile(outFileName.Data());
54  fRun->Init();
55  FairRuntimeDb* rtdb = fRun->GetRuntimeDb();
56 
57  PndOnlineGeometryManager *online_geometry = new PndOnlineGeometryManager(rtdb, parFileName.Data());
58 
59  PndOnlineManager *online = new PndOnlineManager();
60  online->AddGeometryManager(online_geometry);
61 
62  online->LoadStream(PndOnline::kSTT,"STTSortedHits",stthitlifetime);
63  online->LoadStream(PndOnline::kMVDPixel,"MVDHitsPixel",stthitlifetime);
64  online->LoadStream(PndOnline::kMVDStrip,"MVDHitsStrip",stthitlifetime);
65  //online->SetHESRRevolutionDuration(3000);
66 
67  TClonesArray* tubearray = online_geometry->GetDetectorGeometry(PndOnline::kSTT);
68 
69  // load task
70  PndOnlineSttTripletFinder *triplet_finder = new PndOnlineSttTripletFinder(online, tubearray);
71  online->AddTask((FairTask*)triplet_finder);
72 
73  // initialize display
74  TCanvas* c1 = new TCanvas("c1");
75  c1->Range(-42,-42,42,42);
76  c1->SetCanvasSize(1200, 1200);
77  TText* mytext = new TText();
78 
79  // event loop!
80  int current_time = 0;
81  int drawtime = 0;
82  int delta_t = 2050;
83  int drawdeltat = 5;
84  int final_time = maximumTime;
85 
86  online->SetHESRRevolutionDuration(delta_t);
87 
88  online->Init();
89 
90  TObjArray* online_fairhits_stt = 0;
91  TObjArray* online_tracks = 0;
92  TObjArray* online_tracks2 = 0;
93  TObjArray* online_fairhits_mvd = 0;
94 
95  std::vector<double> mctracktime;
96  std::vector<double> mctrackpt;
97  std::vector<double> mctrackphi;
98  std::vector<int> mctrackstatus;
99  std::vector<double> recotracktime;
100  std::vector<double> recotrackpt;
101  std::vector<double> recotrackphi;
102  std::vector<int> recotrackstatus;
103  std::vector<int> recotrackmcid;
104  int goodrecotrack = 0;
105  int badrecotrack = 0;
106  int totalmctracks = 0;
107  int goodmctracks = 0;
108  int totalrecotracks = 0;
109  int goodrecotracks = 0;
110  int norecotracks = 0;
111  int deltaphifailed = 0;
112  int deltaptfailed = 0;
113  int deltabothfailed = 0;
114  int secondaries = 0;
115  int globalmaxmcid = 0;
116  int mctracknonprime = 0;
117  int mctracknonrecomomentum = 0;
118  int lowpthit = 0;
119  int mctracknonhit = 0;
120  int mctrackgrandtotal = 0;
121  int mctrackrecoattempt = 0;
122  int foundrecotrackall = 0;
123  int foundrecotrack90 = 0;
124  int foundrecotrack80 = 0;
125  int foundrecotrackbad = 0;
126  int mctracknotrecoed = 0;
127  int lowbetagamma = 0;
128  int pdgcodetoolarge = 0;
129 
130  PndMCBookkeeper foundmapall;
131  foundmapall.Clear();
132  PndMCBookkeeper foundmap90;
133  foundmap90.Clear();
134  PndMCBookkeeper foundmap80;
135  foundmap80.Clear();
136  PndMCBookkeeper foundmapbad;
137  foundmapbad.Clear();
138 
139  int mctracknullpointers = 0;
140  int recotracknullpointers = 0;
141  int mcidgood90 = 0;
142  int mcidgood80 = 0;
143  int mcidbad = 0;
144 
145  //output histos
146  TH2D* hmcall = new TH2D("hmcall", "hmcall", 1500, -14, 16, 75, 0, 1.5);
147  TH2D* hmcprime = new TH2D("hmcprime", "hmcprime", 1500, -14, 16, 75, 0, 1.5);
148  TH2D* hmcrecoable = new TH2D("hmcrecoable", "hmcrecoable", 1500, -14, 16, 75, 0, 1.5);
149  TH2D* hdafuq = new TH2D("hdafuq", "hdafuq", 1500, -14, 16, 75, 0, 1.5);
150  TH2D* hmcrecoed = new TH2D("hmcrecoed", "hmcrecoed", 1500, -14, 16, 75, 0, 1.5);
151  TH2D* hmcnotrecoed = new TH2D("hmcnotrecoed", "hmcnotrecoed", 1500, -14, 16, 75, 0, 1.5);
152 
153  TH1D* hphires = new TH1D("hphires", "hphires", 701, -3.505, 3.505);
154  TH2D* hptres = new TH2D("hptres", "hptres", 300, -1400, 1600, 200, -1, 1);
155  TH1D* hphiresmcall = new TH1D("hphiresmcall", "hphiresmcall", 701, -3.505, 3.505);
156  TH2D* hptresmcall = new TH2D("hptresmcall", "hptresmcall", 300, -1400, 1600, 200, -1, 1);
157  TH1D* hphiresmc90 = new TH1D("hphiresmc90", "hphiresmc90", 701, -3.505, 3.505);
158  TH2D* hptresmc90 = new TH2D("hptresmc90", "hptresmc90", 300, -1400, 1600, 200, -1, 1);
159  TH1D* hphiresmc80 = new TH1D("hphiresmc80", "hphiresmc80", 701, -3.505, 3.505);
160  TH2D* hptresmc80 = new TH2D("hptresmc80", "hptresmc80", 300, -1400, 1600, 200, -1, 1);
161  TH1D* hphiresmcbad = new TH1D("hphiresmcbad", "hphiresmcbad", 701, -3.505, 3.505);
162  TH2D* hptresmcbad = new TH2D("hptresmcbad", "hptresmcbad", 300, -1400, 1600, 200, -1, 1);
163 
164 
165  for (current_time = 0; current_time < final_time; current_time += delta_t) {
166  online->Clear();
167  //online->LoadHits( PndOnlineManager::kHESRRevolution );
168  online->ClearTracks();
169  online->LoadHits( delta_t );
170  //online->LoadHits( current_time );
171  online->Process();
172  //online->Clear();
173 
174  online_tracks = online->GetTrackObjectList();
175  int maxmcid = 0;
176  int minmcid = 0;
177  int firsti = 0;
178  if (online_tracks->GetEntriesFast() > 0) {
179  PndOnlineTrack* firstTrack = (PndOnlineTrack*) online_tracks->At(firsti);
180  while (firstTrack->DrawMode() != PndOnlineTrack::kVertexMomentumAsCircleHack) {
181  firsti++;
182  if (firsti > online_tracks->GetEntriesFast()) {
183  cout << "No drawable online Track in Array. We have a severe problem." << endl;
184  break;
185  }
186  firstTrack = (PndOnlineTrack*) online_tracks->At(firsti);
187  }
188  maxmcid = firstTrack->MCEntryNumber();
189  minmcid = firstTrack->MCEntryNumber();
190  }
191  for (int i = 0; i < online_tracks->GetEntriesFast(); i++) {
192  PndOnlineTrack* myTrack = (PndOnlineTrack*) online_tracks->At(i);
193  if (myTrack->DrawMode() != PndOnlineTrack::kVertexMomentumAsCircleHack) {
194  continue;
195  }
196  totalrecotracks++;
197 
198  mctree->GetEntry(myTrack->MCEntryNumber());
199  //cout << "MC Entry count: " << mcarray->GetEntriesFast() << endl;
200  PndMCTrack* mctrack = (PndMCTrack*)mcarray->At(myTrack->MCTrackNumber());
201  if (mctrack == 0) {
202  cout << "MC Track NULL pointer in event " << myTrack->MCEntryNumber() << ", track " << myTrack->MCTrackNumber() << endl;
203  cout << "Event Majority: " << myTrack->MCEntryMajority() << ", Track Majority: " << myTrack->MCTrackMajority() << endl;
204  recotracknullpointers++;
205  continue;
206  }
207  if (mctrack->GetMotherID() != -1) {
208  //continue;
209  }
210 
211  TVector2 innerhit = TVector2(myTrack->Vertex().x(), myTrack->Vertex().y());
212  TVector2 outerhit = TVector2(myTrack->Momentum().x(), myTrack->Momentum().y());
213 
214 // TVector2 circleorigin;
215 // Double_t circleradius;
216 // Int_t clockwise;
217 // CalculateCircle(&innerhit, &outerhit, &circleorigin, &circleradius, &clockwise);
218 // Double_t orthophi = circleorigin.Phi();
219 // Double_t trackphi = orthophi + 0.5*TMath::Pi()*clockwise;
220 // if (trackphi < 0) {
221 // trackphi += 2*TMath::Pi();
222 // }
223 // if (trackphi > 2*TMath::Pi()) {
224 // trackphi -= 2*TMath::Pi();
225 // }
226 
227  TVector2 circleorigin = myTrack->CircleOrigin();
228  Double_t orfoundrecotrackallthophi = circleorigin.Phi();
229  Double_t circleradius = myTrack->CircleRadius();
230  Int_t clockwise = myTrack->CircleClockwise();
231  Double_t trackphi = myTrack->CirclePhi();
232  Double_t trackpt = circleradius * 2.0 * 2.99792458;
233 /*
234  cout << "Origin: " << endl;
235  circleorigin.Print();
236  myTrack->CircleOrigin().Print();
237  cout << "Radius1: " << circleradius << endl;
238  cout << "Radius2: " << myTrack->CircleRadius() << endl;
239  cout << "Phi1: " << trackphi << endl;
240  cout << "Phi2: " << myTrack->CirclePhi() << endl;
241  cout << "Clockwise1: " << clockwise << endl;
242  cout << "Clockwise2: " << myTrack->CircleClockwise() << endl;
243 */
244 
245 // cout << "Reco Track MCEntryNumber: " << myTrack->MCEntryNumber() << endl;
246 // cout << "Reco Track MCEntryMajority: " << myTrack->MCEntryMajority() << endl;
247 // cout << "Reco Track MCTrackNumber: " << myTrack->MCTrackNumber() << endl;
248 // cout << "Reco Track MCTrackMajority: " << myTrack->MCTrackMajority() << endl;
249 // cout << "Reco Track HitCount: " << myTrack->HitCount() << endl;
250 
251  //compare reconstructed tracks with the MC input
252  //full check of all MC input data is done afterwards in a separate step
253 
254  Double_t mcphi = mctrack->GetMomentum().Phi();
255  Double_t mcpt = mctrack->GetMomentum().Pt()*1000;
256  if (mcphi < 0) {
257  mcphi += 2*TMath::Pi();
258  }
259  Double_t calcphimc = mcphi;
260  if (calcphimc < 0) calcphimc += 2*TMath::Pi();
261  Double_t calcphireco = trackphi;
262  if (calcphireco < 0) calcphireco += 2*TMath::Pi();
263  Double_t deltaphi = calcphireco - calcphimc;
264  if (deltaphi < 0) {
265  deltaphi *= -1;
266  }
267  if (deltaphi > TMath::Pi()) {
268  deltaphi = 2*TMath::Pi() - deltaphi;
269  }
270  Double_t deltapt = trackpt - mcpt;
271  Int_t chargesign = 1;
272  if ( (mctrack->GetMotherID() == -1) && (pdg->GetParticle(mctrack->GetPdgCode())->Charge() < 0) ) {
273  chargesign = -1;
274  }
275  /*if ( (deltaphi > 2.9) && (deltaphi < 3.2) ) {
276  cout << "Reco Track MCEntryNumber: " << myTrack->MCEntryNumber() << endl;
277  cout << "Reco Track MCEntryMajority: " << myTrack->MCEntryMajority() << endl;
278  cout << "Reco Track MCTrackNumber: " << myTrack->MCTrackNumber() << endl;
279  cout << "Reco Track MCTrackMajority: " << myTrack->MCTrackMajority() << endl;
280  cout << "Reco Track HitCount: " << myTrack->HitCount() << endl;
281  cout << "Reco Raw Phi: " << trackphi << endl;
282  cout << "Reco Phi: " << calcphireco << endl;
283  cout << "MC Comp Phi: " << calcphimc << endl;
284  cout << "MC Raw Phi: " << mcphi << endl;
285  cout << "Orthophi: " << orthophi << endl;
286  mctrack->GetMomentum().Print();
287  cout << "Inner: " << endl;
288  innerhit.Print();
289  cout << "Outer: " << endl;
290  outerhit.Print();
291  } else {
292  //cout << "Delta Phi: " << deltaphi << endl;
293  }*/
294  hphiresmcall->Fill(deltaphi*chargesign);
295  hptresmcall->Fill(mcpt, deltapt/mcpt);
296  int mceventid = myTrack->MCEntryNumber();
297  int mctrackid = myTrack->MCTrackNumber();
298  foundmapall.Add(mceventid, mctrackid);
299  if ( (myTrack->MCEntryMajority() > 0.9) && (myTrack->MCTrackMajority() > 0.9) ) {
300  mcidgood90++;
301  foundmap90.Add(mceventid, mctrackid);
302  hphiresmc90->Fill(deltaphi*chargesign);
303  hptresmc90->Fill(mcpt, deltapt/mcpt);
304  }
305  if ( (myTrack->MCEntryMajority() > 0.8) && (myTrack->MCTrackMajority() > 0.8) ) {
306  mcidgood80++;
307  foundmap80.Add(mceventid, mctrackid);
308  hphiresmc80->Fill(deltaphi*chargesign);
309  hptresmc80->Fill(mcpt, deltapt/mcpt);
310  } else {
311  mcidbad++;
312  foundmapbad.Add(mceventid, mctrackid);
313  hphiresmcbad->Fill(deltaphi*chargesign);
314  hptresmcbad->Fill(mcpt, deltapt/mcpt);
315  }
316 
317  //bool isrecoable = true;
318  //if (mctrack->GetMomentum().Pt() < 0.052) isrecoable = false;
319  if ( (mctrack->GetMomentum().Pt() < 0.052) && (mctrack->GetNPoints(kSTT) > 2) ) {
320  cout << "Low pt hit" << endl;
321  }
322 
323  int mcentrynumber = myTrack->MCEntryNumber();
324  if (mcentrynumber > maxmcid) {
325  maxmcid = mcentrynumber;
326  }
327  if (mcentrynumber < minmcid) {
328  minmcid = mcentrynumber;
329  }
330 
331  //std::cout << "RecoTrack: " << totalrecotracks << ", MCEntryNumber: " << mcentrynumber << " DigiTime: " << myTrack->T0Min() << ", T0time: " << myTrack->T0() << ", Pt " << trackpt << ", Phi: " << trackphi << " TriPhi: " << innerhit.Phi() << ", Radius: " << circleradius << ", DeltaT0: " << myTrack->T0Max()-myTrack->T0Min() << endl;
332 
333  //std::cout << "RecoTrack: " << totalrecotracks << ", MCEntryNumber: " << mcentrynumber << ", Pt " << trackpt << ", Phi: " << trackphi << " TriPhi: " << innerhit.Phi() << ", Radius: " << circleradius << endl;
334 
335  //recotracktime.push_back(myTrack->T0Min());
336  recotrackpt.push_back(trackpt);
337  recotrackphi.push_back(trackphi);
338  recotrackstatus.push_back(0);
339  recotrackmcid.push_back(mcentrynumber);
340 
341 
342  }
343  if (maxmcid > globalmaxmcid) {
344  globalmaxmcid = maxmcid;
345  }
346  cout << "MC ID Range: " << minmcid << ", " << maxmcid << endl;
347  //loop over mc events
348  for (int j = minmcid; j <= maxmcid; j++) {
349  cout << "MC outer loop: " << j << endl;
350  if (j > globalmaxmcid) {
351  break;
352  }
353  mctree->GetEntry(j);
354  cout << "MC Entry count: " << mcarray->GetEntriesFast() << endl;
355  //loop over mc tracks
356  for (int k = 0; k < mcarray->GetEntriesFast(); k++) {
357  mctrackgrandtotal++;
358  PndMCTrack* mctrack = (PndMCTrack*)mcarray->At(k);
359  if (mctrack == 0) {
360  cout << "MC Track has NULL pointer. Oooops." << endl;
361  mctracknullpointers++;
362  continue;
363  }
364  hmcall->Fill(mctrack->GetMomentum().Pz(), mctrack->GetMomentum().Pt());
365  if (mctrack->GetPdgCode() > 1000000000) {
366  cout << "PDG Code too large: " << mctrack->GetPdgCode() << endl;
367  //cout << "Mass: " << pdg->GetParticle(mctrack->GetPdgCode())->Mass() << endl;
368  pdgcodetoolarge++;
369  continue;
370  }
371  if (mctrack->GetMotherID() > 0) {
372  //cout << "Non primary." << endl;
373  mctracknonprime++;
374  continue;
375  }
376  hmcprime->Fill(mctrack->GetMomentum().Pz(), mctrack->GetMomentum().Pt());
377  if (mctrack->GetMomentum().Pt() < 0.06) {
378  mctracknonrecomomentum++;
379  if (mctrack->GetNPoints(kSTT) > 2) {
380  //cout << "Low pt hit" << endl;
381  lowpthit++;
382  }
383  continue;
384  }
385  if ( (mctrack->GetMomentum().Pt()/mctrack->GetMomentum().Pz() < 30.0/110.0) && (mctrack->GetMomentum().Pz() > 0) ) {
386  mctracknonrecomomentum++;
387  if (mctrack->GetNPoints(kSTT) > 0) {
388  //cout << "Low pt hit" << endl;
389  lowpthit++;
390  }
391  continue;
392  }
393  if ( (-(mctrack->GetMomentum().Pt())/mctrack->GetMomentum().Pz() < 30.0/55.0) && (mctrack->GetMomentum().Pz() < 0) ) {
394  mctracknonrecomomentum++;
395  if (mctrack->GetNPoints(kSTT) > 0) {
396  //cout << "Low pt hit" << endl;
397  lowpthit++;
398  }
399  continue;
400  }
401 
402  if ( (mctrack->GetMomentum().Mag())/(pdg->GetParticle(mctrack->GetPdgCode())->Mass()) < 0.4) ) {
403  cout << "Low betagamma particle" << endl;
404  lowbetagamma++;
405  continue;
406  }
407 
408  if (mctrack->GetNPoints(kSTT) < 3) {
409  mctracknonhit++;
410  cout << "Didn't hit." << endl;
411  continue;
412  }
413 
414  mctrackrecoattempt++;
415  hmcrecoable->Fill(mctrack->GetMomentum().Pz(), mctrack->GetMomentum().Pt());
416 
417 
418  if (foundmapall.Get(j, k) > 0) {
419  //cout << "Event " << j << ", Track " << k << " found " << foundmap90.Get(j, k) << " times." << endl;
420  foundrecotrackall++;
421  hmcrecoed->Fill(mctrack->GetMomentum().Pz(), mctrack->GetMomentum().Pt());
422  } else {
423  cout << "Did not find Event " << j << ", Track " << k << endl;
424  mctracknotrecoed++;
425  hmcnotrecoed->Fill(mctrack->GetMomentum().Pz(), mctrack->GetMomentum().Pt());
426  }
427  if (foundmap90.Get(j, k) > 0) {
428  //cout << "Event " << j << ", Track " << k << " found " << foundmap90.Get(j, k) << " times." << endl;
429  foundrecotrack90++;
430  }
431 
432  if (foundmap80.Get(j, k) > 0) {
433  //cout << "Event " << j << ", Track " << k << " found " << foundmap90.Get(j, k) << " times." << endl;
434  foundrecotrack80++;
435  }
436 
437  if (foundmapbad.Get(j, k) > 0) {
438  //cout << "Event " << j << ", Track " << k << " found " << foundmap90.Get(j, k) << " times." << endl;
439  foundrecotrackbad++;
440  }
441 
442  //legacy code below this point
443 // hmcall->Fill(mctrack->GetMomentum().Pz(), mctrack->GetMomentum().Pt());
444 // bool isrecoable = false;
445 // if (mctrack->GetMotherID() == -1) {
446 // hmcprime->Fill(mctrack->GetMomentum().Pz(), mctrack->GetMomentum().Pt());
447 // isrecoable = true;
448 // if (mctrack->GetMomentum().Pz() < -0.4) {
449 // cout << "Critical Track at Event: " << j << " Track " << k << " Particle: " << mctrack->GetPdgCode() << endl;
450 // //cout << "Particle Mass: " << pdg->GetParticle(mctrack->GetPdgCode())->Mass() << " Charge: " << pdg->GetParticle(mctrack->GetPdgCode())->Charge() << endl;
451 // //cout << "MCTrack STT Hits: " << mctrack->GetNPoints(kSTT) << ", Pz: " << mctrack->GetMomentum().Pz() << ", Pt: " << mctrack->GetMomentum().Pt() << ", Phi: " << mctrack->GetMomentum().Phi() << ", Radius: " << mctrack->GetMomentum().Pt() / (3.0*2.0) << endl;
452 // cout << "MCVertex: " << mctrack->GetStartVertex().Mag() << " XYZ: " << mctrack->GetStartVertex().x() << ", " << mctrack->GetStartVertex().y() << ", " << mctrack->GetStartVertex().z() << endl;
453 // }
454 // }
455  //if (mctrack->GetMotherID() == -1 && mctrack->GetNPoints(kSTT) > 0) {
456  //cout << TRho::Instance()->GetPDG()->GetParticle(mctrack->GetPdgCode())->Mass() << endl;
457  //}
458 
459  //bool isrecoable = true;
460 // if (mctrack->GetMomentum().Pt() < 0.052) isrecoable = false;
461 // if ( (mctrack->GetMomentum().Pt() < 0.052) && (mctrack->GetNPoints(kSTT) > 0) ) {
462 // cout << "Low pt hit" << endl;
463 // }
464 // if ( (mctrack->GetMomentum().Pt()/mctrack->GetMomentum().Pz() < 28.0/110.0) && (mctrack->GetMomentum().Pz() > 0) ) isrecoable = false;
465 // if ( (-(mctrack->GetMomentum().Pt())/mctrack->GetMomentum().Pz() < 28.0/55.0) && (mctrack->GetMomentum().Pz() < 0) ) isrecoable = false;
466 // if ( (mctrack->GetMotherID() == -1) && isrecoable && (mctrack->GetNPoints(kSTT) < 3) ) {
467 // if (TMath::Abs(pdg->GetParticle(mctrack->GetPdgCode())->Charge()) < 1) isrecoable = false;
468 // if (pdg->GetParticle(mctrack->GetPdgCode())->Mass() > 1) isrecoable = false;
469 // }
470 // if ( (mctrack->GetMotherID() == -1) && isrecoable && (mctrack->GetMomentum().Mag()/pdg->GetParticle(mctrack->GetPdgCode())->Mass() < 0.4) ) {
471 // cout << "Low betagamma particle" << endl;
472 // isrecoable = false;
473 // }
474 // if ( (mctrack->GetMotherID() == -1) && isrecoable && (mctrack->GetNPoints(kSTT) < 3) ) {
475 // cout << "Meh, reco possibility mismatch." << endl;
476 // hdafuq->Fill(mctrack->GetMomentum().Pz(), mctrack->GetMomentum().Pt());
477 // cout << "Dafuq at Event: " << j << " Track " << k << " Particle: " << mctrack->GetPdgCode() << " Charge: " << pdg->GetParticle(mctrack->GetPdgCode())->Charge() << endl;
478 // mctrack->GetMomentum().Print();
479 // isrecoable = false;
480 // }
481 
482  //if (mctrack->GetMotherID() == -1 && mctrack->GetNPoints(kSTT) > 0) {
483 // if (mctrack->GetMotherID() == -1 && isrecoable) {
484 // hmcrecoable->Fill(mctrack->GetMomentum().Pz(), mctrack->GetMomentum().Pt());
485 // Double_t mcphi = mctrack->GetMomentum().Phi();
486 // Double_t mcpt = mctrack->GetMomentum().Pt()*1000;
487 // if (mcphi < 0) {
488 // mcphi += 2*TMath::Pi();
489 // }
490 // //cout << "MCTrack: " << j << " " << k << ", Pt: " << mcpt << ", Phi: " << mcphi << ", Radius: " << mcpt / (3.0*2.0) << endl;
491 // //cout << "MCVertex: " << mctrack->GetStartVertex().Mag() << " XYZ: " << mctrack->GetStartVertex().x() << ", " << mctrack->GetStartVertex().y() << ", " << mctrack->GetStartVertex().z() << endl;
492 // int trackstatus = 0;
493 // if (mctrack->GetStartVertex().Mag() < 2) {
494 // for (int r = 0; r < recotrackmcid.size(); r++) {
495 // if (recotrackmcid.at(r) == j) {
496 // //Double_t deltaphi = recotrackphi.at(r) - mcphi;
497 // Double_t calcphimc = mcphi;
498 // if (calcphimc < 0) calcphimc += 2*TMath::Pi();
499 // Double_t calcphireco = recotrackphi.at(r);
500 // if (calcphireco < 0) calcphireco += 2*TMath::Pi();
501 // Double_t deltaphi = calcphireco - calcphimc;
502 // if (deltaphi < 0) {
503 // deltaphi *= -1;
504 // }
505 // if (deltaphi > TMath::Pi()) {
506 // deltaphi = 2*TMath::Pi() - deltaphi;
507 // }
508 // Double_t deltapt = recotrackpt.at(r) - mcpt;
509 //
510 // Int_t chargesign = 1;
511 // if (pdg->GetParticle(mctrack->GetPdgCode())->Charge() < 0) {
512 // chargesign = -1;
513 // }
514 //
515 // hphires->Fill(deltaphi*chargesign);
516 // hptres->Fill(mcpt, deltapt/mcpt);
517 //
518 // if ( (deltapt/mcpt < -0.75) && (deltapt/mcpt > -0.8) && (mcpt > 300) && (mcpt < 650) ) {
519 // cout << "Funny resolution at Event: " << j << " Track " << k << " Particle: " << mctrack->GetPdgCode() << " Charge: " << pdg->GetParticle(mctrack->GetPdgCode())->Charge() << endl;
520 // mctrack->GetMomentum().Print();
521 // cout << "MC pt: " << mcpt << endl;
522 // cout << "Re pt: " << recotrackpt.at(r) << endl;
523 // }
524 //
525 // bool phiok = (deltaphi < 0.1);
526 //
527 // if (deltapt < 0) {
528 // deltapt *= -1;
529 // }
530 // bool ptok = ( (deltapt < mcpt * 0.1) || ((mcpt > 450)&&(recotrackpt.at(r) > 450)) );
531 // if (phiok && ptok) {
532 // trackstatus++;
533 // recotrackstatus.at(r) = 1;
534 // //hmcrecoed->Fill(mctrack->GetMomentum().Pz(), mctrack->GetMomentum().Pt());
535 // } else {
536 // //hmcnotrecoed->Fill(mctrack->GetMomentum().Pz(), mctrack->GetMomentum().Pt());
537 // }
538 // if ((!phiok) && ptok) {
539 // deltaphifailed++;
540 // }
541 // if (phiok && (!ptok)) {
542 // deltaptfailed++;
543 // }
544 // if ((!phiok) && (!ptok)) {
545 // deltabothfailed++;
546 // }
547 // }
548 // }
549 // //cout << "Found " << trackstatus << " matching reco tracks." << endl;
550 // if (trackstatus == 0) {
551 // hmcnotrecoed->Fill(mctrack->GetMomentum().Pz(), mctrack->GetMomentum().Pt());
552 // norecotracks++;
553 // } else {
554 // hmcrecoed->Fill(mctrack->GetMomentum().Pz(), mctrack->GetMomentum().Pt());
555 // }
556 // mctrackstatus.push_back(trackstatus);
557 // goodmctracks++;
558 // } else {
559 // secondaries++;
560 // }
561 // }
562  }
563  }
564  }
565 
566  for (int s = 0; s < mctrackstatus.size(); s++) {
567  if (mctrackstatus.at(s) > 0) {
568  goodrecotracks++;
569  }
570  }
571 
572  int recotrackstatusgoodcount = 0;
573  int recotrackstatusbadcount = 0;
574 
575  for (int s = 0; s < recotrackstatus.size(); s++) {
576  if (recotrackmcid.at(s) <= globalmaxmcid) {
577  if (recotrackstatus.at(s) > 0) {
578  recotrackstatusgoodcount++;
579  } else {
580  recotrackstatusbadcount++;
581  }
582  }
583  }
584 
585  //for (std::map<int, std::map<int, int> >::iterator iterevent = foundmap90.begin(); iterevent != foundmap90.end(); iterevent++) {
586  //for (itertrack = *iterevent.begin(); itertrack != *iterevent.end(); itertrack++
587  //}
588 
589  cout << "Total reconstructed track count: " << totalrecotracks << endl;
590  cout << "Reached MC ID: " << globalmaxmcid << endl;
591  cout << "Reconstructable MC track count: " << goodmctracks << endl;
592  cout << "MC Tracks reconstructed: " << goodrecotracks << endl;
593  cout << "MC Tracks NOT reconstructed: " << norecotracks << endl;
594  cout << "Amount of reconstructed MC tracks: " << goodrecotrack << endl;
595  cout << "Amount of not reconstructed MC tracks: " << badrecotrack << endl;
596  cout << "Reco Status Good: " << recotrackstatusgoodcount << endl;
597  cout << "Reco Status Bad: " << recotrackstatusbadcount << endl;
598  cout << "Delta phi failed: " << deltaphifailed << endl;
599  cout << "Delta pt failed: " << deltaptfailed << endl;
600  cout << "Delta both failed: " << deltabothfailed << endl;
601  cout << "Secondary Tracks: " << secondaries << endl;
602 
603  cout << "MC Track Grand Total: " << mctrackgrandtotal << endl;
604  cout << "MC Track PDG Code too large: " << pdgcodetoolarge << endl;
605  cout << "MC Track Non Prime: " << mctracknonprime << endl;
606  cout << "MC Track Low Momentum: " << mctracknonrecomomentum << endl;
607  cout << " MC Track Low Pt Hit: " << lowpthit << endl;
608  cout << "MC Track Low BetaGamma: " << lowbetagamma << endl;
609  cout << "MC Track missed STT: " << mctracknonhit << endl;
610  cout << "MC Track to be recoed: " << mctrackrecoattempt << endl;
611  cout << "MC Track recoed (all): " << foundrecotrackall << endl;
612  cout << "MC Track recoed (90): " << foundrecotrack90 << endl;
613  cout << "MC Track recoed (80): " << foundrecotrack80 << endl;
614  cout << "MC Track recoed (bad): " << foundrecotrackbad << endl;
615  cout << "MC Track not recoed: " << mctracknotrecoed << endl;
616 
617 
618  cout << "MC Match Good 90: " << mcidgood90 << endl;
619  cout << "MC Match Good 80: " << mcidgood80 << endl;
620  cout << "MC Match Bad: " << mcidbad << endl;
621  cout << "MC Tracks with NULL pointers in Reco: " << recotracknullpointers << endl;
622  cout << "MC Tracks with NULL pointers: " << mctracknullpointers << endl;
623 
624  //c1->Print("hitdisplay.gif++");
625 
626  TFile outputstorage(histofilename, "UPDATE");
627  outputstorage.cd();
628  SaveAndUpdateHisto(hmcall, outputstorage);
629  SaveAndUpdateHisto(hmcprime, outputstorage);
630  SaveAndUpdateHisto(hmcrecoable, outputstorage);
631  SaveAndUpdateHisto(hdafuq, outputstorage);
632  SaveAndUpdateHisto(hmcrecoed, outputstorage);
633  SaveAndUpdateHisto(hmcnotrecoed, outputstorage);
634  SaveAndUpdateHisto(hphires, outputstorage);
635  SaveAndUpdateHisto(hptres, outputstorage);
636  SaveAndUpdateHisto(hphiresmcall, outputstorage);
637  SaveAndUpdateHisto(hptresmcall, outputstorage);
638  SaveAndUpdateHisto(hphiresmc90, outputstorage);
639  SaveAndUpdateHisto(hptresmc90, outputstorage);
640  SaveAndUpdateHisto(hphiresmc80, outputstorage);
641  SaveAndUpdateHisto(hptresmc80, outputstorage);
642  SaveAndUpdateHisto(hphiresmcbad, outputstorage);
643  SaveAndUpdateHisto(hptresmcbad, outputstorage);
644 
645  // ----- Finish -------------------------------------------------------
646  timer.Stop();
647  Double_t rtime = timer.RealTime();
648  Double_t ctime = timer.CpuTime();
649  cout << endl << endl;
650  cout << "Macro finished succesfully." << endl;
651  cout << "Real time " << rtime << " s, CPU time " << ctime << " s" << endl;
652  cout << endl;
653  // ------------------------------------------------------------------------
654  return 0;
655 }
mctrack
PndMCTrack * mctrack
Definition: analysis_point_hit_fwendcap.C:250
filereco
TFile filereco("MvdStt_Test_reco.root")
kMVD
Definition: production/tdrct/eta_c/track_check.C:4
kTOF
Definition: production/tdrct/eta_c/track_check.C:4
i
Int_t i
Definition: run_full.C:25
kHYP
Definition: production/tdrct/eta_c/track_check.C:4
PndMCTrack::GetNPoints
Int_t GetNPoints(DetectorId detId) const
Definition: PndMCTrack.cxx:120
PndMCTrack::GetPdgCode
Int_t GetPdgCode() const
Definition: PndMCTrack.h:73
s
TLorentzVector s
Definition: Pnd2DStar.C:50
kLUMI
Definition: production/tdrct/eta_c/track_check.C:4
kTPC
Definition: production/tdrct/eta_c/track_check.C:4
PndMCTrack::GetMomentum
TVector3 GetMomentum() const
Definition: PndMCTrack.h:78
kSTT
Definition: production/tdrct/eta_c/track_check.C:4
fRun
FairRunAna * fRun
Definition: hit_dirc.C:58
kDSK
Definition: production/tdrct/eta_c/track_check.C:4
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
kDRC
Definition: production/tdrct/eta_c/track_check.C:4
DetectorId
DetectorId
Definition: production/tdrct/eta_c/track_check.C:2
treedigi
TTree * treedigi
Definition: checkdedx_helixhit.C:19
timer
TStopwatch timer
Definition: hit_dirc.C:51
filedigi
TFile filedigi("testdigi.root")
kMDT
Definition: production/tdrct/eta_c/track_check.C:4
rtdb
FairRuntimeDb * rtdb
Definition: hit_dirc.C:66
kRPC
Definition: production/tdrct/eta_c/track_check.C:4
rootlogon
rootlogon()
Definition: outdated/mpiTools/macros/emc/rootlogon.C:1
kHYPG
Definition: production/tdrct/eta_c/track_check.C:4
c1
c1
Definition: plot_dirc.C:35
ctime
Double_t ctime
Definition: hit_dirc.C:114
kEMC
Definition: production/tdrct/eta_c/track_check.C:4
iVerbose
Int_t iVerbose
Definition: Lars/runMvdDigi.C:7
kGEM
Definition: production/tdrct/eta_c/track_check.C:4
runOnlineDisplayMCCheckFaster2
int runOnlineDisplayMCCheckFaster2(Int_t maximumTime=2050)
Definition: runOnlineDisplayMCCheckFaster2.C:12
PndMCTrack::GetMotherID
Int_t GetMotherID() const
Definition: PndMCTrack.h:74
rtime
Double_t rtime
Definition: hit_dirc.C:113
Pi
Double_t Pi
Definition: PndGeoHypGeCrystal.cxx:10
kFTS
Definition: production/tdrct/eta_c/track_check.C:4
SaveAndUpdateHisto
void SaveAndUpdateHisto(TH1 *currenthisto, TFile &storagefile)
Definition: runOnlineDisplayMCCheck.C:3
kDCH
Definition: production/tdrct/eta_c/track_check.C:4
Generated on Wed Apr 3 2019 10:03:02 for FairRoot/PandaRoot by   doxygen 1.8.5