FairRoot/PandaRoot
PndSttTrackFitterQATask.cxx
Go to the documentation of this file.
1 // -------------------------------------------------------------------------
2 // ----- PndSttTrackFitterQATask source file -----
3 // -------------------------------------------------------------------------
4 
5 #include "PndSttTrackFitterQATask.h"
6 
7 #include "PndSttHit.h"
8 #include "PndSttHitInfo.h"
9 #include "PndSttPoint.h"
10 //#include "PndSttSingleStraw.h"
11 #include "PndSttMapCreator.h"
12 #include "PndSttTube.h"
13 #include "PndSttTrack.h"
14 #include "PndTrackCand.h"
15 #include "PndTrack.h"
16 #include "PndTrackCandHit.h"
17 #include "PndSttHelixHit.h"
18 #include "PndMCTrack.h"
19 #include "FairRootManager.h"
20 #include "FairRunAna.h"
21 #include "FairRuntimeDb.h"
22 
23 #include "TClonesArray.h"
24 #include "TRandom.h"
25 #include "TH1F.h"
26 #include <iostream>
27 #include <cmath>
28 
29 using std::cout;
30 using std::endl;
31 using std::sqrt;
32 
33 // ----- Default constructor -------------------------------------------
34 PndSttTrackFitterQATask::PndSttTrackFitterQATask() :
35  FairTask("")
36 {
37 
38 }
39 // -------------------------------------------------------------------------
40 
41 PndSttTrackFitterQATask::PndSttTrackFitterQATask(Int_t verbose) :
42  FairTask("")
43 {
44  fVerbose = verbose;
45 }
46 
47 
48 // ----- Destructor ----------------------------------------------------
49 PndSttTrackFitterQATask::~PndSttTrackFitterQATask()
50 {
51 }
52 // -------------------------------------------------------------------------
53 
54 
55 
56 // ----- Public method Init --------------------------------------------
57 InitStatus PndSttTrackFitterQATask::Init()
58 {
59 
60  hptfit = new TH1F("hptfit", "pt fit: reco - mc", 100, -0.3, 0.3);
61  hplfit = new TH1F("hplfit", "pl fit: reco - mc", 100, -0.4, 0.4);
62  hptotfit = new TH1F("hptotfit", "ptot fit: reco - mc", 100, -0.3, 0.3);
63 
64  hptfit_perc = new TH1F("hptfit_perc", "pt fit: (reco - mc) / mc", 100, -0.3, 0.3);
65  hplfit_perc = new TH1F("hplfit_perc", "pl fit: (reco - mc) / mc", 100, -0.4, 0.4);
66  hptotfit_perc = new TH1F("hptotfit_perc", "ptot fit: (reco - mc) / mc", 100, -0.3, 0.3);
67 
68  hDist = new TH1F("hDist", "Dist fit: reco - mc", 100, -1., 1.);
69  hRad = new TH1F("hRad", "Rad fit: reco - mc", 100, -15., 15.);
70  hPhi = new TH1F("hPhi", "Phi fit: reco - mc", 100, -0.1, 0.1);
71  hTanL = new TH1F("hTanL", "TanL fit: reco - mc", 100, -0.2, 0.2);
72  hZ = new TH1F("hZ", "Z fit: reco - mc", 100, -0.3, 0.3);
73 
74  hQ = new TH1F("hQ", "Q fit: mc + reco", 6, -3, 3);
75 
76  hptfound = new TH1F("hptfound", "pt found: reco - mc", 100, -0.3, 0.3);
77  hplfound = new TH1F("hplfound", "pl found: reco - mc", 100, -0.4, 0.4);
78  hptotfound = new TH1F("hptotfound", "ptot found: reco - mc", 100, -0.3, 0.3);
79 
80  hptfound_perc = new TH1F("hptfound_perc", "pt found: (reco - mc) / mc", 100, -0.3, 0.3);
81  hplfound_perc = new TH1F("hplfound_perc", "pl found: (reco - mc) / mc", 100, -0.4, 0.4);
82  hptotfound_perc = new TH1F("hptotfound_perc", "ptot found: (reco - mc) / mc", 100, -0.3, 0.3);
83 
84  hpxfit = new TH1F("hpxfit", "px fit: reco - mc", 100, -0.3, 0.3);
85  hpyfit = new TH1F("hpyfit", "py fit: reco - mc", 100, -0.3, 0.3);
86  hpzfit = new TH1F("hpzfit", "pz fit: reco - mc", 100, -0.4, 0.4);
87 
88  hresx = new TH1F("hresx", "x: reco - mc", 100, -0.2, 0.2);
89  hresy = new TH1F("hresy", "y: reco - mc", 100, -0.2, 0.2);
90  hresz = new TH1F("hresz", "z: reco - mc", 100, -3, 3);
91 
92  hx = new TH1F("hx", "x no skewed: reco - mc", 100, -0.2, 0.2);
93  hy = new TH1F("hy", "y no skewed: reco - mc", 100, -0.2, 0.2);
94  hz = new TH1F("hz", "z no skewed: reco - mc", 100, -3, 3);
95 
96  hxs = new TH1F("hxs", "x skewed: reco - mc", 100, -0.2, 0.2);
97  hys = new TH1F("hys", "y skewed: reco - mc", 100, -0.2, 0.2);
98  hzs = new TH1F("hzs", "z skewed: reco - mc", 100, -3, 3);
99 
100  // Get RootManager
101  FairRootManager* ioman = FairRootManager::Instance();
102  if ( ! ioman ) {
103  cout << "-E- PndSttHelixHitProducer::Init: "
104  << "RootManager not instantiated!" << endl;
105  return kFATAL;
106  }
107 
108  // Get MCTrack array
109  fMCTrackArray = (TClonesArray*) ioman->GetObject("MCTrack");
110  if ( ! fMCTrackArray)
111  {
112  cout << "-E- PndSttTrackFitterQATask::Init: No MCTrack array!"
113  << endl;
114  return kERROR;
115  }
116 
117  // Get SttTrack array
118  fTrackArray = (TClonesArray*) ioman->GetObject("STTTrack");
119  if ( ! fTrackArray)
120  {
121  cout << "-E- PndSttTrackFitterQATask::Init: No SttTrack array!"
122  << endl;
123  return kERROR;
124  }
125 
126  // Get SttFoundTrack array
127  fFoundTrackArray = (TClonesArray*) ioman->GetObject("STTFoundTrack");
128  if ( ! fFoundTrackArray)
129  {
130  cout << "-E- PndSttTrackFitterQATask::Init: No SttFoundTrack array!"
131  << endl;
132  return kERROR;
133  }
134 
135 
136  // Get SttTrackCand array
137  fTrackCandArray = (TClonesArray*) ioman->GetObject("STTTrackCand");
138  if ( ! fTrackCandArray)
139  {
140  cout << "-E- PndSttTrackFitterQATask::Init: No SttTrack Cand array!"
141  << endl;
142  return kERROR;
143  }
144 
145  // MC points
146  fPointArray = (TClonesArray*) ioman->GetObject("STTPoint");
147  if ( ! fPointArray ) {
148  cout << "-W- PndSttTrackFitterQATask::Init: "
149  << "No STTPoint array!" << endl;
150  return kERROR;
151  }
152 
153  // hits
154  fHitArray = (TClonesArray*) ioman->GetObject("STTHit");
155  if ( ! fHitArray ) {
156  cout << "-W- PndSttTrackFitterQATask::Init: "
157  << "No STTHit array!" << endl;
158  return kERROR;
159  }
160 
161  // helix hits
162  fHelixHitArray = (TClonesArray*) ioman->GetObject("SttHelixHit");
163  if ( ! fHelixHitArray ) {
164  cout << "-W- PndSttTrackFitterQATask::Init: "
165  << "No SttHelixHit array!" << endl;
166  return kERROR;
167  }
168 
169  // CHECK added
170  PndSttMapCreator *mapper = new PndSttMapCreator(fSttParameters);
171  fTubeArray = mapper->FillTubeArray();
172 
173  cout << "-I- PndSttTrackFitterQATask: Intialization successfull" << endl;
174 
175  return kSUCCESS;
176 
177 }
178 // -------------------------------------------------------------------------
179 
180 // CHECK added
181 void PndSttTrackFitterQATask::SetParContainers() {
182  FairRuntimeDb* rtdb = FairRunAna::Instance()->GetRuntimeDb();
183  fSttParameters = (PndGeoSttPar*) rtdb->getContainer("PndGeoSttPar");
184 }
185 
186 // ----- Public method Exec --------------------------------------------
187 void PndSttTrackFitterQATask::Exec(Option_t*)
188 {
189  // Declare some variables
190  PndSttTrack* pTrack = NULL;
191  PndTrackCand* pTrackCand = NULL;
192  PndTrack* pFoundTrack = NULL;
193  PndMCTrack * mcTrack = NULL;
194 
195  if ( ! fTrackArray ) Fatal("Exec", "No fTrackArray");
196 
197  // Loop over fitted tracks
198  Int_t nTracks = fTrackArray->GetEntriesFast();
199 
200  for (Int_t j = 0; j < nTracks; j++) {
201  pTrack = (PndSttTrack *) fTrackArray->At(j);
202  if(!pTrack) continue;
203  Int_t trackCandID = pTrack->GetTrackCandIndex();
204  pTrackCand = (PndTrackCand *) fTrackCandArray->At(trackCandID);
205  if(!pTrackCand) continue;
206 
207  bool foundtrack = false;
208  for(int k = 0; k < fFoundTrackArray->GetEntriesFast(); k++) {
209  pFoundTrack = (PndTrack*) fFoundTrackArray->At(k);
210  if(!pFoundTrack) continue;
211  if(pFoundTrack->GetRefIndex() == trackCandID) {
212  foundtrack = true;
213  break;
214  }
215  }
216  if(foundtrack == false) continue;
217 
218  // ================= momentum residual ========================
219  if(pTrack->GetFlag() < 3) continue; // only prefit-fit-zfit CHECK
220  // -------------------------------------------------------------------
221  // FITTED track
222  // xy
223  Int_t h = -(Int_t) pTrack->GetCharge();
224  Double_t d0 = pTrack->GetDist();
225  Double_t phi0 = pTrack->GetPhi();
226  Double_t Rad = pTrack->GetRad();
227  // z
228  //Double_t z0 = pTrack->GetZ(); //[R.K. 01/2017] unused variable?
229  Double_t zslope = pTrack->GetTanL();
230  // center of curvature of helix
231  TVector2 vec((d0+Rad)*cos(phi0), (d0+Rad)*sin(phi0));
232 
233  Double_t ptran = 0.003 * 2 * Rad;
234  Double_t plong = ptran * zslope;
235  Double_t ptot = sqrt(plong*plong + ptran*ptran);
236 
237  // -------------------------------------------------------------------
238  // FOUND track
239  TVector3 foundMom = pFoundTrack->GetParamFirst().GetMomentum();
240  //Double_t momMag = foundMom.Mag(); //[R.K. 01/2017] unused variable?
241 
242  // -------------------------------------------------------------------
243  // MC track
244  Int_t mcIndex = pTrackCand->getMcTrackId();
245  mcTrack = (PndMCTrack *) fMCTrackArray->At(mcIndex);
246  if(!mcTrack) continue;
247  // if(mcTrack->GetMotherID() != -1) continue;
248  TVector3 mcmom = mcTrack->GetMomentum();
249  TVector3 vertex = mcTrack->GetStartVertex();
250 
251  int mcQ = (int)((TDatabasePDG::Instance())->GetParticle(mcTrack->GetPdgCode())->Charge()/3.);
252  hQ->Fill(- h + mcQ);
253 
254 
255  hptfit->Fill(ptran - mcmom.Perp());
256  hplfit->Fill(plong - mcmom.Z());
257  hptotfit->Fill(ptot - mcmom.Mag());
258 
259  hptfound->Fill(foundMom.Perp() - mcmom.Perp());
260  hplfound->Fill(foundMom.Z() - mcmom.Z());
261  hptotfound->Fill(foundMom.Mag() - mcmom.Mag());
262 
263  hptfit_perc->Fill((ptran - mcmom.Perp())/ mcmom.Perp());
264  hplfit_perc->Fill((plong - mcmom.Z())/ mcmom.Z());
265  hptotfit_perc->Fill((ptot - mcmom.Mag())/ mcmom.Mag());
266 
267  hptfound_perc->Fill((foundMom.Perp() - mcmom.Perp())/ mcmom.Perp());
268  hplfound_perc->Fill((foundMom.Z() - mcmom.Z())/ mcmom.Z());
269  hptotfound->Fill((foundMom.Mag() - mcmom.Mag()) / mcmom.Mag());
270 
271 
272 // cout << "PT : " << mcmom.Perp() << " " << foundMom.Perp() << " " << ptran << endl;
273 // cout << "PL : " << mcmom.Z() << " " << foundMom.Z() << " " << plong << endl;
274 // cout << "PTOT : " << mcmom.Mag() << " " << foundMom.Mag() << " " << ptot << endl;
275 // cout << "--------------------------------------------------" << endl;
276 
277  // === momentum coordinates ===
278  TVector3 fitmom(ptran * TMath::Cos(phi0 - h * TMath::Pi()/2.),
279  ptran * TMath::Sin(phi0 - h * TMath::Pi()/2.),
280  plong);
281 
282  hpxfit->Fill(fitmom.X() - mcmom.X());
283  hpyfit->Fill(fitmom.Y() - mcmom.Y());
284  hpzfit->Fill(fitmom.Z() - mcmom.Z());
285 
286 
287  // ================= parameters residual ========================
288  double mcRad = mcmom.Perp()/0.006;
289 
290  // track from tangent ---------------------
291  //double mc_m1 = mcmom.Y() / mcmom.X(); //[R.K. 01/2017] unused variable?
292  //double mc_q1 = vertex.Y() - vertex.X() * mc_m1; //[R.K. 01/2017] unused variable?
293  //double mc_m2 = -1./mc_m1; //[R.K. 01/2017] unused variable?
294  //double mc_q2 = vertex.Y() - vertex.X() * mc_m2; //[R.K. 01/2017] unused variable?
295 
296  double alpha = TMath::ATan2(mcmom.X(), mcmom.Y());
297  double mcX0, mcY0;
298  if(-h > 0) {
299  mcX0 = vertex.X() + mcRad * TMath::Cos(alpha);
300  mcY0 = vertex.Y() - mcRad * TMath::Sin(alpha);
301  }
302  else {
303  mcX0 = vertex.X() - mcRad * TMath::Cos(alpha);
304  mcY0 = vertex.Y() + mcRad * TMath::Sin(alpha);
305  }
306 
307  Double_t mcDist, mcPhi;
308 
309  mcDist = TMath::Sqrt(mcX0 * mcX0 + mcY0 * mcY0) - mcRad;
310  mcPhi = atan2(mcY0, mcX0);
311 
312  Double_t mcTanL; //[R.K. 01/2017] unused variable?//, mcZ;
313  mcTanL = mcmom.Z()/mcmom.Perp();
314 
315  // mcZ = ??
316 
317  if(pTrack->GetFlag() >= 2) hDist->Fill(d0 - mcDist);
318  if(pTrack->GetFlag() >= 2) hRad->Fill(Rad - mcRad);
319  if(pTrack->GetFlag() >= 2) hPhi->Fill(phi0 - mcPhi);
320  if(pTrack->GetFlag() >= 3) hTanL->Fill(zslope - mcTanL);
321  // if(pTrack->GetFlag() >= 3) hZ->Fill(z0 - mcZ); // ??
322 
323 
324  // ========================= helix hits ==============================
325  Int_t hitcounter = pTrack->GetNofHelixHits();
326  TVector2 point; // point
327 
328  for (Int_t k = 0; k < hitcounter; k++) {
329  Int_t iHHit = pTrack->GetHelixHitIndex(k);
330  PndSttHelixHit *helixhit = (PndSttHelixHit*) fHelixHitArray->At(iHHit);
331  Int_t hitindex = helixhit->GetHitIndex();
332  PndSttHit* hit = (PndSttHit*) fHitArray->At(hitindex);
333 
334  // tubeID CHECK added
335  Int_t tubeID = hit->GetTubeID();
336  PndSttTube *tube = (PndSttTube*) fTubeArray->At(tubeID);
337 
338  PndSttPoint *mcpoint = (PndSttPoint*) fPointArray->At(hit->GetRefIndex());
339 
340  hresx->Fill(helixhit->GetX() - mcpoint->GetX());
341  hresy->Fill(helixhit->GetY() - mcpoint->GetY());
342  hresz->Fill(helixhit->GetZ() - mcpoint->GetZ());
343 
344  if(tube->GetWireDirection() == TVector3(0, 0, 1)) {
345  hx->Fill(helixhit->GetX() - mcpoint->GetX());
346  hy->Fill(helixhit->GetY() - mcpoint->GetY());
347  hz->Fill(helixhit->GetZ() - mcpoint->GetZ());
348  }
349  else {
350  hxs->Fill(helixhit->GetX() - mcpoint->GetX());
351  hys->Fill(helixhit->GetY() - mcpoint->GetY());
352  hzs->Fill(helixhit->GetZ() - mcpoint->GetZ());
353  }
354 
355  }
356  }
357 }
358 // -------------------------------------------------------------------------
359 
360 
361 void PndSttTrackFitterQATask::WriteHistograms()
362 {
363  TFile* file = FairRootManager::Instance()->GetOutFile();
364  file->cd();
365  file->mkdir("SttFitterQATask");
366  file->cd("SttFitterQATask");
367 
368  hptfit->Write();
369  delete hptfit;
370  hplfit->Write();
371  delete hplfit;
372  hptotfit->Write();
373  delete hptotfit;
374 
375  hptfit_perc->Write();
376  delete hptfit_perc;
377  hplfit_perc->Write();
378  delete hplfit_perc;
379  hptotfit_perc->Write();
380  delete hptotfit_perc;
381 
382  hpxfit->Write();
383  delete hpxfit;
384  hpyfit->Write();
385  delete hpyfit;
386  hpzfit->Write();
387  delete hpzfit;
388 
389  hDist->Write();
390  delete hDist;
391  hRad->Write();
392  delete hRad;
393  hPhi->Write();
394  delete hPhi;
395  hTanL->Write();
396  delete hTanL;
397  hZ->Write();
398  delete hZ;
399  hQ->Write();
400  delete hQ;
401 
402 
403  hptfound->Write();
404  delete hptfound;
405  hplfound->Write();
406  delete hplfound;
407  hptotfound->Write();
408  delete hptotfound;
409 
410  hptfound_perc->Write();
411  delete hptfound_perc;
412  hplfound_perc->Write();
413  delete hplfound_perc;
414  hptotfound_perc->Write();
415  delete hptotfound_perc;
416 
417  hresx->Write();
418  delete hresx;
419  hresy->Write();
420  delete hresy;
421  hresz->Write();
422  delete hresz;
423 
424  hx->Write();
425  delete hx;
426  hy->Write();
427  delete hy;
428  hz->Write();
429  delete hz;
430 
431  hxs->Write();
432  delete hxs;
433  hys->Write();
434  delete hys;
435  hzs->Write();
436  delete hzs;
437 
438 }
439 
440 
441 
442 ClassImp(PndSttTrackFitterQATask)
PndSttTrackFitterQATask::Exec
virtual void Exec(Option_t *opt)
Definition: PndSttTrackFitterQATask.cxx:187
fVerbose
int fVerbose
Definition: poormantracks.C:24
cos
friend F32vec4 cos(const F32vec4 &a)
Definition: P4_F32vec4.h:112
PndSttTrack::GetTanL
Double_t GetTanL()
Definition: PndSttTrack.h:60
PndSttHelixHit::GetHitIndex
Int_t GetHitIndex()
Definition: PndSttHelixHit.h:54
PndSttTrack::GetNofHelixHits
Int_t GetNofHelixHits() const
Definition: PndSttTrack.h:48
PndTrackCand::getMcTrackId
int getMcTrackId() const
Definition: PndTrackCand.h:60
file
TFile * file
Definition: drc/old/draw_geom.C:12
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
PndSttTrack::GetTrackCandIndex
Int_t GetTrackCandIndex()
Definition: PndSttTrack.h:45
CAMath::Sqrt
static T Sqrt(const T &x)
Definition: PndCAMath.h:37
PndMCTrack::GetPdgCode
Int_t GetPdgCode() const
Definition: PndMCTrack.h:73
verbose
#define verbose
Definition: PndRadMapBoxMesh.cxx:15
sin
friend F32vec4 sin(const F32vec4 &a)
Definition: P4_F32vec4.h:111
CAMath::Sin
static T Sin(const T &x)
Definition: PndCAMath.h:42
ptot
Double_t ptot
Definition: checkhelixhit.C:68
PndMCTrack::GetMomentum
TVector3 GetMomentum() const
Definition: PndMCTrack.h:78
CAMath::Cos
static T Cos(const T &x)
Definition: PndCAMath.h:43
PndSttTrack::GetDist
Double_t GetDist()
Definition: PndSttTrack.h:57
d0
Double_t d0
Definition: checkhelixhit.C:59
PndSttTrack::GetRad
Double_t GetRad()
Definition: PndSttTrack.h:59
PndSttTrack::GetPhi
Double_t GetPhi()
Definition: PndSttTrack.h:58
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
phi0
Double_t phi0
Definition: checkhelixhit.C:60
CAMath::ATan2
static T ATan2(const T &y, const T &x)
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
PndSttTrack::GetCharge
Int_t GetCharge()
Definition: PndSttTrack.h:63
ptran
Double_t ptran
Definition: checkhelixhit.C:65
atan2
friend F32vec4 atan2(const F32vec4 &y, const F32vec4 &x)
Definition: P4_F32vec4.h:117
PndSttTrack::GetHelixHitIndex
Int_t GetHelixHitIndex(Int_t iHit) const
Definition: PndSttTrack.h:49
ClassImp
ClassImp(PndAnaContFact)
hit
PndSdsMCPoint * hit
Definition: anasim.C:70
PndTrack::GetRefIndex
Int_t GetRefIndex() const
Definition: PndTrack.h:36
alpha
double alpha
Definition: f_Init.h:9
PndMCTrack::GetStartVertex
TVector3 GetStartVertex() const
Definition: PndMCTrack.h:76
Pi
Double_t Pi
Definition: PndGeoHypGeCrystal.cxx:10
PndSttTrack::GetFlag
Int_t GetFlag() const
Definition: PndSttTrack.h:51
vec
dble_vec_t vec[12]
Definition: ranlxd.cxx:380
PndSttTube::GetWireDirection
TVector3 GetWireDirection()
Definition: PndSttTube.cxx:107
PndTrack::GetParamFirst
FairTrackParP GetParamFirst()
Definition: PndTrack.h:49
plong
Double_t plong
Definition: checkhelixhit.C:67
point
PndSdsMCPoint * point
Definition: anaLmdCluster.C:72
Generated on Wed Apr 3 2019 10:02:58 for FairRoot/PandaRoot by   doxygen 1.8.5