FairRoot/PandaRoot
PndTrackingQualityBarrelTaskNewLinks.cxx
Go to the documentation of this file.
1 // -------------------------------------------------------------------------
2 // ----- PndTrackingQualityBarrelTaskNewLinks source file -----
3 // ----- Created 18/07/08 by T.Stockmanns -----
4 // -------------------------------------------------------------------------
5 
10 // C++ includes
11 #include <iostream>
12 #include <functional>
13 
14 // ROOT includes
15 #include "TROOT.h"
16 #include "TClonesArray.h"
17 #include "THStack.h"
18 
19 // FairRoot includes
20 #include "FairRootManager.h"
21 #include "FairRun.h"
22 #include "FairRuntimeDb.h"
23 #include "FairHit.h"
24 #include "FairMultiLinkedData.h"
25 
26 // PandaRoot includes
27 #include "PndMCTrack.h"
28 #include "PndTrack.h"
29 #include "PndMCEntry.h"
30 #include "PndSttHit.h"
31 #include "PndSttTube.h"
32 #include "PndSttMapCreator.h"
33 #include "RhoHistogram/RhoTuple.h"
34 
35 // Class includes
38 #include "PndSdsHit.h"
39 #include "TGeoManager.h"
40 #include "PndGeoHandling.h"
41 
42 
43 
44 // ----- Default constructor -------------------------------------------
46  FairTask("Creates PndMC Barrel test"),
47  fBranchNames(),fMCInfoBranchName("MCTrackInfo"), fRecoInfoBranchName("RecoTrackInfo"),
48  fMapLinkData(), fMapEfficiencies(), fMCInfoIdIdealId(),
49  fNGhosts(0),fTrack(NULL), fMCTrack(NULL), fTrackCand(NULL), fSttHitArray(NULL), fMCTrackInfo(NULL), fRecoTrackInfo(NULL),
50  fIdealTrack(NULL), fSttTubeArray(NULL),
51  ioman(NULL),
52  fMapTrackQualifikation(),
53  fTrackBranchName(trackBranchName), fIdealTrackBranchName(idealBranchName), fPndTrackOrTrackCand(pndTrackData),
54  fSttParameters(NULL), fTuple(NULL),
55  fPHisto(NULL),fPRelHisto(NULL),fPtHisto(NULL),fPtRelHisto(NULL),fQualyHisto(NULL),fQualyStack(NULL),
56  fQualyHisto_mc(NULL),fQualyHisto_neg(NULL),fQualyHisto_pos(NULL),fQualyHisto_all(NULL),
57  fEventNr(0) {
58 }
59 // -------------------------------------------------------------------------
60 
61 
62 // ----- Destructor ----------------------------------------------------
64 }
65 
66 // ----- Public method Init --------------------------------------------
69 
70  ioman = FairRootManager::Instance();
71  if (!ioman) {
72  std::cout << "-E- PndTrackingQualityBarrelTaskNewLinks::Init: "
73  << "RootManager not instantiated!" << std::endl;
74  return kFATAL;
75  }
76 
77  fTrack = (TClonesArray*) ioman->GetObject(fTrackBranchName);
78  fMCTrack = (TClonesArray*) ioman->GetObject("MCTrack");
79  // fIdealTrackCand = (TClonesArray*) ioman->GetObject(fIdealTrackBranchName);
80  fIdealTrack = (TClonesArray*) ioman->GetObject(fIdealTrackBranchName);
81  fSttHitArray = (TClonesArray*) ioman->GetObject("STTHit");
82 
83  // MC info for quality
84  fMCTrackInfo = new TClonesArray("PndTrackingQualityMCInfo");
85  ioman->Register(fMCInfoBranchName, "QualityAssurance", fMCTrackInfo, kTRUE); // CHECK
86  fRecoTrackInfo = new TClonesArray("PndTrackingQualityRecoInfo");
87  ioman->Register(fRecoInfoBranchName, "QualityAssurance", fRecoTrackInfo, kTRUE); // CHECK
88 
89  fTuple = new RhoTuple("qaTuple", "QA Rho");
90 
91  if (fBranchNames.size() == 0){
92  AddHitsBranchName("MVDHitsPixel");
93  AddHitsBranchName("MVDHitsStrip");
94  AddHitsBranchName("STTHit");
95  AddHitsBranchName("GEMHit");
96  // std::cout << "PndTrackingQualityAnalysis::Init() CorrectedSkewedHits present: " << FairRootManager::Instance()->GetBranchId("CorrectedSkewedHits") << " ";
97  // if (FairRootManager::Instance()->GetBranchId("CorrectedSkewedHits") > 0){
98  // std::cout << "kTRUE";
99  // AddHitsBranchName("CorrectedSkewedHits");
100  // }
101  // std::cout << std::endl;
102  }
103 
104 
105  for (unsigned int i = 0; i < fBranchNames.size(); i++){
106  fMapEfficiencies[fBranchNames[i]] = new TH2D(fBranchNames[i], fBranchNames[i], 100, 0., 100., 50, 0, 1.1);
107  fMapEfficiencies[fBranchNames[i]]->SetDrawOption("COLz");
108  }
109 
110  // ---------------------------------------- maps of STT tubes
112  fSttTubeArray = mapperStt->FillTubeArray();
113  // ---------------------------------------------------- end map
114 
115  std::cout
116  << "-I- PndTrackingQualityBarrelTaskNewLinks::Init: Initialization successfull"
117  << std::endl;
118 
119  return kSUCCESS;
120 }
121 
123  fPHisto = new TH1D("fPHisto", "Momentum Resolution", 1000, -1, 1);
124  fPHisto->GetXaxis()->SetTitle("p^{RECO} - p^{MC} / GeV");
125  fPHisto->GetYaxis()->SetTitle("counts");
126  fPRelHisto = new TH1D("fPRelHisto", "Relative Momentum Resolution", 1000, -1, 1);
127  fPRelHisto->GetXaxis()->SetTitle("(p^{RECO} - p^{MC}) / p^{MC}");
128  fPRelHisto->GetYaxis()->SetTitle("counts");
129  fPtHisto = new TH1D("fPtHisto", "Transverse Momentum Resolution", 1000, -1, 1);
130  fPtHisto->GetXaxis()->SetTitle("p_{t}^{RECO} - p_{t}^{MC} / GeV");
131  fPtHisto->GetYaxis()->SetTitle("counts");
132  fPtRelHisto = new TH1D("fPtRelHisto", "Relative Transverse Momentum Resolution", 1000, -1, 1);
133  fPtRelHisto->GetXaxis()->SetTitle("(p_{t}^{RECO} - p_{t}^{MC}) / p_{t}^{MC}");
134  fPtRelHisto->GetYaxis()->SetTitle("counts");
135 
136  fQualyHisto = new TH1I("fQualyHisto", "Quality of Trackfinding;;Counts", 26, -15.5, 10.5);
137  fQualyHisto->SetDrawOption("TEXT HIST");
139 
140  fQualyHisto_all = new TH1I("fQualyHisto_all", "Quality of Trackfinding;;Counts", 26, -15.5, 10.5);
141  fQualyHisto_pos = new TH1I("fQualyHisto_pos", "Quality of Trackfinding;;Counts", 26, -15.5, 10.5);
142  fQualyHisto_neg = new TH1I("fQualyHisto_neg", "Quality of Trackfinding;;Counts", 26, -15.5, 10.5);
143  fQualyHisto_mc = new TH1I("fQualyHisto_mc", "Quality of Trackfinding;;Counts", 26, -15.5, 10.5);
144  fQualyStack = new THStack();
145 }
146 
148  hist->GetXaxis()->SetBinLabel(hist->FindFixBin(qualityNumbers::kFullyFound), "Fully found");
149  hist->GetXaxis()->SetBinLabel(hist->FindFixBin(qualityNumbers::kPartiallyFound), "Partially found");
150  hist->GetXaxis()->SetBinLabel(hist->FindFixBin(qualityNumbers::kSpuriousFound), "Spurious found");
151  hist->GetXaxis()->SetBinLabel(hist->FindFixBin(qualityNumbers::kGhost), "Ghosts");
152  hist->GetXaxis()->SetBinLabel(hist->FindFixBin(qualityNumbers::kNotFound), "Total not found");
153  hist->GetXaxis()->SetBinLabel(hist->FindFixBin(qualityNumbers::kFound), "Total found");
154  hist->GetXaxis()->SetBinLabel(hist->FindFixBin(qualityNumbers::kPossibleSec), "Possible, Sec.");
155  hist->GetXaxis()->SetBinLabel(hist->FindFixBin(qualityNumbers::kPossiblePrim), "Possible, Prim.");
156  hist->GetXaxis()->SetBinLabel(hist->FindFixBin(qualityNumbers::kAtLeastThreeSec), ">= 3 Hits, Sec.");
157  hist->GetXaxis()->SetBinLabel(hist->FindFixBin(qualityNumbers::kAtLeastThreePrim), ">= 3 Hits, Prim.");
158  hist->GetXaxis()->SetBinLabel(hist->FindFixBin(qualityNumbers::kLessThanThreePrim), "< 3 Hits, Prim.");
159  hist->GetXaxis()->SetBinLabel(hist->FindFixBin(qualityNumbers::kMcPossibleSec), "MC: Possible, Sec.");
160  hist->GetXaxis()->SetBinLabel(hist->FindFixBin(qualityNumbers::kMcPossiblePrim), "MC: Possible, Prim.");
161  hist->GetXaxis()->SetBinLabel(hist->FindFixBin(qualityNumbers::kMcAtLeastThreeSec), "MC: >= 3 Hits, Sec.");
162  hist->GetXaxis()->SetBinLabel(hist->FindFixBin(qualityNumbers::kMcAtLeastThreePrim), "MC: >= 3 Hits, Prim.");
163  hist->GetXaxis()->SetBinLabel(hist->FindFixBin(qualityNumbers::kMcLessThanThreePrim), "MC: < 3 Hits, Prim.");
164 }
165 
166 
167 // -------------------------------------------------------------------------
169  FairRuntimeDb* rtdb = FairRun::Instance()->GetRuntimeDb();
170  fSttParameters = (PndGeoSttPar*) rtdb->getContainer("PndGeoSttPar");
171 }
172 
173 // ----- Public method Exec --------------------------------------------
175  fMCTrackInfo->Delete();
176  fRecoTrackInfo->Delete();
177 
178  std::cout << "----- Event " << fEventNr << " ------" << std::endl;
179 
180 
182  qaAna.SetVerbose(fVerbose);
184  qaAna.Init();
186 
187  std::map<Int_t, Int_t> qualiMap = qaAna.GetTrackQualification();
188  std::map<Int_t, Int_t> mcStatusMap = qaAna.GetTrackMCStatus();
189  std::map<Int_t, Int_t> mcFoundMap = qaAna.GetMCTrackFound();
190  std::map<Int_t, TVector3> recoPMap = qaAna.GetP();
191 // std::map<Int_t, Double_t> recoPtMap = qaAna.GetPt();
192 
193  FillQualyHisto(qualiMap, qaAna.GetNGhosts());
194  FillMCStatus(mcStatusMap);
196  MapToHist(qaAna.GetPResolution(), fPHisto);
200 
201  // fill MC Track Info ......................................
202  //for (std::map<Int_t, Int_t>::iterator iter = qualiMap.begin(); iter != qualiMap.end(); iter++) {
203  for (std::map<Int_t, Int_t>::iterator iter = mcStatusMap.begin(); iter != mcStatusMap.end(); iter++) {
204  Int_t mcTrackId = iter->first;
205  Int_t idealTrackId = qaAna.GetIdealTrackIdFromMCTrackId(mcTrackId);
206  Int_t trackQuality = iter->second;
207 
208  PndTrack *idealtrack = (PndTrack*) fIdealTrack->At(idealTrackId);
209 
210  PndMCTrack * myMcTrack = (PndMCTrack *) fMCTrack->At(mcTrackId);
211  Int_t pdgId = myMcTrack->GetPdgCode();
212 
213  int size = fMCTrackInfo->GetEntriesFast();
214  PndTrackingQualityMCInfo mctrackinfo = GetMCInfoFromIdealTrack(idealtrack);
215  mctrackinfo.SetMCTrackID(mcTrackId);
216  mctrackinfo.SetQuality(trackQuality);
217  mctrackinfo.SetPDGCode(pdgId);
218  // mctrackinfo.SetReconstructabilityStatus();
219 
220  if(mctrackinfo.GetNofMCPoints() > 0) {
221  new((*fMCTrackInfo)[size]) PndTrackingQualityMCInfo(mctrackinfo);
222  fMCInfoIdIdealId[idealTrackId] = size;
223  }
224  // cout << "MCTRack " << mctrackinfo.GetMCTrackID() << endl;
225  }
226  // ............................................................
227 
228  // loop over reco track info and associate the mc track info
229  for(int itrk = 0; itrk < fRecoTrackInfo->GetEntriesFast(); itrk++) {
231  Int_t idealTrackId = qaAna.GetIdealTrackIdFromRecoTrackId(recoinfo->GetRecoTrackID());
233  recoinfo->SetMCTrackInfo(mctrackinfo);
234  }
235 
236  // associate reconstructed and mc tracks
238 
239 
240  // Save the Tree (/RhoTuple) for some possible additional analysis
241  for (std::map<Int_t, Int_t>::iterator iter = qualiMap.begin(); iter != qualiMap.end(); iter++) {
242  Int_t mcTrackId = iter->first;
243  Int_t trackQuality = iter->second;
244 
245  TVector3 recoMomentum = recoPMap[mcTrackId];
246 
247  PndMCTrack * myMcTrack = (PndMCTrack *) fMCTrack->At(mcTrackId);
248  TVector3 mcMomentum = myMcTrack->GetMomentum();
249  Int_t pdgId = myMcTrack->GetPdgCode();
250 
251  fTuple->Column("EvtNr", (Int_t) fEventNr); // number of the currently processed event
252 
253  fTuple->Column("McTrackId", (Int_t) mcTrackId); // id of the currently processed track
254  fTuple->Column("McTrackFoundNTimes", (Int_t) mcFoundMap[mcTrackId]); // A given MC track was found N times
255  fTuple->Column("TrackQuality", (Int_t) trackQuality); // the quality of the current track
256  fTuple->Column("Mc_TrackQuality", (Int_t) mcStatusMap[mcTrackId] - 6); // the mc quality of the current track, -6 to match the global qualityNumbers
257 
258  fTuple->Column("PdgId", (Int_t) pdgId);
259  fTuple->Column("Mc_px", (Double_t) mcMomentum.Px()); // To compare to the additionally saved pt histogram, call Draw("Reco_pt:Mc_pt","TrackQuality > 0") on the Tuple
260  fTuple->Column("Mc_py", (Double_t) mcMomentum.Py());
261  fTuple->Column("Mc_pz", (Double_t) mcMomentum.Pz());
262  fTuple->Column("Mc_pt", (Double_t) mcMomentum.Pt());
263  fTuple->Column("Reco_px", (Double_t) recoMomentum.Px());
264  fTuple->Column("Reco_py", (Double_t) recoMomentum.Py());
265  fTuple->Column("Reco_pz", (Double_t) recoMomentum.Pz());
266  fTuple->Column("Reco_pt", (Double_t) recoMomentum.Pt());
267 // fTuple->Column("Reco_pt2", (Double_t) recoPtMap[mcTrackId]); // cross check
268 
269  fTuple->DumpData();
270  }
271 
272 // if (fVerbose > 1)
273  qaAna.PrintTrackQualityMap();
274  fEventNr++;
275 }
276 
278 {
279  Int_t result = 0;
280  for (unsigned int branchIndex = 0; branchIndex < fBranchNames.size(); branchIndex++){
281  result += trackData->GetLinksWithType(ioman->GetBranchId(fBranchNames[branchIndex])).GetNLinks();
282  }
283  return result;
284 }
285 
286 
287 void PndTrackingQualityBarrelTaskNewLinks::FillQualyHisto(std::map<Int_t, Int_t> trackQualifikation, Int_t nGhosts)
288 {
289 
290  fQualyHisto->Fill(qualityNumbers::kGhost, nGhosts);
291  for(std::map<Int_t, Int_t>::iterator iter = trackQualifikation.begin(); iter != trackQualifikation.end(); iter++){
292  fQualyHisto->Fill(iter->second);
293  if (iter->second > 0){
295  }
296  else {
298  }
299  }
300 }
301 
302 void PndTrackingQualityBarrelTaskNewLinks::FillMCStatus(std::map<Int_t, Int_t> trackMCStatus)
303 {
305  for(std::map<Int_t, Int_t>::iterator iter = trackMCStatus.begin(); iter != trackMCStatus.end(); iter++){
306  fQualyHisto->Fill(iter->second - mcOffset); // mcOffset = 6
307  }
308 }
309 
310 
311 void PndTrackingQualityBarrelTaskNewLinks::FillEfficiencies(std::map<Int_t, std::map<TString, std::pair<Double_t, Int_t > > > efficiencies)
312 {
313  for (std::map<Int_t, std::map<TString, std::pair<Double_t, Int_t> > >::iterator iterTracks = efficiencies.begin(); iterTracks != efficiencies.end(); iterTracks++ ){
314  std::map<TString, std::pair<Double_t, Int_t> > branchEfficiency = iterTracks->second;
315  for (std::map<TString, std::pair<Double_t, Int_t> >::iterator iterBranch = branchEfficiency.begin(); iterBranch != branchEfficiency.end(); iterBranch++ ){
316  fMapEfficiencies[iterBranch->first]->Fill(iterBranch->second.second, iterBranch->second.first);
317  }
318  }
319 }
320 void PndTrackingQualityBarrelTaskNewLinks::MapToHist (std::map<Int_t, Double_t> map, TH1 * histo) {
321  for (std::map<Int_t, Double_t>::iterator iter = map.begin(); iter != map.end(); iter++ ){
322  histo->Fill(iter->second);
323  }
324 }
325 
327  ColorHistogram();
332  fQualyStack->SetName("fQualyHistoColor");
333  fQualyStack->SetTitle(fQualyHisto->GetTitle());
334  // gROOT->SetBatch(kTRUE);
335  // fQualyHisto->Draw();
336  // LabelQualyHistogram((TH1*) fQualyStack);
337 
338  for (unsigned int i = 0; i < fBranchNames.size(); i++){
339  fMapEfficiencies[fBranchNames[i]]->Write();
340  }
341  fPHisto->Write();
342  fPRelHisto->Write();
343  fPtHisto->Write();
344  fPtRelHisto->Write();
345  fQualyHisto->Write();
346  fQualyStack->Write();
347  fTuple->GetInternalTree()->Write();
348 
349  Int_t allTracksWithHits = 0;
350  Int_t allPossibleTracksWithHits = 0;
351  Int_t allTracksWithHitsNotFound = 0;
352 
353  allTracksWithHits += fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kMcAtLeastThreePrim)); // Todo
354  allTracksWithHits += fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kMcAtLeastThreeSec));
355  allTracksWithHits += fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kMcPossiblePrim));
356  allTracksWithHits += fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kMcPossibleSec));
357 
358  allPossibleTracksWithHits += fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kMcPossiblePrim));
359  allPossibleTracksWithHits += fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kMcPossibleSec));
360 
361 
362 
363  allTracksWithHitsNotFound += fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kAtLeastThreePrim));
364  allTracksWithHitsNotFound += fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kAtLeastThreeSec));
365  allTracksWithHitsNotFound += fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kPossiblePrim));
366  allTracksWithHitsNotFound += fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kPossibleSec));
367 
368  std::cout << "fQualyHisto: All Tracks: " << fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(-11)) + allTracksWithHits << std::endl
369  << " Primary Tracks < 3 hits: " << fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(-11)) << std::endl
370  << " All Tracks with hits: " << allTracksWithHits << ". Not Found: " << allTracksWithHitsNotFound << std::endl
371  << " Primary Tracks with >= 3 hits, but not a possible track: " << fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kMcAtLeastThreePrim)) << ". Not Found: " << fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kAtLeastThreePrim)) << std::endl
372  << " Secondary Tracks with >= 3 hits, but not a possible track: " << fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kMcAtLeastThreeSec)) << ". Not Found: " << fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kAtLeastThreeSec)) << std::endl
373  << " Primary Tracks possible: " << fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kMcPossiblePrim)) << ". Not Found: " << fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kPossiblePrim)) << std::endl
374  << " Secondary Tracks possible: " << fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kMcPossibleSec)) << ". Not Found: " << fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kPossibleSec)) << std::endl
375 
376  << " FullyFound: " << fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kFullyFound)) << " "
377  << (Double_t)fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kFullyFound)) / allTracksWithHits * 100 << "% (of all tracks), "
378  << (Double_t)fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kFullyFound)) / allPossibleTracksWithHits * 100 << "% (of all tracks possible)"
379 
380  << " PartlyFound: " << fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kPartiallyFound)) << " "
381  << (Double_t)fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kPartiallyFound)) / allTracksWithHits * 100 << "% "
382  << (Double_t)fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kPartiallyFound)) / allPossibleTracksWithHits * 100 << "% "
383 
384  << " Spurious: " << fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kSpuriousFound)) << " "
385  << (Double_t)fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kSpuriousFound)) / allTracksWithHits * 100 << "% "
386  << (Double_t)fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kSpuriousFound)) / allPossibleTracksWithHits * 100 << "% "
387 
388  << " Ghosts: " << fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kGhost)) << " "
389  << (Double_t)fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kGhost)) / allTracksWithHits * 100 << "% "
390  << (Double_t)fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kGhost)) / allPossibleTracksWithHits * 100 << "% " << std::endl;
391 }
392 
394  fQualyHisto_pos->SetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kGhost), fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kGhost)));
395  fQualyHisto_pos->SetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kSpuriousFound), fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kSpuriousFound)));
397  fQualyHisto_pos->SetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kFullyFound), fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kFullyFound)));
398  fQualyHisto_pos->SetLineColor(kGreen + 2);
399 
400  fQualyHisto_all->SetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kFound), fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kFound)));
401  fQualyHisto_all->SetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kNotFound), fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kNotFound)));
402  fQualyHisto_all->SetLineColor(kBlack);
403 
404 
405  fQualyHisto_neg->SetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kPossibleSec), fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kPossibleSec)));
406  fQualyHisto_neg->SetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kPossiblePrim), fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kPossiblePrim)));
410  fQualyHisto_neg->SetLineColor(kRed + 2);
411 
412  fQualyHisto_mc->SetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kMcPossibleSec), fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kMcPossibleSec)));
417  fQualyHisto_mc->SetLineColor(kBlue);
418 
419 }
420 
421 
422 //SG!!!
423 
425 {
426  if(!hit) return 0;
427  Int_t sensorID = hit->GetSensorID();
428  gGeoManager->cd(PndGeoHandling::Instance()->GetPath(sensorID));
429  TGeoHMatrix* transMat = gGeoManager->GetCurrentMatrix();
430  if( !transMat ){
431  std::cout<<"\n\nIsBarrelMVD():: NO transition Matrix for MVD sensor "<<sensorID<<"\n\n!!!"<<std::endl;
432  return 0;
433  }
434  Double_t *mmm = transMat->GetRotationMatrix();
435  if( !mmm ){
436  std::cout<<"\n\nIsBarrelMVD():: Can not extract transition Matrix for MVD sensor "<<sensorID<<", something is completely wrong\n\n!!!"<<std::endl;
437  exit(0);//SG!!
438  return 0;
439  }
440  bool forward = ( fabs(mmm[6]) < 0.999 && fabs(mmm[7]) < 0.999 );
441  return ( !forward );
442 }
443 
444 
445 
447 
448  PndTrackCand *idealtrkcand = idealtrack->GetTrackCandPtr();
449 
450  //Int_t nofsttpoint = idealtrkcand->GetNHitsDet(FairRootManager::Instance()->GetBranchId("STTHit"));
451  Int_t nofmvdpixpoint = idealtrkcand->GetNHitsDet(FairRootManager::Instance()->GetBranchId("MVDHitsPixel"));
452  Int_t nofmvdstrpoint = idealtrkcand->GetNHitsDet(FairRootManager::Instance()->GetBranchId("MVDHitsStrip"));
453  //Int_t nofmvdpoint = nofmvdpixpoint + nofmvdstrpoint;
454  Int_t nofgempoint = idealtrkcand->GetNHitsDet(FairRootManager::Instance()->GetBranchId("GEMHit"));
455 
456  int nofsttskewpoint = 0, nofsttparalpoint = 0;
457  // this loop counts skewed (--> parallel) STT/FTS hits
458  for(UInt_t ihit = 0; ihit < idealtrkcand->GetNHits(); ihit++) {
459  PndTrackCandHit idealcandhit = idealtrkcand->GetSortedHit(ihit);
460  Int_t hitID = idealcandhit.GetHitId();
461  Int_t detID = idealcandhit.GetDetId();
462 
463  if(detID != FairRootManager::Instance()->GetBranchId("STTHit")) continue;
464  PndSttHit *stthit = (PndSttHit*) fSttHitArray->At(hitID);
465  Int_t tubeID = stthit->GetTubeID();
466  PndSttTube *tube = (PndSttTube*) fSttTubeArray->At(tubeID);
467  if(tube->IsSkew()) nofsttskewpoint++;
468  else nofsttparalpoint++;
469  }
470 
471  //SG!!!
472  /*
473  nofmvdpixpoint = 0;
474  nofmvdstrpoint = 0;
475  nofgempoint = 0;
476 
477  TClonesArray *mvdPixArray = (TClonesArray*) FairRootManager::Instance()->GetObject("MVDHitsPixel");
478  TClonesArray *mvdStripArray = (TClonesArray*) FairRootManager::Instance()->GetObject("MVDHitsStrip");
479 
480  for(Int_t ihit = 0; ihit < idealtrkcand->GetNHits(); ihit++) {
481  PndTrackCandHit idealcandhit = idealtrkcand->GetSortedHit(ihit);
482  Int_t hitID = idealcandhit.GetHitId();
483  Int_t detID = idealcandhit.GetDetId();
484 
485  if(detID == FairRootManager::Instance()->GetBranchId("MVDHitsPixel")){
486  const PndSdsHit *hit = dynamic_cast<const PndSdsHit*>(mvdPixArray->At(hitID));
487  if(!hit) {
488  std::cout << "GetMCInfoFromIdealTrack() : MVD pix. hit is not found "<<std::endl;
489  exit(0);
490  continue;
491  }
492  nofmvdpixpoint+=IsBarrelMVD(hit);
493  }
494  else if(detID == FairRootManager::Instance()->GetBranchId("MVDHitsStrip")){
495  const PndSdsHit *hit = dynamic_cast<const PndSdsHit*>(mvdStripArray->At(hitID));
496  if(!hit) {
497  std::cout << "GetMCInfoFromIdealTrack() : MVD strip hit is not found "<<std::endl;
498  exit(0);
499  continue;
500  }
501  nofmvdstrpoint+=IsBarrelMVD(hit);
502  }
503  }
504 
505  nofmvdpoint = nofmvdpixpoint + nofmvdstrpoint;
506  */
507  PndTrackingQualityMCInfo info(nofmvdpixpoint, nofmvdstrpoint, nofsttparalpoint, nofsttskewpoint, nofgempoint);
508 
509  // CHECK
510  // Bool_t isreco = Reconstructability(nofmvdpixpoint, nofmvdstrpoint, nofsttparalpoint, nofsttskewpoint, nofgempoint, nofscitilpoint);
511  // info.SetReconstructability(isreco);
512 
513 
514  info.SetPositionFirst(idealtrack->GetParamFirst().GetPosition());
515  info.SetMomentumFirst(idealtrack->GetParamFirst().GetMomentum());
516  info.SetPositionLast(idealtrack->GetParamLast().GetPosition());
517  info.SetMomentumLast(idealtrack->GetParamLast().GetMomentum());
518  info.SetCharge(idealtrack->GetParamFirst().GetQ());
519 
520 
521  return info;
522 
523 }
696  // loop over mc track infos
697  for(int imctrk = 0; imctrk < fMCTrackInfo->GetEntriesFast(); imctrk++) {
699  int mctrackid0 = mcinfo->GetMCTrackID();
700  if(mcinfo->GetAssoRecoTrackID() != -1) continue;
701  double tmpeff = 0., tmppur = 0.;
702  int tmptruerecotrackid = -1;
703 
704  // loop over reco track infos
705  PndTrackingQualityRecoInfo *tmprecoinfo = NULL;
706  for(int itrk = 0; itrk < fRecoTrackInfo->GetEntriesFast(); itrk++) {
708  int mctrackid = recoinfo->GetMCTrackID();
709  if(mctrackid != mctrackid0) continue;
710  mcinfo->SetRecoTrackID(recoinfo->GetRecoTrackID());
711  recoinfo->SetClone();
712  // it must have either the higher efficiency ...
713  if(recoinfo->GetEfficiency() < tmpeff) continue;
714  // ... or, if they are even, the highest purity
715  if(recoinfo->GetEfficiency() == tmpeff && recoinfo->GetPurity() < tmppur) continue;
716 
717  tmpeff = recoinfo->GetEfficiency();
718  tmppur = recoinfo->GetPurity();
719  tmptruerecotrackid = recoinfo->GetRecoTrackID();
720  tmprecoinfo = recoinfo;
721  }
722  if(tmprecoinfo == NULL) continue;
723 
724  tmprecoinfo->SetTrue();
725  mcinfo->SetAssoRecoTrackID(tmptruerecotrackid);
726 
727  }
728 }
729 
730 
731 
732 
UInt_t GetNHitsDet(Int_t detId)
int fVerbose
Definition: poormantracks.C:24
static const int kMcAtLeastThreeSec
Definition: PndTrackingQA.h:52
Int_t i
Definition: run_full.C:25
static const int kGhost
Definition: PndTrackingQA.h:63
void SetMCTrackID(Int_t mctrackid)
PndTransMap * map
Definition: sim_emc_apd.C:99
exit(0)
Int_t GetPdgCode() const
Definition: PndMCTrack.h:73
static const int kLessThanThreePrim
Definition: PndTrackingQA.h:47
PndTrackCandHit GetSortedHit(UInt_t i)
Definition: PndTrackCand.h:54
static const int kPossiblePrim
Definition: PndTrackingQA.h:44
static const int kAtLeastThreeSec
Definition: PndTrackingQA.h:45
TVector3 GetMomentum() const
Definition: PndMCTrack.h:78
void SetRecoTrackID(int recotrkid)
static const int kPartiallyFound
Definition: PndTrackingQA.h:60
static const int kMcPossibleSec
Definition: PndTrackingQA.h:50
void SetMCTrackInfo(PndTrackingQualityMCInfo *info)
TGeoManager * gGeoManager
static const int kPossibleSec
Definition: PndTrackingQA.h:43
Double_t
static const int kAtLeastThreePrim
Definition: PndTrackingQA.h:46
FairTrackParP GetParamLast()
Definition: PndTrack.h:50
FairRuntimeDb * rtdb
Definition: hit_dirc.C:66
Int_t GetTubeID() const
Definition: PndSttHit.h:75
TClonesArray * FillTubeArray()
void Column(const char *label, Bool_t value, Bool_t defval=0, const char *block=0)
Definition: RhoTuple.cxx:56
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
static PndGeoHandling * Instance()
static const int kMcLessThanThreePrim
Definition: PndTrackingQA.h:54
void DumpData()
Definition: RhoTuple.cxx:391
static const int kFullyFound
Definition: PndTrackingQA.h:61
static const int kFound
Definition: PndTrackingQA.h:66
UInt_t GetNHits() const
Definition: PndTrackCand.h:59
static const int kMcAtLeastThreePrim
Definition: PndTrackingQA.h:53
ClassImp(PndAnaContFact)
TTree * GetInternalTree()
Definition: RhoTuple.h:207
PndSdsMCPoint * hit
Definition: anasim.C:70
Int_t GetSensorID() const
Definition: PndSdsHit.h:90
PndTrackCand * GetTrackCandPtr()
Definition: PndTrack.h:48
Int_t GetHitId() const
TH1F * hist
Int_t GetDetId() const
static const int kMcPossiblePrim
Definition: PndTrackingQA.h:51
FairTrackParP GetParamFirst()
Definition: PndTrack.h:49
static const int kSpuriousFound
Definition: PndTrackingQA.h:59
bool IsSkew()
Definition: PndSttTube.h:70
static const int kNotFound
Definition: PndTrackingQA.h:65