STEER/AliTrackFitterRieman.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 ///////////////////////////////////////////////////////////////////////////////
  19 //
  20 // Class to the track points on the Riemann sphere. Inputs are
  21 // the set of id's (volids) of the volumes in which residuals are
  22 // calculated to construct a chi2 function to be minimized during
  23 // the alignment procedures. For the moment the track extrapolation is
  24 // taken at the space-point reference plane. The reference plane is
  25 // found using the covariance matrix of the point
  26 // (assuming sigma(x)=0 at the reference coordinate system.
  27 //
  28 // Internal usage of AliRieman class for minimization
  29 //
  30 //////////////////////////////////////////////////////////////////////////////
  31
  32 #include "TMatrixDSym.h"
  33 #include "TMatrixD.h"
  34 #include "TArrayI.h"
  35 #include "TLinearFitter.h"
  36 #include "AliTrackFitterRieman.h"
  37 #include "AliLog.h"
  38 #include "TTreeStream.h"
  39 #include "AliRieman.h"
  40 #include "AliLog.h"
  41 #include "TRandom.h"
  42
  43 ClassImp(AliTrackFitterRieman)
  44
  45
  46 AliTrackFitterRieman::AliTrackFitterRieman():
  47   AliTrackFitter(),
  48   fBCorrection(kFALSE),
  49   fAlpha(0.),
  50   fNUsed(0),
  51   fConv(kFALSE),
  52   fMaxDelta(3),
  53   fRieman(new AliRieman(10000)),  // allocate rieman
  54   fDebugStream(new TTreeSRedirector("RiemanAlignDebug.root"))
  55 {
  56   //
  57   // default constructor
  58   //
  59 }
  60
  61
  62 AliTrackFitterRieman::AliTrackFitterRieman(AliTrackPointArray *array, Bool_t owner):
  63   AliTrackFitter(array,owner),
  64   fBCorrection(kFALSE),
  65   fAlpha(0.),
  66   fNUsed(0),
  67   fConv(kFALSE),
  68   fMaxDelta(3),
  69   fRieman(new AliRieman(10000)),  //allocate rieman
  70   fDebugStream(0)
  71 {
  72   //
  73   // Constructor
  74   //
  75   if (AliLog::GetDebugLevel("","AliTrackFitterRieman")) fDebugStream = new TTreeSRedirector("RiemanAlignDebug.root");
  76 }
  77
  78 AliTrackFitterRieman::AliTrackFitterRieman(const AliTrackFitterRieman &rieman):
  79   AliTrackFitter(rieman),
  80   fBCorrection(rieman.fBCorrection),
  81   fAlpha(rieman.fAlpha),
  82   fNUsed(rieman.fNUsed),
  83   fConv(rieman.fConv),
  84   fMaxDelta(rieman.fMaxDelta),
  85   fRieman(new AliRieman(*(rieman.fRieman))),
  86   fDebugStream(0)
  87 {
  88   //
  89   // copy constructor
  90   //
  91   if (AliLog::GetDebugLevel("","AliTrackFitterRieman")) fDebugStream = new TTreeSRedirector("RiemanAlignDebug.root");
  92 }
  93
  94 //_____________________________________________________________________________
  95 AliTrackFitterRieman &AliTrackFitterRieman::operator =(const AliTrackFitterRieman& rieman)
  96 {
  97   //
  98   // Assignment operator
  99   //
 100   if(this==&rieman) return *this;
 101   ((AliTrackFitter *)this)->operator=(rieman);
 102
 103   fBCorrection = rieman.fBCorrection;
 104   fAlpha  = rieman.fAlpha;
 105   fNUsed  = rieman.fNUsed;
 106   fConv   = rieman.fConv;
 107   fMaxDelta = rieman.fMaxDelta;
 108   fRieman = new AliRieman(*(rieman.fRieman));
 109   fDebugStream = 0;
 110   if (AliLog::GetDebugLevel("","AliTrackFitterRieman")) fDebugStream = new TTreeSRedirector("RiemanAlignDebug.root");
 111   return *this;
 112 }
 113
 114
 115 AliTrackFitterRieman::~AliTrackFitterRieman(){
 116   //
 117   //
 118   //
 119   delete fRieman;
 120   delete fDebugStream;
 121 }
 122
 123 void AliTrackFitterRieman::Reset()
 124 {
 125   // Reset the track parameters and
 126   // rieman sums
 127   //
 128   AliTrackFitter::Reset();
 129   fRieman->Reset();
 130   fAlpha = 0.;
 131   fNUsed = 0;
 132   fConv =kFALSE;
 133 }
 134
 135 Bool_t AliTrackFitterRieman::Fit(const TArrayI *volIds,const TArrayI *volIdsFit,
 136                                  AliAlignObj::ELayerID layerRangeMin,
 137                          AliAlignObj::ELayerID layerRangeMax)
 138 {
 139   // Fit the track points. The method takes as an input
 140   // the set of id's (volids) of the volumes in which
 141   // one wants to calculate the residuals.
 142   // The following parameters are used to define the
 143   // range of volumes to be used in the fitting
 144   // As a result two AliTrackPointArray's obects are filled.
 145   // The first one contains the space points with
 146   // volume id's from volids list. The second array of points represents
 147   // the track extrapolations corresponding to the space points
 148   // in the first array. The two arrays can be used to find
 149   // the residuals in the volids and consequently construct a
 150   // chi2 function to be minimized during the alignment
 151   // procedures. For the moment the track extrapolation is taken
 152   // at the space-point reference plane. The reference plane is
 153   // found using the covariance matrix of the point
 154   // (assuming sigma(x)=0 at the reference coordinate system.
 155   Int_t debugLevel = AliLog::GetDebugLevel("","AliTrackFitterRieman");
 156   Float_t debugRatio = 1./(1.+debugLevel);
 157
 158   const Int_t kMinPoints =1;
 159   Int_t npoints = fPoints->GetNPoints();
 160   if ( npoints<fMinNPoints) return kFALSE;
 161   //
 162   // fast count points
 163   Int_t countFit = 0;
 164   for (Int_t ifit=0; ifit<volIdsFit->GetSize(); ifit++){
 165     Int_t volIdFit  = volIdsFit->At(ifit);
 166     for (Int_t ipoint = 0; ipoint < npoints; ipoint++)
 167       if (volIdFit==fPoints->GetVolumeID()[ipoint]) countFit++;
 168   }
 169   if (countFit<fMinNPoints) return kFALSE;
 170   //
 171   Int_t countPoint = 0;
 172   for (Int_t jpoint=0; jpoint<volIds->GetSize(); jpoint++){
 173     Int_t volIdPoint  = volIds->At(jpoint);
 174     for (Int_t ipoint = 0; ipoint < npoints; ipoint++)
 175       if (volIdPoint==fPoints->GetVolumeID()[ipoint]) countPoint++;
 176   }
 177   if (countPoint<kMinPoints) return kFALSE;
 178   //
 179   //
 180
 181   Reset();
 182
 183   if (fPoints && AliLog::GetDebugLevel("","AliTrackFitterRieman")>1&& gRandom->Rndm()<debugRatio){
 184     Int_t nVol    = volIds->GetSize();
 185     Int_t nVolFit = volIdsFit->GetSize();
 186     Int_t volId  = volIds->At(0);
 187     (*fDebugStream)<<"PInput"<<
 188       "NPoints="<<npoints<<    // number of points
 189       "VolId="<<volId<<        // first vol ID
 190       "NVol="<<nVol<<          // number of volumes
 191       "NvolFit="<<nVolFit<<    // number of volumes to fit
 192       "fPoints.="<<fPoints<<   // input points
 193       "\n";
 194   }
 195
 196   Bool_t isAlphaCalc = kFALSE;
 197   AliTrackPoint p,plocal;
 198 //   fPoints->GetPoint(p,0);
 199 //   fAlpha = TMath::ATan2(p.GetY(),p.GetX());
 200
 201   Int_t npVolId = 0;
 202   fNUsed = 0;
 203   Int_t *pindex = new Int_t[npoints];
 204   for (Int_t ipoint = 0; ipoint < npoints; ipoint++)
 205     {
 206       fPoints->GetPoint(p,ipoint);
 207       UShort_t iVolId = p.GetVolumeID();
 208       if (FindVolId(volIds,iVolId)) {
 209         pindex[npVolId] = ipoint;
 210         npVolId++;
 211       }
 212       if (volIdsFit != 0x0) {
 213         if (!FindVolId(volIdsFit,iVolId)) continue;
 214       }
 215       else {
 216         if (iVolId < AliAlignObj::LayerToVolUID(layerRangeMin,0) ||
 217             iVolId > AliAlignObj::LayerToVolUID(layerRangeMax,
 218                                                 AliAlignObj::LayerSize(layerRangeMax))) continue;
 219       }
 220       if (!isAlphaCalc) {
 221         fAlpha = p.GetAngle();
 222         isAlphaCalc = kTRUE;
 223       }
 224       plocal = p.Rotate(fAlpha);
 225       if (TMath::Abs(plocal.GetX())>500 || TMath::Abs(plocal.GetX())<2 || plocal.GetCov()[3]<=0 ||plocal.GetCov()[5]<=0 ){
 226         printf("<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<</n");
 227         p.Dump();
 228         plocal.Dump();
 229         printf("Problematic point\n");
 230         printf("<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<</n");
 231       }else{
 232         AddPoint(plocal.GetX(),plocal.GetY(),plocal.GetZ(),
 233                  TMath::Sqrt(plocal.GetCov()[3]),TMath::Sqrt(plocal.GetCov()[5]));
 234       }
 235       // fNUsed++;  AddPoint should be responsible
 236     }
 237
 238   if (npVolId == 0 || fNUsed < fMinNPoints) {
 239     delete [] pindex;
 240     return kFALSE;
 241   }
 242
 243   Update();
 244
 245   if (!fConv) {
 246     delete [] pindex;
 247     return kFALSE;
 248   }
 249
 250   if ((fParams[0] == 0) ||
 251       ((-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1) <= 0)) {
 252     delete [] pindex;
 253     return kFALSE;
 254   }
 255
 256
 257   if (fNUsed < fMinNPoints) {
 258     delete [] pindex;
 259     return kFALSE;
 260   }
 261
 262   fPVolId = new AliTrackPointArray(npVolId);
 263   fPTrack = new AliTrackPointArray(npVolId);
 264   AliTrackPoint p2;
 265   for (Int_t ipoint = 0; ipoint < npVolId; ipoint++)
 266     {
 267       Int_t index = pindex[ipoint];
 268       fPoints->GetPoint(p,index);
 269       if (GetPCA(p,p2) && (
 270           TMath::Abs(p.GetX()-p2.GetX())<fMaxDelta &&
 271           TMath::Abs(p.GetY()-p2.GetY())<fMaxDelta &&
 272           TMath::Abs(p.GetZ()-p2.GetZ())<fMaxDelta
 273           )) {
 274         Float_t xyz[3],xyz2[3];
 275         p.GetXYZ(xyz); p2.GetXYZ(xyz2);
 276         //      printf("residuals %f %d %d %f %f %f %f %f %f\n",fChi2,fNUsed,fConv,xyz[0],xyz[1],xyz[2],xyz2[0]-xyz[0],xyz2[1]-xyz[1],xyz2[2]-xyz[2]);
 277         fPVolId->AddPoint(ipoint,&p);
 278         fPTrack->AddPoint(ipoint,&p2);
 279       }else{
 280         // what should be default bahavior -
 281         delete [] pindex;
 282         delete fPVolId;
 283         delete fPTrack;
 284         fPVolId =0;
 285         fPTrack =0;
 286         return kFALSE;
 287       }
 288     }
 289
 290   if (AliLog::GetDebugLevel("","AliTrackFitterRieman")>0 && gRandom->Rndm()<debugRatio){
 291     //
 292     // Debug Info
 293     //
 294     AliTrackPointArray *lPVolId = new AliTrackPointArray(npVolId);
 295     AliTrackPointArray *lPTrack = new AliTrackPointArray(npVolId);
 296     AliTrackPointArray *lPTrackE = new AliTrackPointArray(npVolId);
 297     AliRieman * residual = fRieman->MakeResiduals();
 298     for (Int_t ipoint = 0; ipoint < npVolId; ipoint++){
 299       AliTrackPoint p0, p0local;
 300       AliTrackPoint pFit, pFitlocal, pFitLocalE;
 301       fPVolId->GetPoint(p0,ipoint);
 302       Float_t lAngle = p0.GetAngle();
 303       p0local= p0.MasterToLocal();
 304       fPTrack->GetPoint(pFit,ipoint);
 305       pFitlocal= pFit.Rotate(lAngle);
 306       //
 307       Float_t xyz[3], cov[3];
 308       xyz[0] = pFitlocal.GetX();
 309       xyz[1] = pFitlocal.GetY();
 310       xyz[2] = pFitlocal.GetZ();
 311       for (Int_t icov=0; icov<6; icov++) cov[icov]=0;
 312       cov[3] = GetErrY2at(xyz[0]);
 313       cov[5] = GetErrZ2at(xyz[0]);
 314       pFitLocalE.SetXYZ(xyz,cov);
 315       //
 316       lPVolId->AddPoint(ipoint,&p0local);
 317       lPTrack->AddPoint(ipoint,&pFitlocal);
 318       lPTrackE->AddPoint(ipoint,&pFitLocalE);
 319     }
 320     //
 321     // debug info
 322     //
 323     Int_t nVol    = volIds->GetSize();
 324     Int_t nVolFit = volIdsFit->GetSize();
 325     Int_t volId   = volIds->At(0);
 326     Int_t modId   =0;
 327     Int_t layer   =  AliAlignObj::VolUIDToLayer(volId,modId);
 328     Int_t volIdFit   = volIdsFit->At(0);
 329     Int_t modIdFit   =0;
 330     Int_t layerFit   =  AliAlignObj::VolUIDToLayer(volIdFit,modIdFit);
 331
 332     (*fDebugStream)<<"Fit"<<
 333       "VolId="<<volId<<        // volume ID
 334       "Layer="<<layer<<        // layer ID
 335       "Module="<<modId<<       // module ID
 336       "LayerFit="<<layerFit<<        // layer ID fit
 337       "ModuleFit="<<modIdFit<<       // module ID fit
 338       "NVol="<<nVol<<          // number of volumes
 339       "NvolFit="<<nVolFit<<    // number of volumes to fit
 340       "Points0.="<<fPVolId<<    // original points
 341       "Points1.="<<fPTrack<<    // fitted points
 342       "LPoints0.="<<lPVolId<<   // original points - local frame
 343       "LPoints1.="<<lPTrack<<   // fitted points   - local frame
 344       "LPointsE.="<<lPTrackE<<   // fitted points with ext error - local frame
 345       "Rieman.="<<this<<        // original rieman fit
 346       "Res.="<<residual<<       // residuals of rieman fit
 347       "\n";
 348     delete lPVolId;
 349     delete lPTrack;
 350     delete residual;
 351   }
 352
 353   delete [] pindex;
 354   return kTRUE;
 355 }
 356
 357
 358 void AliTrackFitterRieman::AddPoint(Float_t x, Float_t y, Float_t z, Float_t sy, Float_t sz)
 359 {
 360   //
 361   // add point to rieman fitter
 362   //
 363   fRieman->AddPoint(x,y,z,sy,sz);
 364   fNUsed = fRieman->GetN();
 365 }
 366
 367
 368
 369 void AliTrackFitterRieman::Update(){
 370   //
 371   //
 372   //
 373   fRieman->Update();
 374   fConv = kFALSE;
 375   if (fRieman->IsValid()){
 376     for (Int_t ipar=0; ipar<6; ipar++){
 377       fParams[ipar] = fRieman->GetParam()[ipar];
 378     }
 379     fChi2  =  fRieman->GetChi2();
 380     fNdf   =  fRieman->GetN()- 2;
 381     fNUsed =  fRieman->GetN();
 382     fConv = kTRUE;
 383   }
 384   //
 385   //
 386   //
 387   TLinearFitter fitY(3,"pol2");
 388   TLinearFitter fitZ(3,"pol2");
 389   for (Int_t ip=0; ip<fRieman->GetN();ip++){
 390     Double_t x = fRieman->GetX()[ip];
 391     fitY.AddPoint(&x,fRieman->GetY()[ip]-fRieman->GetYat(x),1);
 392     fitZ.AddPoint(&x,fRieman->GetZ()[ip]-fRieman->GetZat(x),1);
 393   }
 394   fitY.Eval();
 395   fitZ.Eval();
 396   for (Int_t iparam=0; iparam<3; iparam++){
 397     fCorrY[iparam]=fitY.GetParameter(iparam);
 398     fCorrZ[iparam]=fitZ.GetParameter(iparam);
 399   }
 400   fCorrY[3]=fitY.GetChisquare()/Float_t(fRieman->GetN()-3);
 401   fCorrZ[3]=fitZ.GetChisquare()/Float_t(fRieman->GetN()-3);
 402 }
 403
 404
 405
 406 //_____________________________________________________________________________
 407 Bool_t AliTrackFitterRieman::GetPCA(const AliTrackPoint &p, AliTrackPoint &p2) const
 408 {
 409   //
 410   // Get the closest to a given spacepoint track trajectory point
 411   // Look for details in the description of the Fit() method
 412   //
 413   if (!fConv) return kFALSE;
 414
 415   // First X and Y coordinates
 416   Double_t sin = TMath::Sin(fAlpha);
 417   Double_t cos = TMath::Cos(fAlpha);
 418   //  fParam[0]  = 1/y0
 419   //  fParam[1]  = -x0/y0
 420   //  fParam[2]  = - (R^2 - x0^2 - y0^2)/y0
 421   if (fParams[0] == 0) return kFALSE;
 422   // Track parameters in the global coordinate system
 423   Double_t x0 = -fParams[1]/fParams[0]*cos -         1./fParams[0]*sin;
 424   Double_t y0 =          1./fParams[0]*cos - fParams[1]/fParams[0]*sin;
 425   if ((-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1) <= 0) return kFALSE;
 426   Double_t r  = TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1)/
 427                 fParams[0];
 428
 429   // Define space-point refence plane
 430   Double_t alphap = p.GetAngle();
 431   Double_t sinp = TMath::Sin(alphap);
 432   Double_t cosp = TMath::Cos(alphap);
 433   Double_t x  = p.GetX()*cosp + p.GetY()*sinp;
 434   Double_t y  = p.GetY()*cosp - p.GetX()*sinp;
 435   Double_t x0p= x0*cosp + y0*sinp;
 436   Double_t y0p= y0*cosp - x0*sinp;
 437   if ((r*r - (x-x0p)*(x-x0p))<0) {
 438     AliWarning(Form("Track extrapolation failed ! (Track radius = %f, track circle x = %f, space-point x = %f, reference plane angle = %f\n",r,x0p,x,alphap));
 439     return kFALSE;
 440   }
 441   Double_t temp = TMath::Sqrt(r*r - (x-x0p)*(x-x0p));
 442   Double_t y1 = y0p + temp;
 443   Double_t y2 = y0p - temp;
 444   Double_t yprime = y1;
 445   if(TMath::Abs(y2-y) < TMath::Abs(y1-y)) yprime = y2;
 446
 447   // Back to the global coordinate system
 448   Double_t xsecond = x*cosp - yprime*sinp;
 449   Double_t ysecond = yprime*cosp + x*sinp;
 450
 451   // Now Z coordinate and track angles
 452   Double_t x2 = xsecond*cos + ysecond*sin;
 453   Double_t zsecond = GetZat(x2);
 454   Double_t dydx = GetDYat(x2);
 455   Double_t dzdx = GetDZat(x2);
 456
 457   // Fill the cov matrix of the track extrapolation point
 458   Double_t cov[6] = {0,0,0,0,0,0};
 459   Double_t sigmax = 100*100.;
 460   cov[0] = sigmax;           cov[1] = sigmax*dydx;      cov[2] = sigmax*dzdx;
 461   cov[3] = sigmax*dydx*dydx; cov[4] = sigmax*dydx*dzdx;
 462   cov[5] = sigmax*dzdx*dzdx;
 463
 464   Double_t sigmay2 = GetErrY2at(x2);
 465   Double_t sigmaz2 = GetErrZ2at(x2);
 466   cov[3] += sigmay2;
 467   cov[5] += sigmaz2;
 468
 469
 470   Float_t  newcov[6];
 471   newcov[0] = cov[0]*cos*cos-
 472             2*cov[1]*sin*cos+
 473               cov[3]*sin*sin;
 474   newcov[1] = cov[1]*(cos*cos-sin*sin)-
 475              (cov[3]-cov[0])*sin*cos;
 476   newcov[2] = cov[2]*cos-
 477               cov[4]*sin;
 478   newcov[3] = cov[0]*sin*sin+
 479             2*cov[1]*sin*cos+
 480               cov[3]*cos*cos;
 481   newcov[4] = cov[4]*cos+
 482               cov[2]*sin;
 483   newcov[5] = cov[5];
 484
 485   p2.SetXYZ(xsecond,ysecond,zsecond,newcov);
 486   Int_t debugLevel = AliLog::GetDebugLevel("","AliTrackFitterRieman");
 487   Float_t debugRatio = 1./(1.+debugLevel);
 488   if (AliLog::GetDebugLevel("","AliTrackFitterRieman")>0 && gRandom->Rndm()<debugRatio){
 489     AliTrackPoint lp0(p);
 490     AliTrackPoint lp2(p2);
 491     (*fDebugStream)<<"PCA"<<
 492       "P0.="<<&lp0<<
 493       "P2.="<<&lp2<<
 494       "\n";
 495   }
 496   return kTRUE;
 497 }
 498
 499 Double_t AliTrackFitterRieman::GetYat(Double_t x) const  {
 500   //
 501   // get y position at given point
 502   //
 503   Double_t correction=0;
 504   if (fBCorrection){ // systematic effect correction
 505     correction  =  fCorrY[0]+fCorrY[1]*x +fCorrY[2]*x*x;
 506   }
 507   return fRieman->GetYat(x)+correction;
 508 }
 509
 510 Double_t AliTrackFitterRieman::GetZat(Double_t x) const  {
 511   //
 512   // get z position at given point
 513   //
 514   Double_t correction=0;
 515   if (fBCorrection){ // systematic effect correction
 516     correction  =  fCorrZ[0]+fCorrZ[1]*x +fCorrZ[2]*x*x;
 517   }
 518   return fRieman->GetZat(x)+correction;
 519 }
 520
 521 Double_t AliTrackFitterRieman::GetErrY2at(Double_t x) const {
 522   //
 523   // get estimate of extrapolation error
 524   //
 525   Double_t error = fRieman->GetErrY(x);
 526   Double_t correction=0;
 527   if (fBCorrection){ // everestimate error due systematic effect
 528     error      *=fCorrY[3];
 529     correction  =  fCorrY[0]+fCorrY[1]*x +fCorrY[2]*x*x;
 530     correction *=correction;
 531   }
 532   return TMath::Sqrt(error+correction);
 533 }
 534
 535 Double_t AliTrackFitterRieman::GetErrZ2at(Double_t x) const {
 536   //
 537   // get estimate of extrapolation error
 538   //
 539   Double_t error = fRieman->GetErrZ(x)*fCorrZ[3];
 540   Double_t correction=0;
 541   if (fBCorrection){
 542      error    *=  fCorrZ[3];
 543     correction =  fCorrZ[0]+fCorrZ[1]*x +fCorrZ[2]*x*x;
 544     correction*=  correction;
 545   }
 546   return TMath::Sqrt(error+correction);
 547 }
 548