FairRoot/PandaRoot
Public Member Functions | Private Attributes | List of all members
PndTrkChi2Fits Class Reference

#include <PndTrkChi2Fits.h>

Inheritance diagram for PndTrkChi2Fits:

Public Member Functions

 PndTrkChi2Fits ()
 
 ~PndTrkChi2Fits ()
 
void Calculations_Mvd (bool *InclusionMvd, Double_t *Mvd_DipVar_DipVar, Double_t *Mvd_IndVar_DipVar, Double_t *Mvd_IndVar_IndVar, Short_t nMvdHits, Double_t &Mvd_DipVar_DipVar_Sum, Double_t &Mvd_IndVar_DipVar_Sum, Double_t &Mvd_IndVar_IndVar_Sum)
 
void Calculations_SkewStt_AllLeftRightCombinations (Short_t nSttHits, Double_t *Stt_DriftRad_DipVar, Double_t *Stt_DriftRad_IndVar, Double_t *Stt_DriftRad_DipVar_Sum, Double_t *Stt_DriftRad_IndVar_Sum)
 
Short_t FitHelixCylinder (Short_t nHitsinTrack, Double_t *Xconformal, Double_t *Yconformal, Double_t *DriftRadiusconformal, Double_t *ErrorDriftRadiusconformal, Double_t rotationangle, Double_t trajectory_vertex[2], Short_t NMAX, Double_t *m, Double_t *q, Double_t *pAlfa, Double_t *pBeta, Double_t *pGamma, bool *Type, int istampa, int IVOLTE)
 
Short_t FitSZspace (Short_t nHitsinTrack, Double_t *S, Double_t *Z, Double_t *DriftRadius, Double_t *ErrorDriftRadius, Double_t FInot, Short_t NMAX, Double_t *emme, int IVOLTE)
 
Short_t FitSZspace_Chi2_AnnealingtheMvdOnly (Short_t nHitsinTrack, Double_t *S, Double_t *Z, Double_t *DriftRadius, Double_t *ErrorDriftRadius, Double_t FInot, Short_t NMAX, Double_t *emme, int IVOLTE)
 
void GSumCalculation (Double_t *S, Double_t *Z1, Double_t *Z2, Double_t *Sigma, Double_t FInot, int nHits, Double_t *outSum)
 
void UinvSumCalculation (Double_t *S, Double_t *Sigma, Double_t FInot, int nHits, Double_t *outSum)
 
void ZqSumCalculation (Double_t *Z1, Double_t *Z2, Double_t *Sigma, Double_t FInot, int nHits, Double_t *outSum)
 
 ClassDef (PndTrkChi2Fits, 1)
 

Private Attributes

int fIcounter
 

Detailed Description

Definition at line 8 of file PndTrkChi2Fits.h.

Constructor & Destructor Documentation

PndTrkChi2Fits::PndTrkChi2Fits ( )

Default constructor

Definition at line 16 of file PndTrkChi2Fits.cxx.

16 {};
PndTrkChi2Fits::~PndTrkChi2Fits ( )
inline

Destructor

Definition at line 24 of file PndTrkChi2Fits.h.

24 {};

Member Function Documentation

void PndTrkChi2Fits::Calculations_Mvd ( bool *  InclusionMvd,
Double_t Mvd_DipVar_DipVar,
Double_t Mvd_IndVar_DipVar,
Double_t Mvd_IndVar_IndVar,
Short_t  nMvdHits,
Double_t Mvd_DipVar_DipVar_Sum,
Double_t Mvd_IndVar_DipVar_Sum,
Double_t Mvd_IndVar_IndVar_Sum 
)

Definition at line 20 of file PndTrkChi2Fits.cxx.

References i.

31 {
32  Short_t
33  i;
34  //j; //[R.K. 01/2017] unused variable?
35 
36  Mvd_DipVar_DipVar_Sum = 0. ;
37  Mvd_IndVar_DipVar_Sum = 0. ;
38  Mvd_IndVar_IndVar_Sum = 0. ;
39 
40 
41  for(i=0;i<nMvdHits; i++){
42  // if this Mvd hit has been annihiled, skip;
43  if(!InclusionMvd[i]) continue;
44 
45 
46  Mvd_DipVar_DipVar_Sum += Mvd_DipVar_DipVar[i] ;
47  Mvd_IndVar_DipVar_Sum += Mvd_IndVar_DipVar[i];
48  Mvd_IndVar_IndVar_Sum += Mvd_IndVar_IndVar[i];
49 
50  } // end of for(i=0;i<nMvdHits; i++)
51 
52 
53 }
Int_t i
Definition: run_full.C:25
void PndTrkChi2Fits::Calculations_SkewStt_AllLeftRightCombinations ( Short_t  nSttHits,
Double_t Stt_DriftRad_DipVar,
Double_t Stt_DriftRad_IndVar,
Double_t Stt_DriftRad_DipVar_Sum,
Double_t Stt_DriftRad_IndVar_Sum 
)

Definition at line 61 of file PndTrkChi2Fits.cxx.

References Double_t, and i.

