FairRoot/PandaRoot
PndRiemannTrack.cxx
Go to the documentation of this file.
1 //-----------------------------------------------------------
2 // File and Version Information:
3 // $Id$
4 //
5 // Description:
6 // Implementation of class PndRiemannTrack
7 // see PndRiemannTrack.hh for details
8 //
9 // Environment:
10 // Software developed for the PANDA Detector at FAIR.
11 //
12 // Author List:
13 // Sebastian Neubert TUM (original author)
14 //
15 //
16 //-----------------------------------------------------------
17 
18 // Panda Headers ----------------------
19 
20 // This Class' Header ------------------
21 #include "PndRiemannTrack.h"
22 
23 // C/C++ Headers ----------------------
24 #include <iostream>
25 #include <assert.h>
26 #include <cmath>
27 #include <ostream>
28 
29 // Collaborating Class Headers --------
30 #include "PndRiemannHit.h"
31 #include "FairTrackParP.h"
32 #include "FairRootManager.h"
33 #include "TGraph.h"
34 #include "TMath.h"
35 #include "TMatrixTSym.h"
36 #include "TF1.h"
37 // Class Member definitions -----------
38 
39 void MatrixOutput(TMatrixD mat)
40 {
41  for (int row = 0; row < mat.GetNrows(); row++)
42  {
43  for (int col = 0; col < mat.GetNcols(); col++)
44  {
45  std::cout << mat[row][col] << " ";
46  }
47  std::cout << std::endl;
48  }
49 }
50 
51 ClassImp(PndRiemannTrack);
52 
53 PndRiemannTrack::PndRiemannTrack() :
54  fn(3),fav(3), fc(0), fm(0), ft(0), fmError(0), ftError(0), fChi2(0), fcovPlane(4,4), fjacRXY(3,4), fcovRXY(3,3),
55  fVerbose(0), fFitDone(false), fSZFitDone(false), fErrorCalcDone(false), fweight(0),ftrefit(false),fVertexCut(0.5),
56  fStartAlpha(0), fStopAlpha(0)
57 {
58  fHits.clear();
59 }
60 
61 PndRiemannTrack::PndRiemannTrack(PndTrackCand* trackCand) :
62  fn(3),fav(3), fc(0), fm(0), ft(0), fmError(0), ftError(0), fChi2(0), fcovPlane(4,4), fjacRXY(3,4), fcovRXY(3,3),
63  fVerbose(0), fFitDone(false), fSZFitDone(false), fErrorCalcDone(false), fweight(0),ftrefit(false),fVertexCut(0.5),
64  fStartAlpha(0), fStopAlpha(0)
65 {
66  fHits.clear();
67  addPndTrackCand(trackCand);
68 }
69 
70 void PndRiemannTrack::addPndTrackCand(PndTrackCand* trackCand)
71 {
72  FairRootManager* ioman = FairRootManager::Instance();
73  std::set<FairLink> linksToHits = trackCand->GetLinks();
74  for (std::set<FairLink>::iterator iter = linksToHits.begin(); iter != linksToHits.end(); iter++)
75  {
76  PndRiemannHit riemannHit((FairHit*)ioman->GetCloneOfLinkData(*iter));
77  addHit(riemannHit);
78  }
79 }
80 
81 
82 PndRiemannTrack::~PndRiemannTrack()
83 {
84 }
85 
86 void
87 PndRiemannTrack::init(double x0_, double y0_, double R_,
88  double mydip, double){// z0 //[R.K.03/2017] unused variable(s)
89 
90  fn[1] = TMath::Sqrt(4*y0_*y0_/(4*y0_*y0_+4*x0_*x0_+1));
91  fn[0] = x0_/y0_*fn[1];
92  fn[2] = TMath::Sqrt(1 - fn[0]*fn[0] - fn[1]*fn[1]);
93  double fn2 = fn[2] * fn[2];
94  fc = fn2*4*R_*R_-1+fn2/(-4*fn[2]);
95 
96  fm = tan(acos(mydip));
97 
98 }
99 
100 
101 void
102 PndRiemannTrack::addHit(PndRiemannHit& hit){
103 // SetVerbose(3);
104  //int nbefore=fHits.size(); //[R.K. 01/2017] unused variable?
105  fHits.push_back(hit);
106 // std::cout << "fweight before addHit: " << fweight << std::endl;
107  //ADDED by ME//
108  if ( ftrefit ) {
109  fav *= fweight;
110  ftrefit = false;
111  }
112 
113  if (fVerbose > 1) std::cout << "-I- PndRiemannTrack::addHit " << fHits.size() -1 << ": " << hit.x().X() << " " << hit.x().Y() << " "
114  << hit.sigmaX() << " " << hit.sigmaY() << " " << hit.sigmaXY() << std::endl;
115  //fav*=(double)nbefore;
116  fav[0]+=hit.x().X() /(hit.sigmaXY()*hit.sigmaXY());
117  fav[1]+=hit.x().Y() /(hit.sigmaXY()*hit.sigmaXY());;
118  fav[2]+=hit.x().Z() /(hit.sigmaXY()*hit.sigmaXY());;
119  //fav*=1./(double)(nbefore+1);
120  fweight += 1/(hit.sigmaXY()*hit.sigmaXY());
121 
122 // std::cout << "fweight after addHit: " << fweight << std::endl;
123  fFitDone = false;
124  fSZFitDone = false;
125  fErrorCalcDone = false;
126 }
127 
128 double
129 PndRiemannTrack::dist(PndRiemannHit* hit){
130  double d2=fc;
131  d2+=hit->x().X()*fn[0];
132  d2+=hit->x().Y()*fn[1];
133  d2+=hit->x().Z()*fn[2];
134 
135  if (fVerbose > 1) {
136  std::cout << "PndRiemannTrack::dist: c " << fc << " n: " << fn[0] << "/" << fn[1] << "/" << fn[2] <<
137  " hit: " << hit->x().X() << "/" << hit->x().Y() << "/" << hit->x().Z() << " Dist: " << d2 << std::endl;
138  }
139  return d2;
140 }
141 
142 double
143 PndRiemannTrack::distError(PndRiemannHit* hit){
144  //for the time being the plane errors are ignored!
145  return TMath::Sqrt(TMath::Power(hit->sigmaX(),2) +TMath::Power(hit->sigmaY(),2) +TMath::Power(hit->sigmaW(),2));
146 }
147 
148 double
149 PndRiemannTrack::distCircle(PndRiemannHit* hit){
150  TVector2 pos;
151  TVectorD origin(2);
152  origin = orig();
153  TVector2 diff;
154 
155  pos.Set(hit->x().x(), hit->x().y());
156  diff.Set(pos.X() - origin[0], pos.Y() - origin[1]);
157 
158  return (diff.Mod() - r());
159 }
160 
161 
162 double PndRiemannTrack::ChiSquareDistCircle(){
163  double sum = 0;
164  for (size_t i = 0; i < fHits.size(); i++){
165  double distance = distCircle(&(fHits[i]));
166  sum += TMath::Power(distance / fHits[i].sigmaXY(), 2);
167  }
168  return sum;
169 }
170 
171 void
172 PndRiemannTrack::refit(bool withErrorCalc)
173 {
174 
175  TMatrixT<double> my_av(3,1);
176  if (fVerbose > 1) std::cout << "fweight: " << fweight << std::endl;
177  //fav *= TMath::Power(fweight,-1); //TS
179  if (!ftrefit){
180  fav *= TMath::Power(fweight,-1); //TS
181  }
183  my_av[0][0]=fav[0];// *fweight; //TS *fweight added
184  my_av[1][0]=fav[1];// *fweight;
185  my_av[2][0]=fav[2];// *fweight;
186 
187  if (fVerbose > 1) std::cout << "fav: " << fav[0] << " " << fav[1] << " " << fav[2] << std::endl;
188 
189  TMatrixD sampleCov(3,3);
190 
191  int nh=fHits.size();
192  for(int i=0;i<nh;++i){
193  TMatrixD h(3,1);
194  h[0][0]=fHits[i].x().X();
195  h[1][0]=fHits[i].x().Y();
196  h[2][0]=fHits[i].x().Z();
197  TMatrixD d(3,1);
198  d=h-my_av;
199  TMatrixD dt(TMatrixD::kTransposed,d);
200  TMatrixD ddt(d,TMatrixD::kMult,dt);
201  if (fHits[i].sigmaXY() != 0)
202  ddt *= 1/(fHits[i].sigmaXY()*fHits[i].sigmaXY()); //TS
203  else
204  std::cout << "-E- PndRiemannTrack::refit sigmaXY() == 0" << std::endl;
205  sampleCov+=ddt;
206  }
207 
208  if (nh != 0)
209  sampleCov*=1./(double)nh;
210  else {
211  std::cout << "-E- PndRiemannTrack::refit nh == 0" << std::endl;
212  return;
213  }
214 
215  if (fVerbose > 1) std::cout << "sampleCov: " << std::endl;
216  if (fVerbose > 1) MatrixOutput(sampleCov);
217 
218  TVectorD eigenValues(3);
219  TMatrixD eigenVec=sampleCov.EigenVectors(eigenValues);
220 
221  // find smallest eigenvalue
222 // double val=eigenValues[0]; int g=0;
223 // for(int k=1;k<3;++k){
224 // if(eigenValues[k]<val){
225 // val=eigenValues[k];
226 // g=k;
227 // }
228 // }
229 
230  int g = 0;
231  double smallestDistance = -1;
232  int smallestVec = 0;
233  TVectorD correctN(3);
234  double correctC;
235  double val;
236  for (int i = 0; i < 3; i++){
237  g = i;
238  fn=TMatrixDColumn(eigenVec,g);
239 
240  if (fVerbose > 1) std::cout << "eigenVec: " << std::endl;
241  if (fVerbose > 1) MatrixOutput(eigenVec);
242 
243  if (fVerbose > 1) std::cout << "eigenValues: ";
244  for (int j = 0; j < 3; j++){
245  if (fVerbose > 1) std::cout << eigenValues[j] << " ";
246  }
247  if (fVerbose > 1) std::cout << std::endl;
248 
249  double norm = 1;
250  if (fn.Norm2Sqr() != 0)
251  norm=1./TMath::Sqrt(fn.Norm2Sqr());
252 
253  fn*=norm;
254  fc=-1.*fn*fav;
255  if (smallestDistance < 0){
256  smallestDistance = ChiSquareDistCircle();
257  smallestVec = i;
258  correctN = fn;
259  correctC = fc;
260  val = eigenValues[i];
261  } else if (smallestDistance > ChiSquareDistCircle()) {
262  smallestDistance = ChiSquareDistCircle();
263  smallestVec = i;
264  correctN = fn;
265  correctC = fc;
266  val = eigenValues[i];
267  }
268 
269  if (fVerbose > 1) std::cout << "fn: " << fn[0] << " " << fn[1] << " " << fn[2] << std::endl;
270  if (fVerbose > 1) std::cout << "fc: " << fc << std::endl;
271  if (fVerbose > 1) {
272  TVectorD origin = orig();
273  std::cout << "Origin: " << origin[0] << " " << origin[1] << std::endl;
274  }
275  if (fVerbose > 1) std::cout << "r: " << r() << std::endl;
276  if (fVerbose > 1) std::cout << "ChiSquareCircleDistance: " << ChiSquareDistCircle() << std::endl;
277  }
278 
279  g = smallestVec;
280  fn = correctN;
281  fc = correctC;
282  //------------ Calculation of the full covariance matrix -----------
283 
284  // 1. Calculation of the covarianz matrix of the normal vector
285  if (withErrorCalc){
286  TMatrixD CovNorm(3,3);
287 
288  for (int i = 0; i < 3; i++){
289  if (i != g){
290  TMatrixD res = eigenVec.GetSub(0,2,i,i)*eigenVec.GetSub(0,2,i,i).T();
291  CovNorm += res * (val * eigenValues[i]/TMath::Power(val-eigenValues[i],2));
292  }
293  }
294  // if (fHits.size() > 5)
295  // CovNorm *= TMath::Power(fHits.size()-5,-1); //-5 is very strange. According to the paper this value has to be fixed with simulation???
296  // else
297  CovNorm *= TMath::Power(fHits.size(),-1);
298  if (fVerbose > 1) std::cout << "CovNorm: " << std::endl;
299  if (fVerbose > 1) MatrixOutput(CovNorm);
300 
301 
302  // 2. Calculation of the covarianz matrix of fav
303 
304  TMatrixD covAv(3,3);
305  for (size_t i = 0; i < fHits.size();i++){
306  covAv[0][0] += fHits[i].covX(0,0)/(TMath::Power(fHits[i].sigmaXY(),4));
307  covAv[1][1] += fHits[i].covX(1,1)/(TMath::Power(fHits[i].sigmaXY(),4));
308  covAv[1][0] += fHits[i].covX(1,0)/(TMath::Power(fHits[i].sigmaXY(),4));
309  covAv[2][0] += 2 * (fHits[i].covX(0,0) * fHits[i].x().X() + fHits[i].covX(1,0) * fHits[i].x().Y())
310  /(TMath::Power(fHits[i].sigmaXY(),4));
311  covAv[2][1] += 2 * (fHits[i].covX(1,1) * fHits[i].x().Y() + fHits[i].covX(1,0) * fHits[i].x().X())
312  /(TMath::Power(fHits[i].sigmaXY(),4));
313  covAv[2][2] += 4 * (fHits[i].covX(0,0) * fHits[i].x().X() * fHits[i].x().X() +
314  2 * fHits[i].covX(1,0) * fHits[i].x().X() * fHits[i].x().Y() +
315  fHits[i].covX(1,1) * fHits[i].x().Y() * fHits[i].x().Y())
316  /(TMath::Power(fHits[i].sigmaXY(),4));
317  }
318  covAv[0][1] = covAv[1][0];
319  covAv[0][2] = covAv[2][0];
320  covAv[1][2] = covAv[2][1];
321 
322  covAv *= TMath::Power(fweight,-2);
323 
324 
325  if (fVerbose > 1) std::cout << "covAv: " << std::endl;
326  if (fVerbose > 1) MatrixOutput(covAv);
327 
328 
329  // 3. Calculation of var(fc)
330 
331  double nCrn = (fn * (covAv * fn));
332  double rCnr = (fav * (CovNorm * fav));
333 
334  TMatrixD CnCr(CovNorm,TMatrixD::kMult,covAv);
335  double trCnCr = 0;
336  for (int i = 0; i < CnCr.GetNcols(); i++){
337  trCnCr += CnCr[i][i];
338  }
339  if (fVerbose > 1) std::cout << "nCrn: " << nCrn << " rCnr: " << rCnr << " trCnCr: " << trCnCr << std::endl;
340  double varc = nCrn + rCnr + trCnCr;
341 
342  TVectorD corr_cn = CovNorm * fav;
343  corr_cn *= -1;
344 
345  // 4. Calculation of the covarianz matrix of c,n
346 
347  fcovPlane[0][0] = varc;
348  fcovPlane[1][1] = CovNorm[0][0];
349  fcovPlane[2][1] = CovNorm[1][0];
350  fcovPlane[2][2] = CovNorm[1][1];
351  fcovPlane[3][1] = CovNorm[2][0];
352  fcovPlane[3][2] = CovNorm[2][1];
353  fcovPlane[3][3] = CovNorm[2][2];
354  fcovPlane[1][0] = corr_cn[0];
355  fcovPlane[2][0] = corr_cn[1];
356  fcovPlane[3][0] = corr_cn[2];
358 
359  fcovPlane[1][2] = CovNorm[1][0];
360 
361  fcovPlane[1][3] = CovNorm[2][0];
362  fcovPlane[2][3] = CovNorm[2][1];
363 
364  fcovPlane[0][1] = corr_cn[0];
365  fcovPlane[0][2] = corr_cn[1];
366  fcovPlane[0][3] = corr_cn[2];
368 
369  //5. Convert plane covariance in start parameter(r,x0,y0) covariances
370  calcJacRXY();
371 
372  if (fVerbose > 1) std::cout << "jacRXY: " << std::endl;
373  if (fVerbose > 1) MatrixOutput(fjacRXY);
374 
375  //fcovRXY = fjacRXY * fcovPlane * fjacRXY.T();
376 
377  TMatrixD temp(fjacRXY,TMatrixD::kMult, fcovPlane);
378  if (fVerbose > 1) std::cout << "temp: " << temp.GetNrows() << "x" << temp.GetNcols() << std::endl;
379 
380  if (fVerbose > 1) std::cout << "temp: " << std::endl;
381  if (fVerbose > 1) MatrixOutput(temp);
382 
383  fcovRXY = temp * fjacRXY.T(); // J * fcovPlane * J.T()
384 
385  if (fVerbose > 1) std::cout << "covPlane: " << std::endl;
386  if (fVerbose > 1) MatrixOutput(fcovPlane);
387 
388  if (fVerbose > 1) std::cout << "covRXY: " << std::endl;
389  if (fVerbose > 1) MatrixOutput(fcovRXY);
390  fErrorCalcDone = true;
391  }
392 
393  fFitDone = true;
394  fSZFitDone = false;
395  ftrefit = true;
396 
397  calcSForHits();
398 
399  sortHits();
400  calcStartStopAlpha();
401 }
402 
403 
404 TVectorD
405 PndRiemannTrack::orig() const {
406  TVectorD o(2);
407  if(fn[2]==0) return o; // RK is that good to return here?
408  double den=TMath::Power((2 * fn[2]),-1); // modified for paraboloid
409  o[0]=-fn[0]*den;
410  o[1]=-fn[1]*den;
411  return o;
412 }
413 
414 
415 double
416 PndRiemannTrack::r() const {
417 // if(fc==0)return 0;
418  //if(fabs(fc)>100)return 0;???????????????????????????
419  double a=2.*fn[2]; // modified for paraboloid
420  //if(a==0){
421  // if (fVerbose > 0) std::cout<<"PndRiemannTrack:: a==0 cannot calc r! set r=1E4"<<std::endl;
422  // return 1.E-4;
423  //}
424  //if (fVerbose > 0) std::cout<<fn[0]<<" "<<fn[1]<<" "<<fc<<" "<<a<<std::endl;
425  double nom = 1-fn[2]*fn[2]-4*fc*fn[2]; // modified for paraboloid
426  if(nom<0.0){
427  nom=1.E-4;
428  if (fVerbose > 1) std::cout<<"PndRiemannTrack::nom<0! set r=1E-4!"<<std::endl;
429  }
430  if (TMath::Abs(a) > 0)
431  return sqrt(nom)/TMath::Abs(a);
432  else return 0;
433 }
434 
435 double PndRiemannTrack::dR()
436 {
437  return TMath::Sqrt(fcovRXY[0][0]);
438 }
439 
440 void
441 PndRiemannTrack::szFit(bool withErrorCalc){
442  if (fFitDone == false)
443  refit(withErrorCalc);
444  unsigned int num=getNumHits();
445  if (fVerbose > 1) std::cout << "szFit() for " << num << " Points!" << std::endl;
446  if (r() > 0) {
447 
448  TGraph g(num);
449  // get s'es and zs
450  int j=0;
451  for(unsigned int i=0;i<num;++i){
452  if (fVerbose > 1)
453  if (fHits[i].hit() != NULL)
454  std::cout << "Point: " << i << ": " << fHits[i].hit()->GetEntryNr() << " ";
455  if (fHits[i].hit() != 0 && fHits[i].hit()->GetEntryNr().GetType() == GetBranchId("STTHit")){
456  continue;
457  }
458  fHits[i].calcPosOnTrk(this);
459  if (fVerbose > 1) std::cout << fHits[i].s() << " " << fHits[i].z() << std::endl;
460  g.SetPoint(j,fHits[i].s(),fHits[i].z());
461  j++;
462  }
463  if (j > 1){
464  g.Set(j);
465  g.Fit("pol1","Q0"); // << std::endl;
466  TF1* f = g.GetFunction("pol1");
467  //std::cout << "f: " << f << std::endl;
468  ft = f->GetParameter(0);
469  fm = f->GetParameter(1);
470  ftError = f->GetParError(0);
471  fmError = f->GetParError(1);
472  fChi2 = f->GetChisquare()/f->GetNDF();
473  fSZFitDone = true;
474  } else {
475  ft = 0;
476  fm = 0;
477  ftError = 0;
478  fmError = 0;
479  fChi2 = -1;
480  fSZFitDone = true;
481 
482  if (fVerbose > 0)
483  std::cout << "-W- PndRiemannTrack::calcSZ less than 2 points with z-Info: " << *this << std::endl;
484  }
485  } else {
486  ft = 0;
487  fm = 0;
488  ftError = 0;
489  fmError = 0;
490  fChi2 = -1;
491  fSZFitDone = true;
492 
493  if (fVerbose > 0)
494  std::cout << "-W- PndRiemannTrack::calcSZ r == 0: " << *this << std::endl;
495 
496  }
497  if (fVerbose > 1) std::cout << "t, m: " << ft << " +/- " << ftError << " / " << fm << " +/- " << fmError << " Chi2: " << fChi2 << std::endl;
498  return;
499 }
500 
501 double PndRiemannTrack::calcChi2Plane()
502 {
503  double chiSquare = 0;
504  if (fFitDone == false)
505  refit ();
506  if (r() > 0) {
507  for (size_t i = 0; i < getNumHits(); i++){
508  PndRiemannHit* actualHit = getHit(i);
509 // chiSquare += TMath::Power((dist(actualHit)/actualHit->sigmaXY()),2);
510 // std::cout << "Dist: " << " "<< dist(actualHit) << std::endl;
511  chiSquare += TMath::Power((dist(actualHit)/distError(actualHit)),2);
512  }
513  if (getNumHits() > 3)
514  chiSquare /= (getNumHits() - 3);
515  else
516  chiSquare = 0;
517  }
518  return chiSquare;
519 }
520 
521 double
522 PndRiemannTrack::calcSZChi2(PndRiemannHit* hit){
523  // get s'es and zs
524  TF1* f = 0;
525  if (fFitDone == false)
526  refit ();
527  if (r() > 0) {
528  unsigned int num = getNumHits();
529  if (fVerbose > 1)
530  std::cout << "szChi2Fit(hit) for " << num + 1 << " Points!" << std::endl;
531  TGraph g(num + 1);
532  int j = 0;
533  for (unsigned int i = 0; i < num; ++i) {
534  if (fHits[i].hit() != 0 && fHits[i].hit()->GetEntryNr().GetType() == GetBranchId("STTHit")){
535  continue;
536  }
537  fHits[i].calcPosOnTrk(this);
538  if (fVerbose > 1){
539 
540  std::cout << "Point: " << j << ": ";
541  if (fHits[i].hit() != 0)
542  std::cout << fHits[i].hit()->GetEntryNr() << " ";
543  std::cout << fHits[i].s() << " " << fHits[i].z() << std::endl;
544 
545  }
546  g.SetPoint(j, fHits[i].s(), fHits[i].z());
547  j++;
548  }
549  if (fVerbose > 1)
550  std::cout << "Additional hit: ";
551  hit->calcPosOnTrk(this);
552  if (fVerbose > 1)
553  std::cout << hit->s() << " " << hit->z() << std::endl;
554  g.SetPoint(j, hit->s(), hit->z());
555  g.Set(j+1);
556  g.Fit("pol1", "Q0");
557  f = g.GetFunction("pol1");
558 
559  // fChi2 = chis;
560  if (fVerbose > 1)
561  std::cout << "t, m: " << f->GetParameter(0) << " +/- "
562  << f->GetParError(0) << " / " << f->GetParameter(1)
563  << " +/- " << f->GetParError(1) << "Chi2: "
564  << f->GetChisquare() << std::endl;
565 
566  return f->GetChisquare() / f->GetNDF();
567  // delete(f);
568  } else {
569  std::cout << "-W- PndRiemannTrack::calcSZChi2 r == 0: " << *this << std::endl;
570  return -1;
571  }
572 }
573 
574 double PndRiemannTrack::calcZPosByS(double s)
575 {
576 // TVectorD o=orig();
577 // const PndRiemannHit* firstHit=getHit(0);
578 // assert(firstHit!=NULL);
579 // TVector2 k(firstHit->x().X()-o[0],firstHit->x().Y()-o[1]);
580 // TVector2 l = k.Rotate(s/r());
581  Double_t zCoord = s*m()+t();
582 // TVector3 result(l.X()+o[0], l.Y()+o[1], zCoord);
583 
584 // if (fVerbose > 0) std::cout<< "PosByS: s:" << s << " Vector: " << result.x() << " " << result.y() << " " << result.z() << std::endl;
585 // return result;
586  return zCoord;
587 }
588 
589 TVector3 PndRiemannTrack::calcPosByS(double s)
590 {
591 // calcSForHits();
592  TVectorD o=orig();
593  const PndRiemannHit* firstHit=getHit(0);
594 
595 
596 
597  TVector2 k(firstHit->x().X()-o[0],firstHit->x().Y()-o[1]);
598 
599  Double_t start_phi = k.Phi();
600  Double_t delta_phi = s / r();
601  TVector2 Res2D(r(), 0);
602 
603  TVector2 rotated = Res2D.Rotate(start_phi + delta_phi);
604 
605  TVector3 result(rotated.X() + o[0], rotated.Y() + o[1], calcZPosByS(s));
606 
607  return result;
608 }
609 
610 int PndRiemannTrack::calcIntersection(PndRiemannTrack& track, TVector3& p1, TVector3& p2)
611 {
612  const double SMALL = 1E-12;
613  const double VERTEX_CUT = fVertexCut;
614  if (track.getNumHits() < 3){
615  if (fVerbose > 1) std::cout << "-I- PndRiemannTrack::clacIntersection: less than 3 hits in track!" << std::endl;
616  return 0;
617  }
618 
619  TVector3 n1(fn[0],fn[1],fn[2]); //normal vector of THIS plane
620  TVectorD tempn = track.n();
621  TVector3 n2(tempn[0], tempn[1], tempn[2]); //normal vector of track plane
622  Double_t c1 = c(); //offset of THIS plane
623  Double_t c2 = track.c(); //offset of track plane
624 
625  if(fVerbose > 1) std::cout << "n1: " << n1.X() << " " << n1.Y() << " " << n1.Z() << " c1: " << c1 << std::endl;
626  if(fVerbose > 1) std::cout << "n2: " << n2.X() << " " << n2.Y() << " " << n2.Z() << " c2: " << c2 << std::endl;
627 
628  TVector3 lineNorm = n1.Cross(n2);
629 
630  TVector3 dLineNorm = calcErrorLineNorm(track);
631 
632  if (lineNorm.Mag() < SMALL)
633  return 0;
634  //lineNorm.SetMag(1.); //normalized normal vector of intersection line between the two planes
635  if(fVerbose > 1) std::cout << "n1 x n2: " << lineNorm.X() << "+/-" << dLineNorm.X() << " "
636  << lineNorm.Y() << "+/-" << dLineNorm.Y() << " "
637  << lineNorm.Z() << "+/-" << dLineNorm.Z() << std::endl;
638 
639  double denom = n1[0]*n2[1]-n1[1]*n2[0];
640  double x0 = (n1[1]*c2-c1*n2[1])/denom;
641  double y0 = (c1*n2[0] - n1[0]*c2)/denom;
642  TVector3 lineOffset(x0, y0, 0); //point on the intersection line
643  TVector3 dLineOffset = calcErrorLineOffset(track);
644  if(fVerbose > 1) std::cout << "lineOffset: " << lineOffset.X() << "+/-" << dLineOffset.X() << " "
645  << lineOffset.Y() << "+/-" << dLineOffset.Y() << std::endl;
646 
647  //---------- Calculation of intersection of line with parabola ------------
648  double denom2 = lineNorm[0]*lineNorm[0]+lineNorm[1]*lineNorm[1];
649  double p = (2*(lineNorm[0]*lineOffset[0]+lineNorm[1]*lineOffset[1]) - lineNorm[2])/denom2;
650  double q = (lineOffset[0]*lineOffset[0]+lineOffset[1]*lineOffset[1] - lineOffset[2])/denom2;
651 
652  double rad = (p/2)*(p/2) - q;
653  if (rad < 0){
654  if (fVerbose > 1)
655  std::cout << "PndRiemannTrack::calcIntersection: No match with parabola" << std::endl;
656  return 0;
657  }
658  double l1 = -(p/2)+TMath::Sqrt(rad);
659  double l2 = -(p/2)-TMath::Sqrt(rad);
660 
661  TVectorD dXY12 = calcErrorXY1XY2(lineNorm, dLineNorm, lineOffset, dLineOffset);
662 
663  if (fVerbose > 1)
664  std::cout << "l1: " << l1 << "+/-" << dXY12[4] << " l2: " << l2 << "+/-" << dXY12[5] << std::endl;
665 
666  TVector3 inter1 = lineNorm * l1 + lineOffset; //intersection point1 of line with parabola
667  TVector3 inter2 = lineNorm * l2 + lineOffset; //intersection point2 of line with parabola
668  if (fVerbose > 1){
669  std::cout << "intersection parabola 1: " << inter1.X() << "+/-" << dXY12[0] << " " << inter1.Y() << "+/-" << dXY12[1] << " " << inter1.Z()<< std::endl;
670  std::cout << "intersection parabola 2: " << inter2.X() << "+/-" << dXY12[2] << " " << inter2.Y() << "+/-" << dXY12[3] << " " << inter2.Z()<< std::endl;
671  }
672  PndRiemannHit hit1, hit2;
673  hit1.setXYZ(inter1.X(), inter1.Y(), inter1.Z());
674  hit2.setXYZ(inter2.X(), inter2.Y(), inter2.Z());
675 
676  hit1.calcPosOnTrk(this);
677  double s1 = hit1.s(); //arc length of hit1 for THIS track
678  TVector2 dummy1(inter1.X(), inter2.Y());
679  TVector2 dummy2(dXY12[0], dXY12[1]);
680  double dS1 = calcErrorS(dummy1, dummy2, this);
681  if (fVerbose > 1) std::cout << "S1 for track1: " << s1 << "+/-" << dS1 << std::endl;
682 
683  hit1.calcPosOnTrk(&track);
684  double s2 = hit1.s(); //arc length of hit1 for track
685  double dS2 = calcErrorS(dummy1, dummy2, &track);
686  if (fVerbose > 1) std::cout << "S1 for track2: " << s2 << "+/-" << dS2 << std::endl;
687 
688 
689  TVector3 vertex1;// = calcPosByS(s1); //backcalculated vertex position from s1 with z-Value
690  TVector3 vertex2;// = track.calcPosByS(s2); //backcalculated vertex position from s2 with z-Value
691  TVector3 dVertex1 = calcErrorPosByS(s1, dS1);
692  TVector3 dVertex2 = track.calcErrorPosByS(s2, dS2);
693 
694  vertex1.SetXYZ(inter1.X(),inter1.Y(),calcZPosByS(s1)); // skip calculation error (from (Xv,Yv) to S and then back from S to (Xv,Yv))
695  vertex2.SetXYZ(inter1.X(),inter1.Y(),track.calcZPosByS(s2));
696 
697  if (fVerbose > -1){
698  std::cout << "vertex candidate1 for hit1: " << vertex1.X() << " " << vertex1.Y() << " " << vertex1.Z()<< std::endl;
699  std::cout << "vertex candidate2 for hit1: " << vertex2.X() << " " << vertex2.Y() << " " << vertex2.Z()<< std::endl;
700  }
701 
702  if ((vertex1-vertex2).Mag() < VERTEX_CUT){
703  if (fVerbose > -1) std::cout << "Vertex Found!" << std::endl;
704  p1 = vertex1;
705  p2 = vertex2;
706  return 1;
707  }
708 
709  hit2.calcPosOnTrk(this);
710  double s1b = hit2.s();
711  dummy1.Set(inter2.X(),inter2.Y());
712  dummy2.Set(dXY12[2],dXY12[3]);
713  double dS1b = calcErrorS(dummy1, dummy2, this);
714  if (fVerbose > 1) std::cout << "S1b for track1: " << s1b << "+/-" << dS1b << std::endl;
715 
716  hit2.calcPosOnTrk(&track);
717  double s2b = hit2.s();
718  double dS2b = calcErrorS(dummy1, dummy2, &track);
719  if (fVerbose > 1) std::cout << "S1b for track2: " << s2b << "+/-" << dS2b << std::endl;
720 
721 // vertex1 = calcPosByS(s1b);
722 // vertex2 = track.calcPosByS(s2b);
723 
724  vertex1.SetXYZ(inter2.X(),inter2.Y(),calcZPosByS(s1b)); // skip calculation error (from (Xv,Yv) to S and then back from S to (Xv,Yv))
725  vertex2.SetXYZ(inter2.X(),inter2.Y(),track.calcZPosByS(s2b));
726 
727  dVertex1 = calcErrorPosByS(s1b, dS1b);
728  dVertex2 = track.calcErrorPosByS(s2b, dS2b);
729 
730  if (fVerbose > -1){
731  std::cout << "vertex candidate1 for hit2: " << vertex1.X() << " " << vertex1.Y() << " " << vertex1.Z()<< std::endl;
732  std::cout << "vertex candidate2 for hit2: " << vertex2.X() << " " << vertex2.Y() << " " << vertex2.Z()<< std::endl;
733  }
734  if ((vertex1-vertex2).Mag() < VERTEX_CUT){
735  if (fVerbose > -1) std::cout << "Vertex Found!" << std::endl;
736  p1 = vertex1;
737  p2 = vertex2;
738  return 1;
739  }
740  return 0;
741 }
742 
743 TVector3 PndRiemannTrack::calcErrorLineNorm(PndRiemannTrack& track)
744 {
745  TVectorD n1=n();
746  TVectorD n2=track.n();
747 
748  TMatrixD covN1 = covPlane();
749  TMatrixD covN2 = track.covPlane();
750 
751  TVector3 dN1(TMath::Sqrt(covN1[1][1]), TMath::Sqrt(covN1[2][2]), TMath::Sqrt(covN1[3][3]));
752  TVector3 dN2(TMath::Sqrt(covN2[1][1]), TMath::Sqrt(covN2[2][2]), TMath::Sqrt(covN2[3][3]));
753 
754  Double_t dLineNorm1 = TMath::Sqrt(TMath::Power(dN1[1]*n2[2],2) + TMath::Power(n1[1]*dN2[2],2) +
755  TMath::Power(dN1[2]*n2[1],2) + TMath::Power(n1[2]*dN2[1],2));
756 
757  Double_t dLineNorm2 = TMath::Sqrt(TMath::Power(dN1[2]*n2[0],2) + TMath::Power(n1[2]*dN2[0],2) +
758  TMath::Power(dN1[0]*n2[2],2) + TMath::Power(n1[0]*dN2[2],2));
759 
760  Double_t dLineNorm3 = TMath::Sqrt(TMath::Power(dN1[0]*n2[1],2) + TMath::Power(n1[0]*dN2[1],2) +
761  TMath::Power(dN1[1]*n2[0],2) + TMath::Power(n1[1]*dN2[0],2));
762 
763  return TVector3(dLineNorm1, dLineNorm2, dLineNorm3);
764 
765 }
766 
767 TVector3 PndRiemannTrack::calcErrorLineOffset(PndRiemannTrack& track)
768 {
769  TVectorD n1=n();
770  Double_t c1=c();
771  TVectorD n2=track.n();
772  Double_t c2=track.c();
773 
774  TMatrixD covN1 = covPlane();
775  TMatrixD covN2 = track.covPlane();
776 
777  TVector3 dN1(TMath::Sqrt(covN1[1][1]), TMath::Sqrt(covN1[2][2]), TMath::Sqrt(covN1[3][3]));
778  TVector3 dN2(TMath::Sqrt(covN2[1][1]), TMath::Sqrt(covN2[2][2]), TMath::Sqrt(covN2[3][3]));
779  Double_t dC1(TMath::Sqrt(covN1[0][0]));
780  Double_t dC2(TMath::Sqrt(covN2[0][0]));
781 
782  Double_t denom = n1[0]*n2[1] - n1[1]*n2[0];
783 // Double_t dDenom2 = (TMath::Power(dN1[0]*n2[1],2) + TMath::Power(n1[0]*dN2[1],2) +
784 // TMath::Power(dN1[1]*n2[0],2) + TMath::Power(n1[1]*dN2[0],2));
785 
786  Double_t qX = (n1[1]*c2-n2[1]*c1)/(TMath::Power(denom,2));
787  Double_t qY = (n2[0]*c1-n1[0]*c2)/(TMath::Power(denom,2));
788 
789  Double_t dx = TMath::Sqrt(TMath::Power(n2[1]*dC1/denom,2) + TMath::Power(n1[1]*dC2/denom,2) +
790  TMath::Power(qX*n2[1]*dN1[0],2) + TMath::Power(qX*n1[1]*dN2[0],2) +
791  TMath::Power((c1/denom + qX*n2[0])*dN1[1],2) +
792  TMath::Power((c1/denom + qX*n1[0])*dN2[1],2));
793 
794  Double_t dy = TMath::Sqrt(TMath::Power(n2[0]*dC1/denom,2) + TMath::Power(n1[0]*dC2/denom,2) +
795  TMath::Power(qX*n2[0]*dN1[1],2) + TMath::Power(qX*n1[0]*dN2[1],2) +
796  TMath::Power((c1/denom + qY*n2[1])*dN1[0],2) +
797  TMath::Power((c1/denom + qY*n1[1])*dN2[0],2));
798 
799  return TVector3(dx, dy, 0);
800 
801 }
802 
803 TVectorD PndRiemannTrack::calcErrorXY1XY2(TVector3& line, TVector3& dLine, TVector3& offset, TVector3& dOffset)
804 {
805  const double SMALL = 1E-12;
806  Double_t denom = (TMath::Power(line[0],2) + TMath::Power(line[1],2));
807  Double_t p = (2*(offset[0]*line[0] + offset[1]*line[1]) - line[2]) / denom;
808  Double_t q = (TMath::Power(offset[0],2)+TMath::Power(offset[1],2) - offset[2]) / denom;
809  Double_t S1 = -(p/2) + TMath::Sqrt(TMath::Power(p/2,2) - q);
810  Double_t S2 = -(p/2) - TMath::Sqrt(TMath::Power(p/2,2) - q);
811 
812 // std::cout << "S1: " << S1 << " S2: " << S2 << std::endl;
813  Double_t prod = (2*(offset[0]*line[0] + offset[1]*line[1]) - line[2]) / TMath::Power(denom,2);
814 // std::cout << "denom: " << denom << " prod1: " << prod * line[0] * dLine[0] << " prod2: " << prod * line[1] * dLine[1]<< std::endl;
815 
816  if ((TMath::Power(p/2,2)-q) < SMALL)
817  return TVectorD(6);
818 
819  Double_t sqrterm = TMath::Sqrt(TMath::Power(p/2,2)-q);
820 
821  Double_t dP = TMath::Sqrt(TMath::Power(2*line[0]*dOffset[0]/denom,2) + TMath::Power(2*line[1]*dOffset[1]/denom,2) +
822  TMath::Power(dLine[2]/denom,2) +
823  TMath::Power((2*offset[0]/denom - 2*prod*line[0])*dLine(0),2) +
824  TMath::Power((2*offset[1]/denom - 2*prod*line[1])*dLine(1),2));
825 // std::cout << "dP1: " << TMath::Power(2*line[0]*dOffset[0]/denom,2) << " dP2: " << TMath::Power(2*line[1]*dOffset[1]/denom,2) << std::endl;
826 // std::cout << "dP3: " << TMath::Power(dLine[2]/denom,2) << std::endl;
827 // std::cout << "dP4: " << TMath::Power((2*offset[0]/denom - 2*prod*line[0])*dLine(0),2) << " dP5: " << TMath::Power((2*offset[1]/denom - 2*prod*line[1])*dLine(1),2) << std::endl;
828 
829 // std::cout << "P: " << p << "+/-" << dP << std::endl;
830  Double_t dQ = TMath::Sqrt(TMath::Power(2*offset[0]/denom*dOffset[0],2) +
831  TMath::Power(2*offset[1]/denom*dOffset[1],2) +
832  TMath::Power(dOffset[2]/denom,2));
833 // std::cout << "Q: " << q << "+/-" << dQ << std::endl;
834  Double_t dS1 = TMath::Sqrt(TMath::Power((-1/2 + p/(4*sqrterm)) * dP,2) + TMath::Power(1/(2*sqrterm)*dQ, 2));
835  Double_t dS2 = TMath::Sqrt(TMath::Power((-1/2 - p/(4*sqrterm)) * dP,2) + TMath::Power(1/(2*sqrterm)*dQ, 2));
836 
837  Double_t dX1 = TMath::Sqrt(TMath::Power(dOffset[0],2) + TMath::Power(line[0]*dS1,2) + TMath::Power(dLine[0]*S1,2));
838  Double_t dY1 = TMath::Sqrt(TMath::Power(dOffset[1],2) + TMath::Power(line[1]*dS1,2) + TMath::Power(dLine[1]*S1,2));
839 
840 // std::cout << "dX1: " << dX1 << " dY1: " << dY1 << std::endl;
841 
842  Double_t dX2 = TMath::Sqrt(TMath::Power(dOffset[0],2) + TMath::Power(line[0]*dS2,2) + TMath::Power(dLine[0]*S2,2));
843  Double_t dY2 = TMath::Sqrt(TMath::Power(dOffset[1],2) + TMath::Power(line[1]*dS2,2) + TMath::Power(dLine[1]*S2,2));
844 
845 // std::cout << "dX2: " << dX2 << " dY2: " << dY2 << std::endl;
846  TVectorD result(6);
847  result[0] = dX1;
848  result[1] = dY1;
849  result[2] = dX2;
850  result[3] = dY2;
851  result[4] = dS1;
852  result[5] = dS2;
853  return result;
854 
855 }
856 
857 double PndRiemannTrack::calcErrorS(TVector2& XY, TVector2& dXY, PndRiemannTrack* track)
858 {
859  assert(track!=NULL);
860  TVectorD o=track->orig();
861  double R=track->r();
862  double dr = track->dR();
863 
864  const PndRiemannHit* firstHit=track->getHit(0);
865  assert(firstHit!=NULL);
866  TVector2 k(firstHit->x().X()-o[0],firstHit->x().Y()-o[1]);
867  TVector2 l(XY.X()-o[0],XY.Y()-o[1]);
868 
869  double alpha=l.DeltaPhi(k);
870  double dAlpha = TMath::Sqrt(TMath::Power(XY.X()*dXY.Y(),2) + TMath::Power(XY.Y()*dXY.X(),2))/(TMath::Power(XY.X(),2) + TMath::Power(XY.Y(),2));
871 
872  double dS = TMath::Sqrt(TMath::Power(dAlpha*R,2) + TMath::Power(alpha*dr,2));
873 
874  return dS;
875 
876 }
877 
878 TVector3 PndRiemannTrack::calcErrorPosByS(Double_t s, Double_t dS)
879 {
880  TVectorD o=orig();
881  const PndRiemannHit* firstHit=getHit(0);
882  assert(firstHit!=NULL);
883  TVector2 k(firstHit->x().X()-o[0],firstHit->x().Y()-o[1]);
884  Double_t dPhi = TMath::Sqrt(TMath::Power(dS/r(),2) + TMath::Power(s/(r()*r())*dR(),2));
885  Double_t phi = s/r();
886  Double_t dL1 = TMath::Sqrt(TMath::Power((-k.X()*TMath::Sin(phi)-k.Y()*TMath::Cos(phi))*dPhi,2));
887  Double_t dL2 = TMath::Sqrt(TMath::Power((k.X()*TMath::Cos(phi)-k.Y()*TMath::Sin(phi))*dPhi,2));
888  Double_t dZCoord = TMath::Sqrt(TMath::Power(m()*dS,2) + TMath::Power(s*mError(),2) + TMath::Power(tError(),2));
889 
890  Double_t resX = TMath::Sqrt(TMath::Power(dL1,2) + TMath::Power(dX(),2));
891  Double_t resY = TMath::Sqrt(TMath::Power(dL2,2) + TMath::Power(dY(),2));
892 
893  TVector3 result(resX, resY, dZCoord);
894 
895  if (fVerbose > 1) std::cout<< "PosBySError: s:" << dS << " Vector: " << result.x() << " " << result.y() << " " << result.z() << std::endl;
896  return result;
897 }
898 
899 double
900 PndRiemannTrack::szDist(PndRiemannHit* hit){
901 
902  if (fSZFitDone == false)
903  szFit();
904  hit->calcPosOnTrk(this);
905  double hits=hit->s();
906  double predz=hits*fm+ft;
907  return predz-hit->z();
908 }
909 
910 double PndRiemannTrack::szError(PndRiemannHit* hit){
911  if (fSZFitDone == false)
912  szFit();
913  hit->calcPosOnTrk(this);
914  double hits=hit->s();
915  double result = TMath::Sqrt(TMath::Power((hits*fmError),2)+TMath::Power(ftError,2));
916  //if (fVerbose > 0) std::cout << "Error on szDist is: " << result;
917  return result;
918 }
919 
920 void PndRiemannTrack::calcStartStopAlpha()
921 {
922  PndRiemannHit* first = getHit(0);
923  PndRiemannHit* last = getLastHit();
924 
925  fStartAlpha = calcAlpha(first);
926  fStopAlpha = calcAlpha(last);
927 
928 }
929 
930 double PndRiemannTrack::calcAlpha(PndRiemannHit* myHit)
931 {
932  TVectorD origin = orig();
933  TVector2 myVector(myHit->x().x() - origin[0], myHit->x().y() - origin[1]);
934  return myVector.Phi();
935 }
936 
937 void PndRiemannTrack::calcSForHits(){
938  for (size_t i = 0; i < fHits.size(); i++){
939  fHits[i].calcPosOnTrk(this);
940  }
941  sortHits();
942 }
943 
944 // only after szFit!
945 double
946 PndRiemannTrack::dip() {
947  if (fSZFitDone == false)
948  szFit();
949  return cos(atan(fm));
950 }
951 
952 // only after szFit!
953 double
954 PndRiemannTrack::dipangle() {
955  if (fSZFitDone == false)
956  szFit();
957  return (atan(fm));
958 }
959 
960 double PndRiemannTrack::dDip()
961 {
962  if (fSZFitDone == false)
963  szFit();
964  //return (fabs(sin(1/(1+fm*fm))) * fmError);
965  //WRONG derivative!!!!!!!
966  return (fabs(sin(atan(fm))/(1+fm*fm)) * fmError);
967 }
968 
969 double PndRiemannTrack::Pt(double B)
970 {
971  double result = 0.3 * B * r() * 0.01;
972  if (fVerbose > 0) std::cout << "Pt: " << result << std::endl;
973 
974  return result;
975 }
976 
977 double PndRiemannTrack::Pl(double B)
978 {
979  double pl = Pt(B) * m();
980 
981  if (getNumHits() > 2){
982  PndRiemannHit* startHit = getHit(0);
983  PndRiemannHit* nextHit = getHit(1);
984  if (nextHit->z() - startHit->z() > 0){
985  if (pl < 0)
986  pl *= -1;
987  }
988  else{
989  if (pl > 0){
990  pl *= -1;
991  }
992  }
993  }
994  return pl;
995 }
996 
997 double PndRiemannTrack::P(double B)
998 {
999  double result = TMath::Sqrt(Pl(B) * Pl(B) + Pt(B) * Pt(B));
1000  return result;
1001 }
1002 
1003 TVector3 PndRiemannTrack::getPforHit(int i, double B)
1004 {
1005  double pt = Pt(B);
1006  //double p = P(B); //[R.K. 01/2017] unused variable?
1007  double pl = Pl(B);
1008  TVectorD origin = orig();
1009  TVector3 result;
1010 // PndRiemannHit* startHit = getHit(0);
1011 // PndRiemannHit* lastHit = getLastHit();
1012 // if (lastHit->z() < startHit->z()){
1013 // pl *= -1;
1014 // }
1015 
1016 
1017  if (i < (int)getNumHits())
1018  {
1019  PndRiemannHit* myHit = getHit(i);
1020  //std::cout << "MyHit: " << myHit->x().X() << " " << myHit->x().Y() << std::endl;
1021  Double_t arclength = myHit->s()/r();
1022 
1023  Double_t phi = fStartAlpha + arclength;
1024  TVector2 ptVec(0,pt);
1025  ptVec = ptVec.Rotate(phi);
1026  if (i > 0){
1027  PndRiemannHit* myHitBefore = getHit(i-1);
1028  //std::cout << "MyHitBefore: " << myHitBefore->x().X() << " " << myHitBefore->x().Y() << std::endl;
1029  TVector2 difVec(myHit->x().x() - myHitBefore->x().x(), myHit->x().y() - myHitBefore->x().y());
1030  //std::cout <<"DifVec: " << difVec.X() << " " << difVec.Y() << std::endl;
1031  Double_t length = ptVec * difVec;
1032  //std::cout << "Length: " << length << std::endl;
1033  if (length < 0){
1034  ptVec *= -1;
1035  }
1036  }
1037  else if(getNumHits() > 0){
1038  PndRiemannHit* myHitAfter = getHit(1);
1039  //std::cout << "MyHitAfter: " << myHitAfter->x().X() << " " << myHitAfter->x().Y() << std::endl;
1040  TVector2 difVec2(myHitAfter->x().x() - myHit->x().x(), myHitAfter->x().y() - myHit->x().y());
1041  //std::cout <<"DifVec: " << difVec2.X() << " " << difVec2.Y() << std::endl;
1042  Double_t length2 = ptVec * difVec2;
1043  //std::cout << "Length: " << length2 << std::endl;
1044  if (length2 < 0){
1045  ptVec *= -1;
1046  }
1047  }
1048  result.SetXYZ(ptVec.X(), ptVec.Y(), pl);
1049  }
1050 
1051  if (fVerbose > 0) std::cout << "P-Vector for point " << i << " : " << result.X() << " " << result.Y() << " " << result.Z() << std::endl;
1052  return result;
1053 }
1054 
1055 Int_t PndRiemannTrack::getCharge(Double_t B)
1056 {
1057  TVector3 p = getPforHit(0, B);
1058  Double_t s = getHit(0)->s();
1059  TVector3 hitPos = calcPosByS(s);
1060  TVector2 hitPos2D(hitPos.x(), hitPos.y());
1061 
1062  TVector2 pt(p.x(), p.y());
1063  TVectorD origin = orig();
1064  TVector2 orig2(origin[0], origin[1]);
1065 
1066  orig2 -= hitPos2D;
1067 
1068  TVector2 origRotated = orig2.Rotate(-pt.Phi());
1069 
1070  //std::cout << "OrigRotated: " << origRotated.X() << " " << origRotated.Y() << std::endl;
1071 
1072  if (origRotated.Y() < 0)
1073  return 1;
1074  else
1075  return -1;
1076 }
1077 
1078 FairTrackParP PndRiemannTrack::getTrackParPForHit(Int_t i, Double_t B)
1079 {
1080  calcSForHits();
1081 
1082 
1083  TVector3 hitPos;
1084  TVector3 hitPosError;
1085  TVector3 momError(2, 2, 2);
1086  TVector3 dj(1,0,0);
1087  TVector3 dk(0,1,0);
1088  TVector3 origin(0, 0, 1);
1089 
1090  if (i < (int)getNumHits()){
1091  Double_t s = getHit(i)->s();
1092  getHit(i)->hit()->Position(hitPos);
1093  hitPos = calcPosByS(s);
1094  // getHit(i)->hit()->Position(hitPos);
1095  getHit(i)->hit()->PositionError(hitPosError);
1096  FairTrackParP result(hitPos, getPforHit(i, B), hitPosError, momError, getCharge(B), origin, dj, dk);
1097 
1098  return result;
1099  }
1100  else {
1101  return FairTrackParP();
1102  }
1103 }
1104 
1105 PndTrack PndRiemannTrack::getPndTrack(Double_t B)
1106 {
1107  FairTrackParP first, last;
1108  PndTrackCand myCand;
1109  if (getNumHits() > 0){
1110  first = getTrackParPForHit(0, B);
1111  last = getTrackParPForHit(getNumHits()-1, B);
1112  //FairTrackParP last;
1113  }
1114  return PndTrack(first, last, myCand);
1115 }
1116 
1117 // Todo: check for backward going tracks!!!
1118 double
1119 PndRiemannTrack::sign() const {
1120  double s=fHits[3].s();
1121  return s<0 ? -1. : 1;
1122 }
1123 
1124 void PndRiemannTrack::calcJacRXY()
1125 {
1126  double val = -1.;
1127  if ((1-fn[2]*fn[2]-4*fc*fn[2]) > 0)
1128  val = (TMath::Sqrt(1-fn[2]*fn[2]-4*fc*fn[2]));
1129  else{
1130  std::cout << "-E- PndRiemannTrack::calcJacRXY val is imaginary: Sqrt of " << (1-fn[2]*fn[2]-4*fc*fn[2]) << std::endl;
1131  return;
1132  }
1133  fjacRXY.Clear();
1134  fjacRXY.ResizeTo(3,4);
1135 /* fjacRXY[0][0] = -1/val;
1136  fjacRXY[0][3] = -1/(2*fn[2]) * (fn[2]+2*fc)/val - val/(2*fn[2]*fn[2]);
1137  fjacRXY[1][1] = -fn[0]/(2*fn[2]);
1138  fjacRXY[1][3] = fn[0]/(2*fn[2]*fn[2]);
1139  fjacRXY[2][2] = -fn[1]/(2*fn[2]);
1140  fjacRXY[2][3] = fn[1]/(2*fn[2]*fn[2]);*/
1141 
1142  if (fn[2] != 0.){
1143  fjacRXY[0][0] = -1/val;
1144  fjacRXY[0][3] = -1/(2*fn[2]) * (fn[2]+2*fc)/val - val/(2*fn[2]*fn[2]);
1145  fjacRXY[1][1] = -1/(2*fn[2]);//<---------
1146  fjacRXY[1][3] = fn[0]/(2*fn[2]*fn[2]);
1147  fjacRXY[2][2] = -1/(2*fn[2]);//<---------
1148  fjacRXY[2][3] = fn[1]/(2*fn[2]*fn[2]);
1149  }
1150  else
1151  std::cout << "-E- PndRiemannTrack::calcJacRXY fn[2] is zero!" << std::endl;
1152 }
1153 
1154 void PndRiemannTrack::PrintHits()
1155 {
1156  std::cout << "-I- PndRiemannTrack::PrintHits: " << fHits.size() << std::endl;
1157  //bool first = false; //[R.K. 01/2017] unused variable?
1158  for (size_t i = 0; i < fHits.size(); i++){
1159  std::cout << i << ": ";
1160 
1161  if (fHits[i].hit() != 0) {
1162 // std::cout << fHits[i].hit()->GetEntryNr() << " : ";
1163  }
1164  std::cout << fHits[i].x().X() << " +/- " << fHits[i].sigmaX()
1165  << "\t" << fHits[i].x().Y()<< " +/- " << fHits[i].sigmaY()
1166  << "\t" << fHits[i].z() << " s: " << fHits[i].s() << "\t" << " dist: " << dist(&fHits[i]);
1167  // calcPosByS(fHits[i].s()).Print();
1168  std::cout << std::endl;
1169 
1170  }
1171 }
1172 
1173 
1174 PndRiemannHit PndRiemannTrack::correctSttHit(PndSttHit* mySttHit)
1175 {
1176 
1177  TVectorD normalVector = n();
1178  if (fVerbose > 0)
1179  std::cout << "MySttHit: " << mySttHit << std::endl;
1180  Double_t x_val = mySttHit->GetX();
1181  Double_t y_val = mySttHit->GetY();
1182 
1183  Double_t phi = TMath::ATan2((normalVector[1]+2*normalVector[2]*y_val), (normalVector[0]+2*normalVector[2]*x_val));
1184 
1185  Double_t correctedX = x_val+TMath::Cos(phi)*mySttHit->GetIsochrone();
1186  Double_t correctedY = y_val+TMath::Sin(phi)*mySttHit->GetIsochrone();
1187 
1188  Double_t correctedX1 = x_val+TMath::Cos(phi+TMath::Pi())*mySttHit->GetIsochrone();
1189  Double_t correctedY1 = y_val+TMath::Sin(phi+TMath::Pi())*mySttHit->GetIsochrone();
1190 
1191  PndRiemannHit myRiemannHit(mySttHit);
1192  PndRiemannHit myRiemannHit1(mySttHit);
1193  myRiemannHit1.setXYZ(correctedX, correctedY, 0);
1194  PndRiemannHit myRiemannHit2(mySttHit);
1195  myRiemannHit2.setXYZ(correctedX1, correctedY1, 0);
1196 
1197  Double_t finalX, finalY;
1198 
1199  if (fabs(dist(&myRiemannHit)) < fabs(dist(&myRiemannHit1))){
1200  if (fabs(dist(&myRiemannHit)) < fabs(dist(&myRiemannHit2))){
1201  finalX = x_val;
1202  finalY = y_val;
1203  }
1204  else {
1205  finalX = correctedX1;
1206  finalY = correctedY1;
1207  }
1208  } else {
1209  if (fabs(dist(&myRiemannHit1)) < fabs(dist(&myRiemannHit2))){
1210  finalX = correctedX;
1211  finalY = correctedY;
1212  } else {
1213  finalX = correctedX1;
1214  finalY = correctedY1;
1215  }
1216 
1217  }
1218 
1219  if (fVerbose > 1){
1220  std::cout << "TrackNormal + Length: " << normalVector[0] << "/" << normalVector[1] << "/" << normalVector[2] << " o: " << c() << std::endl;
1221  std::cout << "CalculateCorrectedSTTHit: TubeCenter(" << x_val << "/" << y_val << ") Radius: " << mySttHit->GetIsochrone() << " Distance: " << dist(&myRiemannHit) << std::endl;
1222  std::cout << "CalculateCorrectedSTTHit: Corrected1(" << correctedX << "/" << correctedY << ") Distance: " << dist(&myRiemannHit1) << " phi: " << phi << std::endl;
1223  std::cout << "CalculateCorrectedSTTHit: Corrected2(" << correctedX1 << "/" << correctedY1 << ") Distance: " << dist(&myRiemannHit2) << " phi: " << phi << std::endl;
1224 
1225  std::cout << "Result: " << finalX << "/" << finalY << std::endl;
1226  }
1227  PndRiemannHit result(mySttHit);
1228  if (fVerbose > 1){
1229  std::cout << "result: " << mySttHit << " " << result.x().X() << "/" << result.x().Y() << "/" << result.x().Z() << " " << result.z() << std::endl;
1230  }
1231  result.setXYZ(finalX, finalY, 0);
1232  result.setDXYZ(0.01,0.01,200);
1233 
1234  if (fVerbose > 1){
1235  std::cout << "result: " << result.x().X() << "/" << result.x().Y() << "/" << result.x().Z() << " " << result.z() << std::endl;
1236  }
1237  return result;
1238 }
1239 
1240 
1241 PndRiemannHit PndRiemannTrack::correctSttSkewedHit(PndSttHit* mySttHit, PndSttTube* myTube)
1242 {
1243  TVectorD nVec = n();
1244  if (fVerbose > 0)
1245  std::cout << "MySttHit: " << mySttHit << std::endl;
1246 
1247  TVector3 tubeDir = myTube->GetWireDirection();
1248  TVector3 tubeMidPos = myTube->GetPosition();
1249 
1250  Double_t nominator = nVec[0]*tubeDir.x() + nVec[1]*tubeDir.y() + 2*nVec[2]*(tubeDir.x() + tubeDir.y());
1251  Double_t denominator = 2*nVec[2]*myTube->GetHalfLength()*(tubeDir.x() * tubeDir.x() + tubeDir.y()*tubeDir.y());
1252 
1253  Double_t posOnTube = nominator/denominator;
1254  if (posOnTube > 1){
1255  posOnTube = 1;
1256  } else if (posOnTube < -1) {
1257  posOnTube = -1;
1258  }
1259 
1260  PndRiemannHit result(mySttHit);
1261  result.setXYZ(tubeMidPos.x()+posOnTube*myTube->GetHalfLength()*tubeDir.x(),
1262  tubeMidPos.y()+posOnTube*myTube->GetHalfLength()*tubeDir.y(),
1263  tubeMidPos.z()+posOnTube*myTube->GetHalfLength()*tubeDir.z());
1264  return result;
1265 }
1266 
1267 
1268 void PndRiemannTrack::correctSttHits()
1269 {
1270 
1271  for (size_t i = 0; i < fHits.size(); i++){
1272  if (fHits[i].hit()->GetEntryNr().GetType() == GetBranchId("STTHit")){
1273 
1274  PndRiemannHit myHit = fHits[i];
1275  PndSttHit* mySttHit = (PndSttHit*)fHits[i].hit();
1276 
1277  PndRiemannHit correctedHit = correctSttHit(mySttHit);
1278 // TVectorD normalVector = n();
1279 // Double_t phi = TMath::ATan2((normalVector[1]+2*normalVector[2]*myHit.x().Y()), (normalVector[0]+2*normalVector[2]*myHit.x().X()));
1280 //
1281 // Double_t correctedX = myHit.x().X()+TMath::Cos(phi)*mySttHit->GetIsochrone();
1282 // Double_t correctedY = myHit.x().Y()+TMath::Sin(phi)*mySttHit->GetIsochrone();
1283 //
1284 // Double_t correctedX1 = myHit.x().X()+TMath::Cos(phi+TMath::Pi())*mySttHit->GetIsochrone();
1285 // Double_t correctedY1 = myHit.x().Y()+TMath::Sin(phi+TMath::Pi())*mySttHit->GetIsochrone();
1286 //
1287 // PndRiemannHit myRiemannHit(myHit);
1288 // PndRiemannHit myRiemannHit1(myHit);
1289 // myRiemannHit1.setXYZ(correctedX, correctedY, 0);
1290 // PndRiemannHit myRiemannHit2(myHit);
1291 // myRiemannHit2.setXYZ(correctedX1, correctedY1, 0);
1292 //
1293 // Double_t finalX, finalY;
1294 //
1295 // if (fabs(dist(&myRiemannHit)) < fabs(dist(&myRiemannHit1))){
1296 // if (fabs(dist(&myRiemannHit)) < fabs(dist(&myRiemannHit2))){
1297 // finalX = myHit.x().X();
1298 // finalY = myHit.x().Y();
1299 // }
1300 // else {
1301 // finalX = correctedX1;
1302 // finalY = correctedY1;
1303 // }
1304 // } else {
1305 // if (fabs(dist(&myRiemannHit1)) < fabs(dist(&myRiemannHit2))){
1306 // finalX = correctedX;
1307 // finalY = correctedY;
1308 // } else {
1309 // finalX = correctedX1;
1310 // finalY = correctedY1;
1311 // }
1312 //
1313 // }
1314 //
1315 // if (fVerbose > 1){
1316 // std::cout << "TrackNormal + Length: " << normalVector[0] << "/" << normalVector[1] << "/" << normalVector[2] << " o: " << c() << std::endl;
1317 // std::cout << "CalculateCorrectedSTTHit: TubeCenter(" << myHit.x().X() << "/" << myHit.x().Y() << ") Radius: " << mySttHit->GetIsochrone() << " Distance: " << dist(&myRiemannHit) << std::endl;
1318 // std::cout << "CalculateCorrectedSTTHit: Corrected1(" << correctedX << "/" << correctedY << ") Distance: " << dist(&myRiemannHit1) << " phi: " << phi << std::endl;
1319 // std::cout << "CalculateCorrectedSTTHit: Corrected2(" << correctedX1 << "/" << correctedY1 << ") Distance: " << dist(&myRiemannHit2) << " phi: " << phi << std::endl;
1320 //
1321 // std::cout << "Result: " << finalX << "/" << finalY << std::endl;
1322 // }
1323  fHits[i].setXYZ(correctedHit.x().X(), correctedHit.x().Y(), 0);
1324  fHits[i].setDXYZ(0.001,0.001,200);
1325  }
1326  }
1327 
1328 // return result;
1329 }
1330 
1331 
1332 Int_t PndRiemannTrack::GetBranchId(TString branchName)
1333 {
1334  if (fBranchNameMap.count(branchName) == 0){
1335  fBranchNameMap[branchName] = FairRootManager::Instance()->GetBranchId(branchName);
1336  }
1337  return fBranchNameMap[branchName];
1338 }
row
int row
Definition: anaLmdDigi.C:67
pos
TVector3 pos
Definition: GammaSpectraAnalysis.C:43
acos
friend F32vec4 acos(const F32vec4 &a)
Definition: P4_F32vec4.h:113
x0
Double_t x0
Definition: checkhelixhit.C:70
PndRiemannTrack::getLastHit
PndRiemannHit * getLastHit()
Definition: PndRiemannTrack.h:76
PndRiemannTrack::getNumHits
unsigned int getNumHits()
Definition: PndRiemannTrack.h:74
p
Double_t p
Definition: anasim.C:58
fVerbose
int fVerbose
Definition: poormantracks.C:24
cos
friend F32vec4 cos(const F32vec4 &a)
Definition: P4_F32vec4.h:112
dy
double dy
Definition: runGeaneTestSimBox.cxx:45
PndRiemannTrack::fweight
double fweight
sum over all weights (1/(sigmaXY*sigmaXY))
Definition: PndRiemannTrack.h:167
PndRiemannTrack::calcSZChi2
double calcSZChi2(PndRiemannHit *hit)
Definition: PndRiemannTrack.cxx:522
d
TObjArray * d
Definition: Pnd_Hc_etaee7G.C:78
PndRiemannTrack::distError
double distError(PndRiemannHit *hit)
Definition: PndRiemannTrack.cxx:143
PndRiemannTrack::correctSttSkewedHit
PndRiemannHit correctSttSkewedHit(PndSttHit *mySttHit, PndSttTube *myTube)
Definition: PndRiemannTrack.cxx:1241
PndRiemannTrack::init
void init(double x0, double y0, double R, double dip, double z0)
Definition: PndRiemannTrack.cxx:87
res
Int_t res
Definition: anadigi.C:166
PndRiemannTrack::calcZPosByS
double calcZPosByS(double s)
Definition: PndRiemannTrack.cxx:574
num
int num[96]
Definition: ranlxd.cxx:381
i
Int_t i
Definition: run_full.C:25
PndRiemannTrack::correctSttHit
PndRiemannHit correctSttHit(PndSttHit *mySttHit)
Definition: PndRiemannTrack.cxx:1174
PndSttTube::GetHalfLength
Double_t GetHalfLength()
Definition: PndSttTube.cxx:99
PndRiemannHit::sigmaXY
double sigmaXY() const
Definition: PndRiemannHit.cxx:109
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
PndRiemannTrack::fHits
std::vector< PndRiemannHit > fHits
Definition: PndRiemannTrack.h:171
val
Double_t val[nBoxes][nFEBox]
Definition: createCalib.C:11
CAMath::Sqrt
static T Sqrt(const T &x)
Definition: PndCAMath.h:37
PndRiemannTrack::distCircle
double distCircle(PndRiemannHit *hit)
Definition: PndRiemannTrack.cxx:149
PndRiemannHit::hit
const FairHit * hit() const
Definition: PndRiemannHit.h:72
sin
friend F32vec4 sin(const F32vec4 &a)
Definition: P4_F32vec4.h:111
s
TLorentzVector s
Definition: Pnd2DStar.C:50
CAMath::Sin
static T Sin(const T &x)
Definition: PndCAMath.h:42
col
int col
Definition: anaLmdDigi.C:67
PndRiemannTrack::szDist
double szDist(PndRiemannHit *hit)
Definition: PndRiemannTrack.cxx:900
g
TFile * g
Definition: GammaSpectraAnalysis.C:22
c2
c2
Definition: plot_dirc.C:39
PndRiemannTrack::P
double P(double B)
Definition: PndRiemannTrack.cxx:997
PndRiemannTrack::dist
double dist(PndRiemannHit *hit)
Definition: PndRiemannTrack.cxx:129
offset
TVector3 offset(2, 0, 0)
PndRiemannTrack::GetBranchId
Int_t GetBranchId(TString branchName)
Definition: PndRiemannTrack.cxx:1332
CAMath::Cos
static T Cos(const T &x)
Definition: PndCAMath.h:43
PndRiemannHit::s
double s() const
Definition: PndRiemannHit.h:74
PndRiemannTrack::fmError
double fmError
Error of fit.
Definition: PndRiemannTrack.h:155
PndRiemannHit::setDXYZ
void setDXYZ(double dx, double dy, double dz)
Definition: PndRiemannHit.cxx:70
PndRiemannTrack::fBranchNameMap
std::map< TString, Int_t > fBranchNameMap
Definition: PndRiemannTrack.h:175
PndRiemannTrack::fChi2
double fChi2
Chisquare of sz fit.
Definition: PndRiemannTrack.h:157
PndRiemannTrack::calcJacRXY
void calcJacRXY()
calcualtes fjacRXY
Definition: PndRiemannTrack.cxx:1124
PndRiemannTrack::calcErrorPosByS
TVector3 calcErrorPosByS(Double_t s, Double_t dS)
Definition: PndRiemannTrack.cxx:878
PndRiemannTrack::sign
double sign() const
Definition: PndRiemannTrack.cxx:1119
PndRiemannTrack::calcErrorLineOffset
TVector3 calcErrorLineOffset(PndRiemannTrack &track)
Definition: PndRiemannTrack.cxx:767
pt
TString pt(TString pts, TString exts="px py pz")
Definition: invexp.C:133
CAMath::Abs
static T Abs(const T &x)
Definition: PndCAMath.h:39
PndRiemannTrack::n
const TVectorD & n() const
Definition: PndRiemannTrack.h:56
PndRiemannHit::sigmaW
double sigmaW() const
Definition: PndRiemannHit.h:80
PndRiemannTrack::calcErrorLineNorm
TVector3 calcErrorLineNorm(PndRiemannTrack &track)
Definition: PndRiemannTrack.cxx:743
a
Int_t a
Definition: anaLmdDigi.C:126
PndRiemannTrack::getTrackParPForHit
FairTrackParP getTrackParPForHit(Int_t i, Double_t B)
Definition: PndRiemannTrack.cxx:1078
PndRiemannTrack::m
double m() const
Definition: PndRiemannTrack.h:98
PndRiemannTrack::t
double t() const
Definition: PndRiemannTrack.h:100
PndRiemannTrack::tError
double tError() const
Definition: PndRiemannTrack.h:101
PndRiemannTrack::getCharge
Int_t getCharge(Double_t B)
Definition: PndRiemannTrack.cxx:1055
PndSttHit::GetIsochrone
Double_t GetIsochrone() const
Definition: PndSttHit.h:62
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
PndRiemannTrack::addPndTrackCand
void addPndTrackCand(PndTrackCand *trackCand)
Definition: PndRiemannTrack.cxx:70
PndRiemannTrack::fn
TVectorD fn
normal vector to plane;
Definition: PndRiemannTrack.h:149
y0
Double_t y0
Definition: checkhelixhit.C:71
CAMath::ATan2
static T ATan2(const T &y, const T &x)
PndRiemannHit::setXYZ
void setXYZ(double x, double y, double z)
Definition: PndRiemannHit.cxx:57
track
PndMCTrack * track
Definition: anaLmdCluster.C:89
PndSttTube::GetPosition
TVector3 GetPosition()
Definition: PndSttTube.cxx:87
PndRiemannHit::sigmaX
double sigmaX() const
Definition: PndRiemannHit.h:78
PndRiemannTrack::calcErrorXY1XY2
TVectorD calcErrorXY1XY2(TVector3 &line, TVector3 &dLine, TVector3 &offset, TVector3 &dOffset)
Definition: PndRiemannTrack.cxx:803
PndRiemannTrack::refit
void refit(bool withErrorCalc=true)
Definition: PndRiemannTrack.cxx:172
f
TFile * f
Definition: bump_analys.C:12
z
Double_t z
Definition: createRootFscGeometryFile.C:146
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
c1
c1
Definition: plot_dirc.C:35
PndRiemannTrack::Pt
double Pt(double B)
Definition: PndRiemannTrack.cxx:969
p2
TPad * p2
Definition: hist-t7.C:117
PndRiemannTrack::ChiSquareDistCircle
double ChiSquareDistCircle()
Definition: PndRiemannTrack.cxx:162
dx
double dx
Definition: runGeaneTestSimBox.cxx:44
phi
double phi
Definition: macro/detectors/emc/dedicated/reco_analys.C:24
PndRiemannTrack::r
double r() const
Definition: PndRiemannTrack.cxx:416
PndRiemannHit::z
double z() const
Definition: PndRiemannHit.cxx:87
PndRiemannTrack::calcIntersection
int calcIntersection(PndRiemannTrack &track, TVector3 &p1, TVector3 &p2)
Definition: PndRiemannTrack.cxx:610
PndRiemannTrack::covPlane
TMatrixD covPlane() const
Definition: PndRiemannTrack.h:95
PndRiemannTrack::mError
double mError() const
Definition: PndRiemannTrack.h:99
PndRiemannTrack::c
double c() const
Definition: PndRiemannTrack.h:57
PndRiemannTrack::orig
TVectorD orig() const
Definition: PndRiemannTrack.cxx:405
PndRiemannHit::x
const TVector3 & x() const
Definition: PndRiemannHit.h:71
ClassImp
ClassImp(PndAnaContFact)
p1
TPad * p1
Definition: hist-t7.C:116
PndRiemannTrack::fjacRXY
TMatrixD fjacRXY
jacobian matrix to transform from c,n1,n2,n2 to r,x,y
Definition: PndRiemannTrack.h:159
PndRiemannTrack::szFit
void szFit(bool withErrorCalc=true)
Definition: PndRiemannTrack.cxx:441
PndRiemannTrack::addHit
void addHit(PndRiemannHit &hit)
Definition: PndRiemannTrack.cxx:102
PndRiemannTrack::calcErrorS
Double_t calcErrorS(TVector2 &XY, TVector2 &dXY, PndRiemannTrack *track)
Definition: PndRiemannTrack.cxx:857
PndRiemannTrack::fc
double fc
distance of plane to origin
Definition: PndRiemannTrack.h:151
hit
PndSdsMCPoint * hit
Definition: anasim.C:70
MatrixOutput
void MatrixOutput(TMatrixD mat)
Definition: PndRiemannTrack.cxx:39
hits
CbmHit * hits[nHits]
Definition: RiemannTest.C:19
alpha
double alpha
Definition: f_Init.h:9
PndRiemannTrack::fcovPlane
TMatrixD fcovPlane
full covarince matrix of the plane;
Definition: PndRiemannTrack.h:158
PndRiemannHit::calcPosOnTrk
void calcPosOnTrk(PndRiemannTrack *trk)
Definition: PndRiemannHit.cxx:92
PndRiemannTrack::getPforHit
TVector3 getPforHit(int i, double B)
Definition: PndRiemannTrack.cxx:1003
dPhi
double dPhi
Definition: RiemannTest.C:17
PndRiemannTrack::fm
double fm
parameters of sz-fit
Definition: PndRiemannTrack.h:153
Pi
Double_t Pi
Definition: PndGeoHypGeCrystal.cxx:10
PndRiemannTrack::fav
TVectorD fav
average over all hits
Definition: PndRiemannTrack.h:150
PndRiemannHit::sigmaY
double sigmaY() const
Definition: PndRiemannHit.h:79
R
Double_t R
Definition: checkhelixhit.C:61
PndRiemannTrack::getPndTrack
PndTrack getPndTrack(Double_t B)
Definition: PndRiemannTrack.cxx:1105
PndSttTube::GetWireDirection
TVector3 GetWireDirection()
Definition: PndSttTube.cxx:107
PndRiemannTrack::calcAlpha
double calcAlpha(PndRiemannHit *myHit)
Definition: PndRiemannTrack.cxx:930
PndRiemannTrack::calcPosByS
TVector3 calcPosByS(double s)
Definition: PndRiemannTrack.cxx:589
PndRiemannTrack::getHit
PndRiemannHit * getHit(unsigned int i)
Definition: PndRiemannTrack.h:75
TMatrixD
TMatrixT< double > TMatrixD
Definition: PndLmdDim.h:52
PndRiemannTrack::ftError
double ftError
Error of fit.
Definition: PndRiemannTrack.h:156
PndRiemannTrack::szError
double szError(PndRiemannHit *hit)
Definition: PndRiemannTrack.cxx:910
PndRiemannTrack::Pl
double Pl(double B)
Definition: PndRiemannTrack.cxx:977
Generated on Wed Apr 3 2019 10:02:56 for FairRoot/PandaRoot by   doxygen 1.8.5