HLT/src/AliL3ConfMapFit.cxx

   1 #include <math.h>
   2
   3 #include "AliL3Logging.h"
   4 #include "AliL3ConfMapFit.h"
   5 #include "AliL3Vertex.h"
   6 #include "AliL3ConfMapTrack.h"
   7 #include "AliL3ConfMapPoint.h"
   8
   9 //_______________________________
  10 // AliL3ConfMapFit
  11 //
  12 // Fit class for conformal mapping tracking
  13
  14 ClassImp(AliL3ConfMapFit)
  15
  16 Double_t AliL3ConfMapFit::pi=3.14159265358979323846;
  17
  18 AliL3ConfMapFit::AliL3ConfMapFit(AliL3ConfMapTrack *track,AliL3Vertex *vertex)
  19 {
  20   //constructor
  21   fTrack = track;
  22   fVertex = vertex;
  23   BFACT = 0.0029980;
  24   bField = 0.2;
  25
  26 }
  27
  28 Int_t AliL3ConfMapFit::FitHelix()
  29 {
  30   if(FitCircle())
  31     {
  32       LOG(AliL3Log::kError,"AliL3ConfMapFit::FitHelix","TrackFit")<<AliL3Log::kDec<<
  33         "Problems during circle fit"<<ENDLOG;
  34       return 1;
  35     }
  36   if(FitLine())
  37     {
  38       LOG(AliL3Log::kError,"AliL3ConfMapFit::FitHelix","TrackFit")<<AliL3Log::kDec<<
  39         "Problems during line fit"<<ENDLOG;
  40       return 1;
  41     }
  42   return 0;
  43 }
  44
  45 Int_t AliL3ConfMapFit::FitCircle()
  46 {
  47   //-----------------------------------------------------------------
  48   //Fits circle parameters using algorithm
  49   //described by ChErnov and Oskov in Computer Physics
  50   //Communications.
  51   //
  52   //Written in FORTRAN by Jawluen Tang, Physics department , UT-Austin
  53   //Moved to C by Pablo Yepes
  54   //Moved to AliROOT by ASV.
  55   //------------------------------------------------------------------
  56
  57   Double_t wsum  = 0.0 ;
  58   Double_t xav   = 0.0 ;
  59   Double_t yav   = 0.0 ;
  60
  61   Int_t num_of_hits = fTrack->GetNumberOfPoints();
  62   //
  63   //     Loop over hits calculating average
  64   Int_t co=0;
  65
  66   for(fTrack->StartLoop(); fTrack->LoopDone(); fTrack->GetNextHit())
  67     {
  68       co++;
  69       AliL3ConfMapPoint *cHit = (AliL3ConfMapPoint*)fTrack->currentHit;
  70       cHit->SetXYWeight( 1./ (Double_t)(cHit->GetXerr()*cHit->GetXerr() + cHit->GetYerr()*cHit->GetYerr()) );
  71       wsum      += cHit->GetXYWeight() ;
  72       xav       += cHit->GetXYWeight() * cHit->GetX() ;
  73       yav       += cHit->GetXYWeight() * cHit->GetY() ;
  74     }
  75   if(co!=num_of_hits)
  76     LOG(AliL3Log::kError,"AliL3ConfMapFit::FitCircle","TrackFit")<<AliL3Log::kDec<<
  77       "Mismatch of hits. Counter: "<<co<<" nHits: "<<num_of_hits<<ENDLOG;
  78   if (fTrack->ComesFromMainVertex() == true)
  79     {
  80       wsum += fVertex->GetXYWeight() ;
  81       xav  += fVertex->GetX() ;
  82       yav  += fVertex->GetY() ;
  83     }
  84
  85   xav = xav / wsum ;
  86   yav = yav / wsum ;
  87 //
  88 //  CALCULATE <X**2>, <XY>, AND <Y**2> WITH <X> = 0, & <Y> = 0
  89 //
  90   Double_t xxav  = 0.0 ;
  91   Double_t xyav  = 0.0 ;
  92   Double_t yyav  = 0.0 ;
  93   Double_t xi, yi ;
  94
  95   for(fTrack->StartLoop(); fTrack->LoopDone(); fTrack->GetNextHit())
  96     {
  97       //AliL3ConfMapPoint *cHit = (AliL3ConfMapPoint *)hits->At(hit_counter);
  98       AliL3ConfMapPoint *cHit = (AliL3ConfMapPoint*)fTrack->currentHit;
  99       xi        = cHit->GetX() - xav ;
 100       yi        = cHit->GetY() - yav ;
 101       xxav     += xi * xi * cHit->GetXYWeight() ;
 102       xyav     += xi * yi * cHit->GetXYWeight() ;
 103       yyav     += yi * yi * cHit->GetXYWeight() ;
 104     }
 105
 106   if (fTrack->ComesFromMainVertex() == true)
 107     {
 108       xi        = fVertex->GetX() - xav ;
 109       yi        = fVertex->GetY() - yav ;
 110       xxav     += xi * xi * fVertex->GetXYWeight() ;
 111       xyav     += xi * yi * fVertex->GetXYWeight() ;
 112       yyav     += yi * yi * fVertex->GetXYWeight() ;
 113     }
 114   xxav = xxav / wsum ;
 115   xyav = xyav / wsum ;
 116   yyav = yyav / wsum ;
 117 //
 118 //-->  ROTATE COORDINATES SO THAT <XY> = 0
 119 //
 120 //-->  SIGN(C**2 - S**2) = SIGN(XXAV - YYAV) >
 121 //-->  &                                     > ==> NEW : (XXAV-YYAV) > 0
 122 //-->  SIGN(S) = SIGN(XYAV)                  >
 123
 124   Double_t a = fabs( xxav - yyav ) ;
 125   Double_t b = 4.0 * xyav * xyav ;
 126
 127   Double_t asqpb  = a * a + b  ;
 128   Double_t rasqpb = sqrt ( asqpb) ;
 129
 130   Double_t splus  = 1.0 + a / rasqpb ;
 131   Double_t sminus = b / (asqpb * splus) ;
 132
 133   splus  = sqrt (0.5 * splus ) ;
 134   sminus = sqrt (0.5 * sminus) ;
 135 //
 136 //->  FIRST REQUIRE : SIGN(C**2 - S**2) = SIGN(XXAV - YYAV)
 137 //
 138   Double_t sinrot, cosrot ;
 139   if ( xxav <= yyav ) {
 140          cosrot = sminus ;
 141          sinrot = splus  ;
 142   }
 143   else {
 144           cosrot = splus ;
 145           sinrot = sminus ;
 146   }
 147 //
 148 //->  REQUIRE : SIGN(S) = SIGN(XYAV) * SIGN(C) (ASSUMING SIGN(C) > 0)
 149 //
 150   if ( xyav < 0.0 ) sinrot = - sinrot ;
 151 //
 152 //-->  WE NOW HAVE THE SMALLEST ANGLE THAT GUARANTEES <X**2> > <Y**2>
 153 //-->  TO GET THE SIGN OF THE CHARGE RIGHT, THE NEW X-AXIS MUST POINT
 154 //-->  OUTWARD FROM THE ORGIN.  WE ARE FREE TO CHANGE SIGNS OF BOTH
 155 //-->  COSROT AND SINROT SIMULTANEOUSLY TO ACCOMPLISH THIS.
 156 //
 157 //-->  CHOOSE SIGN OF C WISELY TO BE ABLE TO GET THE SIGN OF THE CHARGE
 158 //
 159   if ( cosrot*xav+sinrot*yav < 0.0 ) {
 160           cosrot = -cosrot ;
 161           sinrot = -sinrot ;
 162   }
 163 //
 164 //->  NOW GET <R**2> AND RSCALE= SQRT(<R**2>)
 165 //
 166   Double_t rrav   = xxav + yyav ;
 167   Double_t rscale = sqrt(rrav) ;
 168
 169   xxav   = 0.0 ;
 170   yyav   = 0.0 ;
 171   xyav   = 0.0 ;
 172   Double_t xrrav = 0.0 ;
 173   Double_t yrrav = 0.0 ;
 174   Double_t rrrrav  = 0.0 ;
 175
 176   Double_t xixi, yiyi, riri, wiriri, xold, yold ;
 177
 178   //for (hit_counter=0; hit_counter<num_of_hits; hit_counter++)
 179   for(fTrack->StartLoop(); fTrack->LoopDone(); fTrack->GetNextHit())
 180     {
 181       //AliL3ConfMapPoint *cHit = (AliL3ConfMapPoint*)hits->At(hit_counter);
 182       AliL3ConfMapPoint* cHit = (AliL3ConfMapPoint*)fTrack->currentHit;
 183
 184       xold = cHit->GetX() - xav ;
 185       yold = cHit->GetY() - yav ;
 186       //
 187       //-->  ROTATE SO THAT <XY> = 0 & DIVIDE BY RSCALE SO THAT <R**2> = 1
 188       //
 189       xi = (  cosrot * xold + sinrot * yold ) / rscale ;
 190       yi = ( -sinrot * xold + cosrot * yold ) / rscale ;
 191
 192       xixi   = xi * xi ;
 193       yiyi   = yi * yi ;
 194       riri   = xixi + yiyi ;
 195       wiriri = cHit->GetXYWeight() * riri ;
 196
 197       xyav   += cHit->GetXYWeight() * xi * yi ;
 198       xxav   += cHit->GetXYWeight() * xixi ;
 199       yyav   += cHit->GetXYWeight() * yiyi ;
 200
 201       xrrav  += wiriri * xi ;
 202       yrrav  += wiriri * yi ;
 203       rrrrav += wiriri * riri ;
 204     }
 205   //
 206 //   Include vertex if required
 207 //
 208   if (fTrack->ComesFromMainVertex() == true)
 209     {
 210         xold = fVertex->GetX() - xav ;
 211         yold = fVertex->GetY() - yav ;
 212         //
 213         //-->  ROTATE SO THAT <XY> = 0 & DIVIDE BY RSCALE SO THAT <R**2> = 1
 214         //
 215         xi = (  cosrot * xold + sinrot * yold ) / rscale ;
 216         yi = ( -sinrot * xold + cosrot * yold ) / rscale ;
 217
 218         xixi   = xi * xi ;
 219         yiyi   = yi * yi ;
 220         riri   = xixi + yiyi ;
 221         wiriri = fVertex->GetXYWeight() * riri ;
 222
 223         xyav   += fVertex->GetXYWeight() * xi * yi ;
 224         xxav   += fVertex->GetXYWeight() * xixi ;
 225         yyav   += fVertex->GetXYWeight() * yiyi ;
 226
 227         xrrav  += wiriri * xi ;
 228         yrrav  += wiriri * yi ;
 229         rrrrav += wiriri * riri ;
 230   }
 231   //
 232   //
 233   //
 234   //-->  DIVIDE BY WSUM TO MAKE AVERAGES
 235   //
 236   xxav    = xxav   / wsum ;
 237   yyav    = yyav   / wsum ;
 238   xrrav   = xrrav  / wsum ;
 239   yrrav   = yrrav  / wsum ;
 240   rrrrav  = rrrrav / wsum ;
 241   xyav    = xyav   / wsum ;
 242
 243   Int_t const ntry = 5 ;
 244 //
 245 //-->  USE THESE TO GET THE COEFFICIENTS OF THE 4-TH ORDER POLYNIMIAL
 246 //-->  DON'T PANIC - THE THIRD ORDER TERM IS ZERO !
 247 //
 248   Double_t xrrxrr = xrrav * xrrav ;
 249   Double_t yrryrr = yrrav * yrrav ;
 250   Double_t rrrrm1 = rrrrav - 1.0  ;
 251   Double_t xxyy   = xxav  * yyav  ;
 252
 253   Double_t c0  =          rrrrm1*xxyy - xrrxrr*yyav - yrryrr*xxav ;
 254   Double_t c1  =        - rrrrm1      + xrrxrr      + yrryrr   - 4.0*xxyy ;
 255   Double_t c2  =   4.0  + rrrrm1                               - 4.0*xxyy ;
 256   Double_t c4  = - 4.0  ;
 257 //
 258 //-->  COEFFICIENTS OF THE DERIVATIVE - USED IN NEWTON-RAPHSON ITERATIONS
 259 //
 260   Double_t c2d =   2.0 * c2 ;
 261   Double_t c4d =   4.0 * c4 ;
 262 //
 263 //-->  0'TH VALUE OF LAMDA - LINEAR INTERPOLATION BETWEEN P(0) & P(YYAV)
 264 //
 265 //   LAMDA = YYAV * C0 / (C0 + YRRSQ * (XXAV-YYAV))
 266   Double_t lamda  = 0.0 ;
 267   Double_t dlamda = 0.0 ;
 268 //
 269   Double_t chiscl = wsum * rscale * rscale ;
 270   Double_t dlamax = 0.001 / chiscl ;
 271
 272   Double_t p, pd ;
 273   for ( int itry = 1 ; itry <= ntry ; itry++ ) {
 274      p      = c0 + lamda * (c1 + lamda * (c2 + lamda * lamda * c4 )) ;
 275      pd     = (c1 + lamda * (c2d + lamda * lamda * c4d)) ;
 276      dlamda = -p / pd ;
 277      lamda  = lamda + dlamda ;
 278      if (fabs(dlamda)<   dlamax) break ;
 279   }
 280
 281   Double_t chi2 = (Double_t)(chiscl * lamda) ;
 282
 283   fTrack->SetChiSq1(chi2);
 284   // Double_t dchisq = chiscl * dlamda ;
 285 //
 286 //-->  NOW CALCULATE THE MATRIX ELEMENTS FOR ALPHA, BETA & KAPPA
 287 //
 288   Double_t h11   = xxav  -     lamda ;
 289   Double_t h14   = xrrav ;
 290   Double_t h22   = yyav  -     lamda ;
 291   Double_t h24   = yrrav ;
 292   Double_t h34   = 1.0   + 2.0*lamda ;
 293   if ( h11 == 0.0 || h22 == 0.0 ){
 294     LOG(AliL3Log::kError,"AliL3ConfMapFit::FitCircle","TrackFit")<<AliL3Log::kDec<<
 295       "Problems fitting circle"<<ENDLOG;
 296     return 1 ;
 297   }
 298   Double_t rootsq = (h14*h14)/(h11*h11) + 4.0*h34 ;
 299
 300   Double_t ratio, kappa, beta ;
 301   if ( fabs(h22) > fabs(h24) ) {
 302      ratio  = h24 / h22 ;
 303      rootsq = ratio * ratio + rootsq ;
 304      kappa = 1.0 / sqrt(rootsq) ;
 305      beta  = - ratio * kappa ;
 306   }
 307   else {
 308      ratio  = h22 / h24 ;
 309      rootsq = 1.0 + ratio * ratio * rootsq ;
 310      beta  = 1.0 / sqrt(rootsq) ;
 311      if ( h24 > 0 ) beta = - beta ;
 312      kappa = -ratio * beta ;
 313   }
 314   Double_t alpha = - (h14/h11) * kappa ;
 315 //
 316 //-->  transform these into the lab coordinate system
 317 //-->  first get kappa and back to real dimensions
 318 //
 319   Double_t kappa1 = kappa / rscale ;
 320   Double_t dbro   = 0.5   / kappa1 ;
 321 //
 322 //-->  next rotate alpha and beta and scale
 323 //
 324   Double_t alphar = (cosrot * alpha - sinrot * beta)* dbro ;
 325   Double_t betar  = (sinrot * alpha + cosrot * beta)* dbro ;
 326 //
 327 //-->  then translate by (xav,yav)
 328 //
 329   Double_t acent  = (double)(xav - alphar) ;
 330   Double_t bcent  = (double)(yav - betar ) ;
 331   Double_t radius = (double)dbro ;
 332 //
 333 //   Get charge
 334 //
 335   Int_t q = ( ( yrrav < 0 ) ? 1 : -1 ) ;
 336
 337   fTrack->SetCharge(q);
 338
 339 //
 340 //    Get other track parameters
 341 //
 342   Double_t x0, y0,phi0,r0,psi,pt ;
 343   if ( fTrack->ComesFromMainVertex() == true )
 344     {
 345       //flag = 1 ; // primary track flag
 346       x0   = fVertex->GetX() ;
 347       y0   = fVertex->GetY() ;
 348       phi0 = fVertex->GetPhi() ;
 349       r0   = fVertex->GetR() ;
 350       fTrack->SetPhi0(phi0);
 351       fTrack->SetR0(r0);
 352     }
 353   else
 354     {
 355       //AliL3ConfMapPoint *lHit = (AliL3ConfMapPoint*)hits->Last();
 356       AliL3ConfMapPoint *lHit = (AliL3ConfMapPoint*)fTrack->lastHit;
 357       //flag =  0 ; // primary track flag
 358       x0   =  lHit->GetX()  ;
 359       y0   =  lHit->GetY()  ;
 360       phi0 =  atan2(lHit->GetY(),lHit->GetX());
 361       if ( phi0 < 0 ) phi0 += 2*pi;
 362       r0   =  sqrt ( lHit->GetX() * lHit->GetX() + lHit->GetY() * lHit->GetY() )  ;
 363       fTrack->SetPhi0(phi0);
 364       fTrack->SetR0(r0);
 365     }
 366   //
 367   psi  = (Double_t)atan2(bcent-y0,acent-x0) ;
 368   psi  = psi + q * 0.5F * pi ;
 369   if ( psi < 0 ) psi = psi + 2*pi;
 370
 371   pt   = (Double_t)(BFACT * bField * radius ) ;
 372   fTrack->SetPsi(psi);
 373   fTrack->SetPt(pt);
 374
 375   //
 376 //    Get errors from fast fit
 377 //
 378   //if ( getPara()->getErrors ) getErrorsCircleFit ( acent, bcent, radius ) ;
 379 //
 380   return 0 ;
 381
 382 }
 383
 384 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 385 //    Fit Line in s-z plane
 386 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 387 Int_t AliL3ConfMapFit::FitLine ( )
 388 {
 389   //
 390   //Initialization
 391   //
 392   Double_t sum = 0.F ;
 393   Double_t ss  = 0.F ;
 394   Double_t sz  = 0.F ;
 395   Double_t sss = 0.F ;
 396   Double_t ssz = 0.F ;
 397   //
 398   //find sum , sums ,sumz, sumss
 399   //
 400   Double_t dx, dy ;
 401   Double_t radius = (Double_t)(fTrack->GetPt() / ( BFACT * bField ) ) ;
 402
 403   //TObjArray *hits = fTrack->GetHits();
 404   //Int_t num_of_hits = fTrack->GetNumberOfPoints();
 405
 406   if ( fTrack->ComesFromMainVertex() == true )
 407     {
 408       dx = ((AliL3ConfMapPoint*)fTrack->firstHit)->GetX() - fVertex->GetX();
 409       dy = ((AliL3ConfMapPoint*)fTrack->firstHit)->GetY() - fVertex->GetY() ;
 410     }
 411   else
 412     {
 413       dx = ((AliL3ConfMapPoint *)fTrack->firstHit)->GetX() - ((AliL3ConfMapPoint *)fTrack->lastHit)->GetX() ;
 414       dy = ((AliL3ConfMapPoint *)fTrack->firstHit)->GetY() - ((AliL3ConfMapPoint *)fTrack->lastHit)->GetY() ;
 415       //dx = ((AliL3ConfMapPoint *)hits->First())->GetX() - ((AliL3ConfMapPoint *)hits->Last())->GetX() ;
 416       //dy = ((AliL3ConfMapPoint *)hits->First())->GetY() - ((AliL3ConfMapPoint *)hits->Last())->GetY() ;
 417     }
 418
 419   Double_t localPsi = 0.5F * sqrt ( dx*dx + dy*dy ) / radius ;
 420   Double_t total_s ;
 421
 422   if ( fabs(localPsi) < 1. )
 423     {
 424       total_s = 2.0 * radius * asin ( localPsi ) ;
 425     }
 426   else
 427     {
 428       total_s = 2.0 * radius * pi ;
 429     }
 430
 431   AliL3ConfMapPoint *previousHit = NULL;
 432
 433   // FtfBaseHit *previousHit = 0  ;
 434
 435   //for ( startLoop() ; done() ; nextHit() ) {
 436   Double_t dpsi,s;
 437
 438   //  for(hit_counter=0; hit_counter<num_of_hits; hit_counter++)
 439   for(fTrack->StartLoop(); fTrack->LoopDone(); fTrack->GetNextHit())
 440     {
 441       // AliL3ConfMapPoint *cHit = (AliL3ConfMapPoint*)hits->At(hit_counter);
 442       AliL3ConfMapPoint *cHit = (AliL3ConfMapPoint*)fTrack->currentHit;
 443       // if ( currentHit != firstHit )
 444       if(cHit != fTrack->firstHit)//  hits->First())
 445         {
 446           dx   = cHit->GetX() - previousHit->GetX() ;
 447           dy   = cHit->GetY() - previousHit->GetY() ;
 448           dpsi = 0.5 * (Double_t)sqrt ( dx*dx + dy*dy ) / radius ;
 449           fTrack->SetPsierr(dpsi);
 450           s = previousHit->GetS() - 2.0 * radius * (Double_t)asin ( dpsi ) ;
 451           cHit->SetS(s);
 452         }
 453       else
 454         cHit->SetS(total_s);
 455       //        cHit->s = total_s ;
 456
 457       sum += cHit->GetZWeight() ;
 458       ss  += cHit->GetZWeight() * cHit->GetS() ;
 459       sz  += cHit->GetZWeight() * cHit->GetZ() ;
 460       sss += cHit->GetZWeight() * cHit->GetS() * cHit->GetS() ;
 461       ssz += cHit->GetZWeight() * cHit->GetS() * cHit->GetZ() ;
 462       previousHit = cHit ;
 463     }
 464
 465   Double_t chi2,det = sum * sss - ss * ss;
 466   if ( fabs(det) < 1e-20)
 467     {
 468       chi2 = 99999.F ;
 469       fTrack->SetChiSq2(chi2);
 470       return 0 ;
 471     }
 472
 473   //Compute the best fitted parameters A,B
 474   Double_t tanl,z0,dtanl,dz0;
 475
 476   tanl = (Double_t)((sum * ssz - ss * sz ) / det );
 477   z0   = (Double_t)((sz * sss - ssz * ss ) / det );
 478
 479   fTrack->SetTgl(tanl);
 480   fTrack->SetZ0(z0);
 481
 482   //     calculate chi-square
 483
 484   chi2 = 0.;
 485   Double_t r1 ;
 486
 487   //for(hit_counter=0; hit_counter<num_of_hits; hit_counter++)
 488   for(fTrack->StartLoop(); fTrack->LoopDone(); fTrack->GetNextHit())
 489     {
 490       //AliL3ConfMapPoint *cHit = (AliL3ConfMapPoint*)hits->At(hit_counter);
 491       AliL3ConfMapPoint *cHit = (AliL3ConfMapPoint*)fTrack->currentHit;
 492       r1   = cHit->GetZ() - tanl * cHit->GetS() - z0 ;
 493       chi2 += (Double_t) ( (Double_t)cHit->GetZWeight() * (r1 * r1) );
 494     }
 495   fTrack->SetChiSq2(chi2);
 496   //
 497   //     calculate estimated variance
 498   //      varsq=chi/(double(n)-2.)
 499   //     calculate covariance matrix
 500   //      siga=sqrt(varsq*sxx/det)
 501   //      sigb=sqrt(varsq*sum/det)
 502   //
 503   dtanl = (Double_t) ( sum / det );
 504   dz0   = (Double_t) ( sss / det );
 505
 506   fTrack->SetTglerr(dtanl);
 507   fTrack->SetZ0err(dz0);
 508
 509   return 0 ;
 510 }