69 {
70 
71  Short_t
72  current_number,
73  i,
74  j;
75 
76  Double_t
77  aux;
78 
79  current_number = 1;
80  Stt_DriftRad_IndVar_Sum[0] = 0.;
81  Stt_DriftRad_DipVar_Sum[0] = 0.;
82 
83  // calculate the quantities for the chi**2 minimization for all the other combinations;
84  // this part was written previously in an iterative way;
85  for(i=0; i<nSttHits; i++){
86  for(j=0; j<current_number; j++){ // current_number --> number of combinations formed so far;
87 
88  aux = Stt_DriftRad_IndVar_Sum[i];
89  Stt_DriftRad_IndVar_Sum[j] += Stt_DriftRad_IndVar[i] ;
90  Stt_DriftRad_IndVar_Sum[current_number + j] = aux - Stt_DriftRad_IndVar[i] ;
91 
92  aux = Stt_DriftRad_DipVar_Sum[i];
93  Stt_DriftRad_DipVar_Sum[j] += Stt_DriftRad_DipVar[i] ;
94  Stt_DriftRad_DipVar_Sum[current_number + j] = aux - Stt_DriftRad_DipVar[i] ;
95 
96  } // end of for(j=0; j<current_number; j++)
97 
98  current_number *= 2;
99 
100  } // end of for(i=0; i<nHitsinTrack; i++)
101 
102 
103 }
Int_t i
Definition: run_full.C:25
Double_t
PndTrkChi2Fits::ClassDef ( PndTrkChi2Fits  ,
 
)
Short_t PndTrkChi2Fits::FitHelixCylinder ( Short_t  nHitsinTrack,
Double_t Xconformal,
Double_t Yconformal,
Double_t DriftRadiusconformal,
Double_t ErrorDriftRadiusconformal,
Double_t  rotationangle,
Double_t  trajectory_vertex[2],
Short_t  NMAX,
Double_t m,
Double_t q,
Double_t pAlfa,
Double_t pBeta,
Double_t pGamma,
bool *  Type,
int  istampa,
int  IVOLTE 
)

Definition at line 111 of file PndTrkChi2Fits.cxx.

References cos(), Double_t, fabs(), i, mm, sin(), sqrt(), v1, and v2.

Referenced by PndTrkCTFindTrackInXY2::FindTrackInXYProjection(), and PndTrkTracking2::RefitMvdStt().

