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 "AliTrackFitterRieman.h"
  36 #include "AliLog.h"
  37 #include "TTreeStream.h"
  38 #include "AliRieman.h"
  39 #include "AliLog.h"
  40
  41 ClassImp(AliTrackFitterRieman)
  42
  43
  44 AliTrackFitterRieman::AliTrackFitterRieman():
  45   AliTrackFitter(),
  46   fAlpha(0.),
  47   fNUsed(0),
  48   fConv(kFALSE),
  49   fMaxDelta(3),
  50   fRieman(new AliRieman(10000)),  // allocate rieman
  51   fDebugStream(new TTreeSRedirector("RiemanAlignDebug.root"))
  52 {
  53   //
  54   // default constructor
  55   //
  56 }
  57
  58
  59 AliTrackFitterRieman::AliTrackFitterRieman(AliTrackPointArray *array, Bool_t owner):
  60   AliTrackFitter(array,owner),
  61   fAlpha(0.),
  62   fNUsed(0),
  63   fConv(kFALSE),
  64   fMaxDelta(3),
  65   fRieman(new AliRieman(10000)),  //allocate rieman
  66   fDebugStream(0)
  67 {
  68   //
  69   // Constructor
  70   //
  71   if (AliLog::GetDebugLevel("","AliTrackFitterRieman")) fDebugStream = new TTreeSRedirector("RiemanAlignDebug.root");
  72 }
  73
  74 AliTrackFitterRieman::AliTrackFitterRieman(const AliTrackFitterRieman &rieman):
  75   AliTrackFitter(rieman),
  76   fAlpha(rieman.fAlpha),
  77   fNUsed(rieman.fNUsed),
  78   fConv(rieman.fConv),
  79   fMaxDelta(rieman.fMaxDelta),
  80   fRieman(new AliRieman(*(rieman.fRieman))),
  81   fDebugStream(0)
  82 {
  83   //
  84   // copy constructor
  85   //
  86   if (AliLog::GetDebugLevel("","AliTrackFitterRieman")) fDebugStream = new TTreeSRedirector("RiemanAlignDebug.root");
  87 }
  88
  89 //_____________________________________________________________________________
  90 AliTrackFitterRieman &AliTrackFitterRieman::operator =(const AliTrackFitterRieman& rieman)
  91 {
  92   //
  93   // Assignment operator
  94   //
  95   if(this==&rieman) return *this;
  96   ((AliTrackFitter *)this)->operator=(rieman);
  97
  98   fAlpha  = rieman.fAlpha;
  99   fNUsed  = rieman.fNUsed;
 100   fConv   = rieman.fConv;
 101   fMaxDelta = rieman.fMaxDelta;
 102   fRieman = new AliRieman(*(rieman.fRieman));
 103   fDebugStream = 0;
 104   if (AliLog::GetDebugLevel("","AliTrackFitterRieman")) fDebugStream = new TTreeSRedirector("RiemanAlignDebug.root");
 105   return *this;
 106 }
 107
 108
 109 AliTrackFitterRieman::~AliTrackFitterRieman(){
 110   //
 111   //
 112   //
 113   delete fRieman;
 114   delete fDebugStream;
 115 }
 116
 117 void AliTrackFitterRieman::Reset()
 118 {
 119   // Reset the track parameters and
 120   // rieman sums
 121   //
 122   AliTrackFitter::Reset();
 123   fRieman->Reset();
 124   fAlpha = 0.;
 125   fNUsed = 0;
 126   fConv =kFALSE;
 127 }
 128
 129 Bool_t AliTrackFitterRieman::Fit(const TArrayI *volIds,const TArrayI *volIdsFit,
 130                                  AliAlignObj::ELayerID layerRangeMin,
 131                          AliAlignObj::ELayerID layerRangeMax)
 132 {
 133   // Fit the track points. The method takes as an input
 134   // the set of id's (volids) of the volumes in which
 135   // one wants to calculate the residuals.
 136   // The following parameters are used to define the
 137   // range of volumes to be used in the fitting
 138   // As a result two AliTrackPointArray's obects are filled.
 139   // The first one contains the space points with
 140   // volume id's from volids list. The second array of points represents
 141   // the track extrapolations corresponding to the space points
 142   // in the first array. The two arrays can be used to find
 143   // the residuals in the volids and consequently construct a
 144   // chi2 function to be minimized during the alignment
 145   // procedures. For the moment the track extrapolation is taken
 146   // at the space-point reference plane. The reference plane is
 147   // found using the covariance matrix of the point
 148   // (assuming sigma(x)=0 at the reference coordinate system.
 149
 150   Reset();
 151   Int_t npoints = fPoints->GetNPoints();
 152   if (fPoints && AliLog::GetDebugLevel("","AliTrackFitterRieman")>1){
 153     Int_t nVol    = volIds->GetSize();
 154     Int_t nVolFit = volIdsFit->GetSize();
 155     Int_t volId  = volIds->At(0);
 156     (*fDebugStream)<<"PInput"<<
 157       "NPoints="<<npoints<<    // number of points
 158       "VolId="<<volId<<        // first vol ID
 159       "NVol="<<nVol<<          // number of volumes
 160       "NvolFit="<<nVolFit<<    // number of volumes to fit
 161       "fPoints.="<<fPoints<<   // input points
 162       "\n";
 163   }
 164   if (npoints < 3) return kFALSE;
 165
 166   Bool_t isAlphaCalc = kFALSE;
 167   AliTrackPoint p,plocal;
 168 //   fPoints->GetPoint(p,0);
 169 //   fAlpha = TMath::ATan2(p.GetY(),p.GetX());
 170
 171   Int_t npVolId = 0;
 172   fNUsed = 0;
 173   Int_t *pindex = new Int_t[npoints];
 174   for (Int_t ipoint = 0; ipoint < npoints; ipoint++)
 175     {
 176       fPoints->GetPoint(p,ipoint);
 177       UShort_t iVolId = p.GetVolumeID();
 178       if (FindVolId(volIds,iVolId)) {
 179         pindex[npVolId] = ipoint;
 180         npVolId++;
 181       }
 182       if (volIdsFit != 0x0) {
 183         if (!FindVolId(volIdsFit,iVolId)) continue;
 184       }
 185       else {
 186         if (iVolId < AliAlignObj::LayerToVolUID(layerRangeMin,0) ||
 187             iVolId > AliAlignObj::LayerToVolUID(layerRangeMax,
 188                                                 AliAlignObj::LayerSize(layerRangeMax))) continue;
 189       }
 190       if (!isAlphaCalc) {
 191         fAlpha = p.GetAngle();
 192         isAlphaCalc = kTRUE;
 193       }
 194       plocal = p.Rotate(fAlpha);
 195       if (TMath::Abs(plocal.GetX())>500 || TMath::Abs(plocal.GetX())<2 || plocal.GetCov()[3]<=0 ||plocal.GetCov()[5]<=0 ){
 196         printf("<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<</n");
 197         p.Dump();
 198         plocal.Dump();
 199         printf("Problematic point\n");
 200         printf("<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<</n");
 201       }else{
 202         AddPoint(plocal.GetX(),plocal.GetY(),plocal.GetZ(),
 203                  TMath::Sqrt(plocal.GetCov()[3]),TMath::Sqrt(plocal.GetCov()[5]));
 204       }
 205       // fNUsed++;  AddPoint should be responsible
 206     }
 207
 208   if (npVolId == 0 || fNUsed < fMinNPoints) {
 209     delete [] pindex;
 210     return kFALSE;
 211   }
 212
 213   Update();
 214
 215   if (!fConv) {
 216     delete [] pindex;
 217     return kFALSE;
 218   }
 219
 220   if ((fParams[0] == 0) ||
 221       ((-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1) <= 0)) {
 222     delete [] pindex;
 223     return kFALSE;
 224   }
 225
 226
 227   if (fNUsed < fMinNPoints) {
 228     delete [] pindex;
 229     return kFALSE;
 230   }
 231
 232   fPVolId = new AliTrackPointArray(npVolId);
 233   fPTrack = new AliTrackPointArray(npVolId);
 234   AliTrackPoint p2;
 235   for (Int_t ipoint = 0; ipoint < npVolId; ipoint++)
 236     {
 237       Int_t index = pindex[ipoint];
 238       fPoints->GetPoint(p,index);
 239       if (GetPCA(p,p2) && (
 240           TMath::Abs(p.GetX()-p2.GetX())<fMaxDelta &&
 241           TMath::Abs(p.GetY()-p2.GetY())<fMaxDelta &&
 242           TMath::Abs(p.GetZ()-p2.GetZ())<fMaxDelta
 243           )) {
 244         Float_t xyz[3],xyz2[3];
 245         p.GetXYZ(xyz); p2.GetXYZ(xyz2);
 246         //      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]);
 247         fPVolId->AddPoint(ipoint,&p);
 248         fPTrack->AddPoint(ipoint,&p2);
 249       }else{
 250         // what should be default bahavior -
 251         delete [] pindex;
 252         delete fPVolId;
 253         delete fPTrack;
 254         fPVolId =0;
 255         fPTrack =0;
 256         return kFALSE;
 257       }
 258     }
 259
 260   if (AliLog::GetDebugLevel("","AliTrackFitterRieman")>0){
 261     //
 262     // Debug Info
 263     //
 264     AliTrackPointArray *lPVolId = new AliTrackPointArray(npVolId);
 265     AliTrackPointArray *lPTrack = new AliTrackPointArray(npVolId);
 266     AliRieman * residual = fRieman->MakeResiduals();
 267     AliTrackPoint p2local;
 268     for (Int_t ipoint = 0; ipoint < npVolId; ipoint++){
 269       Int_t index = pindex[ipoint];
 270       fPoints->GetPoint(p,index);
 271       if (GetPCA(p,p2)) {
 272         plocal = p.Rotate(fAlpha);
 273         //      plocal.Rotate(fAlpha);
 274         Float_t xyz[3],xyz2[3];
 275         p2local = plocal;
 276         plocal.GetXYZ(xyz);
 277         xyz2[0] = xyz[0];
 278         xyz2[1] = GetYat(xyz[0]);
 279         xyz2[2] = GetZat(xyz[0]);
 280         p2local.SetXYZ(xyz2);
 281         lPVolId->AddPoint(ipoint,&plocal);
 282         lPTrack->AddPoint(ipoint,&p2local);
 283       }
 284     }
 285     //
 286     // debug info
 287     //
 288     Int_t nVol    = volIds->GetSize();
 289     Int_t nVolFit = volIdsFit->GetSize();
 290     Int_t volId   = volIds->At(0);
 291     Int_t modId   =0;
 292     Int_t layer   =  AliAlignObj::VolUIDToLayer(volId,modId);
 293
 294     (*fDebugStream)<<"Fit"<<
 295       "VolId="<<volId<<        // volume ID
 296       "Layer="<<layer<<        // layer ID
 297       "Module="<<modId<<       // module ID
 298       "NVol="<<nVol<<          // number of volumes
 299       "NvolFit="<<nVolFit<<    // number of volumes to fit
 300       "Points0.="<<fPVolId<<    // original points
 301       "Points1.="<<fPTrack<<    // fitted points
 302       "LPoints0.="<<lPVolId<<   // original points - local frame
 303       "LPoints1.="<<lPTrack<<   // fitted points   - local frame
 304       "Rieman.="<<this<<        // original rieman fit
 305       "Res.="<<residual<<       // residuals of rieman fit
 306       "\n";
 307     delete lPVolId;
 308     delete lPTrack;
 309     delete residual;
 310   }
 311
 312   delete [] pindex;
 313   return kTRUE;
 314 }
 315
 316
 317 void AliTrackFitterRieman::AddPoint(Float_t x, Float_t y, Float_t z, Float_t sy, Float_t sz)
 318 {
 319   //
 320   // add point to rieman fitter
 321   //
 322   fRieman->AddPoint(x,y,z,sy,sz);
 323   fNUsed = fRieman->GetN();
 324 }
 325
 326
 327
 328 void AliTrackFitterRieman::Update(){
 329   //
 330   //
 331   //
 332   fRieman->Update();
 333   fConv = kFALSE;
 334   if (fRieman->IsValid()){
 335     for (Int_t ipar=0; ipar<6; ipar++){
 336       fParams[ipar] = fRieman->GetParam()[ipar];
 337     }
 338     fChi2  =  fRieman->GetChi2();
 339     fNdf   =  fRieman->GetN()- 2;
 340     fNUsed =  fRieman->GetN();
 341     fConv = kTRUE;
 342   }
 343 }
 344
 345 //_____________________________________________________________________________
 346 Bool_t AliTrackFitterRieman::GetPCA(const AliTrackPoint &p, AliTrackPoint &p2) const
 347 {
 348   //
 349   // Get the closest to a given spacepoint track trajectory point
 350   // Look for details in the description of the Fit() method
 351   //
 352   if (!fConv) return kFALSE;
 353
 354   // First X and Y coordinates
 355   Double_t sin = TMath::Sin(fAlpha);
 356   Double_t cos = TMath::Cos(fAlpha);
 357   //  fParam[0]  = 1/y0
 358   //  fParam[1]  = -x0/y0
 359   //  fParam[2]  = - (R^2 - x0^2 - y0^2)/y0
 360   if (fParams[0] == 0) return kFALSE;
 361   // Track parameters in the global coordinate system
 362   Double_t x0 = -fParams[1]/fParams[0]*cos -         1./fParams[0]*sin;
 363   Double_t y0 =          1./fParams[0]*cos - fParams[1]/fParams[0]*sin;
 364   if ((-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1) <= 0) return kFALSE;
 365   Double_t r  = TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1)/
 366                 fParams[0];
 367
 368   // Define space-point refence plane
 369   Double_t alphap = p.GetAngle();
 370   Double_t sinp = TMath::Sin(alphap);
 371   Double_t cosp = TMath::Cos(alphap);
 372   Double_t x  = p.GetX()*cosp + p.GetY()*sinp;
 373   Double_t y  = p.GetY()*cosp - p.GetX()*sinp;
 374   Double_t x0p= x0*cosp + y0*sinp;
 375   Double_t y0p= y0*cosp - x0*sinp;
 376   if ((r*r - (x-x0p)*(x-x0p))<0) {
 377     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));
 378     return kFALSE;
 379   }
 380   Double_t temp = TMath::Sqrt(r*r - (x-x0p)*(x-x0p));
 381   Double_t y1 = y0p + temp;
 382   Double_t y2 = y0p - temp;
 383   Double_t yprime = y1;
 384   if(TMath::Abs(y2-y) < TMath::Abs(y1-y)) yprime = y2;
 385
 386   // Back to the global coordinate system
 387   Double_t xsecond = x*cosp - yprime*sinp;
 388   Double_t ysecond = yprime*cosp + x*sinp;
 389
 390   // Now Z coordinate and track angles
 391   Double_t x2 = xsecond*cos + ysecond*sin;
 392   Double_t zsecond = GetZat(x2);
 393   Double_t dydx = GetDYat(x2);
 394   Double_t dzdx = GetDZat(x2);
 395
 396   // Fill the cov matrix of the track extrapolation point
 397   Double_t cov[6] = {0,0,0,0,0,0};
 398   Double_t sigmax = 100*100.;
 399   cov[0] = sigmax;           cov[1] = sigmax*dydx;      cov[2] = sigmax*dzdx;
 400   cov[3] = sigmax*dydx*dydx; cov[4] = sigmax*dydx*dzdx;
 401   cov[5] = sigmax*dzdx*dzdx;
 402
 403   Float_t  newcov[6];
 404   newcov[0] = cov[0]*cos*cos-
 405             2*cov[1]*sin*cos+
 406               cov[3]*sin*sin;
 407   newcov[1] = cov[1]*(cos*cos-sin*sin)-
 408              (cov[3]-cov[0])*sin*cos;
 409   newcov[2] = cov[2]*cos-
 410               cov[4]*sin;
 411   newcov[3] = cov[0]*sin*sin+
 412             2*cov[1]*sin*cos+
 413               cov[3]*cos*cos;
 414   newcov[4] = cov[4]*cos+
 415               cov[2]*sin;
 416   newcov[5] = cov[5];
 417
 418   p2.SetXYZ(xsecond,ysecond,zsecond,newcov);
 419   if (AliLog::GetDebugLevel("","AliTrackFitterRieman")>1){
 420     AliTrackPoint lp0(p);
 421     AliTrackPoint lp2(p2);
 422     if (0)(*fDebugStream)<<"PCA"<<
 423       "P0.="<<&lp0<<
 424       "P2.="<<&lp2<<
 425       "\n";
 426   }
 427   return kTRUE;
 428 }