FairRoot/PandaRoot
PndMvdEventAnaTask.cxx
Go to the documentation of this file.
1 // -------------------------------------------------------------------------
2 // ----- PndMvdEventAnaTask source file -----
3 // ----- Created 18/07/08 by T.Stockmanns -----
4 // -------------------------------------------------------------------------
5 // libc includes
6 #include <iostream>
7 
8 // Root includes
9 #include "TROOT.h"
10 #include "TClonesArray.h"
11 #include "TVector3.h"
12 
13 
14 // framework includes
15 #include "FairRootManager.h"
16 #include "PndMvdEventAnaTask.h"
17 #include "FairRun.h"
18 #include "FairRuntimeDb.h"
19 #include "FairHit.h"
20 #include "PndMCTrack.h"
21 // PndMvd includes
22 #include "PndSdsMCPoint.h"
23 #include "PndSdsDigiPixel.h"
24 #include "PndSdsDigiStrip.h"
25 #include "PndSdsHit.h"
26 #include "PndSdsCluster.h"
27 #include "GFTrackCand.h"
28 #include "PndRiemannTrack.h"
29 
30 #include <vector>
31 #include <map>
32 
33 
34 // ----- Default constructor -------------------------------------------
35 PndMvdEventAnaTask::PndMvdEventAnaTask() : FairTask("Event Display task for PANDA PndMvd"),
36  fPrintTrack(true), fPrintMCHit(true), fPrintCluster(true), fPrintPixDigis(true), fPrintPixHit(true),
37  fPrintStripCluster(true), fPrintStripDigis(true), fPrintStripHit(true), fPrintTrackMatch(true), fPrintGhosts(true),
38  fNTracks(0), fNPossibleTracks(0), fNCompleteTracks(0), fNPartTracks(0), fNNotFoundPossibleTracks(0), fNNotFoundTracks(0),
39  fNGhostTracks(0), fEventNr(0)
40 {
41 }
42 // -------------------------------------------------------------------------
43 
44 
45 // ----- Destructor ----------------------------------------------------
46 PndMvdEventAnaTask::~PndMvdEventAnaTask()
47 {
48 }
49 
50 // ----- Public method Init --------------------------------------------
51 InitStatus PndMvdEventAnaTask::Init()
52 {
53  // Get RootManager
54  FairRootManager* ioman = FairRootManager::Instance();
55  if ( !ioman){
56  std::cout << "-E- PndMvdEventAnaTask::Init: "<< "RootManager not instantiated!" << std::endl;
57  return kFATAL;
58  }
59 
60  // Get input arrays
61  fMCTracks = (TClonesArray*) ioman->GetObject("MCTrack");
62  if (!fMCTracks){
63  std::cout << "-W- PndMvdEventAnaTask::Init: "<< "No MCTrack" << " array!" << std::endl;
64  return kERROR;
65  }
66 
67  fMCHits = (TClonesArray*) ioman->GetObject("MVDPoint");
68  if ( !fMCHits) {
69  std::cout << "-W- PndMvdEventAnaTask::Init: "<< "No MVDPoint"<<" array!" << std::endl;
70  return kERROR;
71  }
72 
73  fStripDigis = (TClonesArray*) ioman->GetObject("MVDStripDigis");
74  if ( !fStripDigis){
75  std::cout << "-W- PndMvdEventAnaTask::Init: "<< "No MVDStripDigis"<<" array!" << std::endl;
76  }
77 
78  fPixDigis = (TClonesArray*) ioman->GetObject("MVDPixelDigis");
79  if ( !fPixDigis){
80  std::cout << "-W- PndMvdEventAnaTask::Init: "<< "No MVDPixDigis"<<" array!" << std::endl;
81  }
82 
83  fPixReco = (TClonesArray*) ioman->GetObject("MVDHitsPixel");
84  if ( !fPixReco) {
85  std::cout << "-W- PndMvdEventAnaTask::Init: "<< "No MVDHitsPixel"<<" array!" << std::endl;
86  }
87 
88  fStripReco = (TClonesArray*) ioman->GetObject("MVDHitsStrip");
89  if ( !fStripReco) {
90  std::cout << "-W- PndMvdEventAnaTask::Init: "<< "No MVDHitsStrip" << " array!" << std::endl;
91  }
92 
93  fPixCluster = (TClonesArray*) ioman->GetObject("MVDClusterCand");
94  if ( !fPixCluster){
95  std::cout << "-W- PndMvdEventAnaTask::Init: "<< "No MVDClusterCand" <<" array!" << std::endl;
96  }
97 
98  fStripCluster = (TClonesArray*) ioman->GetObject("MVDStripClusterCand");
99  if ( !fStripCluster) {
100  std::cout << "-W- PndMvdEventAnaTask::Init: " << "No MVDStripClusterCand" <<" array!" << std::endl;
101  }
102 
103  fTrackCand = (TClonesArray*) ioman->GetObject("MVDRiemannTrackCand");
104  if ( !fTrackCand) {
105  std::cout << "-W- PndMvdEventAnaTask::Init: " << "No MVDRiemannTrackCand" <<" array!" << std::endl;
106  fTrackCand = (TClonesArray*) ioman->GetObject("MVDIdealTrackCand");
107  if (!fTrackCand)
108  std::cout << "-W- PndMvdEventAnaTask::Init: " << "No MVDIdealTrackCand" <<" array!" << std::endl;
109 
110  }
111 
112  fHTracksPerEvent = new TH1I("HTracksPerEvent", "Tracks per Event", 21, -0.5, 20.5);
113  fHHitsPerTrack = new TH1I("HHitsPerTrack", "Hits per Track", 21,0.5,20.5);
114  fHEnergyPerHit = new TH1D("HEnergyPerHit","Energy per Hit", 200, 0, 10000000);
115 
116  fHPointRes = new TH1D("HPointRes", "Point Resolution", 1000, 0,0.1);
117  fHPointRes->SetLineColor(1);
118  fHPointResS = new TH1D("HPointResS", "Point Resolution", 1000, 0,0.1);
119  fHPointResS->SetLineColor(3);
120  fHPointResD = new TH1D("HPointResD", "Point Resolution", 1000, 0,0.1);
121  fHPointResD->SetLineColor(4);
122  fHPointResM = new TH1D("HPointResM", "Point Resolution", 1000, 0,0.1);
123  fHPointResM->SetLineColor(5);
124 
125  fHDigisPerCluster = new TH1I("HDigisPerCluster","Digis per Cluster", 20, 0, 20);
126  fHEnergyRes = new TH1D("HEnergyRes", "Energy Resolution", 1000, -100000, 100000);
127  fHPRes = new TH1D("HPRes","Momentum Resolution", 1000, -5,5);
128  fHPtRes = new TH1D("HPtRes","Transversal Momentum Resolution", 1000, -5,5);
129 
130  fHPointResStrip = new TH1D("HPointResStrip", "Point Resolution Strips", 1000, 0,0.1);
131  fHPointResStrip->SetLineColor(6);
132  fHPointResSStrip = new TH1D("HPointResSStrip", "Point Resolution Strip", 1000, 0,0.1);
133  fHPointResSStrip->SetLineColor(7);
134  fHPointResDStrip = new TH1D("HPointResDStrip", "Point Resolution Strip", 1000, 0,0.1);
135  fHPointResDStrip->SetLineColor(8);
136  fHPointResMStrip = new TH1D("HPointResMStrip", "Point Resolution Strip", 1000, 0,0.1);
137  fHPointResMStrip->SetLineColor(9);
138 
139  fHDigisPerClusterStrip = new TH1I("HDigisPerClusterStrip","Digis per Cluster Strips", 20, 0, 20);
140  fHDigisPerClusterStrip->SetLineColor(2);
141  fHEnergyResStrip = new TH1D("HEnergyResStrip", "Energy Resolution Strips", 1000, -100000, 100000);
142  fHEnergyResStrip->SetLineColor(2);
143 
144  fHRiemannRes = new TH1I("HRiemannRes", "Quality of the Riemann Trackfinding", 16,-5.5,10.5);
145  fHRiemannFakes = new TH1I("HRiemannFakes", "Wrong assign hits in Track", 16, -5.5,10.5);
146  fHRiemannFakes->SetLineColor(2);
147  fHRiemannTracksPerTrack = new TH1I("HRiemannTracksPerTrack", "Found Tracks per MC Track", 11, -0.5, 10.5);
148  fHRiemannTracksPerTrackAdd = new TH1I("HRiemannTracksPerTrackAdd", "Found Tracks per MC Track", 11, -0.5, 10.5);
149  fHRiemannTracksPerTrackAdd->SetLineColor(2);
150 
151  fHRiemannVertexResolutionX = new TH1D("HRiemannVertexResolutionX","Difference between reco vertex and MC vertex", 1000, -1, 1);
152  fHRiemannVertexResolutionY = new TH1D("HRiemannVertexResolutionY","Difference between reco vertex and MC vertex", 1000, -1, 1);
153  fHRiemannVertexResolutionZ = new TH1D("HRiemannVertexResolutionZ","Difference between reco vertex and MC vertex", 1000, -1, 1);
154 
155  return kSUCCESS;
156 }
157 // -------------------------------------------------------------------------
158 void PndMvdEventAnaTask::SetParContainers()
159 {
160  // Get Base Container
161 // FairRun* ana = FairRun::Instance();
162 // FairRuntimeDb* rtdb=ana->GetRuntimeDb();
163 
164 }
165 
166 
167 // ----- Public method Exec --------------------------------------------
168 void PndMvdEventAnaTask::Exec(Option_t*)
169 {
170  std::map<int, std::vector<int> > mcHitMap; //Track -> MCHits
171  std::map<int, std::vector<int> > trackToTrackCandMap; //Track -> TrackCand
172 
173  std::vector<PndRiemannTrack> riemannTracks;
174  std::vector<int> MCTrackOrderRiemann; //information which riemannTrack belongs to which MCTrack;
175 
176  TVector3 MCPos, RecoPos;
177  double MCEnergy;
178  double TrackP, TrackPt;
179 
180  mcHitMap = AssignHitsToTracks();
181  fGhostCand.clear();
182  fGhostCand.resize(fTrackCand->GetEntriesFast(), 0);
183  fTrackPixHitIdMap.clear();
184  fTrackStripHitIdMap.clear();
185 
186  std::cout << "------------Event " << fEventNr << "-------------" << std::endl;
187  fEventNr++;
188  fHTracksPerEvent->Fill(mcHitMap.size());
189  fNTracks += mcHitMap.size();
190 
191  for (std::map<int, std::vector<int> >::const_iterator kIt = mcHitMap.begin(); kIt!= mcHitMap.end(); kIt++){ //go through all tracks
192  PndMCTrack* myTrack = (PndMCTrack*)(fMCTracks->At(kIt->first));
193  std::vector<int> MChits = kIt->second;
194  TrackPt = myTrack->GetMomentum().Pt();
195  TrackP = myTrack->GetMomentum().Mag();
196  fHHitsPerTrack->Fill(MChits.size());
197 
198  if (fPrintTrack){
199  std::cout << "<<<<<<<<<<< MCTrack >>>>>>>>>> " << std::endl;
200  myTrack->Print(kIt->first);
201  TVector3 startVertex = myTrack->GetStartVertex();
202  std::cout << "Pt: " << TrackPt << " GeV/c; P: " << TrackP << " GeV/c" << std::endl;
203  std::cout << "StartVertex: " << startVertex.X() << " " << startVertex.Y() << " " << startVertex.Z() << std::endl;
204  }
205  for (unsigned int p = 0; p < MChits.size(); p++){ //go through all hits in track
206  PndSdsMCPoint* myPoint = (PndSdsMCPoint*)(fMCHits->At(MChits[p]));
207  if (fPrintMCHit){
208  std::cout << "-------------------------------" << std::endl;
209  myPoint->Print();
210  }
211  MCPos = (myPoint->GetPosition() + myPoint->GetPositionOut())*0.5;
212  MCEnergy = myPoint->GetEnergyLoss()*10E9;
213  fHEnergyPerHit->Fill(MCEnergy);
214 
215  std::vector<int> pixCluster = GetClusters(MChits[p], true); //Indices of all clusters belonging to one hit
216 
217  for (unsigned int clInd = 0; clInd < pixCluster.size(); clInd++) {
218 
219  PndSdsCluster* myCluster = (PndSdsCluster*)(fPixCluster->At(pixCluster[clInd]));
220  std::vector<Int_t> digiInd = myCluster->GetClusterList(); //gets the list of pixel Digis belonging to the cluster
221  int recoHit = GetRecoHit(pixCluster[clInd], true);
222  fHDigisPerCluster->Fill(digiInd.size());
223  if (recoHit > -1)
224  fTrackPixHitIdMap[kIt->first].push_back(recoHit);
225 
226  PrintClusterDigiInfo(pixCluster[clInd], digiInd, true);
227  if (recoHit < 0)
228  std::cout << "-W- No Reco Hit found!" << std::endl;
229  else
230  PrintRecoHitInfo(recoHit, digiInd.size(), MCPos, MCEnergy, true);
231  }
232 
233  std::vector<int> stripCluster = GetClusters(MChits[p], false);
234 
235  int hitCount = 0;
236  int oldRecoHit = -1;
237  for (unsigned int clInd = 0; clInd < stripCluster.size(); clInd++) {
238  int recoHit = -1;
239 
240  PndSdsCluster* myCluster = (PndSdsCluster*)(fStripCluster->At(stripCluster[clInd]));
241  std::vector<Int_t> digiInd = myCluster->GetClusterList(); //gets the list of pixel Digis belonging to the cluster
242  recoHit = GetRecoHit(stripCluster[clInd], false);
243 
244  fHDigisPerClusterStrip->Fill(digiInd.size());
245  std::cout << "RecoHit: " << recoHit << std::endl;
246  if (recoHit > -1){
247  hitCount++;
248  if (oldRecoHit != recoHit){
249  fTrackStripHitIdMap[kIt->first].push_back(recoHit);
250  oldRecoHit = recoHit;
251  }
252  }
253  PrintClusterDigiInfo(stripCluster[clInd], digiInd, false);
254  //std::cout << " hitCount: " << hitCount << " oldReco " << oldRecoHit << " recoHit: " << recoHit << std::endl;
255  if (recoHit > -1){
256  if (hitCount > 1){
257  PrintRecoHitInfo(recoHit, digiInd.size(), MCPos, MCEnergy, false);
258  if ((oldRecoHit > -1) && (oldRecoHit != recoHit))
259  PrintRecoHitInfo(oldRecoHit, digiInd.size(), MCPos, MCEnergy, false);
260  }
261  }
262  }
263  }
264  std::vector<int> pixHits = fTrackPixHitIdMap[kIt->first];
265  std::vector<int> stripHits = fTrackStripHitIdMap[kIt->first];
266 
267  std::cout << std::endl;
268  std::cout << "TrackID " << kIt->first << ": ";
269  for (unsigned int testInd = 0; testInd < pixHits.size(); testInd++)
270  std::cout << " 5/" << pixHits[testInd];
271 
272  for (unsigned int testInd = 0; testInd < stripHits.size(); testInd++)
273  std::cout << " 4/" << stripHits[testInd];
274  std::cout << std::endl;
275 
276  std::cout << "MCHitMap: ";
277  for (unsigned int testInd = 0; testInd < mcHitMap[kIt->first].size(); testInd++)
278  std::cout << mcHitMap[kIt->first].at(testInd) << " ";
279  std::cout << std::endl;
280 
281  std::vector<int> matches;
282  std::vector<int> candidates;
283  GetTrackCandsForMCTrack(pixHits, stripHits, matches, candidates);
284  std::cout << "TrackCands for MCTrack: ";
285  for (unsigned int i = 0; i < candidates.size(); i++) std::cout << candidates[i];
286  std::cout << std::endl;
287  trackToTrackCandMap[kIt->first] = candidates;
288 
289  int TrackMatch = 0;
290  int oldmatches = 0;
291  int highestMatch=-1;
292  for (unsigned int i = 0; i < matches.size(); i++){
293  if (oldmatches < matches[i]){
294  oldmatches = matches[i];
295  highestMatch = i;
296  }
297  GFTrackCand* myCand = (GFTrackCand*)fTrackCand->At(candidates[i]);
298  PndRiemannTrack myRiemannTrack;
299  MCTrackOrderRiemann.push_back(kIt->first);
300  myRiemannTrack.SetVerbose(1);
301  for (UInt_t j = 0; j < myCand->getNHits(); j++){
302  unsigned int detId, hitId;
303  myCand->getHit(j, detId, hitId);
304  PndRiemannHit hit(GetFairHit(detId, hitId));
305  myRiemannTrack.addHit(hit);
306  }
307 
308  myRiemannTrack.refit();
309  myRiemannTrack.szFit();
310  TVectorD n = myRiemannTrack.n();
311  //std::cout << "RiemannNormal: " << n[0] << " " << n[1] << " " << n[2] << std::endl;
312  riemannTracks.push_back(myRiemannTrack);
313 
314  double curv = myCand->getCurv();
315  double dip = myCand->getDip();
316  double pt = 1/curv;
317  pt*=2;
318  pt*=0.00299792;
319  unsigned int detId, hitId;
320  if (matches[i]>2){
321  TrackMatch++;
322  myCand->getHit(0,detId, hitId);
323  if (detId > 0)
324  fHRiemannFakes->Fill(myCand->getNHits()- matches[i]);
325  else fHRiemannFakes->Fill(myCand->getNHits()- matches[i] -1);
326  }
327  if (fPrintTrackMatch){
328  myCand->getHit(0,detId, hitId);
329  std::cout << std::endl;
330  std::cout << "HitMatch: " << matches[i] << std::endl;
331  std::cout << "Hits in Track: " << myCand->getNHits() << " Difference: ";
332  if (detId > 0) std::cout << myCand->getNHits() - matches[i] << std::endl;
333  else std::cout << myCand->getNHits() - matches [i] -1 << std::endl;
334 
335  PrintTrackCand(myCand);
336  std::cout << "Pt: " << pt << " GeV/c " << " P: " << pt/dip << " GeV/c" << std::endl;
337  std::cout << "Pt error: " << TrackPt - pt << " P error: " << TrackP - pt/dip << std::endl;
338  std::cout << "TrackID: " << kIt->first << " mcHitMap.size: " << mcHitMap[kIt->first].size() << " matches: " << matches[i] << std::endl;
339  std::cout << "Missing Hits: " << mcHitMap[kIt->first].size() - matches[i] << std::endl;
340  if (mcHitMap[kIt->first].size() < (myCand->getNHits()-1))
341  fHRiemannRes->Fill(-2); //more hits than expected
342  }
343  }
344 
345  if (mcHitMap[kIt->first].size() > 2){
346  if (mcHitMap[kIt->first].size() > 3){
347  fNPossibleTracks++;
348  if (oldmatches == 0)
349  fNNotFoundPossibleTracks++;
350  }
351  else if (myTrack->GetMotherID() < 0){
352  fNPossibleTracks++;
353  if (oldmatches == 0)
354  fNNotFoundPossibleTracks++;
355  }
356  }
357  if (oldmatches == 0)
358  fNNotFoundTracks++;
359 
360 
361  if (matches.size() > 0){ // there are riemann tracks to a MC Track
362  GFTrackCand* myCand = (GFTrackCand*)fTrackCand->At(candidates[highestMatch]);
363  double curv = myCand->getCurv();
364  double dip = myCand->getDip();
365  double pt = 1/curv;
366  pt*=2;
367  pt*=0.00299792;
368  fHPtRes->Fill(TrackPt - pt);
369  fHPRes->Fill(TrackP - pt/dip);
370  fHRiemannRes->Fill(mcHitMap[kIt->first].size() - matches[highestMatch]);
371  if (mcHitMap[kIt->first].size() - matches[highestMatch] == 0)
372  fNCompleteTracks++;
373  else fNPartTracks++;
374  }
375  std::cout << "Found TracksPerTrack: " << TrackMatch;
376  if (mcHitMap[kIt->first].size() > 2) std::cout << " 3Hits+" << std::endl;
377  if (mcHitMap[kIt->first].size() > 2) fHRiemannTracksPerTrackAdd->Fill(TrackMatch);
378  fHRiemannTracksPerTrack->Fill(TrackMatch);
379  }
380  if (riemannTracks.size() > 1){
381  std::cout << riemannTracks.size() << " Riemann Tracks used!" << std::endl;
382  for (unsigned int i = 0; i < riemannTracks.size() - 1; i++){
383  for (unsigned int j = i + 1; j < riemannTracks.size(); j++) {
384  TVector3 p1, p2;
385  //riemannTracks[i].refit();
386  //riemannTracks[i].szFit();
387  //std::cout << "Track: " << i << std::endl;
388  //riemannTracks[i].PrintHits();
389 
390  //riemannTracks[j].refit();
391  //riemannTracks[j].szFit();
392  //std::cout << "Track: " << j << std::endl;
393  //riemannTracks[j].PrintHits();
394 
395  std::cout << "Vertex Test for: " << MCTrackOrderRiemann[i] << " and " << MCTrackOrderRiemann[j] << std::endl;
396  int result = riemannTracks[i].calcIntersection(riemannTracks[j], p1, p2);
397 
398  if (result == 1){
399  std::cout << "Vertex1: " << p1.X() << " " << p1.Y() << " " << p1.Z() << std::endl;
400  std::cout << "Vertex2: " << p2.X() << " " << p2.Y() << " " << p2.Z() << std::endl;
401  PndMCTrack* myMCTrack = (PndMCTrack*)(fMCTracks->At(MCTrackOrderRiemann[i]));
402  TVector3 MCVertex = myMCTrack->GetStartVertex();
403  std::cout << "MCVertex " << MCTrackOrderRiemann[i] << ": "<< MCVertex.X() << " " << MCVertex.Y() << " " << MCVertex.Z() << std::endl;
404  MCVertex -= (p1+p2)*0.5;
405  std::cout << "Difference: " << MCVertex.X() << " " << MCVertex.Y() << " " << MCVertex.Z() << std::endl;
406  fHRiemannVertexResolutionX->Fill(MCVertex.X());
407  fHRiemannVertexResolutionY->Fill(MCVertex.Y());
408  fHRiemannVertexResolutionZ->Fill(MCVertex.Z());
409 
410 
411  }
412  }
413  }
414  }
415  if (fPrintGhosts){
416  std::cout << std::endl;
417  std::cout << "----Ghost TrackCandidates:-------" << std::endl;
418  for (unsigned int tcIndex = 0; tcIndex < fGhostCand.size(); tcIndex++){
419  if (fGhostCand[tcIndex] < 3){
420  GFTrackCand* myTrackCand = (GFTrackCand*)fTrackCand->At(tcIndex);
421  PrintTrackCand(myTrackCand);
422  fHRiemannRes->Fill(-1); //Ghosts
423  fNGhostTracks++;
424  }
425  }
426  std::cout << std::endl;
427  }
428 }
429 
430 void PndMvdEventAnaTask::Finish()
431 {
432  std::cout << std::endl;
433  std::cout << "--------- Summary Track Fitting -----------" << std::endl;
434  std::cout << "Number of Tracks: " << fNTracks << std::endl;
435  std::cout << "Number of Tracks with more than 3 hits: " << fNPossibleTracks << std::endl;
436  std::cout << (double)fNCompleteTracks / fNTracks * 100<< " % found completely" << std::endl;
437  std::cout << (double)fNCompleteTracks / fNPossibleTracks * 100<< " % of possible tracks found completely" << std::endl;
438  std::cout << (double)fNPartTracks / fNTracks * 100<< " % found partially" << std::endl;
439  std::cout << (double)fNPartTracks / fNPossibleTracks * 100<< " % of possible tracks found partially" << std::endl;
440  std::cout << (double)fNNotFoundPossibleTracks / fNPossibleTracks * 100 << " % of possible tracks not found" << std::endl;
441  std::cout << (double)fNNotFoundTracks / fNTracks * 100<< " % of Tracks not found in total" << std::endl;
442  std::cout << (double)fNGhostTracks / fNTracks * 100<< " % of Ghost Tracks" << std::endl;
443 }
444 
445 
446 FairHit* PndMvdEventAnaTask::GetFairHit(Int_t detId, Int_t hitId)
447 {
448  FairHit* p;
449 
450  if (detId == 5){
451  return (FairHit *) fPixReco->At(hitId);
452  }
453  else if (detId == 4){
454  return p=(FairHit *) fStripReco->At(hitId);
455  }
456 
457  std::cout << "-E- FairRiemannTrackCandDraw::GetFairHit : Unknown Detector with ID: " << detId << std::endl;
458  return 0;
459 }
460 
461 // -------------------------------------------------------------------------
462 std::vector<int> PndMvdEventAnaTask::GetClusters(int MCHit, bool pixel)
463 {
464  TClonesArray* cluster;
465 
466  std::vector<int> result;
467 
468  if (pixel == true){
469  cluster = fPixCluster;
470  }
471  else{
472  cluster = fStripCluster;
473  }
474 
475  if (cluster != 0) {
476  for (int clIndex = 0; clIndex < cluster->GetEntriesFast(); clIndex++){ //go through all pixel clusters
477  PndSdsCluster* myCluster = (PndSdsCluster*)(cluster->At(clIndex));
478  if (MCHitBelongsToCluster(MCHit, myCluster, pixel)){
479  result.push_back(clIndex);
480  }
481  }
482  }
483  return result;
484 }
485 
486 int PndMvdEventAnaTask::GetRecoHit(int clIndex, bool pixel) const
487 {
488  TClonesArray* reco;
489 
490  if (pixel == true){
491  reco = fPixReco;
492  }
493  else{
494  reco = fStripReco;
495  }
496  if (reco != 0){
497  for (int hitIndex = 0; hitIndex < reco->GetEntriesFast(); hitIndex++){
498  PndSdsHit* myHit = (PndSdsHit*)reco->At(hitIndex);
499  if ((myHit->GetClusterIndex() == clIndex) || myHit->GetBotIndex() == clIndex){ //test if RecoHit belongs to cluster
500  return hitIndex;
501  }
502  }
503  }
504  return -1;
505 }
506 // -------------------------------------------------------------------------
507 void PndMvdEventAnaTask::PrintClusterDigiInfo(int clIndex,
508  std::vector<Int_t> digiInd, bool pixel)
509 {
510  TClonesArray* digis;
511 
512  if (pixel){
513  digis = fPixDigis;
514  }
515  else{
516  digis = fStripDigis;
517  }
518 
519  if (digis == 0)
520  fPrintPixDigis = false;
521 
522  if (fPrintCluster)
523  std::cout << "Hit belongs to cluster " << clIndex << std::endl;
524  if (fPrintPixDigis) {
525  for (unsigned int s = 0; s < digiInd.size(); s++) {
526  PndSdsDigi* myDigi = (PndSdsDigi*)(digis->At(digiInd[s]));
527  if (pixel) ((PndSdsDigiPixel*)myDigi)->Print();
528  else ((PndSdsDigiStrip*)myDigi)->Print();
529  }
530  }
531 }
532 
533 void PndMvdEventAnaTask::PrintRecoHitInfo(int hitInd, int digiSize, TVector3 MCPos, double MCEnergy, bool pixel) const
534 {
535  TClonesArray* reco;
536  TH1* hPointRes;
537  TH1* hPointResS;
538  TH1* hPointResD;
539  TH1* hPointResM;
540  TH1* hEnergyRes;
541 
542  if (pixel){
543  reco = fPixReco;
544 
545  hPointRes = fHPointRes;
546  hPointResS = fHPointResS;
547  hPointResD = fHPointResD;
548  hPointResM = fHPointResM;
549  hEnergyRes = fHEnergyRes;
550  }
551  else{
552  reco = fStripReco;
553 
554  hPointRes = fHPointResStrip;
555  hPointResS = fHPointResSStrip;
556  hPointResD = fHPointResDStrip;
557  hPointResM = fHPointResMStrip;
558  hEnergyRes = fHEnergyResStrip;
559  }
560 
561  PndSdsHit* myHit = (PndSdsHit*)reco->At(hitInd);
562  if (fPrintPixHit)
563  {
564  if (pixel) std::cout << "PixHit: ";
565  else std::cout << "StripHit: ";
566  std::cout << hitInd << std::endl; //write out RecoHit
567  myHit->Print();
568  }
569  TVector3 RecoPos = myHit->GetPosition();
570  double RecoEnergy = myHit->GetEloss()*10E9;
571  TVector3 result = MCPos - RecoPos;
572  ((TH1D*)hPointRes)->Fill(result.Mag());
573  if (digiSize == 1)
574  hPointResS->Fill(result.Mag());
575  if (digiSize == 2)
576  hPointResD->Fill(result.Mag());
577  if (digiSize > 2)
578  hPointResM->Fill(result.Mag());
579  ((TH1D*)hEnergyRes)->Fill(MCEnergy-RecoEnergy);
580 }
581 
582 
583 bool PndMvdEventAnaTask::MCHitBelongsToCluster(int HitIndex, PndSdsCluster* cluster, bool pixCluster)
584 {
585  bool result = false;
586  std::vector<Int_t> clusterList = cluster->GetClusterList();
587  for (unsigned int i = 0; i < clusterList.size() && result == false; i++){
588  PndSdsDigi* myDigi;
589  if (pixCluster)
590  myDigi = (PndSdsDigi*)(fPixDigis->At(clusterList[i]));
591  else
592  myDigi = (PndSdsDigi*)(fStripDigis->At(clusterList[i]));
593  for (int j = 0; j < myDigi->GetNIndices(); j++)
594  if (myDigi->GetIndex(j) == HitIndex)
595  result = true;
596  }
597  return result;
598 }
599 // -------------------------------------------------------------------------
600 
601 void PndMvdEventAnaTask::GetTrackCandsForMCTrack(std::vector<int> pixHitId, std::vector<int> stripHitId,
602  std::vector<int>& matches, std::vector<int>& result)
603 {
604  unsigned int detId, hitId;
605  int hitMatch = 0;
606 
607  if (fTrackCand != 0){
608  for (int i = 0; i < fTrackCand->GetEntriesFast(); i++){
609  GFTrackCand* myTrackCand = (GFTrackCand*)fTrackCand->At(i);
610  hitMatch = 0;
611  for (unsigned int j = 0; j < myTrackCand->getNHits(); j++){
612  myTrackCand->getHit(j, detId, hitId);
613  //std::cout << "DetId: " << detId << " HitId: " << hitId << std::endl;
614  if (detId == 5){
615  for (unsigned int k = 0; k < pixHitId.size(); k++){
616  //std::cout << "PixHitId: " << pixHitId[k] << std::endl;
617  if ((int)hitId == pixHitId[k])
618  hitMatch++;
619  }
620  }
621  else if (detId == 4){
622  for (unsigned int k = 0; k < stripHitId.size(); k++){
623  //std::cout << "StripHitId: " << stripHitId[k] << std::endl;
624  if ((int)hitId == stripHitId[k])
625  hitMatch++;
626  }
627  }
628  }
629  //std::cout << "internal HitMatch: " << hitMatch << std::endl;
630  if (hitMatch > 2){
631  matches.push_back(hitMatch);
632  result.push_back(i);
633  }
634  if (fGhostCand[i] < hitMatch)
635  fGhostCand[i] = hitMatch;
636  }
637  }
638 }
639 
640 std::map<int, std::vector<int> > PndMvdEventAnaTask::AssignHitsToTracks()
641 {
642  std::map<int, std::vector<int> > result;
643  for (int i = 0; i < fMCHits->GetEntriesFast(); i++){ //get all MC Hits
644  PndSdsMCPoint* myPoint = (PndSdsMCPoint*)(fMCHits->At(i)); //sort MCHits with Tracks
645  //PndMCTrack* myTrack = (PndMCTrack*)(fMCTracks->At(myPoint->GetTrackID()));
646  result[myPoint->GetTrackID()].push_back(i);
647 
648  }
649  return result;
650 }
651 
652 void PndMvdEventAnaTask::PrintTrackCand(GFTrackCand* cand) const
653 {
654  unsigned int det, hit;
655  //std::cout << "TrackCand: " << cand->getCurv() << " curv, " << cand->getDip() << " dip, " << (int)(cand->inverted()) << " inverted." << "\n";
656  for (unsigned int i = 0; i < cand->getNHits(); i++){
657  cand->getHit(i, det, hit);
658  std::cout << det << "/" << hit << " ";
659  }
660  std::cout << "\n";
661 }
662 
663 ClassImp(PndMvdEventAnaTask);
PndMvdEventAnaTask::SetParContainers
virtual void SetParContainers()
Definition: PndMvdEventAnaTask.cxx:158
p
Double_t p
Definition: anasim.C:58
PndSdsCluster::GetClusterList
std::vector< Int_t > GetClusterList() const
Definition: PndSdsCluster.h:38
PndMvdEventAnaTask::fStripDigis
TClonesArray * fStripDigis
Definition: PndMvdEventAnaTask.h:87
PndSdsHit::Print
virtual void Print(const Option_t *opt=0) const
Definition: PndSdsHit.cxx:66
PndSdsDigi
Base class for Digi information.
Definition: PndSdsDigi.h:29
PndMvdEventAnaTask::PrintRecoHitInfo
void PrintRecoHitInfo(int hitInd, int digiSize, TVector3 MCPos, double MCEnergy, bool pixel) const
Definition: PndMvdEventAnaTask.cxx:533
i
Int_t i
Definition: run_full.C:25
PndSdsHit::GetPosition
TVector3 GetPosition() const
Definition: PndSdsHit.h:93
PndSdsCluster
Class to store the Digis which belong to one cluster This class holds the information which Digi belo...
Definition: PndSdsCluster.h:19
GFTrackCand::getDip
double getDip() const
Definition: GFTrackCand.h:115
map
PndTransMap * map
Definition: sim_emc_apd.C:99
PndMvdEventAnaTask::GetFairHit
FairHit * GetFairHit(Int_t detId, Int_t hitId)
Definition: PndMvdEventAnaTask.cxx:446
PndSdsDigi::GetIndex
Int_t GetIndex(int i=0) const
Definition: PndSdsDigi.h:63
PndMvdEventAnaTask::AssignHitsToTracks
std::map< int, std::vector< int > > AssignHitsToTracks()
Definition: PndMvdEventAnaTask.cxx:640
s
TLorentzVector s
Definition: Pnd2DStar.C:50
PndSdsDigiStrip
Class for digitised strip hits.
Definition: PndSdsDigiStrip.h:29
GFTrackCand::getNHits
unsigned int getNHits() const
Definition: GFTrackCand.h:113
n
int n
Definition: HypStatDecay.cxx:21
PndMCTrack::GetMomentum
TVector3 GetMomentum() const
Definition: PndMCTrack.h:78
PndSdsMCPoint::GetPositionOut
TVector3 GetPositionOut() const
Definition: PndSdsMCPoint.h:91
PndSdsHit::GetEloss
Double_t GetEloss() const
Definition: PndSdsHit.h:97
PndSdsDigi::Print
virtual void Print()
Definition: PndSdsDigi.h:87
GFTrackCand::getCurv
double getCurv() const
Definition: GFTrackCand.h:114
PndMvdEventAnaTask::PrintTrackCand
void PrintTrackCand(GFTrackCand *cand) const
Definition: PndMvdEventAnaTask.cxx:652
pt
TString pt(TString pts, TString exts="px py pz")
Definition: invexp.C:133
PndRiemannTrack::n
const TVectorD & n() const
Definition: PndRiemannTrack.h:56
PndMvdEventAnaTask::Init
virtual InitStatus Init()
Definition: PndMvdEventAnaTask.cxx:51
PndMvdEventAnaTask::fGhostCand
std::vector< int > fGhostCand
Definition: PndMvdEventAnaTask.h:94
PndMvdEventAnaTask::fTrackPixHitIdMap
std::map< int, std::vector< int > > fTrackPixHitIdMap
Definition: PndMvdEventAnaTask.h:147
PndMvdEventAnaTask::fStripCluster
TClonesArray * fStripCluster
Definition: PndMvdEventAnaTask.h:91
PndMvdEventAnaTask::fStripReco
TClonesArray * fStripReco
Definition: PndMvdEventAnaTask.h:89
GFTrackCand::getHit
void getHit(unsigned int i, unsigned int &detId, unsigned int &hitId) const
Get detector ID and cluster index (hitId) for hit number i.
Definition: GFTrackCand.h:84
PndMvdEventAnaTask::MCHitBelongsToCluster
bool MCHitBelongsToCluster(int HitIndex, PndSdsCluster *cluster, bool pixCluster)
Definition: PndMvdEventAnaTask.cxx:583
PndMvdEventAnaTask::fTrackStripHitIdMap
std::map< int, std::vector< int > > fTrackStripHitIdMap
Definition: PndMvdEventAnaTask.h:148
PndMvdEventAnaTask::Exec
virtual void Exec(Option_t *opt)
Definition: PndMvdEventAnaTask.cxx:168
PndRiemannTrack::refit
void refit(bool withErrorCalc=true)
Definition: PndRiemannTrack.cxx:172
PndMCTrack::Print
void Print(Int_t iTrack=0) const
Definition: PndMCTrack.cxx:95
GFTrackCand
Track candidate – a list of cluster indices.
Definition: GFTrackCand.h:55
PndSdsHit::GetBotIndex
Int_t GetBotIndex() const
Definition: PndSdsHit.h:96
PndMvdEventAnaTask::PrintClusterDigiInfo
void PrintClusterDigiInfo(int clIndex, std::vector< Int_t > digiInd, bool pixel)
Definition: PndMvdEventAnaTask.cxx:507
PndRiemannTrack::SetVerbose
void SetVerbose(int i)
Definition: PndRiemannTrack.h:122
PndMvdEventAnaTask::fTrackCand
TClonesArray * fTrackCand
Definition: PndMvdEventAnaTask.h:92
PndMvdEventAnaTask::GetClusters
std::vector< int > GetClusters(int MCHit, bool pixel)
Definition: PndMvdEventAnaTask.cxx:462
p2
TPad * p2
Definition: hist-t7.C:117
PndSdsMCPoint::GetPosition
TVector3 GetPosition() const
Definition: PndSdsMCPoint.h:90
reco
int reco()
Definition: detectors/disc/reco.C:1
PndMvdEventAnaTask::GetTrackCandsForMCTrack
void GetTrackCandsForMCTrack(std::vector< int > pixHitId, std::vector< int > stripHitId, std::vector< int > &matches, std::vector< int > &result)
Definition: PndMvdEventAnaTask.cxx:601
ClassImp
ClassImp(PndAnaContFact)
p1
TPad * p1
Definition: hist-t7.C:116
PndRiemannTrack::szFit
void szFit(bool withErrorCalc=true)
Definition: PndRiemannTrack.cxx:441
PndRiemannTrack::addHit
void addHit(PndRiemannHit &hit)
Definition: PndRiemannTrack.cxx:102
PndSdsDigiPixel
Data class to store the digi output of a pixel module.
Definition: PndSdsDigiPixel.h:25
PndSdsHit::GetClusterIndex
Int_t GetClusterIndex() const
Definition: PndSdsHit.h:94
hit
PndSdsMCPoint * hit
Definition: anasim.C:70
PndSdsDigi::GetNIndices
Int_t GetNIndices() const
Definition: PndSdsDigi.h:64
PndMCTrack::GetStartVertex
TVector3 GetStartVertex() const
Definition: PndMCTrack.h:76
PndMvdEventAnaTask::fPixReco
TClonesArray * fPixReco
Definition: PndMvdEventAnaTask.h:88
PndMCTrack::GetMotherID
Int_t GetMotherID() const
Definition: PndMCTrack.h:74
PndMvdEventAnaTask::GetRecoHit
int GetRecoHit(int clIndex, bool pixel) const
Definition: PndMvdEventAnaTask.cxx:486
PndMvdEventAnaTask::fPixCluster
TClonesArray * fPixCluster
Definition: PndMvdEventAnaTask.h:90
PndSdsMCPoint::Print
virtual void Print(const Option_t *opt=0) const
Definition: PndSdsMCPoint.cxx:52
PndMvdEventAnaTask::fMCTracks
TClonesArray * fMCTracks
Definition: PndMvdEventAnaTask.h:84
PndMvdEventAnaTask::fMCHits
TClonesArray * fMCHits
Definition: PndMvdEventAnaTask.h:85
PndMvdEventAnaTask::fPixDigis
TClonesArray * fPixDigis
Definition: PndMvdEventAnaTask.h:86
Generated on Wed Apr 3 2019 08:43:59 for FairRoot/PandaRoot by   doxygen 1.8.5