129 {
130 
131  int
132  i;
133 
134  Double_t
135  Alfa,
136  alfetta,
137  Beta,
138  cose = cos(rotationangle),
139  dete,
140  mm,
141  sigma2,
142  sine = sin(rotationangle),
143  Su,
144  Suu,
145  Suv,
146  Sv,
147  S1,
148  ui,
149  u1,
150  u2,
151  vi,
152  v1,
153  v2;
154 
155  // the parameters of the track trajectory in conformal space is a straight line;
156  // here it is parametrized as : (*pBeta) * V + (*pAlfa)*U +1 =0
157  // where U, V are the conformal coordinates : U = X/( X**2 + Y**2 );
158  // In this method *pAlfa, *pBeta are inputs from the fit to the candidate track
159  // performed earlier; also, assume *pGamma = 0.
160 
161  // perform the required translation of coordinates in new center trajectory_vertex[2]
162  // in XY space !!
163  // The arrays Xconformal and Yconformal are SUPPOSED TO BE CALCULATED ALREADY TAKING
164  // INTO ACCOUNT THE TRANSLATION, therefore no need to modify them here for that. However
165  // they still need to be rotated.
166 
167  // modify the input values for the translation :
168 
169  Alfa = (*pAlfa) + 2.*trajectory_vertex[0];
170  Beta = (*pBeta) + 2.*trajectory_vertex[1];
171 
172 
173  // now take into account the rotation by rotationangle;
174 
175  alfetta = Alfa;
176  Alfa = Alfa*cose + Beta*sine;
177  Beta = -alfetta*sine + Beta*cose;
178 
179 
180  Double_t
181  Xp[nHitsinTrack],
182  Yp[nHitsinTrack];
183 
184  // rotate the points;
185  for(i=0;i<nHitsinTrack; i++){
186  Xp[i] = Xconformal[ i ] *cose + Yconformal[ i ]*sine;
187  Yp[i] = -Xconformal[ i ] *sine + Yconformal[ i ]*cose;
188 
189 
190  }
191 //--------------------
192 
193  // starts fitting procedure by finding the POCA first and then doing the Chi**2 minimization;
194  // use the starting value of Alfa and Beta for finding the POCA
195 
196 
197  mm = sqrt( Alfa*Alfa + Beta*Beta );
198 
199 
200  Su = 0.;
201  Sv = 0.;
202  Suv = 0.;
203  Suu = 0.;
204  S1 = 0.;
205  for(i=0; i<nHitsinTrack; i++){
206 
207  if( DriftRadiusconformal[i]<0){ // this is a Mvd hit;
208  ui = Xp[i]; // U
209  vi = Yp[i]; // V
210  }else { // this is a Stt axial;
211  // find the POCA of the Stt hit to the original trajectory in conformal space;
212  u1 = Xp[i] - DriftRadiusconformal[i] * Alfa/mm;
213  v1 = Yp[i] - DriftRadiusconformal[i] * Beta/mm;
214  u2 = Xp[i] + DriftRadiusconformal[i] * Alfa/mm;
215  v2 = Yp[i] + DriftRadiusconformal[i] * Beta/mm;
216  // poca is the point closest to the straight line;
217  if(fabs(Beta*v1+Alfa*u1+1) < fabs(Beta*v2+Alfa*u2+1) ){
218  ui = u1;
219  vi = v1;
220  } else {
221  ui = u2;
222  vi = v2;
223  }
224  }
225 
226  // now the minimization with the chi2 formula;
227  sigma2 = ErrorDriftRadiusconformal[i]*ErrorDriftRadiusconformal[i];
228  Su += ui/sigma2;
229  Sv += vi/sigma2;
230 
231  Suv += ui * vi/sigma2;
232  Suu += ui* ui/sigma2;
233 
234  S1 += 1./sigma2;
235 
236  } // end of for(i=0; i<nHitsinTrack; i++)
237 
238 
239  dete = Suu * S1 - Su * Su;
240  if(fabs(dete) < 1e-10) { *Type = false; return -5;} // fit fails;
241 
242  *emme = (Suv * S1 - Su * Sv)/dete;
243  *qu = (Suu * Sv - Su * Suv)/dete;
244  // protect the extreme case of q=0 (in principle not possible, it would
245  // correspond to a straight line trajectory y = emme*x in the XY plane;
246  if( fabs(*qu) < 1.e-9 ) {
247 // if ((*qu) > 0. ) *qu = 1.e-9 ; else *qu = -1.e-9;
248  *Type = false;
249  return -4 ; // fit fails;
250  } else {
251  *Type = true ; // normal case;
252  }
253 
254  // calculate the output coefficients of the circumference in XY plane;
255 
256  Alfa = (*emme)/(*qu);
257  Beta = -1./(*qu);
258 
259 // now take into account the rotation and correct back; the affected quantities are ALFA and BETA,
260 // (*emme and *qu also but later);
261  alfetta = Alfa;
262  Alfa = Alfa*cose - Beta*sine;
263  Beta = alfetta*sine + Beta*cose;
264 
265 // calculate the coefficients taking into account the translation;
266 // now take into account the displacement and calculate Gamma;
267 
268  *pGamma = trajectory_vertex[0]*trajectory_vertex[0]+ // *pGamma is assumed to be 0 at the beginning;
269  trajectory_vertex[1]*trajectory_vertex[1]
270  -Alfa*trajectory_vertex[0]-Beta*trajectory_vertex[1];
271 
272  // calculate the coefficients taking into account the translation;
273  Alfa -= 2.*trajectory_vertex[0];
274  Beta -= 2.*trajectory_vertex[1];
275 
276 // calculate *emme and *qu using the newly calculated *pAlfa, *pBeta, *pGamma of the
277 // circular trajectory in XY. Assuming also that *pGamma ~ 0, that is, the circumference
278 // goes thru the origin.
279 
280  if( abs(Beta)>1e-9) {
281  // normal case of a straight line in UV that can be put in the V = m*U +q form;
282  *emme = -Alfa/Beta;
283  *qu = -1./Beta;
284  *pAlfa = Alfa;
285  *pBeta = Beta;
286  return 1;
287  } else {
288  *emme = -Alfa/1e-9;
289  *qu = -1./1e-9;
290  *pAlfa = Alfa;
291  *pBeta = Beta;
292  return 99;
293  }
294 
295 }
friend F32vec4 cos(const F32vec4 &a)
Definition: P4_F32vec4.h:112
Int_t i
Definition: run_full.C:25
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
friend F32vec4 sin(const F32vec4 &a)
Definition: P4_F32vec4.h:111
Double_t
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
Double_t const mm
TVector3 v1
Definition: bump_analys.C:40
TVector3 v2
Definition: bump_analys.C:40
Short_t PndTrkChi2Fits::FitSZspace ( Short_t  nHitsinTrack,
Double_t S,
Double_t Z,
Double_t DriftRadius,
Double_t ErrorDriftRadius,
Double_t  FInot,
Short_t  NMAX,
Double_t emme,
int  IVOLTE 
)

Definition at line 301 of file PndTrkChi2Fits.cxx.

References Double_t, exit(), fabs(), and i.

