FairRoot/PandaRoot
PndTrackingQATask.cxx
Go to the documentation of this file.
1 // -------------------------------------------------------------------------
2 // ----- PndTrackingQATask source file -----
3 // ----- Created 18/07/08 by T.Stockmanns -----
4 // -------------------------------------------------------------------------
5 
10 // C++ includes
11 #include <PndTrackingQA.h>
12 #include <PndTrackingQATask.h>
13 #include <iostream>
14 #include <functional>
15 
16 // ROOT includes
17 #include "TROOT.h"
18 #include "TClonesArray.h"
19 #include "THStack.h"
20 
21 // FairRoot includes
22 #include "FairRootManager.h"
23 #include "FairRun.h"
24 #include "FairRuntimeDb.h"
25 #include "FairHit.h"
26 #include "FairMultiLinkedData.h"
27 
28 // PandaRoot includes
29 #include "PndMCTrack.h"
30 #include "PndTrack.h"
31 #include "PndSttHit.h"
32 #include "PndSttTube.h"
33 #include "PndSttMapCreator.h"
34 #include "RhoHistogram/RhoTuple.h"
35 
36 // Class includes
37 
38 // ----- Default constructor -------------------------------------------
39 PndTrackingQATask::PndTrackingQATask(TString trackBranchName, TString idealBranchName, Bool_t pndTrackData) :
40  FairTask("Creates PndMC test"), fMCInfoBranchName("MCTrackInfo"), fRecoInfoBranchName("RecoTrackInfo"), fTrackBranchName(trackBranchName), fIdealTrackBranchName(idealBranchName), fPndTrackOrTrackCand(pndTrackData), fEventNr(0) {
41 }
42 // -------------------------------------------------------------------------
43 
44 
45 // ----- Destructor ----------------------------------------------------
47 }
48 
49 // ----- Public method Init --------------------------------------------
52 
53  ioman = FairRootManager::Instance();
54  if (!ioman) {
55  std::cout << "-E- PndTrackingQATask::Init: "
56  << "RootManager not instantiated!" << std::endl;
57  return kFATAL;
58  }
59 
60  fTrack = (TClonesArray*) ioman->GetObject(fTrackBranchName);
61  fMCTrack = (TClonesArray*) ioman->GetObject("MCTrack");
62  // fIdealTrackCand = (TClonesArray*) ioman->GetObject(fIdealTrackBranchName);
63  fIdealTrack = (TClonesArray*) ioman->GetObject(fIdealTrackBranchName);
64  fSttHitArray = (TClonesArray*) ioman->GetObject("STTHit");
65 
66  if (fTrack == nullptr){
67  std::cout << "-E- PndTrackingQATask::Init " << "no track branch " << fTrackBranchName << std::endl;
68  return kFATAL;
69  }
70  if (fMCTrack == nullptr){
71  std::cout << "-E- PndTrackingQATask::Init " << "no MC track branch " << std::endl;
72  return kFATAL;
73  }
74 
75  if (fIdealTrack == nullptr){
76  std::cout << "-E- PndTrackingQATask::Init " << "no ideal track branch " << fIdealTrackBranchName << std::endl;
77  return kFATAL;
78  }
79 
80  // MC info for quality
81  fMCTrackInfo = new TClonesArray("PndTrackingQualityMCInfo");
82  ioman->Register(fMCInfoBranchName, "QualityAssurance", fMCTrackInfo, kTRUE); // CHECK
83  fRecoTrackInfo = new TClonesArray("PndTrackingQualityRecoInfo");
84  ioman->Register(fRecoInfoBranchName, "QualityAssurance", fRecoTrackInfo, kTRUE); // CHECK
85 
86  fTuple = new RhoTuple("qaTuple", "QA Rho");
87 
88  if (fBranchNames.size() == 0){
89  AddHitsBranchName("MVDHitsPixel");
90  AddHitsBranchName("MVDHitsStrip");
91  AddHitsBranchName("STTHit");
92  AddHitsBranchName("GEMHit");
93  AddHitsBranchName("FTSHit");
94  // std::cout << "PndTrackingQualityAnalysis::Init() CorrectedSkewedHits present: " << FairRootManager::Instance()->GetBranchId("CorrectedSkewedHits") << " ";
95  // if (FairRootManager::Instance()->GetBranchId("CorrectedSkewedHits") > 0){
96  // std::cout << "kTRUE";
97  // AddHitsBranchName("CorrectedSkewedHits");
98  // }
99  // std::cout << std::endl;
100  }
101 
102  if (fPossibleTrackFunctorName.Length() == 0)
103  fPossibleTrackFunctorName = "StandardTrackFunctor";
104 
105  SetFunctor();
106  for (size_t i = 0; i < fBranchNames.size(); i++){
107  fMapEfficiencies[fBranchNames[i]] = new TH2D(fBranchNames[i], fBranchNames[i], 100, 0., 100., 50, 0, 1.1);
108  fMapEfficiencies[fBranchNames[i]]->SetDrawOption("COLz");
109  }
110 
111  // ---------------------------------------- maps of STT tubes
113  fSttTubeArray = mapperStt->FillTubeArray();
114  // ---------------------------------------------------- end map
115 
116 // std::cout << "-I- PndTrackingQATask::Init: Initialization successfull" << std::endl;
117 
118  return kSUCCESS;
119 }
120 
122  fIdealTracksPerEvent = new TH1I("fIdealTracksPerEvent", "Ideal Tracks per Event", 1000, -0.5,999.5);
123  fIdealTracksPerEvent->GetXaxis()->SetTitle("Tracks/Event");
124  fIdealPHisto = new TH1I("fIdealPHisto", "Ideal total Momentum", 1500, -0.5,14.5);
125  fIdealPHisto->GetXaxis()->SetTitle("p [GeV/c]");
126  fIdealPtHisto = new TH1I("fIdealPtHisto", "Ideal Tansversal Momentum", 1500, -0.5,14.5);
127  fIdealPtHisto->GetXaxis()->SetTitle("p_{t} [GeV/c]");
128  fIdealPlHisto = new TH1I("fIdealPlHisto", "Ideal Longitudinal Momentum", 1500, -0.5,14.5);
129  fIdealPlHisto->GetXaxis()->SetTitle("p_{l} [GeV/c]");
130  fPHisto = new TH1D("fPHisto", "Momentum Resolution", 1000, -1, 1);
131  fPHisto->GetXaxis()->SetTitle("p^{RECO} - p^{MC} / GeV");
132  fPHisto->GetYaxis()->SetTitle("counts");
133  fPRelHisto = new TH1D("fPRelHisto", "Relative Momentum Resolution", 1000, -1, 1);
134  fPRelHisto->GetXaxis()->SetTitle("(p^{RECO} - p^{MC}) / p^{MC}");
135  fPRelHisto->GetYaxis()->SetTitle("counts");
136  fPtHisto = new TH1D("fPtHisto", "Transverse Momentum Resolution", 1000, -1, 1);
137  fPtHisto->GetXaxis()->SetTitle("p_{t}^{RECO} - p_{t}^{MC} / GeV");
138  fPtHisto->GetYaxis()->SetTitle("counts");
139  fPtRelHisto = new TH1D("fPtRelHisto", "Relative Transverse Momentum Resolution", 1000, -1, 1);
140  fPtRelHisto->GetXaxis()->SetTitle("(p_{t}^{RECO} - p_{t}^{MC}) / p_{t}^{MC}");
141  fPtRelHisto->GetYaxis()->SetTitle("counts");
142  fPlHisto = new TH1D("fPlHisto", "Longitudinal Momentum Resolution", 1000, -1, 1);
143  fPlHisto->GetXaxis()->SetTitle("p_{l}^{RECO} - p_{l}^{MC} / GeV");
144  fPlHisto->GetYaxis()->SetTitle("counts");
145  fPlRelHisto = new TH1D("fPlRelHisto", "Relative Longitudinal Momentum Resolution", 1000, -1, 1);
146  fPlRelHisto->GetXaxis()->SetTitle("(p_{l}^{RECO} - p_{l}^{MC}) / p_{l}^{MC}");
147  fPlRelHisto->GetYaxis()->SetTitle("counts");
148 
149  fQualyHisto = new TH1I("fQualyHisto", "Quality of Trackfinding;;Counts", 26, -15.5, 10.5);
150  fQualyHisto->SetDrawOption("TEXT HIST");
152 
153  fQualyHisto_all = new TH1I("fQualyHisto_all", "Quality of Trackfinding;;Counts", 26, -15.5, 10.5);
154  fQualyHisto_pos = new TH1I("fQualyHisto_pos", "Quality of Trackfinding;;Counts", 26, -15.5, 10.5);
155  fQualyHisto_neg = new TH1I("fQualyHisto_neg", "Quality of Trackfinding;;Counts", 26, -15.5, 10.5);
156  fQualyHisto_mc = new TH1I("fQualyHisto_mc", "Quality of Trackfinding;;Counts", 26, -15.5, 10.5);
157  fQualyStack = new THStack();
158 
159  fQualyHisto_rel_all = new TH1D("fQualyHisto_rel_all", "Quality of Trackfinding;;Relative", 26, -15.5, 10.5);
160  fQualyHisto_rel_all->SetBarWidth(0.45);
161  fQualyHisto_rel_all->SetBarOffset(0.1);
162  fQualyHisto_rel_all->SetFillColor(kBlue);
163 
164  fQualyHisto_rel_possible = new TH1D("fQualyHisto_rel_possible", "Quality of Trackfinding;;Relative", 26, -15.5, 10.5);
165  fQualyHisto_rel_possible->SetBarWidth(0.4);
166  fQualyHisto_rel_possible->SetBarOffset(0.55);
167  fQualyHisto_rel_possible->SetFillColor(kRed);
168 
169 
172 }
173 
175  hist->GetXaxis()->SetBinLabel(hist->FindFixBin(qualityNumbers::kFullyFound), "Fully found");
176  hist->GetXaxis()->SetBinLabel(hist->FindFixBin(qualityNumbers::kPartiallyFound), "Partially found");
177  hist->GetXaxis()->SetBinLabel(hist->FindFixBin(qualityNumbers::kSpuriousFound), "Spurious found");
178  hist->GetXaxis()->SetBinLabel(hist->FindFixBin(qualityNumbers::kGhost), "Ghosts");
179  hist->GetXaxis()->SetBinLabel(hist->FindFixBin(qualityNumbers::kNotFound), "Total not found");
180  hist->GetXaxis()->SetBinLabel(hist->FindFixBin(qualityNumbers::kFound), "Total found");
181  hist->GetXaxis()->SetBinLabel(hist->FindFixBin(qualityNumbers::kPossibleSec), "Possible, Sec.");
182  hist->GetXaxis()->SetBinLabel(hist->FindFixBin(qualityNumbers::kPossiblePrim), "Possible, Prim.");
183  hist->GetXaxis()->SetBinLabel(hist->FindFixBin(qualityNumbers::kAtLeastThreeSec), ">= 3 Hits, Sec.");
184  hist->GetXaxis()->SetBinLabel(hist->FindFixBin(qualityNumbers::kAtLeastThreePrim), ">= 3 Hits, Prim.");
185  hist->GetXaxis()->SetBinLabel(hist->FindFixBin(qualityNumbers::kLessThanThreePrim), "< 3 Hits, Prim.");
186  hist->GetXaxis()->SetBinLabel(hist->FindFixBin(qualityNumbers::kMcPossibleSec), "MC: Possible, Sec.");
187  hist->GetXaxis()->SetBinLabel(hist->FindFixBin(qualityNumbers::kMcPossiblePrim), "MC: Possible, Prim.");
188  hist->GetXaxis()->SetBinLabel(hist->FindFixBin(qualityNumbers::kMcAtLeastThreeSec), "MC: >= 3 Hits, Sec.");
189  hist->GetXaxis()->SetBinLabel(hist->FindFixBin(qualityNumbers::kMcAtLeastThreePrim), "MC: >= 3 Hits, Prim.");
190  hist->GetXaxis()->SetBinLabel(hist->FindFixBin(qualityNumbers::kMcLessThanThreePrim), "MC: < 3 Hits, Prim.");
191 }
192 
193 
194 // -------------------------------------------------------------------------
196  FairRuntimeDb* rtdb = FairRun::Instance()->GetRuntimeDb();
197  fSttParameters = (PndGeoSttPar*) rtdb->getContainer("PndGeoSttPar");
198 }
199 
201 {
203 }
204 
205 // ----- Public method Exec --------------------------------------------
206 void PndTrackingQATask::Exec(Option_t*) {
207  fMCTrackInfo->Delete();
208  fRecoTrackInfo->Delete();
209 
210  if (fVerbose > 0)
211  std::cout << "----- Event " << fEventNr << " ------" << std::endl;
212 
213 
215  qaAna.SetVerbose(fVerbose);
217  qaAna.Init();
219 
220  std::map<Int_t, Int_t> qualiMap = qaAna.GetTrackQualification();
221  std::map<Int_t, Int_t> mcStatusMap = qaAna.GetTrackMCStatus();
222  std::map<Int_t, Int_t> mcFoundMap = qaAna.GetMCTrackFound();
223  std::map<Int_t, TVector3> recoPMap = qaAna.GetP();
224 // std::map<Int_t, Double_t> recoPtMap = qaAna.GetPt();
225 
226 
227  FillQualyHisto(qualiMap, qaAna.GetNGhosts());
228  FillMCStatus(mcStatusMap);
229  fIdealTracksPerEvent->Fill(fIdealTrack->GetEntries());
230  for (int i = 0; i < fIdealTrack->GetEntries(); i++){
231  PndTrack* myTrack = (PndTrack*)fIdealTrack->At(i);
232  fIdealPHisto->Fill(myTrack->GetParamFirst().GetMomentum().Mag());
233  fIdealPtHisto->Fill(myTrack->GetParamFirst().GetMomentum().Pt());
234  fIdealPlHisto->Fill(myTrack->GetParamFirst().GetMomentum().Pz());
235  }
236 
238  MapToHist(qaAna.GetPResolution(), fPHisto);
244 
245  // fill MC Track Info ......................................
246  for (std::map<Int_t, Int_t>::iterator iter = qualiMap.begin(); iter != qualiMap.end(); iter++)
247  {
248  Int_t mcTrackId = iter->first;
249 
250  Int_t idealTrackId = qaAna.GetIdealTrackIdFromMCTrackId(mcTrackId);
251  if (idealTrackId < 0){
252  std::cout << "-W- PndTrackingQATask::Exec no idealTrack for mcTrack " << mcTrackId << std::endl;
253  continue;
254  }
255  Int_t trackQuality = iter->second;
256 
257  PndTrack *idealtrack = (PndTrack*) fIdealTrack->At(idealTrackId);
258 
259  if (idealtrack == 0){
260  std::cout << "-E- No ideal track found for idealTrackId " << idealTrackId << std::endl;
261  continue;
262  }
263  if (mcTrackId == -1){
264  std::cout << "-W- PndTrackingQATask::Exec mcTrackId == -1" << std::endl;
265  continue;
266  }
267 
268  PndMCTrack * myMcTrack = (PndMCTrack *) fMCTrack->At(mcTrackId);
269 
270  if (myMcTrack == 0){
271  std::cout << "-E- PndTrackingQATask::Exec mcMyTrack == 0" << std::endl;
272  continue;
273  }
274  Int_t pdgId = myMcTrack->GetPdgCode();
275 
276  int size = fMCTrackInfo->GetEntriesFast();
277  PndTrackingQualityMCInfo mctrackinfo = GetMCInfoFromIdealTrack(idealtrack);
278  mctrackinfo.SetMCTrackID(mcTrackId);
279  mctrackinfo.SetQuality(trackQuality);
280  mctrackinfo.SetPDGCode(pdgId);
281  mctrackinfo.SetMomentum(myMcTrack->GetMomentum());
282  mctrackinfo.SetMCQuality(mcStatusMap[mcTrackId]);
283  // mctrackinfo.SetReconstructabilityStatus();
284 
285  if(mctrackinfo.GetNofMCPoints() > 0) {
286  new((*fMCTrackInfo)[size]) PndTrackingQualityMCInfo(mctrackinfo);
287  fMCInfoIdIdealId[idealTrackId] = size;
288  }
289  // cout << "MCTRack " << mctrackinfo.GetMCTrackID() << endl;
290  }
291  // ............................................................
292 
293  // loop over reco track info and associate the mc track info
294  for(int itrk = 0; itrk < fRecoTrackInfo->GetEntriesFast(); itrk++) {
296  Int_t idealTrackId = qaAna.GetIdealTrackIdFromRecoTrackId(recoinfo->GetRecoTrackID());
297  recoinfo->SetIdealTrackId(idealTrackId);
298  if (idealTrackId < 0){
299  std::cout << "-W- PndTrackingQATask::Exec no idealTrack for recoTrack " << itrk << std::endl;
300  continue;
301  }
302 
304  recoinfo->SetMCTrackInfo(mctrackinfo);
305  }
306  // associate reconstructed and mc tracks
308 
309  // Save the Tree (/RhoTuple) for some possible additional analysis
310  for (std::map<Int_t, Int_t>::iterator iter = qualiMap.begin(); iter != qualiMap.end(); iter++) {
311  Int_t mcTrackId = iter->first;
312  Int_t trackQuality = iter->second;
313 
314  if (mcTrackId == -1){
315  std::cout << "-W- PndTrackingQATask::Exec mcTrackId == -1" << std::endl;
316  continue;
317  }
318 
319  TVector3 recoMomentum = recoPMap[mcTrackId];
320 
321  PndMCTrack * myMcTrack = (PndMCTrack *) fMCTrack->At(mcTrackId);
322 
323  if (myMcTrack == 0){
324  std::cout << "-E- PndTrackingQATask::Exec mcMyTrack == 0" << std::endl;
325  continue;
326  }
327  TVector3 mcMomentum = myMcTrack->GetMomentum();
328  Int_t pdgId = myMcTrack->GetPdgCode();
329 
330  fTuple->Column("EvtNr", (Int_t) fEventNr); // number of the currently processed event
331 
332  fTuple->Column("McTrackId", (Int_t) mcTrackId); // id of the currently processed track
333  fTuple->Column("McTrackFoundNTimes", (Int_t) mcFoundMap[mcTrackId]); // A given MC track was found N times
334  fTuple->Column("TrackQuality", (Int_t) trackQuality); // the quality of the current track
335  fTuple->Column("Mc_TrackQuality", (Int_t) mcStatusMap[mcTrackId] - 6); // the mc quality of the current track, -6 to match the global qualityNumbers
336 
337  fTuple->Column("PdgId", (Int_t) pdgId);
338  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
339  fTuple->Column("Mc_py", (Double_t) mcMomentum.Py());
340  fTuple->Column("Mc_pz", (Double_t) mcMomentum.Pz());
341  fTuple->Column("Mc_pt", (Double_t) mcMomentum.Pt());
342  fTuple->Column("Reco_px", (Double_t) recoMomentum.Px());
343  fTuple->Column("Reco_py", (Double_t) recoMomentum.Py());
344  fTuple->Column("Reco_pz", (Double_t) recoMomentum.Pz());
345  fTuple->Column("Reco_pt", (Double_t) recoMomentum.Pt());
346 // fTuple->Column("Reco_pt2", (Double_t) recoPtMap[mcTrackId]); // cross check
347 
348  fTuple->DumpData();
349  }
350 
351  if (fVerbose > 0)
352  qaAna.PrintTrackQualityMap();
353  fEventNr++;
354 }
355 
356 Int_t PndTrackingQATask::GetSumOfAllValidMCHits(FairMultiLinkedData* trackData)
357 {
358  Int_t result = 0;
359  for (size_t branchIndex = 0; branchIndex < fBranchNames.size(); branchIndex++){
360  result += trackData->GetLinksWithType(ioman->GetBranchId(fBranchNames[branchIndex])).GetNLinks();
361  }
362  return result;
363 }
364 
365 
366 void PndTrackingQATask::FillQualyHisto(std::map<Int_t, Int_t> trackQualifikation, Int_t nGhosts)
367 {
368 
369  fQualyHisto->Fill(qualityNumbers::kGhost, nGhosts);
370  for(std::map<Int_t, Int_t>::iterator iter = trackQualifikation.begin(); iter != trackQualifikation.end(); iter++){
371  fQualyHisto->Fill(iter->second);
372  if (iter->second > 0){
374  }
375  else {
377  }
378  }
379 }
380 
381 void PndTrackingQATask::FillMCStatus(std::map<Int_t, Int_t> trackMCStatus)
382 {
384  for(std::map<Int_t, Int_t>::iterator iter = trackMCStatus.begin(); iter != trackMCStatus.end(); iter++){
385  fQualyHisto->Fill(iter->second - mcOffset); // mcOffset = 6
386  }
387 }
388 
389 
390 void PndTrackingQATask::FillEfficiencies(std::map<Int_t, std::map<TString, std::pair<Double_t, Int_t > > > efficiencies)
391 {
392  for (std::map<Int_t, std::map<TString, std::pair<Double_t, Int_t> > >::iterator iterTracks = efficiencies.begin(); iterTracks != efficiencies.end(); iterTracks++ ){
393  std::map<TString, std::pair<Double_t, Int_t> > branchEfficiency = iterTracks->second;
394  for (std::map<TString, std::pair<Double_t, Int_t> >::iterator iterBranch = branchEfficiency.begin(); iterBranch != branchEfficiency.end(); iterBranch++ ){
395  fMapEfficiencies[iterBranch->first]->Fill(iterBranch->second.second, iterBranch->second.first);
396  }
397  }
398 }
399 void PndTrackingQATask::MapToHist (std::map<Int_t, Double_t> map, TH1 * histo) {
400  for (std::map<Int_t, Double_t>::iterator iter = map.begin(); iter != map.end(); iter++ ){
401  histo->Fill(iter->second);
402  }
403 }
404 
405 void PndTrackingQATask::SetQualyHisto(TH1* histo, Bool_t relative, Int_t base)
406 {
407  Int_t allTracks = 0;
408  Int_t allTracksWithHits = 0;
409  Int_t allPossibleTracksWithHits = 0;
410  Int_t allTracksWithHitsNotFound = 0;
411 
412  Int_t mcLessThanThreePrim = fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kMcLessThanThreePrim));
413  Int_t mcAtLeastThreePrim = fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kMcAtLeastThreePrim));
414  Int_t mcAtLeastThreeSec = fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kMcAtLeastThreeSec));
415  Int_t mcPossiblePrim = fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kMcPossiblePrim));
416  Int_t mcPossibleSec = fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kMcPossibleSec));
417 
418  Int_t lessThanThreePrim = fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kLessThanThreePrim));
419  Int_t atLeastThreePrim = fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kAtLeastThreePrim));
420  Int_t atLeastThreeSec = fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kAtLeastThreeSec));
421  Int_t possiblePrim = fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kPossiblePrim));
422  Int_t possibleSec = fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kPossibleSec));
423 
424  Double_t fullyFound = fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kFullyFound));
425  Double_t partiallyFound = fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kPartiallyFound));
426  Double_t spuriousFound = fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kSpuriousFound));
427  Double_t allFound = fullyFound + partiallyFound + spuriousFound;
428 
429  Double_t ghosts = fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kGhost));
430 
431  allTracksWithHits += mcAtLeastThreePrim; // Todo
432  allTracksWithHits += mcAtLeastThreeSec;
433  allTracksWithHits += mcPossiblePrim;
434  allTracksWithHits += mcPossibleSec;
435 
436  allTracks = allTracksWithHits + lessThanThreePrim;
437 
438  allPossibleTracksWithHits += mcPossiblePrim;
439  allPossibleTracksWithHits += mcPossibleSec;
440 
441  allTracksWithHitsNotFound += atLeastThreePrim;
442  allTracksWithHitsNotFound += atLeastThreeSec;
443  allTracksWithHitsNotFound += possiblePrim;
444  allTracksWithHitsNotFound += possibleSec;
445 
446  Double_t divisor = 1.0;
447  if (relative == kTRUE){
448  divisor = allTracks / 100.0;
449  }
450 
451  histo->Fill(qualityNumbers::kMcAllTracksWithHits, (Double_t)allTracksWithHits / divisor);
452  histo->Fill(qualityNumbers::kMcLessThanThreePrim, (Double_t)mcLessThanThreePrim / divisor);
453  histo->Fill(qualityNumbers::kMcAtLeastThreePrim, (Double_t)mcAtLeastThreePrim / divisor);
454  histo->Fill(qualityNumbers::kMcAtLeastThreeSec, (Double_t)mcAtLeastThreeSec / divisor);
455  histo->Fill(qualityNumbers::kMcPossiblePrim, (Double_t)mcPossiblePrim / divisor);
456  histo->Fill(qualityNumbers::kMcPossibleSec, (Double_t)mcPossibleSec / divisor);
457 
458  divisor = 1.0;
459 
460  if (mcLessThanThreePrim > 0){
461  if (relative == kTRUE){
462  divisor = mcLessThanThreePrim / 100.0;
463  }
464  histo->Fill(qualityNumbers::kLessThanThreePrim, (Double_t)(mcLessThanThreePrim - lessThanThreePrim) /divisor);
465  }
466 
467  divisor = 1.0;
468 
469  if (mcAtLeastThreePrim > 0) {
470  if (relative == kTRUE){
471  divisor = mcAtLeastThreePrim / 100.0;
472  }
473  histo->Fill(qualityNumbers::kAtLeastThreePrim, (Double_t)(mcAtLeastThreePrim - atLeastThreePrim) / divisor);
474  }
475 
476  divisor = 1.0;
477 
478  if (mcAtLeastThreeSec > 0){
479  if (relative == kTRUE){
480  divisor = mcAtLeastThreeSec / 100.0;
481  }
482  histo->Fill(qualityNumbers::kAtLeastThreeSec, (Double_t)(Double_t)(mcAtLeastThreeSec - atLeastThreeSec) / divisor);
483  }
484 
485  divisor = 1.0;
486 
487  if (mcPossiblePrim > 0){
488  if (relative == kTRUE){
489  divisor = mcPossiblePrim / 100.0;
490  }
491  histo->Fill(qualityNumbers::kPossiblePrim, (Double_t)(mcPossiblePrim - possiblePrim) / divisor);
492  }
493 
494  divisor = 1.0;
495 
496  if (mcPossibleSec > 0){
497  if (relative == kTRUE){
498  divisor = mcPossibleSec / 100.0;
499  }
500  histo->Fill(qualityNumbers::kPossibleSec, (Double_t)(mcPossibleSec - possibleSec) / divisor);
501  }
502 
503  Double_t baseDouble = base;
504  if (relative == kTRUE){
505  baseDouble /= 100.0;
506  }
507 
508  if (base == 0)
509  return;
510 
511  histo->Fill(qualityNumbers::kFullyFound, fullyFound / baseDouble);
512  histo->Fill(qualityNumbers::kPartiallyFound, partiallyFound / baseDouble);
513  histo->Fill(qualityNumbers::kSpuriousFound, spuriousFound / baseDouble);
514  histo->Fill(qualityNumbers::kGhost, ghosts / baseDouble);
515  histo->Fill(qualityNumbers::kFound, allFound / baseDouble);
516  histo->Fill(qualityNumbers::kNotFound, (base - allFound) / baseDouble);
517 }
518 
520  ColorHistogram();
525  fQualyStack->SetName("fQualyHistoColor");
526  fQualyStack->SetTitle(fQualyHisto->GetTitle());
527 
528  // gROOT->SetBatch(kTRUE);
529  // fQualyHisto->Draw();
530  // LabelQualyHistogram((TH1*) fQualyStack);
531 
532  for (size_t i = 0; i < fBranchNames.size(); i++){
533  fMapEfficiencies[fBranchNames[i]]->Write();
534  }
535  fIdealTracksPerEvent->Write();
536  fIdealPHisto->Write();
537  fIdealPtHisto->Write();
538  fIdealPlHisto->Write();
539 
540  fPHisto->Write();
541  fPRelHisto->Write();
542  fPtHisto->Write();
543  fPtRelHisto->Write();
544  fPlHisto->Write();
545  fPlRelHisto->Write();
546  fQualyHisto->Write();
547  fQualyStack->Write();
548 
549  Int_t allTracks = 0;
550  Int_t allTracksWithHits = 0;
551  Int_t allPossibleTracksWithHits = 0;
552  Int_t allTracksWithHitsNotFound = 0;
553 
554  Int_t mcLessThanThreePrim = fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kMcLessThanThreePrim));
555  Int_t mcAtLeastThreePrim = fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kMcAtLeastThreePrim));
556  Int_t mcAtLeastThreeSec = fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kMcAtLeastThreeSec));
557  Int_t mcPossiblePrim = fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kMcPossiblePrim));
558  Int_t mcPossibleSec = fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kMcPossibleSec));
559 
560  Int_t lessThanThreePrim = fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kLessThanThreePrim));
561  Int_t atLeastThreePrim = fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kAtLeastThreePrim));
562  Int_t atLeastThreeSec = fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kAtLeastThreeSec));
563  Int_t possiblePrim = fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kPossiblePrim));
564  Int_t possibleSec = fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kPossibleSec));
565 
566  allTracksWithHits += mcAtLeastThreePrim; // Todo
567  allTracksWithHits += mcAtLeastThreeSec;
568  allTracksWithHits += mcPossiblePrim;
569  allTracksWithHits += mcPossibleSec;
570 
571  allTracks = allTracksWithHits + lessThanThreePrim;
572 
573  allPossibleTracksWithHits += mcPossiblePrim;
574  allPossibleTracksWithHits += mcPossibleSec;
575 
576  allTracksWithHitsNotFound += atLeastThreePrim;
577  allTracksWithHitsNotFound += atLeastThreeSec;
578  allTracksWithHitsNotFound += possiblePrim;
579  allTracksWithHitsNotFound += possibleSec;
580 
581  Double_t fullyFound = fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kFullyFound));
582  Double_t partiallyFound = fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kPartiallyFound));
583  Double_t spuriousFound = fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kSpuriousFound));
584  //Double_t allFound = fullyFound + partiallyFound + spuriousFound; //[R.K.04/2017] unused variable
585  //Double_t notFound = fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kNotFound)); //[R.K.03/2017] unused variable
586 
587  Double_t ghosts = fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kGhost));
588 
589  std::cout << "fQualyHisto: All Tracks: " << allTracks << std::endl
590  << " Primary Tracks < 3 hits: " << mcLessThanThreePrim << ": Not Found: " << lessThanThreePrim << std::endl
591  << " All Tracks with hits: " << allTracksWithHits << ". Not Found: " << allTracksWithHitsNotFound << std::endl
592  << " Primary Tracks with >= 3 hits, but not a possible track: " << mcAtLeastThreePrim << ". Not Found: " << atLeastThreePrim << std::endl
593  << " Secondary Tracks with >= 3 hits, but not a possible track: " << mcAtLeastThreeSec << ". Not Found: " << atLeastThreeSec << std::endl
594  << " Primary Tracks possible: " << mcPossiblePrim << ". Not Found: " << possiblePrim << std::endl
595  << " Secondary Tracks possible: " << mcPossibleSec << ". Not Found: " << possibleSec << std::endl
596  << " All Possible Tracks with hits: " << allPossibleTracksWithHits << std::endl
597 
598  << " FullyFound: " << fullyFound << " "
599  << fullyFound / allTracksWithHits * 100.0 << "% (of all tracks), "
600  << fullyFound / allPossibleTracksWithHits * 100.0 << "% (of all tracks possible)"
601 
602  << " PartlyFound: " << partiallyFound << " "
603  << partiallyFound / allTracksWithHits * 100.0 << "% "
604  << partiallyFound / allPossibleTracksWithHits * 100.0 << "% "
605 
606  << " Spurious: " << spuriousFound << " "
607  << spuriousFound / allTracksWithHits * 100.0 << "% "
608  << spuriousFound / allPossibleTracksWithHits * 100.0 << "% "
609 
610  << " Ghosts: " << ghosts << " "
611  << ghosts / allTracksWithHits * 100.0 << "% "
612  << ghosts / allPossibleTracksWithHits * 100.0 << "% " << std::endl;
613 
614  SetQualyHisto(fQualyHisto_rel_all, kTRUE, allTracksWithHits);
615  SetQualyHisto(fQualyHisto_rel_possible, kTRUE, allPossibleTracksWithHits);
616 
617  fQualyHisto->SetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kMcAllTracksWithHits), allTracksWithHits);
618  fQualyHisto->SetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kLessThanThreePrim), mcLessThanThreePrim - lessThanThreePrim);
619  fQualyHisto->SetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kAtLeastThreePrim), mcAtLeastThreePrim - atLeastThreePrim);
620  fQualyHisto->SetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kAtLeastThreeSec), mcAtLeastThreeSec - atLeastThreeSec);
621  fQualyHisto->SetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kPossiblePrim), mcPossiblePrim - possiblePrim);
622  fQualyHisto->SetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kPossibleSec), mcPossibleSec - possibleSec);
623 
624  fQualyHisto_rel_all->Write();
625  fQualyHisto_rel_possible->Write();
626  fQualyHisto->Write();
627  fTuple->GetInternalTree()->Write();
628 
629  std::cout << "Finish finished!" << std::endl;
630 }
631 
633  fQualyHisto_pos->SetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kGhost), fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kGhost)));
634  fQualyHisto_pos->SetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kSpuriousFound), fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kSpuriousFound)));
636  fQualyHisto_pos->SetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kFullyFound), fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kFullyFound)));
637  fQualyHisto_pos->SetLineColor(kGreen + 2);
638 
639  fQualyHisto_all->SetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kFound), fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kFound)));
640  fQualyHisto_all->SetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kNotFound), fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kNotFound)));
641  fQualyHisto_all->SetLineColor(kBlack);
642 
643 
644  fQualyHisto_neg->SetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kPossibleSec), fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kPossibleSec)));
645  fQualyHisto_neg->SetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kPossiblePrim), fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kPossiblePrim)));
649  fQualyHisto_neg->SetLineColor(kRed + 2);
650 
651  fQualyHisto_mc->SetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kMcPossibleSec), fQualyHisto->GetBinContent(fQualyHisto->FindFixBin(qualityNumbers::kMcPossibleSec)));
656  fQualyHisto_mc->SetLineColor(kBlue);
657 
658 }
659 
661 
662  PndTrackCand *idealtrkcand = idealtrack->GetTrackCandPtr();
663 
664  //Int_t nofsttpoint = idealtrkcand->GetNHitsDet(FairRootManager::Instance()->GetBranchId("STTHit")); //[R.K. 01/2017] unused variable
665  Int_t nofmvdpixpoint = idealtrkcand->GetNHitsDet(FairRootManager::Instance()->GetBranchId("MVDHitsPixel"));
666  Int_t nofmvdstrpoint = idealtrkcand->GetNHitsDet(FairRootManager::Instance()->GetBranchId("MVDHitsStrip"));
667  //Int_t nofmvdpoint = nofmvdpixpoint + nofmvdstrpoint; //[R.K. 01/2017] unused variable
668  Int_t nofgempoint = idealtrkcand->GetNHitsDet(FairRootManager::Instance()->GetBranchId("GEMHit"));
669  Int_t nofftspoint = idealtrkcand->GetNHitsDet(FairRootManager::Instance()->GetBranchId("FTSHit"));
670 
671  int nofsttskewpoint = 0, nofsttparalpoint = 0;
672  // this loop counts skewed (--> parallel) STT/FTS hits
673  for(size_t ihit = 0; ihit < idealtrkcand->GetNHits(); ihit++) {
674  PndTrackCandHit idealcandhit = idealtrkcand->GetSortedHit(ihit);
675  Int_t hitID = idealcandhit.GetHitId();
676  Int_t detID = idealcandhit.GetDetId();
677 
678  if(detID != FairRootManager::Instance()->GetBranchId("STTHit")) continue;
679  PndSttHit *stthit = (PndSttHit*) fSttHitArray->At(hitID); //todo: For time based sim this has to be replaced with FariRootManager::GetCloneOfLinkData
680  Int_t tubeID = stthit->GetTubeID();
681  PndSttTube *tube = (PndSttTube*) fSttTubeArray->At(tubeID);
682  if(tube->IsSkew()) nofsttskewpoint++;
683  else nofsttparalpoint++;
684  }
685 
686  PndTrackingQualityMCInfo info(nofmvdpixpoint, nofmvdstrpoint, nofsttparalpoint, nofsttskewpoint, nofgempoint, nofftspoint);
687  std::vector<FairLink> mcTracks = idealtrack->GetSortedMCTracks();
688  if (mcTracks.size() > 0){
689  PndMCTrack* myMCTrack = (PndMCTrack*)FairRootManager::Instance()->GetCloneOfLinkData(mcTracks[0]);
690  if (myMCTrack != nullptr){
691  info.SetVertex(myMCTrack->GetStartVertex());
692  if (myMCTrack->GetMotherID() < 0){
693  info.SetIsPrimary(kTRUE);
694  }
695  else{
696  info.SetIsPrimary(kFALSE);
697  }
698  }
699  }
700  // CHECK
701  // Bool_t isreco = Reconstructability(nofmvdpixpoint, nofmvdstrpoint, nofsttparalpoint, nofsttskewpoint, nofgempoint, nofscitilpoint);
702  // info.SetReconstructability(isreco);
703 
704 
705  info.SetPositionFirst(idealtrack->GetParamFirst().GetPosition());
706  info.SetMomentumFirst(idealtrack->GetParamFirst().GetMomentum());
707  info.SetPositionLast(idealtrack->GetParamLast().GetPosition());
708  info.SetMomentumLast(idealtrack->GetParamLast().GetMomentum());
709  info.SetCharge(idealtrack->GetParamFirst().GetQ());
710 
711 
712  return info;
713 
714 }
887  // loop over mc track infos
888  for(int imctrk = 0; imctrk < fMCTrackInfo->GetEntriesFast(); imctrk++) {
890  int mctrackid0 = mcinfo->GetMCTrackID();
891  if(mcinfo->GetAssoRecoTrackID() != -1) continue;
892  double tmpeff = 0., tmppur = 0.;
893  int tmptruerecotrackid = -1;
894 
895  // loop over reco track infos
896  PndTrackingQualityRecoInfo *tmprecoinfo = NULL;
897  for(int itrk = 0; itrk < fRecoTrackInfo->GetEntriesFast(); itrk++) {
899  int mctrackid = recoinfo->GetMCTrackID();
900  if(mctrackid != mctrackid0) continue;
901  mcinfo->SetRecoTrackID(recoinfo->GetRecoTrackID());
902  recoinfo->SetClone();
903  // it must have either the higher efficiency ...
904  if(recoinfo->GetEfficiency() < tmpeff) continue;
905  // ... or, if they are even, the highest purity
906  if(recoinfo->GetEfficiency() == tmpeff && recoinfo->GetPurity() < tmppur) continue;
907 
908  tmpeff = recoinfo->GetEfficiency();
909  tmppur = recoinfo->GetPurity();
910  tmptruerecotrackid = recoinfo->GetRecoTrackID();
911  tmprecoinfo = recoinfo;
912  }
913  if(tmprecoinfo == NULL) continue;
914 
915  tmprecoinfo->SetTrue();
916  mcinfo->SetAssoRecoTrackID(tmptruerecotrackid);
917 
918  }
919 }
920 
921 
922 
923 
virtual void FillQualyHisto(std::map< Int_t, Int_t > trackQualifikation, Int_t nGhosts)
TClonesArray * fRecoTrackInfo
UInt_t GetNHitsDet(Int_t detId)
std::map< Int_t, Int_t > GetMCTrackFound()
int fVerbose
Definition: poormantracks.C:24
void LabelQualyHistogram(TH1 *)
Int_t GetNGhosts()
static const int kMcAtLeastThreeSec
Definition: PndTrackingQA.h:52
Int_t i
Definition: run_full.C:25
std::map< Int_t, Double_t > GetPlResolution()
void AnalyseEvent(TClonesArray *recoTrackInfo)
static const int kGhost
Definition: PndTrackingQA.h:63
void SetMCTrackID(Int_t mctrackid)
PndTransMap * map
Definition: sim_emc_apd.C:99
TClonesArray * fMCTrackInfo
virtual void MapToHist(std::map< Int_t, Double_t >, TH1 *)
Int_t GetPdgCode() const
Definition: PndMCTrack.h:73
virtual void Init()
static const int kLessThanThreePrim
Definition: PndTrackingQA.h:47
void SetHitsBranchNames(std::vector< TString > names)
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
PndTrackingQATask(TString trackBranchName, TString idealBranchName, Bool_t pndTrackData=kTRUE)
std::map< Int_t, Double_t > GetPResolution()
std::map< Int_t, Double_t > GetPResolutionRel()
TVector3 GetMomentum() const
Definition: PndMCTrack.h:78
void SetRecoTrackID(int recotrkid)
FairRootManager * ioman
TClonesArray * fIdealTrack
static const int kPartiallyFound
Definition: PndTrackingQA.h:60
static const int kMcPossibleSec
Definition: PndTrackingQA.h:50
void SetMCTrackInfo(PndTrackingQualityMCInfo *info)
TClonesArray * fSttTubeArray
virtual void FillEfficiencies(std::map< Int_t, std::map< TString, std::pair< Double_t, Int_t > > > efficiencies)
virtual void SetQualyHisto(TH1 *histo, Bool_t relative, Int_t base=1)
TClonesArray * fMCTrack
Int_t GetMCInfoIdFromIdealTrackId(int idealtrackid)
TClonesArray * fSttHitArray
static const int kPossibleSec
Definition: PndTrackingQA.h:43
std::map< Int_t, Int_t > GetTrackMCStatus()
Double_t
PndTrackingQualityMCInfo GetMCInfoFromIdealTrack(PndTrack *idealtrack)
virtual void SetParContainers()
TClonesArray * fTrack
static const int kAtLeastThreePrim
Definition: PndTrackingQA.h:46
static PndTrackFunctor * make_PndTrackFunctor(std::string functorName)
TString fPossibleTrackFunctorName
std::map< Int_t, TVector3 > GetP()
FairTrackParP GetParamLast()
Definition: PndTrack.h:50
FairRuntimeDb * rtdb
Definition: hit_dirc.C:66
Int_t GetTubeID() const
Definition: PndSttHit.h:75
TClonesArray * FillTubeArray()
std::map< Int_t, std::map< TString, std::pair< Double_t, Int_t > > > GetEfficiencies()
void Column(const char *label, Bool_t value, Bool_t defval=0, const char *block=0)
Definition: RhoTuple.cxx:56
std::map< Int_t, Double_t > GetPlResolutionRel()
std::map< Int_t, Double_t > GetPtResolution()
void SetVerbose(Int_t val)
Definition: PndTrackingQA.h:99
static const int kMcLessThanThreePrim
Definition: PndTrackingQA.h:54
std::map< int, int > fMCInfoIdIdealId
ClassImp(PndTrackingQATask)
void DumpData()
Definition: RhoTuple.cxx:391
static const int kMcAllTracksWithHits
Definition: PndTrackingQA.h:55
std::map< Int_t, Int_t > GetTrackQualification()
PndTrackFunctor * fPossibleTrackFunctor
static const int kFullyFound
Definition: PndTrackingQA.h:61
std::vector< TString > fBranchNames
static const int kFound
Definition: PndTrackingQA.h:66
std::map< TString, TH2 * > fMapEfficiencies
Int_t GetIdealTrackIdFromMCTrackId(int mctrackid)
UInt_t GetNHits() const
Definition: PndTrackCand.h:59
static const int kMcAtLeastThreePrim
Definition: PndTrackingQA.h:53
TTree * GetInternalTree()
Definition: RhoTuple.h:207
virtual InitStatus Init()
std::map< Int_t, Double_t > GetPtResolutionRel()
TVector3 GetStartVertex() const
Definition: PndMCTrack.h:76
virtual void Exec(Option_t *opt)
PndTrackCand * GetTrackCandPtr()
Definition: PndTrack.h:48
Int_t GetMotherID() const
Definition: PndMCTrack.h:74
Int_t GetHitId() const
virtual Int_t GetSumOfAllValidMCHits(FairMultiLinkedData *trackData)
TH1F * hist
void AddHitsBranchName(TString name)
Adds branch names of detector data which should be taken into account in the analysis.
void PrintTrackQualityMap(Bool_t detailedInfo=kFALSE)
Int_t GetDetId() const
virtual void FillMCStatus(std::map< Int_t, Int_t > trackMCStatus)
static const int kMcPossiblePrim
Definition: PndTrackingQA.h:51
PndGeoSttPar * fSttParameters
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
Int_t GetIdealTrackIdFromRecoTrackId(int trackid)