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