312 {
313 
314  if(nHitsinTrack<0)
315  {
316  cout<<"from PndTrkChi2Fits::FitSZspace error, n hits to fit = 0, exit(-1).";
317  exit(-1);
318  }
319 
320 
321  int
322  Combinations,
323  i,
324  icomb,
325  nSttHits;
326 
327  Double_t
328  e2,
329  mvdGSum,
330  mvdUinvSum,
331  mvdZqSum,
332  sttS[nHitsinTrack], // S of only the STT hits;
333  // this is a little overdimensioned, but easy to do;
334  sttSigma[nHitsinTrack], // as above; first possibility of Z;
335  sttZ1[nHitsinTrack], // as above; first possibility of Z;
336  sttZ2[nHitsinTrack] // as above; second possibility of Z;
337  ;
338 
339 
340 
341  nSttHits = 0 ;
342  mvdGSum=0.;
343  mvdUinvSum=0.;
344  mvdZqSum=0.;
345 
346  for(i=0; i<nHitsinTrack; i++){
347  if( DriftRadius[i]<0. )
348  {
349  // Mvd hit;
350  // relevant quantities for the chi2 calculation and minimization;
351  e2 = ErrorDriftRadius[i]*ErrorDriftRadius[i];
352  mvdGSum += Z[i]*(S[i]-FInot)/e2;
353  mvdUinvSum += (S[i]-FInot)*(S[i]-FInot)/e2;
354  mvdZqSum += Z[i]*Z[i]/e2 ;
355 
356 
357  } else {
358  // Stt hit;
359  sttS[nSttHits] = S[i];
360  sttSigma[nSttHits] = ErrorDriftRadius[i];
361  sttZ1[nSttHits] = Z[i] + DriftRadius[i];
362  sttZ2[nSttHits] = Z[i] - DriftRadius[i];
363  nSttHits++;
364  }
365  } // end of for(i=0; i<nHitsinTrack; i++)
366 
367 // limit the number of Skew hits in fit, otherwise the number of combinations
368 // becomes too large;
369  if(nSttHits>NMAX) nSttHits=NMAX;
370 
371 
372 //---- begin the fit procedure. Here S is the INDEPENDENT-like variable and Z the
373 // dependent-like one.
374 // For each Stt (Skew) hit the two possibilities are taken into account; therefore the
375 // number f chi**2 are 2**nSttHits ; at the end of the procedure the Stt (Skew) hit pattern
376 // is chosen with the least chi**2;
377 // use a recursive algorithm to calculate all the chi**2 and minimize them because the
378 // a priori is NOT know how may Stt (Skew) hits are there in the track candidate.
379 
380 
381  Combinations = (int) (pow(2.,nSttHits) + 0.1) ;// +0.1 only for beeing absolutely sure agains rounding errors;
382 
383  Double_t
384  GSum[Combinations],
385  UinvSum[Combinations],
386  ZqSum[Combinations];
387 
388 
389 
390  Double_t
391  chi2,
392  Chi2min,
393  Kappa,
394  TotalGSum,
395  TotalUinvSum,
396  TotalZqSum;
397 
398 
399  if(nSttHits>0){
400 // the following is a recursive function;
401 
403  sttS, // input; the independent-like variable;
404  sttZ1, // input; the dependent-like variable;
405  sttZ2, // input; the dependent-like variable;
406  sttSigma, // input; the errors on sttZ;
407  FInot,
408  nSttHits,
409  GSum // the output; this must be an array of 2**nSttHits elements;
410  );
411 
412 
413 
414 // the following is a recursive function;
415 
417  sttS, // input; the independent-like variable;
418  sttSigma, // input; the errors on sttZ;
419  FInot,
420  nSttHits,
421  UinvSum // the output; this must be an array of 2**nSttHits elements;
422  );
423 
424 // the following is a recursive function;
425 
427  sttZ1, // input; the dependent-like variable;
428  sttZ2, // input; the dependent-like variable;
429  sttSigma, // input; the errors on sttZ;
430  FInot,
431  nSttHits,
432  ZqSum // the output; this must be an array of 2**nSttHits elements;
433  );
434 
435 // calculation of the Kappa parameter (the parametrization is : Z = (S-FInot)/Kappa );
436 // and the all chi**2;
437 
438  // the first combination;
439 
440  TotalGSum = mvdGSum + GSum[0];
441  TotalUinvSum = mvdUinvSum + UinvSum[0];
442  TotalZqSum = mvdZqSum + ZqSum[0];
443  // calculation of the fit parameter obtained by minimization of the chi**2;
444  if(fabs(TotalGSum)>1e-9)
445  {
446  Kappa = TotalUinvSum / TotalGSum;
447  }else {
448  Kappa = 1e9;
449  }
450  // calculation of the corresponding chi**2;
451  Chi2min = TotalZqSum + TotalUinvSum/Kappa -2.*TotalGSum/Kappa;
452  // for the moment this is the solution;
453  *emme = Kappa;
454 
455  // all the other combinations;
456  for(icomb=1; icomb<Combinations; icomb++){
457  TotalGSum = mvdGSum + GSum[icomb];
458  TotalUinvSum = mvdUinvSum + UinvSum[icomb];
459  TotalZqSum = mvdZqSum + ZqSum[icomb];
460  // calculation of the fit parameter obtained by minimization of the chi**2;
461  if(fabs(TotalGSum) > 1.e-8){
462  Kappa = TotalUinvSum / TotalGSum;
463  // calculation of the corresponding chi**2;
464  chi2 = TotalZqSum + TotalUinvSum/Kappa -2.*TotalGSum/Kappa;
465  // choose the solution with best chi**2;
466  if( chi2 <Chi2min) *emme = Kappa;
467  } else {
468  if(TotalGSum>0.) Kappa = 1.e8 ; else Kappa = -1.e8 ;
469  } // end of if(fabs(TotalGSum) > 1.e-8)
470 
471  } // end of for(icomb=1; icomb<Combinations; i++)
472 
473 
474  } else { // continuation of if(nSttHits>0)
475  // here nSttHits is 0;
476  TotalGSum = mvdGSum;
477  TotalUinvSum = mvdUinvSum;
478  TotalZqSum = mvdZqSum ;
479 
480  if(fabs(TotalGSum) > 1.e-8){
481  Kappa = TotalUinvSum / TotalGSum;
482  // choose the solution with best chi**2;
483  *emme = Kappa;
484  } else {
485  if(TotalGSum>0.) Kappa = 1.e8 ; else Kappa = -1.e8 ;
486  } // end of if(fabs(TotalGSum) > 1.e-8)
487 
488  } // end of if(nSttHits>0)
489 
490 
491 
492 return 1;
493 }
Int_t i
Definition: run_full.C:25
exit(0)
void ZqSumCalculation(Double_t *Z1, Double_t *Z2, Double_t *Sigma, Double_t FInot, int nHits, Double_t *outSum)
Double_t
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
void GSumCalculation(Double_t *S, Double_t *Z1, Double_t *Z2, Double_t *Sigma, Double_t FInot, int nHits, Double_t *outSum)
void UinvSumCalculation(Double_t *S, Double_t *Sigma, Double_t FInot, int nHits, Double_t *outSum)
double Z
Definition: anaLmdDigi.C:68
Short_t PndTrkChi2Fits::FitSZspace_Chi2_AnnealingtheMvdOnly ( Short_t  nHitsinTrack,
Double_t S,
Double_t Z,
Double_t DriftRadius,
Double_t ErrorDriftRadius,
Double_t  FInot,
Short_t  NMAX,
Double_t emme,
int  IVOLTE 
)

