FairRoot/PandaRoot
Functions
runOnlineDisplayMCCheckFaster.C File Reference
#include "drawhelpers.hxx"

Go to the source code of this file.

Functions

void SaveAndUpdateHisto (TH1 *currenthisto, TFile &storagefile)
 
int runOnlineDisplayMCCheckFaster (Int_t maximumTime=2050)
 

Function Documentation

int runOnlineDisplayMCCheckFaster ( Int_t  maximumTime = 2050)

Definition at line 12 of file runOnlineDisplayMCCheckFaster.C.

References CAMath::Abs(), c1, ctime, Double_t, filedigi(), filereco(), fRun, PndMCTrack::GetMomentum(), PndMCTrack::GetMotherID(), PndMCTrack::GetNPoints(), PndMCTrack::GetPdgCode(), PndMCTrack::GetStartVertex(), i, iVerbose, kDCH, kDRC, kDSK, kEMC, kFTS, kGEM, kHYP, kHYPG, kLUMI, kMDT, kMVD, kRPC, kSTT, kTOF, kTPC, Mass, mctrack, Pi, r, rootlogon(), rtdb, rtime, s, SaveAndUpdateHisto(), timer, treedigi, and TString.

12  {
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 totalmctracks = 0;
105  int goodmctracks = 0;
106  int totalrecotracks = 0;
107  int goodrecotracks = 0;
108  int norecotracks = 0;
109  int deltaphifailed = 0;
110  int deltaptfailed = 0;
111  int deltabothfailed = 0;
112  int secondaries = 0;
113  int globalmaxmcid = 0;
114 
115  //std::map<int, std::map<int, int> > foundmap90;
116  //std::map<std::string, int> foundmap90;
117  //std::map<int, std::map<int, int> > foundmap80;
118  //std::map<int, std::map<int, int> > foundmapbad;
119 
120  int mctracknullpointers = 0;
121  int recotracknullpointers = 0;
122  int mcidgood90 = 0;
123  int mcidgood80 = 0;
124  int mcidbad = 0;
125 
126  //output histos
127  TH2D* hmcall = new TH2D("hmcall", "hmcall", 1500, -14, 16, 75, 0, 1.5);
128  TH2D* hmcprime = new TH2D("hmcprime", "hmcprime", 1500, -14, 16, 75, 0, 1.5);
129  TH2D* hmcrecoable = new TH2D("hmcrecoable", "hmcrecoable", 1500, -14, 16, 75, 0, 1.5);
130  TH2D* hdafuq = new TH2D("hdafuq", "hdafuq", 1500, -14, 16, 75, 0, 1.5);
131  TH2D* hmcrecoed = new TH2D("hmcrecoed", "hmcrecoed", 1500, -14, 16, 75, 0, 1.5);
132  TH2D* hmcnotrecoed = new TH2D("hmcnotrecoed", "hmcnotrecoed", 1500, -14, 16, 75, 0, 1.5);
133 
134  TH1D* hphires = new TH1D("hphires", "hphires", 701, -3.505, 3.505);
135  TH2D* hptres = new TH2D("hptres", "hptres", 300, -1400, 1600, 200, -1, 1);
136  TH1D* hphiresmc = new TH1D("hphiresmc", "hphiresmc", 701, -3.505, 3.505);
137  TH2D* hptresmc = new TH2D("hptresmc", "hptresmc", 300, -1400, 1600, 200, -1, 1);
138  TH1D* hphiresmc90 = new TH1D("hphiresmc90", "hphiresmc90", 701, -3.505, 3.505);
139  TH2D* hptresmc90 = new TH2D("hptresmc90", "hptresmc90", 300, -1400, 1600, 200, -1, 1);
140  TH1D* hphiresmc80 = new TH1D("hphiresmc80", "hphiresmc80", 701, -3.505, 3.505);
141  TH2D* hptresmc80 = new TH2D("hptresmc80", "hptresmc80", 300, -1400, 1600, 200, -1, 1);
142  TH1D* hphiresmcbad = new TH1D("hphiresmcbad", "hphiresmcbad", 701, -3.505, 3.505);
143  TH2D* hptresmcbad = new TH2D("hptresmcbad", "hptresmcbad", 300, -1400, 1600, 200, -1, 1);
144 
145 
146  for (current_time = 0; current_time < final_time; current_time += delta_t) {
147  online->Clear();
148  //online->LoadHits( PndOnlineManager::kHESRRevolution );
149  online->ClearTracks();
150  online->LoadHits( delta_t );
151  //online->LoadHits( current_time );
152  online->Process();
153  //online->Clear();
154 
155  online_tracks = online->GetTrackObjectList();
156  int maxmcid = 0;
157  int minmcid = 0;
158  int firsti = 0;
159  if (online_tracks->GetEntriesFast() > 0) {
160  PndOnlineTrack* firstTrack = (PndOnlineTrack*) online_tracks->At(firsti);
161  while (firstTrack->DrawMode() != PndOnlineTrack::kVertexMomentumAsCircleHack) {
162  firsti++;
163  if (firsti > online_tracks->GetEntriesFast()) {
164  cout << "No drawable online Track in Array. We have a severe problem." << endl;
165  break;
166  }
167  firstTrack = (PndOnlineTrack*) online_tracks->At(firsti);
168  }
169  maxmcid = firstTrack->MCEntryNumber();
170  minmcid = firstTrack->MCEntryNumber();
171  }
172  for (int i = 0; i < online_tracks->GetEntriesFast(); i++) {
173  PndOnlineTrack* myTrack = (PndOnlineTrack*) online_tracks->At(i);
174  if (myTrack->DrawMode() != PndOnlineTrack::kVertexMomentumAsCircleHack) {
175  continue;
176  }
177  totalrecotracks++;
178 
179  mctree->GetEntry(myTrack->MCEntryNumber());
180  //cout << "MC Entry count: " << mcarray->GetEntriesFast() << endl;
181  PndMCTrack* mctrack = (PndMCTrack*)mcarray->At(myTrack->MCTrackNumber());
182  if (mctrack == 0) {
183  cout << "MC Track NULL pointer in event " << myTrack->MCEntryNumber() << ", track " << myTrack->MCTrackNumber() << endl;
184  cout << "Event Majority: " << myTrack->MCEntryMajority() << ", Track Majority: " << myTrack->MCTrackMajority() << endl;
185  recotracknullpointers++;
186  continue;
187  }
188  if (mctrack->GetMotherID() != -1) {
189  //continue;
190  }
191 
192  TVector2 innerhit = TVector2(myTrack->Vertex().x(), myTrack->Vertex().y());
193  TVector2 outerhit = TVector2(myTrack->Momentum().x(), myTrack->Momentum().y());
194 
195 // TVector2 circleorigin;
196 // Double_t circleradius;
197 // Int_t clockwise;
198 // CalculateCircle(&innerhit, &outerhit, &circleorigin, &circleradius, &clockwise);
199 // Double_t orthophi = circleorigin.Phi();
200 // Double_t trackphi = orthophi + 0.5*TMath::Pi()*clockwise;
201 // if (trackphi < 0) {
202 // trackphi += 2*TMath::Pi();
203 // }
204 // if (trackphi > 2*TMath::Pi()) {
205 // trackphi -= 2*TMath::Pi();
206 // }
207 
208  TVector2 circleorigin = myTrack->CircleOrigin();
209  Double_t orthophi = circleorigin.Phi();
210  Double_t circleradius = myTrack->CircleRadius();
211  Int_t clockwise = myTrack->CircleClockwise();
212  Double_t trackphi = myTrack->CirclePhi();
213  Double_t trackpt = circleradius * 2.0 * 2.99792458;
214 /*
215  cout << "Origin: " << endl;
216  circleorigin.Print();
217  myTrack->CircleOrigin().Print();
218  cout << "Radius1: " << circleradius << endl;
219  cout << "Radius2: " << myTrack->CircleRadius() << endl;
220  cout << "Phi1: " << trackphi << endl;
221  cout << "Phi2: " << myTrack->CirclePhi() << endl;
222  cout << "Clockwise1: " << clockwise << endl;
223  cout << "Clockwise2: " << myTrack->CircleClockwise() << endl;
224 */
225 
226 // cout << "Reco Track MCEntryNumber: " << myTrack->MCEntryNumber() << endl;
227 // cout << "Reco Track MCEntryMajority: " << myTrack->MCEntryMajority() << endl;
228 // cout << "Reco Track MCTrackNumber: " << myTrack->MCTrackNumber() << endl;
229 // cout << "Reco Track MCTrackMajority: " << myTrack->MCTrackMajority() << endl;
230 // cout << "Reco Track HitCount: " << myTrack->HitCount() << endl;
231 
232  //compare reconstructed tracks with the MC input
233  //full check of all MC input data is done afterwards in a separate step
234 
235  Double_t mcphi = mctrack->GetMomentum().Phi();
236  Double_t mcpt = mctrack->GetMomentum().Pt()*1000;
237  if (mcphi < 0) {
238  mcphi += 2*TMath::Pi();
239  }
240  Double_t calcphimc = mcphi;
241  if (calcphimc < 0) calcphimc += 2*TMath::Pi();
242  Double_t calcphireco = trackphi;
243  if (calcphireco < 0) calcphireco += 2*TMath::Pi();
244  Double_t deltaphi = calcphireco - calcphimc;
245  if (deltaphi < 0) {
246  deltaphi *= -1;
247  }
248  if (deltaphi > TMath::Pi()) {
249  deltaphi = 2*TMath::Pi() - deltaphi;
250  }
251  Double_t deltapt = trackpt - mcpt;
252  Int_t chargesign = 1;
253  if ( (mctrack->GetMotherID() == -1) && (pdg->GetParticle(mctrack->GetPdgCode())->Charge() < 0) ) {
254  chargesign = -1;
255  }
256  /*if ( (deltaphi > 2.9) && (deltaphi < 3.2) ) {
257  cout << "Reco Track MCEntryNumber: " << myTrack->MCEntryNumber() << endl;
258  cout << "Reco Track MCEntryMajority: " << myTrack->MCEntryMajority() << endl;
259  cout << "Reco Track MCTrackNumber: " << myTrack->MCTrackNumber() << endl;
260  cout << "Reco Track MCTrackMajority: " << myTrack->MCTrackMajority() << endl;
261  cout << "Reco Track HitCount: " << myTrack->HitCount() << endl;
262  cout << "Reco Raw Phi: " << trackphi << endl;
263  cout << "Reco Phi: " << calcphireco << endl;
264  cout << "MC Comp Phi: " << calcphimc << endl;
265  cout << "MC Raw Phi: " << mcphi << endl;
266  cout << "Orthophi: " << orthophi << endl;
267  mctrack->GetMomentum().Print();
268  cout << "Inner: " << endl;
269  innerhit.Print();
270  cout << "Outer: " << endl;
271  outerhit.Print();
272  } else {
273  //cout << "Delta Phi: " << deltaphi << endl;
274  }*/
275  hphiresmc->Fill(deltaphi*chargesign);
276  hptresmc->Fill(mcpt, deltapt/mcpt);
277  int mceventid = myTrack->MCEntryNumber();
278  int mctrackid = myTrack->MCTrackNumber();
279  if ( (myTrack->MCEntryMajority() > 0.9) && (myTrack->MCTrackMajority() > 0.9) ) {
280  mcidgood90++;
281  hphiresmc90->Fill(deltaphi*chargesign);
282  hptresmc90->Fill(mcpt, deltapt/mcpt);
283  }
284  if ( (myTrack->MCEntryMajority() > 0.8) && (myTrack->MCTrackMajority() > 0.8) ) {
285  mcidgood80++;
286  hphiresmc80->Fill(deltaphi*chargesign);
287  hptresmc80->Fill(mcpt, deltapt/mcpt);
288  } else {
289  mcidbad++;
290  hphiresmcbad->Fill(deltaphi*chargesign);
291  hptresmcbad->Fill(mcpt, deltapt/mcpt);
292  }
293 
294  //bool isrecoable = true;
295  //if (mctrack->GetMomentum().Pt() < 0.052) isrecoable = false;
296  if ( (mctrack->GetMomentum().Pt() < 0.052) && (mctrack->GetNPoints(kSTT) > 0) ) {
297  cout << "Low pt hit" << endl;
298  }
299 
300  int mcentrynumber = myTrack->MCEntryNumber();
301  if (mcentrynumber > maxmcid) {
302  maxmcid = mcentrynumber;
303  }
304  if (mcentrynumber < minmcid) {
305  minmcid = mcentrynumber;
306  }
307 
308  //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;
309 
310  //std::cout << "RecoTrack: " << totalrecotracks << ", MCEntryNumber: " << mcentrynumber << ", Pt " << trackpt << ", Phi: " << trackphi << " TriPhi: " << innerhit.Phi() << ", Radius: " << circleradius << endl;
311 
312  //recotracktime.push_back(myTrack->T0Min());
313  recotrackpt.push_back(trackpt);
314  recotrackphi.push_back(trackphi);
315  recotrackstatus.push_back(0);
316  recotrackmcid.push_back(mcentrynumber);
317 
318 
319  }
320  if (maxmcid > globalmaxmcid) {
321  globalmaxmcid = maxmcid;
322  }
323  cout << "MC ID Range: " << minmcid << ", " << maxmcid << endl;
324  for (int j = minmcid; j <= maxmcid; j++) {
325  cout << "MC outer loop: " << j << endl;
326  if (j > globalmaxmcid) {
327  break;
328  }
329  mctree->GetEntry(j);
330  cout << "MC Entry count: " << mcarray->GetEntriesFast() << endl;
331  for (int k = 0; k < mcarray->GetEntriesFast(); k++) {
332  PndMCTrack* mctrack = (PndMCTrack*)mcarray->At(k);
333  if (mctrack == 0) {
334  cout << "MC Track has NULL pointer. Oooops." << endl;
335  mctracknullpointers++;
336  continue;
337  }
338  hmcall->Fill(mctrack->GetMomentum().Pz(), mctrack->GetMomentum().Pt());
339  bool isrecoable = false;
340  if (mctrack->GetMotherID() == -1) {
341  hmcprime->Fill(mctrack->GetMomentum().Pz(), mctrack->GetMomentum().Pt());
342  isrecoable = true;
343  if (mctrack->GetMomentum().Pz() < -0.4) {
344  cout << "Critical Track at Event: " << j << " Track " << k << " Particle: " << mctrack->GetPdgCode() << endl;
345  //cout << "Particle Mass: " << pdg->GetParticle(mctrack->GetPdgCode())->Mass() << " Charge: " << pdg->GetParticle(mctrack->GetPdgCode())->Charge() << endl;
346  //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;
347  cout << "MCVertex: " << mctrack->GetStartVertex().Mag() << " XYZ: " << mctrack->GetStartVertex().x() << ", " << mctrack->GetStartVertex().y() << ", " << mctrack->GetStartVertex().z() << endl;
348  }
349  }
350  //if (mctrack->GetMotherID() == -1 && mctrack->GetNPoints(kSTT) > 0) {
351  //cout << TRho::Instance()->GetPDG()->GetParticle(mctrack->GetPdgCode())->Mass() << endl;
352  //}
353 
354  //bool isrecoable = true;
355  if (mctrack->GetMomentum().Pt() < 0.052) isrecoable = false;
356  if ( (mctrack->GetMomentum().Pt() < 0.052) && (mctrack->GetNPoints(kSTT) > 0) ) {
357  cout << "Low pt hit" << endl;
358  }
359  if ( (mctrack->GetMomentum().Pt()/mctrack->GetMomentum().Pz() < 28.0/110.0) && (mctrack->GetMomentum().Pz() > 0) ) isrecoable = false;
360  if ( (-(mctrack->GetMomentum().Pt())/mctrack->GetMomentum().Pz() < 28.0/55.0) && (mctrack->GetMomentum().Pz() < 0) ) isrecoable = false;
361  if ( (mctrack->GetMotherID() == -1) && isrecoable && (mctrack->GetNPoints(kSTT) < 3) ) {
362  if (TMath::Abs(pdg->GetParticle(mctrack->GetPdgCode())->Charge()) < 1) isrecoable = false;
363  if (pdg->GetParticle(mctrack->GetPdgCode())->Mass() > 1) isrecoable = false;
364  }
365  if ( (mctrack->GetMotherID() == -1) && isrecoable && (mctrack->GetMomentum().Mag()/pdg->GetParticle(mctrack->GetPdgCode())->Mass() < 0.4) ) {
366  cout << "Low betagamma particle" << endl;
367  isrecoable = false;
368  }
369  if ( (mctrack->GetMotherID() == -1) && isrecoable && (mctrack->GetNPoints(kSTT) < 3) ) {
370  cout << "Meh, reco possibility mismatch." << endl;
371  hdafuq->Fill(mctrack->GetMomentum().Pz(), mctrack->GetMomentum().Pt());
372  cout << "Dafuq at Event: " << j << " Track " << k << " Particle: " << mctrack->GetPdgCode() << " Charge: " << pdg->GetParticle(mctrack->GetPdgCode())->Charge() << endl;
373  mctrack->GetMomentum().Print();
374  isrecoable = false;
375  }
376 
377  //if (mctrack->GetMotherID() == -1 && mctrack->GetNPoints(kSTT) > 0) {
378  if (mctrack->GetMotherID() == -1 && isrecoable) {
379  hmcrecoable->Fill(mctrack->GetMomentum().Pz(), mctrack->GetMomentum().Pt());
380  Double_t mcphi = mctrack->GetMomentum().Phi();
381  Double_t mcpt = mctrack->GetMomentum().Pt()*1000;
382  if (mcphi < 0) {
383  mcphi += 2*TMath::Pi();
384  }
385  //cout << "MCTrack: " << j << " " << k << ", Pt: " << mcpt << ", Phi: " << mcphi << ", Radius: " << mcpt / (3.0*2.0) << endl;
386  //cout << "MCVertex: " << mctrack->GetStartVertex().Mag() << " XYZ: " << mctrack->GetStartVertex().x() << ", " << mctrack->GetStartVertex().y() << ", " << mctrack->GetStartVertex().z() << endl;
387  int trackstatus = 0;
388  if (mctrack->GetStartVertex().Mag() < 2) {
389  for (int r = 0; r < recotrackmcid.size(); r++) {
390  if (recotrackmcid.at(r) == j) {
391  //Double_t deltaphi = recotrackphi.at(r) - mcphi;
392  Double_t calcphimc = mcphi;
393  if (calcphimc < 0) calcphimc += 2*TMath::Pi();
394  Double_t calcphireco = recotrackphi.at(r);
395  if (calcphireco < 0) calcphireco += 2*TMath::Pi();
396  Double_t deltaphi = calcphireco - calcphimc;
397  if (deltaphi < 0) {
398  deltaphi *= -1;
399  }
400  if (deltaphi > TMath::Pi()) {
401  deltaphi = 2*TMath::Pi() - deltaphi;
402  }
403  Double_t deltapt = recotrackpt.at(r) - mcpt;
404 
405  Int_t chargesign = 1;
406  if (pdg->GetParticle(mctrack->GetPdgCode())->Charge() < 0) {
407  chargesign = -1;
408  }
409 
410  hphires->Fill(deltaphi*chargesign);
411  hptres->Fill(mcpt, deltapt/mcpt);
412 
413  if ( (deltapt/mcpt < -0.75) && (deltapt/mcpt > -0.8) && (mcpt > 300) && (mcpt < 650) ) {
414  cout << "Funny resolution at Event: " << j << " Track " << k << " Particle: " << mctrack->GetPdgCode() << " Charge: " << pdg->GetParticle(mctrack->GetPdgCode())->Charge() << endl;
415  mctrack->GetMomentum().Print();
416  cout << "MC pt: " << mcpt << endl;
417  cout << "Re pt: " << recotrackpt.at(r) << endl;
418  }
419 
420  bool phiok = (deltaphi < 0.1);
421 
422  if (deltapt < 0) {
423  deltapt *= -1;
424  }
425  bool ptok = ( (deltapt < mcpt * 0.1) || ((mcpt > 450)&&(recotrackpt.at(r) > 450)) );
426  if (phiok && ptok) {
427  trackstatus++;
428  recotrackstatus.at(r) = 1;
429  //hmcrecoed->Fill(mctrack->GetMomentum().Pz(), mctrack->GetMomentum().Pt());
430  } else {
431  //hmcnotrecoed->Fill(mctrack->GetMomentum().Pz(), mctrack->GetMomentum().Pt());
432  }
433  if ((!phiok) && ptok) {
434  deltaphifailed++;
435  }
436  if (phiok && (!ptok)) {
437  deltaptfailed++;
438  }
439  if ((!phiok) && (!ptok)) {
440  deltabothfailed++;
441  }
442  }
443  }
444  //cout << "Found " << trackstatus << " matching reco tracks." << endl;
445  if (trackstatus == 0) {
446  hmcnotrecoed->Fill(mctrack->GetMomentum().Pz(), mctrack->GetMomentum().Pt());
447  norecotracks++;
448  } else {
449  hmcrecoed->Fill(mctrack->GetMomentum().Pz(), mctrack->GetMomentum().Pt());
450  }
451  mctrackstatus.push_back(trackstatus);
452  goodmctracks++;
453  } else {
454  secondaries++;
455  }
456  }
457  }
458  }
459  }
460 
461  for (int s = 0; s < mctrackstatus.size(); s++) {
462  if (mctrackstatus.at(s) > 0) {
463  goodrecotracks++;
464  }
465  }
466 
467  int recotrackstatusgoodcount = 0;
468  int recotrackstatusbadcount = 0;
469 
470  for (int s = 0; s < recotrackstatus.size(); s++) {
471  if (recotrackmcid.at(s) <= globalmaxmcid) {
472  if (recotrackstatus.at(s) > 0) {
473  recotrackstatusgoodcount++;
474  } else {
475  recotrackstatusbadcount++;
476  }
477  }
478  }
479 
480  //for (std::map<int, std::map<int, int> >::iterator iterevent = foundmap90.begin(); iterevent != foundmap90.end(); iterevent++) {
481  //for (itertrack = *iterevent.begin(); itertrack != *iterevent.end(); itertrack++
482  //}
483 
484  cout << "Total reconstructed track count: " << totalrecotracks << endl;
485  cout << "Reached MC ID: " << globalmaxmcid << endl;
486  cout << "Reconstructable MC track count: " << goodmctracks << endl;
487  cout << "MC Tracks reconstructed: " << goodrecotracks << endl;
488  cout << "MC Tracks NOT reconstructed: " << norecotracks << endl;
489  cout << "Reco Status Good: " << recotrackstatusgoodcount << endl;
490  cout << "Reco Status Bad: " << recotrackstatusbadcount << endl;
491  cout << "Delta phi failed: " << deltaphifailed << endl;
492  cout << "Delta pt failed: " << deltaptfailed << endl;
493  cout << "Delta both failed: " << deltabothfailed << endl;
494  cout << "Secondary Tracks: " << secondaries << endl;
495 
496  cout << "MC Match Good 90: " << mcidgood90 << endl;
497  cout << "MC Match Good 80: " << mcidgood80 << endl;
498  cout << "MC Match Bad: " << mcidbad << endl;
499  cout << "MC Tracks with NULL pointers in Reco: " << recotracknullpointers << endl;
500  cout << "MC Tracks with NULL pointers: " << mctracknullpointers << endl;
501 
502  //c1->Print("hitdisplay.gif++");
503 
504  TFile outputstorage(histofilename, "UPDATE");
505  outputstorage.cd();
506  SaveAndUpdateHisto(hmcall, outputstorage);
507  SaveAndUpdateHisto(hmcprime, outputstorage);
508  SaveAndUpdateHisto(hmcrecoable, outputstorage);
509  SaveAndUpdateHisto(hdafuq, outputstorage);
510  SaveAndUpdateHisto(hmcrecoed, outputstorage);
511  SaveAndUpdateHisto(hmcnotrecoed, outputstorage);
512  SaveAndUpdateHisto(hphires, outputstorage);
513  SaveAndUpdateHisto(hptres, outputstorage);
514  SaveAndUpdateHisto(hphiresmc, outputstorage);
515  SaveAndUpdateHisto(hptresmc, outputstorage);
516  SaveAndUpdateHisto(hphiresmc90, outputstorage);
517  SaveAndUpdateHisto(hptresmc90, outputstorage);
518  SaveAndUpdateHisto(hphiresmc80, outputstorage);
519  SaveAndUpdateHisto(hptresmc80, outputstorage);
520  SaveAndUpdateHisto(hphiresmcbad, outputstorage);
521  SaveAndUpdateHisto(hptresmcbad, outputstorage);
522 
523  // ----- Finish -------------------------------------------------------
524  timer.Stop();
525  Double_t rtime = timer.RealTime();
526  Double_t ctime = timer.CpuTime();
527  cout << endl << endl;
528  cout << "Macro finished succesfully." << endl;
529  cout << "Real time " << rtime << " s, CPU time " << ctime << " s" << endl;
530  cout << endl;
531  // ------------------------------------------------------------------------
532  return 0;
533 }
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
r
double r
Definition: RiemannTest.C:14
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
CAMath::Abs
static T Abs(const T &x)
Definition: PndCAMath.h:39
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
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
SaveAndUpdateHisto
void SaveAndUpdateHisto(TH1 *currenthisto, TFile &storagefile)
Definition: runOnlineDisplayMCCheck.C:3
kDCH
Definition: production/tdrct/eta_c/track_check.C:4
void SaveAndUpdateHisto ( TH1 *  currenthisto,
TFile &  storagefile 
)

Definition at line 3 of file runOnlineDisplayMCCheckFaster.C.

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 }
Generated on Wed Apr 3 2019 08:44:23 for FairRoot/PandaRoot by   doxygen 1.8.5