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