Definition at line 503 of file PndTrkChi2Fits.cxx.

References Double_t, fabs(), i, and status.

Referenced by PndTrkTracking2::Exec().

514 {
515  bool
516  InclusionMvd[nHitsinTrack];
517 
518  Short_t
519  i,
520  j,
521  //nHitsinFit, //[R.K. 01/2017] unused variable?
522  nMvdHits,
523  nSttHits,
524  status;
525 
526  Int_t
527  Combinations;
528 
529  Double_t
530  A,
531  chi2,
532  chi2_fixed,
533  chi2_best,
534  e2,
535  M,
536 
537  Mvd_DipVar_DipVar[nHitsinTrack],
538  Mvd_DipVar_DipVar_Sum,
539  Mvd_IndVar_DipVar[nHitsinTrack],
540  Mvd_IndVar_DipVar_Sum,
541  Mvd_IndVar_IndVar[nHitsinTrack],
542  Mvd_IndVar_IndVar_Sum,
543  //Mvd_invError2[nHitsinTrack], //[R.K.02/2017] Unused variable?
544 
545  Penalty,
546 
547  //Stt_DipVar_Sum, //[R.K. 01/2017] unused variable?
548  Stt_DipVar_DipVar_Sum,
549  Stt_DriftRad_DipVar[nHitsinTrack],
550  Stt_DriftRad_DriftRad_Sum,
551  Stt_DriftRad_IndVar[nHitsinTrack],
552  //Stt_IndVar_Sum, //[R.K. 01/2017] unused variable?
553  Stt_IndVar_DipVar_Sum,
554  Stt_IndVar_IndVar_Sum;
555 
556 
557 
558 
559 
560  // nHitsinTrack = n. hits in this track candidate;
561  // *S = array of the independent variable of the fit, namely the R*fi variable on
562  // the side surface of the Helix cylinder;
563  // *Z = array of the dependent variable of the fit; namely the Z position of the hit
564  // projected on the side surface of the Helix cylinder; (Z coordinate for a
565  // Mvd hit, intersection between wire and Helix cylinder for Skew Stt);
566  // *DriftRadius = array of Radius of drift for the Skew Sraw PROJECTED onto the lateral surface
567  // of the Helix, namely the major axis of the ellipsis projection
568  // of the Skew Straw; for the Mvd hits it is < 0;
569  // *ErrorDriftRadius = array of Errors associated to DriftRadius; presently (8 Aug 2013) they
570  // are set to 0.5 cm for the Mvd hits, and to 2 times DriftRadius for the
571  // Skew Straws;
572  // FInot = this is a fixed input from the previous fits in XY; it is NOT a output variable
573  // in this fit;
574  // NMAX = maximum number of Skew Straws allowed to partecipate in this fit;
575  // *emme = output of the fit;
576  // IVOLTE = service variable indicating the current event;
577 
578 
579  // This method calculates the Chi2 for :
580  // 1) all Mvd hits plus all left/right combinations of the Skew Stt;
581  // 2) like the above but excluding one Mvd hit and replacing it with a penalty term;
582  // at the end the minimum Chi2 among all these is chosen and the fit parameter accordingly.
583 
584 
585 
586  // Step 1 : calculation of all the Mvd hits plus all the combinations left/right of the Skew Stt;
587  // here Mvd hits means Mvd (these have DriftRadius = -1) or SciTil (these have DriftRadius = -2);
588  nMvdHits=0;
589  nSttHits=0;
590  Stt_DipVar_DipVar_Sum = 0. ;
591  Stt_DriftRad_DriftRad_Sum = 0. ;
592  Stt_IndVar_IndVar_Sum = 0. ;
593  Stt_IndVar_DipVar_Sum = 0. ;
594 
595  for(i=0; i<nHitsinTrack; i++){
596 
597  e2 = 1./(ErrorDriftRadius[i]*ErrorDriftRadius[i]);
598 
599  if( DriftRadius[i]<0. )
600  {
601  // Mvd hit;
602  // relevant quantities for the chi2 calculation and minimization;
603  // Here S is the independent variable (== Mvd_IndVar), Z is the
604  // dependent variable (== Mvd_DepVar);
605  // the fit function is : Z = (1/Kappa)*(S-FInot);
606  Mvd_DipVar_DipVar[nMvdHits] = Z[i]*Z[i]*e2;
607  Mvd_IndVar_IndVar[nMvdHits] =(S[i]-FInot)*(S[i]-FInot)*e2;
608  Mvd_IndVar_DipVar[nMvdHits] =(S[i]-FInot)*Z[i]*e2;
609  // the following is useful only for adding the penalty term later;
610  //Mvd_invError2[nMvdHits] = e2; //[R.K.02/2017] Unused variable?
611  nMvdHits++;
612  } else {
613  // Skew Straw hit;
614  // relevant quantities for the chi2 calculation and minimization;
615  // Here S is the independent variable (== Stt_IndVar), Z is the
616  // dependent variable (== Stt_DepVar);
617  // the fit function is : Z = (1/Kappa)*(S-FInot);
618 
619  if( nSttHits < NMAX){
620 
621 
622  Stt_DipVar_DipVar_Sum += Z[i]*Z[i]*e2;
623  Stt_DriftRad_DriftRad_Sum += DriftRadius[i]*DriftRadius[i]*e2;
624  Stt_DriftRad_IndVar[nSttHits] = DriftRadius[i]*(S[i]-FInot)*e2;
625  Stt_DriftRad_DipVar[nSttHits] = DriftRadius[i]*Z[i]*e2;
626 
627  Stt_IndVar_IndVar_Sum += (S[i]-FInot)*(S[i]-FInot)*e2;
628  Stt_IndVar_DipVar_Sum += (S[i]-FInot)*Z[i]*e2;
629  }
630  nSttHits++;
631  }
632  } // end of for(i=0; i<nHitsinTrack; i++)
633 
634 
635  // limit the number of Skew hits in fit, otherwise the number of combinations
636  // becomes too large;
637  if(nSttHits>NMAX) nSttHits=NMAX;
638 
639  //---- begin the fit procedure. Here S is the INDEPENDENT-like variable and Z the
640  // dependent-like one.
641 
642 
643  status = 0; // failure in fitting (just initializing!);
644 
645  // For the Skew Stt hit calculate the contribution of all possible left/right combinations once and for all
646  // since no 'annealing' mechanism is used for them in this algorithm;
647  // only the terms containing DriftRadius (linearly) change depending on the combinations;
648 
649  // Combinations == number of all possible left/right combinations given the Skew Stt in the track;
650  Combinations = (Int_t) (pow(2., (double) nSttHits) + 0.1) ;// +0.1 only for being absolutely
651  // sure agains rounding errors;
652 
653  Double_t
654  Stt_DriftRad_DipVar_Sum[Combinations],
655  Stt_DriftRad_IndVar_Sum[Combinations];
656 
657 
659  nSttHits, // input;
660  Stt_DriftRad_DipVar, // input;
661  Stt_DriftRad_IndVar, // input;
662 
663  Stt_DriftRad_DipVar_Sum, // output;
664  Stt_DriftRad_IndVar_Sum // output;
665  );
666 
667  // For the Mvd hits, calculate the contribution of all Mvd, and all Mvd except 1 (a penalty term instead)
668 
669  // all Mvd contribution;
670  memset(InclusionMvd, true, sizeof(InclusionMvd));
672  InclusionMvd, // input;
673  Mvd_DipVar_DipVar, // input;
674  Mvd_IndVar_DipVar, // input;
675  Mvd_IndVar_IndVar, // input;
676  nMvdHits, // input;
677 
678  Mvd_DipVar_DipVar_Sum, // output;
679  Mvd_IndVar_DipVar_Sum, // output;
680  Mvd_IndVar_IndVar_Sum // output;
681  );
682 
683  // calculation with all the Mvd hits considered;
684 
685  // formula used here : emme * (Sum_i (x_i **2 / sigma_i**2)) = Sum_i ( x_i*y_i/sigma_i**2) - qu *(Sum_i x_i/sigma_i**2);
686  chi2_best = 9999999.;
687  A = Mvd_IndVar_DipVar_Sum + Stt_IndVar_DipVar_Sum ;
688  chi2_fixed =
689  Mvd_DipVar_DipVar_Sum + Stt_DipVar_DipVar_Sum
690  + Stt_DriftRad_DriftRad_Sum;
691 
692  for(i=0;i<Combinations; i++){
693  M = (A + Stt_DriftRad_IndVar_Sum[i])/(Stt_IndVar_IndVar_Sum + Mvd_IndVar_IndVar_Sum) ;
694  chi2 = chi2_fixed +
695  M*M*(Mvd_IndVar_IndVar_Sum + Stt_IndVar_IndVar_Sum)
696  + 2.*Stt_DriftRad_DipVar_Sum[i]
697  - 2.*M*(Mvd_IndVar_DipVar_Sum + Stt_IndVar_DipVar_Sum)
698  - 2.*Stt_DriftRad_IndVar_Sum[i]*M;
699  if(chi2<chi2_best){
700  chi2_best = chi2;
701  *emme = M;
702  status = 1; // success in fitting;
703  }
704  } // end of for(i=0;i<nCombinations; i++)
705 
706 
707  // calculate the chi2/degrees of freedom; in case there is only 1 hit than just keep the chi2_best as it is;
708  if(nSttHits+nMvdHits > 1) { chi2_best = chi2_best/( nSttHits+nMvdHits-1); }
709 
710 
711 
712 
713  if( nSttHits+nMvdHits >1){
714 
715  // exclude 1 Mvd hit, add penalty term to chi2; DON'T DO IT in the situation when
716  // there is only one hit associated to this track (because when such a hit is a Mvd hit,
717  // the program crashes; in fact the track remains without hits and the chi**2 cannot be
718  // calculated);
719 
720 // // Penalty = 9.*Mvd_invError2[j];
721 // Penalty = 36.;
722  Penalty = 0.;
723  for( j =0; j<nMvdHits;j++){
724  InclusionMvd[j] = false;
725  // the following is necessary to include again the Mvd hit excluded in the previous loop;
726  if(j>0) InclusionMvd[j-1] = true;
727 
728  // calculate the new Mvd contributions with 1 hit excluded;
730  InclusionMvd, // input;
731  Mvd_DipVar_DipVar, // input;
732  Mvd_IndVar_DipVar, // input;
733  Mvd_IndVar_IndVar, // input;
734  nMvdHits, // input;
735 
736  Mvd_DipVar_DipVar_Sum, // output;
737  Mvd_IndVar_DipVar_Sum, // output;
738  Mvd_IndVar_IndVar_Sum // output;
739  );
740 
741  // redo the chi2 minimization and best chi2 choice;
742  A = Mvd_IndVar_DipVar_Sum + Stt_IndVar_DipVar_Sum ;
743  chi2_fixed = Penalty +
744  Mvd_DipVar_DipVar_Sum + Stt_DipVar_DipVar_Sum
745  + Stt_DriftRad_DriftRad_Sum;
746  for(i=0;i<Combinations; i++){
747  M = (A + Stt_DriftRad_IndVar_Sum[i])/(Stt_IndVar_IndVar_Sum + Mvd_IndVar_IndVar_Sum) ;
748  chi2 = chi2_fixed +
749  M*M*(Mvd_IndVar_IndVar_Sum + Stt_IndVar_IndVar_Sum)
750  + 2.*Stt_DriftRad_DipVar_Sum[i]
751  - 2.*M*(Mvd_IndVar_DipVar_Sum + Stt_IndVar_DipVar_Sum)
752  - 2.*Stt_DriftRad_IndVar_Sum[i]*M;
753 
754  // comparison now has to be done between chi2/degrees of freedom; in the case
755  // when nSttHits+nMvdHits = 2, just leave the chi2 as it is;
756 
757  if( nSttHits+nMvdHits-2 >0) chi2 = chi2/(nSttHits+nMvdHits-2);
758 
759  if(chi2<chi2_best){
760  chi2_best = chi2;
761  *emme = M;
762  status = 1; // success in fitting;
763  }
764  } // end of for(i=0;i<nCombinations; i++)
765 
766  } // end of for( j =0; j<nMvdHits;j++)
767 
768 
769 
770  } // end of if(nSttHits+nMvdHits > 1)
771 
772 
773  // at this moment *emme corresponds to 1/KAPPA so now it is inverted;
774  if( status>0){
775  if( fabs(*emme) > 1.e-10 ) { *emme = 1./(*emme); } else { *emme = 1.e10; }
776  }
777 
778 
779 
780 
781 
782  return status;
783 
784 }
Int_t i
Definition: run_full.C:25
Double_t
void Calculations_SkewStt_AllLeftRightCombinations(Short_t nSttHits, Double_t *Stt_DriftRad_DipVar, Double_t *Stt_DriftRad_IndVar, Double_t *Stt_DriftRad_DipVar_Sum, Double_t *Stt_DriftRad_IndVar_Sum)
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
void Calculations_Mvd(bool *InclusionMvd, Double_t *Mvd_DipVar_DipVar, Double_t *Mvd_IndVar_DipVar, Double_t *Mvd_IndVar_IndVar, Short_t nMvdHits, Double_t &Mvd_DipVar_DipVar_Sum, Double_t &Mvd_IndVar_DipVar_Sum, Double_t &Mvd_IndVar_IndVar_Sum)
double Z
Definition: anaLmdDigi.C:68
int status[10]
Definition: f_Init.h:28
void PndTrkChi2Fits::GSumCalculation ( Double_t S,
Double_t Z1,
Double_t Z2,
Double_t Sigma,
Double_t  FInot,
int  nHits,
Double_t outSum 
)

