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