TPC/AliTPCkalmanFit.cxx

   1 /*
   2   marian.ivanov@cern.ch
   3
   4   AliTPCkalmanFit:
   5
   6   Kalman filter(s) for fitting of the tracks together with calibration/transformation
   7   parameters.
   8
   9   Correction/Transformation are currently described by set (TObjArray) of primitive
  10   correction/transformatio - AliTPCTransformation.  Currently we assume that transformation
  11   comute (in first order). AliTPCTransformation describe general non linear transformation.
  12
  13   Current calibration parameters and covariance stored (fCalibParam, fCalibCovar).
  14
  15   Currenly only linear track model implemented.
  16   Fits to be implemented:
  17    0. Plane fitting (Laser CE)
  18    1. Primary vertex fitting.
  19    2. Propagation in magnetic field and fit of planes
  20
  21
  22
  23   How to use it - see  AliTPCkalmanFit::Test function
  24
  25   Simple test (see AliTPCkalmanFit::Test)
  26   AliTPCTransformation::BuildBasicFormulas();
  27   AliTPCkalmanFit *kalmanFit0 = AliTPCkalmanFit::Test(2000);
  28   TFile f("kalmanfitTest.root");
  29
  30   Transformation visualization:
  31   Transforamtion can be visualized using TFn (TF1,TF2 ...) and using tree->Draw()
  32   e.g:
  33   kalmanFit0->SetInstance(kalmanFit0);   //
  34   kalmanFit0->InitTransformation();      //
  35   TF2 fxRZdz("fxRZdz","AliTPCkalmanFit::SGetTPCDeltaXYZ(0,-1,x,0,y)-AliTPCkalmanFit::SGetTPCDeltaXYZ(0,-1,x,0,y-1)",85,245,-250,250);
  36   fxRZdz->Draw("");
  37
  38
  39
  40 */
  41
  42 #include "TRandom.h"
  43 #include "TMath.h"
  44 #include "TBits.h"
  45 #include "TFormula.h"
  46 #include "TF1.h"
  47 #include "TLinearFitter.h"
  48 #include "TFile.h"
  49 #include "THnSparse.h"
  50 #include "TAxis.h"
  51
  52
  53 #include "TTreeStream.h"
  54 #include "AliTrackPointArray.h"
  55 #include "AliLog.h"
  56 #include "AliTPCTransformation.h"
  57 #include "AliTPCkalmanFit.h"
  58
  59 ClassImp(AliTPCkalmanFit)
  60
  61 AliTPCkalmanFit* AliTPCkalmanFit::fgInstance = 0;
  62
  63 AliTPCkalmanFit::AliTPCkalmanFit():
  64   TNamed(),
  65   fCalibration(0),
  66   fCalibParam(0),
  67   fCalibCovar(0),
  68   fLinearParam(0),
  69   fLinearCovar(0),
  70   fLastTimeStamp(-1),
  71   fCurrentAlpha(0),  //! current rotation frame
  72   fCA(0),            //! cosine of current angle
  73   fSA(0)            //! sinus of current angle
  74 {
  75   //
  76   // Default constructor
  77   //
  78   for (Int_t ihis=0; ihis<8; ihis++){
  79     fLinearTrackDelta[ihis]=0;
  80     fLinearTrackPull[ihis]=0;
  81   }
  82 }
  83
  84 void AliTPCkalmanFit::InitTransformation(){
  85   //
  86   // Initialize pointers to the transforamtion functions
  87   //
  88   //
  89   Int_t ncalibs = fCalibration->GetEntries();
  90   for (Int_t icalib=0;icalib<ncalibs; icalib++){
  91     AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(icalib);
  92     transform->Init();
  93   }
  94 }
  95
  96 void AliTPCkalmanFit::Init(){
  97   //
  98   // Initialize parameter vector and covariance matrix
  99   // To be called after initialization of all of the transformations
 100   //
 101   //
 102   Int_t ncalibs = fCalibration->GetEntries();
 103   fCalibParam = new TMatrixD(ncalibs,1);
 104   fCalibCovar = new TMatrixD(ncalibs,ncalibs);
 105   for (Int_t icalib=0;icalib<ncalibs; icalib++){
 106     AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(icalib);
 107     (*fCalibParam)(icalib,0) = transform->fParam;
 108     for (Int_t jcalib=0;jcalib<ncalibs; jcalib++){
 109       if (icalib!=jcalib) (*fCalibCovar)(icalib,jcalib)= 0;
 110       if (icalib==jcalib) (*fCalibCovar)(icalib,jcalib) = transform->fSigma*transform->fSigma;
 111     }
 112   }
 113   //
 114   // Build QA histograms
 115   //
 116   Double_t mpi = TMath::Pi();
 117   //
 118   //                    delta  alpha   y0     z0    ky   kz
 119   Int_t binsQA[6]    = {300,   36*4,   30,    25,   20,  20};
 120   Double_t xminQA[6] = {-0.5,  -mpi,  -120, -250,   -1,  -1.};
 121   Double_t xmaxQA[6] = { 0.5,   mpi,   120,  250,    1,   1.};
 122   TString  axisName[6]={"#Delta",
 123                         "#alpha",
 124                         "y_{0}",
 125                         "z_{0}",
 126                         "tan(#phi)",
 127                         "tan(theta)"};
 128   //
 129   TString deltaName[8]={"#Delta_{y}(cm)",
 130                         "100*#Delta_{#phi}(cm)",
 131                         "100^{2}dy0^{2}/dx0^{2}(cm)",
 132                         "100^{2}dy1^{2}/dx1^{2}(cm)",
 133                         "#Delta_{z}(cm)",
 134                         "100*#Delta_{#theta}(cm)",
 135                         "100^{2}*dz0^{2}/dx0^{2}(cm)",
 136                         "100^{2}*dz1^{2}/dx1^{2}(cm)"};
 137   TString pullName[8]={"#Delta_{y}(unit)",
 138                        "100*#Delta_{#phi}(unit)",
 139                        "100^{2}dy0^{2}/dx0^{2}(unit)",
 140                        "100^{2}dy1^{2}/dx1^{2}(unit)",
 141                        "#Delta_{z}(unit)",
 142                        "100*#Delta_{#theta}(unit)",
 143                        "100^{2}*dz0^{2}/dx0^{2}(unit)",
 144                        "100^{2}*dz1^{2}/dx1^{2}(unit)"};
 145   for (Int_t ihis=0; ihis<8; ihis++){
 146     fLinearTrackDelta[ihis]=0;
 147     fLinearTrackPull[ihis]=0;
 148     xminQA[0]=-0.5; xmaxQA[0] = 0.5;
 149     fLinearTrackDelta[ihis]  = new THnSparseS(deltaName[ihis],deltaName[ihis], 6, binsQA,xminQA, xmaxQA);
 150     xminQA[0]=-10; xmaxQA[0]  = 10;
 151     fLinearTrackPull[ihis]   = new THnSparseS(pullName[ihis],pullName[ihis],   6, binsQA,xminQA, xmaxQA);
 152     for (Int_t iaxis=1; iaxis<6; iaxis++){
 153       fLinearTrackDelta[ihis]->GetAxis(iaxis)->SetName(axisName[iaxis]);
 154       fLinearTrackDelta[ihis]->GetAxis(iaxis)->SetTitle(axisName[iaxis]);
 155       fLinearTrackPull[ihis]->GetAxis(iaxis)->SetName(axisName[iaxis]);
 156       fLinearTrackPull[ihis]->GetAxis(iaxis)->SetTitle(axisName[iaxis]);
 157     }
 158     fLinearTrackDelta[ihis]->GetAxis(0)->SetName(deltaName[ihis]);
 159     fLinearTrackDelta[ihis]->GetAxis(0)->SetTitle(deltaName[ihis]);
 160     fLinearTrackPull[ihis]->GetAxis(0)->SetName(deltaName[ihis]);
 161     fLinearTrackPull[ihis]->GetAxis(0)->SetTitle(deltaName[ihis]);
 162   }
 163
 164
 165 }
 166
 167
 168 void AliTPCkalmanFit::FitTrackLinear(AliTrackPointArray &points, Int_t step, TTreeSRedirector *debug){
 169   //
 170   //
 171
 172   if (!fCalibParam) {
 173     AliError("Kalman Fit not initialized");
 174     return;
 175   }
 176   //
 177   // 1. Make initial track hypothesis
 178   //
 179   TLinearFitter lfitY(2,"pol1");
 180   TLinearFitter lfitZ(2,"pol1");
 181   TVectorD vecZ(2);
 182   TVectorD vecY(2);
 183   //
 184   lfitY.StoreData(kTRUE);
 185   lfitZ.StoreData(kTRUE);
 186   lfitY.ClearPoints();
 187   lfitZ.ClearPoints();
 188   Int_t npoints = points.GetNPoints();
 189   if (npoints<2) return;
 190   const Double_t kFac=npoints*npoints*100;
 191   const Double_t kOff0=0.01*0.01;
 192   const Double_t kOff1=kOff0/(250.*250.);
 193   //
 194   // 0. Make seed
 195   //
 196   //  AliTrackPointArray pointsc(points);
 197   //ApplyCalibration(&pointsc, -1.);
 198   //
 199   // 1.a Choosing reference rotation alpha
 200   //
 201   fCurrentAlpha = TMath::ATan2(points.GetY()[npoints-1]-points.GetY()[0], points.GetX()[npoints-1]-points.GetX()[0]);
 202   fCA = TMath::Cos(fCurrentAlpha);
 203   fSA = TMath::Sin(fCurrentAlpha);
 204   //
 205   // 1.b Fit the track in the rotated frame - MakeSeed
 206   //
 207   for (Int_t ipoint=0; ipoint<npoints-1; ipoint+=step){
 208     Double_t rx =   fCA*points.GetX()[ipoint]+fSA*points.GetY()[ipoint];
 209     Double_t ry =  -fSA*points.GetX()[ipoint]+fCA*points.GetY()[ipoint];
 210     Double_t rz =  points.GetZ()[ipoint];
 211     lfitY.AddPoint(&rx,ry,1);
 212     lfitZ.AddPoint(&rx,rz,1);
 213   }
 214   lfitY.Eval();
 215   lfitZ.Eval();
 216   lfitY.GetParameters(vecY);
 217   lfitZ.GetParameters(vecZ);
 218   //
 219   // 2. Set initial parameters and the covariance matrix
 220   //
 221   Int_t ncalibs = fCalibration->GetEntries();
 222   if (!fLinearParam) {
 223     fLinearParam = new TMatrixD(ncalibs+4,1);
 224     fLinearCovar = new TMatrixD(ncalibs+4,ncalibs+4);
 225   }
 226   //
 227   for (Int_t i=0; i<ncalibs+4;i++){
 228     (*fLinearParam)(i,0)=0;
 229     if (i<ncalibs) (*fLinearParam)(i,0) = (*fCalibParam)(i,0);
 230     for (Int_t j=0; j<ncalibs+4;j++){
 231       (*fLinearCovar)(i,j)=0;
 232       if (i<ncalibs&&j<ncalibs) (*fLinearCovar)(i,j) = (*fCalibCovar)(i,j);
 233     }
 234   }
 235   Double_t chi2Y = lfitY.GetChisquare()/lfitY.GetNpoints();
 236   Double_t chi2Z = lfitZ.GetChisquare()/lfitZ.GetNpoints();
 237   (*fLinearParam)(ncalibs+0,0) =  lfitY.GetParameter(0);
 238   (*fLinearCovar)(ncalibs+0,ncalibs+0)= lfitY.GetCovarianceMatrixElement(0,0)*chi2Y*kFac+kOff0;
 239   (*fLinearParam)(ncalibs+1,0) =  lfitY.GetParameter(1);
 240   (*fLinearCovar)(ncalibs+1,ncalibs+1)= lfitY.GetCovarianceMatrixElement(1,1)*chi2Y*kFac+kOff1;
 241   //
 242   (*fLinearParam)(ncalibs+2,0) =  lfitZ.GetParameter(0);
 243   (*fLinearCovar)(ncalibs+2,ncalibs+2)= lfitZ.GetCovarianceMatrixElement(0,0)*chi2Z*kFac+kOff0;
 244   (*fLinearParam)(ncalibs+3,0) =  lfitZ.GetParameter(1);
 245   (*fLinearCovar)(ncalibs+3,ncalibs+3)= lfitZ.GetCovarianceMatrixElement(1,1)*chi2Z*kFac+kOff1;
 246   //
 247   // Fit thetrack together with correction
 248   //
 249   AliTrackPoint point;
 250   for (Int_t ipoint=0; ipoint<npoints-1; ipoint+=step){
 251     //
 252     if (!points.GetPoint(point,ipoint)) continue;
 253     Double_t erry2 = chi2Y;
 254     Double_t errz2 = chi2Z;
 255     //set error - it is hack
 256     Float_t *cov = (Float_t*) point.GetCov();
 257     cov[1] = erry2+kOff0;
 258     cov[2] = errz2+kOff0;
 259     UpdateLinear(point,debug);
 260     if (!points.GetPoint(point,npoints-1-ipoint)) continue;
 261     //set error - it is hack
 262     cov = (Float_t*) point.GetCov();
 263     cov[1] = erry2+kOff0;
 264     cov[2] = errz2+kOff0;
 265     UpdateLinear(point,debug);
 266   }
 267   //
 268   // save current param and covariance
 269   for (Int_t i=0; i<ncalibs;i++){
 270     AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(i);
 271     transform->fParam= (*fLinearParam)(i,0);
 272     (*fCalibParam)(i,0) = (*fLinearParam)(i,0);
 273     for (Int_t j=0; j<ncalibs;j++){
 274       (*fCalibCovar)(i,j) = (*fLinearCovar)(i,j);
 275     }
 276   }
 277   if (debug){ // dump debug info
 278     (*debug)<<"fitLinear"<<
 279       "vecY.="<<&vecY<<
 280       "vecZ.="<<&vecZ<<
 281       "chi2Y="<<chi2Y<<
 282       "chi2Z="<<chi2Z<<
 283       "lP.="<<fLinearParam<<
 284       "cP.="<<fCalibParam<<
 285       "lC.="<<fLinearCovar<<
 286       "cC.="<<fCalibCovar<<
 287       "\n";
 288   }
 289 }
 290
 291 void AliTPCkalmanFit::DumpTrackLinear(AliTrackPointArray &points, TTreeSRedirector *debug){
 292   //
 293   // Dump the track parameters before and after current calibration
 294   //
 295   // Track is divided to two parts -
 296   // X mean is defined as middle point
 297   //
 298
 299   if (!fCalibParam) {
 300     AliError("Kalman Fit not initialized");
 301     return;
 302   }
 303   //
 304   //
 305   //
 306   TLinearFitter *fitters[16];
 307   TVectorD      *params[16];
 308   TVectorD      *errs[16];
 309   TVectorD      chi2N(16);
 310   Int_t npoints = points.GetNPoints();
 311   AliTrackPointArray pointsTrans(points);
 312   ApplyCalibration(&pointsTrans,-1.);
 313   for (Int_t ifit=0; ifit<8;ifit++){
 314     fitters[ifit]  = new TLinearFitter(2,"pol1");
 315     params[ifit]   = new TVectorD(2);
 316     fitters[ifit+8]= new TLinearFitter(3,"pol2");
 317     params[ifit+8] = new TVectorD(3);
 318     errs[ifit]     = new TVectorD(2);
 319     errs[ifit+8]   = new TVectorD(3);
 320   }
 321   //
 322   // calculate mean x point  and corrdinate frame
 323   //
 324   fCurrentAlpha = TMath::ATan2(points.GetY()[npoints-1]-points.GetY()[0], points.GetX()[npoints-1]-points.GetX()[0]);
 325   fCA = TMath::Cos(fCurrentAlpha);
 326   fSA = TMath::Sin(fCurrentAlpha);
 327   Double_t xmean=0, sum=0;
 328   for (Int_t ipoint=0; ipoint<npoints-1; ipoint++){
 329     Double_t rx  =   fCA*points.GetX()[ipoint]+fSA*points.GetY()[ipoint];
 330     xmean+=rx;
 331     sum++;
 332   }
 333   xmean/=sum;
 334   //
 335   for (Int_t ipoint=0; ipoint<npoints-1; ipoint++){
 336     Double_t rx  =   fCA*points.GetX()[ipoint]+fSA*points.GetY()[ipoint];
 337     Double_t ry  =  -fSA*points.GetX()[ipoint]+fCA*points.GetY()[ipoint];
 338     Double_t rz  =  points.GetZ()[ipoint];
 339     //
 340     Double_t rxT =   fCA*pointsTrans.GetX()[ipoint]+fSA*pointsTrans.GetY()[ipoint];
 341     Double_t ryT =  -fSA*pointsTrans.GetX()[ipoint]+fCA*pointsTrans.GetY()[ipoint];
 342     Double_t rzT =  pointsTrans.GetZ()[ipoint];
 343     if (rx>xmean){
 344       fitters[0]->AddPoint(&rxT, ryT,1);
 345       fitters[2]->AddPoint(&rxT, rzT,1);
 346       fitters[4]->AddPoint(&rx, ry,1);
 347       fitters[6]->AddPoint(&rx, rz,1);
 348       fitters[8]->AddPoint(&rxT, ryT,1);
 349       fitters[10]->AddPoint(&rxT, rzT,1);
 350       fitters[12]->AddPoint(&rx, ry,1);
 351       fitters[14]->AddPoint(&rx, rz,1);
 352     }else{
 353       fitters[1]->AddPoint(&rxT, ryT,1);
 354       fitters[3]->AddPoint(&rxT, rzT,1);
 355       fitters[5]->AddPoint(&rx, ry,1);
 356       fitters[7]->AddPoint(&rx, rz,1);
 357       fitters[9]->AddPoint(&rxT, ryT,1);
 358       fitters[11]->AddPoint(&rxT, rzT,1);
 359       fitters[13]->AddPoint(&rx, ry,1);
 360       fitters[15]->AddPoint(&rx, rz,1);
 361     }
 362   }
 363   for (Int_t ifit=0;ifit<16;ifit++){
 364     fitters[ifit]->Eval();
 365     fitters[ifit]->GetParameters(*(params[ifit]));
 366     fitters[ifit]->GetErrors(*(errs[ifit]));
 367     chi2N[ifit] = TMath::Sqrt(fitters[ifit]->GetChisquare()/(fitters[ifit]->GetNpoints()-1));
 368     (*(errs[ifit]))[0]*=chi2N[ifit];
 369     (*(errs[ifit]))[1]*=chi2N[ifit];
 370     if (ifit>=8) (*(errs[ifit]))[2]*=chi2N[ifit];  // second derivative
 371   }
 372   if (debug){
 373     (*debug)<<"dumpLinear"<<
 374       "alpha="<<fCurrentAlpha<<
 375       "xmean="<<xmean<<
 376       "y0T.="<<params[0]<<
 377       "y1T.="<<params[1]<<
 378       "z0T.="<<params[2]<<
 379       "z1T.="<<params[3]<<
 380       "y0O.="<<params[4]<<
 381       "y1O.="<<params[5]<<
 382       "z0O.="<<params[6]<<
 383       "z1O.="<<params[7]<<
 384       "y0T2.="<<params[8]<<
 385       "y1T2.="<<params[9]<<
 386       "z0T2.="<<params[10]<<
 387       "z1T2.="<<params[11]<<
 388       "y0O2.="<<params[12]<<
 389       "y1O2.="<<params[13]<<
 390       "z0O2.="<<params[14]<<
 391       "z1O2.="<<params[15]<<
 392       "y0TErr.="<<errs[0]<<
 393       "y1TErr.="<<errs[1]<<
 394       "z0TErr.="<<errs[2]<<
 395       "z1TErr.="<<errs[3]<<
 396       "y0OErr.="<<errs[4]<<
 397       "y1OErr.="<<errs[5]<<
 398       "z0OErr.="<<errs[6]<<
 399       "z1OErr.="<<errs[7]<<
 400       "y0T2Err.="<<errs[8]<<
 401       "y1T2Err.="<<errs[9]<<
 402       "z0T2Err.="<<errs[10]<<
 403       "z1T2Err.="<<errs[11]<<
 404       "y0O2Err.="<<errs[12]<<
 405       "y1O2Err.="<<errs[13]<<
 406       "z0O2Err.="<<errs[14]<<
 407       "z1O2Err.="<<errs[15]<<
 408       "chi2.="<<&chi2N<<
 409       "\n";
 410   }
 411   //
 412   //                    delta  alpha   y0     z0    ky   kz
 413   Double_t x[6]={0,0,0,0,0,0};
 414   x[1]=fCurrentAlpha;
 415   x[2]=(*params[0])[0];
 416   x[3]=(*params[2])[0];
 417   x[4]=(*params[0])[1];
 418   x[5]=(*params[2])[1];
 419   //
 420   //Delta y
 421   //
 422   x[0]= (*params[1])[0]-(*params[0])[0];
 423   fLinearTrackDelta[0]->Fill(x);
 424   x[0]/=TMath::Sqrt((*errs[1])[0]*(*errs[1])[0]+(*errs[0])[0]*(*errs[0])[0]);
 425   fLinearTrackPull[0]->Fill(x);
 426   //
 427   //Delta ky
 428   //
 429   x[0]= 100*((*params[1])[1]-(*params[0])[1]);
 430   fLinearTrackDelta[1]->Fill(x);
 431   x[0]/=100*TMath::Sqrt((*errs[1])[1]*(*errs[1])[1]+(*errs[0])[1]*(*errs[0])[1]);
 432   fLinearTrackPull[1]->Fill(x);
 433   //
 434   // ky2_0
 435   //
 436   x[0]= 100*100*((*params[8])[2]);
 437   fLinearTrackDelta[2]->Fill(x);
 438   x[0]/=100*100*TMath::Sqrt((*errs[8])[2]*(*errs[8])[2]);
 439   fLinearTrackPull[2]->Fill(x);
 440   //
 441   // ky2_1
 442   //
 443   x[0]= 100*100*((*params[9])[2]);
 444   fLinearTrackDelta[3]->Fill(x);
 445   x[0]/=100*100*TMath::Sqrt((*errs[9])[2]*(*errs[9])[2]);
 446   fLinearTrackPull[3]->Fill(x);
 447   //
 448   //
 449   //Delta z
 450   //
 451   x[0]= (*params[3])[0]-(*params[2])[0];
 452   fLinearTrackDelta[4]->Fill(x);
 453   x[0]/=TMath::Sqrt((*errs[3])[0]*(*errs[3])[0]+(*errs[2])[0]*(*errs[2])[0]);
 454   fLinearTrackPull[4]->Fill(x);
 455   //
 456   //Delta kz
 457   //
 458   x[0]= 100*((*params[3])[1]-(*params[2])[1]);
 459   fLinearTrackDelta[5]->Fill(x);
 460   x[0]/=100*TMath::Sqrt((*errs[3])[1]*(*errs[3])[1]+(*errs[2])[1]*(*errs[2])[1]);
 461   fLinearTrackPull[5]->Fill(x);
 462   //
 463   // kz2_0
 464   //
 465   x[0]= 100*100*((*params[10])[2]);
 466   fLinearTrackDelta[6]->Fill(x);
 467   x[0]/=100*100*TMath::Sqrt((*errs[10])[2]*(*errs[10])[2]);
 468   fLinearTrackPull[6]->Fill(x);
 469   //
 470   // kz2_1
 471   //
 472   x[0]= 100*100*((*params[11])[2]);
 473   fLinearTrackDelta[7]->Fill(x);
 474   x[0]/=100*100*TMath::Sqrt((*errs[11])[2]*(*errs[11])[2]);
 475   fLinearTrackPull[7]->Fill(x);
 476
 477
 478
 479   for (Int_t ifit=0; ifit<8;ifit++){
 480     delete fitters[ifit];
 481     delete params[ifit];
 482     delete fitters[ifit+8];
 483     delete params[ifit+8];
 484     delete errs[ifit];
 485     delete errs[ifit+8];
 486   }
 487 }
 488
 489
 490 void AliTPCkalmanFit::Propagate(TTreeSRedirector */*debug*/){
 491   //
 492   // Propagate the Kalman
 493   //
 494 }
 495
 496
 497 void AliTPCkalmanFit::PropagateTime(Int_t time){
 498   //
 499   // Propagate the calibration in time
 500   // - Increase covariance matrix
 501   //
 502   if (fCalibCovar) return;
 503   Int_t ncalibs = fCalibration->GetEntries();
 504   Double_t deltaT = 0;
 505   if (fLastTimeStamp>0) deltaT = (fLastTimeStamp-time)/(60*60.);  // delta T in hours
 506   fLastTimeStamp = time;
 507   for (Int_t icalib=0;icalib<ncalibs; icalib++){
 508     AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(icalib);
 509     (*fCalibCovar)(icalib,icalib)+=  transform->fSigma2Time*deltaT;
 510   }
 511 }
 512
 513
 514 void  AliTPCkalmanFit::UpdateLinear(AliTrackPoint &point, TTreeSRedirector *debug){
 515   //
 516   // Update Kalman
 517   //
 518   //
 519   Int_t ncalibs = fCalibration->GetEntries();
 520   Int_t kNmeas  = 2;
 521   Int_t nelem   = ncalibs+4;
 522   TMatrixD &vecXk=*fLinearParam;     // X vector
 523   TMatrixD &covXk=*fLinearCovar;     // X covariance
 524   //
 525   TMatrixD matHk(kNmeas,nelem);     // vector to mesurement
 526   TMatrixD vecYk(kNmeas,1);         // Innovation or measurement residual
 527   TMatrixD vecZk(kNmeas,1);         // Innovation or measurement residual
 528   TMatrixD measR(kNmeas,kNmeas);
 529   TMatrixD matHkT(nelem,kNmeas);    // helper matrix Hk transpose
 530   TMatrixD matSk(kNmeas,kNmeas);    // Innovation (or residual) covariance
 531   TMatrixD matKk(nelem,kNmeas);     // Optimal Kalman gain
 532   TMatrixD covXk2(nelem,nelem);     // helper matrix
 533   TMatrixD covXk3(nelem,nelem);     // helper matrix
 534   TMatrixD mat1(nelem,nelem);
 535   //
 536   // reset matHk
 537   for (Int_t iel=0;iel<nelem;iel++){
 538     for (Int_t ip=0;ip<kNmeas;ip++) {
 539       matHk(ip,iel)=0;
 540     }
 541   }
 542   for (Int_t iel0=0;iel0<nelem;iel0++)
 543     for (Int_t iel1=0;iel1<nelem;iel1++){
 544       if (iel0!=iel1) mat1(iel0,iel1)=0;
 545       if (iel0==iel1) mat1(iel0,iel1)=1;
 546     }
 547   //
 548   // fill matrix Hk
 549   //
 550   Int_t volId = point.GetVolumeID();
 551   Double_t gxyz[3]={point.GetX(), point.GetY(),point.GetZ()};
 552   Double_t rxyz[3]={fCA*gxyz[0]+fSA*gxyz[1],-fSA*gxyz[0]+fCA*gxyz[1] ,point.GetZ()};
 553   //
 554   Double_t dxdydz[3]={0,0,0};
 555   Double_t rdxdydz[3]={0,0,0};
 556
 557   for (Int_t icalib=0;icalib<ncalibs;icalib++){
 558     AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(icalib);
 559     dxdydz[0] = transform->GetDeltaXYZ(0,volId, 1, gxyz[0], gxyz[1], gxyz[2]);
 560     dxdydz[1] = transform->GetDeltaXYZ(1,volId, 1, gxyz[0], gxyz[1], gxyz[2]);
 561     dxdydz[2] = transform->GetDeltaXYZ(2,volId, 1, gxyz[0], gxyz[1], gxyz[2]);
 562     rdxdydz[0] =  fCA*dxdydz[0]+fSA*dxdydz[1];
 563     rdxdydz[1] = -fSA*dxdydz[0]+fCA*dxdydz[1];
 564     rdxdydz[2] =  dxdydz[2];
 565     //
 566     matHk(0,icalib)= rdxdydz[1]-rdxdydz[0]*(*fLinearParam)(ncalibs+1,0);  // shift y + shift x * angleY
 567     matHk(1,icalib)= rdxdydz[2]-rdxdydz[0]*(*fLinearParam)(ncalibs+3,0);  // shift z + shift x * angleZ
 568   }
 569   matHk(0,ncalibs+0)=1;
 570   matHk(0,ncalibs+1)=rxyz[0];
 571   matHk(1,ncalibs+2)=1;
 572   matHk(1,ncalibs+3)=rxyz[0];
 573   //
 574   //
 575   //
 576   vecZk(0,0) =  rxyz[1];
 577   vecZk(1,0) =  rxyz[2];
 578   measR(0,0) = point.GetCov()[1]; measR(0,1)=0;
 579   measR(1,1) = point.GetCov()[2]; measR(1,0)=0;
 580   vecYk = vecZk-matHk*vecXk;               // Innovation or measurement residual
 581   matHkT=matHk.T(); matHk.T();
 582   matSk = (matHk*(covXk*matHkT))+measR;    // Innovation (or residual) covariance
 583   matSk.Invert();
 584   matKk = (covXk*matHkT)*matSk;            //  Optimal Kalman gain
 585   //
 586   covXk2= (mat1-(matKk*matHk));
 587   covXk3 =  covXk2*covXk;
 588   //
 589   CheckCovariance(covXk3,100);
 590   vecXk += matKk*vecYk;                    //  updated vector
 591   covXk = covXk3;
 592   if (debug){ // dump debug info
 593     (*debug)<<"updateLinear"<<
 594       //
 595       "point.="<<&point<<
 596       "vecYk.="<<&vecYk<<
 597       "vecZk.="<<&vecZk<<
 598       "measR.="<<&measR<<
 599       "matHk.="<<&matHk<<
 600       "matHkT.="<<&matHkT<<
 601       "matSk.="<<&matSk<<
 602       "matKk.="<<&matKk<<
 603       "covXk2.="<<&covXk2<<
 604       "covXk.="<<&covXk<<
 605       "vecXk.="<<&vecXk<<
 606       "\n";
 607   }
 608 }
 609
 610
 611 AliTrackPointArray * AliTPCkalmanFit::MakePointArrayLinear(Double_t alpha, Double_t y0, Double_t z0, Double_t ky, Double_t kz, Double_t err){
 612   //
 613   //
 614   //
 615   AliTrackPointArray array(500);
 616   Float_t cov[10];  // dummy covariance
 617   Int_t npoints=0;
 618   for (Int_t i=0;i<6;i++) cov[i]=0.001;
 619   for (Int_t i=0;i<500;i++){
 620     AliTrackPoint point(0, 0, 0, cov, -1,0,0);
 621     array.AddPoint(npoints, &point);
 622     npoints++;
 623   }
 624   npoints=0;
 625   for (Float_t ir = -245; ir<245; ir++){
 626     //
 627     //
 628     if (TMath::Abs(ir)<80) continue;
 629     Double_t ca = TMath::Cos(alpha);
 630     Double_t sa = TMath::Sin(alpha);
 631     Double_t lx = ir;
 632     Double_t ly = y0+lx*ky+gRandom->Gaus(0,err);
 633     Double_t lz = z0+lx*kz+gRandom->Gaus(0,err);
 634     Double_t gx = lx*ca-ly*sa;
 635     Double_t gy = lx*sa+ly*ca;
 636     Double_t gz = lz;
 637     Double_t galpha = TMath::ATan2(gy,gx);
 638     Int_t isec = TMath::Nint((9*galpha/TMath::Pi()+0.5));
 639     if (isec<0) isec+=18;
 640     if (gz<0) isec+=18;
 641     if (ir>130) isec+=36;
 642     //
 643     AliTrackPoint point(gx, gy, gz, cov, isec,0,0);
 644     array.AddPoint(npoints, &point);
 645     npoints++;
 646   }
 647   AliTrackPointArray *parray = new AliTrackPointArray(npoints);
 648   AliTrackPoint point;
 649   for (Int_t i=0;i<npoints;i++){
 650     array.GetPoint(point,i);
 651     parray->AddPoint(i,&point);
 652   }
 653
 654   return parray;
 655 }
 656
 657 void AliTPCkalmanFit::AddCalibration(AliTPCTransformation * calib){
 658   //
 659   // Add calibration
 660   //
 661   if (!fCalibration) fCalibration = new TObjArray(2000);
 662   fCalibration->AddLast(calib);
 663 }
 664
 665 void AliTPCkalmanFit::ApplyCalibration(AliTrackPointArray *array, Double_t csign){
 666   //
 667   //
 668   //
 669   if (!fCalibration) return;
 670   Int_t ncalibs = fCalibration->GetEntries();
 671   if (ncalibs==0) return;
 672   Int_t npoints = array->GetNPoints();
 673   for (Int_t ipoint=0; ipoint<npoints; ipoint++){
 674     Int_t volId = array->GetVolumeID()[ipoint];
 675     Double_t xyz[3]={array->GetX()[ipoint], array->GetY()[ipoint],array->GetZ()[ipoint]};
 676     Double_t dxdydz[3]={0,0,0};
 677     for (Int_t icalib=0; icalib<ncalibs; icalib++){
 678       AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(icalib);
 679       dxdydz[0] += transform->GetDeltaXYZ(0,volId, transform->fParam, xyz[0], xyz[1], xyz[2]);
 680       dxdydz[1] += transform->GetDeltaXYZ(1,volId, transform->fParam, xyz[0], xyz[1], xyz[2]);
 681       dxdydz[2] += transform->GetDeltaXYZ(2,volId, transform->fParam, xyz[0], xyz[1], xyz[2]);
 682     }
 683     ((Float_t*)array->GetX())[ipoint]+=csign*dxdydz[0];
 684     ((Float_t*)array->GetY())[ipoint]+=csign*dxdydz[1];
 685     ((Float_t*)array->GetZ())[ipoint]+=csign*dxdydz[2];
 686   }
 687 }
 688
 689 Bool_t AliTPCkalmanFit::DumpCorelation(Double_t threshold){
 690   //
 691   //
 692   //
 693   TMatrixD &mat = *fCalibCovar;
 694   Int_t nrow= mat.GetNrows();
 695   for (Int_t irow=0; irow<nrow; irow++){
 696     for (Int_t icol=irow+1; icol<nrow; icol++){
 697       Double_t diag = TMath::Sqrt(TMath::Abs(mat(irow,irow)*mat(icol,icol)));
 698       if (diag<=0){
 699         printf("Negative covariance\t%d\t%d\t%f\n",irow,icol, mat(irow,icol));
 700         continue;
 701       }
 702       Double_t corr0 = mat(irow,icol)/diag;
 703       if (TMath::Abs(corr0)>threshold){
 704         AliTPCTransformation * trans0 = GetTransformation(irow);
 705         AliTPCTransformation * trans1 = GetTransformation(icol);
 706         printf("%d\t%d\t%s\t%s\t%f\t%f\t%f\n", irow,icol, trans0->GetName(), trans1->GetName(),
 707                TMath::Sqrt(mat(irow,irow)), TMath::Sqrt(mat(icol,icol)), corr0);
 708       }
 709     }
 710   }
 711   return (nrow>0);
 712 }
 713
 714 Bool_t AliTPCkalmanFit::DumpCalib(){
 715   //
 716   // Print calibration entries - name, value, error
 717   //
 718   TMatrixD &mat = *fCalibCovar;
 719   Int_t nrow= mat.GetNrows();
 720   for (Int_t irow=0; irow<nrow; irow++){
 721     AliTPCTransformation * trans0 = GetTransformation(irow);
 722     printf("%d\t%s\t%f\t%f\n",
 723            irow,
 724            trans0->GetName(),
 725            (*fCalibParam)(irow,0),
 726            TMath::Sqrt(mat(irow,irow)));
 727   }
 728   return (nrow>0);
 729 }
 730
 731
 732 Bool_t  AliTPCkalmanFit::CheckCovariance(TMatrixD &mat, Float_t /*maxEl*/){
 733   //
 734   // Check consistency of covariance matrix
 735   // + symetrize coavariance matrix
 736   Bool_t isOK=kTRUE;
 737   Int_t nrow= mat.GetNrows();
 738   for (Int_t irow=0; irow<nrow; irow++){
 739     if (mat(irow,irow)<=0){
 740       printf("Negative covariance\t%d\t%f\n",irow,mat(irow,irow));
 741       isOK=kFALSE;
 742     }
 743 //     if (mat(irow,irow)>maxEl){
 744 //       printf("Too big  covariance\t%d\t%f\n",irow,mat(irow,irow));
 745 //       isOK=kFALSE;
 746 //     }
 747
 748     for (Int_t icol=0; icol<nrow; icol++){
 749       //      if (mat(irow,irow)
 750       Double_t diag = TMath::Sqrt(TMath::Abs(mat(irow,irow)*mat(icol,icol)));
 751       if (diag<=0){
 752         printf("Negative covariance\t%d\t%d\t%f\n",irow,icol, mat(irow,icol));
 753         isOK=kFALSE;
 754         continue;
 755       }
 756       Double_t cov0 = mat(irow,icol)/diag;
 757       Double_t cov1 = mat(icol,irow)/diag;
 758       if (TMath::Abs(cov0)>1 || TMath::Abs(cov1)>1 ){
 759         printf("Covariance Problem %d\t%d\t%f\t%f\n",irow,icol, cov0, cov1);
 760         isOK=kFALSE;
 761       }
 762       if (TMath::Abs(cov0-cov1)>0.0000001){
 763         printf("Asymetry problem %d\t%d\t%f\t%f\n",irow,icol, cov0, cov1);
 764         isOK=kFALSE;
 765       }
 766       //
 767       // symetrize the covariance matrix
 768       Double_t mean = (mat(irow,icol)+ mat(icol,irow))*0.5;
 769       mat(irow,icol)=mean;
 770       mat(icol,irow)=mean;
 771     }
 772   }
 773   return isOK;
 774 }
 775
 776
 777 AliTPCkalmanFit *  AliTPCkalmanFit::Test(Int_t ntracks){
 778   //
 779   // This is test example
 780   //
 781
 782   //
 783   // 1. Setup transformation
 784   //
 785
 786
 787   TVectorD fpar(10);
 788   AliTPCTransformation * transformation=0;
 789   AliTPCkalmanFit * kalmanFit0 =  new AliTPCkalmanFit;
 790   AliTPCkalmanFit * kalmanFit2 =  new AliTPCkalmanFit;
 791   //
 792   // Radial scaling
 793   //
 794   for (Int_t iside=0; iside<=1; iside++)
 795     for (Int_t ipar0=0; ipar0<3; ipar0++)
 796       for (Int_t ipar1=0; ipar1<3; ipar1++){
 797         fpar[0]=ipar0;
 798         fpar[1]=ipar1;
 799         if (ipar0+ipar1==0) continue;
 800         Double_t param = (gRandom->Rndm()-0.5)*0.5;  // generate random parameters
 801         char tname[100];
 802         sprintf(tname,"tscalingR%d%dSide%d",ipar0,ipar1,iside);
 803         transformation = new AliTPCTransformation(tname,AliTPCTransformation::BitsSide(iside),"TPCscalingRPol",0,0,1);
 804         transformation->SetParams(0,5*0.25,0,&fpar);
 805         kalmanFit0->AddCalibration(transformation);
 806         transformation = new AliTPCTransformation(tname,AliTPCTransformation::BitsSide(iside),"TPCscalingRPol",0,0,1);
 807         transformation->SetParams(param,0.25,0,&fpar);
 808         kalmanFit2->AddCalibration(transformation);
 809       }
 810
 811
 812   //
 813   // 2. Init transformation
 814   //
 815   kalmanFit2->Init();
 816   kalmanFit0->Init();
 817
 818   //
 819   // 3. Run kalman filter
 820   //
 821   Int_t ncalibs = kalmanFit0->fCalibration->GetEntries();
 822   TVectorD err(ncalibs);
 823   TTreeSRedirector *pcstream = new TTreeSRedirector("kalmanfitTest.root");
 824   for (Int_t i=0;i<ntracks;i++){
 825     if (i%100==0) printf("%d\n",i);
 826      Double_t alpha = gRandom->Rndm()*TMath::TwoPi();
 827      Double_t   y0  = (gRandom->Rndm()-0.5)*180;
 828      Double_t   z0  = (gRandom->Rndm()-0.5)*250*2;
 829      Double_t   ky  = (gRandom->Rndm()-0.5)*1;
 830      Double_t   kz  = (gRandom->Rndm()-0.5)*1;
 831      //generate random TPC track
 832      AliTrackPointArray * array  = AliTPCkalmanFit::MakePointArrayLinear(alpha,y0,z0, ky, kz,0.04);
 833      AliTrackPointArray * arrayB = new AliTrackPointArray(*array);  // backup
 834      kalmanFit2->ApplyCalibration(array,1.);  // misalign ideal track
 835      for (Int_t icalib=0; icalib<ncalibs; icalib++){
 836        err[icalib] = TMath::Sqrt((*kalmanFit0->fCalibCovar)(icalib,icalib));
 837      }
 838      //
 839      (*pcstream)<<"dump0"<<
 840        "alpha="<<alpha<<
 841        "y0="<<y0<<
 842        "z0="<<z0<<
 843        "ky="<<ky<<
 844        "lz="<<kz<<
 845        "p.="<<array<<
 846        "pB.="<<arrayB<<
 847        "cparam.="<<kalmanFit0->fCalibParam<<
 848        "ccovar.="<<kalmanFit0->fCalibCovar<<
 849        "err.="<<&err<<
 850        "gparam.="<<kalmanFit2->fCalibParam<<
 851        "gcovar.="<<kalmanFit2->fCalibCovar<<
 852        "\n";
 853      if (i%20==0) {
 854        kalmanFit0->FitTrackLinear(*array,1,pcstream); // fit track - dump intermediate results
 855      }else{
 856        kalmanFit0->FitTrackLinear(*array,1,0);        // fit track + calibration
 857      }
 858      kalmanFit0->DumpTrackLinear(*array,pcstream);    // dump track residuals to the tree + fill histograms
 859   }
 860   pcstream->GetFile()->cd();
 861   kalmanFit0->Write("kalmanFit");
 862   delete pcstream;
 863   return kalmanFit0;
 864 }
 865
 866
 867 Double_t AliTPCkalmanFit::GetTPCDeltaXYZ(Int_t coord, Int_t volID, Double_t x, Double_t y, Double_t z){
 868   //
 869   // function for visualization purposes
 870   //
 871   if (!fCalibration) return 0;
 872   Int_t ncalibs = fCalibration->GetEntries();
 873   if (ncalibs==0) return 0;
 874   Double_t dxdydz[3]={0,0,0};
 875   //
 876   if (volID<0){
 877     Double_t alpha       = TMath::ATan2(y,x);
 878     Double_t r           = TMath::Sqrt(y*y+x*x);
 879     volID                = TMath::Nint(9*alpha/TMath::Pi()-0.5);
 880     if (volID<0) volID+=18;
 881     if (z<0) volID+=18;
 882     if (r>120) volID+=36;
 883   }
 884   for (Int_t icalib=0; icalib<ncalibs; icalib++){
 885     AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(icalib);
 886     Double_t param = (*fCalibParam)(icalib,0);
 887     dxdydz[coord] += transform->GetDeltaXYZ(coord,volID, param, x, y,z);
 888   }
 889   return dxdydz[coord];
 890 }
 891
 892 Double_t AliTPCkalmanFit::SGetTPCDeltaXYZ(Int_t coord, Int_t volID, Double_t x, Double_t y, Double_t z){
 893   //
 894   //
 895   //
 896   if (AliTPCkalmanFit::fgInstance==0) return 0;
 897   return AliTPCkalmanFit::fgInstance->GetTPCDeltaXYZ(coord, volID,x,y,z);
 898 }
 899