Definition at line 796 of file PndTrkChi2Fits.cxx.

References Double_t, and i.

805 {
806 
807  int i,
808  Combinations,
809  Combinations2;
810 
811  double e2;
812 
813 
814  // this method calculates the sum calleg g in the PDG minimization formula;
815  // it is a recursive method;
816 
817 
818  Combinations = (int) (pow(2., nHits) + 0.1) ; // +0.1 to be absolutely sure agains rounding errors;
819  Combinations2 = (int) (pow(2., nHits-1) + 0.1);
820 
821  if(nHits==1)
822  {
823  // two possible combinations;
824  outSum[0] = Z1[0]*(S[0]-FInot)/(Sigma[0]*Sigma[0]);
825  outSum[1] = Z2[0]*(S[0]-FInot)/(Sigma[0]*Sigma[0]);
826  return;
827  }
828 
829  Double_t aux[Combinations];
830 
831  // the following is the recursion;
832 
834  S, // input; the independent-like variable;
835  Z1, // input; the dependent-like variable;
836  Z2, // input; the dependent-like variable;
837  Sigma, // input; the errors on sttZ;
838  FInot,
839  nHits-1,
840  outSum // the output; this must be an array of 2**nSttHits elements;
841  );
842 
843  // for each of the Combinations2 = 2**(nHits-1) already calculated combinations
844  // of Stt Skew hits, calculate two combinations by using the last SttSkew hit
845  // (its Z1, Z2, S are Z1[nHits-1], Z2[nHits-1] and so on);
846  // in this way a total of Combinations = 2**nHits combinations exist now.
847 
848  for(i=0; i<Combinations2; i++){
849  e2 = Sigma[nHits-1]*Sigma[nHits-1];
850  aux[i] = outSum[i] + Z1[nHits-1]*(S[nHits-1]-FInot)/e2;
851  aux[i+Combinations2] = outSum[i] + Z2[nHits-1]*(S[nHits-1]-FInot)/e2;
852  } // end of for(i=0; i<Combinations; i++)
853 
854  // reload the output array (that has now more elements);
855  for(i=0; i<Combinations; i++){
856  outSum[i] = aux[i];
857  }
858 
859 
860  return;
861 }
Int_t i
Definition: run_full.C:25
int nHits
Definition: RiemannTest.C:16
Double_t
void GSumCalculation(Double_t *S, Double_t *Z1, Double_t *Z2, Double_t *Sigma, Double_t FInot, int nHits, Double_t *outSum)
void PndTrkChi2Fits::UinvSumCalculation ( Double_t S,
Double_t Sigma,
Double_t  FInot,
int  nHits,
Double_t outSum 
)

