FairRoot/PandaRoot
PndSttCellTrackFinderAnalysisTask.cxx
Go to the documentation of this file.
1 /*
2  * PndSttCellTrackFinderAnalysisTask.cpp
3  *
4  * Created on: Jun 26, 2013
5  * Author: schumann
6  */
7 
8 // Root includes
9 #include "TROOT.h"
10 #include "TClonesArray.h"
11 
12 // framework includes
13 #include "FairRunAna.h"
14 #include "FairRootManager.h"
15 #include "FairRun.h"
16 #include "FairRuntimeDb.h"
17 #include "FairHit.h"
18 #include "FairEventHeader.h"
19 
20 #include "PndSttCellTrackFinderAnalysisTask.h"
21 #include "PndTrackCand.h"
22 #include "PndTrackCandHit.h"
23 #include "PndRiemannTrack.h"
24 #include "PndMCMatch.h"
25 #include "PndMCEntry.h"
26 #include "PndMCResult.h"
27 #include "PndMCTrack.h"
28 #include "PndSttMapCreator.h"
29 #include "PndSttHit.h"
30 #include "PndGeoSttPar.h"
31 #include "PndSttTube.h"
32 
33 #include "TGraph.h"
34 #include "TArc.h"
35 #include "TLine.h"
36 #include "TMath.h"
37 
38 #include <iostream>
39 #include <vector>
40 #include <set>
41 #include <string>
42 
43 ClassImp(PndSttCellTrackFinderAnalysisTask)
44 ;
45 
46 using namespace std;
47 
48 void PndSttCellTrackFinderAnalysisTask::SetParContainers() {
49  FairRuntimeDb* rtdb = FairRunAna::Instance()->GetRuntimeDb();
50  fSttParameters = (PndGeoSttPar*) rtdb->getContainer("PndGeoSttPar");
51 }
52 
53 InitStatus PndSttCellTrackFinderAnalysisTask::Init() {
54 
55  FairRootManager* ioman = FairRootManager::Instance();
56 
57  if (!ioman) {
58  cout << "-E- PndSttCellTrackFinderAnalysisTask::Init: "
59  << "RootManager not instantiated!" << endl;
60  return kFATAL;
61  }
62 
63  fMCMatch = (PndMCMatch*) ioman->GetObject("MCMatch");
64 
65  if (!fMCMatch) {
66  cout
67  << "-W- PndSttCellTrackFinderAnalysisTask::Init: No MCMatch array! Needed for MC Truth"
68  << endl;
69  return kERROR;
70  }
71 
72  fEventHeader = (TClonesArray*) ioman->GetObject("EventHeader.");
73  if (!fEventHeader) {
74  cout
75  << "-W- PndSttCellTrackFinderAnalysisTask::Init: No EventHeader array! Needed for EventNumber"
76  << endl;
77  return kERROR;
78  }
79 
80  // get data-object
81  fMCTrack = (TClonesArray*) ioman->GetObject("MCTrack");
82  fFirstTrackCand = (TClonesArray*) ioman->GetObject("FirstTrackCand");
83  fFirstRiemannTrack = (TClonesArray*) ioman->GetObject("FirstRiemannTrack");
84  fCombiTrackCand = (TClonesArray*) ioman->GetObject("CombiTrackCand");
85  fCombiRiemannTrack = (TClonesArray*) ioman->GetObject("CombiRiemannTrack");
86  fSTTHits = (TClonesArray*) ioman->GetObject("STTHit");
87 
88  // get tube-array
89  PndSttMapCreator *mapper = new PndSttMapCreator(fSttParameters);
90  fTubeArray = mapper->FillTubeArray();
91 
92  // Initializing histograms
93  fHistoNumberOfLinks1 =
94  new TH1I("fHistoNumberOfLinks1",
95  "Number Of Links Per Track Candidate ; # Links ; # Track Candidates ",
96  5, -0.5, 4.5);
97 
98  fHistoNumberOfTracklets1 = new TH1I("fHistoNumberOfTracklets1",
99  "Number Of Tracklets Per MCTrack; #Tracklets ; # MCTracks ", 20,
100  -0.5, 19.5);
101 
102  fHistoNumberOfAssignedHits1 = new TH1I("fHistoNumberOfAssignedHits1",
103  "Ratio Of Assigned And Unassigned STTHits ; ; # Hits ", 2, -0.5,
104  1.5);
105 
106  fHistoNumberOfHits1 = new TH1I("fHistoNumberOfHits1",
107  "Number Of Hits Per Tracklet ; # Hits ; # Tracklets ", 40, -0.5,
108  39.5);
109 
110  fHistoNumberOfErrors1 = new TH1I("fHistoNumberOfErrors1",
111  "Number Of Errors Per Event ; # Errors ; # Events ", 11, -1.5,
112  9.5);
113 
114  fHistoNumberOfLinks2 =
115  new TH1I("fHistoNumberOfLinks2",
116  "Number Of Links Per Track Candidate ; # Links ; # Track Candidates ",
117  5, -0.5, 4.5);
118 
119  fHistoNumberOfTracklets2 = new TH1I("fHistoNumberOfTracklets2",
120  "Number Of Tracklets Per MCTrack; # Tracklets ; # MCTracks ", 25,
121  -0.5, 24.5);
122 
123  fHistoNumberOfHits2 = new TH1I("fHistoNumberOfHits2",
124  "Number Of Hits Per Tracklet ; # Hits ; # Tracklets ", 40, -0.5,
125  39.5);
126 
127  fHistoNumberOfErrors2 = new TH1I("fHistoNumberOfErrors2",
128  "Number Of Errors Per Event ; # Errors ; # Events ", 11, -1.5,
129  9.5);
130 
131  fHistoQualityCombi = new TH1I("fHistoQualityCombi",
132  "Quality Of Trackfinding; ; # Tracks", 8, -0.5, 7.5);
133 
134  fHistoQualityFirstStep = new TH1I("fHistoQualityFirstStep",
135  "Quality Of Trackfinding; ; # Tracks", 8, -0.5, 7.5);
136 
137  fHistoMaxHitsOfMCTrack1 = new TH1I("fHistoMaxHitsOfMCTrack1",
138  "Max Number Of Hits At A Stretch Per MCTrack; # Hits; # MCTracks",
139  40, -0.5, 39.5);
140 
141  fHistoMaxHitsOfMCTrack2 = new TH1I("fHistoMaxHitsOfMCTrack2",
142  "Max Number Of Hits At A Stretch Per MCTrack; # Hits; # MCTracks",
143  40, -0.5, 39.5);
144 
145  fHistoNumberOfMissingHits1 =
146  new TH1I("fHistoNumberOfMissingHits1",
147  "Amount Of Missing Hits Of Spurious Tracklets; Missing Hits in %; # Tracklets",
148  4, -0.5, 3.5);
149 
150  fHistoNumberOfMissingHits2 =
151  new TH1I("fHistoNumberOfMissingHits2",
152  "Amount Of Missing Hits Of Spurious Tracklets; Missing Hits in %; # Tracklets",
153  4, -0.5, 3.5);
154 
155  // set labels
156  fHistoQualityFirstStep->GetXaxis()->SetBinLabel(1, "#MCTracks");
157  fHistoQualityFirstStep->GetXaxis()->SetBinLabel(2, "#Found in pure");
158  fHistoQualityFirstStep->GetXaxis()->SetBinLabel(3, "#Found in faulty");
159  fHistoQualityFirstStep->GetXaxis()->SetBinLabel(4, "#Not Found");
160  fHistoQualityFirstStep->GetXaxis()->SetBinLabel(5, "#TrackCands");
161  fHistoQualityFirstStep->GetXaxis()->SetBinLabel(6, "#Spuriuous Cands");
162  fHistoQualityFirstStep->GetXaxis()->SetBinLabel(7, "#FullyReco Cands");
163  fHistoQualityFirstStep->GetXaxis()->SetBinLabel(8, "#Ghost Cands");
164 
165  fHistoQualityCombi->GetXaxis()->SetBinLabel(1, "#MCTracks");
166  fHistoQualityCombi->GetXaxis()->SetBinLabel(2, "#Found in pure");
167  fHistoQualityCombi->GetXaxis()->SetBinLabel(3, "#Found in faulty");
168  fHistoQualityCombi->GetXaxis()->SetBinLabel(4, "#Not Found");
169  fHistoQualityCombi->GetXaxis()->SetBinLabel(5, "#TrackCands");
170  fHistoQualityCombi->GetXaxis()->SetBinLabel(6, "#Spuriuous Cands");
171  fHistoQualityCombi->GetXaxis()->SetBinLabel(7, "#FullyReco Cands");
172  fHistoQualityCombi->GetXaxis()->SetBinLabel(8, "#Ghost Cands");
173 
174  fHistoNumberOfMissingHits1->GetXaxis()->SetBinLabel(1, ">0");
175  fHistoNumberOfMissingHits1->GetXaxis()->SetBinLabel(2, ">25");
176  fHistoNumberOfMissingHits1->GetXaxis()->SetBinLabel(3, ">50");
177  fHistoNumberOfMissingHits1->GetXaxis()->SetBinLabel(4, ">75");
178 
179  fHistoNumberOfMissingHits2->GetXaxis()->SetBinLabel(1, ">0");
180  fHistoNumberOfMissingHits2->GetXaxis()->SetBinLabel(2, ">25");
181  fHistoNumberOfMissingHits2->GetXaxis()->SetBinLabel(3, ">50");
182  fHistoNumberOfMissingHits2->GetXaxis()->SetBinLabel(4, ">75");
183 
184  fHistoNumberOfHitsMCTrack = new TH1I("fHistoNumberOfHitsMCTrack",
185  "Number Of Hits Per MCTrack; # Hits; # MCTracks", 70, -0.5, 69.5);
186 
187  fCanvas = new TCanvas();
188  fCanvas->SetName("fCanvas");
189  fCanvas->SetGrid();
190  fCanvas->SetCanvasSize(700, 700);
191 
192  if (fVerbose == 2) {
193  cout
194  << "-I- PndSttCellTrackFinderAnalysisTask::Init: Initialization successfull"
195  << endl;
196  }
197 
198  return kSUCCESS;
199 
200 }
201 
202 void PndSttCellTrackFinderAnalysisTask::Exec(Option_t *opt) {
203 
204  if (fVerbose > 0) {
205  cout << "============= Begin PndSttCellTrackFinderAnalysisTask::Exec"
206  << endl;
207  cout << endl;
208  }
209 
210  // init fMapTubeIDToHits
211  PndSttHit* hit;
212  for (int j = 0; j < fSTTHits->GetEntriesFast(); ++j) {
213  hit = (PndSttHit*) fSTTHits->At(j);
214  fMapTubeIDToHits[hit->GetTubeID()].push_back(
215  hit->GetLink(0).GetIndex());
216  fMapHitIndexToTubeID[hit->GetLink(0).GetIndex()] = hit->GetTubeID();
217  }
218 
219  fMCMatch->CreateArtificialStage("MCTrack");
220 
221  CheckFirstTracklets();
222  TestRecoQualityFirstStep();
223  CountMaxHitsFirstStep();
224 
225  CheckTrackletCombinations();
226  TestRecoQualityCombi();
227  CountMaxHitsCombi();
228 
229 }
230 
231 void PndSttCellTrackFinderAnalysisTask::CheckFirstTracklets() {
232 
233  if (fVerbose > 1)
234  cout << "=== CheckFirstTracklets() ===" << endl;
235 
236  //PndMCTrack* myMCTrack; //[R.K. 01/2017] unused variable?
237  PndTrackCand* myTrackCand;
238  vector<PndTrackCandHit> hits;
239  int eventNumber;
240 
241  FairRootManager* ioman = FairRootManager::Instance();
242  FairEventHeader* myEventHeader = (FairEventHeader*) fEventHeader;
243 
244  // get FairLinks of stored data
245  PndMCResult myResultCandToMCTrack = fMCMatch->GetMCInfo("FirstTrackCand",
246  "MCTrack");
247  PndMCResult myResultHitsToMCTrack = fMCMatch->GetMCInfo("STTHit",
248  "MCTrack");
249  PndMCResult myResultHitToPoint = fMCMatch->GetMCInfo("STTHit", "STTPoint");
250 
251  eventNumber = myEventHeader->GetMCEntryNumber();
252 
253  if (fVerbose > 2) {
254  cout << "EventNumber: " << eventNumber << endl;
255  cout << "# TrackCand: " << myResultCandToMCTrack.GetNEntries() << endl;
256  }
257 
258  // variables for histograms
259  map<int, int> trackletsCounter; // map<MCTrack-Index, #Tracklets>
260  int assignedHits = 0;
261  int errorCounter = 0;
262 
263  int numberOfLinks;
264 
265  // check tracklets
266  for (int i = 0; i < fFirstTrackCand->GetEntriesFast(); ++i) {
267 
268  // increase tracklet-counter for actual MCTrack-Index
269  for (int j = 0; j < myResultCandToMCTrack.GetEntry(i).GetNLinks();
270  ++j) {
271  trackletsCounter[myResultCandToMCTrack.GetEntry(i).GetLink(j).GetIndex()] +=
272  1;
273  }
274 
275  // get current TrackCand
276  myTrackCand = (PndTrackCand*) fFirstTrackCand->At(i);
277 
278  // get all FairLinks of STTHits
279  FairMultiLinkedData hitLinks = myTrackCand->GetLinksWithType(
280  ioman->GetBranchId("STTHit"));
281 
282  // add number of STTHits of the tracklet to the assigned hits
283  assignedHits += hitLinks.GetNLinks();
284  fHistoNumberOfHits1->Fill(hitLinks.GetNLinks());
285 
286  // clean tracklet?
287  if (myResultCandToMCTrack.GetEntry(i).GetNLinks() > 1) {
288  // hits of tracklet point to more than one MCTrack
289 
290  // there could be more than one hit per tube -> get "real" number of different links
291  numberOfLinks = GetNumLinksOfHits(hitLinks, myResultHitsToMCTrack,
292  myResultHitToPoint);
293 
294  // fill histo with number of elements in the set (only different ones were stored)
295  fHistoNumberOfLinks1->Fill(numberOfLinks);
296 
297  if (numberOfLinks > 1) {
298  // faulty tracklet
299  ++errorCounter;
300 
301  if (fVerbose > 2) {
302  cout << "Warning: Found more than one TrackID!" << endl;
303  cout << "Nlinks:"
304  << myResultCandToMCTrack.GetEntry(i).GetNLinks()
305  << endl;
306  }
307  }
308 
309  } else {
310  // only one link to MCTrack
311  fHistoNumberOfLinks1->Fill(1);
312  }
313  }
314 
315  // store data of tracklet-counter in histogram
316  map<int, int>::iterator it;
317  for (it = trackletsCounter.begin(); it != trackletsCounter.end(); ++it) {
318  fHistoNumberOfTracklets1->Fill(it->second);
319  }
320 
321  // store data of assigned and unassigned hits in histogram
322  for (int i = 0; i < assignedHits; ++i) {
323  fHistoNumberOfAssignedHits1->Fill("Assigned", 1);
324  }
325  for (int i = 0; i < fSTTHits->GetEntriesFast() - assignedHits; ++i) {
326  fHistoNumberOfAssignedHits1->Fill("Unassigned", 1);
327  }
328 
329  // store errcorCounter of this event in histogram
330  fHistoNumberOfErrors1->Fill(errorCounter);
331 
332 // // plot STTHits + Riemann-circle for the first 100 events
333 // if (eventNumber < 100) {
334 // DrawFirstRiemannPlots(eventNumber);
335 // }
336 }
337 
338 void PndSttCellTrackFinderAnalysisTask::CheckTrackletCombinations() {
339 
340  if (fVerbose > 1)
341  cout << "=== CheckTrackletCombinations() ===" << endl;
342  //PndMCTrack* myMCTrack; //[R.K. 01/2017] unused variable?
343  PndTrackCand* myTrackCand;
344  vector<PndTrackCandHit> hits;
345  int eventNumber;
346 
347  FairRootManager* ioman = FairRootManager::Instance();
348  FairEventHeader* myEventHeader = (FairEventHeader*) fEventHeader;
349 
350  // get FairLinks of stored data
351  PndMCResult myResultCandToMCTrack = fMCMatch->GetMCInfo("CombiTrackCand",
352  "MCTrack");
353  PndMCResult myResultHitsToMCTrack = fMCMatch->GetMCInfo("STTHit",
354  "MCTrack");
355  PndMCResult myResultHitToPoint = fMCMatch->GetMCInfo("STTHit", "STTPoint");
356 
357  eventNumber = myEventHeader->GetMCEntryNumber();
358 
359  if (fVerbose > 2) {
360  cout << "EventNumber: " << eventNumber << endl;
361  cout << "# TrackCand: " << myResultCandToMCTrack.GetNEntries() << endl;
362  }
363 
364  // variables for histogram
365  map<int, int> trackletsCounter; // map<MCTrack-Index, #Tracklets>
366  int errorCounter = 0;
367 
368  int numberOfLinks;
369 
370  // check each tracklet
371  for (int i = 0; i < fCombiTrackCand->GetEntriesFast(); ++i) {
372 
373  // increase tracklet-counter for actual MCTrack-Index
374  for (int j = 0; j < myResultCandToMCTrack.GetEntry(i).GetNLinks();
375  ++j) {
376  trackletsCounter[myResultCandToMCTrack.GetEntry(i).GetLink(j).GetIndex()] +=
377  1;
378  }
379 
380  // get current TrackCand
381  myTrackCand = (PndTrackCand*) fCombiTrackCand->At(i);
382 
383  // get all FairLinks of STTHits
384  FairMultiLinkedData hitLinks = myTrackCand->GetLinksWithType(
385  ioman->GetBranchId("STTHit"));
386 
387  // fill histo with number of hits of this trackCand
388  fHistoNumberOfHits2->Fill(hitLinks.GetNLinks());
389 
390  // clean tracklet?
391  if (myResultCandToMCTrack.GetEntry(i).GetNLinks() > 1) {
392 
393  // there could be more than one hit per tube -> get "real" number of different links
394  numberOfLinks = GetNumLinksOfHits(hitLinks, myResultHitsToMCTrack,
395  myResultHitToPoint);
396 
397  // fill histogram with number of elements in the set (only different ones were stored)
398  fHistoNumberOfLinks2->Fill(numberOfLinks);
399 
400  if (numberOfLinks > 1) {
401  // faulty tracklet
402  ++errorCounter;
403 
404  if (fVerbose > 2) {
405  cout << "Warning: Found more than one TrackID!" << endl;
406  cout << "Nlinks:"
407  << myResultCandToMCTrack.GetEntry(i).GetNLinks()
408  << endl;
409  }
410  }
411 
412  } else {
413  // only one link to MCTrack
414  fHistoNumberOfLinks2->Fill(1);
415  }
416  }
417 
418  // store data of tracklet-counter in histogram
419  map<int, int>::iterator it;
420  for (it = trackletsCounter.begin(); it != trackletsCounter.end(); ++it) {
421  fHistoNumberOfTracklets2->Fill(it->second);
422  }
423 
424  // store errorCounter ot this event in histogramm
425  fHistoNumberOfErrors2->Fill(errorCounter);
426 
427 // // draw STTHits + Riemann-circle for the first 100 events
428 // if (eventNumber < 100) {
429 // DrawCombiRiemannPlots(eventNumber);
430 // }
431 }
432 
433 int PndSttCellTrackFinderAnalysisTask::GetNumLinksOfHits(
434  FairMultiLinkedData& hitLinks, PndMCResult& sttHitsToMCTrack,
435  PndMCResult& sttHitToPoint) {
436 
437  if (fVerbose > 3)
438  cout << "GetNumLinksOfHits()" << endl;
439 
440  /* Check hits of a tracklet in hitLinks. Approach: check MCTrack-index of consecutive hits.
441  * If they are not equal: tracklet might be faulty -> search for another hit in the same tube with
442  * the faulty hit (if available) and get index of MCTrack. Count the different MCTrack-Indices.*/
443 
444  /* Set containing MCTrack-Indices of STTHits of one PndTrackCand.
445  * If it contains more than one Index, the tracklet ist faulty. */
446  set<int> trackIndices;
447 
448  /* PriorTrackIndex and curTrackIndex of prior and current hit in tracklet should be the same.*/
449  int priorTrackIndex;
450  int curTrackIndex;
451  int priorHitIndex;
452  int curHitIndex;
453 
454  int tmpHitIndex;
455  int tubeID;
456  FairLink curLink;
457 
458  // get MCTrackIndex the first hit points to
459  priorHitIndex = hitLinks.GetLink(0).GetIndex();
460  priorTrackIndex =
461  sttHitsToMCTrack.GetEntry(priorHitIndex).GetLink(0).GetIndex();
462 
463  //only one hit?
464  if (hitLinks.GetNLinks() < 2) {
465  return sttHitsToMCTrack.GetEntry(priorHitIndex).GetNLinks();
466  }
467 
468  if (fVerbose > 3)
469  cout << "TubeID: " << fMapHitIndexToTubeID[priorHitIndex]
470  << ", TrackID: " << priorTrackIndex << endl;
471 
472  // check each hit of myTrackCand: compare MCTrackIndex of prior and current hit
473  for (int j = 1; j < hitLinks.GetNLinks(); ++j) {
474 
475  // get MCTrackIndex the current Hit points to
476  curLink = hitLinks.GetLink(j);
477  // !! index of STTPoint and STTHit must not be the same !!
478  curHitIndex =
479  sttHitToPoint.GetEntry(curLink.GetIndex()).GetLink(0).GetIndex();
480  curTrackIndex =
481  sttHitsToMCTrack.GetEntry(curHitIndex).GetLink(0).GetIndex();
482 
483  if (fVerbose > 3)
484  cout << "TubeID: " << fMapHitIndexToTubeID[curHitIndex]
485  << ", TrackID: " << curTrackIndex << endl;
486 
487  // consecutive hits doesn't point to same MCTrack?
488  if (priorTrackIndex != curTrackIndex) {
489 
490  // more than one hit per tube?
491  if (fMapTubeIDToHits[fMapHitIndexToTubeID[curHitIndex]].size()
492  > 1) {
493  // tube of current hit has more than one hit
494 
495  tubeID = fMapHitIndexToTubeID[curHitIndex];
496 
497  // get hitIndex of the other hit
498  if (curHitIndex == fMapTubeIDToHits[tubeID].at(0)) {
499  tmpHitIndex = fMapTubeIDToHits[tubeID].at(1);
500  } else {
501  tmpHitIndex = fMapTubeIDToHits[tubeID].at(0);
502  }
503 
504  // insert trackIndex of the other hit into set
505  trackIndices.insert(
506  sttHitsToMCTrack.GetEntry(tmpHitIndex).GetLink(0).GetIndex());
507 
508  } else if (fMapTubeIDToHits[fMapHitIndexToTubeID[priorHitIndex]].size()
509  > 1) {
510  // tube of prior hit has more than one hit
511 
512  tubeID = fMapHitIndexToTubeID[priorHitIndex];
513 
514  // get hitIndex of the other hit
515  if (priorHitIndex == fMapTubeIDToHits[tubeID].at(0)) {
516  tmpHitIndex = fMapTubeIDToHits[tubeID].at(1);
517  } else {
518  tmpHitIndex = fMapTubeIDToHits[tubeID].at(0);
519  }
520 
521  // insert trackIndex of the other prior hit into set
522  trackIndices.insert(
523  sttHitsToMCTrack.GetEntry(tmpHitIndex).GetLink(0).GetIndex());
524  // inset trackIndex of current hit
525  trackIndices.insert(curTrackIndex);
526 
527  } else {
528  // only one hit per tube -> insert current trackIndex into set
529  trackIndices.insert(curTrackIndex);
530 
531  if (fVerbose > 3)
532  cout << "---------> TubeID of faulty hit: "
533  << fMapHitIndexToTubeID[curHitIndex] << endl;
534  }
535 
536  } else {
537  // same trackIndices
538  trackIndices.insert(curTrackIndex);
539  }
540 
541  // update data of prior hit
542  priorTrackIndex = curTrackIndex;
543  priorHitIndex = curHitIndex;
544  }
545 
546  return trackIndices.size();
547 }
548 
549 void PndSttCellTrackFinderAnalysisTask::DrawFirstRiemannPlots(int eventNumber) {
550 
551  PndMCResult myResultHitToPoint = fMCMatch->GetMCInfo("STTHit", "STTPoint");
552  PndTrackCand* myTrackCand;
553  int curHitIndex;
554 
555  fCanvas->Clear();
556  fCanvas->SetTitle(
557  "STTHits and circle fits for event #"
558  + TString::Itoa(eventNumber, 10));
559 
560  TGraph* graph = new TGraph();
561  int hitCounter = 0; // for counting points of TGraph
562 
563  /* Draw points of all STTHits of each TrackCand.
564  * Hits of skewed tubes were ignored.*/
565  for (int i = 0; i < fFirstTrackCand->GetEntriesFast(); ++i) {
566 
567  // get current TrackCand
568  myTrackCand = (PndTrackCand*) fFirstTrackCand->At(i);
569  // get all FairLinks of STTHits
570  FairMultiLinkedData hitLinks = myTrackCand->GetLinksWithType(
571  FairRootManager::Instance()->GetBranchId("STTHit"));
572 
573  // for each STTHit
574  for (int j = 0; j < hitLinks.GetNLinks(); ++j) {
575 
576  // Get HitIndex -> get Tube -> Get Position
577  curHitIndex = hitLinks.GetLink(j).GetIndex();
578  PndSttTube* tube = (PndSttTube*) fTubeArray->At(
579  fMapHitIndexToTubeID[curHitIndex]);
580  TVector3 pos = tube->GetPosition();
581 
582  //ignore skewed straw tubes
583  if (!tube->IsSkew()) {
584  graph->SetPoint(hitCounter, pos[0], pos[1]);
585  ++hitCounter;
586  }
587 
588  }
589 
590  }
591 
592  // draw graph
593  graph->Draw("SAMEA*");
594  graph->GetXaxis()->SetLimits(-50, 50);
595  graph->GetYaxis()->SetRangeUser(-50, 50);
596 
597  // draw circles of riemannTracks
598 
599  PndRiemannTrack* myRiemannTrack;
600  double r, x, y;
601 
602  // for each riemannTrack
603  for (int i = 0; i < fFirstRiemannTrack->GetEntriesFast(); ++i) {
604  myRiemannTrack = (PndRiemannTrack*) fFirstRiemannTrack->At(i);
605 
606  // get radius + position of center
607  r = myRiemannTrack->r();
608  x = myRiemannTrack->orig()[0];
609  y = myRiemannTrack->orig()[1];
610 
611  // draw circle
612  TArc* circle = new TArc(x, y, r);
613  circle->SetFillStyle(0);
614  circle->Draw("SAME");
615  }
616 
617  // mark center of STT
618  TLine* xLine = new TLine(-50, 0, 50, 0);
619  TLine* yLine = new TLine(0, 50, 0, -50);
620  xLine->Draw("SAME");
621  yLine->Draw("SAME");
622 
623  // save plot
624  fCanvas->Print(
625  "Plots/FirstTracklets/Event" + TString::Itoa(eventNumber, 10)
626  + ".pdf", "pdf");
627 }
628 
629 void PndSttCellTrackFinderAnalysisTask::DrawCombiRiemannPlots(int eventNumber) {
630 
631  PndMCResult myResultHitToPoint = fMCMatch->GetMCInfo("STTHit", "STTPoint");
632  PndTrackCand* myTrackCand;
633  int curHitIndex;
634 
635  fCanvas->Clear();
636  fCanvas->SetTitle(
637  "STTHits and circle fits for event #"
638  + TString::Itoa(eventNumber, 10));
639 
640  TGraph* graph = new TGraph();
641  int hitCounter = 0; // for counting the points of TGraph
642 
643  /* Draw points of all STTHits of each TrackCand.
644  * Hits of skewed tubes were ignored.*/
645  for (int i = 0; i < fCombiTrackCand->GetEntriesFast(); ++i) {
646 
647  // get current TrackCand
648  myTrackCand = (PndTrackCand*) fCombiTrackCand->At(i);
649  // get all FairLinks of STTHits
650  FairMultiLinkedData hitLinks = myTrackCand->GetLinksWithType(
651  FairRootManager::Instance()->GetBranchId("STTHit"));
652 
653  // for each STTHit
654  for (int j = 0; j < hitLinks.GetNLinks(); ++j) {
655 
656  // get HitIndex -> get Tube -> get Position
657  curHitIndex = hitLinks.GetLink(j).GetIndex();
658  PndSttTube* tube = (PndSttTube*) fTubeArray->At(
659  fMapHitIndexToTubeID[curHitIndex]);
660  TVector3 pos = tube->GetPosition();
661 
662  //ignore skewed straw tubes
663  if (!tube->IsSkew()) {
664  graph->SetPoint(hitCounter, pos[0], pos[1]);
665  ++hitCounter;
666  }
667 
668  }
669 
670  }
671 
672  // draw graph
673  graph->Draw("SAMEA*");
674  graph->GetXaxis()->SetLimits(-50, 50);
675  graph->GetYaxis()->SetRangeUser(-50, 50);
676 
677  // draw circles of riemannTracks
678  PndRiemannTrack* myRiemannTrack;
679  double r, x, y;
680 
681  // for each riemannTrack
682  for (int i = 0; i < fCombiRiemannTrack->GetEntriesFast(); ++i) {
683  myRiemannTrack = (PndRiemannTrack*) fCombiRiemannTrack->At(i);
684 
685  // get radius + position of center
686  r = myRiemannTrack->r();
687  x = myRiemannTrack->orig()[0];
688  y = myRiemannTrack->orig()[1];
689 
690  // draw circle
691  TArc* circle = new TArc(x, y, r);
692  circle->SetFillStyle(0);
693  circle->Draw("SAME");
694  }
695 
696  // marl center of STT
697  TLine* xLine = new TLine(-50, 0, 50, 0);
698  TLine* yLine = new TLine(0, 50, 0, -50);
699  xLine->Draw("SAME");
700  yLine->Draw("SAME");
701 
702  // save plot
703  fCanvas->Print(
704  "Plots/Combinations/EventCombi" + TString::Itoa(eventNumber, 10)
705  + ".pdf", "pdf");
706 }
707 
708 void PndSttCellTrackFinderAnalysisTask::CheckConformalMap() {
709 
710  TCanvas* test = new TCanvas("test", "Circle");
711  test->SetGrid();
712  TGraph* circlePoints = new TGraph();
713  TGraph* linePoints = new TGraph();
714  int pointCounter = 0;
715 
716  double strawHits[6][2] = { { -1, 2 }, { 0, 3 }, { 1, 2 }, { 2, 3 },
717  { 3, 2 }, { 4, 3 } };
718  double mappedHits[6][2];
719  double div;
720 
721  double refHitX = strawHits[0][0];
722  double refHitY = strawHits[0][1];
723 
724  for (int i = 0; i < 6; ++i) {
725 
726  div = (strawHits[i][0] - refHitX) * (strawHits[i][0] - refHitX)
727  + (strawHits[i][1] - refHitY) * (strawHits[i][1] - refHitY);
728  mappedHits[i][0] = (strawHits[i][0] - refHitX) / div;
729  mappedHits[i][1] = -(strawHits[i][1] - refHitY) / div;
730  circlePoints->SetPoint(i, strawHits[i][0], strawHits[i][1]);
731  linePoints->SetPoint(i, mappedHits[i][0], mappedHits[i][1]);
732  }
733 
734  circlePoints->Draw("SAMEA*");
735  circlePoints->GetXaxis()->SetLimits(-2, 5);
736  circlePoints->GetYaxis()->SetRangeUser(-2, 5);
737  linePoints->Draw("*");
738  linePoints->SetMarkerColor(2);
739 
740  test->SetCanvasSize(700, 700);
741  test->Print("Plots/CircleLine1.pdf", "pdf");
742 
743  //Second example
744  test->Clear();
745 
746  TGraph* pointsOnCircle = new TGraph();
747  TGraph* pointsOnLine = new TGraph();
748 
749  double radius = 2;
750  double centerX = 1;
751  double centerY = 1;
752 
753  double oldPoints[9][2];
754  double mappedPoints[9][2];
755 
756  pointCounter = 0;
757 
758  //calc circle-points
759  for (double y = -1; y <= 3; y += 0.5) {
760 
761  oldPoints[pointCounter][0] = TMath::Sqrt(
762  radius * radius - ((y - centerY) * (y - centerY))) + centerX;
763  oldPoints[pointCounter][1] = y;
764  pointsOnCircle->SetPoint(pointCounter, oldPoints[pointCounter][0],
765  oldPoints[pointCounter][1]);
766  ++pointCounter;
767  }
768 
769  //calc with error
770  double errPoint[2] = { -0.5, 3 };
771  pointsOnCircle->SetPoint(pointCounter, errPoint[0], errPoint[1]);
772 
773  //calc mapped points, take one old point as reference point
774  double refX = oldPoints[0][0];
775  double refY = oldPoints[0][1];
776  pointCounter = 0;
777 
778  for (int i = 0; i < 9; ++i) {
779  div = (oldPoints[pointCounter][0] - refX)
780  * (oldPoints[pointCounter][0] - refX)
781  + (oldPoints[pointCounter][1] - refY)
782  * (oldPoints[pointCounter][1] - refY);
783  mappedPoints[pointCounter][0] = (oldPoints[pointCounter][0] - refX)
784  / div;
785  mappedPoints[pointCounter][1] = -(oldPoints[pointCounter][1] - refY)
786  / div;
787  pointsOnLine->SetPoint(pointCounter, mappedPoints[pointCounter][0],
788  mappedPoints[pointCounter][1]);
789  ++pointCounter;
790  }
791 
792  double mappedErrPoint[2];
793  mappedErrPoint[0] = (errPoint[0] - refX)
794  / ((errPoint[0] - refX) * (errPoint[0] - refX)
795  + (errPoint[1] - refY) * (errPoint[1] - refY));
796  mappedErrPoint[1] = -(errPoint[1] - refY)
797  / ((errPoint[0] - refX) * (errPoint[0] - refX)
798  + (errPoint[1] - refY) * (errPoint[1] - refY));
799  pointsOnLine->SetPoint(pointCounter, mappedErrPoint[0], mappedErrPoint[1]);
800 
801  pointsOnCircle->Draw("A*");
802  pointsOnCircle->GetXaxis()->SetLimits(-2, 4);
803  pointsOnCircle->GetYaxis()->SetRangeUser(-2, 4);
804 
805  TArc* circle2 = new TArc(centerX, centerY, radius);
806  circle2->SetFillStyle(0);
807  circle2->Draw("SAME");
808 
809  pointsOnLine->Draw("SAME*");
810 
811  test->Print("Plots/CircleLine2.pdf", "pdf");
812 }
813 
814 void PndSttCellTrackFinderAnalysisTask::TestRecoQualityFirstStep() {
815 
816  if (fVerbose > 0)
817  cout << "===TestRecoQualityFirstStep===" << endl;
818 
819  FairRootManager* ioman = FairRootManager::Instance();
820 
821  PndMCResult trackCandToMCTrack = fMCMatch->GetMCInfo("FirstTrackCand",
822  "MCTrack");
823  PndMCResult MCTrackToSttHits = fMCMatch->GetMCInfo("MCTrack", "STTHit");
824 
825  // for checking number of links of sttHits
826  PndMCResult myResultHitsToMCTrack = fMCMatch->GetMCInfo("STTHit",
827  "MCTrack");
828  PndMCResult myResultHitToPoint = fMCMatch->GetMCInfo("STTHit", "STTPoint");
829 
830  int index;
831 
832  for (int i = 0;
833  i < fMCTrack->GetEntries() && i < MCTrackToSttHits.GetNEntries();
834  i++) {
835 
836  FairMultiLinkedData STTHitsOfMCTrack =
837  MCTrackToSttHits.GetEntry(i).GetLinksWithType(
838  ioman->GetBranchId("STTHit"));
839 
840  // MCTrack produced more than one STTHits?
841  if (STTHitsOfMCTrack.GetNLinks() > 1) {
842 
843  // increase number of MCTracks in histogram
844  fHistoQualityFirstStep->Fill(0);
845  fFoundMC[i] = false;
846 
847  // search for current MCTrackIndex in trackCands
848  for (int j = 0; j < trackCandToMCTrack.GetNEntries(); j++) {
849 
850  PndTrackCand* myTrackCand = (PndTrackCand*) fFirstTrackCand->At(
851  j);
852  FairMultiLinkedData STTHitsOfTrackCand =
853  myTrackCand->GetLinksWithType(
854  ioman->GetBranchId("STTHit"));
855 
856  if (trackCandToMCTrack.GetEntry(j).GetNLinks() == 1
857  || GetNumLinksOfHits(STTHitsOfTrackCand,
858  myResultHitsToMCTrack, myResultHitToPoint)
859  == 1) {
860  // hits of trackCand points to one MCTrack
861 
862  if (trackCandToMCTrack.GetEntry(j).GetLink(0).GetIndex()
863  == i) {
864  // found current MCTrackIndex
865  fHistoQualityFirstStep->Fill(1);
866  fFoundMC[i] = true;
867  goto nextMCTrack;
868  }
869 
870  }
871  }
872 
873  // search for current MCTrackIndex in faulty trackCands
874  for (int j = 0; j < trackCandToMCTrack.GetNEntries(); j++) {
875 
876  PndTrackCand* myTrackCand = (PndTrackCand*) fFirstTrackCand->At(
877  j);
878  FairMultiLinkedData STTHitsOfTrackCand =
879  myTrackCand->GetLinksWithType(
880  ioman->GetBranchId("STTHit"));
881 
882  if (STTHitsOfTrackCand.GetNLinks() > 1) {
883  // hits of trackCand points to more than one MCTrack
884 
885  // check all links
886  for (int k = 0;
887  k < trackCandToMCTrack.GetEntry(j).GetNLinks();
888  ++k) {
889  if (trackCandToMCTrack.GetEntry(j).GetLink(k).GetIndex()
890  == i) {
891  // found current MCTrackIndex
892  fHistoQualityFirstStep->Fill(2);
893  fFoundMC[i] = true;
894  goto nextMCTrack;
895  }
896  }
897 
898  }
899  }
900 
901  // MCTrack was not found
902  fHistoQualityFirstStep->Fill(3);
903 
904  if (fVerbose > 0) {
905  index = STTHitsOfMCTrack.GetLink(0).GetIndex();
906  cout << "!!! MCTrack not found, tube: "
907  << fMapHitIndexToTubeID[index] << endl;
908 
909  }
910 
911  nextMCTrack: ;
912  }
913 
914  }
915 
916  int candTrackIndex, numMCHits, numTrackHits,
917  numSkewedTubesMC, numSkewedTubesCand, hitIndex; //numberOfLinks, //[R.K. 01/2017] unused variable?
918  double amount;
919 
920  for (int i = 0; i < trackCandToMCTrack.GetNEntries(); i++) {
921  // increase number of TrackCands in histogram
922  fHistoQualityFirstStep->Fill(4);
923 
924  PndTrackCand* myTrackCand = (PndTrackCand*) fFirstTrackCand->At(i);
925  FairMultiLinkedData hitLinks = myTrackCand->GetLinksWithType(
926  ioman->GetBranchId("STTHit"));
927 
928  numTrackHits = hitLinks.GetNLinks();
929 
930  if (trackCandToMCTrack.GetEntry(i).GetNLinks() == 1
931  || GetNumLinksOfHits(hitLinks, myResultHitsToMCTrack,
932  myResultHitToPoint) == 1) {
933  // clean trackCand
934 
935  // get MCTrackIndex
936  candTrackIndex =
937  trackCandToMCTrack.GetEntry(i).GetLink(0).GetIndex();
938 
939  // get number of hits of MCTrack -> compare
940 
941  FairMultiLinkedData STTHitsOfMCTrack = MCTrackToSttHits.GetEntry(
942  candTrackIndex).GetLinksWithType(
943  ioman->GetBranchId("STTHit"));
944 
945  numMCHits = STTHitsOfMCTrack.GetNLinks();
946  numSkewedTubesMC = 0;
947  numSkewedTubesCand = 0;
948 
949  // get number of hits of skewed tubes of MCTrack
950  for (int j = 0; j < STTHitsOfMCTrack.GetNLinks(); ++j) {
951 
952  // Get HitIndex -> get Tube -> check if skewed
953  hitIndex = STTHitsOfMCTrack.GetLink(j).GetIndex();
954 
955  PndSttTube* tube = (PndSttTube*) fTubeArray->At(
956  fMapHitIndexToTubeID[hitIndex]);
957  // count skewed tubes in MCTrack
958  if (tube->IsSkew())
959  ++numSkewedTubesMC;
960 
961  }
962 
963  // get number of hits of skewed tubes ob trackCand
964  for (int j = 0; j < hitLinks.GetNLinks(); ++j) {
965 
966  // Get HitIndex -> get Tube -> check if skewed
967  hitIndex = hitLinks.GetLink(j).GetIndex();
968 
969  PndSttTube* tube = (PndSttTube*) fTubeArray->At(
970  fMapHitIndexToTubeID[hitIndex]);
971  // count skewed tubes in MCTrack
972  if (tube->IsSkew())
973  ++numSkewedTubesCand;
974  }
975 
976  /* TrackCand may contain hits of skewedTubes but not necessarily.
977  * It must not contain more hits than MCTrack.*/
978  if ((numMCHits - numSkewedTubesMC)
979  == (numTrackHits - numSkewedTubesCand)) {
980 
981  fHistoQualityFirstStep->Fill(6);
982  if (fVerbose > 0) {
983  index = hitLinks.GetLink(0).GetIndex();
984  cout << "FULLY RECO: " << fMapHitIndexToTubeID[index]
985  << endl;
986  }
987 
988  } else {
989  fHistoQualityFirstStep->Fill(5);
990 
991  amount = 1
992  - (double) (numTrackHits - numSkewedTubesCand)
993  / (numMCHits - numSkewedTubesMC);
994 
995  if (amount < 0.25) {
996  fHistoNumberOfMissingHits1->Fill(0);
997  } else if (amount < 0.5) {
998  fHistoNumberOfMissingHits1->Fill(1);
999  } else if (amount < 0.75) {
1000  fHistoNumberOfMissingHits1->Fill(2);
1001  } else {
1002  fHistoNumberOfMissingHits1->Fill(3);
1003  }
1004  }
1005 
1006  } else {
1007  //faulty trackCand -> ghost
1008  fHistoQualityFirstStep->Fill(7);
1009  }
1010  }
1011 
1012 }
1013 
1014 void PndSttCellTrackFinderAnalysisTask::TestRecoQualityCombi() {
1015 
1016  if (fVerbose > 0)
1017  cout << "=== TestRecoQualityCombi===" << endl;
1018 
1019  FairRootManager* ioman = FairRootManager::Instance();
1020 
1021  PndMCResult trackCandToMCTrack = fMCMatch->GetMCInfo("CombiTrackCand",
1022  "MCTrack");
1023  PndMCResult MCTrackToSttHits = fMCMatch->GetMCInfo("MCTrack", "STTHit");
1024 
1025  // for checking number of links of sttHits
1026  PndMCResult myResultHitsToMCTrack = fMCMatch->GetMCInfo("STTHit",
1027  "MCTrack");
1028  PndMCResult myResultHitToPoint = fMCMatch->GetMCInfo("STTHit", "STTPoint");
1029 
1030  int index;
1031  // for each MCTrack
1032  for (int i = 0;
1033  i < fMCTrack->GetEntries() && i < MCTrackToSttHits.GetNEntries();
1034  i++) {
1035 
1036  // get STTHit-Links
1037  FairMultiLinkedData STTHitsOfMCTrack =
1038  MCTrackToSttHits.GetEntry(i).GetLinksWithType(
1039  ioman->GetBranchId("STTHit"));
1040 
1041  // MCTrack produced more than one STTHits?
1042  if (STTHitsOfMCTrack.GetNLinks() > 1) {
1043 
1044  // increase number of MCTracks in histogram
1045  fHistoQualityCombi->Fill(0);
1046  index = STTHitsOfMCTrack.GetLink(0).GetIndex();
1047 
1048  // fill histo with number of hits per MCTrack
1049  fHistoNumberOfHitsMCTrack->Fill(STTHitsOfMCTrack.GetNLinks());
1050 
1051  // search for current MCTrackIndex in pure trackCands
1052  for (int j = 0; j < trackCandToMCTrack.GetNEntries(); j++) {
1053 
1054  PndTrackCand* myTrackCand = (PndTrackCand*) fCombiTrackCand->At(
1055  j);
1056 
1057  // get links of trackCand to STTHits
1058  FairMultiLinkedData STTHitsOfTrackCand =
1059  myTrackCand->GetLinksWithType(
1060  ioman->GetBranchId("STTHit"));
1061 
1062  // only one link to MCTrack?
1063  if (trackCandToMCTrack.GetEntry(j).GetNLinks() == 1
1064  || GetNumLinksOfHits(STTHitsOfTrackCand,
1065  myResultHitsToMCTrack, myResultHitToPoint)
1066  == 1) {
1067  // hits of trackCand points to one MCTrack
1068 
1069  if (trackCandToMCTrack.GetEntry(j).GetLink(0).GetIndex()
1070  == i) {
1071  // found current MCTrackIndex in pure tracklet
1072  fHistoQualityCombi->Fill(1);
1073 
1074  if (fVerbose > 0) {
1075  if (fFoundMC[i] != true)
1076  cout << "FOUND MORE: tube "
1077  << fMapHitIndexToTubeID[index] << endl;
1078  }
1079  goto nextMCTrack;
1080  }
1081 
1082  }
1083  }
1084 
1085  // search for current MCTrackIndex in faulty trackCands
1086  for (int j = 0; j < trackCandToMCTrack.GetNEntries(); j++) {
1087 
1088  PndTrackCand* myTrackCand = (PndTrackCand*) fCombiTrackCand->At(
1089  j);
1090  FairMultiLinkedData STTHitsOfTrackCand =
1091  myTrackCand->GetLinksWithType(
1092  ioman->GetBranchId("STTHit"));
1093 
1094  if (STTHitsOfTrackCand.GetNLinks() > 1) {
1095  // hits of trackCand points to more than one MCTrack
1096 
1097  // check all links
1098  for (int k = 0;
1099  k < trackCandToMCTrack.GetEntry(j).GetNLinks();
1100  ++k) {
1101  if (trackCandToMCTrack.GetEntry(j).GetLink(k).GetIndex()
1102  == i) {
1103  // found current MCTrackIndex
1104 
1105  fHistoQualityCombi->Fill(2);
1106 
1107  if (fVerbose > 0) {
1108  if (fFoundMC[i] != true)
1109  cout << "FOUND MORE: tube "
1110  << fMapHitIndexToTubeID[index]
1111  << endl;
1112  }
1113 
1114  goto nextMCTrack;
1115  }
1116  }
1117 
1118  }
1119  }
1120 
1121  //MCTrack was not found
1122  fHistoQualityCombi->Fill(3);
1123 
1124  if (fVerbose > 0) {
1125  cout << "!!! MCTrack not found, tube: "
1126  << fMapHitIndexToTubeID[index] << endl;
1127 
1128  if (fFoundMC[i] == true)
1129  cout << "FOUND LESS: tube " << fMapHitIndexToTubeID[index]
1130  << endl;
1131  }
1132 
1133  nextMCTrack: ;
1134  }
1135 
1136  }
1137 
1138  int candTrackIndex, numMCHits, numTrackHits,
1139  numSkewedTubesMC, numSkewedTubesCand, hitIndex;// numberOfLinks, //[R.K. 01/2017] unused variable?
1140  double amount;
1141 
1142  for (int i = 0; i < trackCandToMCTrack.GetNEntries(); i++) {
1143  // increase number of TrackCands in histogram
1144  fHistoQualityCombi->Fill(4);
1145 
1146  PndTrackCand* myTrackCand = (PndTrackCand*) fCombiTrackCand->At(i);
1147  FairMultiLinkedData hitLinks = myTrackCand->GetLinksWithType(
1148  ioman->GetBranchId("STTHit"));
1149 
1150  numTrackHits = hitLinks.GetNLinks();
1151 
1152  if (trackCandToMCTrack.GetEntry(i).GetNLinks() == 1
1153  || GetNumLinksOfHits(hitLinks, myResultHitsToMCTrack,
1154  myResultHitToPoint) == 1) {
1155  // clean trackCand
1156 
1157  // get MCTrackIndex
1158  candTrackIndex =
1159  trackCandToMCTrack.GetEntry(i).GetLink(0).GetIndex();
1160 
1161  // get number of hits of MCTrack -> compare
1162 
1163  FairMultiLinkedData STTHitsOfMCTrack = MCTrackToSttHits.GetEntry(
1164  candTrackIndex).GetLinksWithType(
1165  ioman->GetBranchId("STTHit"));
1166 
1167  numMCHits = STTHitsOfMCTrack.GetNLinks();
1168  numSkewedTubesMC = 0;
1169  numSkewedTubesCand = 0;
1170 
1171  // get number of hits of skewed tubes of MCTrack
1172  for (int j = 0; j < STTHitsOfMCTrack.GetNLinks(); ++j) {
1173 
1174  // Get HitIndex -> get Tube -> check if skewed
1175  hitIndex = STTHitsOfMCTrack.GetLink(j).GetIndex();
1176 
1177  PndSttTube* tube = (PndSttTube*) fTubeArray->At(
1178  fMapHitIndexToTubeID[hitIndex]);
1179  // count skewed tubes in MCTrack
1180  if (tube->IsSkew())
1181  ++numSkewedTubesMC;
1182 
1183  }
1184 
1185  // get number of hits of skewed tubes of trackCand
1186  for (int j = 0; j < hitLinks.GetNLinks(); ++j) {
1187 
1188  // Get HitIndex -> get Tube -> check if skewed
1189  hitIndex = hitLinks.GetLink(j).GetIndex();
1190 
1191  PndSttTube* tube = (PndSttTube*) fTubeArray->At(
1192  fMapHitIndexToTubeID[hitIndex]);
1193  // count skewed tubes in MCTrack
1194  if (tube->IsSkew())
1195  ++numSkewedTubesCand;
1196  }
1197 
1198  /* TrackCand may contain hits of skewedTubes but not necessarily.
1199  * It must not contain more hits than MCTrack.*/
1200  if ((numMCHits - numSkewedTubesMC)
1201  == (numTrackHits - numSkewedTubesCand)) {
1202  fHistoQualityCombi->Fill(6);
1203 
1204  if (fVerbose > 0) {
1205  index = hitLinks.GetLink(0).GetIndex();
1206  cout << "FULLY RECO: " << fMapHitIndexToTubeID[index]
1207  << endl;
1208  }
1209 
1210  } else {
1211  fHistoQualityCombi->Fill(5);
1212 
1213  amount = 1
1214  - (double) (numTrackHits - numSkewedTubesCand)
1215  / (numMCHits - numSkewedTubesMC);
1216 
1217  if (amount < 0.25) {
1218  fHistoNumberOfMissingHits2->Fill(0);
1219  } else if (amount < 0.5) {
1220  fHistoNumberOfMissingHits2->Fill(1);
1221  } else if (amount < 0.75) {
1222  fHistoNumberOfMissingHits2->Fill(2);
1223  } else {
1224  fHistoNumberOfMissingHits2->Fill(3);
1225  }
1226  }
1227 
1228  } else {
1229  //faulty trackCand -> ghost
1230  fHistoQualityCombi->Fill(7);
1231  }
1232  }
1233 
1234 }
1235 
1236 void PndSttCellTrackFinderAnalysisTask::CountMaxHitsFirstStep() {
1237 
1238  if (fVerbose > 1)
1239  cout << "=== CountMaxHitsFirstStep===" << endl;
1240 
1241  FairRootManager* ioman = FairRootManager::Instance();
1242 
1243  PndMCResult trackCandToMCTrack = fMCMatch->GetMCInfo("FirstTrackCand",
1244  "MCTrack");
1245  PndMCResult MCTrackToSttHits = fMCMatch->GetMCInfo("MCTrack", "STTHit");
1246 
1247 // for checking number of links of sttHits
1248  PndMCResult myResultHitsToMCTrack = fMCMatch->GetMCInfo("STTHit",
1249  "MCTrack");
1250  PndMCResult myResultHitToPoint = fMCMatch->GetMCInfo("STTHit", "STTPoint");
1251 
1252  int maxHits;
1253 
1254  for (int i = 0;
1255  i < fMCTrack->GetEntries() && i < MCTrackToSttHits.GetNEntries();
1256  i++) {
1257 
1258  FairMultiLinkedData STTHitsOfMCTrack =
1259  MCTrackToSttHits.GetEntry(i).GetLinksWithType(
1260  ioman->GetBranchId("STTHit"));
1261 
1262 // MCTrack produced STTHits?
1263  if (STTHitsOfMCTrack.GetNLinks() > 0) {
1264 
1265  maxHits = 0;
1266 
1267  // search for appropriate trackCand with max number of hits
1268  for (int j = 0; j < trackCandToMCTrack.GetNEntries(); j++) {
1269 
1270  // get trackCand and links of hits
1271  PndTrackCand* myTrackCand = (PndTrackCand*) fFirstTrackCand->At(
1272  j);
1273  FairMultiLinkedData STTHitsOfTrackCand =
1274  myTrackCand->GetLinksWithType(
1275  ioman->GetBranchId("STTHit"));
1276 
1277  // if pure tracklet, check number of hits
1278  if (trackCandToMCTrack.GetEntry(j).GetNLinks() == 1
1279  || GetNumLinksOfHits(STTHitsOfTrackCand,
1280  myResultHitsToMCTrack, myResultHitToPoint)
1281  == 1) {
1282 
1283  // hits of trackCand points to one MCTrack
1284  if (trackCandToMCTrack.GetEntry(j).GetLink(0).GetIndex()
1285  == i) {
1286  // found current MCTrackIndex
1287 
1288  maxHits = TMath::Max(maxHits,
1289  STTHitsOfTrackCand.GetNLinks());
1290  }
1291 
1292  }
1293  }
1294 
1295  // maxHits could be 0, if no trackCand was found
1296  fHistoMaxHitsOfMCTrack1->Fill(maxHits);
1297  }
1298  }
1299 }
1300 
1301 void PndSttCellTrackFinderAnalysisTask::CountMaxHitsCombi() {
1302 
1303  if (fVerbose > 1)
1304  cout << "=== CountMaxHitsCombi===" << endl;
1305 
1306  FairRootManager* ioman = FairRootManager::Instance();
1307 
1308  PndMCResult trackCandToMCTrack = fMCMatch->GetMCInfo("CombiTrackCand",
1309  "MCTrack");
1310  PndMCResult MCTrackToSttHits = fMCMatch->GetMCInfo("MCTrack", "STTHit");
1311 
1312 // for checking number of links of sttHits
1313  PndMCResult myResultHitsToMCTrack = fMCMatch->GetMCInfo("STTHit",
1314  "MCTrack");
1315  PndMCResult myResultHitToPoint = fMCMatch->GetMCInfo("STTHit", "STTPoint");
1316 
1317  int maxHits;
1318 
1319  for (int i = 0;
1320  i < fMCTrack->GetEntries() && i < MCTrackToSttHits.GetNEntries();
1321  i++) {
1322 
1323  FairMultiLinkedData STTHitsOfMCTrack =
1324  MCTrackToSttHits.GetEntry(i).GetLinksWithType(
1325  ioman->GetBranchId("STTHit"));
1326 
1327 // MCTrack produced STTHits?
1328  if (STTHitsOfMCTrack.GetNLinks() > 0) {
1329 
1330  maxHits = 0;
1331 
1332  // search for appropriate trackCand with max number of hits
1333  for (int j = 0; j < trackCandToMCTrack.GetNEntries(); j++) {
1334 
1335  // get trackCand and links of hits
1336  PndTrackCand* myTrackCand = (PndTrackCand*) fCombiTrackCand->At(
1337  j);
1338  FairMultiLinkedData STTHitsOfTrackCand =
1339  myTrackCand->GetLinksWithType(
1340  ioman->GetBranchId("STTHit"));
1341 
1342  // if pure tracklet, check number of hits
1343  if (trackCandToMCTrack.GetEntry(j).GetNLinks() == 1
1344  || GetNumLinksOfHits(STTHitsOfTrackCand,
1345  myResultHitsToMCTrack, myResultHitToPoint)
1346  == 1) {
1347 
1348  // hits of trackCand points to one MCTrack
1349  if (trackCandToMCTrack.GetEntry(j).GetLink(0).GetIndex()
1350  == i) {
1351  // found current MCTrackIndex
1352 
1353  maxHits = TMath::Max(maxHits,
1354  STTHitsOfTrackCand.GetNLinks());
1355  }
1356 
1357  }
1358  }
1359 
1360  // maxHits could be 0 if no trackCand was found
1361  fHistoMaxHitsOfMCTrack2->Fill(maxHits);
1362  }
1363  }
1364 }
1365 
1366 void PndSttCellTrackFinderAnalysisTask::FinishEvent() {
1367 
1368  fMapTubeIDToHits.clear();
1369  fMapHitIndexToTubeID.clear();
1370  fFoundMC.clear();
1371 
1372 }
1373 
1374 void PndSttCellTrackFinderAnalysisTask::Finish() {
1375  fHistoNumberOfLinks1->Write();
1376  fHistoNumberOfTracklets1->Write();
1377 
1378  fHistoNumberOfAssignedHits1->Write();
1379  fHistoNumberOfHits1->Write();
1380 
1381  fHistoNumberOfHits2->Write();
1382  fHistoNumberOfLinks2->Write();
1383  fHistoNumberOfTracklets2->Write();
1384 
1385  fHistoNumberOfErrors1->Write();
1386  fHistoNumberOfErrors2->Write();
1387 
1388  fHistoNumberOfHitsMCTrack->Write();
1389  fHistoMaxHitsOfMCTrack1->Write();
1390  fHistoMaxHitsOfMCTrack2->Write();
1391 
1392  fHistoQualityFirstStep->Write();
1393  fHistoQualityCombi->Write();
1394 
1395  fHistoNumberOfMissingHits1->Write();
1396  fHistoNumberOfMissingHits2->Write();
1397 
1398  if (fVerbose > 0) {
1399 
1400  std::cout << "fHistoQualityFirstStep: #Possible Tracks "
1401  << fHistoQualityFirstStep->GetBinContent(1)
1402  << " #Found in clean: "
1403  << fHistoQualityFirstStep->GetBinContent(2)
1404  << " #Found in faulty: "
1405  << fHistoQualityFirstStep->GetBinContent(3) << " #Not Found: "
1406  << fHistoQualityFirstStep->GetBinContent(4) << " #TrackCands: "
1407  << fHistoQualityFirstStep->GetBinContent(5) << " #Spurious: "
1408  << fHistoQualityFirstStep->GetBinContent(6) << " #FullyReco: "
1409  << fHistoQualityFirstStep->GetBinContent(7) << " #Ghosts: "
1410  << fHistoQualityFirstStep->GetBinContent(8) << std::endl;
1411 
1412  std::cout << "fHistoQualityCombi: #Possible Tracks "
1413  << fHistoQualityCombi->GetBinContent(1) << " #Found in clean: "
1414  << fHistoQualityCombi->GetBinContent(2) << " #Found in faulty: "
1415  << fHistoQualityCombi->GetBinContent(3) << " #Not Found: "
1416  << fHistoQualityCombi->GetBinContent(4) << " #TrackCands: "
1417  << fHistoQualityCombi->GetBinContent(5) << " #Spurious: "
1418  << fHistoQualityCombi->GetBinContent(6) << " #FullyReco: "
1419  << fHistoQualityCombi->GetBinContent(7) << " #Ghosts: "
1420  << fHistoQualityCombi->GetBinContent(8) << std::endl;
1421  }
1422 
1423 }
1424 
pos
TVector3 pos
Definition: GammaSpectraAnalysis.C:43
fVerbose
int fVerbose
Definition: poormantracks.C:24
r
double r
Definition: RiemannTest.C:14
i
Int_t i
Definition: run_full.C:25
CAMath::Sqrt
static T Sqrt(const T &x)
Definition: PndCAMath.h:37
PndSttCellTrackFinderAnalysisTask::GetNumLinksOfHits
int GetNumLinksOfHits(FairMultiLinkedData &hitLinks, PndMCResult &sttHitsToMCTrack, PndMCResult &sttHitToPoint)
Definition: PndSttCellTrackFinderAnalysisTask.cxx:433
PndSttTube::GetPosition
TVector3 GetPosition()
Definition: PndSttTube.cxx:87
rtdb
FairRuntimeDb * rtdb
Definition: hit_dirc.C:66
PndSttHit::GetTubeID
Int_t GetTubeID() const
Definition: PndSttHit.h:75
PndSttMapCreator::FillTubeArray
TClonesArray * FillTubeArray()
Definition: PndSttMapCreator.cxx:128
PndRiemannTrack::r
double r() const
Definition: PndRiemannTrack.cxx:416
CAMath::Max
static T Max(const T &x, const T &y)
Definition: PndCAMath.h:36
x
Double_t x
Definition: createRootFscGeometryFile.C:159
PndRiemannTrack::orig
TVectorD orig() const
Definition: PndRiemannTrack.cxx:405
ClassImp
ClassImp(PndAnaContFact)
hit
PndSdsMCPoint * hit
Definition: anasim.C:70
hits
CbmHit * hits[nHits]
Definition: RiemannTest.C:19
y
Double_t y
Definition: createRootFscGeometryFile.C:159
PndSttTube::IsSkew
bool IsSkew()
Definition: PndSttTube.h:70
Generated on Wed Apr 3 2019 08:44:01 for FairRoot/PandaRoot by   doxygen 1.8.5