FairRoot/PandaRoot
PndSttHelixTrackFitter.cxx
Go to the documentation of this file.
1 //
2 
3 #include "PndSttHelixTrackFitter.h"
4 
5 #include "PndSttTrackMatch.h"
6 #include "PndSttHoughDefines.h"
7 #include "PndSttGeomPoint.h"
8 #include "PndSttHit.h"
9 #include "PndSttPoint.h"
10 #include "PndSttTube.h"
11 
12 #include "PndTrackCandHit.h"
13 #include "PndTrackCand.h"
14 
15 #include "FairRootManager.h"
16 #include "FairTask.h"
17 
18 #include "TArc.h"
19 #include "TH2.h"
20 #include "TClonesArray.h"
21 
22 #include "TMatrixD.h"
23 #include "TMarker.h"
24 #include "TLine.h"
25 #include "TPolyLine.h"
26 #include "TMinuit.h"
27 #include "TMath.h"
28 
29 #include <string>
30 #include <sstream>
31 #include <iostream>
32 //#include <cmath>
33 
34 using std::cout;
35 using std::cin;
36 using std::endl;
37 using std::map;
38 using TMath::ATan2;
39 using TMath::ATan;
40 using TMath::Cos;
41 using TMath::Sin;
42 using TMath::Sqrt;
43 
44 Int_t h = 0; // negative charge = 1; positive charge = -1
45 Double_t vote[201][1001], votecon[201];
46 TArrayD marray(200);
47 
48 PndSttHelixTrackFitter::PndSttHelixTrackFitter()
49 {
50  fEventCounter = 0;
51  fVerbose = 1;
52  fConstraint = 0;
53  fDisplayLevel = 0;
54 }
55 
56 PndSttHelixTrackFitter::PndSttHelixTrackFitter(Int_t verbose)
57 {
58  fVerbose = verbose;
59  fEventCounter = 0;
60  fConstraint = 0;
61  fDisplayLevel = 0;
62 }
63 
64 PndSttHelixTrackFitter::~PndSttHelixTrackFitter()
65 {
66 }
67 
68 void PndSttHelixTrackFitter::Init()
69 {
70  fEventCounter = 0;
71 
72  // Get and check FairRootManager
73  FairRootManager
74  *ioman = FairRootManager::Instance();
75 
76  if (! ioman)
77  {
78  cout << "-E- PndSttHelixTrackFitter::Init: "
79  << "RootManager not instantised!" << endl;
80  // return kFATAL;
81  }
82 
83  // Get hit Array
84  fHitArray = (TClonesArray*) ioman->GetObject("STTHit");
85  if ( ! fHitArray)
86  {
87  cout << "-W- PndSttHelixTrackFitter::Init: No Hit array!"
88  << endl;
89  }
90 
91  // Get point Array
92  fPointArray = (TClonesArray*) ioman->GetObject("STTPoint");
93  if ( ! fPointArray)
94  {
95  cout << "-W- PndSttHelixTrackFitter::Init: No Point array!"
96  << endl;
97  }
98 
99 
100  if( fDisplayLevel > 0) InitEventDisplay();
101 
102 }
103 
104 
105 Int_t PndSttHelixTrackFitter::DoFit(PndTrackCand* pTrackCand, PndSttTrack* pTrack, Int_t pidHypo)
106 {
107 
108  if(fDisplayLevel > 0) RunEventDisplay(pTrackCand);
109 
110 
111  if(fConstraint == 0) return DoFitPlain(pTrackCand, pTrack, pidHypo);
112  else if (fConstraint == 1) {
113  cout << "PndSttHelixTrackFitter::DoFit, Constraint 1 not usable now. A big change is needed to use start point from PndTrack and not PndTrackCand anymore" << endl;
114  return 1;
115  // return DoFitThroughOrigin(pTrackCand, pTrack, pidHypo);
116  }
117  else {
118  cout << "PndSttHelixTrackFitter::DoFit, Constraint " << fConstraint << " not implemented" << endl;
119  return 1;
120  }
121 }
122 
123 // ===================================================================================
124 // CONSTRAINT 0
125 Int_t PndSttHelixTrackFitter::DoFitPlain(PndTrackCand* pTrackCand, PndSttTrack* pTrack, Int_t )//pidHypo //[R.K.03/2017] unused variable(s)
126 {
127  // cout << "track fitting event # " << fEventCounter << endl;
128 
129  if(!pTrackCand) return 0;
130  fTrack = pTrack; // CHECK canc
131  fTrackCand = pTrackCand;
132  if(fDisplayLevel > 0) {
133  RunEventDisplay(pTrackCand);
134  char goOnChar;
135  cout << "press any key to continue: " << endl;
136  cin >> goOnChar;
137  }
138 
139  Int_t fit = 0;
140 
141  fit = XYFit(pTrackCand, 1);
142 
143  if(fit == 0 || fTrack->GetRad() == 0 || !(fTrack->GetRad()) || fTrack->GetRad() > 3000) {
144  if(fVerbose == 2) cout << "-E- pre prefit FAILED " << fit << " " << fTrack->GetRad() << endl;
145  fTrack->SetRad(-999);
146  return 0;
147  }
148 
149  // cout << "prefit-------------------------------" << endl;
150  fit = 0;
151  fit = MinuitFit(pTrackCand, 1);
152 
153  // if the fit fails
154  if(fit == 0 || fTrack->GetRad() == 0 || !(fTrack->GetRad()) || fTrack->GetRad() > 3000) {
155  fTrack->SetRad(-999);
156  if(fVerbose == 2) cout << "-E- prefit FAILED" << endl;
157  return 0;
158  }
159  else {
160  pTrack->SetFlag(1); // prefit done
161  Bool_t Rint = IntersectionFinder(pTrackCand);
162  // cout << "refit" << endl;
163  // if refit is OK
164  if(Rint == kTRUE) {
165  fit = XYFit(pTrackCand, 2);
166  Rint = IntersectionFinder(pTrackCand);
167  if(Rint == kTRUE) fit = MinuitFit(pTrackCand, 2);
168  Rint = IntersectionFinder(pTrackCand);
169  if(Rint == kTRUE) fit = MinuitFit(pTrackCand, 2);
170 
171  }
172  else {
173  return 0;
174  }
175  if(fit == 1 && fTrack->GetRad()>0&& fTrack->GetRad() < 3000) {
176 
177  pTrack->SetFlag(2); // refit done
178  Bool_t zint = ZFinder(pTrackCand, 1);
179 
180  if(zint == kTRUE) {
181  Int_t zfit = ZFit(pTrackCand, 1);
182  if(zfit == 1) {
183  pTrack->SetFlag(3); // z fit done
184  }
185  }
186  else if(fVerbose == 2) cout << "-E- zfinder FAILED" << endl;
187  }
188  }
189 
190  if(fVerbose == 2) {
191  cout << "param last d : " << fTrack->GetDist() << endl;
192  cout << "param last phi : " << fTrack->GetPhi() << endl;
193  cout << "param last R : " << fTrack->GetRad() << endl;
194  cout << "param last tanL : " << fTrack->GetTanL() << endl;
195  cout << "param last q : " << fTrack->GetCharge() << endl;
196  double pt, pl;
197  pt = fTrack->GetRad() * 0.006;
198  pl = fTrack->GetRad() * fTrack->GetTanL() * 0.006;
199  cout << "pT : " << pt << endl;
200  cout << "pL : " << pl << endl;
201  cout << "px, py, pz : " << pt * cos(-h * TMath::Pi()/2. + fTrack->GetPhi()) << " " << pt * sin(-h * TMath::Pi()/2. + fTrack->GetPhi()) << " " << pl << endl;
202  }
203 
204  if( fDisplayLevel > 0) FinishEventDisplay(fTrack);
205 
206  return 0;
207 }
208 
209 // XYFit was Fit4b
210 Int_t PndSttHelixTrackFitter::XYFit(PndTrackCand* pTrackCand, Int_t whatToFit) {
211 
212  if(!pTrackCand) return 0;
213 
214  Int_t hitcounter = pTrackCand->GetNHits();
215  Bool_t first = kFALSE;
216  if(hitcounter == 0) return 0;
217  if(hitcounter < 5) { // hitcounter > 50
218  if(fVerbose == 2) cout << "Bad No of hits in STT " << hitcounter << endl;
219  return 0;
220  }
221  if(fVerbose == 2) cout << "FIT xy ********************" << endl;
222 
223  // traslation and rotation -----------
224  // traslation near the first point
225  Double_t s = 0.001;
226  Double_t trasl[2];
227  // rotation
228  Double_t alpha;
229  if(fVerbose == 2) cout << "hitcounter: " << hitcounter << endl;
230  for(Int_t k = 0; k <hitcounter; k++) {
231 
232  PndTrackCandHit candhit = pTrackCand->GetSortedHit(k);
233  Int_t iHit = candhit.GetHitId();
234  PndSttHit *currenthit = (PndSttHit*) fHitArray->At(iHit);
235  if(!currenthit) { cout << "PndSttHelixTrackFitter::XYFit: no hit at " << iHit << endl; continue; }
236 
237  if(currenthit->GetX() == -999 || currenthit->GetY() == -999) continue;
238  Int_t refindex = currenthit->GetRefIndex();
239  // get point
240  PndSttPoint *iPoint = (PndSttPoint*) fPointArray->At(refindex);
241  if(!iPoint) { cout << "PndSttHelixTrackFitter::XYFit: no point at " << refindex << " associated to hit " << iHit << endl; continue; }
242 
243  // tubeID CHECK added
244  Int_t tubeID = currenthit->GetTubeID();
245  PndSttTube *tube = (PndSttTube*) fTubeArray->At(tubeID);
246 
247  PndSttHit *hitfirst, *hitlast;
248 
249  TVector3 wiredirection = tube->GetWireDirection();
250  if(wiredirection != TVector3(0.,0.,1.)) continue;
251  else if(first == kFALSE)
252  {
253  hitfirst = (PndSttHit*) fHitArray->At(iHit);
254  first = kTRUE;
255  trasl[0] = hitfirst->GetX();
256  trasl[1] = hitfirst->GetY();
257  }
258  else{
259  hitlast = (PndSttHit*) fHitArray->At(iHit);
260  alpha = TMath::ATan2(hitlast->GetY() - hitfirst->GetY(), hitlast->GetX() - hitfirst->GetX());
261  }
262  }
263 
264  // error <--> resolution
265  Double_t sigr, sigx, sigy; // = 0.14; //sigxy, //[R.K. 01/2017] unused variable?
266 
267  // FITTING IN X-Y PLANE:
268  // v = a + bu + cu^2
269 
270 
271  Double_t Suu, Su, Sv, Suv, S1, Suuu, Suuv, Suuuu;
272 
273  Su = 0.;
274  Sv = 0.;
275  Suu = 0.;
276  Suv = 0.;
277  Suuu = 0.;
278  S1 = 0.;
279  Suuv = 0.;
280  Suuuu = 0.;
281 
282  Int_t digicounter = 0;
283  TArrayD uarray(hitcounter);
284  TArrayD varray(hitcounter);
285  TArrayD sigv2array(hitcounter);
286  TArrayD sigu2array(hitcounter);
287 
288  for(Int_t i=0; i < hitcounter; i++){
289  uarray.AddAt(-999, i);
290  varray.AddAt(-999, i);
291  sigv2array.AddAt(-999, i);
292  sigu2array.AddAt(-999, i);
293  }
294  for(Int_t i=0; i < hitcounter; i++){
295 
296  PndTrackCandHit candhit = pTrackCand->GetSortedHit(i);
297  Int_t iHit = candhit.GetHitId();
298  PndSttHit *currenthit = (PndSttHit*) fHitArray->At(iHit);
299  if(!currenthit) { cout << "PndSttHelixTrackFitter::XYFit: no hit at " << iHit << endl; continue; }
300  if(currenthit->GetX() == -999 || currenthit->GetY() == -999) continue;
301  Int_t refindex = currenthit->GetRefIndex();
302  // get point
303  PndSttPoint *iPoint = (PndSttPoint*) fPointArray->At(refindex);
304  if(!iPoint) { cout << "PndSttHelixTrackFitter::XYFit: no point at " << refindex << " associated to hit " << iHit << endl; continue; }
305 
306  // tubeID CHECK added
307  Int_t tubeID = currenthit->GetTubeID();
308  PndSttTube *tube = (PndSttTube*) fTubeArray->At(tubeID);
309 
310  TVector3 wiredirection = tube->GetWireDirection();
311  if(wiredirection != TVector3(0.,0.,1.)) continue;
312 
313 
314  if(whatToFit == 2) {
315  //Double_t resx = iPoint->GetX() - currenthit->GetX(); //[R.K. 01/2017] unused variable?
316  //Double_t resy = iPoint->GetY() - currenthit->GetY(); //[R.K. 01/2017] unused variable?
317  //Double_t resdist = TMath::Sqrt((iPoint->GetY() - currenthit->GetY())*(iPoint->GetY() - currenthit->GetY()) + (iPoint->GetX() - currenthit->GetX())*(iPoint->GetX() - currenthit->GetX())); //[R.K. 01/2017] unused variable?
318 
319  }
320 
321  Double_t xtrasl, ytrasl;
322  // traslation
323  xtrasl = currenthit->GetX() - trasl[0];
324  ytrasl = currenthit->GetY() - trasl[1];
325 
326  Double_t xrot, yrot;
327  // rotation
328  xrot = TMath::Cos(alpha)*xtrasl + TMath::Sin(alpha)*ytrasl;
329  yrot = -TMath::Sin(alpha)*xtrasl + TMath::Cos(alpha)*ytrasl;
330 
331  // re-traslation
332  xtrasl = xrot + s;
333  ytrasl = yrot;
334 
335  if( fDisplayLevel >= 3) {
336  eventCanvas->cd(1);
337  TMarker *pt = new TMarker(xrot, yrot, 6);
338  pt->SetMarkerColor(3);
339  // if(whatToFit == 2)
340  // pt->Draw("SAME");
341  }
342 
343  // change coordinate
344  Double_t u, v, sigv2, sigu2;
345  u = xtrasl / (xtrasl*xtrasl + ytrasl*ytrasl);
346  v = ytrasl / (xtrasl*xtrasl + ytrasl*ytrasl);
347 
348  if( fDisplayLevel >= 4) {
349  eventDetails->cd(1);
350  TMarker *uv = new TMarker(u, v, 6);
351  // uv->SetMarkerColor(3);
352  if(whatToFit == 2) uv->Draw("SAME");
353  eventDetails->Update();
354  eventDetails->Modified();
355  }
356 
357  if(whatToFit == 1) {
358  if(currenthit->GetIsochrone() == 0) sigr = sqrt(2.)/sqrt(12.);
359  else sigr = sqrt(2.) * currenthit->GetIsochrone()/sqrt(12.);
360  sigx = sigr;
361  sigy = sigr;
362  }
363  else {
364  if(currenthit->GetIsochrone() == 0) {
365  sigr = sqrt(2.)/sqrt(12.);
366  sigx = sigr;
367  sigy = sigr;
368  }
369  else {
370  sigr = 0.0150;
371  sigx = fabs(sigr * TMath::Cos(TMath::ATan(marray.At(i))));
372  sigy = fabs(sigr * TMath::Sin(TMath::ATan(marray.At(i))));
373  //sigxy = sigr * TMath::Sqrt(fabs(TMath::Cos(TMath::ATan(marray.At(i))) * TMath::Sin(TMath::ATan(marray.At(i)))));; //[R.K. 01/2017] unused variable?
374  }
375  }
376 
377  Double_t dvdx = (-2 * xtrasl * ytrasl)/pow((xtrasl*xtrasl + ytrasl*ytrasl),2);
378  Double_t dvdy = (xtrasl*xtrasl - ytrasl*ytrasl) / pow((xtrasl*xtrasl + ytrasl*ytrasl),2);
379  Double_t dudx = (ytrasl*ytrasl - xtrasl*xtrasl) / pow((xtrasl*xtrasl + ytrasl*ytrasl),2);
380  Double_t dudy = (-2 * xtrasl * ytrasl)/pow((xtrasl*xtrasl + ytrasl*ytrasl),2);
381 
382  sigu2 = dudx * dudx * sigx * sigx + dudy * dudy * sigy * sigy + 2 * dudx * dudy * sigx * sigy;
383  sigv2 = dvdx * dvdx * sigx * sigx + dvdy * dvdy * sigy * sigy + 2 * dvdx * dvdy * sigx * sigy;
384 
385  uarray.AddAt(u, digicounter);
386  varray.AddAt(v, digicounter);
387  sigu2array.AddAt(sigu2, digicounter);
388  sigv2array.AddAt(sigv2, digicounter);
389 
390  Su = Su + (u/sigv2);
391  Sv = Sv + (v/sigv2);
392 
393  Suv = Suv + ((u*v)/sigv2);
394  Suu = Suu + ((u*u)/sigv2);
395 
396  Suuu = Suuu + ((u*u*u)/sigv2);
397  Suuv = Suuv + ((u*u*v)/sigv2);
398 
399  Suuuu = Suuuu + ((u*u*u*u)/sigv2);
400 
401  S1 = S1 + 1/sigv2;
402 
403  digicounter++;
404  }
405 
406  TMatrixD matrix(3,3);
407  matrix[0][0] = S1;
408  matrix[0][1] = Su;
409  matrix[0][2] = Suu;
410 
411  matrix[1][0] = Su;
412  matrix[1][1] = Suu;
413  matrix[1][2] = Suuu;
414 
415  matrix[2][0] = Suu;
416  matrix[2][1] = Suuu;
417  matrix[2][2] = Suuuu;
418 
419  Double_t determ;
420 
421  determ = matrix.Determinant();
422 
423  if (determ != 0) {
424  matrix.Invert();
425  }
426  else {
427  return 0;
428  if(fVerbose == 2) cout << "DET 0" << endl;
429  }
430 
431  TMatrixD column(3,1);
432  column[0][0] = Sv;
433  column[1][0] = Suv;
434  column[2][0] = Suuv;
435 
436  TMatrixD column2(3,1);
437  column2.Mult(matrix, column);
438 
439  Double_t a, b, c;
440  a = column2[0][0];
441  b = column2[1][0];
442  c = column2[2][0];
443 
444  if(fVerbose == 2) {
445  std::cout << "1) parabolic parameters:\n";
446  std::cout << "a = " << column2[0][0] << "\n";
447  std::cout << "b = " << column2[1][0] << "\n";
448  std::cout << "c = " << column2[2][0] << "\n";
449  }
450  if(fabs(a)<0.000001) return 0;
451 
452  Double_t chi2;
453  chi2 = 0.;
454  for(Int_t i=0; i < digicounter; i++){
455  if(uarray.At(i) == -999 || varray.At(i) == -999 || sigv2array.At(i) == -999) continue;
456  chi2 = chi2 + pow(((varray.At(i) - (a + b*uarray.At(i) + c*uarray.At(i)*uarray.At(i))) /sqrt(sigv2array.At(i))),2) ;
457  }
458 
459  if(fVerbose == 2) cout << "digicounter: " << digicounter << endl;
460 
461  //------------------ with right errors
462  Su = 0.;
463  Sv = 0.;
464  Suu = 0.;
465  Suv = 0.;
466  Suuu = 0.;
467  S1 = 0.;
468  Suuv = 0.;
469  Suuuu = 0.;
470 
471  TArrayD sigE2array(digicounter);
472  Double_t sigE2;
473  for(Int_t i=0; i< digicounter; i++){
474  if(uarray.At(i) == -999 || varray.At(i) == -999 || sigv2array.At(i) == -999) continue;
475  sigE2 = sigv2array.At(i) + sigu2array.At(i) * pow((b + 2*c*uarray.At(i)),2);
476  sigE2array.AddAt(sigE2, i);
477 
478  Su = Su + (uarray.At(i)/sigE2);
479  Sv = Sv + (varray.At(i)/sigE2);
480 
481  Suv = Suv + ((uarray.At(i)*varray.At(i))/sigE2);
482  Suu = Suu + ((uarray.At(i)*uarray.At(i))/sigE2);
483 
484  Suuu = Suuu + ((uarray.At(i)*uarray.At(i)*uarray.At(i))/sigE2);
485  Suuv = Suuv + ((uarray.At(i)*uarray.At(i)*varray.At(i))/sigE2);
486 
487  Suuuu = Suuuu + ((uarray.At(i)*uarray.At(i)*uarray.At(i)*uarray.At(i))/sigE2);
488 
489  S1 = S1 + 1/sigE2;
490  }
491 
492  TMatrixD matrixb(3,3);
493  matrixb[0][0] = S1;
494  matrixb[0][1] = Su;
495  matrixb[0][2] = Suu;
496 
497  matrixb[1][0] = Su;
498  matrixb[1][1] = Suu;
499  matrixb[1][2] = Suuu;
500 
501  matrixb[2][0] = Suu;
502  matrixb[2][1] = Suuu;
503  matrixb[2][2] = Suuuu;
504 
505  determ = matrixb.Determinant();
506 
507  if (determ != 0) {
508  matrixb.Invert();
509  }
510  else {
511  return 0;
512  if(fVerbose == 2) cout << "DET 0" << endl;
513  }
514 
515  TMatrixD columnb(3,1);
516  columnb[0][0] = Sv;
517  columnb[1][0] = Suv;
518  columnb[2][0] = Suuv;
519 
520  TMatrixD column2b(3,1);
521  column2b.Mult(matrixb, columnb);
522 
523  a = column2b[0][0];
524  b = column2b[1][0];
525  c = column2b[2][0];
526 
527  // std::cout << "*************\n";
528  // std::cout << "2) parabolic parameters:\n";
529  // std::cout << "a = " << column2b[0][0] << "\n";
530  // std::cout << "b = " << column2b[1][0] << "\n";
531  // std::cout << "c = " << column2b[2][0] << "\n";
532 
533  //v = a + bu + cu^2
534  chi2 = 0.;
535  for(Int_t i=0; i<digicounter; i++){
536  if(uarray.At(i) == -999 || varray.At(i) == -999 || sigv2array.At(i) == -999) continue;
537  chi2 = chi2 + pow(((varray.At(i) - (a + b*uarray.At(i) + c*uarray.At(i)*uarray.At(i)))/sqrt(sigE2array.At(i))), 2);
538  }
539 
540  if( fDisplayLevel >= 4) {
541  eventDetails->cd(1);
542  Double_t uu[100];
543  Double_t vv[100];
544  for(Int_t p = 0; p<100; p++){
545  uu[p] = -1.5 + p*3*0.01;
546  vv[p] = a + b*uu[p] + c*uu[p]*uu[p];
547  }
548  TPolyLine *uvline = new TPolyLine(100, uu, vv);
549  if(whatToFit == 2) uvline->Draw("SAME");
550  eventDetails->Update();
551  eventDetails->Modified();
552  }
553 
554  // center and radius
555  Double_t xcrot, ycrot, xc, yc, epsilon, R;
556  ycrot = 1/(2*a);
557  xcrot = -b/(2*a);
558  epsilon = -c*pow((1+(b*b)), -3/2);
559  R = epsilon + sqrt((xcrot*xcrot)+(ycrot*ycrot));
560 
561  // // errors on parameters -------------------
562  // // a, b, c
563  // Double_t Da[MAXNUMSTTHITS], Db[MAXNUMSTTHITS], Dc[MAXNUMSTTHITS];
564  // Double_t p1[MAXNUMSTTHITS], pu[MAXNUMSTTHITS], puu[MAXNUMSTTHITS];
565  // Double_t erra2, errb2, errc2;
566  // erra2 = 0.;
567  // errb2 = 0.;
568  // errc2 = 0.;
569 
570  // for(Int_t i=0; i<fTrackNumSttHits[fitTrack]; i++){
571  // p1[i] = 1/sigE2[i];
572  // pu[i] = u[i]*p1[i];
573  // puu[i] = u[i]*pu[i];
574 
575  // Da[i] = matrixb[0][0] * p1[i] + matrixb[0][1] * pu[i] + matrixb[0][2] * puu[i];
576  // Db[i] = matrixb[1][0] * p1[i] + matrixb[1][1] * pu[i] + matrixb[1][2] * puu[i];
577  // Dc[i] = matrixb[2][0] * p1[i] + matrixb[2][1] * pu[i] + matrixb[2][2] * puu[i];
578 
579  // erra2 = erra2 + (Da[i]*Da[i]*sigE2[i]);
580  // errb2 = errb2 + (Db[i]*Db[i]*sigE2[i]);
581  // errc2 = errc2 + (Dc[i]*Dc[i]*sigE2[i]);
582  // }
583 
584  // // std::cout << "**************\n";
585  // // std::cout << "erra = " << sqrt(erra2) << "\n";
586  // // std::cout << "errb = " << sqrt(errb2)<< "\n";
587  // // std::cout << "errc = " << sqrt(errc2)<< "\n";
588 
589  // // errors on xc, yc, R
590  // Double_t errxcrot2, errycrot2, errR2;
591  // errxcrot2 = (1./(4*pow(a,4))) * (b*b*erra2 + a*a*errb2 - 2*a*b*sqrt(erra2*errb2)) ;
592  // errycrot2 = (1./4)*(erra2/pow(a,4)) ;
593  // errR2 = (1./(R*R)) * (ycrot*ycrot*errxcrot2 + xcrot*xcrot*errycrot2 + 2*xcrot*ycrot*sqrt(errxcrot2*errycrot2)) ;
594 
595  // re-rotation and re-traslation of xc and yc
596  // translation
597  xcrot = xcrot - s;
598  // rotation
599  xc = TMath::Cos(alpha)*xcrot - TMath::Sin(alpha)*ycrot;
600  yc = TMath::Sin(alpha)*xcrot + TMath::Cos(alpha)*ycrot;
601  // traslation
602  xc = xc + trasl[0];
603  yc = yc + trasl[1];
604  Double_t phi = TMath::ATan2(yc, xc);
605  Double_t d;
606  d = ((xc + yc) - R*(TMath::Cos(phi) + TMath::Sin(phi)))/(TMath::Cos(phi) + TMath::Sin(phi));
607 
608  fTrack->SetDist(d);
609  fTrack->SetPhi(phi);
610  // Double_t newZ = -999.; // CHECK da cambiare
611  // fTrack->SetZ(newZ);
612  fTrack->SetRad(R);
613  // Double_t newTheta = -999.; // CHECK da cambiare
614  // fTrack->SetTanL(newTheta);
615  h = -GetCharge(d, phi, R);
616  fTrack->SetCharge(-h);
617 
618  // cout << "FIT: " << xc << " " << yc << endl;
619  // cout << "RAGGIO: " << R << endl;
620 
621  if( fDisplayLevel >= 3) {
622  eventCanvas->cd(1);
623  TArc *fitarc = new TArc(((fTrack->GetDist() + fTrack->GetRad()) * cos(fTrack->GetPhi())), ((fTrack->GetDist() + fTrack->GetRad()) * sin(fTrack->GetPhi())), fTrack->GetRad());
624  if(whatToFit == 2) fitarc->SetLineColor(kRed - 9);
625  fitarc->SetFillStyle(0);
626  fitarc->Draw("SAME");
627  eventCanvas->Update();
628  eventCanvas->Modified();
629  }
630  return 1;
631 }
632 
633 
634 // -------------- IntersectionFinder --------------------------------------
635 Bool_t PndSttHelixTrackFitter::IntersectionFinder(PndTrackCand *pTrackCand)
636 {
637 
638  ResetMArray();
639 
640  // calculation of the curvature from the helix prefit
641  if(fTrack->GetRad() == 0) return kFALSE;
642 
643  // vec = (xc, yc)
644  TVector2 vec((fTrack->GetDist() + fTrack->GetRad()) * cos(fTrack->GetPhi()), (fTrack->GetDist() + fTrack->GetRad()) * sin(fTrack->GetPhi()));
645 
646  //==========
647  // POINT ----------------------------------------------------
648  // 1. find the cooordinates of the point fired wire of the track
649  TVector2 point; // point
650  Double_t radius;
651 
652  Int_t counter = 0;
653  // loop over input points
654  Int_t hitcounter = pTrackCand->GetNHits();
655  for(Int_t k = 0; k < hitcounter; k++){
656  // get hit
657  PndTrackCandHit candhit = pTrackCand->GetSortedHit(k);
658  Int_t iHit = candhit.GetHitId();
659 
660  PndSttHit *pMhit = (PndSttHit*) fHitArray->At(iHit);
661  if (!pMhit ) { cout << "PndSttHelixTrackFitter::IntersectionFinder: no hit at " << iHit << endl; continue; }
662 
663  Int_t refindex = pMhit->GetRefIndex();
664  // get point
665  PndSttPoint *iPoint = (PndSttPoint*) fPointArray->At(refindex);
666  if (!iPoint ) { cout << "PndSttHelixTrackFitter::IntersectionFinder: no point at " << refindex << " associated to hit " << iHit << endl; continue; }
667 
668 
669  // tubeID CHECK added
670  Int_t tubeID = pMhit->GetTubeID();
671  PndSttTube *tube = (PndSttTube*) fTubeArray->At(tubeID);
672  TVector3 wiredirection = tube->GetWireDirection();
673 
674  if(wiredirection != TVector3(0.,0.,1.)) continue;
675 
676  // [xp, yp] point = coordinates xy of the centre of the firing tube
677  point.Set(tube->GetPosition().X(), tube->GetPosition().Y());
678  radius = pMhit->GetIsochrone();
679 
680  // the coordinates of the point are taken from the intersection
681  // between the circumference from the drift time and the R radius of
682  // curvature. -------------------------------------------------------
683  // 2. find the intersection between the little circle and the line // R
684  // 2.a
685  // find the line passing throught [xc, yc] (centre of curvature) and [xp, yp] (first wire)
686  // y = mx + q
687  Double_t m = (point.Y() - vec.Y())/(point.X() - vec.X());
688  Double_t q = point.Y() - m*point.X();
689 
690  if( fDisplayLevel >= 3) {
691  eventCanvas->cd(1);
692 // TArc *archetto = new TArc(point.X(), point.Y(), radius);
693 // archetto->SetFillStyle(0);// archetto->Draw("SAME");
694  TLine *line = new TLine(-50, -50*m + q, 50, 50*m + q);
695  line->SetLineColor(5);
696  // line->Draw("SAME");
697  eventCanvas->Update();
698  eventCanvas->Modified();
699  }
700 
701  // cut on radius CHECK
702  // if the simulated radius is too small, the pMhit
703  // is not used for the fit
704  if(radius < 0.1) {
705  marray.AddAt(-999, k);
706  pMhit->SetX(-999);
707  pMhit->SetY(-999);
708  continue; // CHECK
709  }
710 
711  Double_t x1 = 0, y1 = 0,
712  x2 = 0, y2 = 0,
713  xb1 = 0, yb1 = 0,
714  xb2 = 0, yb2 = 0;
715 
716  // CHECK the vertical track
717  if(fabs(point.X() - vec.X()) < 1e-6) {
718 
719  // 2.b
720  // intersection little circle and line --> [x1, y1]
721  // + and - refer to the 2 possible intersections
722  // +
723  x1 = point.X();
724  y1 = point.Y() + sqrt(radius * radius - (x1 - point.X()) * (x1 - point.X()));
725  // -
726  x2 = x1;
727  y2 = point.Y() - sqrt(radius * radius - (x2 - point.X()) * (x2 - point.X()));
728 
729  // 2.c intersection between line and circle
730  // +
731  xb1 = vec.X();
732  yb1 = vec.Y() + sqrt(fTrack->GetRad() * fTrack->GetRad() - (xb1 - vec.X()) * (xb1 - vec.X()));
733  // -
734  xb2 = xb1;
735  yb2 = vec.Y() - sqrt(fTrack->GetRad() * fTrack->GetRad() - (xb2 - vec.X()) * (xb2 - vec.X()));
736 
737  } // END CHECK
738  else {
739 
740  // 2.b
741  // intersection little circle and line --> [x1, y1]
742  // + and - refer to the 2 possible intersections
743  // +
744  x1 = (-(m*(q - point.Y()) - point.X()) + sqrt((m*(q - point.Y()) - point.X())*(m*(q - point.Y()) - point.X()) - (m*m + 1)*((q - point.Y())*(q - point.Y()) + point.X()*point.X() - radius*radius))) / (m*m + 1);
745  y1 = m*x1 + q;
746  // -
747  x2 = (-(m*(q - point.Y()) - point.X()) - sqrt((m*(q - point.Y()) - point.X())*(m*(q - point.Y()) - point.X()) - (m*m + 1)*((q - point.Y())*(q - point.Y()) + point.X()*point.X() - radius*radius))) / (m*m + 1);
748  y2 = m*x2 + q;
749 
750  // 2.c intersection between line and circle
751  // +
752  xb1 = (-(m*(q - vec.Y()) - vec.X()) + sqrt((m*(q - vec.Y()) - vec.X())*(m*(q - vec.Y()) - vec.X()) - (m*m + 1)*((q - vec.Y())*(q - vec.Y()) + vec.X()*vec.X() - (fTrack->GetRad()) *(fTrack->GetRad()) ))) / (m*m + 1);
753  yb1 = m*xb1 + q;
754  // -
755  xb2 = (-(m*(q - vec.Y()) - vec.X()) - sqrt((m*(q - vec.Y()) - vec.X())*(m*(q - vec.Y()) - vec.X()) - (m*m + 1)*((q - vec.Y())*(q - vec.Y()) + vec.X()*vec.X() - (fTrack->GetRad()) *(fTrack->GetRad())))) / (m*m + 1);
756  yb2 = m*xb2 + q;
757  }
758 
759  // calculation of the distance between [xb, yb] and [xp, yp]
760  Double_t distb1 = sqrt((yb1 - point.Y())*(yb1 - point.Y()) + (xb1 - point.X())*(xb1 - point.X()));
761  Double_t distb2 = sqrt((yb2 - point.Y())*(yb2 - point.Y()) + (xb2 - point.X())*(xb2 - point.X()));
762 
763  // choice of [xb, yb]
764  TVector2 xyb;
765  if(distb1 > distb2) xyb.Set(xb2, yb2);
766  else xyb.Set(xb1, yb1);
767 
768  // calculation of the distance between [x, y] and [xb. yb]
769  Double_t dist1 = sqrt((xyb.Y() - y1)*(xyb.Y() - y1) + (xyb.X() - x1)*(xyb.X() - x1));
770  Double_t dist2 = sqrt((xyb.Y() - y2)*(xyb.Y() - y2) + (xyb.X() - x2)*(xyb.X() - x2));
771 
772  // choice of [x, y]
773  TVector2 *xy;
774  if(dist1 > dist2) xy = new TVector2(x2, y2);
775  else xy = new TVector2(x1, y1); // <========= THIS IS THE NEW POINT to be used for the fit
776 
777  // SET AS DEBUG
778  // if (TMath::Sqrt((xy->X() - point.X())*(xy->X() - point.X()) + (xy->Y() - point.Y())*(xy->Y() - point.Y())) - radius > 0.000001) {
779  // cout << "ATTENZIONE: " << "differenza = " << TMath::Sqrt((xy->X() - point.X())*(xy->X() - point.X()) + (xy->Y() - point.Y())*(xy->Y() - point.Y())) - radius << endl;
780  // }
781 
782  // cout << "x hit prima: " << tube->GetPosition().X() << " " << pMhit->GetX() << endl;
783  // cout << prima: " << pMhit->GetX() << " " << pMhit->GetY() << " " << pMhit->GetZ() << endl;
784 
785  marray.AddAt(m, k);
786  pMhit->SetX(xy->X());
787  pMhit->SetY(xy->Y());
788 
789  if( fDisplayLevel >= 3) {
790  eventCanvas->cd(1);
791  TMarker *np = new TMarker(pMhit->GetX(), pMhit->GetY(), 6);
792  np->SetMarkerColor(3);
793  // np->Draw("SAME");
794  eventCanvas->Update();
795  eventCanvas->Modified();
796  }
797  delete xy;
798  counter++;
799  }
800 
801  if(counter==0) return kFALSE;
802  else return kTRUE;
803 }
804 
805 // ZFinder was ZFinderbb3
806 // -------- ZFinder --------------------------------------------
807 Bool_t PndSttHelixTrackFitter::ZFinder(PndTrackCand* pTrackCand, Int_t ) { //whatToFit //[R.K.03/2017] unused variable(s)
808  // the z finding procedure uses the hough transform to find the line
809  // in the plane z - track length on which the correct points lie.
810 
811  if(fVerbose == 2) cout << "ZFINDER" << endl;
812 
813  if(!pTrackCand) return 0;
814 
815  Int_t hitcounter = pTrackCand->GetNHits();
816 
817  // cut on number of hits
818  // if(hitcounter > 50) {
819  // cout << "more than 50 hits " << hitcounter << endl;
820  // return 0;
821  // }
822  if(hitcounter < 5) {
823  if(fVerbose == 2) cout << "less than 5 hits" << endl;
824  return 0;
825  }
826 
827  ZPointsArray = new TObjArray(2*hitcounter);
828 
829  // ZFIT =======
830  if(hitcounter == 0) return kFALSE;
831 
832  //Double_t Sxx, Sx, Sz, Sxz, S1z; //[R.K.01/2017]unused variable?
833  //Double_t Detz = 0.; //[R.K. 01/2017] unused variable?
834  //Double_t fitm, fitp; //[R.K. 01/2017] unused variable?
835  Double_t sigz = 1.; // CHECK
836 
837  //Sx = 0.; //[R.K.01/2017]unused variable?
838  //Sz = 0.; //[R.K.01/2017]unused variable?
839  //Sxx = 0.; //[R.K.01/2017]unused variable?
840  //Sxz = 0.; //[R.K.01/2017]unused variable?
841  //S1z = 0.; //[R.K.01/2017]unused variable?
843 
844  TVector3 *tofit, *tofit2;
845 
846  // centre of curvature
847  Double_t x_0 = (fTrack->GetDist() + fTrack->GetRad()) * cos(fTrack->GetPhi());
848  Double_t y_0 = (fTrack->GetDist() + fTrack->GetRad()) * sin(fTrack->GetPhi());
849 
850  // radius of curvature
851  Double_t R = fTrack->GetRad();
852  Int_t wireOk = 0;
853  //Bool_t first = kTRUE; //[R.K.01/2017]unused variable?
854 
855  for (Int_t i = 0; i < hitcounter; i++) {
856 
857  // get index of hit
858  PndTrackCandHit candhit = pTrackCand->GetSortedHit(i);
859  Int_t iHit = candhit.GetHitId();
860 
861  // get hit
862  PndSttHit *pMhit = (PndSttHit*) fHitArray->At(iHit);
863  if (!pMhit ) { cout << "PndSttHelixTrackFitter::ZFinder: no hit at " << iHit << endl; continue; }
864 
865  Int_t refindex = pMhit->GetRefIndex();
866  // get point
867  PndSttPoint *iPoint = (PndSttPoint*) fPointArray->At(refindex);
868  if (!iPoint ) { cout << "PndSttHelixTrackFitter::ZFinder: no point at " << refindex << " associated to hit " << iHit << endl; continue; }
869 
870  // tubeID CHECK added
871  Int_t tubeID = pMhit->GetTubeID();
872  PndSttTube *tube = (PndSttTube*) fTubeArray->At(tubeID);
873 
874  TVector3 wiredirection = tube->GetWireDirection();
875  TVector3 wiredirection2;
876 
877  wiredirection2 = tube->GetHalfLength() * wiredirection;
878  TVector3 cenposition = tube->GetPosition();
879  // if(tube->GetPosition().Z() != 35) continue; // to throw away short tubes
880 
881  TVector3 min, max;
882  min = cenposition - wiredirection2;
883  max = cenposition + wiredirection2;
884 
885  // first extremity
886  Double_t x_1= min.X();
887  Double_t y_1= min.Y();
888  Double_t z_1= min.Z();
889 
890  // second extremity
891  Double_t x_2= max.X();
892  Double_t y_2= max.Y();
893  Double_t z_2= max.Z();
894 
895  Double_t rcur = pMhit->GetIsochrone();
896 
897  Double_t x1 = -9999.;
898  Double_t y1 = -9999.;
899  Double_t x2 = -9999.;
900  Double_t y2 = -9999.;
901 
902  // from xy plane fit
903  Double_t phi0 = fTrack->GetPhi();
904  Double_t d0 = fTrack->GetDist();
905  Double_t x0 = d0*TMath::Cos(phi0);
906  Double_t y0 = d0*TMath::Sin(phi0);
907  // in xy plane: angle of the PCA to the origin
908  // with respect to the curvature center
909  Double_t Phi0 = TMath::ATan2((y0 - y_0),(x0 - x_0));
910 
911  if(wiredirection != TVector3(0.,0.,1.))
912  {
913  Double_t a = -999;
914  Double_t b = -999;
915 
916  // intersection point between the reconstructed
917  // circumference and the line joining the centres
918  // of the reconstructed circle and the i_th drift circle
919  if(fabs(x_2-x_1)>0.0001) {
920  a =(y_2-y_1)/(x_2-x_1);
921  b =(y_1-a*x_1);
922  Double_t A = a*a+1;
923  Double_t B = x_0+a*y_0-a*b;
924  Double_t C = x_0*x_0+y_0*y_0+b*b-R*R-2*y_0*b;
925  if((B*B-A*C)>0) {
926  x1= (B+TMath::Sqrt(B*B-A*C))/A;
927  x2= (B-TMath::Sqrt(B*B-A*C))/A;
928  y1=a*x1+b;
929  y2=a*x2+b;
930  }
931  }
932  else if(fabs(y_2-y_1)>0.0001) {
933  Double_t A = 1;
934  Double_t B = y_0;
935  Double_t C = y_0*y_0 +(x_1-x_0)*(x_1-x_0) -R*R;
936 
937  if((B*B-A*C)>0) {
938  y1= (B+TMath::Sqrt(B*B-A*C))/A;
939  y2= (B-TMath::Sqrt(B*B-A*C))/A;
940  x1=x2=x_1;
941  }
942  }
943  else {
944  if(fVerbose == 2) cout << "-E- intersection point not found" << endl;
945  continue;
946  }
947 
948  //x1 and x2 are the 2 intersection points
949  Double_t d1=TMath::Sqrt((x1-cenposition.X())*(x1-cenposition.X())+
950  (y1-cenposition.Y())*(y1-cenposition.Y()));
951  Double_t d2=TMath::Sqrt((x2-cenposition.X())*(x2-cenposition.X())+
952  (y2-cenposition.Y())*(y2-cenposition.Y()));
953 
954  Double_t x_ = x1;
955  Double_t y_ = y1;
956 
957  // the intersection point nearest to the drift circle's centre is taken
958  if(d2<d1) {x_=x2;y_=y2;}
959 
960  if( fDisplayLevel >= 3) {
961  eventCanvas->cd(1);
962  TMarker *pt = new TMarker(x_, y_, 6);
963  pt->SetMarkerColor(6);
964  pt->Draw("SAME");
965 
966  // intersection lines
967  TLine *intline = new TLine(x_, y_, x_0, y_0);
968  intline->SetLineColor(5);
969  intline->Draw("SAME");
970  }
971 
972  // now we need to find the actual centre of the drift circle,
973  // by translating the drift circle until it becomes tangent
974  // to the reconstructed circle.
975  // Two solutions are possible (left rigth abiguity),
976  // they are both kept, only the following zed fit will discard the wrong ones.
977  // Using the parametric equation of the 3d-straigth line and taking the
978  // x points just obtained, the zed coordinate of the skewed tube centre is calculated.
979 
980  if(x_1 == x_2)
981  {
982  // if the skewed layer lies in yz plane
983  x1 = x_;
984  y1 = y_ + rcur;
985  x2 = x_;
986  y2 = y_ - rcur;
987  }
988  else {
989  //solving the equation to find out the centre of the tangent circle
990  Double_t A = a*a+1;
991  Double_t B = -(a*b-a*y_-x_);
992  Double_t C = x_*x_+ y_*y_+b*b-2*b*y_-rcur*rcur;
993  if((B*B-A*C)>0) {
994  x1= (B+TMath::Sqrt(B*B-A*C))/A;
995  x2= (B-TMath::Sqrt(B*B-A*C))/A;
996  y1=a*x1+b;
997  y2=a*x2+b;
998  }
999  else if((B*B-A*C)==0){
1000  x1= B/A;
1001  x2 = x1;
1002  y1=a*x1+b;
1003  y2=a*x2+b;
1004  }
1005  else {
1006  if(fVerbose == 2) cout << "NO WAY2" << endl;
1007  continue;
1008  }
1009  }
1010  if( fDisplayLevel >= 3) {
1011  eventCanvas->cd(1);
1012  TArc *drftarc = new TArc(x1, y1, rcur);
1013  drftarc->SetLineColor(7);
1014  drftarc->SetFillStyle(0);
1015  drftarc->Draw("SAME");
1016  TArc *drftarc2 = new TArc(x2, y2, rcur);
1017  drftarc2->SetLineColor(7);
1018  drftarc2->SetFillStyle(0);
1019  drftarc2->Draw("SAME");
1020  }
1021 
1022  d1=TMath::Sqrt((x1-cenposition.X())*(x1-cenposition.X())+(y1-cenposition.Y())*(y1-cenposition.Y()));
1023  d2=TMath::Sqrt((x2-cenposition.X())*(x2-cenposition.X())+(y2-cenposition.Y())*(y2-cenposition.Y()));
1024 
1025  Double_t xcen0=x1;
1026  Double_t xcen1=x2;
1027  Double_t ycen0=y1;
1028  Double_t ycen1=y2;
1029 
1030  if(d2<d1) { // z2 contains the points nearest (in x-y) to the initial centre of the skewed tube
1031  xcen0=x2;
1032  xcen1=x1;
1033  ycen0=y2;
1034  ycen1=y1;
1035 
1036  }
1037 
1038  // zed association
1039  Double_t t_, z_, t_bis, z_bis;
1040  if(fabs(x_2-x_1)<0.001) {
1041  t_ =(ycen0-y_1)/(y_2-y_1); // k= a_y*t + y_1 [t=1 y=y_2]
1042  z_ =(z_2-z_1)*t_ +z_1; // z= a_z*t + z_1 [t=1 z=z_2]
1043  t_bis =(ycen1-y_1)/(y_2-y_1); // from y_'s (the 2 solutions of the 2nd order equation)
1044  z_bis =(z_2-z_1)*t_bis +z_1; // and the 2 parametric equations the z coord. are obtained
1045  }
1046  else{
1047  t_ =(xcen0-x_1)/(x_2-x_1); // x= a_x*t + x_1 [t=1 x=x_2]
1048  z_ =(z_2-z_1)*t_ +z_1; // z= a_z*t + z_1 [t=1 z=z_2]
1049  t_bis =(xcen1-x_1)/(x_2-x_1); // from x_'s (the 2 solutions of the 2nd order equation)
1050  z_bis =(z_2-z_1)*t_bis +z_1; // and the 2 parametric equations the z coord. are obtained
1051  }
1052 
1053  // tofit = new TVector3(xcen0,ycen0,z_);
1054  tofit = new TVector3(x_,y_,z_);
1055  ZPointsArray->Add(tofit);
1056  // tofit2 = new TVector3(xcen1,ycen1,z_bis);
1057  tofit2 = new TVector3(x_,y_,z_bis);
1058  ZPointsArray->Add(tofit2);
1059 
1060  if( fDisplayLevel >= 4) {
1061  eventDetails->cd(3);
1062 
1063  // xz plane
1064  // found point 1
1065  TMarker *mrkt2 = new TMarker(x_, z_, 6);
1066  mrkt2->SetMarkerColor(wireOk);
1067  mrkt2->Draw("SAME");
1068  // found point2
1069  TMarker *mrkt2bis = new TMarker(x_, z_bis, 6);
1070  mrkt2bis->SetMarkerColor(wireOk);
1071  mrkt2bis->Draw("SAME");
1072 
1073  eventDetails->cd(4);
1074  // yz plane
1075  // found point 1
1076  TMarker *mrkt2b = new TMarker(y_, z_, 6);
1077  mrkt2b->SetMarkerColor(wireOk);
1078  mrkt2b->Draw("SAME");
1079  // found point 2
1080  TMarker *mrkt2bbis = new TMarker(y_, z_bis, 6);
1081  mrkt2bbis->SetMarkerColor(wireOk);
1082  mrkt2bbis->Draw("SAME");
1083  }
1084 
1085  if( fDisplayLevel >= 2) {
1086  eventCanvas->cd(2);
1087  // MC point in z - track length plane
1088  TMarker *MCmrk = new TMarker(h* R*TMath::ATan2((iPoint->GetY() - y0)*TMath::Cos(Phi0) - (iPoint->GetX() - x0)*TMath::Sin(Phi0) , R + (iPoint->GetX() - x0) * TMath::Cos(Phi0) + (iPoint->GetY()-y0) * TMath::Sin(Phi0)), iPoint->GetZ(), 6);
1089  MCmrk->SetMarkerColor(4);
1090  MCmrk->Draw("SAME");
1091 
1092  }
1093 
1094  // ------- HOUGH TRANSFORM ------------
1095  // FIRST CHOICE
1096  if(tofit) Hough(tofit, Phi0, x0, y0, R);
1097 
1098  // SECOND CHOICE
1099  if(tofit2) Hough(tofit2, Phi0, x0, y0, R);
1100  wireOk++;
1101  }
1102 
1103  }
1104 
1105  // cout << "skewed: " << wireOk << endl;
1106 
1107 
1108  TVector3 outz = GetHoughResponse();
1109 
1110  if( fDisplayLevel >= 3) {
1111  eventCanvas->cd(2);
1112  // found line in z track length plane
1113  TLine *line = new TLine(-40, -40*outz.X() + outz.Y(), 40, (40*outz.X() + outz.Y()));
1114  line->Draw("SAME");
1115  }
1116 
1117  // pTrack->GetParamFirst()->SetX(outz.Z()); // CHECK what is this?!!!!!
1118 
1119  // if votes < 6 consider the fit failed CHECK
1120  // if(outz.Z() < 6) return kFALSE;
1121 
1122  Int_t counter = 0;
1123 
1124  if(outz.X() == 0. && outz.Y() == 0.) return kFALSE;
1125 
1126  if( wireOk < 2){
1127  return kFALSE;
1128  }
1129 
1130  //first = kTRUE; //[R.K.01/2017]unused variable?
1131  wireOk = 0;
1132  Int_t okcounter = 0;
1133  for(Int_t i = 0; i < (2*hitcounter); i+=2) {
1134  // get index of hit
1135  PndTrackCandHit candhit = pTrackCand->GetSortedHit(counter);
1136  Int_t iHit = candhit.GetHitId();
1137  counter++;
1138 
1139  // get hit
1140  PndSttHit *pMhit = (PndSttHit*) fHitArray->At(iHit);
1141  if (!pMhit ) { cout << "PndSttHelixTrackFitter::ZFinder: no hit at " << iHit << endl; continue; }
1142 
1143  Int_t refindex = pMhit->GetRefIndex();
1144  // get point
1145  PndSttPoint *iPoint = (PndSttPoint*) fPointArray->At(refindex);
1146  if (!iPoint ) { cout << "PndSttHelixTrackFitter::ZFinder: no point at " << refindex << " associated to hit " << iHit << endl; continue; }
1147 
1148  // tubeID CHECK added
1149  Int_t tubeID = pMhit->GetTubeID();
1150  PndSttTube *tube = (PndSttTube*) fTubeArray->At(tubeID);
1151  TVector3 wiredirection = tube->GetWireDirection();
1152 
1153  if(wiredirection == TVector3(0.,0.,1.)) continue;
1154  wireOk++;
1155 
1156  sigz = 1.; // CHECK
1157 
1158  // FIRST CHOICE .................
1159  TVector3 *vi = (TVector3*) ZPointsArray->At(okcounter);
1160 
1161  if(vi == NULL) continue;
1162  Double_t scos = fTrack->CalculateScosl(vi->X(), vi->Y());
1163 
1164  if( fDisplayLevel >= 3) {
1165  eventCanvas->cd(2);
1166  // first point z track lenght plane
1167  TMarker *mrk = new TMarker(scos, vi->Z(), 6);
1168  mrk->Draw("SAME");
1169  }
1170 
1171  Bool_t fitdone = kFALSE;
1172 
1173  // distance between the found z and the predicted from the found line one
1174  Double_t distfirst = pow(((vi->Z() - (outz.Y() + outz.X() * scos))/sigz), 2);
1175  if(distfirst < 10) {
1176  pMhit->SetX(vi->X());
1177  pMhit->SetY(vi->Y());
1178  pMhit->SetZ(vi->Z());
1179  fitdone = kTRUE;
1180  }
1181 
1182 
1183  // SECOND CHOICE ...................
1184  vi = (TVector3*) ZPointsArray->At(okcounter+1);
1185  if(vi == NULL) continue;
1186  scos = fTrack->CalculateScosl(vi->X(), vi->Y());
1187 
1188  Double_t distsecond = pow(((vi->Z() - (outz.Y() + outz.X() * scos))/sigz), 2);
1189 
1190  if( fDisplayLevel >= 3) {
1191  eventCanvas->cd(2);
1192  TMarker *mrk2 = new TMarker(scos, vi->Z(), 6);
1193  mrk2->Draw("SAME");
1194  }
1195 
1196  // Xint Yin Zint set
1197  if(fitdone == kTRUE){
1198  if(distsecond < 10 && distsecond < distfirst) {
1199  pMhit->SetX(vi->X());
1200  pMhit->SetY(vi->Y());
1201  pMhit->SetZ(vi->Z());
1202  }
1203  }
1204  else {
1205  if (distsecond < 10) {
1206  pMhit->SetX(vi->X());
1207  pMhit->SetY(vi->Y());
1208  pMhit->SetZ(vi->Z());
1209  fitdone = kTRUE;
1210  }
1211  else {
1212  pMhit->SetX(-999);
1213  pMhit->SetY(-999);
1214  pMhit->SetZ(-999);
1215  }
1216  }
1217 
1218  if( fDisplayLevel >= 3) {
1219  if(pMhit->GetZ()!= -999) {
1220  eventCanvas->cd(2);
1221  TMarker *mrk3 = new TMarker(scos, pMhit->GetZ(), 6);
1222  mrk3->SetMarkerColor(3);
1223  mrk3->Draw("SAME");
1224  }
1225  }
1226 
1227  okcounter = okcounter + 2;
1228  }
1229 
1230  return kTRUE;
1231  // fiterrp = sqrt(fiterrp2);
1232  // fiterrm = sqrt(fiterrm2);
1233 }
1234 
1235 
1236 
1237 // ZFit was Zfitbb2
1238 // ----- Zfit ----------------------------------------
1239 Int_t PndSttHelixTrackFitter::ZFit(PndTrackCand* pTrackCand, Int_t ) { //whatToFit //[R.K.03/2017] unused variable(s)
1240 
1241  if(fVerbose == 2) cout << "ZFIT" << endl;
1242 
1243  if(!pTrackCand) return 0;
1244 
1245  Int_t hitcounter = pTrackCand->GetNHits();
1246 
1247  // cut on number of hits
1248  // if(hitcounter > 50) return 0;
1249  if(hitcounter < 5) return 0;
1250 
1251  // CHECK: ATTENTION
1252  // refit for error correction missing
1253 
1254 
1255  // SCOSL ======
1256  // get 1st hit
1257  //PndSttHit *fMhit = (PndSttHit*) fHitArray->At(pTrackCand->GetSortedHit(0).GetHitId()); //[R.K. 01/2017] unused variable?
1258 
1259 
1260  Double_t Sxx, Sx, Sz, Sxz, S1z;
1261  Double_t Detz = 0.;
1262  Double_t fitm, fitp;
1263  Double_t sigz = 1.; // CHECK
1264 
1265  Sx = 0.;
1266  Sz = 0.;
1267  Sxx = 0.;
1268  Sxz = 0.;
1269  S1z = 0.;
1270 
1271  // centre of curvature
1272  //Double_t x_0 = (fTrack->GetDist() + fTrack->GetRad()) * cos(fTrack->GetPhi()); //[R.K. 01/2017] unused variable?
1273  //Double_t y_0 = (fTrack->GetDist() + fTrack->GetRad()) * sin(fTrack->GetPhi()); //[R.K. 01/2017] unused variable?
1274  // radius of curvature
1275  //Double_t R = fTrack->GetRad(); //[R.K. 01/2017] unused variable?
1276  Int_t counter = 0;
1277  Int_t wireOk = 0;
1278  //Bool_t first = kTRUE; //[R.K. 01/2017] unused variable?
1279  for (Int_t i = 0; i < hitcounter; i++) {
1280 
1281  // get index of hit
1282  PndTrackCandHit candhit = fTrackCand->GetSortedHit(i);
1283  Int_t iHit = candhit.GetHitId();
1284 
1285  // get hit
1286  PndSttHit *pMhit = (PndSttHit*) fHitArray->At(iHit);
1287  if (!pMhit ) { cout << "PndSttHelixTrackFitter::ZFit: no hit at " << iHit << endl; continue; }
1288 
1289  if(pMhit->GetX() == -999 || pMhit->GetY() == -999 || pMhit->GetZ() == -999) continue; // CHECK
1290 
1291  Int_t refindex = pMhit->GetRefIndex();
1292  // get point
1293  PndSttPoint *iPoint = (PndSttPoint*) fPointArray->At(refindex);
1294  if (!iPoint ) { cout << "PndSttHelixTrackFitter::ZFit: no point at " << refindex << " associated to hit " << iHit << endl; continue; }
1295 
1296  // tubeID CHECK added
1297  Int_t tubeID = pMhit->GetTubeID();
1298  PndSttTube *tube = (PndSttTube*) fTubeArray->At(tubeID);
1299 
1300  TVector3 wiredirection = tube->GetWireDirection();
1301 
1302  if(wiredirection == TVector3(0.,0.,1.)) continue;
1303  // if(tube->GetPosition().Z() != 35) continue; // for the moment I throw away short tubes
1304 
1305  TVector3 *vi = new TVector3(pMhit->GetX(), pMhit->GetY(), pMhit->GetZ());
1306 
1307  // if the found z is > 75 cm or < -75 cm continue: this has to be fixed
1308  if(pMhit->GetZ() < (tube->GetPosition().Z() - tube->GetHalfLength()) || pMhit->GetZ() > (tube->GetPosition().Z() + tube->GetHalfLength())) continue; // CHECK
1309 
1310  Double_t scos = fTrack->CalculateScosl(pMhit->GetX(), pMhit->GetY());
1311 
1312  Sx = Sx + (scos /(sigz * sigz));
1313  Sz = Sz + (vi->Z()/(sigz * sigz));
1314  Sxz = Sxz + ((scos *vi->Z())/(sigz * sigz));
1315  Sxx = Sxx + ((scos * scos)/(sigz * sigz));
1316  S1z = S1z + 1/(sigz * sigz);
1317  counter++;
1318  wireOk++;
1319 
1320  }
1321 
1322  if(counter == 0) return 0;
1323 
1324  Detz = S1z*Sxx - Sx*Sx;
1325  if(Detz == 0) {
1326  return 0;
1327  }
1328 
1329  fitp = (1/Detz)*(Sxx*Sz - Sx*Sxz);
1330  fitm = (1/Detz)*(S1z*Sxz - Sx*Sz);
1331 
1332  // y = m*x + p
1333  TVector2 outz(fitm, fitp);
1334 
1335  // fill in parameters
1336  fTrack->SetTanL(fitm); // tan(lambda)
1337  fTrack->SetZ(fitp); // z0
1338 
1339  if( fDisplayLevel >= 3) {
1340  eventCanvas->cd(2);
1341 
1342  // fit line
1343  TLine *line2 = new TLine(-40, -40*fitm + fitp, 40, (40*fitm + fitp));
1344  line2->SetLineColor(kGreen - 9);
1345  line2->Draw("SAME");
1346  }
1347  return 1;
1348  // fiterrp = sqrt(fiterrp2);
1349  // fiterrm = sqrt(fiterrm2);
1350 }
1351 // -------------------------------------------------------------------------------
1352 
1353 void PndSttHelixTrackFitter::Hough(TVector3* choice, Double_t Phi0, Double_t x0, Double_t y0, Double_t R)
1354 {
1355 
1356  // since the points stay on a line in z track length plane
1357  // I can use the hough transform to transform the points in
1358  // z track length plane to lines in the parameters plane ab
1359  // (scos = a * z + b) and find their intersections.
1360  //
1361  // grid:
1362  // a [-pi/2 , pi/2]; 200 steps
1363  // b [-150, 150]; 1000 steps // CHECK deve essere -75, 75
1364 
1365  Double_t y = choice->Z();
1366  Double_t x = h* R*TMath::ATan2((choice->Y() - y0)*TMath::Cos(Phi0) - (choice->X() - x0)*TMath::Sin(Phi0) , R + (choice->X() - x0) * TMath::Cos(Phi0) + (choice->Y()-y0) * TMath::Sin(Phi0));
1367 
1368 
1369  // Hough -------------------------- angle [-90, 90]
1370  Double_t a = -TMath::Pi()/2.;
1371  Double_t tga;
1372  Double_t b = -150;
1373 
1374  // 200 step per -pi/2 < a < pi/2
1375  for(Int_t i = 0; i < 200; i++) {
1376  b = -150;
1377  tga = TMath::Tan(a);
1378  // 1000 step per -150 < b < 150
1379  for(Int_t j = 0; j <= 1000; j++) {
1380  Double_t ycalc = x * tga + b;
1381 
1382  // if(fabs(ycalc - y) < (0.3 + TMath::Tan(TMath::Pi()/200.) * x)) {
1383  if(fabs(ycalc - y) < 0.2) {
1384  vote[i][j] += 1;
1385  }
1386 
1387  b += 0.3;
1388  }
1389  a += (TMath::Pi()/200.);
1390  }
1391 }
1392 
1393 TVector3 PndSttHelixTrackFitter::GetHoughResponse()
1394 {
1395  Double_t a = -TMath::Pi()/2.;
1396  Double_t tga;
1397  Double_t b = -150;
1398  TMarker *mrk2;
1399 
1400  Double_t aref, bref, vref = 0;
1401  //Double_t aref2, bref2, vref2 = 0, cref2 = 0; //[R.K. 01/2017] unused variable?
1402  for(Int_t i=0; i <= 200; i++)
1403  {
1404  tga = TMath::Tan(a);
1405  b = -150;
1406  for(Int_t j = 0; j <= 1000; j++) {
1407  if(vote[i][j] != 0)
1408  {
1409  if( fDisplayLevel >= 4) {
1410  eventDetails->cd(2);
1411  mrk2 = new TMarker(tga, b, 6);
1412  mrk2->SetMarkerColor(3);
1413  mrk2->Draw("SAME");
1414  }
1415 
1416  if(vote[i][j] >= vref) {
1417  aref = tga;
1418  bref = b;
1419  vref = vote[i][j];
1420 
1421  // if(vote[i][j] > vref) {
1422  // aref2 = tga;
1423  // bref2 = b;
1424  // vref2 = vote[i][j];
1425  // cref2 = 1;
1426  // }
1427  // else if(vote[i][j] > vref) {
1428  // aref2 += tga;
1429  // bref2 += b;
1430  // cref2++;
1431  // }
1432  }
1433  }
1434  b += 0.3;
1435  }
1436  a += (TMath::Pi()/200.);
1437  }
1438  // if(vref < vref2)
1439  // {
1440  // aref = aref2/cref2;
1441  // bref = bref2/cref2;
1442  // vref = vref2;
1443  // }
1444 
1445  TVector3 hough(aref, bref, vref);
1446 
1447  // reset votes
1448  for(Int_t i=0; i <= 200; i++) for(Int_t j = 0; j <= 1000; j++) vote[i][j] = 0;
1449  return hough;
1450 
1451 }
1452 
1453 
1454 // ------ Extrapolate ------------------------------------------------------------
1455 void PndSttHelixTrackFitter::Extrapolate(PndSttTrack* , Double_t , FairTrackParam * ) //track, r, param //[R.K.03/2017] unused variable(s)
1456 {
1457  cout << "-W- PndSttMinuitTrackFitter::Extrapolate: Not yet implemented, sorry!"
1458  << endl;
1459 }
1460 
1461 
1462 void PndSttHelixTrackFitter::ResetMArray() {
1463  for(Int_t i = 0; i < 100; i++) marray.AddAt(-999, i);
1464 }
1465 
1466 
1467 Int_t PndSttHelixTrackFitter::SetUpFitVector(PndTrackCand* pTrackCand, TMatrixT<Double_t> &fitvect)
1468 {
1469  int nhits = pTrackCand->GetNHits();
1470  fitvect.ResizeTo(nhits, 4); // x y r errr
1471 
1472  PndSttHit *currenthit = NULL;
1473  int counter = 0;
1474  for(int i = 0; i < nhits; i++)
1475  {
1476  PndTrackCandHit candhit = fTrackCand->GetSortedHit(i);
1477  Int_t iHit = candhit.GetHitId();
1478 
1479  currenthit = (PndSttHit*) fHitArray->At(iHit);
1480  if(!currenthit) { cout << "PndSttHelixTrackFitter::SetUpFitVector: no hit at " << iHit << endl; continue; }
1481 
1482 
1483  // tubeID CHECK added
1484  Int_t tubeID = currenthit->GetTubeID();
1485  PndSttTube *tube = (PndSttTube*) fTubeArray->At(tubeID);
1486 
1487  if(tube->GetWireDirection() != TVector3(0.,0.,1.)) continue;
1488  if(currenthit->GetX() == -999 || currenthit->GetY() == -999) continue;
1489 
1490  fitvect[counter][0] = currenthit->GetX();
1491  fitvect[counter][1] = currenthit->GetY();
1492  fitvect[counter][2] = currenthit->GetIsochrone();
1493  fitvect[counter][3] = currenthit->GetIsochroneError();
1494 
1495  counter++;
1496  }
1497  if(nhits != counter) {
1498  fitvect.ResizeTo(counter, 4); // x y r errr
1499  }
1500 
1501  return counter;
1502 }
1503 
1504 Int_t PndSttHelixTrackFitter::MinuitFit(PndTrackCand* pTrackCand, Int_t whatToFit)
1505 {
1506  if(fVerbose == 2) cout << "MINUIT FIT " << pTrackCand->GetNHits() << endl;
1507 
1508  fEventCounter++;
1509  //Double_t hitcounter = pTrackCand->GetNHits(); //[R.K. 01/2017] unused variable?
1510 
1511  Double_t rstart = fTrack->GetRad();
1512  Double_t xcstart = (fTrack->GetDist() + fTrack->GetRad()) * cos(fTrack->GetPhi());
1513  Double_t ycstart = (fTrack->GetDist() + fTrack->GetRad()) * sin(fTrack->GetPhi());
1514 
1515 
1516  TMinuit minimizer(3);
1517 
1518  // set to quiet
1519  if(fVerbose < 2) minimizer.SetPrintLevel(-1);
1520 
1521  if(fVerbose == 2) {
1522  cout << "*******MINUIT********" << endl;
1523  cout << "D SEED: " << fTrack->GetDist() << endl;
1524  cout << "PHI SEED: " << fTrack->GetPhi() << endl;
1525  cout << "R SEED: " << fTrack->GetRad() << endl;
1526  cout << "********************" << endl;
1527  }
1528 
1529  // set the object to be fitted:
1530  // TMatrixT<Double_t> [x][y][r][err_r]
1531  TMatrixT<Double_t> fitvect;
1532  SetUpFitVector(pTrackCand, fitvect);
1533 
1534  if(whatToFit == 1) minimizer.SetFCN(fcnHelix);
1535  else minimizer.SetFCN(fcnHelix2);
1536  // minimizer.SetErrorDef(1); // ???
1537 
1538  minimizer.DefineParameter(0, "xc", xcstart, 0.1, -3000., 3000.); // ???
1539  minimizer.DefineParameter(1, "yc", ycstart, 0.1, -3000., 3000.); // ??? LIMITS ???
1540  minimizer.DefineParameter(2, "r", rstart, 0.1, 0., 3000.); // ???
1541  if(fVerbose == 2) cout << "xcstart: " << xcstart << " ycxtart: " << ycstart << " rstart: " << rstart << endl;
1542  minimizer.SetObjectFit(&fitvect);
1543 
1544  minimizer.SetPrintLevel(-1);
1545 
1546  minimizer.SetMaxIterations(500);
1547 
1548  minimizer.Migrad();
1549 
1550  //Double_t chisquare; //[R.K. 01/2017] unused variable?
1551  Double_t resultsRadial[3], errorsRadial[3];
1552 
1553  minimizer.GetParameter(0, resultsRadial[0], errorsRadial[0]);
1554  minimizer.GetParameter(1, resultsRadial[1], errorsRadial[1]);
1555  minimizer.GetParameter(2, resultsRadial[2], errorsRadial[2]);
1556 
1557  // minimizer.Eval(3, NULL, chisquare, resultsRadial, 0); // ???
1558 
1559  if(fVerbose == 2) {
1560  cout << "xc: " << resultsRadial[0] << endl;
1561  cout << "yc: " << resultsRadial[1] << endl;
1562  cout << "R: " << resultsRadial[2] << endl;
1563  }
1564 
1565  Double_t phi = TMath::ATan2(resultsRadial[1], resultsRadial[0]); // CHECK
1566  Double_t d;
1567  d = ((resultsRadial[0] + resultsRadial[1]) - resultsRadial[2] *(TMath::Cos(phi) + TMath::Sin(phi)))/(TMath::Cos(phi) + TMath::Sin(phi)); // CHECK
1568 
1569 // pTrack->SetChi2Rad(chisquare); // CHECK how to handle this
1570 // pTrack->SetNDF(fTrackCand->GetNHits()); // CHECK how to handle this
1571 
1572  fTrack->SetDist(d);
1573  fTrack->SetPhi(phi);
1574  // fTrack->SetZ(z0);
1575  fTrack->SetRad(resultsRadial[2]);
1576  // fTrack->SetTanL(alpha);
1577  h = -GetCharge(d, phi, resultsRadial[2]);
1578  fTrack->SetCharge(-h);
1579 
1580  if( fDisplayLevel >= 3) {
1581  eventCanvas->cd(1);
1582  TArc *fitarc = new TArc(((fTrack->GetDist() + fTrack->GetRad()) * cos(fTrack->GetPhi())), ((fTrack->GetDist() + fTrack->GetRad()) * sin(fTrack->GetPhi())), fTrack->GetRad());
1583  fitarc->SetLineColor(kGreen - 9);
1584  fitarc->SetFillStyle(0);
1585  fitarc->Draw("SAME");
1586  eventCanvas->Update();
1587  eventCanvas->Modified();
1588  }
1589 
1590  return 1;
1591 }
1592 
1593 void fcnHelix(Int_t &, Double_t *, Double_t &f, Double_t *par, Int_t ) // npar, gin, iflag //[R.K.03/2017] unused variable(s)
1594 {
1595 
1596  TMatrixT<Double_t> *mama = (TMatrixT<Double_t> *)gMinuit->GetObjectFit();
1597 
1598  Double_t chisq = 0;
1599  Double_t delta = 0;
1600  Int_t hitcounter = mama->GetNrows();
1601  for (Int_t i = 0; i < hitcounter; i++)
1602  {
1603  delta =sqrt((mama[0][i][0]-par[0])*(mama[0][i][0]-par[0])+(mama[0][i][1]-par[1])*(mama[0][i][1]-par[1])) -par[2] ;
1604  if(mama[0][i][2] == 0) chisq += (delta * delta * 12.);
1605  else chisq += (delta*delta)/(mama[0][i][2] * mama[0][i][2] / 12.);
1606  }
1607  f = chisq;
1608 }
1609 
1610 void fcnHelix2(Int_t &, Double_t *, Double_t &f, Double_t *par, Int_t ) //npar, gin, iflag //[R.K.03/2017] unused variable(s)
1611 {
1612  TMatrixT<Double_t> *mama = (TMatrixT<Double_t> *)gMinuit->GetObjectFit();
1613 
1614  Double_t chisq = 0;
1615  Double_t delta = 0;
1616  Int_t hitcounter = mama->GetNrows();
1617  for (Int_t i = 0; i < hitcounter; i++)
1618  {
1619  delta =sqrt((mama[0][i][0]-par[0])*(mama[0][i][0]-par[0])+(mama[0][i][1]-par[1])*(mama[0][i][1]-par[1])) -par[2] ;
1620  if(mama[0][i][2] == 0) chisq += (delta * delta * 12.);
1621  else chisq += (delta*delta)/(pow(mama[0][i][3],2));
1622 
1623  }
1624 
1625  f = chisq;
1626 }
1627 
1628 
1629 
1630 
1631 PndSttHit* PndSttHelixTrackFitter::GetHitFromCollections(Int_t hitCounter) const
1632 {
1633  PndSttHit
1634  *retval = NULL;
1635 
1636  Int_t
1637  relativeCounter = hitCounter;
1638 
1639  for (Int_t collectionCounter = 0; collectionCounter < fHitCollectionList.GetEntries(); collectionCounter++)
1640  {
1641  Int_t
1642  size = ((TClonesArray *)fHitCollectionList.At(collectionCounter))->GetEntriesFast();
1643 
1644  if (relativeCounter < size)
1645  {
1646  retval = (PndSttHit*) ((TClonesArray *)fHitCollectionList.At(collectionCounter))->At(relativeCounter);
1647  break;
1648  }
1649  else
1650  {
1651  relativeCounter -= size;
1652  }
1653  }
1654 
1655  return retval;
1656 }
1657 
1658 Double_t PndSttHelixTrackFitter::GetHitAngle(Int_t hitNo, Double_t dCenter,
1659  Double_t phiCenter, Double_t radius)
1660 {
1661  PndSttHit
1662  *pMhit = NULL;
1663 
1664  pMhit = GetHitFromCollections(hitNo);
1665 
1666  // CHECK added
1667  Int_t tubeID = pMhit->GetTubeID();
1668  PndSttTube *tube = (PndSttTube*) fTubeArray->At(tubeID);
1669 
1670  Double_t
1671  xCenter = (dCenter + radius) * cos(phiCenter),
1672  yCenter = (dCenter + radius) * sin(phiCenter);
1673 
1674  // get XY of wire
1675  Double_t
1676  deltaX = tube->GetPosition().X() - xCenter,
1677  deltaY = tube->GetPosition().Y() - yCenter;
1678 
1679  // CHECK: the old way is method 1, but I thinkl method 2 is better. Check what
1680  // kind of angle exactly is used and choose one of the two methods (anyway, up to now
1681  // both works, so it is not urgent).
1682 
1683  // method 1
1684  Double_t angle = 0;
1685  if(deltaX != 0) angle = atan(deltaY / deltaX);
1686  else angle = TMath::Pi()/2.;
1687  if(deltaY < 0.) angle += dPi;
1688 
1689  // bring angle into the usual frame of reference
1690  if (deltaX < 0.)
1691  {
1692  if (deltaY < 0.)
1693  angle -= dPi;
1694  else
1695  angle += dPi;
1696  }
1697 
1698  // method 2
1699  // Double_t
1700  // angle = TMath::ATan2(deltaY , deltaX); // CHECK
1701  // if (deltaY < 0.) angle += (2 * dPi); // CHECK
1702 
1703  return angle;
1704 }
1705 
1706 void PndSttHelixTrackFitter::OrderHitsByR(map<Double_t, Int_t> &hitMap)
1707 {
1708  // sort the hits by r position
1709  PndSttHit
1710  *pMhit = NULL;
1711 
1712  if(fVerbose == 2) cout << "n of hits: " << fTrackCand->GetNHits() << endl;
1713 
1714  for (size_t i = 0; i < fTrackCand->GetNHits(); i++)
1715  {
1716  // get index of hit
1717  PndTrackCandHit candhit = fTrackCand->GetSortedHit(i);
1718  Int_t iHit = candhit.GetHitId();
1719 
1720  // get hit
1721  pMhit = GetHitFromCollections(iHit);
1722  if ( ! pMhit )
1723  continue;
1724 
1725  // tubeID CHECK added
1726  Int_t tubeID = pMhit->GetTubeID();
1727  PndSttTube *tube = (PndSttTube*) fTubeArray->At(tubeID);
1728  if(tube->GetWireDirection() != TVector3(0., 0., 1.)) continue;
1729  // hitMap[sqrt(pMhit->GetX() * pMhit->GetX() + pMhit->GetY() * pMhit->GetY())] = iHit;
1730  hitMap[tube->GetPosition().Perp()] = iHit;
1731  }
1732 }
1733 
1734 Int_t PndSttHelixTrackFitter::GetCharge(Double_t dCenter, Double_t phiCenter, Double_t radius)
1735 {
1736  map<Double_t, Int_t> hitMap;
1737  Int_t charge = 0;
1738 
1739  OrderHitsByR(hitMap);
1740 
1741  map<Double_t, Int_t>::iterator
1742  hitMapStart = hitMap.begin(),
1743  hitMapEnd = hitMap.end();
1744 
1745  hitMapEnd--;
1746 
1747  PndSttHit
1748  *startHit = GetHitFromCollections(hitMapStart->second),
1749  *endHit = GetHitFromCollections(hitMapEnd->second);
1750 
1751  // tubeID CHECK added
1752  Int_t starttubeID = startHit->GetTubeID();
1753  Int_t endtubeID = endHit->GetTubeID();
1754 
1755  PndSttTube
1756  *startTube = (PndSttTube*) fTubeArray->At(starttubeID),
1757  *endTube = (PndSttTube*) fTubeArray->At(endtubeID);
1758 
1759  PndSttGeomPoint
1760  vertex(0., 0., 0.),
1761  firstPoint(startTube->GetPosition().X(),
1762  startTube->GetPosition().Y(),
1763  startTube->GetPosition().Z()),
1764  lastPoint(endTube->GetPosition().X(),
1765  endTube->GetPosition().Y(),
1766  endTube->GetPosition().Z());
1767 
1768  //Bool_t
1769  //retval = kTRUE; //[R.K. 01/2017] unused variable?
1770 
1771  Double_t
1772  angle,
1773  oldAngle,
1774  difAngle;
1775 
1776  Int_t
1777  votesForPositive = 0,
1778  votesForNegative = 0,
1779  hitNo = 0;
1780 
1781  map<Double_t, Int_t>::iterator
1782  hitMapIter = hitMap.begin();
1783 
1784  while (hitMapIter != hitMap.end())
1785  {
1786  // get hit's angle
1787  angle = GetHitAngle(hitMapIter->second, dCenter, phiCenter, radius);
1788 
1789  // store the first hit's angle as a reference (hit with smallest z coordinate)
1790  if (hitNo == 0)
1791  {
1792  refAngle = angle;
1793  }
1794 
1795  // calculate the angle difference between the reference hit and this hit
1796  difAngle = angle - refAngle;
1797 
1798  while (difAngle < 0)
1799  {
1800  difAngle += 2 * dPi;
1801  }
1802 
1803  // cast votes for criterium 1
1804  if (hitNo > 0)
1805  {
1806  if (difAngle > oldAngle)
1807  votesForPositive++;
1808  else
1809  votesForNegative++;
1810  }
1811 
1812  oldAngle = difAngle;
1813 
1814  hitNo++;
1815  hitMapIter++;
1816  }
1817 
1818  // majority is enough to decide on track
1819  if (votesForPositive > votesForNegative) //2*
1820  {
1821 
1822  // when swimming the track from -z to z we get a positive rotation
1823  if (firstPoint.DistanceTo(vertex) < lastPoint.DistanceTo(vertex))
1824  {
1825  charge = -1;
1826  }
1827  else
1828  {
1829  charge = 1;
1830  }
1831  }
1832  else if (votesForNegative > votesForPositive) //2*
1833  {
1834  // when swimming the track from -z to z we get a positive rotation
1835  if (firstPoint.DistanceTo(vertex) < lastPoint.DistanceTo(vertex))
1836  {
1837  charge = 1;
1838  }
1839  else
1840  {
1841  charge = -1;
1842  }
1843  }
1844  else
1845  {
1846  charge = -1; // most likely an electron which makes a spiral
1847  }
1848 
1849  return charge;
1850 }
1851 
1852 
1853 // ===================================================================================
1854 // CONSTRAINT 1
1855 
1856 
1857 // XYFit spposing the track passes through (0, 0)
1858 Int_t PndSttHelixTrackFitter::XYFitThroughOrigin(PndTrackCand* pTrackCand, Int_t whatToFit) {
1859 
1860  if(!pTrackCand) return 0;
1861 
1862  Int_t hitcounter = pTrackCand->GetNHits();
1863  Bool_t first = kFALSE;
1864  if(hitcounter == 0) return 0;
1865  if(hitcounter < 5) { // hitcounter > 50
1866  if(fVerbose == 2) cout << "Bad No of hits in STT " << hitcounter << endl;
1867  return 0;
1868  }
1869  if(fVerbose == 2) cout << "FIT xy ********************" << endl;
1870 
1871 
1872  // traslation and rotation -----------
1873  // traslation near the first point
1874  Double_t trasl[2];
1875  // rotation
1876  Double_t alpha;
1877  if(fVerbose == 2) cout << "hitcounter: " << hitcounter << endl;
1878  for(Int_t k = 0; k <hitcounter; k++) {
1879 
1880  PndTrackCandHit candhit = pTrackCand->GetSortedHit(k);
1881  Int_t iHit = candhit.GetHitId();
1882  PndSttHit *currenthit = (PndSttHit*) fHitArray->At(iHit);
1883  if(!currenthit) continue;
1884  if(currenthit->GetX() == -999 || currenthit->GetY() == -999) continue;
1885 
1886  // tubeID CHECK added
1887  Int_t tubeID = currenthit->GetTubeID();
1888  PndSttTube *tube = (PndSttTube*) fTubeArray->At(tubeID);
1889 
1890  PndSttHit *hitfirst, *hitlast;
1891 
1892  TVector3 wiredirection = tube->GetWireDirection();
1893  if(wiredirection != TVector3(0.,0.,1.)) continue;
1894  else if(first == kFALSE)
1895  {
1896  hitfirst = (PndSttHit*) fHitArray->At(iHit);
1897  first = kTRUE;
1898  trasl[0] = hitfirst->GetX();
1899  trasl[1] = hitfirst->GetY();
1900  }
1901  else{
1902  hitlast = (PndSttHit*) fHitArray->At(iHit);
1903  alpha = TMath::ATan2(hitlast->GetY() - hitfirst->GetY(), hitlast->GetX() - hitfirst->GetX());
1904  }
1905  }
1906 
1907  // error <--> resolution
1908  Double_t sigr, sigx, sigy; // = 0.14; sigxy, //[R.K.01/2017]unused variable?
1909 
1910  // FITTING IN X-Y PLANE:
1911  // v = a + bu
1912 
1913  Double_t Suu, Su, Sv, Suv, S1;
1914  Su = 0.;
1915  Sv = 0.;
1916  Suu = 0.;
1917  Suv = 0.;
1918  S1 = 0.;
1919 
1920  Int_t digicounter = 0;
1921  TArrayD uarray(hitcounter);
1922  TArrayD varray(hitcounter);
1923  TArrayD sigv2array(hitcounter);
1924  TArrayD sigu2array(hitcounter);
1925 
1926  for(Int_t i=0; i < hitcounter; i++){
1927  uarray.AddAt(-999, i);
1928  varray.AddAt(-999, i);
1929  sigv2array.AddAt(-999, i);
1930  sigu2array.AddAt(-999, i);
1931  }
1932  for(Int_t i=0; i < hitcounter; i++){
1933 
1934  PndTrackCandHit candhit = pTrackCand->GetSortedHit(i);
1935  Int_t iHit = candhit.GetHitId();
1936  PndSttHit *currenthit = (PndSttHit*) fHitArray->At(iHit);
1937  if(!currenthit) { cout << "PndSttHelixTrackFitter::XYFit: no hit at " << iHit << endl; continue; }
1938  if(currenthit->GetX() == -999 || currenthit->GetY() == -999) continue;
1939 
1940  // tubeID CHECK added
1941  Int_t tubeID = currenthit->GetTubeID();
1942  PndSttTube *tube = (PndSttTube*) fTubeArray->At(tubeID);
1943 
1944  Int_t refindex = currenthit->GetRefIndex();
1945  // get point
1946  PndSttPoint *iPoint = (PndSttPoint*) fPointArray->At(refindex);
1947  if(!iPoint) { cout << "PndSttHelixTrackFitter::XYFit: no point at " << refindex << " associated to hit " << iHit << endl; continue; }
1948  TVector3 wiredirection = tube->GetWireDirection();
1949  if(wiredirection != TVector3(0.,0.,1.)) continue;
1950 
1951 
1952  if(whatToFit == 2) {
1953  //Double_t resx = iPoint->GetX() - currenthit->GetX(); //[R.K. 01/2017] unused variable?
1954  //Double_t resy = iPoint->GetY() - currenthit->GetY(); //[R.K. 01/2017] unused variable?
1955  //Double_t resdist = TMath::Sqrt((iPoint->GetY() - currenthit->GetY())*(iPoint->GetY() - currenthit->GetY()) + (iPoint->GetX() - currenthit->GetX())*(iPoint->GetX() - currenthit->GetX())); //[R.K. 01/2017] unused variable?
1956 
1957  }
1958 
1959  Double_t xtrasl, ytrasl;
1960  //Double_t xi, yi; //[R.K. 01/2017] unused variable?
1961  // traslation
1962  xtrasl = currenthit->GetX() - trasl[0];
1963  ytrasl = currenthit->GetY() - trasl[1];
1964 
1965  if(xtrasl == 0 && ytrasl == 0) continue;
1966 
1967  Double_t xrot, yrot;
1968  // rotation
1969  xrot = TMath::Cos(alpha)*xtrasl + TMath::Sin(alpha)*ytrasl;
1970  yrot = -TMath::Sin(alpha)*xtrasl + TMath::Cos(alpha)*ytrasl;
1971 
1972  if( fDisplayLevel >= 3) {
1973  eventCanvas->cd(1);
1974  TMarker *pt = new TMarker(xrot, yrot, 6);
1975  pt->SetMarkerColor(3);
1976  // if(whatToFit == 2)
1977  // pt->Draw("SAME");
1978  }
1979 
1980  // change coordinate
1981  Double_t u, v, sigv2, sigu2;
1982  u = xrot / (xrot*xrot + yrot*yrot);
1983  v = yrot / (xrot*xrot + yrot*yrot);
1984 
1985  if( fDisplayLevel >= 4) {
1986  eventDetails->cd(1);
1987  TMarker *uv = new TMarker(u, v, 6);
1988  uv->SetMarkerColor(3);
1989  if(whatToFit == 2) uv->Draw("SAME");
1990  eventDetails->Update();
1991  eventDetails->Modified();
1992  }
1993 
1994  if(whatToFit == 1) {
1995  if(currenthit->GetIsochrone() == 0) sigr = sqrt(2.)/sqrt(12.);
1996  else sigr = sqrt(2.) * currenthit->GetIsochrone()/sqrt(12.);
1997  sigx = sigr;
1998  sigy = sigr;
1999  }
2000  else {
2001  if(currenthit->GetIsochrone() == 0) {
2002  sigr = sqrt(2.)/sqrt(12.);
2003  sigx = sigr;
2004  sigy = sigr;
2005  }
2006  else {
2007  sigr = 0.0150;
2008  sigx = fabs(sigr * TMath::Cos(TMath::ATan(marray.At(i))));
2009  sigy = fabs(sigr * TMath::Sin(TMath::ATan(marray.At(i))));
2010  //sigxy = sigr * TMath::Sqrt(fabs(TMath::Cos(TMath::ATan(marray.At(i))) * TMath::Sin(TMath::ATan(marray.At(i)))));; //[R.K.01/2017]unused variable?
2011  }
2012  }
2013 
2014  Double_t dvdx = (-2 * xrot * yrot)/pow((xrot*xrot + yrot*yrot),2);
2015  Double_t dvdy = (xrot*xrot - yrot*yrot) / pow((xrot*xrot + yrot*yrot),2);
2016  Double_t dudx = (yrot*yrot - xrot*xrot) / pow((xrot*xrot + yrot*yrot),2);
2017  Double_t dudy = (-2 * xrot * yrot)/pow((xrot*xrot + yrot*yrot),2);
2018 
2019  sigu2 = dudx * dudx * sigx * sigx + dudy * dudy * sigy * sigy + 2 * dudx * dudy * sigx * sigy;
2020  sigv2 = dvdx * dvdx * sigx * sigx + dvdy * dvdy * sigy * sigy + 2 * dvdx * dvdy * sigx * sigy;
2021 
2022  uarray.AddAt(u, digicounter);
2023  varray.AddAt(v, digicounter);
2024  sigu2array.AddAt(sigu2, digicounter);
2025  sigv2array.AddAt(sigv2, digicounter);
2026 
2027  Su = Su + (u/sigv2);
2028  Sv = Sv + (v/sigv2);
2029 
2030  Suv = Suv + ((u*v)/sigv2);
2031  Suu = Suu + ((u*u)/sigv2);
2032 
2033  S1 = S1 + 1/sigv2;
2034 
2035  digicounter++;
2036  }
2037 
2038 
2039  Double_t determ = S1*Suu - Su*Su;
2040  if(determ == 0) {
2041  return 0;
2042  }
2043 
2044  // v = a + bu
2045  double fita = (1/determ)*(Suu*Sv - Su*Suv);
2046  double fitb = (1/determ)*(S1*Suv - Su*Sv);
2047 
2048  if(fVerbose == 2) {
2049  std::cout << "1) linear parameters:\n";
2050  std::cout << "a = " << fita << "\n";
2051  std::cout << "b = " << fitb << "\n";
2052  }
2053 
2054  Double_t chi2;
2055  chi2 = 0.;
2056  for(Int_t i=0; i < digicounter; i++){
2057  if(uarray.At(i) == -999 || varray.At(i) == -999 || sigv2array.At(i) == -999) continue;
2058  chi2 = chi2 + pow(((varray.At(i) - (fita + fitb*uarray.At(i)))/sqrt(sigv2array.At(i))),2) ;
2059  }
2060 
2061  if(fVerbose == 2) cout << "digicounter: " << digicounter << endl;
2062 
2111  if( fDisplayLevel >= 4) {
2112  eventDetails->cd(1);
2113  TLine *uvline = new TLine(-0.1, fita - 0.1 * fitb, 0.1, fita + 0.1 * fitb);
2114  if(whatToFit == 2) uvline->Draw("SAME");
2115  eventDetails->Update();
2116  eventDetails->Modified();
2117  }
2118 
2119  // center and radius
2120  Double_t xc, yc, xcrot, ycrot, R;
2121  // if(fabs(a)<0.000001) return 0; // CHECK proteggi dalla divisione per 0
2122  ycrot = 1 / (2 * fita);
2123  xcrot = - fitb * ycrot;
2124  R = sqrt(xcrot * xcrot + ycrot * ycrot);
2125 
2126  // re-rotation and re-traslation of xc and yc
2127  // rotation
2128  xc = TMath::Cos(alpha)*xcrot - TMath::Sin(alpha)*ycrot;
2129  yc = TMath::Sin(alpha)*xcrot + TMath::Cos(alpha)*ycrot;
2130  // traslation
2131  xc = xc + trasl[0];
2132  yc = yc + trasl[1];
2133 
2134  // // errors on parameters -------------------
2135  // MISSING
2136  // =============================================
2137 
2138  Double_t phi = TMath::ATan2(yc, xc);
2139  Double_t d;
2140  d = ((xc + yc) - R*(TMath::Cos(phi) + TMath::Sin(phi)))/(TMath::Cos(phi) + TMath::Sin(phi));
2141 
2142  fTrack->SetDist(d);
2143  fTrack->SetPhi(phi);
2144  // Double_t newZ = -999.; // CHECK da cambiare
2145  // fTrack->SetZ(newZ);
2146  fTrack->SetRad(R);
2147  // Double_t newTheta = -999.; // CHECK da cambiare
2148  // fTrack->SetTanL(newTheta);
2149  h = -GetCharge(d, phi, R);
2150  fTrack->SetCharge(-h);
2151 
2152  // cout << "FIT: " << xc << " " << yc << endl;
2153  // cout << "RAGGIO: " << R << endl;
2154 
2155  if(fDisplayLevel >= 3) {
2156  eventCanvas->cd(1);
2157  TArc *fitarc = new TArc(((fTrack->GetDist() + fTrack->GetRad()) * cos(fTrack->GetPhi())), ((fTrack->GetDist() + fTrack->GetRad()) * sin(fTrack->GetPhi())), fTrack->GetRad());
2158  if(whatToFit == 2) fitarc->SetLineColor(kRed - 9);
2159  fitarc->SetFillStyle(0);
2160  fitarc->Draw("SAME");
2161  eventCanvas->Update();
2162  eventCanvas->Modified();
2163  }
2164  return 1;
2165 }
2166 
2167 void PndSttHelixTrackFitter::HoughThroughOrigin(TVector3* choice, Double_t Phi0, Double_t x0, Double_t y0, Double_t R)
2168 {
2169 
2170  // since the points stay on a line in z track length plane
2171  // I can use the hough transform to transform the points in
2172  // z track length plane to lines in the parameters plane ab
2173  // (scos = a * z) and find their intersections.
2174  //
2175  // grid:
2176  // a [-pi/2 , pi/2]; 200 steps
2177  // b = 0
2178 
2179  Double_t y = choice->Z();
2180  Double_t x = h* R*TMath::ATan2((choice->Y() - y0)*TMath::Cos(Phi0) - (choice->X() - x0)*TMath::Sin(Phi0) , R + (choice->X() - x0) * TMath::Cos(Phi0) + (choice->Y()-y0) * TMath::Sin(Phi0));
2181 
2182 
2183  // Hough -------------------------- angle [-90, 90]
2184  Double_t a = -TMath::Pi()/2.;
2185  Double_t tga;
2186 
2187  // 200 step per -pi/2 < a < pi/2
2188  for(Int_t i = 0; i < 200; i++) {
2189  tga = TMath::Tan(a);
2190  Double_t ycalc = x * tga;
2191 
2192  // if(fabs(ycalc - y) < (0.3 + TMath::Tan(TMath::Pi()/200.) * x)) {
2193  if(fabs(ycalc - y) < 0.2) {
2194  votecon[i] += 1;
2195  }
2196  a += (TMath::Pi()/200.);
2197  }
2198 }
2199 
2200 TVector3 PndSttHelixTrackFitter::GetHoughResponseThroughOrigin()
2201 {
2202  Double_t a = -TMath::Pi()/2.;
2203  Double_t tga;
2204 
2205  //TMarker *mrk2; //[R.K. 01/2017] unused variable?
2206 
2207  Double_t aref, /*bref,*/ vref = 0; //[R.K. 01/2017] unused variable?
2208  //Double_t aref2, bref2, vref2 = 0, cref2 = 0; //[R.K. 01/2017] unused variable?
2209  for(Int_t i=0; i <= 200; i++)
2210  {
2211  tga = TMath::Tan(a);
2212 
2213  if(votecon[i] != 0)
2214  {
2215  if(fDisplayLevel >= 3) {
2216  hougcon->Fill(tga);
2217  }
2218 
2219  if(votecon[i]>= vref) {
2220  aref = tga;
2221  vref = votecon[i];
2222  }
2223  }
2224  a += (TMath::Pi()/200.);
2225  }
2226 
2227  if(fDisplayLevel >= 4) {
2228  eventDetails->cd(2);
2229  hougcon->Draw();
2230  }
2231  TVector3 hough(aref, 0, vref);
2232 
2233  // reset votes
2234  for(Int_t i=0; i <= 200; i++) votecon[i] = 0;
2235  return hough;
2236 
2237 }
2238 
2239 // ZFinder was ZFinderbb3
2240 // -------- ZFinder --------------------------------------------
2241 Bool_t PndSttHelixTrackFitter::ZFinderThroughOrigin(PndTrackCand* pTrackCand, Int_t ) { //whatToFit //[R.K.03/2017] unused variable(s)
2242  // the z finding procedure uses the hough transform to find the line
2243  // in the plane z - track length on which the correct points lie.
2244 
2245  if(fVerbose == 2) cout << "ZFINDER" << endl;
2246 
2247  if(!pTrackCand) return 0;
2248 
2249  Int_t hitcounter = pTrackCand->GetNHits();
2250 
2251  // cut on number of hits
2252  // if(hitcounter > 50) {
2253  // cout << "more than 50 hits " << hitcounter << endl;
2254  // return 0;
2255  // }
2256  if(hitcounter < 5) {
2257  if(fVerbose == 2) cout << "less than 5 hits" << endl;
2258  return 0;
2259  }
2260 
2261  ZPointsArray = new TObjArray(2*hitcounter);
2262 
2263  // ZFIT =======
2264  if(hitcounter == 0) return kFALSE;
2265 
2266  //Double_t Sxx, Sxz; //[R.K. 01/2017] unused variable?
2267  //Double_t fitm, fitp; //[R.K. 01/2017] unused variable?
2268  Double_t sigz = 1.; // CHECK
2269 
2270  //Sxx = 0.; //[R.K. 01/2017] unused variable?
2271  //Sxz = 0.; //[R.K. 01/2017] unused variable?
2272  // ============
2273 
2274  TVector3 *tofit, *tofit2;
2275 
2276  // from xy plane fit
2277  Double_t phi0 = fTrack->GetPhi();
2278  Double_t d0 = fTrack->GetDist();
2279  Double_t x0 = d0*TMath::Cos(phi0);
2280  Double_t y0 = d0*TMath::Sin(phi0);
2281 
2282  // centre of curvature
2283  Double_t x_0 = (fTrack->GetDist() + fTrack->GetRad()) * cos(fTrack->GetPhi());
2284  Double_t y_0 = (fTrack->GetDist() + fTrack->GetRad()) * sin(fTrack->GetPhi());
2285 
2286  // radius of curvature
2287  Double_t R = fTrack->GetRad();
2288  Int_t wireOk = 0;
2289  //Bool_t first = kTRUE; //[R.K. 01/2017] unused variable?
2290 
2291  for (Int_t i = 0; i < hitcounter; i++) {
2292 
2293  // get index of hit
2294  PndTrackCandHit candhit = pTrackCand->GetSortedHit(i);
2295  Int_t iHit = candhit.GetHitId();
2296 
2297  // get hit
2298  PndSttHit *pMhit = (PndSttHit*) fHitArray->At(iHit);
2299  if (!pMhit ) { cout << "PndSttHelixTrackFitter::ZFinder: no hit at " << iHit << endl; continue; }
2300  // tubeID CHECK added
2301  Int_t tubeID = pMhit->GetTubeID();
2302  PndSttTube *tube = (PndSttTube*) fTubeArray->At(tubeID);
2303 
2304  Int_t refindex = pMhit->GetRefIndex();
2305  // get point
2306  PndSttPoint *iPoint = (PndSttPoint*) fPointArray->At(refindex);
2307  if (!iPoint ) { cout << "PndSttHelixTrackFitter::ZFinder: no point at " << refindex << " associated to hit " << iHit << endl; continue; }
2308 
2309  TVector3 wiredirection = tube->GetWireDirection();
2310  TVector3 wiredirection2;
2311 
2312  wiredirection2 = tube->GetHalfLength() * wiredirection;
2313  TVector3 cenposition = tube->GetPosition();
2314  // if(cenposition.Z() != 35) continue; // to throw away short tubes
2315 
2316  TVector3 min, max;
2317  min = cenposition - wiredirection2;
2318  max = cenposition + wiredirection2;
2319 
2320  // first extremity
2321  Double_t x_1= min.X();
2322  Double_t y_1= min.Y();
2323  Double_t z_1= min.Z();
2324 
2325  // second extremity
2326  Double_t x_2= max.X();
2327  Double_t y_2= max.Y();
2328  Double_t z_2= max.Z();
2329 
2330  Double_t rcur = pMhit->GetIsochrone();
2331 
2332  Double_t x1 = -9999.;
2333  Double_t y1 = -9999.;
2334  Double_t x2 = -9999.;
2335  Double_t y2 = -9999.;
2336 
2337  // in xy plane: angle of the PCA to the origin
2338  // with respect to the curvature center
2339  Double_t Phi0 = TMath::ATan2((y0 - y_0),(x0 - x_0));
2340 
2341  if(wiredirection != TVector3(0.,0.,1.))
2342  {
2343  Double_t a = -999;
2344  Double_t b = -999;
2345 
2346  // intersection point between the reconstructed
2347  // circumference and the line joining the centres
2348  // of the reconstructed circle and the i_th drift circle
2349  if(fabs(x_2-x_1)>0.0001) {
2350  a =(y_2-y_1)/(x_2-x_1);
2351  b =(y_1-a*x_1);
2352  Double_t A = a*a+1;
2353  Double_t B = x_0+a*y_0-a*b;
2354  Double_t C = x_0*x_0+y_0*y_0+b*b-R*R-2*y_0*b;
2355  if((B*B-A*C)>0) {
2356  x1= (B+TMath::Sqrt(B*B-A*C))/A;
2357  x2= (B-TMath::Sqrt(B*B-A*C))/A;
2358  y1=a*x1+b;
2359  y2=a*x2+b;
2360  }
2361  }
2362  else if(fabs(y_2-y_1)>0.0001) {
2363  Double_t A = 1;
2364  Double_t B = y_0;
2365  Double_t C = y_0*y_0 +(x_1-x_0)*(x_1-x_0) -R*R;
2366 
2367  if((B*B-A*C)>0) {
2368  y1= (B+TMath::Sqrt(B*B-A*C))/A;
2369  y2= (B-TMath::Sqrt(B*B-A*C))/A;
2370  x1=x2=x_1;
2371  }
2372  }
2373  else {
2374  if(fVerbose == 2) cout << "-E- intersection point not found" << endl;
2375  continue;
2376  }
2377 
2378 
2379 
2380  //x1 and x2 are the 2 intersection points
2381  Double_t d1=TMath::Sqrt((x1-cenposition.X())*(x1-cenposition.X())+
2382  (y1-cenposition.Y())*(y1-cenposition.Y()));
2383  Double_t d2=TMath::Sqrt((x2-cenposition.X())*(x2-cenposition.X())+
2384  (y2-cenposition.Y())*(y2-cenposition.Y()));
2385 
2386  Double_t x_ = x1;
2387  Double_t y_ = y1;
2388 
2389  // the intersection point nearest to the drift circle's centre is taken
2390  if(d2<d1) {x_=x2;y_=y2;}
2391 
2392  if(fDisplayLevel >= 3) {
2393  eventCanvas->cd(1);
2394  TMarker *pt = new TMarker(x_, y_, 6);
2395  pt->SetMarkerColor(6);
2396  pt->Draw("SAME");
2397 
2398  // intersection lines
2399  TLine *intline = new TLine(x_, y_, x_0, y_0);
2400  intline->SetLineColor(5);
2401  intline->Draw("SAME");
2402  }
2403 
2404  // now we need to find the actual centre of the drift circle,
2405  // by translating the drift circle until it becomes tangent
2406  // to the reconstructed circle.
2407  // Two solutions are possible (left rigth abiguity),
2408  // they are both kept, only the following zed fit will discard the wrong ones.
2409  // Using the parametric equation of the 3d-straigth line and taking the
2410  // x points just obtained, the zed coordinate of the skewed tube centre is calculated.
2411 
2412  if(x_1 == x_2)
2413  {
2414  // if the skewed layer lies in yz plane
2415  x1 = x_;
2416  y1 = y_ + rcur;
2417  x2 = x_;
2418  y2 = y_ - rcur;
2419  }
2420  else {
2421  //solving the equation to find out the centre of the tangent circle
2422  Double_t A = a*a+1;
2423  Double_t B = -(a*b-a*y_-x_);
2424  Double_t C = x_*x_+ y_*y_+b*b-2*b*y_-rcur*rcur;
2425  if((B*B-A*C)>0) {
2426  x1= (B+TMath::Sqrt(B*B-A*C))/A;
2427  x2= (B-TMath::Sqrt(B*B-A*C))/A;
2428  y1=a*x1+b;
2429  y2=a*x2+b;
2430  }
2431  else if((B*B-A*C)==0){
2432  x1= B/A;
2433  x2 = x1;
2434  y1=a*x1+b;
2435  y2=a*x2+b;
2436  }
2437  else {
2438  if(fVerbose == 2) cout << "NO WAY2" << endl;
2439  continue;
2440  }
2441  }
2442  if(fDisplayLevel >= 3) {
2443  eventCanvas->cd(1);
2444  TArc *drftarc = new TArc(x1, y1, rcur);
2445  drftarc->SetFillStyle(0);
2446  drftarc->SetLineColor(3);
2447  drftarc->Draw("SAME");
2448  TArc *drftarc2 = new TArc(x2, y2, rcur);
2449  drftarc2->SetFillStyle(0);
2450  drftarc2->SetLineColor(3);
2451  drftarc2->Draw("SAME");
2452  }
2453 
2454  d1=TMath::Sqrt((x1-cenposition.X())*(x1-cenposition.X())+(y1-cenposition.Y())*(y1-cenposition.Y()));
2455  d2=TMath::Sqrt((x2-cenposition.X())*(x2-cenposition.X())+(y2-cenposition.Y())*(y2-cenposition.Y()));
2456 
2457  Double_t xcen0=x1;
2458  Double_t xcen1=x2;
2459  Double_t ycen0=y1;
2460  Double_t ycen1=y2;
2461 
2462  if(d2<d1) { // z2 contains the points nearest (in x-y) to the initial centre of the skewed tube
2463  xcen0=x2;
2464  xcen1=x1;
2465  ycen0=y2;
2466  ycen1=y1;
2467 
2468  }
2469 
2470  // zed association
2471  Double_t t_, z_, t_bis, z_bis;
2472  if(fabs(x_2-x_1)<0.001) {
2473  t_ =(ycen0-y_1)/(y_2-y_1); // k= a_y*t + y_1 [t=1 y=y_2]
2474  z_ =(z_2-z_1)*t_ +z_1; // z= a_z*t + z_1 [t=1 z=z_2]
2475  t_bis =(ycen1-y_1)/(y_2-y_1); // from y_'s (the 2 solutions of the 2nd order equation)
2476  z_bis =(z_2-z_1)*t_bis +z_1; // and the 2 parametric equations the z coord. are obtained
2477  }
2478  else{
2479  t_ =(xcen0-x_1)/(x_2-x_1); // x= a_x*t + x_1 [t=1 x=x_2]
2480  z_ =(z_2-z_1)*t_ +z_1; // z= a_z*t + z_1 [t=1 z=z_2]
2481  t_bis =(xcen1-x_1)/(x_2-x_1); // from x_'s (the 2 solutions of the 2nd order equation)
2482  z_bis =(z_2-z_1)*t_bis +z_1; // and the 2 parametric equations the z coord. are obtained
2483  }
2484 
2485  // tofit = new TVector3(xcen0,ycen0,z_);
2486  tofit = new TVector3(x_,y_,z_);
2487  ZPointsArray->Add(tofit);
2488  // tofit2 = new TVector3(xcen1,ycen1,z_bis);
2489  tofit2 = new TVector3(x_,y_,z_bis);
2490  ZPointsArray->Add(tofit2);
2491 
2492  if(fDisplayLevel >= 4) {
2493  eventDetails->cd(3);
2494  // xz plane
2495  // found point 1
2496  TMarker *mrkt2 = new TMarker(x_, z_, 6);
2497  mrkt2->SetMarkerColor(wireOk);
2498  mrkt2->Draw("SAME");
2499  // found point2
2500  TMarker *mrkt2bis = new TMarker(x_, z_bis, 6);
2501  mrkt2bis->SetMarkerColor(wireOk);
2502  mrkt2bis->Draw("SAME");
2503 
2504 
2505  eventDetails->cd(4);
2506  // yz plane
2507  // found point 1
2508  TMarker *mrkt2b = new TMarker(y_, z_, 6);
2509  mrkt2b->SetMarkerColor(wireOk);
2510  mrkt2b->Draw("SAME");
2511  // found point 2
2512  TMarker *mrkt2bbis = new TMarker(y_, z_bis, 6);
2513  mrkt2bbis->SetMarkerColor(wireOk);
2514  mrkt2bbis->Draw("SAME");
2515  }
2516 
2517  if( fDisplayLevel >= 2) {
2518  eventCanvas->cd(2);
2519  // MC point in z - track length plane
2520  TMarker *MCmrk = new TMarker(h* R*TMath::ATan2((iPoint->GetY() - y0)*TMath::Cos(Phi0) - (iPoint->GetX() - x0)*TMath::Sin(Phi0) , R + (iPoint->GetX() - x0) * TMath::Cos(Phi0) + (iPoint->GetY()-y0) * TMath::Sin(Phi0)), iPoint->GetZ(), 6);
2521  MCmrk->SetMarkerColor(4);
2522  MCmrk->Draw("SAME");
2523  }
2524 
2525  // ------- HOUGH TRANSFORM ------------
2526 // // FIRST CHOICE
2527 // if(tofit) HoughThroughOrigin(tofit, Phi0, x0, y0, R); // CHECK
2528 
2529 // // SECOND CHOICE
2530 // if(tofit2) HoughThroughOrigin(tofit2, Phi0, x0, y0, R); // CHECK
2531 
2532  wireOk++;
2533  }
2534  }
2535 
2536  // cout << "skewed: " << wireOk << endl;
2537 
2538 
2539 // TVector3 outz = GetHoughResponseThroughOrigin(); // CHECK
2540  TVector3 outz = FindCorrectZ(ZPointsArray, x_0, y_0, x0, y0, R); // CHEC
2541  if(fDisplayLevel >= 3) {
2542  eventCanvas->cd(2);
2543 
2544  // found line in z track length plane
2545  TLine *line = new TLine(-40, -40*outz.X() + outz.Y(), 40, (40*outz.X() + outz.Y()));
2546  line->Draw("SAME");
2547  }
2548 
2549  // pTrack->GetParamFirst()->SetX(outz.Z()); // CHECK what is this?!!!!!
2550 
2551  // if votes < 6 consider the fit failed CHECK
2552  // if(outz.Z() < 6) return kFALSE;
2553 
2554  Int_t counter = 0;
2555 
2556  if(outz.X() == 0. && outz.Y() == 0.) return kFALSE;
2557 
2558  if( wireOk < 2){
2559  return kFALSE;
2560  }
2561 
2562  //first = kTRUE; //[R.K. 01/2017] unused variable?
2563  wireOk = 0;
2564  Int_t okcounter = 0;
2565  for(Int_t i = 0; i < (2*hitcounter); i+=2) {
2566  // get index of hit
2567  PndTrackCandHit candhit = pTrackCand->GetSortedHit(counter);
2568  Int_t iHit = candhit.GetHitId();
2569  counter++;
2570 
2571  // get hit
2572  PndSttHit *pMhit = (PndSttHit*) fHitArray->At(iHit);
2573  if (!pMhit ) { cout << "PndSttHelixTrackFitter::ZFinder: no hit at " << iHit << endl; continue; }
2574  // tubeID CHECK added
2575  Int_t tubeID = pMhit->GetTubeID();
2576  PndSttTube *tube = (PndSttTube*) fTubeArray->At(tubeID);
2577 
2578  Int_t refindex = pMhit->GetRefIndex();
2579  // get point
2580  PndSttPoint *iPoint = (PndSttPoint*) fPointArray->At(refindex);
2581  if (!iPoint ) { cout << "PndSttHelixTrackFitter::ZFinder: no point at " << refindex << " associated to hit " << iHit << endl; continue; }
2582 
2583 
2584  TVector3 wiredirection = tube->GetWireDirection();
2585 
2586 
2587  if(wiredirection == TVector3(0.,0.,1.)) continue;
2588  wireOk++;
2589 
2590  sigz = 1.; // CHECK
2591 
2592  // FIRST CHOICE .................
2593  TVector3 *vi = (TVector3*) ZPointsArray->At(okcounter);
2594 
2595  if(vi == NULL) continue;
2596  Double_t scos = fTrack->CalculateScosl(vi->X(), vi->Y());
2597 
2598  if(fDisplayLevel >= 3) {
2599  eventCanvas->cd(2);
2600  // first point z track lenght plane
2601  TMarker *mrk = new TMarker(scos, vi->Z(), 6);
2602  mrk->Draw("SAME");
2603  }
2604 
2605  Bool_t fitdone = kFALSE;
2606 
2607  // distance between the found z and the predicted from the found line one
2608  Double_t distfirst = pow(((vi->Z() - (outz.Y() + outz.X() * scos))/sigz), 2);
2609  if(distfirst < 10) {
2610  pMhit->SetX(vi->X());
2611  pMhit->SetY(vi->Y());
2612  pMhit->SetZ(vi->Z());
2613  fitdone = kTRUE;
2614  }
2615 
2616 
2617  // SECOND CHOICE ...................
2618  vi = (TVector3*) ZPointsArray->At(okcounter+1);
2619  if(vi == NULL) continue;
2620  scos = fTrack->CalculateScosl(vi->X(), vi->Y());
2621 
2622  Double_t distsecond = pow(((vi->Z() - (outz.Y() + outz.X() * scos))/sigz), 2);
2623 
2624  if(fDisplayLevel >= 3) {
2625  eventCanvas->cd(2);
2626  TMarker *mrk2 = new TMarker(scos, vi->Z(), 6);
2627  mrk2->Draw("SAME");
2628  }
2629 
2630  // Xint Yin Zint set
2631  if(fitdone == kTRUE){
2632  if(distsecond < 10 && distsecond < distfirst) {
2633  pMhit->SetX(vi->X());
2634  pMhit->SetY(vi->Y());
2635  pMhit->SetZ(vi->Z());
2636  }
2637  }
2638  else {
2639  if (distsecond < 10) {
2640  pMhit->SetX(vi->X());
2641  pMhit->SetY(vi->Y());
2642  pMhit->SetZ(vi->Z());
2643  fitdone = kTRUE;
2644  }
2645  else {
2646  pMhit->SetX(-999);
2647  pMhit->SetY(-999);
2648  pMhit->SetZ(-999);
2649  }
2650  }
2651 
2652  if(fDisplayLevel >= 3) {
2653  if(pMhit->GetZ()!= -999) {
2654  eventCanvas->cd(2);
2655  TMarker *mrk3 = new TMarker(scos, pMhit->GetZ(), 6);
2656  mrk3->SetMarkerColor(3);
2657  mrk3->Draw("SAME");
2658  }
2659  }
2660 
2661  okcounter = okcounter + 2;
2662  }
2663 
2664  return kTRUE;
2665  // fiterrp = sqrt(fiterrp2);
2666  // fiterrm = sqrt(fiterrm2);
2667 }
2668 
2669 
2670 
2671 // ZFit was Zfitbb2
2672 // ----- Zfit ----------------------------------------
2673 Int_t PndSttHelixTrackFitter::ZFitThroughOrigin(PndTrackCand* pTrackCand, Int_t ) {// whatToFit //[R.K.03/2017] unused variable(s)
2674 
2675  if(fVerbose == 2) cout << "ZFIT" << endl;
2676 
2677  if(!pTrackCand) return 0;
2678 
2679  Int_t hitcounter = pTrackCand->GetNHits();
2680 
2681  // cut on number of hits
2682  // if(hitcounter > 50) return 0;
2683  if(hitcounter < 5) return 0;
2684 
2685  // CHECK: ATTENTION
2686  // refit for error correction missing
2687 
2688 
2689  // SCOSL ======
2690  // get 1st hit
2691  //PndSttHit *fMhit = (PndSttHit*) fHitArray->At(pTrackCand->GetSortedHit(0).GetHitId()); //[R.K. 01/2017] unused variable?
2692 
2693 
2694  Double_t Sxx, Sxz;
2695  Double_t fitm, fitp;
2696  Double_t sigz = 1.; // CHECK
2697 
2698  Sxx = 0.;
2699  Sxz = 0.;
2700 
2701  // centre of curvature
2702  //Double_t x_0 = (fTrack->GetDist() + fTrack->GetRad()) * cos(fTrack->GetPhi()); //[R.K. 01/2017] unused variable?
2703  //Double_t y_0 = (fTrack->GetDist() + fTrack->GetRad()) * sin(fTrack->GetPhi()); //[R.K. 01/2017] unused variable?
2704  // radius of curvature
2705  //Double_t R = fTrack->GetRad(); //[R.K. 01/2017] unused variable?
2706  Int_t counter = 0;
2707  Int_t wireOk = 0;
2708  //Bool_t first = kTRUE; //[R.K. 01/2017] unused variable?
2709  for (Int_t i = 0; i < hitcounter; i++) {
2710 
2711  // get index of hit
2712  PndTrackCandHit candhit = fTrackCand->GetSortedHit(i);
2713  Int_t iHit = candhit.GetHitId();
2714 
2715  // get hit
2716  PndSttHit *pMhit = (PndSttHit*) fHitArray->At(iHit);
2717  if (!pMhit ) { cout << "PndSttHelixTrackFitter::ZFit: no hit at " << iHit << endl; continue; }
2718  // tubeID CHECK added
2719  Int_t tubeID = pMhit->GetTubeID();
2720  PndSttTube *tube = (PndSttTube*) fTubeArray->At(tubeID);
2721 
2722  if(pMhit->GetX() == -999 || pMhit->GetY() == -999 || pMhit->GetZ() == -999) continue; // CHECK
2723 
2724  Int_t refindex = pMhit->GetRefIndex();
2725  // get point
2726  PndSttPoint *iPoint = (PndSttPoint*) fPointArray->At(refindex);
2727  if (!iPoint ) { cout << "PndSttHelixTrackFitter::ZFit: no point at " << refindex << " associated to hit " << iHit << endl; continue; }
2728 
2729  TVector3 wiredirection = tube->GetWireDirection();
2730 
2731  if(wiredirection == TVector3(0.,0.,1.)) continue;
2732 
2733  // if(tube->GetPosition().Z() != 35) continue; // for the moment I throw away short tubes
2734 
2735  TVector3 *vi = new TVector3(pMhit->GetX(), pMhit->GetY(), pMhit->GetZ());
2736 
2737  // if the found z is > 75 cm or < -75 cm continue: this has to be fixed
2738  if(pMhit->GetZ() < (tube->GetPosition().Z() - tube->GetHalfLength()) || pMhit->GetZ() > (tube->GetPosition().Z() + tube->GetHalfLength())) continue; // CHECK
2739  Double_t scos = fTrack->CalculateScosl(pMhit->GetX(), pMhit->GetY());
2740 
2741  Sxz = Sxz + ((scos *vi->Z())/(sigz * sigz));
2742  Sxx = Sxx + ((scos * scos)/(sigz * sigz));
2743 
2744  counter++;
2745  wireOk++;
2746  }
2747 
2748  if(counter == 0) return 0;
2749 
2750  fitp = 0.;
2751  fitm = Sxz/Sxx;
2752 
2753  // y = m*x + p
2754  TVector2 outz(fitm, fitp);
2755 
2756  // fill in parameters
2757  fTrack->SetTanL(fitm); // tan(lambda)
2758  fTrack->SetZ(fitp); // z0
2759 
2760  if(fDisplayLevel >= 3) {
2761  eventCanvas->cd(2);
2762 
2763  // fit line
2764  TLine *line2 = new TLine(-40, -40*fitm + fitp, 40, (40*fitm + fitp));
2765  line2->SetLineColor(kGreen - 9);
2766  line2->Draw("SAME");
2767  }
2768  return 1;
2769  // fiterrp = sqrt(fiterrp2);
2770  // fiterrm = sqrt(fiterrm2);
2771 }
2772 // -------------------------------------------------------------------------------
2773 // ***************** CHECK *******************
2774 // TO USE THIS A BIG CHANGE IS NEEDED TO RETRIEVE THE START
2775 // POINT FROM PndTrack
2912 void PndSttHelixTrackFitter::InitEventDisplay()
2913 {
2914  // 0 no display
2915  // 1 only final result
2916  // 2 final + mc points
2917  // 3 also details
2918 
2919  h1 = new TH2F("h1","xy plane", 100, -50, 50, 100, -50, 50);
2920  h2 = new TH2F("h2","z cos#{lambda} plane", 100, -75, 75, 100, -45, 115);
2921  h3 = new TH2F("h3","conformal plane", 100, -1.5, 1.5, 100, -1.5, 1.5);
2922  h4 = new TH2F("h4","tubes", 100, -45, 45, 100, -45, 115);
2923 
2924  houg = new TH2F("houg","houg", 100,-1.001,1.001,100,-150.001,151.001);
2925  hougcon = new TH1F("hougcon","hougcon", 200,-1.001,1.001);
2926 
2927  switch(fDisplayLevel) {
2928  case 0:
2929  break;
2930  case 1:
2931  case 2:
2932  case 3:
2933  eventCanvas = new TCanvas("eventCanvas", "eventcanvas", 900, 500);
2934  eventCanvas->Divide(2, 1);
2935  break;
2936  case 4:
2937  eventCanvas = new TCanvas("eventCanvas", "eventcanvas", 0, 0, 400, 600);
2938  eventCanvas->Divide(1, 2);
2939  eventDetails = new TCanvas("eventDetails", "eventdetails", 400, 0, 600, 600);
2940  eventDetails->Divide(2, 2);
2941  break;
2942  }
2943 
2944 }
2945 
2946 Bool_t PndSttHelixTrackFitter::RunEventDisplay(PndTrackCand* pTrackCand)
2947 {
2948  // 0 no display
2949  // 1 only final result
2950  // 2 final + mc points
2951  // 3 also details
2952 
2953  eventCanvas->cd(1);
2954  h1->Draw();
2955  eventCanvas->cd(2);
2956  h2->Draw();
2957  if(fDisplayLevel >= 4) {
2958  eventDetails->cd(1);
2959  h3->Draw();
2960  eventDetails->cd(3);
2961  h4->Draw();
2962  eventDetails->cd(4);
2963  h4->Draw();
2964  eventDetails->cd(2);
2965  houg->Draw();
2966  }
2967  //.....................
2968  if(!pTrackCand) return kFALSE;
2969  Int_t hitcounter = pTrackCand->GetNHits();
2970 
2971  for(Int_t i = 0; i < hitcounter; i++){
2972  PndTrackCandHit candhit = pTrackCand->GetSortedHit(i);
2973  Int_t iHit = candhit.GetHitId();
2974  PndSttHit *currenthit = (PndSttHit*) fHitArray->At(iHit);
2975  if(!currenthit) continue;
2976  if(currenthit->GetX() == -999 || currenthit->GetY() == -999) continue;
2977  Int_t refindex = currenthit->GetRefIndex();
2978  // get point
2979  PndSttPoint *iPoint = (PndSttPoint*) fPointArray->At(refindex);
2980  if(!iPoint) continue;
2981  Int_t tubeID = currenthit->GetTubeID();
2982  PndSttTube *tube = (PndSttTube*) fTubeArray->At(tubeID);
2983  TVector3 wiredirection = tube->GetWireDirection();
2984 
2985  // draw MC points (BLUE)
2986  if(fDisplayLevel >= 2) {
2987  eventCanvas->cd(1);
2988  TMarker *cir0 = new TMarker(iPoint->GetX(), iPoint->GetY(), 6);
2989  cir0->SetMarkerColor(4);
2990  cir0->Draw("SAME");
2991 
2992  if(fDisplayLevel >= 4) {
2993  eventDetails->cd(3);
2994  TMarker *mcmrkt2bis = new TMarker(iPoint->GetX(), iPoint->GetZ(), 6);
2995  mcmrkt2bis->SetMarkerColor(4);
2996  mcmrkt2bis->Draw("SAME");
2997 
2998  eventDetails->cd(4);
2999  TMarker *mcmrkt2bbis = new TMarker(iPoint->GetY(), iPoint->GetZ(), 6);
3000  mcmrkt2bbis->SetMarkerColor(4);
3001  mcmrkt2bbis->Draw("SAME");
3002  }
3003  }
3004 
3005  // draw hit (GREEN = NON SKEWED / YELLOW = SKEWED)
3006  eventCanvas->cd(1);
3007  TArc *cir1 = new TArc(tube->GetPosition().X(), tube->GetPosition().Y(), currenthit->GetIsochrone());
3008  cir1->SetLineColor(3);
3009  if(wiredirection != TVector3(0.,0.,1.)) cir1->SetLineColor(5);
3010  cir1->SetFillStyle(0);
3011  cir1->Draw("SAME");
3012 
3013  if(fDisplayLevel >= 4) {
3014 
3015  // draw tube (BLACK LINE / PURPLE EXTREMITIES)
3016  TVector3 wiredirection2;
3017  wiredirection2 = tube->GetHalfLength() * wiredirection;
3018  TVector3 cenposition = tube->GetPosition();
3019 
3020  TVector3 min, max;
3021  min = cenposition - wiredirection2;
3022  max = cenposition + wiredirection2;
3023 
3024  // first extremity
3025  Double_t x_1= min.X();
3026  Double_t y_1= min.Y();
3027  Double_t z_1= min.Z();
3028 
3029  // second extremity
3030  Double_t x_2= max.X();
3031  Double_t y_2= max.Y();
3032  Double_t z_2= max.Z();
3033 
3034  if(wiredirection != TVector3(0.,0.,1.)) {
3035  eventDetails->cd(3);
3036 
3037  // first wire extremity
3038  TMarker *pt1z = new TMarker(x_1, z_1, 6);
3039  pt1z->SetMarkerColor(6);
3040  pt1z->Draw("SAME");
3041  // last wire extremity
3042  TMarker *pt2z = new TMarker(x_2, z_2, 6);
3043  pt2z->SetMarkerColor(6);
3044  pt2z->Draw("SAME");
3045  // tube line
3046  TLine *ztubeline = new TLine(x_1, z_1, x_2, z_2); // wire position
3047  ztubeline->Draw("SAME");
3048 
3049  eventDetails->cd(4);
3050 
3051  // first wire extremity
3052  TMarker *pt1z2 = new TMarker(y_1, z_1, 6);
3053  pt1z2->SetMarkerColor(6);
3054  pt1z2->Draw("SAME");
3055  // last wire extremity
3056  TMarker *pt2z2 = new TMarker(y_2, z_2, 6);
3057  pt2z2->SetMarkerColor(6);
3058  pt2z2->Draw("SAME");
3059  // tube line
3060  TLine *ztubeline2 = new TLine(y_1, z_1, y_2, z_2); // wire position
3061  ztubeline2->Draw("SAME");
3062 
3063  // xy
3064  eventCanvas->cd(1);
3065  TMarker *pt1 = new TMarker(x_1, y_1, 6);
3066  pt1->SetMarkerColor(6);
3067  pt1->Draw("SAME");
3068  TMarker *pt2 = new TMarker(x_2, y_2, 6);
3069  pt2->SetMarkerColor(6);
3070  pt2->Draw("SAME");
3071  TLine *tubeline = new TLine(x_1, y_1, x_2, y_2);
3072  tubeline->Draw("SAME");
3073  }
3074  }
3075  }
3076 
3077  return kTRUE;
3078 }
3079 
3080 void PndSttHelixTrackFitter::FinishEventDisplay(PndSttTrack *track)
3081 {
3082  if(!track) return;
3083 
3084 // cout << "param last x: " << track->GetDist() << endl;
3085 // cout << "param last y: " << track->GetPhi() << endl;
3086 // cout << "param last tx: " << track->GetRad() << endl;
3087 // cout << "param last ty: " << track->GetTanL() << endl;
3088 // cout << "param last qp: " << track->GetCharge() << endl;
3089 
3090  eventCanvas->cd(1);
3091  TArc *fitarc = new TArc(((track->GetDist() + track->GetRad()) * cos(track->GetPhi())),
3092  ((track->GetDist() + track->GetRad()) * sin(track->GetPhi())),
3093  track->GetRad());
3094  fitarc->SetLineColor(2);
3095  fitarc->SetFillStyle(0);
3096  fitarc->Draw("SAME");
3097 
3098  eventCanvas->cd(2);
3099 
3100  TLine *line = new TLine(-40, -40 * track->GetTanL() + track->GetZ(),
3101  40, 40 * track->GetTanL() + track->GetZ());
3102  line->SetLineColor(2);
3103  if(track->GetFlag() >= 3) line->Draw("SAME");
3104 
3105  eventCanvas->Update();
3106  eventCanvas->Modified();
3107  if( fDisplayLevel >= 4) {
3108  eventDetails->Update();
3109  eventDetails->Modified();
3110  }
3111 
3112 
3113 }
3114 
3115 
3116 TVector3 PndSttHelixTrackFitter::FindCorrectZ(TObjArray* choices, Double_t x_0, Double_t y_0, Double_t x0, Double_t y0, Double_t R)
3117 {
3118 
3119  int hitcounter = choices->GetEntriesFast();
3120  TH1F hlocal("hlocal", "", 200, -0.001, 6.281);
3121  Double_t Phi0 = TMath::ATan2((y0 - y_0),(x0 - x_0));
3122 
3123  TMatrixT<double> found_atan(hitcounter, 1);
3124  for(int i = 0; i < hitcounter; i++) {
3125 
3126  TVector3 *choice = (TVector3*) choices->At(i);
3127 
3128  Double_t y = choice->Z();
3129 
3130  Double_t x = h * R * TMath::ATan2((choice->Y() - y0)*TMath::Cos(Phi0) - (choice->X() - x0)*TMath::Sin(Phi0) , R + (choice->X() - x0) * TMath::Cos(Phi0) + (choice->Y()-y0) * TMath::Sin(Phi0));
3131 
3132 // TVector2 pos(choice->X() - x_0, choice->Y() - y_0);
3133 // TVector2 piv(x0 - x_0, y0 - y_0);
3134 // double alpha = TMath::ACos((pos.X() * piv.X() + pos.Y() * piv.Y()) / (pos.Mod() * piv.Mod()));
3135 // Double_t x = R * alpha;
3136 
3137  // the straingt line starting from (0, 0) is y = m * x;
3138  double atang = TMath::ATan2(y, x);
3139  if(y < 0) atang += (2 * TMath::Pi());
3140  found_atan[i][0] = atang;
3141 
3142  if(fDisplayLevel >= 4) hougcon->Fill(found_atan[i][0]);
3143  hlocal.Fill(found_atan[i][0]);
3144  // cout << "m " << i << " " << found_atan[i][0] * TMath::RadToDeg() << endl;
3145  }
3146 
3147  if(fDisplayLevel >= 4) {
3148  eventDetails->cd(2);
3149  hougcon->Draw();
3150  }
3151 
3152  double foundatan = hlocal.GetBinCenter(hlocal.GetMaximumBin());
3153  // cout << "found atan " << foundatan << " " << hlocal.GetMaximumBin() << endl;
3154 
3155  TVector3 foundz(TMath::Tan(foundatan), 0, hlocal.GetMaximumBin());
3156 
3157  return foundz;
3158 
3159 }
3160 
3161 ClassImp(PndSttTrackFitter)
PndSttHelixTrackFitter::XYFit
Int_t XYFit(PndTrackCand *pTrackCand, Int_t whatToFit)
Definition: PndSttHelixTrackFitter.cxx:210
x0
Double_t x0
Definition: checkhelixhit.C:70
PndSttHelixTrackFitter::GetHitAngle
Double_t GetHitAngle(Int_t hitNo, Double_t dCenter, Double_t phiCenter, Double_t radius)
Definition: PndSttHelixTrackFitter.cxx:1658
p
Double_t p
Definition: anasim.C:58
cos
friend F32vec4 cos(const F32vec4 &a)
Definition: P4_F32vec4.h:112
PndSttTrack::GetTanL
Double_t GetTanL()
Definition: PndSttTrack.h:60
PndSttHelixTrackFitter::Hough
void Hough(TVector3 *choice, Double_t Phi0, Double_t x0, Double_t y0, Double_t R)
Definition: PndSttHelixTrackFitter.cxx:1353
fcnHelix2
void fcnHelix2(Int_t &, Double_t *, Double_t &f, Double_t *par, Int_t)
Definition: PndSttHelixTrackFitter.cxx:1610
d
TObjArray * d
Definition: Pnd_Hc_etaee7G.C:78
choice
#define choice(c1, c2, c3)
Definition: PndCAMath.h:75
PndSttTrack::GetZ
Double_t GetZ()
Definition: PndSttTrack.h:61
PndSttHelixTrackFitter::Extrapolate
virtual void Extrapolate(PndSttTrack *track, Double_t r, FairTrackParam *param)
Definition: PndSttHelixTrackFitter.cxx:1455
i
Int_t i
Definition: run_full.C:25
m
__m128 m
Definition: P4_F32vec4.h:28
b
TTree * b
Definition: GammaSpectraAnalysis.C:24
PndSttTube::GetHalfLength
Double_t GetHalfLength()
Definition: PndSttTube.cxx:99
map
PndTransMap * map
Definition: sim_emc_apd.C:99
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
CAMath::Sqrt
static T Sqrt(const T &x)
Definition: PndCAMath.h:37
votecon
Double_t votecon[201]
Definition: PndSttHelixTrackFitter.cxx:45
verbose
#define verbose
Definition: PndRadMapBoxMesh.cxx:15
PndSttHelixTrackFitter::HoughThroughOrigin
void HoughThroughOrigin(TVector3 *choice, Double_t Phi0, Double_t x0, Double_t y0, Double_t R)
Definition: PndSttHelixTrackFitter.cxx:2167
nhits
Int_t nhits
Definition: analysis_point_hit_fwendcap.C:256
sin
friend F32vec4 sin(const F32vec4 &a)
Definition: P4_F32vec4.h:111
s
TLorentzVector s
Definition: Pnd2DStar.C:50
marray
TArrayD marray(200)
PndTrackCand::GetSortedHit
PndTrackCandHit GetSortedHit(UInt_t i)
Definition: PndTrackCand.h:54
CAMath::Sin
static T Sin(const T &x)
Definition: PndCAMath.h:42
PndSttHelixTrackFitter::ZFit
Int_t ZFit(PndTrackCand *pTrackCand, Int_t whatToFit)
Definition: PndSttHelixTrackFitter.cxx:1239
par
Double_t par[3]
Definition: anaLmdCluster.C:146
PndSttHelixTrackFitter::RunEventDisplay
Bool_t RunEventDisplay(PndTrackCand *trackCand)
Definition: PndSttHelixTrackFitter.cxx:2946
CAMath::Tan
float Tan(float x)
Definition: PndCAMath.h:165
GetEntriesFast
cout<< "POINTs for new FwEndCap == "<< tsim-> GetEntriesFast()
Definition: anal_point_fwendcap.C:90
PndSttHelixTrackFitter::ZFitThroughOrigin
Int_t ZFitThroughOrigin(PndTrackCand *pTrackCand, Int_t whatToFit)
Definition: PndSttHelixTrackFitter.cxx:2673
max
friend F32vec4 max(const F32vec4 &a, const F32vec4 &b)
Definition: P4_F32vec4.h:26
PndSttHelixTrackFitter::FindCorrectZ
TVector3 FindCorrectZ(TObjArray *choices, Double_t x_0, Double_t y_0, Double_t x0, Double_t y0, Double_t R)
Definition: PndSttHelixTrackFitter.cxx:3116
PndSttHelixTrackFitter::IntersectionFinder
Bool_t IntersectionFinder(PndTrackCand *pTrackCand)
Definition: PndSttHelixTrackFitter.cxx:635
PndSttHelixTrackFitter::FinishEventDisplay
void FinishEventDisplay(PndSttTrack *track)
Definition: PndSttHelixTrackFitter.cxx:3080
CAMath::Cos
static T Cos(const T &x)
Definition: PndCAMath.h:43
v
__m128 v
Definition: P4_F32vec4.h:4
PndSttTrack::GetDist
Double_t GetDist()
Definition: PndSttTrack.h:57
pt
TString pt(TString pts, TString exts="px py pz")
Definition: invexp.C:133
d0
Double_t d0
Definition: checkhelixhit.C:59
C
int Pic_FED Eff_lEE C()
Definition: Pic_FEP-Eff_lE2.C:1
PndSttTrack::SetCharge
void SetCharge(Int_t charge)
Definition: PndSttTrack.h:92
PndSttTrack::GetRad
Double_t GetRad()
Definition: PndSttTrack.h:59
PndSttTrack::SetZ
void SetZ(Double_t z)
Definition: PndSttTrack.h:91
a
Int_t a
Definition: anaLmdDigi.C:126
PndSttTrack::GetPhi
Double_t GetPhi()
Definition: PndSttTrack.h:58
counter
int counter
Definition: ZeeAnalysis.C:59
PndSttTrack::SetPhi
void SetPhi(Double_t phi)
Definition: PndSttTrack.h:88
PndSttTrack::SetDist
void SetDist(Double_t dist)
Definition: PndSttTrack.h:87
PndSttHit::GetIsochrone
Double_t GetIsochrone() const
Definition: PndSttHit.h:62
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
phi0
Double_t phi0
Definition: checkhelixhit.C:60
y0
Double_t y0
Definition: checkhelixhit.C:71
CAMath::ATan2
static T ATan2(const T &y, const T &x)
track
PndMCTrack * track
Definition: anaLmdCluster.C:89
PndSttTube::GetPosition
TVector3 GetPosition()
Definition: PndSttTube.cxx:87
f
TFile * f
Definition: bump_analys.C:12
PndSttHelixTrackFitter::DoFit
Int_t DoFit(PndTrackCand *pTrackCand, PndSttTrack *pTrack, Int_t pidHypo=211)
Definition: PndSttHelixTrackFitter.cxx:105
PndSttTrack::SetRad
void SetRad(Double_t r)
Definition: PndSttTrack.h:89
PndSttHit::GetTubeID
Int_t GetTubeID() const
Definition: PndSttHit.h:75
PndSttTrack::GetCharge
Int_t GetCharge()
Definition: PndSttTrack.h:63
PndSttHit::GetIsochroneError
Double_t GetIsochroneError() const
Definition: PndSttHit.h:63
min
friend F32vec4 min(const F32vec4 &a, const F32vec4 &b)
Definition: P4_F32vec4.h:25
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
PndSttHelixTrackFitter::OrderHitsByR
void OrderHitsByR(std::map< Double_t, Int_t > &hitMap)
Definition: PndSttHelixTrackFitter.cxx:1706
c
GeV c
Definition: ReadMesDecCorrRes_PiT.C:28
vote
Double_t vote[201][1001]
Definition: PndSttHelixTrackFitter.cxx:45
PndSttHelixTrackFitter::GetHitFromCollections
PndSttHit * GetHitFromCollections(Int_t hitCounter) const
Definition: PndSttHelixTrackFitter.cxx:1631
phi
double phi
Definition: macro/detectors/emc/dedicated/reco_analys.C:24
x
Double_t x
Definition: createRootFscGeometryFile.C:159
PndSttTrack::CalculateScosl
Double_t CalculateScosl(Double_t x, Double_t y)
Definition: PndSttTrack.cxx:129
PndTrackCand::GetNHits
UInt_t GetNHits() const
Definition: PndTrackCand.h:59
PndSttHelixTrackFitter::DoFitPlain
Int_t DoFitPlain(PndTrackCand *pTrackCand, PndSttTrack *pTrack, Int_t pidHypo=211)
Definition: PndSttHelixTrackFitter.cxx:125
dPi
#define dPi
Definition: PndSttHoughDefines.h:4
ClassImp
ClassImp(PndAnaContFact)
PndSttHelixTrackFitter::ZFinderThroughOrigin
Bool_t ZFinderThroughOrigin(PndTrackCand *pTrackCand, Int_t whatToFit)
Definition: PndSttHelixTrackFitter.cxx:2241
PndSttHelixTrackFitter::XYFitThroughOrigin
Int_t XYFitThroughOrigin(PndTrackCand *pTrackCand, Int_t whatToFit)
Definition: PndSttHelixTrackFitter.cxx:1858
y
Double_t y
Definition: createRootFscGeometryFile.C:159
alpha
double alpha
Definition: f_Init.h:9
angle
Double_t angle
Definition: createROOTEnDrcGeo.C:13
PndSttHelixTrackFitter::ZFinder
Bool_t ZFinder(PndTrackCand *pTrackCand, Int_t whatToFit)
Definition: PndSttHelixTrackFitter.cxx:807
PndSttHelixTrackFitter::MinuitFit
Int_t MinuitFit(PndTrackCand *pTrackCand, Int_t whatToFit)
Definition: PndSttHelixTrackFitter.cxx:1504
PndTrackCandHit::GetHitId
Int_t GetHitId() const
Definition: PndTrackCandHit.h:55
Pi
Double_t Pi
Definition: PndGeoHypGeCrystal.cxx:10
PndSttTrack::GetFlag
Int_t GetFlag() const
Definition: PndSttTrack.h:51
At
cout<<"the Event No is "<< i<< endl;{{if(hit_array->GetEntriesFast()!=mc_array->GetEntriesFast()) continue;PndSdsHit *hit=(PndSdsHit *) hit_array-> At(j)
Definition: anaLmdCluster.C:71
PndSttTrack::SetFlag
void SetFlag(Int_t flag)
Definition: PndSttTrack.h:95
R
Double_t R
Definition: checkhelixhit.C:61
PndSttHelixTrackFitter::GetCharge
Int_t GetCharge(Double_t dCenter, Double_t phiCenter, Double_t radius)
Definition: PndSttHelixTrackFitter.cxx:1734
PndSttHelixTrackFitter::SetUpFitVector
Int_t SetUpFitVector(PndTrackCand *pTrackCand, TMatrixT< Double32_t > &fitvect)
Definition: PndSttHelixTrackFitter.cxx:1467
vec
dble_vec_t vec[12]
Definition: ranlxd.cxx:380
PndSttTube::GetWireDirection
TVector3 GetWireDirection()
Definition: PndSttTube.cxx:107
TMatrixD
TMatrixT< double > TMatrixD
Definition: PndLmdDim.h:52
fcnHelix
void fcnHelix(Int_t &, Double_t *, Double_t &f, Double_t *par, Int_t)
Definition: PndSttHelixTrackFitter.cxx:1593
point
PndSdsMCPoint * point
Definition: anaLmdCluster.C:72
PndSttTrack::SetTanL
void SetTanL(Double_t tanl)
Definition: PndSttTrack.h:90
Generated on Wed Apr 3 2019 10:02:57 for FairRoot/PandaRoot by   doxygen 1.8.5