Definition at line 871 of file PndTrkChi2Fits.cxx.

References Double_t, and i.

878 {
879 
880  int i,
881  Combinations,
882  Combinations2;
883 
884  double e2;
885 
886  // this method calculates the sum calleg g in the PDG minimization formula;
887  // it is a recursive method;
888 
889 
890  Combinations = (int) (pow(2., nHits)+0.1); // +0.1 to be absolutely sure agains rounding errors;
891  Combinations2 = (int) (pow(2., nHits-1)+0.1);
892 
893  if(nHits==1)
894  {
895  // two possible combinations; they have the same Uinv;
896  outSum[0] = (S[0]-FInot)*(S[0]-FInot)/(Sigma[0]*Sigma[0]);
897  outSum[1] = outSum[0];
898  return;
899  }
900 
901  Double_t aux[Combinations];
902 
903  // the following is the recursion;
904 
906  S, // input; the independent-like variable;
907  Sigma, // input; the errors on sttZ;
908  FInot,
909  nHits-1,
910  outSum // the output; this must be an array of 2**nSttHits elements;
911  );
912 
913  // for each of the Combinations2 = 2**(nHits-1) already calculated combinations
914  // of Stt Skew hits, calculate two combinations by using the last SttSkew hit
915  // (its Z1, Z2, S are Z1[nHits-1], Z2[nHits-1] and so on);
916  // in this way a total of Combinations = 2**nHits combinations exist now.
917 
918  for(i=0; i<Combinations2; i++){
919  e2 = Sigma[nHits-1]*Sigma[nHits-1];
920  // two combinations but with the same U inv;
921  aux[i] = outSum[i] + (S[nHits-1]-FInot)*(S[nHits-1]-FInot)/e2;
922  aux[i+Combinations2] = aux[i] ;
923  } // end of for(i=0; i<Combinations; i++)
924 
925  // reload the output array (that has now more elements);
926  for(i=0; i<Combinations; i++){
927  outSum[i] = aux[i];
928  }
929 
930  return;
931 }
Int_t i
Definition: run_full.C:25
int nHits
Definition: RiemannTest.C:16
Double_t
void UinvSumCalculation(Double_t *S, Double_t *Sigma, Double_t FInot, int nHits, Double_t *outSum)
void PndTrkChi2Fits::ZqSumCalculation ( Double_t Z1,
Double_t Z2,
Double_t Sigma,
Double_t  FInot,
int  nHits,
Double_t outSum 
)

