FairRoot/PandaRoot
TtCracowTask.cxx
Go to the documentation of this file.
1 // This Class' Header ------------------
2 #include "TtCracowTask.h"
3 #include "TtFitRes.h"
4 
5 // C/C++ Headers ----------------------
6 #include <iostream>
7 
8 // Collaborating Class Headers --------
9 #include "FairRootManager.h"
10 #include "TClonesArray.h"
11 #include "TGraph2DErrors.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 
22 // Fit Classes -----------
23 #include <TMath.h>
24 #include <TVector3.h>
25 #include <TRandom2.h>
26 //#include <TStyle.h>
27 #include <TCanvas.h>
28 #include <TF2.h>
29 #include <TH1.h>
30 #include <TVirtualFitter.h>
31 #include <TPolyLine3D.h>
32 #include <Math/Vector3D.h>
33 #include <TGraphErrors.h>
34 #include <TCanvas.h>
35 
36 #include <fstream>
37 
38 
39 //using namespace ROOT::Math;
40 using namespace std;
41 
42 TtCracowTask::TtCracowTask()
43  : FairTask("Cracow"),
44  fTCandArray(NULL),
45  fTCandBranchName("MVDHitsStrip"),
46  fTrackArray(NULL),
47  fTrackcount(0),
48  fEvent(0),
49  fEloss()
50 {
51  for(int i=0;i<6;i++)fEloss[i]=0.;
52 }
53 
54 
55 TtCracowTask::~TtCracowTask()
56 {
57 
58 }
59 
60 InitStatus TtCracowTask::Init()
61 {
62 
63 
64  //Get ROOT Manager
65  FairRootManager* ioman= FairRootManager::Instance();
66 
67  if(ioman==0)
68  {
69  Error("TtCracowTask::Init","RootManager not instantiated!");
70  return kERROR;
71  }
72 
73  //Get input collection
74 
75  fTCandArray=(TClonesArray*) ioman->GetObject(fTCandBranchName);
76 
77  if(fTCandArray==0)
78  {
79  Error("TtCracowTask::Init","trackcand-array not found!");
80  return kERROR;
81  }
82 
83  std::cout << "Input: " << fTCandArray->GetClass()->GetName() << std::endl;
84 
85 
86  fTrackArray = new TClonesArray("TtFitRes");
87  ioman->Register("TTFit", "PndMvd", fTrackArray, kTRUE);
88 
89  std::cout << "-I- TtCracowTask: Initialisation successfull" << std::endl;
90 
91  fEvent = 0;
92 
93  return kSUCCESS;
94 }
95 
96 
97 // line parametrization
98 // x = par[1]*t +par[0]
99 // y = par[3]*t +par[2]
100 // z = t
101 
102 Double_t TtCracowTask::distance2(double x,double y,double z,double *par) {
103  // distance line point is D= | (xp-x0) cross ux |
104  // where ux is direction of line and x0 is a point in the line (like t = 0)
105  TVector3 xp(x,y,z);
106  TVector3 x0(par[0],par[2],0);
107  TVector3 x1(par[0]+par[1],par[2]+par[3],1);
108  TVector3 u = (x1-x0).Unit();
109  double d2 = ((xp-x0).Cross(u)) .Mag2();
110 
111  //std::cout << "Root of the distance: " << TMath::Sqrt(d2) << std::endl;
112 
113  return d2;
114 }
115 
116 Double_t TtCracowTask::distance2Single(double x,double y,double z,double ex,double ey,double ,double *par) {// ez //[R.K.03/2017] unused variable(s)
117  // distance line point is D= | (xp-x0) cross ux |
118  // where ux is direction of line and x0 is a point in the line (like t = 0)
119 
120  TVector3 ddd(0,0,0);
121  double d2=-999999.;
122 
123  if (ex < 10000.) // this means only x
124  {
125  ddd.SetXYZ(x-(par[1]*z +par[0]),0.,0.);
126  d2 = ddd.Mag2();
127  }
128 
129  if (ey < 10000.) // this means only y
130  {
131  ddd.SetXYZ(0.,y-(par[3]*z +par[2]),0.);
132  d2 = ddd.Mag2();
133  }
134 
135  //std::cout << "Root of the distance: " << TMath::Sqrt(d2) << std::endl;
136 
137  return d2;
138 }
139 
140 void TtCracowTask::SumDistance2(int &, double *, double & sum, double * par, int ) {
141 
142  TGraph2D * gr = dynamic_cast<TGraph2D*>( (TVirtualFitter::GetFitter())->GetObjectFit() );
143  assert(gr != 0);
144  double * x = gr->GetX();
145  double * y = gr->GetY();
146  double * z = gr->GetZ();
147  double * Ex = gr->GetEX();
148  double * Ey = gr->GetEY();
149  double * Ez = gr->GetEZ();
150 
151  int npoints = gr->GetN();
152  sum = 0;
153  double d = 0.;
154 
155  //for (int i = 0; i < 3 ; ++i) {
156  for (int i = 0; i < npoints ; ++i) {
157  //std::cout << "Sensor: " << i;
158  if ((Ex[i] < 10000.) && (Ey[i] < 10000.))
159  {
160  //std::cout << " - Double Sided" << std::endl;
161  d = distance2(x[i],y[i],z[i],par);
162  }
163 
164  if ((Ex[i] > 10000.) && (Ey[i] > 10000.))
165  {
166  //std::cout << " - Wrong Errors!" << std::endl;
167  d = 0.;
168  }
169 
170  if ((Ex[i] > 10000.) || (Ey[i] > 10000.))
171  {
172  //std::cout << " - Single Sided" << std::endl;
173  d = distance2Single(x[i],y[i],z[i],Ex[i],Ey[i],Ez[i],par);
174  }
175 
176 
177 
178  sum += d;
179  }
180 }
181 
182 
183 
184 void TtCracowTask::Exec(Option_t*)
185 {
186 
187 
188 
189  fTrackcount = 0;
190 
191  Bool_t s0=kFALSE;
192  Bool_t s1=kFALSE;
193  Bool_t s2=kFALSE;
194  Bool_t s3=kFALSE;
195  Bool_t s4=kFALSE;
196  Bool_t s5=kFALSE;
197 
198  // resetting energy losses
199 
200  for (Int_t qq = 0 ; qq < 6 ; qq++)
201  {
202  fEloss[qq]=0;
203  }
204 
205  Int_t ntcand=fTCandArray->GetEntriesFast();
206 
207  if (ntcand == 6)
208  {
209 
210  // std::cout << "Entries: " << ntcand << std::endl;
211 
212 
213 
214  std::map<Double_t,Int_t> SensorsPos;
215 
216  std::cout << "Event: " << fEvent << std::endl;
217  std::cout << "Initial map size: " << SensorsPos.size() << std::endl;
218 
219 // std::cout<<"TtCracowTask::Exec"<<std::endl;
220  // Reset output Array
221  if(fTrackArray==0) Fatal("TtCracowTask::Exec","No TrackArray");
222 
223  Int_t name;
224 
225 
226  for(Int_t itr=0;itr<ntcand;++itr){
227 
228  PndSdsHit *theHit = (PndSdsHit*) fTCandArray->At(itr);
229 
230  name = theHit->GetSensorID();
231 
232  if (SensorsPos.size() < 6) SensorsPos[theHit->GetZ()] = name;
233 
234  }
235 
236  const Int_t sizeMap = (Int_t) SensorsPos.size();
237 
238  Int_t DetNames[sizeMap];
239  //Double_t Pos[sizeMap]; //[R.K.02/2017] Unused variable?
240 
241 
242  Int_t jj = 0;
243 
244  if (SensorsPos.size()>=6)
245  {
246 
247  std::cout << "Map size after looping: " << SensorsPos.size() << std::endl;
248  //std::cout << "Number of Sensors:" << SensorsPos.size() << std::endl;
249  for (std::map<Double_t,Int_t>::iterator it=SensorsPos.begin();it!=SensorsPos.end();++it)
250  {
251 
252  std::cout << "position: " << it->first << " name: " << (it->second) << std::endl;
253 
254  DetNames[jj] = it->second;
255  //Pos[jj] = it -> first; //[R.K.02/2017] Unused variable?
256  jj++;
257  }
258 
259  }
260 
261 
262 
263 
264  if (fEvent < 5) // just to check
265  {
266 
267  for (Int_t l = 0 ; l < sizeMap ; l++)
268  {
269  //std::cout << DetNames[l] << std::endl;
270  }
271  }
272 
273  if (SensorsPos.size()==6)
274  {
275 
276  Double_t x[6];
277  Double_t y[6];
278  Double_t z[6];
279 
280  Double_t Erx[6];
281  Double_t Ery[6];
282  //Double_t Erz[6]; //[R.K.02/2017] Unused variable?
283 
284  //Int_t track = 0; //[R.K. 01/2017] unused variable?
285 
286 
287  //Double_t buffErrX = -0.099999; //[R.K. 01/2017] unused variable?
288  //Double_t buffErrY = -0.099999; //[R.K. 01/2017] unused variable?
289 
290 
291  // std::cout << "Event: " << fEvent << std::endl;
292 
293 
294  Int_t ihit = 0;
295 
296  TGraph2DErrors * gr = new TGraph2DErrors();
297 
298  for (Int_t it1 = 0 ; it1 < ntcand ; it1++)
299  {
300 
301  PndSdsHit *theHit = (PndSdsHit*) fTCandArray->At(it1);
302 
303  if ((theHit->GetSensorID()) == DetNames[0])
304  {
305  x[0] = theHit->GetX();
306  y[0] = theHit->GetY();
307  z[0] = theHit->GetZ();
308  Erx[0] = theHit->GetDx();
309  Ery[0] = theHit->GetDy();
310  //Erz[0] = theHit->GetDz(); //[R.K.02/2017] Unused variable?
311  fEloss[0] = theHit->GetEloss();
312  s0 = kTRUE;
313  }
314 
315  if ((theHit->GetSensorID()) == DetNames[1])
316  {
317  x[1] = theHit->GetX();
318  y[1] = theHit->GetY();
319  z[1] = theHit->GetZ();
320  Erx[1] = theHit->GetDx();
321  Ery[1] = theHit->GetDy();
322  //Erz[1] = theHit->GetDz(); //[R.K.02/2017] Unused variable?
323  fEloss[1] = theHit->GetEloss();
324  s1 = kTRUE;
325  }
326  if ((theHit->GetSensorID()) == DetNames[2])
327  {
328  x[2] = theHit->GetX();
329  y[2] = theHit->GetY();
330  z[2] = theHit->GetZ();
331  Erx[2] = theHit->GetDx();
332  Ery[2] = theHit->GetDy();
333  //Erz[2] = theHit->GetDz(); //[R.K.02/2017] Unused variable?
334  fEloss[2] = theHit->GetEloss();
335  s2 = kTRUE;
336  }
337 
338  if ((theHit->GetSensorID()) == DetNames[3])
339  {
340  x[3] = theHit->GetX();
341  y[3] = theHit->GetY();
342  z[3] = theHit->GetZ();
343  Erx[3] = theHit->GetDx();
344  Ery[3] = theHit->GetDy();
345  //Erz[3] = theHit->GetDz(); //[R.K.02/2017] Unused variable?
346  fEloss[3] = theHit->GetEloss();
347  s3 = kTRUE;
348  }
349  if ((theHit->GetSensorID()) == DetNames[4])
350  {
351  x[4] = theHit->GetX();
352  y[4] = theHit->GetY();
353  z[4] = theHit->GetZ();
354  Erx[4] = theHit->GetDx();
355  Ery[4] = theHit->GetDy();
356  //Erz[4] = theHit->GetDz(); //[R.K.02/2017] Unused variable?
357  fEloss[4] = theHit->GetEloss();
358  s4 = kTRUE;
359  }
360 
361  if ((theHit->GetSensorID()) == DetNames[5])
362  {
363  x[5] = theHit->GetX();
364  y[5] = theHit->GetY();
365  z[5] = theHit->GetZ();
366  Erx[5] = theHit->GetDx();
367  Ery[5] = theHit->GetDy();
368  //Erz[5] = theHit->GetDz(); //[R.K.02/2017] Unused variable?
369  fEloss[5] = theHit->GetEloss();
370  s5 = kTRUE;
371  }
372 
373  ihit++;
374 
375 
376  } // end loop on points
377 
378 
379  //Int_t uu = 0; //[R.K. 01/2017] unused variable?
380 
381 
382  for (Int_t ww = 0 ; ww < 6 ; ww++)
383  { // setting big errors for the bottom side
384  {
385  if (TMath::Abs(Erx[ww]) > 0.5) Erx[ww] = 10000.;
386  if (TMath::Abs(Ery[ww]) > 0.5) Ery[ww] = 10000.;
387  }
388 
389  // std::cout << "(" << x[ww] << "," << y[ww] << "," << z[ww] << ")" << std::endl;
390 
391  }
392 
393 
394  //Double_t parFit[4]; //fit-parameter //[R.K. 01/2017] unused variable?
395 
396 
397  if (s0 && s1 && s2 && s3 && s4 && s5)
398  {
399  for (Int_t ss = 0 ; ss < 6 ; ss++)
400  {
401  gr->SetPoint(ss,x[ss],y[ss],z[ss]);
402  gr->SetPointError(ss,x[ss],y[ss],z[ss]);
403  }
404  }
405 
406  Int_t Npoint = gr->GetN();
407 
408  Double_t *xx = gr->GetX();
409  Double_t *yy = gr->GetY();
410  Double_t *zz = gr->GetZ();
411 
412  TVector3 dir;
413  dir.SetXYZ(xx[Npoint-1]-xx[0],yy[Npoint-1]-yy[0],zz[Npoint-1]-zz[0]);
414 
415  TVirtualFitter::SetDefaultFitter("Minuit");
416  TVirtualFitter *min = TVirtualFitter::Fitter(0,4);
417  min->SetObjectFit(gr);
418  min->SetFCN(*SumDistance2);//using local coordinate in FCN
419 
420  Double_t arglist[100];
421 
422  arglist[0] = 0;
423  min->ExecuteCommand("SET NOWarnings",arglist,1);
424  arglist[0] = -1;
425  min->ExecuteCommand("SET PRINT",arglist,1);
426 
427  //double pStart[4] = {xx[0],(dir.X())/(dir.Y()),zz[0],(dir.Z())/(dir.Y())}; //[R.K. 01/2017] unused variable?
428  //double pStartErr[4] = {gr->GetErrorX(0),0.1*(dir.X())/(dir.Y()),gr->GetErrorZ(0),0.1*(dir.Z())/(dir.Y())}; //[R.K. 01/2017] unused variable?
429 
430  /*min->SetParameter(0,"y0",pStart[0],pStartErr[0],0,0);
431  min->SetParameter(1,"A",pStart[1],pStartErr[1],0,0);
432  min->SetParameter(2,"z0",pStart[2],pStartErr[2],0,0);
433  min->SetParameter(3,"B",pStart[3],pStartErr[3],0,0);*/
434 
435  min->SetParameter(0,"y0",0,1,0,0);
436  min->SetParameter(1,"A",0,1,0,0);
437  min->SetParameter(2,"z0",0,1,0,0);
438  min->SetParameter(3,"B",0,1,0,0);
439 
440  arglist[0] = 1000; //number of functiona calls
441  arglist[1] = 0.001; //tolerance
442  min->ExecuteCommand("MIGRAD", arglist ,2);
443 
444  int nvpar,nparx;
445  double amin,edm, errdef;
446  min->GetStats(amin,edm,errdef,nvpar,nparx);
447  if(fVerbose>1) min->PrintResults(1,amin);
448 
449  //Double_t fitPar[4]; //[R.K. 01/2017] unused variable?
450  //Double_t fitParErr[4]; //[R.K. 01/2017] unused variable?
451 
452  Double_t chi2 = amin/(2.*Npoint-4);
453 
454  Double_t SumEn;
455 
456  for (Int_t gg = 0 ; gg < 6 ; gg++) SumEn += fEloss[gg];
457 
458  std::cout << "Scrivo" << endl;
459  new ((*fTrackArray)[0]) TtFitRes(min->GetParameter(0), min->GetParameter(1), min->GetParameter(2), min->GetParameter(3), SumEn, chi2, chi2, Npoint) ;
460 // new ((*fTrackArray)[0]) TtFitRes(min->GetParameter(0), min->GetParameter(1), min->GetParameter(2), min->GetParameter(3), fEloss, chi2, chi2, 3) ;
461  delete gr;
462 
463 
464  }// if 6 hits
465 
466 
467  }
468 
469  fEvent++;
470 
471 // delete gr;
472 
473  return;
474 
475 }
476 
477 
478 void TtCracowTask::MyFit(Double_t *x,Double_t *y,Double_t *z,Double_t *Erx,Double_t *Ery,Double_t *Erz,Double_t *par, Double_t &chiX, Double_t &chiY)
479 {
480 
481  TGraphErrors grX;
482  TGraphErrors grY;
483 
484  Int_t ix=0,iy=0;
485 
486  for (Int_t j = 0 ; j < 6 ; j++)
487  {
488  //std::cout << j << "@(" << Erx[j] << "," << Ery[j] << "," << Erz[j] << ")" << std::endl;
489 
490  if (TMath::Abs(Erx[j])<0.5)
491  {
492  grX.SetPoint(ix,z[j],x[j]);
493  grX.SetPointError(ix,TMath::Abs(Erz[j]),TMath::Abs(Erx[j]));
494  //cout << "erz : " << Erz[j] << " erx : " << Erx[j] << endl;
495  ix++;
496  }
497 
498  if (TMath::Abs(Ery[j])<0.5)
499  {
500  grY.SetPoint(iy,z[j],y[j]);
501  grY.SetPointError(iy,TMath::Abs(Erz[j]),TMath::Abs(Ery[j]));
502  //cout << "erz : " << Erz[j] << " ery : " << Ery[j] << endl;
503  iy++;
504  }
505 
506  }
507 
508  // std::cout << "POINTS " << ix << " " << iy << std::endl;
509 
510 
511  Double_t mx,nx,my,ny;
512 
513  grX.Fit("pol1","");
514 
515  nx = (grX.GetFunction("pol1"))->GetParameter(0);
516  mx = (grX.GetFunction("pol1"))->GetParameter(1);
517 
518 
519  grY.Fit("pol1","");//"Q"
520 
521  ny = (grY.GetFunction("pol1"))->GetParameter(0);
522  my = (grY.GetFunction("pol1"))->GetParameter(1);
523 
524 
525  par[0] = nx;
526  par[1] = mx;
527  par[2] = ny;
528  par[3] = my;
529 
530  chiX = (grX.GetFunction("pol1"))->GetChisquare();
531 
532  chiY = (grY.GetFunction("pol1"))->GetChisquare();
533 
534 }
535 
536 void TtCracowTask::FinishEvent()
537 {
538  fTrackArray->Clear();
539  FinishEvents();
540 }
541 
542 
543 
544 ClassImp(TtCracowTask);
x0
Double_t x0
Definition: checkhelixhit.C:70
fVerbose
int fVerbose
Definition: poormantracks.C:24
d
TObjArray * d
Definition: Pnd_Hc_etaee7G.C:78
i
Int_t i
Definition: run_full.C:25
TtCracowTask::fEloss
Double_t fEloss[6]
Definition: TtCracowTask.h:58
par
Double_t par[3]
Definition: anaLmdCluster.C:146
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)<< "
PndSdsHit::GetEloss
Double_t GetEloss() const
Definition: PndSdsHit.h:97
jj
Int_t jj
Definition: createRootGeoFileBwEnd_FwEndCapGeo.C:434
CAMath::Abs
static T Abs(const T &x)
Definition: PndCAMath.h:39
TtCracowTask::fTrackcount
Int_t fTrackcount
Definition: TtCracowTask.h:54
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
TtCracowTask::MyFit
void MyFit(Double_t *x, Double_t *y, Double_t *z, Double_t *Erx, Double_t *Ery, Double_t *Erz, Double_t *par, Double_t &chiX, Double_t &chiY)
Definition: TtCracowTask.cxx:478
TtCracowTask::~TtCracowTask
virtual ~TtCracowTask()
Definition: TtCracowTask.cxx:55
TtCracowTask::fTrackArray
TClonesArray * fTrackArray
Definition: TtCracowTask.h:50
z
Double_t z
Definition: createRootFscGeometryFile.C:146
min
friend F32vec4 min(const F32vec4 &a, const F32vec4 &b)
Definition: P4_F32vec4.h:25
name
TString name
Definition: createRootGeoFile.C:82
TtCracowTask::fTCandBranchName
TString fTCandBranchName
Definition: TtCracowTask.h:47
x
Double_t x
Definition: createRootFscGeometryFile.C:159
TtCracowTask::distance2
static double distance2(double x, double y, double z, double *p)
Definition: TtCracowTask.cxx:102
ClassImp
ClassImp(PndAnaContFact)
TtCracowTask::SumDistance2
static void SumDistance2(int &, double *, double &sum, double *par, int)
Definition: TtCracowTask.cxx:140
y
Double_t y
Definition: createRootFscGeometryFile.C:159
TtCracowTask::fTCandArray
TClonesArray * fTCandArray
Definition: TtCracowTask.h:46
PndSdsHit::GetSensorID
Int_t GetSensorID() const
Definition: PndSdsHit.h:90
TtCracowTask::Exec
virtual void Exec(Option_t *opt)
Definition: TtCracowTask.cxx:184
TtCracowTask::FinishEvent
virtual void FinishEvent()
Definition: TtCracowTask.cxx:536
TtCracowTask::distance2Single
static double distance2Single(double x, double y, double z, double ex, double ey, double ez, double *p)
Definition: TtCracowTask.cxx:116
TtCracowTask::Init
virtual InitStatus Init()
Definition: TtCracowTask.cxx:60
Generated on Wed Apr 3 2019 08:44:14 for FairRoot/PandaRoot by   doxygen 1.8.5