HLT/TPCLib/AliHLTTPCConfMapFit.cxx

   1 // @(#) $Id$
   2 // Original: AliHLTConfMapFit.cxx,v 1.14 2005/06/14 10:55:21 cvetan
   3
   4 /**************************************************************************
   5  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   6  *                                                                        *
   7  * Authors: Anders Vestbo                                                 *
   8  *          Matthias Richter <Matthias.Richter@ift.uib.no>                *
   9  *          for The ALICE Off-line Project.                               *
  10  *                                                                        *
  11  * Permission to use, copy, modify and distribute this software and its   *
  12  * documentation strictly for non-commercial purposes is hereby granted   *
  13  * without fee, provided that the above copyright notice appears in all   *
  14  * copies and that both the copyright notice and this permission notice   *
  15  * appear in the supporting documentation. The authors make no claims     *
  16  * about the suitability of this software for any purpose. It is          *
  17  * provided "as is" without express or implied warranty.                  *
  18  **************************************************************************/
  19
  20 /** @file   AliHLTTPCConfMapFit.cxx
  21     @author Anders Vestbo, maintained by Matthias Richter
  22     @date
  23     @brief  Fit class for conformal mapping tracking.
  24 */
  25
  26 #include "AliHLTTPCRootTypes.h"
  27 #include "AliHLTTPCLogging.h"
  28 #include "AliHLTTPCVertex.h"
  29 #include "AliHLTTPCConfMapTrack.h"
  30 #include "AliHLTTPCConfMapPoint.h"
  31 #include "AliHLTTPCTransform.h"
  32 #include "AliHLTTPCConfMapFit.h"
  33
  34 #if __GNUC__ >= 3
  35 using namespace std;
  36 #endif
  37
  38 ClassImp(AliHLTTPCConfMapFit)
  39
  40
  41 AliHLTTPCConfMapFit::AliHLTTPCConfMapFit()
  42   :
  43   fTrack(NULL),
  44   fVertex(NULL)
  45 {
  46   //constructor
  47 }
  48
  49 AliHLTTPCConfMapFit::AliHLTTPCConfMapFit(AliHLTTPCConfMapTrack *track,AliHLTTPCVertex *vertex)
  50   :
  51   fTrack(track),
  52   fVertex(vertex)
  53
  54 {
  55   //constructor
  56 }
  57
  58 AliHLTTPCConfMapFit::AliHLTTPCConfMapFit(const AliHLTTPCConfMapFit&)
  59   :
  60   fTrack(NULL),
  61   fVertex(NULL)
  62 {
  63   // dummy copy constructor
  64   //HLTFatal("copy constructor untested");
  65 }
  66
  67 AliHLTTPCConfMapFit& AliHLTTPCConfMapFit::operator=(const AliHLTTPCConfMapFit&)
  68 {
  69   // dummy assignment operator
  70   //HLTFatal("assignment operator untested");
  71   return *this;
  72 }
  73
  74 AliHLTTPCConfMapFit::~AliHLTTPCConfMapFit()
  75 {
  76   // destructor
  77 }
  78
  79 Int_t AliHLTTPCConfMapFit::FitHelix()
  80 {
  81   //fit the helix
  82   if(FitCircle())
  83     {
  84       LOG(AliHLTTPCLog::kError,"AliHLTTPCConfMapFit::FitHelix","TrackFit")<<AliHLTTPCLog::kDec<<
  85         "Problems during circle fit"<<ENDLOG;
  86       return 1;
  87     }
  88   if(FitLine())
  89     {
  90       LOG(AliHLTTPCLog::kError,"AliHLTTPCConfMapFit::FitHelix","TrackFit")<<AliHLTTPCLog::kDec<<
  91         "Problems during line fit"<<ENDLOG;
  92       return 1;
  93     }
  94   return 0;
  95 }
  96
  97 // #### -B0-CHANGE-START == JMT
  98
  99 Int_t AliHLTTPCConfMapFit::FitStraightLine() {
 100     //fit the straight line
 101     if(FitLineXY()) {
 102         LOG(AliHLTTPCLog::kError,"AliHLTTPCConfMapFit::FitStraightLine","TrackFit")<<AliHLTTPCLog::kDec<<
 103             "Problems during stright line fit in XY plane"<<ENDLOG;
 104         return 1;
 105     }
 106     if(FitLineSZ()){
 107         LOG(AliHLTTPCLog::kError,"AliHLTTPCConfMapFit::FitStraightLine","TrackFit")<<AliHLTTPCLog::kDec<<
 108             "Problems during stright line fit in SZ plane"<<ENDLOG;
 109         return 1;
 110     }
 111     return 0;
 112 }
 113
 114 // #### -B0-CHANGE-END == JMT
 115
 116 Int_t AliHLTTPCConfMapFit::FitCircle()
 117 {
 118   //-----------------------------------------------------------------
 119   //Fits circle parameters using algorithm
 120   //described by ChErnov and Oskov in Computer Physics
 121   //Communications.
 122   //
 123   //Written in FORTRAN by Jawluen Tang, Physics department , UT-Austin
 124   //Moved to C by Pablo Yepes
 125   //Moved to AliROOT by ASV.
 126   //------------------------------------------------------------------
 127
 128   Double_t wsum  = 0.0 ;
 129   Double_t xav   = 0.0 ;
 130   Double_t yav   = 0.0 ;
 131
 132   Int_t num_of_hits = fTrack->GetNumberOfPoints();
 133   //
 134   //     Loop over hits calculating average
 135   Int_t co=0;
 136
 137   for(fTrack->StartLoop(); fTrack->LoopDone(); fTrack->GetNextHit())
 138     {
 139       co++;
 140       AliHLTTPCConfMapPoint *cHit = (AliHLTTPCConfMapPoint*)fTrack->GetCurrentHit();
 141       cHit->SetXYWeight( 1./ (Double_t)(cHit->GetXerr()*cHit->GetXerr() + cHit->GetYerr()*cHit->GetYerr()) );
 142       wsum      += cHit->GetXYWeight() ;
 143       xav       += cHit->GetXYWeight() * cHit->GetX() ;
 144       yav       += cHit->GetXYWeight() * cHit->GetY() ;
 145     }
 146   if(co!=num_of_hits)
 147     LOG(AliHLTTPCLog::kError,"AliHLTTPCConfMapFit::FitCircle","TrackFit")<<AliHLTTPCLog::kDec<<
 148       "Mismatch of hits. Counter: "<<co<<" nHits: "<<num_of_hits<<ENDLOG;
 149   if (fTrack->ComesFromMainVertex() == true)
 150     {
 151       wsum += fVertex->GetXYWeight() ;
 152       xav  += fVertex->GetX() ;
 153       yav  += fVertex->GetY() ;
 154     }
 155
 156   xav = xav / wsum ;
 157   yav = yav / wsum ;
 158 //
 159 //  CALCULATE <X**2>, <XY>, AND <Y**2> WITH <X> = 0, & <Y> = 0
 160 //
 161   Double_t xxav  = 0.0 ;
 162   Double_t xyav  = 0.0 ;
 163   Double_t yyav  = 0.0 ;
 164   Double_t xi, yi ;
 165
 166   for(fTrack->StartLoop(); fTrack->LoopDone(); fTrack->GetNextHit())
 167     {
 168       //AliHLTTPCConfMapPoint *cHit = (AliHLTTPCConfMapPoint *)hits->At(hit_counter);
 169       AliHLTTPCConfMapPoint *cHit = (AliHLTTPCConfMapPoint*)fTrack->GetCurrentHit();
 170       xi        = cHit->GetX() - xav ;
 171       yi        = cHit->GetY() - yav ;
 172       xxav     += xi * xi * cHit->GetXYWeight() ;
 173       xyav     += xi * yi * cHit->GetXYWeight() ;
 174       yyav     += yi * yi * cHit->GetXYWeight() ;
 175     }
 176
 177   if (fTrack->ComesFromMainVertex() == true)
 178     {
 179       xi        = fVertex->GetX() - xav ;
 180       yi        = fVertex->GetY() - yav ;
 181       xxav     += xi * xi * fVertex->GetXYWeight() ;
 182       xyav     += xi * yi * fVertex->GetXYWeight() ;
 183       yyav     += yi * yi * fVertex->GetXYWeight() ;
 184     }
 185   xxav = xxav / wsum ;
 186   xyav = xyav / wsum ;
 187   yyav = yyav / wsum ;
 188 //
 189 //-->  ROTATE COORDINATES SO THAT <XY> = 0
 190 //
 191 //-->  SIGN(C**2 - S**2) = SIGN(XXAV - YYAV) >
 192 //-->  &                                     > ==> NEW : (XXAV-YYAV) > 0
 193 //-->  SIGN(S) = SIGN(XYAV)                  >
 194
 195   Double_t a = fabs( xxav - yyav ) ;
 196   Double_t b = 4.0 * xyav * xyav ;
 197
 198   Double_t asqpb  = a * a + b  ;
 199   Double_t rasqpb = sqrt ( asqpb) ;
 200
 201   Double_t splus  = 1.0 + a / rasqpb ;
 202   Double_t sminus = b / (asqpb * splus) ;
 203
 204   splus  = sqrt (0.5 * splus ) ;
 205   sminus = sqrt (0.5 * sminus) ;
 206 //
 207 //->  FIRST REQUIRE : SIGN(C**2 - S**2) = SIGN(XXAV - YYAV)
 208 //
 209   Double_t sinrot, cosrot ;
 210   if ( xxav <= yyav ) {
 211          cosrot = sminus ;
 212          sinrot = splus  ;
 213   }
 214   else {
 215           cosrot = splus ;
 216           sinrot = sminus ;
 217   }
 218 //
 219 //->  REQUIRE : SIGN(S) = SIGN(XYAV) * SIGN(C) (ASSUMING SIGN(C) > 0)
 220 //
 221   if ( xyav < 0.0 ) sinrot = - sinrot ;
 222 //
 223 //-->  WE NOW HAVE THE SMALLEST ANGLE THAT GUARANTEES <X**2> > <Y**2>
 224 //-->  TO GET THE SIGN OF THE CHARGE RIGHT, THE NEW X-AXIS MUST POINT
 225 //-->  OUTWARD FROM THE ORGIN.  WE ARE FREE TO CHANGE SIGNS OF BOTH
 226 //-->  COSROT AND SINROT SIMULTANEOUSLY TO ACCOMPLISH THIS.
 227 //
 228 //-->  CHOOSE SIGN OF C WISELY TO BE ABLE TO GET THE SIGN OF THE CHARGE
 229 //
 230   if ( cosrot*xav+sinrot*yav < 0.0 ) {
 231           cosrot = -cosrot ;
 232           sinrot = -sinrot ;
 233   }
 234 //
 235 //->  NOW GET <R**2> AND RSCALE= SQRT(<R**2>)
 236 //
 237   Double_t rrav   = xxav + yyav ;
 238   Double_t rscale = sqrt(rrav) ;
 239
 240   xxav   = 0.0 ;
 241   yyav   = 0.0 ;
 242   xyav   = 0.0 ;
 243   Double_t xrrav = 0.0 ;
 244   Double_t yrrav = 0.0 ;
 245   Double_t rrrrav  = 0.0 ;
 246
 247   Double_t xixi, yiyi, riri, wiriri, xold, yold ;
 248
 249   //for (hit_counter=0; hit_counter<num_of_hits; hit_counter++)
 250   for(fTrack->StartLoop(); fTrack->LoopDone(); fTrack->GetNextHit())
 251     {
 252       //AliHLTTPCConfMapPoint *cHit = (AliHLTTPCConfMapPoint*)hits->At(hit_counter);
 253       AliHLTTPCConfMapPoint* cHit = (AliHLTTPCConfMapPoint*)fTrack->GetCurrentHit();
 254
 255       xold = cHit->GetX() - xav ;
 256       yold = cHit->GetY() - yav ;
 257       //
 258       //-->  ROTATE SO THAT <XY> = 0 & DIVIDE BY RSCALE SO THAT <R**2> = 1
 259       //
 260       xi = (  cosrot * xold + sinrot * yold ) / rscale ;
 261       yi = ( -sinrot * xold + cosrot * yold ) / rscale ;
 262
 263       xixi   = xi * xi ;
 264       yiyi   = yi * yi ;
 265       riri   = xixi + yiyi ;
 266       wiriri = cHit->GetXYWeight() * riri ;
 267
 268       xyav   += cHit->GetXYWeight() * xi * yi ;
 269       xxav   += cHit->GetXYWeight() * xixi ;
 270       yyav   += cHit->GetXYWeight() * yiyi ;
 271
 272       xrrav  += wiriri * xi ;
 273       yrrav  += wiriri * yi ;
 274       rrrrav += wiriri * riri ;
 275     }
 276   //
 277 //   Include vertex if required
 278 //
 279   if (fTrack->ComesFromMainVertex() == true)
 280     {
 281         xold = fVertex->GetX() - xav ;
 282         yold = fVertex->GetY() - yav ;
 283         //
 284         //-->  ROTATE SO THAT <XY> = 0 & DIVIDE BY RSCALE SO THAT <R**2> = 1
 285         //
 286         xi = (  cosrot * xold + sinrot * yold ) / rscale ;
 287         yi = ( -sinrot * xold + cosrot * yold ) / rscale ;
 288
 289         xixi   = xi * xi ;
 290         yiyi   = yi * yi ;
 291         riri   = xixi + yiyi ;
 292         wiriri = fVertex->GetXYWeight() * riri ;
 293
 294         xyav   += fVertex->GetXYWeight() * xi * yi ;
 295         xxav   += fVertex->GetXYWeight() * xixi ;
 296         yyav   += fVertex->GetXYWeight() * yiyi ;
 297
 298         xrrav  += wiriri * xi ;
 299         yrrav  += wiriri * yi ;
 300         rrrrav += wiriri * riri ;
 301   }
 302   //
 303   //
 304   //
 305   //-->  DIVIDE BY WSUM TO MAKE AVERAGES
 306   //
 307   xxav    = xxav   / wsum ;
 308   yyav    = yyav   / wsum ;
 309   xrrav   = xrrav  / wsum ;
 310   yrrav   = yrrav  / wsum ;
 311   rrrrav  = rrrrav / wsum ;
 312   xyav    = xyav   / wsum ;
 313
 314   Int_t const ntry = 5 ;
 315 //
 316 //-->  USE THESE TO GET THE COEFFICIENTS OF THE 4-TH ORDER POLYNIMIAL
 317 //-->  DON'T PANIC - THE THIRD ORDER TERM IS ZERO !
 318 //
 319   Double_t xrrxrr = xrrav * xrrav ;
 320   Double_t yrryrr = yrrav * yrrav ;
 321   Double_t rrrrm1 = rrrrav - 1.0  ;
 322   Double_t xxyy   = xxav  * yyav  ;
 323
 324   Double_t c0  =          rrrrm1*xxyy - xrrxrr*yyav - yrryrr*xxav ;
 325   Double_t c1  =        - rrrrm1      + xrrxrr      + yrryrr   - 4.0*xxyy ;
 326   Double_t c2  =   4.0  + rrrrm1                               - 4.0*xxyy ;
 327   Double_t c4  = - 4.0  ;
 328 //
 329 //-->  COEFFICIENTS OF THE DERIVATIVE - USED IN NEWTON-RAPHSON ITERATIONS
 330 //
 331   Double_t c2d =   2.0 * c2 ;
 332   Double_t c4d =   4.0 * c4 ;
 333 //
 334 //-->  0'TH VALUE OF LAMDA - LINEAR INTERPOLATION BETWEEN P(0) & P(YYAV)
 335 //
 336 //   LAMDA = YYAV * C0 / (C0 + YRRSQ * (XXAV-YYAV))
 337   Double_t lamda  = 0.0 ;
 338   Double_t dlamda = 0.0 ;
 339 //
 340   Double_t chiscl = wsum * rscale * rscale ;
 341   Double_t dlamax = 0.001 / chiscl ;
 342
 343   Double_t p, pd ;
 344   for ( int itry = 1 ; itry <= ntry ; itry++ ) {
 345      p      = c0 + lamda * (c1 + lamda * (c2 + lamda * lamda * c4 )) ;
 346      pd     = (c1 + lamda * (c2d + lamda * lamda * c4d)) ;
 347      dlamda = -p / pd ;
 348      lamda  = lamda + dlamda ;
 349      if (fabs(dlamda)<   dlamax) break ;
 350   }
 351
 352   Double_t chi2 = (Double_t)(chiscl * lamda) ;
 353
 354   fTrack->SetChiSq1(chi2);
 355   // Double_t dchisq = chiscl * dlamda ;
 356 //
 357 //-->  NOW CALCULATE THE MATRIX ELEMENTS FOR ALPHA, BETA & KAPPA
 358 //
 359   Double_t h11   = xxav  -     lamda ;
 360   Double_t h14   = xrrav ;
 361   Double_t h22   = yyav  -     lamda ;
 362   Double_t h24   = yrrav ;
 363   Double_t h34   = 1.0   + 2.0*lamda ;
 364   if ( h11 == 0.0 || h22 == 0.0 ){
 365     LOG(AliHLTTPCLog::kError,"AliHLTTPCConfMapFit::FitCircle","TrackFit")<<AliHLTTPCLog::kDec<<
 366       "Problems fitting circle"<<ENDLOG;
 367     return 1 ;
 368   }
 369   Double_t rootsq = (h14*h14)/(h11*h11) + 4.0*h34 ;
 370
 371   Double_t ratio, kappa, beta ;
 372   if ( fabs(h22) > fabs(h24) ) {
 373      ratio  = h24 / h22 ;
 374      rootsq = ratio * ratio + rootsq ;
 375      kappa = 1.0 / sqrt(rootsq) ;
 376      beta  = - ratio * kappa ;
 377   }
 378   else {
 379      ratio  = h22 / h24 ;
 380      rootsq = 1.0 + ratio * ratio * rootsq ;
 381      beta  = 1.0 / sqrt(rootsq) ;
 382      if ( h24 > 0 ) beta = - beta ;
 383      kappa = -ratio * beta ;
 384   }
 385   Double_t alpha = - (h14/h11) * kappa ;
 386 //
 387 //-->  transform these into the lab coordinate system
 388 //-->  first get kappa and back to real dimensions
 389 //
 390   Double_t kappa1 = kappa / rscale ;
 391   Double_t dbro   = 0.5   / kappa1 ;
 392 //
 393 //-->  next rotate alpha and beta and scale
 394 //
 395   Double_t alphar = (cosrot * alpha - sinrot * beta)* dbro ;
 396   Double_t betar  = (sinrot * alpha + cosrot * beta)* dbro ;
 397 //
 398 //-->  then translate by (xav,yav)
 399 //
 400   Double_t acent  = (double)(xav - alphar) ;
 401   Double_t bcent  = (double)(yav - betar ) ;
 402   Double_t radius = (double)dbro ;
 403 //
 404 //   Get charge
 405 //
 406   Int_t q = ( ( yrrav < 0 ) ? 1 : -1 ) ;
 407
 408   fTrack->SetCharge(q);
 409
 410
 411   //Set the first point on the track to the space point coordinates of the innermost track
 412   //This will be updated to lie on the fit later on (AliHLTTPCTrack::UpdateToFirstPoint).
 413   Double_t x0,y0,psi,pt ;
 414   AliHLTTPCConfMapPoint *lHit = (AliHLTTPCConfMapPoint*)fTrack->GetLastHit();
 415   x0 = lHit->GetX();
 416   y0 = lHit->GetY();
 417   fTrack->SetFirstPoint(x0,y0,0); //Z-value is set in FitLine
 418
 419   psi  = (Double_t)atan2(bcent-y0,acent-x0) ;
 420   psi  = psi + q * AliHLTTPCTransform::PiHalf();
 421   if ( psi < 0 ) psi = psi + AliHLTTPCTransform::TwoPi();
 422   pt   = (Double_t)(AliHLTTPCTransform::GetBFieldValue() * radius ) ;
 423
 424   //Update the track parameters with the parameters from this fit:
 425   fTrack->SetPsi(psi);
 426   fTrack->SetPt(pt);
 427   fTrack->SetRadius(radius);
 428   fTrack->SetCenterX(acent);
 429   fTrack->SetCenterY(bcent);
 430
 431   //
 432 //    Get errors from fast fit
 433 //
 434   //if ( getPara()->getErrors ) getErrorsCircleFit ( acent, bcent, radius ) ;
 435 //
 436   return 0 ;
 437
 438 }
 439
 440 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 441 //    Fit Line in s-z plane
 442 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 443 Int_t AliHLTTPCConfMapFit::FitLine ( )
 444 {
 445   //
 446   //Initialization
 447   //
 448   Double_t sum = 0.F ;
 449   Double_t ss  = 0.F ;
 450   Double_t sz  = 0.F ;
 451   Double_t sss = 0.F ;
 452   Double_t ssz = 0.F ;
 453   //
 454   //find sum , sums ,sumz, sumss
 455   //
 456   Double_t dx, dy ;
 457   Double_t radius = (Double_t)(fTrack->GetPt() / AliHLTTPCTransform::GetBFieldValue() ) ;
 458
 459   //TObjArray *hits = fTrack->GetHits();
 460   //Int_t num_of_hits = fTrack->GetNumberOfPoints();
 461
 462   if (0)// fTrack->ComesFromMainVertex() == true )
 463     {
 464       dx = ((AliHLTTPCConfMapPoint*)fTrack->GetFirstHit())->GetX() - fVertex->GetX();
 465       dy = ((AliHLTTPCConfMapPoint*)fTrack->GetFirstHit())->GetY() - fVertex->GetY() ;
 466     }
 467   else
 468     {
 469       dx = ((AliHLTTPCConfMapPoint *)fTrack->GetFirstHit())->GetX() - ((AliHLTTPCConfMapPoint *)fTrack->GetLastHit())->GetX() ;
 470       dy = ((AliHLTTPCConfMapPoint *)fTrack->GetFirstHit())->GetY() - ((AliHLTTPCConfMapPoint *)fTrack->GetLastHit())->GetY() ;
 471       //dx = ((AliHLTTPCConfMapPoint *)hits->First())->GetX() - ((AliHLTTPCConfMapPoint *)hits->Last())->GetX() ;
 472       //dy = ((AliHLTTPCConfMapPoint *)hits->First())->GetY() - ((AliHLTTPCConfMapPoint *)hits->Last())->GetY() ;
 473     }
 474
 475   Double_t localPsi = 0.5F * sqrt ( dx*dx + dy*dy ) / radius ;
 476   Double_t total_s ;
 477
 478   if ( fabs(localPsi) < 1. )
 479     {
 480       total_s = 2.0 * radius * asin ( localPsi ) ;
 481     }
 482   else
 483     {
 484       total_s = 2.0 * radius * AliHLTTPCTransform::Pi() ;
 485     }
 486
 487   AliHLTTPCConfMapPoint *previousHit = NULL;
 488
 489   // FtfBaseHit *previousHit = 0  ;
 490
 491   //for ( startLoop() ; done() ; nextHit() ) {
 492   Double_t dpsi,s;
 493
 494   //  for(hit_counter=0; hit_counter<num_of_hits; hit_counter++)
 495   for(fTrack->StartLoop(); fTrack->LoopDone(); fTrack->GetNextHit())
 496     {
 497       // AliHLTTPCConfMapPoint *cHit = (AliHLTTPCConfMapPoint*)hits->At(hit_counter);
 498       AliHLTTPCConfMapPoint *cHit = (AliHLTTPCConfMapPoint*)fTrack->GetCurrentHit();
 499       // if ( GetCurrentHit() != GetFirstHit() )
 500       if(cHit != fTrack->GetFirstHit())//  hits->First())
 501         {
 502           dx   = cHit->GetX() - previousHit->GetX() ;
 503           dy   = cHit->GetY() - previousHit->GetY() ;
 504           dpsi = 0.5 * (Double_t)sqrt ( dx*dx + dy*dy ) / radius ;
 505           fTrack->SetPsierr(dpsi);
 506           s = previousHit->GetS() - 2.0 * radius * (Double_t)asin ( dpsi ) ;
 507           cHit->SetS(s);
 508         }
 509       else
 510         cHit->SetS(total_s);
 511       //        cHit->s = total_s ;
 512
 513       sum += cHit->GetZWeight() ;
 514       ss  += cHit->GetZWeight() * cHit->GetS() ;
 515       sz  += cHit->GetZWeight() * cHit->GetZ() ;
 516       sss += cHit->GetZWeight() * cHit->GetS() * cHit->GetS() ;
 517       ssz += cHit->GetZWeight() * cHit->GetS() * cHit->GetZ() ;
 518       previousHit = cHit ;
 519     }
 520
 521   Double_t chi2,det = sum * sss - ss * ss;
 522   if ( fabs(det) < 1e-20)
 523     {
 524       chi2 = 99999.F ;
 525       fTrack->SetChiSq2(chi2);
 526       return 0 ;
 527     }
 528
 529   //Compute the best fitted parameters A,B
 530   Double_t tanl,z0,dtanl,dz0;
 531
 532   tanl = (Double_t)((sum * ssz - ss * sz ) / det );
 533   z0   = (Double_t)((sz * sss - ssz * ss ) / det );
 534
 535   fTrack->SetTgl(tanl);
 536   fTrack->SetZ0(z0);
 537
 538   //     calculate chi-square
 539
 540   chi2 = 0.;
 541   Double_t r1 ;
 542
 543   //for(hit_counter=0; hit_counter<num_of_hits; hit_counter++)
 544   for(fTrack->StartLoop(); fTrack->LoopDone(); fTrack->GetNextHit())
 545     {
 546       //AliHLTTPCConfMapPoint *cHit = (AliHLTTPCConfMapPoint*)hits->At(hit_counter);
 547       AliHLTTPCConfMapPoint *cHit = (AliHLTTPCConfMapPoint*)fTrack->GetCurrentHit();
 548       r1   = cHit->GetZ() - tanl * cHit->GetS() - z0 ;
 549       chi2 += (Double_t) ( (Double_t)cHit->GetZWeight() * (r1 * r1) );
 550     }
 551   fTrack->SetChiSq2(chi2);
 552   //
 553   //     calculate estimated variance
 554   //      varsq=chi/(double(n)-2.)
 555   //     calculate covariance matrix
 556   //      siga=sqrt(varsq*sxx/det)
 557   //      sigb=sqrt(varsq*sum/det)
 558   //
 559   dtanl = (Double_t) ( sum / det );
 560   dz0   = (Double_t) ( sss / det );
 561
 562   fTrack->SetTglerr(dtanl);
 563   fTrack->SetZ0err(dz0);
 564
 565   return 0 ;
 566 }
 567
 568
 569 // #### -B0-CHANGE-START == JMT
 570
 571 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 572 //    Straight Line Fit  in x-y plane
 573 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 574 Int_t AliHLTTPCConfMapFit::FitLineXY ( ){
 575     // -----------------------------------------------------------------------------
 576     // Implementation after Numerical Recipes in C, 2nd Edtion, chapter 15.2, p. 661
 577     // with y = b*x + a
 578     // and Data Analysis for Physical Science Students, Luis Lyons, chapter 2.4 p.51
 579     // with y = a' + bx' , x' = x - <x>
 580     // -----------------------------------------------------------------------------
 581
 582     Double_t S = 0.;
 583     Double_t Sx = 0.;
 584
 585     Double_t SPrime = 0.;
 586     Double_t SxPrime = 0.;
 587     Double_t SxxPrime = 0.;
 588     Double_t SyPrime = 0.;
 589     Double_t SxyPrime = 0.;
 590
 591     Double_t chi2 = 0.;
 592
 593     Int_t num_of_hits = fTrack->GetNumberOfPoints();
 594
 595     Int_t co=0;
 596
 597     // - Loop over hits calculating average : xav
 598     for(fTrack->StartLoop(); fTrack->LoopDone(); fTrack->GetNextHit()) {
 599         co++;
 600         AliHLTTPCConfMapPoint *currentHit = (AliHLTTPCConfMapPoint*)fTrack->GetCurrentHit();
 601         // ** maybe not necessary, already done in ConfMapPoint
 602         currentHit->SetXYWeight( 1./ (Double_t)(currentHit->GetXerr()*currentHit->GetXerr() + currentHit->GetYerr()*currentHit->GetYerr()) );
 603         // **
 604         S   += currentHit->GetXYWeight();
 605         Sx  += currentHit->GetXYWeight() * currentHit->GetX();
 606     }
 607
 608     if(co!=num_of_hits)
 609         LOG(AliHLTTPCLog::kError,"AliHLTTPCConfMapFit::FitLineXY","TrackFit") << "Mismatch of hits. Counter: "<<co<<" nHits: "<<num_of_hits<<ENDLOG;
 610
 611     Double_t xav = (Double_t)Sx / S;
 612
 613     // Calculate weighted means
 614     for(fTrack->StartLoop(); fTrack->LoopDone(); fTrack->GetNextHit()) {
 615         AliHLTTPCConfMapPoint *currentHit = (AliHLTTPCConfMapPoint*)fTrack->GetCurrentHit();
 616
 617         Double_t xPrime =  currentHit->GetX() - xav;
 618         SPrime   += currentHit->GetXYWeight();
 619         SxPrime  += currentHit->GetXYWeight() * xPrime;
 620         SxxPrime += currentHit->GetXYWeight() * xPrime * xPrime;
 621         SyPrime  += currentHit->GetXYWeight() * currentHit->GetY();
 622         SxyPrime += currentHit->GetXYWeight() * xPrime * currentHit->GetY();
 623     }
 624
 625     Double_t det = SPrime*SxxPrime + SxPrime*SxPrime;
 626
 627     if (fabs(det) < 1e-20) {
 628         LOG(AliHLTTPCLog::kDebug,"AliHLTTPCConfMapFit::FitLineXY","TrackFit") << "Determinant == 0" << ENDLOG;
 629         chi2 = 99999.F ;
 630         fTrack->SetChiSq1(chi2);
 631         return -1 ;
 632     }
 633
 634     Double_t b   = (Double_t)(SPrime*SxyPrime - SxPrime*SyPrime) / det;        // line parameter b
 635     Double_t aPrime   = (Double_t)(SxxPrime*SyPrime - SxPrime*SxyPrime) / det; // line parameter a
 636
 637     Double_t sigma2b = (Double_t)1. / SxxPrime;
 638     //-- Double_t sigma2aprime = (Double_t)1. /SPrime;
 639
 640     // Get gradient angle psi of line in xy plane
 641     Double_t psi  = (Double_t) atan(b) ;
 642
 643     // Calculate chi2
 644     for(fTrack->StartLoop(); fTrack->LoopDone(); fTrack->GetNextHit()) {
 645         AliHLTTPCConfMapPoint *currentHit = (AliHLTTPCConfMapPoint*)fTrack->GetCurrentHit();
 646         Double_t tempchi = currentHit->GetY() - aPrime - b*(currentHit->GetX() - xav);
 647         chi2 += tempchi*tempchi*currentHit->GetXYWeight() ;
 648     }
 649
 650     Double_t a = aPrime - b*xav;
 651
 652
 653     // Set TrackParameter
 654     fTrack->SetChiSq1(chi2);
 655     fTrack->SetPsi(psi);
 656     fTrack->SetPsierr(sigma2b);
 657     fTrack->SetCenterX(0.);    // Set to point on the track (for UpdateToFirstPoint)
 658     fTrack->SetCenterY(a);     // Set to point on the track (for UpdateToFirstPoint)
 659
 660     //Set the first point on the track to the space point coordinates of the innermost track
 661     //This will be updated to lie on the fit later on (AliHLTTPCTrack::UpdateToFirstPoint).
 662     AliHLTTPCConfMapPoint *lastHit = (AliHLTTPCConfMapPoint*)fTrack->GetLastHit();
 663     Double_t x0 = lastHit->GetX();
 664     Double_t y0 = lastHit->GetY();
 665     fTrack->SetFirstPoint(x0,y0,0); //Z-value is set in FitLineSZ
 666
 667
 668     //Set Defaults
 669     fTrack->SetRadius(-1.);
 670     fTrack->SetCharge(1);
 671     fTrack->SetPt(-1.);
 672
 673
 674     return 0;
 675 }
 676
 677
 678 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 679 //    Straight Line Fit  in s-z plane
 680 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 681 Int_t AliHLTTPCConfMapFit::FitLineSZ ( ){
 682     // -----------------------------------------------------------------------------
 683     // Implementation after Numerical Recipes in C, 2nd Edtion, chapter 15.2, p. 661
 684     // with z = b*s + a
 685     // and Data Analysis for Physical Science Students, Luis Lyons, chapter 2.4 p.51
 686     // with z = a' + bs' , s' = s - <s>
 687     // -----------------------------------------------------------------------------
 688
 689     Double_t S = 0.;
 690     Double_t Ss = 0.;
 691
 692     Double_t SPrime = 0.;
 693     Double_t SsPrime = 0.;
 694     Double_t SssPrime = 0.;
 695     Double_t SzPrime = 0.;
 696     Double_t SszPrime = 0.;
 697
 698     Double_t chi2 = 0.;
 699     Double_t s = 0.;
 700
 701     AliHLTTPCConfMapPoint *previousHit = NULL;
 702
 703     // - Loop over hits calculating length in xy-plane: s
 704     // - Loop over hits calculating average : sav
 705     for(fTrack->StartLoop(); fTrack->LoopDone(); fTrack->GetNextHit()) {
 706         AliHLTTPCConfMapPoint *currentHit = (AliHLTTPCConfMapPoint*)fTrack->GetCurrentHit();
 707         if(currentHit != fTrack->GetFirstHit()) {
 708             Double_t dx = currentHit->GetX() - previousHit->GetX() ;
 709             Double_t dy = currentHit->GetY() - previousHit->GetY() ;
 710             s = previousHit->GetS() - (Double_t)sqrt ( dx*dx + dy*dy );
 711         }
 712         else{
 713             Double_t dx = ((AliHLTTPCConfMapPoint *)fTrack->GetFirstHit())->GetX() - ((AliHLTTPCConfMapPoint *)fTrack->GetLastHit())->GetX();
 714             Double_t dy = ((AliHLTTPCConfMapPoint *)fTrack->GetFirstHit())->GetY() - ((AliHLTTPCConfMapPoint *)fTrack->GetLastHit())->GetY();
 715             s = (Double_t)sqrt ( dx*dx + dy*dy );
 716         }
 717
 718         currentHit->SetS(s);
 719
 720         S   += currentHit->GetZWeight();
 721         Ss  += currentHit->GetZWeight() * currentHit->GetS();
 722     }
 723
 724     Double_t sav = (Double_t)Ss / S;
 725
 726     // Calculate weighted means
 727     for(fTrack->StartLoop(); fTrack->LoopDone(); fTrack->GetNextHit()) {
 728         AliHLTTPCConfMapPoint *currentHit = (AliHLTTPCConfMapPoint*)fTrack->GetCurrentHit();
 729
 730         Double_t sPrime =  currentHit->GetS() - sav;
 731         SPrime   += currentHit->GetZWeight();
 732         SsPrime  += currentHit->GetZWeight() * sPrime;
 733         SssPrime += currentHit->GetZWeight() * sPrime * sPrime;
 734         SzPrime  += currentHit->GetZWeight() * currentHit->GetZ();
 735         SszPrime += currentHit->GetZWeight() * sPrime * currentHit->GetZ();
 736     }
 737
 738     Double_t det = SPrime*SssPrime + SsPrime*SsPrime;
 739
 740     if (fabs(det) < 1e-20) {
 741         LOG(AliHLTTPCLog::kDebug,"AliHLTTPCConfMapFit::FitLineSZ","TrackFit") << "Determinant == 0" << ENDLOG;
 742         chi2 = 99999.F ;
 743         fTrack->SetChiSq2(chi2);
 744         return -1 ;
 745     }
 746
 747     Double_t b   = (Double_t)(SPrime*SszPrime - SsPrime*SzPrime) / det;        // line parameter b
 748     Double_t aPrime   = (Double_t)(SssPrime*SzPrime - SsPrime*SszPrime) / det; // line parameter a
 749
 750     Double_t a = aPrime - b*sav;
 751
 752     Double_t sigma2b = (Double_t) 1. / SssPrime;
 753     Double_t sigma2aprime = (Double_t) 1. /SPrime;
 754
 755     Double_t sigma2a = sigma2aprime + sav*sav * sigma2b*sigma2b;
 756
 757     // Calculate chi2
 758     for(fTrack->StartLoop(); fTrack->LoopDone(); fTrack->GetNextHit()) {
 759         AliHLTTPCConfMapPoint *currentHit = (AliHLTTPCConfMapPoint*)fTrack->GetCurrentHit();
 760         Double_t tempchi = currentHit->GetZ() - aPrime - b*(currentHit->GetS() - sav);
 761         chi2 += tempchi*tempchi*currentHit->GetZWeight() ;
 762     }
 763
 764     // Set TrackParameter
 765     fTrack->SetChiSq2(chi2);
 766     fTrack->SetTgl(b);
 767     fTrack->SetZ0(a);
 768     fTrack->SetTglerr(sigma2b);
 769 //  fTrack->SetZ0err(sigma2aprime);   // maybe subject to check
 770     fTrack->SetZ0err(sigma2a);        // maybe subject to check
 771     return 0;
 772 }
 773
 774 // #### -B0-CHANGE-END == JMT