Definition at line 940 of file PndTrkChi2Fits.cxx.

References Double_t, and i.

948 {
949 
950  int i,
951  Combinations,
952  Combinations2;
953 
954  double e2;
955 
956  // this method calculates the sum calleg g in the PDG minimization formula;
957  // it is a recursive method;
958 
959 
960  Combinations = pow(2., nHits)+0.1; // +0.1 to be absolutely sure against rounding errors;
961  Combinations2 = pow(2., nHits-1)+0.1;
962 
963  if(nHits==1)
964  {
965  // two possible combinations;
966  outSum[0] = Z1[0]*Z1[0]/(Sigma[0]*Sigma[0]);
967  outSum[1] = Z2[0]*Z2[0]/(Sigma[0]*Sigma[0]);
968  return;
969  }
970 
971  Double_t aux[Combinations];
972 
973  // the following is the recursion;
974 
976  Z1, // input; the dependent-like variable;
977  Z2, // input; the dependent-like variable;
978  Sigma, // input; the errors on sttZ;
979  FInot,
980  nHits-1,
981  outSum // the output; this must be an array of 2**nSttHits elements;
982  );
983 
984  // for each of the Combinations2 = 2**(nHits-1) already calculated combinations
985  // of Stt Skew hits, calculate two combinations by using the last SttSkew hit
986  // (its Z1, Z2, S are Z1[nHits-1], Z2[nHits-1] and so on);
987  // in this way a total of Combinations = 2**nHits combinations exist now.
988 
989  for(i=0; i<Combinations2; i++){
990  e2 = Sigma[nHits-1]*Sigma[nHits-1];
991  aux[i] = outSum[i] + Z1[nHits-1]*Z1[nHits-1]/e2;
992  aux[i+Combinations2] = outSum[i] + Z2[nHits-1]*Z2[nHits-1]/e2;
993  } // end of for(i=0; i<Combinations; i++)
994 
995  // reload the output array (that has now more elements);
996  for(i=0; i<Combinations; i++){
997  outSum[i] = aux[i];
998  }
999 
1000 
1001  return;
1002 }
Int_t i
Definition: run_full.C:25
int nHits
Definition: RiemannTest.C:16
void ZqSumCalculation(Double_t *Z1, Double_t *Z2, Double_t *Sigma, Double_t FInot, int nHits, Double_t *outSum)
Double_t

Member Data Documentation

int PndTrkChi2Fits::fIcounter
private

Definition at line 13 of file PndTrkChi2Fits.h.


The documentation for this class was generated from the following files: