HLT/src/AliL3ConfMapFit.cxx

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