FairRoot/PandaRoot
TtAliTask.cxx
Go to the documentation of this file.
1 // Simone Bianco 14/07/2010
2 
3 // This Class' Header ------------------
4 #include "TtAliTask.h"
5 
6 // C/C++ Headers ----------------------
7 #include <iostream>
8 
9 // Collaborating Class Headers --------
10 #include "FairRootManager.h"
11 #include "TClonesArray.h"
12 
13 #include "PndTrackCand.h"
14 #include "PndSdsHit.h"
15 #include "PndTrackCandHit.h"
16 
17 #include "TFile.h"
18 #include "TGeoTrack.h"
19 #include "TGeoManager.h"
20 #include "TLorentzVector.h"
21 #include "TVector2.h"
22 
23 // Fit Classes -----------
24 #include <TMath.h>
25 #include <TVector3.h>
26 #include <TRandom2.h>
27 //#include <TStyle.h>
28 #include <TCanvas.h>
29 #include <TF2.h>
30 #include <TH1.h>
31 #include <TVirtualFitter.h>
32 #include <TPolyLine3D.h>
33 #include <Math/Vector3D.h>
34 #include <TGraphErrors.h>
35 #include <TCanvas.h>
36 
37 #include <fstream>
38 
39 using namespace ROOT::Math;
40 using namespace std;
41 
42 TtAliTask::TtAliTask()
43  : FairTask("Alignment") ,
44  fTCandArray(NULL),
45  fTCandBranchName("MVDHitsStrip"),
46  fTrackcount(),
47  fEvent(),
48  fExclBox(),
49  fPrint(0),
50  sX(),
51  sY(),
52  sigX(),
53  sigY(),
54  m_X(),
55  m_Y(),
56  hx(NULL),
57  hy(NULL)
58 {
59  for (Int_t gg = 0 ; gg < 4 ; gg++)
60  {
61  sX[gg] = 0.;
62  sY[gg] = 0.;
63  sigX[gg] = 0.;
64  sigY[gg] = 0.;
65  m_X[gg] = 0.;
66  m_Y[gg] = 0.;
67  }
68 }
69 
70 TtAliTask::TtAliTask(Int_t ExcludedBox)
71  : FairTask("3D-Straight-Line-Fit"),
72  fTCandArray(NULL),
73  fTCandBranchName("MVDHitsStrip"),
74  fTrackcount(),
75  fEvent(),
76  fExclBox(ExcludedBox),
77  fPrint(0),
78  sX(),
79  sY(),
80  sigX(),
81  sigY(),
82  m_X(),
83  m_Y(),
84  hx(NULL),
85  hy(NULL)
86 {
87  if (fExclBox < 1 || fExclBox > 6)
88  {
89  std::cout << "Excluded box: Wrong value, setting as default 2!" << std::endl;
90  fExclBox = 2;
91  }
92 
93  for (Int_t gg = 0 ; gg < 6 ; gg++)
94  {
95  sX[gg] = 0.;
96  sY[gg] = 0.;
97  sigX[gg] = 0.;
98  sigY[gg] = 0.;
99  m_X[gg] = 0.;
100  m_Y[gg] = 0.;
101  }
102 }
103 
104 TtAliTask::~TtAliTask()
105 {
106 
107 }
108 
109 InitStatus TtAliTask::Init()
110 {
111 
112  //Get ROOT Manager
113  FairRootManager* ioman= FairRootManager::Instance();
114 
115  if(ioman==0)
116  {
117  Error("TtAliTask::Init","RootManager not instantiated!");
118  return kERROR;
119  }
120 
121  //Get input collection
122 
123  fTCandArray=(TClonesArray*) ioman->GetObject(fTCandBranchName);
124 
125  if(fTCandArray==0)
126  {
127  Error("TtAliTask::Init","trackcand-array not found!");
128  return kERROR;
129  }
130 
131  std::cout << "-I- TtAliTask: Initialisation successfull" << std::endl;
132 
133  fEvent = 0;
134 
135  hx = new TH1F("hx","hx",50000,-5.,+5.);
136  hy = new TH1F("hy","hy",50000,-5.,+5.);
137 
138  return kSUCCESS;
139 }
140 
141 
142 void TtAliTask::Exec(Option_t*)
143 {
144 
145  // std::cout << "Event: " << fEvent << std::endl;
146 
147  Int_t ntcand=fTCandArray->GetEntriesFast();
148 
149  if (ntcand == 6)
150  {
151 
152  // std::cout << "Entries: " << ntcand << std::endl;
153 
154 
155  std::map<Double_t,Int_t> SensorsPos;
156 
157 
158  // std::cout<<"TtAliTask::Exec"<<std::endl;
159  // Reset output Array
160  if(fTCandArray==0) Fatal("TtAliTask::Exec","No Points");
161 
162 
163 
164 
165  Int_t name;
166 
167 
168  // Detailed output
169  // std::cout<<" -I- TtAliTask: contains "<<ntcand<<" TCandidates"<<std::endl;
170 
171  //std::cout<< " Detailed Debug info on the candidates:"<<std::endl;
172 
173  for(Int_t itr=0;itr<ntcand;++itr){
174 
175  PndSdsHit *theHit = (PndSdsHit*) fTCandArray->At(itr);
176 
177  name = theHit->GetSensorID();
178 
179 
180  if (SensorsPos.size() < 6) SensorsPos[theHit->GetZ()] = name;
181 
182 
183  }
184 
185 
186  // const Int_t sizeMap = Sensors.size();
187 
188  const Int_t sizeMap = SensorsPos.size();
189 
190  Int_t DetNames[sizeMap];
191  //Double_t Pos[sizeMap]; //[R.K.03/2017] unused variable
192  Int_t jj = 0;
193 
194  if (SensorsPos.size()>=6)
195  {
196  //std::cout << "Number of Sensors:" << SensorsPos.size() << std::endl;
197  for (std::map<Double_t,Int_t>::iterator it=SensorsPos.begin();it!=SensorsPos.end();++it)
198  {
199  if (fEvent < 5) std::cout << "position: " << it->first << " name: " << (it->second) << std::endl;
200 
201  DetNames[jj] = it->second;
202  //Pos[jj] = it -> first; //[R.K.03/2017] unused variable
203  jj++;
204  }
205 
206  }
207 
208 
209 
210 
211  // if (fEvent < 5)
212 // {
213 
214 // for (Int_t l = 0 ; l < sizeMap ; l++)
215 // {
216 // std::cout << DetNames[l] << std::endl;
217 // }
218 // }
219 
220 
222 
223 
224  if (SensorsPos.size()==6)
225  {
226  //Int_t DetRem[sizeMap-1]; //[R.K.02/2017] Unused variable?
227  //for (Int_t ss = 0 ; ss < (sizeMap-1) ; ss++) //[R.K.02/2017] Unused variable?
228  //{ //[R.K.02/2017] Unused variable?
229 
230  //DetRem[ss] = 0; //[R.K.02/2017] Unused variable?
231 
232  //} //[R.K.02/2017] Unused variable?
233 
234 
235  //Int_t counter = 0; //[R.K.02/2017] Unused variable?
236 
237  //Double_t BuffZ = -9999.;
238 
239  //for (Int_t kk = 0 ; kk < (sizeMap) ; kk++) //[R.K.02/2017] Unused variable?
240  //{ //[R.K.02/2017] Unused variable?
241 
242  //if ((fExclBox-1) == kk) //[R.K.02/2017] Unused variable?
243  //{ //[R.K.02/2017] Unused variable?
244  //BuffZ = Pos[kk]; //[R.K.02/2017] Unused variable?
245  //continue; //[R.K.02/2017] Unused variable?
246  //} //[R.K.02/2017] Unused variable?
247  //else //[R.K.02/2017] Unused variable?
248  //{ //[R.K.02/2017] Unused variable?
249  //DetRem[counter] = kk; //[R.K.02/2017] Unused variable?
250  //counter++; //[R.K.02/2017] Unused variable?
251  //} //[R.K.02/2017] Unused variable?
252 
253  //} //[R.K.02/2017] Unused variable?
254 
255 
256  Double_t x[6];
257  Double_t y[6];
258  Double_t z[6];
259 
260  Double_t Erx[6];
261  Double_t Ery[6];
262  Double_t Erz[6];
263 
264  //Int_t track = 0; //[R.K. 01/2017] unused variable?
265 
266 
267  //Double_t RealX = -1999999., RealY = -1999999;
268 
269 
270  //Double_t buffErrX = -0.099999; //[R.K. 01/2017] unused variable?
271  //Double_t buffErrY = -0.099999; //[R.K. 01/2017] unused variable?
272 
273 
274  Int_t ihit = 0;
275 
276  for (Int_t it1 = 0 ; it1 < ntcand ; it1++)
277 
278  {
279 
280  PndSdsHit *theHit = (PndSdsHit*) fTCandArray->At(it1);
281 
282  if ((theHit->GetSensorID()) == DetNames[0])
283  {
284  x[0] = theHit->GetX();
285  y[0] = theHit->GetY();
286  z[0] = theHit->GetZ();
287  Erx[0] = theHit->GetDx();
288  Ery[0] = theHit->GetDy();
289  Erz[0] = theHit->GetDz();
290  }
291 
292  if ((theHit->GetSensorID()) == DetNames[1])
293  {
294  x[1] = theHit->GetX();
295  y[1] = theHit->GetY();
296  z[1] = theHit->GetZ();
297  Erx[1] = theHit->GetDx();
298  Ery[1] = theHit->GetDy();
299  Erz[1] = theHit->GetDz();
300  }
301  if ((theHit->GetSensorID()) == DetNames[2])
302  {
303  x[2] = theHit->GetX();
304  y[2] = theHit->GetY();
305  z[2] = theHit->GetZ();
306  Erx[2] = theHit->GetDx();
307  Ery[2] = theHit->GetDy();
308  Erz[2] = theHit->GetDz();
309  }
310 
311  if ((theHit->GetSensorID()) == DetNames[3])
312  {
313  x[3] = theHit->GetX();
314  y[3] = theHit->GetY();
315  z[3] = theHit->GetZ();
316  Erx[3] = theHit->GetDx();
317  Ery[3] = theHit->GetDy();
318  Erz[3] = theHit->GetDz();
319  }
320  if ((theHit->GetSensorID()) == DetNames[4])
321  {
322  x[4] = theHit->GetX();
323  y[4] = theHit->GetY();
324  z[4] = theHit->GetZ();
325  Erx[4] = theHit->GetDx();
326  Ery[4] = theHit->GetDy();
327  Erz[4] = theHit->GetDz();
328  }
329 
330  if ((theHit->GetSensorID()) == DetNames[5])
331  {
332  x[5] = theHit->GetX();
333  y[5] = theHit->GetY();
334  z[5] = theHit->GetZ();
335  Erx[5] = theHit->GetDx();
336  Ery[5] = theHit->GetDy();
337  Erz[5] = theHit->GetDz();
338  }
339 
340  ihit++;
341 
342 
343  } // end loop on points
344 
345 
346  // setting the points
347 
348 
349 
350 
351  Int_t uu = 0;
352 
353  // to check alignment...
354 
355  for (Int_t ww = 0 ; ww < 6 ; ww++)
356  {
357  if (TMath::Abs(Erx[ww]) < 0.5)
358  {
359  x[ww]=x[ww]+sX[ww];
360  }
361  if (TMath::Abs(Ery[ww]) < 0.5)
362  {
363  y[ww]=y[ww]+sY[ww];
364  }
365  }
366 
367 
368  for (Int_t ww = 0 ; ww < 6 ; ww++)
369  {
370  if ((fExclBox-1)!=ww)
371  {
372  if (TMath::Abs(Erx[ww]) > 0.5) Erx[ww] = 1000000;
373  if (TMath::Abs(Ery[ww]) > 0.5) Ery[ww] = 1000000;
374  uu++;
375  }
376 
377  // std::cout << "(" << x[ww] << "," << y[ww] << "," << z[ww] << ")" << std::endl;
378 
379  }
380 
381 
382  Double_t DX,DY;
383 
384 
385  MyFit(x,y,z,Erx,Ery,Erz,x[fExclBox-1],y[fExclBox-1],z[fExclBox-1],DX,DY);
386 
387  //Double_t pointX,pointY,pointZ; //[R.K. 01/2017] unused variable?
388 
389  hx->Fill(DX);
390  hy->Fill(DY);
391 
392  }// if 6 hits
393 
394  fEvent++;
395 
396  }
397 
398  return;
399 
400 }
401 
402 
403 void TtAliTask::MyFit(Double_t *x,Double_t *y,Double_t *z,Double_t *Erx,Double_t *Ery,Double_t *,Double_t realX, Double_t realY, Double_t realZ, Double_t &DELTAX, Double_t &DELTAY) //Erz //[R.K.03/2017] unused variable(s)
404 {
405 
406  TGraphErrors grX;
407  TGraphErrors grY;
408 
409  Int_t ix=0,iy=0;
410 
411  for (Int_t j = 0 ; j < 6 ; j++)
412  {
413  //std::cout << j << "@(" << Erx[j] << "," << Ery[j] << "," << Erz[j] << ")" << std::endl;
414  if (j == fExclBox-1) continue;
415 
416  if (TMath::Abs(Erx[j])<0.5)
417  //if (j==0 || j==1 || j==3 || j==5)
418  {
419  grX.SetPoint(ix,z[j],x[j]);
420  ix++;
421  }
422  //if(j==0 || j==2 || j==4 || j==5)
423  if (TMath::Abs(Ery[j])<0.5)
424  {
425  grY.SetPoint(iy,z[j],y[j]);
426  iy++;
427  }
428 
429  }
430 
431  // std::cout << "POINTS " << ix << " " << iy << std::endl;
432 
433 
434  Double_t mx,nx,my,ny;
435 
436  grX.Fit("pol1","Q");
437 
438  nx = (grX.GetFunction("pol1"))->GetParameter(0);
439  mx = (grX.GetFunction("pol1"))->GetParameter(1);
440 
441 
442  grY.Fit("pol1","Q");
443 
444  ny = (grY.GetFunction("pol1"))->GetParameter(0);
445  my = (grY.GetFunction("pol1"))->GetParameter(1);
446 
447  Double_t buffX, buffY;
448 
449  buffX = mx * realZ + nx;
450  buffY = my * realZ + ny;
451 
452  DELTAX = buffX - realX;
453  DELTAY = buffY - realY;
454 
455 
456 
457 }
458 
459 
460 
461 
462 void TtAliTask::FinishTask()
463 {
464 
465 
466  Double_t mX,mY,siX,siY;
467 
468  TF1 *fun = new TF1("fun","gaus");
469 
470  hx->Fit(fun,"Q","",(hx->GetMean() - 2*(hx->GetRMS())),(hx->GetMean() + 2*(hx->GetRMS())));
471 
472  //hx->Fit(fun,"Q");
473 
474  mX = fun->GetParameter(1);
475  siX = fun->GetParameter(2);
476 
477  cout << "RMS: " << hx->GetRMS() << " sigma " << siX << endl;
478 
479  hy->Fit(fun,"Q","",(hy->GetMean() - 2*(hy->GetRMS())),(hy->GetMean() + 2*(hy->GetRMS())));
480 
481  //hy->Fit(fun,"Q");
482 
483  mY = fun->GetParameter(1);
484  siY = fun->GetParameter(2);
485 
486  cout << "X, mean: " << mX << ", sig: " << siX << endl;
487  cout << "Y, mean: " << mY << ", sig: " << siY << endl;
488 
489  if ((TMath::Abs(siX)) < 1.)
490  {
491  sX[fExclBox-1] += mX;
492  m_X[fExclBox-1] = mX;
493  sigX[fExclBox-1] = siX;
494  }
495  else
496  {
497  sX[fExclBox-1] = 0.;
498  m_X[fExclBox-1] = 0.;
499  sigX[fExclBox-1] = 0.;
500  }
501 
502  if ((TMath::Abs(siY)) < 1.)
503  {
504  sY[fExclBox-1] += mY;
505  m_Y[fExclBox-1] = mY;
506  sigY[fExclBox-1] = siY;
507  }
508  else
509  {
510  sY[fExclBox-1] = 0.;
511  m_Y[fExclBox-1] = 0.;
512  sigY[fExclBox-1] = 0.;
513  }
514 
515 
516  /* switch(fExclBox)
517  {
518  case 1:
519  sX[0] += mX;
520  sY[0] += mY;
521  m_X[0] = mX;
522  m_Y[0] = mY;
523  sigX[0] = siX;
524  sigY[0] = siY;
525  break;
526 
527  case 2:
528  sX[1] += mX;
529  sigX[1] = siX;
530  m_X[1] = mX;
531 
532  break;
533 
534  case 3:
535  sY[1] += mY;
536  sigY[1] = siY;
537  m_Y[1] = mY;
538  break;
539 
540  case 4:
541  sX[2] += mX;
542  sigX[2] = siX;
543  m_X[2] = mX;
544  break;
545 
546  case 5:
547  sY[2] += mY;
548  sigY[2] = siY;
549  m_Y[2] = mY;
550  break;
551 
552  case 6:
553  sX[3] += mX;
554  sY[3] += mY;
555  sigX[3] = siX;
556  sigY[3] = siY;
557  m_X[3] = mX;
558  m_Y[3] = mY;
559  break;
560  }
561 
562  */
563  if (fPrint == 1)
564  {
565  PrintHistos();
566  }
567 
568  hx->Reset();
569  hy->Reset();
570 
571  // fEvent = 0;
572 }
573 
574 
575 void TtAliTask::PrintVal()
576 {
577 
578  cout << "SHIFTS" << endl;
579  for (Int_t k = 0 ; k < 6 ; k++)
580  {
581  cout << "X: " << sX[k] << " Y: " << sY[k] << endl;
582  }
583 }
584 
585 
586 
587 
588 void TtAliTask::PrintMeanResiduals()
589 {
590 
591  cout << "RESIDUALS" << endl;
592  for (Int_t k = 0 ; k < 6 ; k++)
593  {
594  cout << "X: " << m_X[k] << " Y: " << m_Y[k] << endl;
595  }
596 }
597 
598 void TtAliTask::PrintSigmaResiduals()
599 {
600 
601  cout << "SIGMA-RESIDUALS" << endl;
602  for (Int_t k = 0 ; k < 6 ; k++)
603  {
604  cout << "sigX: " << sigX[k] << " sigY: " << sigY[k] << endl;
605  }
606 }
607 
608 TVector2 TtAliTask::GetRes()
609 {
610 
611  Double_t resx = TMath::Power((sigX[0]*sigX[1]*sigX[2]*sigX[3]),1./4.);
612  Double_t resy = TMath::Power((sigY[0]*sigY[1]*sigY[2]*sigY[3]),1./4.);
613 
614  return TVector2(resx,resy);
615 
616 }
617 
618 void TtAliTask::PrintHistos()
619 {
620 
621  TCanvas *can = new TCanvas();
622 
623  can->cd();
624 
625  hx->Draw();
626 
627  hx->GetXaxis()->SetRangeUser((hx->GetMean())-8*(hx->GetRMS()),(hx->GetMean())+8*(hx->GetRMS()));
628 
629  std::string nameX = Form("HistResidualsX_%d.png",fExclBox);
630 
631  can->Print(nameX.c_str(),"png");
632 
633 
634  hy->Draw();
635 
636  hy->GetXaxis()->SetRangeUser((hy->GetMean())-8*(hy->GetRMS()),(hy->GetMean())+8*(hy->GetRMS()));
637 
638  std::string nameY = Form("HistResidualsY_%d.png",fExclBox);
639 
640  can->Print(nameY.c_str(),"png");
641 
642 
643 }
644 
645 
646 
647 ClassImp(TtAliTask);
TtAliTask::Exec
virtual void Exec(Option_t *opt)
Definition: TtAliTask.cxx:142
TtAliTask::Init
virtual InitStatus Init()
Definition: TtAliTask.cxx:109
TtAliTask::sY
Double_t sY[6]
Definition: TtAliTask.h:66
TtAliTask::fExclBox
Int_t fExclBox
Definition: TtAliTask.h:62
TtAliTask::hx
TH1F * hx
Definition: TtAliTask.h:70
TtAliTask::m_Y
Double_t m_Y[6]
Definition: TtAliTask.h:68
GetParameter
cout<< "ifile "<< ifile<< endl;cout<< " momentum sampled over "<< nsteps<< " with step width "<< 1.5/nsteps<< endl;cout<< endl;cout<< "MEAN DEDX PARAMETRIZATION"<< endl;cout<< "mom limits "<< mean_inf<< " "<< mean_sup<< endl;cout<< "mu: ";for(int param=0;param< npardedx;param++) cout<< fdedx-> GetParameter(param)<< "
TtAliTask::sX
Double_t sX[6]
Definition: TtAliTask.h:66
TtAliTask::m_X
Double_t m_X[6]
Definition: TtAliTask.h:68
jj
Int_t jj
Definition: createRootGeoFileBwEnd_FwEndCapGeo.C:434
TtAliTask::sigX
Double_t sigX[6]
Definition: TtAliTask.h:67
CAMath::Abs
static T Abs(const T &x)
Definition: PndCAMath.h:39
TtAliTask::FinishTask
virtual void FinishTask()
Definition: TtAliTask.cxx:462
TtAliTask::GetRes
TVector2 GetRes()
Definition: TtAliTask.cxx:608
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
TtAliTask::sigY
Double_t sigY[6]
Definition: TtAliTask.h:67
TtAliTask::PrintHistos
void PrintHistos()
Definition: TtAliTask.cxx:618
z
Double_t z
Definition: createRootFscGeometryFile.C:146
TtAliTask::fTCandArray
TClonesArray * fTCandArray
Definition: TtAliTask.h:50
name
TString name
Definition: createRootGeoFile.C:82
TtAliTask::~TtAliTask
virtual ~TtAliTask()
Definition: TtAliTask.cxx:104
TtAliTask::PrintSigmaResiduals
void PrintSigmaResiduals()
Definition: TtAliTask.cxx:598
TtAliTask::hy
TH1F * hy
Definition: TtAliTask.h:71
x
Double_t x
Definition: createRootFscGeometryFile.C:159
TtAliTask::fPrint
Int_t fPrint
Definition: TtAliTask.h:64
ClassImp
ClassImp(TtAliTask)
TtAliTask::PrintVal
void PrintVal()
Definition: TtAliTask.cxx:575
y
Double_t y
Definition: createRootFscGeometryFile.C:159
TtAliTask::MyFit
void MyFit(Double_t *x, Double_t *y, Double_t *z, Double_t *Erx, Double_t *Ery, Double_t *Erz, Double_t realX, Double_t realY, Double_t realZ, Double_t &DELTAX, Double_t &DELTAY)
Definition: TtAliTask.cxx:403
PndSdsHit::GetSensorID
Int_t GetSensorID() const
Definition: PndSdsHit.h:90
TtAliTask::fEvent
Int_t fEvent
Definition: TtAliTask.h:61
TtAliTask::fTCandBranchName
TString fTCandBranchName
Definition: TtAliTask.h:56
TtAliTask::PrintMeanResiduals
void PrintMeanResiduals()
Definition: TtAliTask.cxx:588
Generated on Sun Aug 18 2019 03:02:19 for FairRoot/PandaRoot by   doxygen 1.8.8