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