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   //
  36   TF2 fxRZdz("fxRZdz","sign(y)*1000*(AliTPCkalmanFit::SGetTPCDeltaXYZ(0,-1,0,x,0,y)-AliTPCkalmanFit::SGetTPCDeltaXYZ(0,-1,0,x,0,y-1))",85,245,-250,250);
  37   fxRZdz->Draw("");
  38
  39   TF2 fxRZ("fxRZ","sign(y)*10*(AliTPCkalmanFit::SGetTPCDeltaXYZ(0,-1,x,0,y))",85,245,-250,250);
  40   fxRZ->Draw("");
  41
  42
  43
  44 */
  45
  46 #include "TRandom.h"
  47 #include "TMath.h"
  48 #include "TBits.h"
  49 #include "TFormula.h"
  50 #include "TF1.h"
  51 #include "TLinearFitter.h"
  52 #include "TFile.h"
  53 #include "THnSparse.h"
  54 #include "TAxis.h"
  55
  56
  57 #include "TTreeStream.h"
  58 #include "AliTrackPointArray.h"
  59 #include "AliLog.h"
  60 #include "AliTPCTransformation.h"
  61 #include "AliTPCkalmanFit.h"
  62
  63 ClassImp(AliTPCkalmanFit)
  64
  65 AliTPCkalmanFit* AliTPCkalmanFit::fgInstance = 0;
  66
  67 AliTPCkalmanFit::AliTPCkalmanFit():
  68   TNamed(),
  69   fCalibration(0),
  70   fCalibParam(0),
  71   fCalibCovar(0),
  72   fLinearParam(0),
  73   fLinearCovar(0),
  74   fLastTimeStamp(-1),
  75   fCurrentAlpha(0),  //! current rotation frame
  76   fCA(0),            //! cosine of current angle
  77   fSA(0)            //! sinus of current angle
  78 {
  79   //
  80   // Default constructor
  81   //
  82   for (Int_t ihis=0; ihis<12; ihis++){
  83     fLinearTrackDelta[ihis]=0;
  84     fLinearTrackPull[ihis]=0;
  85   }
  86 }
  87
  88 void AliTPCkalmanFit::InitTransformation(){
  89   //
  90   // Initialize pointers to the transforamtion functions
  91   //
  92   //
  93   Int_t ncalibs = fCalibration->GetEntries();
  94   for (Int_t icalib=0;icalib<ncalibs; icalib++){
  95     AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(icalib);
  96     transform->Init();
  97   }
  98 }
  99
 100 void AliTPCkalmanFit::Add(const AliTPCkalmanFit * kalman){
 101   //
 102   //
 103   //
 104   Update(kalman);
 105   for (Int_t i=0;i<12;i++){
 106     if (fLinearTrackDelta[i] && kalman->fLinearTrackDelta[i]){
 107       fLinearTrackDelta[i]->Add(kalman->fLinearTrackDelta[i]);
 108     }
 109     if (fLinearTrackPull[i] && kalman->fLinearTrackPull[i]){
 110       fLinearTrackPull[i]->Add(kalman->fLinearTrackPull[i]);
 111     }
 112   }
 113
 114 }
 115
 116
 117 void AliTPCkalmanFit::Init(){
 118   //
 119   // Initialize parameter vector and covariance matrix
 120   // To be called after initialization of all of the transformations
 121   //
 122   //
 123   Int_t ncalibs = fCalibration->GetEntries();
 124   fCalibParam = new TMatrixD(ncalibs,1);
 125   fCalibCovar = new TMatrixD(ncalibs,ncalibs);
 126   for (Int_t icalib=0;icalib<ncalibs; icalib++){
 127     AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(icalib);
 128     (*fCalibParam)(icalib,0) = transform->GetParam();
 129     for (Int_t jcalib=0;jcalib<ncalibs; jcalib++){
 130       if (icalib!=jcalib) (*fCalibCovar)(icalib,jcalib)= 0;
 131       if (icalib==jcalib) (*fCalibCovar)(icalib,jcalib) = transform->GetSigma()*transform->GetSigma();
 132     }
 133   }
 134   //
 135   // Build QA histograms
 136   //
 137   Double_t mpi = TMath::Pi();
 138   //
 139   //                    delta  alpha   y0     z0    ky   kz
 140   Int_t binsQA[6]    = {300,   36*4,   30,    25,   20,  20};
 141   Double_t xminQA[6] = {-0.5,  -mpi,  -120, -250,   -1,  -1.};
 142   Double_t xmaxQA[6] = { 0.5,   mpi,   120,  250,    1,   1.};
 143   TString  axisName[6]={"#Delta",
 144                         "#alpha",
 145                         "y_{0}",
 146                         "z_{0}",
 147                         "tan(#phi)",
 148                         "tan(theta)"};
 149   //
 150   TString deltaName[12]={"#Delta_{y}(cm)",
 151                          "100*#Delta_{#phi}(cm)",
 152                          "100^{2}dy0^{2}/dx0^{2}(cm)",
 153                          "100^{2}dy1^{2}/dx1^{2}(cm)",
 154                          "#Delta_{z}(cm)",
 155                          "100*#Delta_{#theta}(cm)",
 156                          "100^{2}*dz0^{2}/dx0^{2}(cm)",
 157                          "100^{2}*dz1^{2}/dx1^{2}(cm)",
 158                          "RMSy_{0} (cm)",
 159                          "RMSy_{1} (cm)",
 160                          "RMSz_{0} (cm)",
 161                          "RMSz_{1} (cm)"};
 162
 163   TString pullName[12]={"#Delta_{y}(unit)",
 164                         "100*#Delta_{#phi}(unit)",
 165                         "100^{2}dy0^{2}/dx0^{2}(unit)",
 166                         "100^{2}dy1^{2}/dx1^{2}(unit)",
 167                         "#Delta_{z}(unit)",
 168                         "100*#Delta_{#theta}(unit)",
 169                         "100^{2}*dz0^{2}/dx0^{2}(unit)",
 170                         "100^{2}*dz1^{2}/dx1^{2}(unit)"
 171                         "RMSy_{0} (cm)",
 172                         "RMSy_{1} (cm)",
 173                         "RMSz_{0} (cm)",
 174                         "RMSz_{1} (cm)"};
 175   //
 176   //
 177   //
 178   for (Int_t ihis=0; ihis<12; ihis++){
 179     fLinearTrackDelta[ihis]=0;
 180     fLinearTrackPull[ihis]=0;
 181     xminQA[0]=-0.5; xmaxQA[0] = 0.5;
 182     fLinearTrackDelta[ihis]  = new THnSparseS(deltaName[ihis],deltaName[ihis], 6, binsQA,xminQA, xmaxQA);
 183     xminQA[0]=-10; xmaxQA[0]  = 10;
 184     fLinearTrackPull[ihis]   = new THnSparseS(pullName[ihis],pullName[ihis],   6, binsQA,xminQA, xmaxQA);
 185     for (Int_t iaxis=1; iaxis<6; iaxis++){
 186       fLinearTrackDelta[ihis]->GetAxis(iaxis)->SetName(axisName[iaxis]);
 187       fLinearTrackDelta[ihis]->GetAxis(iaxis)->SetTitle(axisName[iaxis]);
 188       fLinearTrackPull[ihis]->GetAxis(iaxis)->SetName(axisName[iaxis]);
 189       fLinearTrackPull[ihis]->GetAxis(iaxis)->SetTitle(axisName[iaxis]);
 190     }
 191     fLinearTrackDelta[ihis]->GetAxis(0)->SetName(deltaName[ihis]);
 192     fLinearTrackDelta[ihis]->GetAxis(0)->SetTitle(deltaName[ihis]);
 193     fLinearTrackPull[ihis]->GetAxis(0)->SetName(deltaName[ihis]);
 194     fLinearTrackPull[ihis]->GetAxis(0)->SetTitle(deltaName[ihis]);
 195   }
 196
 197
 198 }
 199
 200 void AliTPCkalmanFit::SetStatus(const char * mask, Bool_t setOn, Bool_t isOr){
 201   //
 202   // 0. To activate all transforamtion call SetStatus(0,kTRUE)
 203   // 1. To disable everything               SetStatus(0,kFALSE)
 204   // 2. To activate/desactivate             SetStatus("xxx",kTRUE/kFALSE,kFALSE)
 205   Int_t ncalibs = fCalibration->GetEntries();
 206   if (mask==0) {
 207     for (Int_t i=0; i<ncalibs;i++){
 208       AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(i);
 209       transform->SetActive(setOn);
 210     }
 211   }
 212   else{
 213     for (Int_t i=0; i<ncalibs;i++){
 214       AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(i);
 215       TString  strName(transform->GetName());
 216       if (strName.Contains(mask)){
 217         if (!isOr) transform->SetActive( transform->IsActive() && setOn);
 218         if (isOr)  transform->SetActive( transform->IsActive() || setOn);
 219       }
 220     }
 221   }
 222 }
 223
 224
 225 void AliTPCkalmanFit::Update(const AliTPCkalmanFit * kalman){
 226   //
 227   // Update Kalman filter
 228   //
 229   Int_t ncalibs = fCalibration->GetEntries();
 230   TMatrixD vecXk=*fCalibParam;       // X vector
 231   TMatrixD covXk=*fCalibCovar;       // X covariance
 232   TMatrixD &vecZk = *(kalman->fCalibParam);
 233   TMatrixD &measR = *(kalman->fCalibCovar);
 234
 235   TMatrixD matHk(ncalibs,ncalibs);   // vector to mesurement
 236   TMatrixD vecYk(ncalibs,1);         // Innovation or measurement residual
 237   TMatrixD matHkT(ncalibs,ncalibs);  // helper matrix Hk transpose
 238   TMatrixD matSk(ncalibs,ncalibs);   // Innovation (or residual) covariance
 239   TMatrixD matKk(ncalibs,ncalibs);   // Optimal Kalman gain
 240   TMatrixD covXk2(ncalibs,ncalibs);  // helper matrix
 241   TMatrixD covXk3(ncalibs,ncalibs);  // helper matrix
 242   //
 243   for (Int_t i=0;i<ncalibs;i++){
 244     for (Int_t j=0;j<ncalibs;j++) matHk(i,j)=0;
 245     matHk(i,i)=1;
 246   }
 247    vecYk = vecZk-matHk*vecXk;               // Innovation or measurement residual
 248    matHkT=matHk.T(); matHk.T();
 249    matSk = (matHk*(covXk*matHkT))+measR;    // Innovation (or residual) covariance
 250    matSk.Invert();
 251    matKk = (covXk*matHkT)*matSk;            //  Optimal Kalman gain
 252    vecXk += matKk*vecYk;                    //  updated vector
 253    covXk2= (matHk-(matKk*matHk));
 254    covXk3 =  covXk2*covXk;
 255    covXk = covXk3;
 256    (*fCalibParam) = vecXk;
 257    (*fCalibCovar) = covXk;
 258 }
 259
 260
 261
 262
 263
 264 void AliTPCkalmanFit::FitTrackLinear(AliTrackPointArray &points, TTreeSRedirector *debug,  Float_t scalingRMSY, Float_t scalingRMSZ){
 265   //
 266   //
 267
 268   if (!fCalibParam) {
 269     AliError("Kalman Fit not initialized");
 270     return;
 271   }
 272   //
 273   // 1. Make initial track hypothesis
 274   //
 275   TLinearFitter lfitY(2,"pol1");
 276   TLinearFitter lfitZ(2,"pol1");
 277   TVectorD vecZ(2);
 278   TVectorD vecY(2);
 279   //
 280   lfitY.StoreData(kTRUE);
 281   lfitZ.StoreData(kTRUE);
 282   lfitY.ClearPoints();
 283   lfitZ.ClearPoints();
 284   Int_t npoints = points.GetNPoints();
 285   if (npoints<2) return;
 286   const Double_t kFac=npoints*npoints*100;
 287   const Double_t kOff0=0.01*0.01;
 288   const Double_t kOff1=kOff0/(250.*250.);
 289   //
 290   // 0. Make seed
 291   //
 292   //  AliTrackPointArray pointsc(points);
 293   //ApplyCalibration(&pointsc, -1.);
 294   //
 295   // 1.a Choosing reference rotation alpha
 296   //
 297   fCurrentAlpha = TMath::ATan2(points.GetY()[npoints-1]-points.GetY()[0], points.GetX()[npoints-1]-points.GetX()[0]);
 298   fCA = TMath::Cos(fCurrentAlpha);
 299   fSA = TMath::Sin(fCurrentAlpha);
 300   //
 301   // 1.b Fit the track in the rotated frame - MakeSeed
 302   //
 303   for (Int_t ipoint=0; ipoint<npoints-1; ipoint++){
 304     Double_t rx =   fCA*points.GetX()[ipoint]+fSA*points.GetY()[ipoint];
 305     Double_t ry =  -fSA*points.GetX()[ipoint]+fCA*points.GetY()[ipoint];
 306     Double_t rz =  points.GetZ()[ipoint];
 307     lfitY.AddPoint(&rx,ry,1);
 308     lfitZ.AddPoint(&rx,rz,1);
 309   }
 310   lfitY.Eval();
 311   lfitZ.Eval();
 312   lfitY.GetParameters(vecY);
 313   lfitZ.GetParameters(vecZ);
 314   //
 315   // 2. Set initial parameters and the covariance matrix
 316   //
 317   Int_t ncalibs = fCalibration->GetEntries();
 318   if (!fLinearParam) {
 319     fLinearParam = new TMatrixD(ncalibs+4,1);
 320     fLinearCovar = new TMatrixD(ncalibs+4,ncalibs+4);
 321   }
 322   //
 323   for (Int_t i=0; i<ncalibs+4;i++){
 324     (*fLinearParam)(i,0)=0;
 325     if (i<ncalibs) (*fLinearParam)(i,0) = (*fCalibParam)(i,0);
 326     for (Int_t j=0; j<ncalibs+4;j++){
 327       (*fLinearCovar)(i,j)=0;
 328       if (i<ncalibs&&j<ncalibs) (*fLinearCovar)(i,j) = (*fCalibCovar)(i,j);
 329     }
 330   }
 331   Double_t chi2Y = lfitY.GetChisquare()/lfitY.GetNpoints();
 332   Double_t chi2Z = lfitZ.GetChisquare()/lfitZ.GetNpoints();
 333   (*fLinearParam)(ncalibs+0,0) =  lfitY.GetParameter(0);
 334   (*fLinearCovar)(ncalibs+0,ncalibs+0)= lfitY.GetCovarianceMatrixElement(0,0)*chi2Y*kFac+kOff0;
 335   (*fLinearParam)(ncalibs+1,0) =  lfitY.GetParameter(1);
 336   (*fLinearCovar)(ncalibs+1,ncalibs+1)= lfitY.GetCovarianceMatrixElement(1,1)*chi2Y*kFac+kOff1;
 337   //
 338   (*fLinearParam)(ncalibs+2,0) =  lfitZ.GetParameter(0);
 339   (*fLinearCovar)(ncalibs+2,ncalibs+2)= lfitZ.GetCovarianceMatrixElement(0,0)*chi2Z*kFac+kOff0;
 340   (*fLinearParam)(ncalibs+3,0) =  lfitZ.GetParameter(1);
 341   (*fLinearCovar)(ncalibs+3,ncalibs+3)= lfitZ.GetCovarianceMatrixElement(1,1)*chi2Z*kFac+kOff1;
 342   //
 343   // Fit thetrack together with correction
 344   //
 345   AliTrackPoint point;
 346   for (Int_t ipoint=0; ipoint<npoints-1; ipoint++){
 347     //
 348     if (!points.GetPoint(point,ipoint)) continue;
 349     Double_t erry2 = chi2Y;
 350     Double_t errz2 = chi2Z;
 351     //set error - it is hack
 352     Float_t *cov = (Float_t*) point.GetCov();
 353     cov[1] = (erry2+kOff0)*scalingRMSY;
 354     cov[2] = (errz2+kOff0)*scalingRMSZ;
 355     UpdateLinear(point,debug);
 356     if (!points.GetPoint(point,npoints-1-ipoint)) continue;
 357     //set error - it is hack
 358     cov = (Float_t*) point.GetCov();
 359     cov[1] = (erry2+kOff0)*scalingRMSY;
 360     cov[2] = (errz2+kOff0)*scalingRMSZ;
 361     UpdateLinear(point,debug);
 362   }
 363   //
 364   // save current param and covariance
 365   for (Int_t i=0; i<ncalibs;i++){
 366     AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(i);
 367     transform->SetParam( (*fLinearParam)(i,0));
 368     (*fCalibParam)(i,0) = (*fLinearParam)(i,0);
 369     for (Int_t j=0; j<ncalibs;j++){
 370       (*fCalibCovar)(i,j) = (*fLinearCovar)(i,j);
 371     }
 372   }
 373   if (debug){ // dump debug info
 374     (*debug)<<"fitLinear"<<
 375       "vecY.="<<&vecY<<
 376       "vecZ.="<<&vecZ<<
 377       "chi2Y="<<chi2Y<<
 378       "chi2Z="<<chi2Z<<
 379       "lP.="<<fLinearParam<<
 380       "cP.="<<fCalibParam<<
 381       "lC.="<<fLinearCovar<<
 382       "cC.="<<fCalibCovar<<
 383       "\n";
 384   }
 385 }
 386
 387 void AliTPCkalmanFit::DumpTrackLinear(AliTrackPointArray &points, TTreeSRedirector *debug){
 388   //
 389   // Dump the track parameters before and after current calibration
 390   //
 391   // Track is divided to two parts -
 392   // X mean is defined as middle point
 393   //
 394
 395   if (!fCalibParam) {
 396     AliError("Kalman Fit not initialized");
 397     return;
 398   }
 399   //
 400   //
 401   //
 402   TLinearFitter *fitters[16];
 403   TVectorD      *params[16];
 404   TVectorD      *errs[16];
 405   TVectorD      chi2N(16);
 406   Int_t npoints = points.GetNPoints();
 407   AliTrackPointArray pointsTrans(points);
 408   ApplyCalibration(&pointsTrans,-1.);
 409   for (Int_t ifit=0; ifit<8;ifit++){
 410     fitters[ifit]  = new TLinearFitter(2,"pol1");
 411     params[ifit]   = new TVectorD(2);
 412     fitters[ifit+8]= new TLinearFitter(3,"pol2");
 413     params[ifit+8] = new TVectorD(3);
 414     errs[ifit]     = new TVectorD(2);
 415     errs[ifit+8]   = new TVectorD(3);
 416   }
 417   //
 418   // calculate mean x point  and corrdinate frame
 419   //
 420   fCurrentAlpha = TMath::ATan2(points.GetY()[npoints-1]-points.GetY()[0], points.GetX()[npoints-1]-points.GetX()[0]);
 421   fCA = TMath::Cos(fCurrentAlpha);
 422   fSA = TMath::Sin(fCurrentAlpha);
 423   Double_t xmean=0, sum=0;
 424   for (Int_t ipoint=0; ipoint<npoints-1; ipoint++){
 425     Double_t rx  =   fCA*points.GetX()[ipoint]+fSA*points.GetY()[ipoint];
 426     xmean+=rx;
 427     sum++;
 428   }
 429   xmean/=sum;
 430   //
 431   for (Int_t ipoint=0; ipoint<npoints-1; ipoint++){
 432     Double_t rx  =   fCA*points.GetX()[ipoint]+fSA*points.GetY()[ipoint];
 433     Double_t ry  =  -fSA*points.GetX()[ipoint]+fCA*points.GetY()[ipoint];
 434     Double_t rz  =  points.GetZ()[ipoint];
 435     //
 436     Double_t rxT =   fCA*pointsTrans.GetX()[ipoint]+fSA*pointsTrans.GetY()[ipoint];
 437     Double_t ryT =  -fSA*pointsTrans.GetX()[ipoint]+fCA*pointsTrans.GetY()[ipoint];
 438     Double_t rzT =  pointsTrans.GetZ()[ipoint];
 439     if (rx>xmean){
 440       fitters[0]->AddPoint(&rxT, ryT,1);
 441       fitters[2]->AddPoint(&rxT, rzT,1);
 442       fitters[4]->AddPoint(&rx, ry,1);
 443       fitters[6]->AddPoint(&rx, rz,1);
 444       fitters[8]->AddPoint(&rxT, ryT,1);
 445       fitters[10]->AddPoint(&rxT, rzT,1);
 446       fitters[12]->AddPoint(&rx, ry,1);
 447       fitters[14]->AddPoint(&rx, rz,1);
 448     }else{
 449       fitters[1]->AddPoint(&rxT, ryT,1);
 450       fitters[3]->AddPoint(&rxT, rzT,1);
 451       fitters[5]->AddPoint(&rx, ry,1);
 452       fitters[7]->AddPoint(&rx, rz,1);
 453       fitters[9]->AddPoint(&rxT, ryT,1);
 454       fitters[11]->AddPoint(&rxT, rzT,1);
 455       fitters[13]->AddPoint(&rx, ry,1);
 456       fitters[15]->AddPoint(&rx, rz,1);
 457     }
 458   }
 459   for (Int_t ifit=0;ifit<16;ifit++){
 460     fitters[ifit]->Eval();
 461     fitters[ifit]->GetParameters(*(params[ifit]));
 462     fitters[ifit]->GetErrors(*(errs[ifit]));
 463     chi2N[ifit] = TMath::Sqrt(fitters[ifit]->GetChisquare()/(fitters[ifit]->GetNpoints()-1));
 464     (*(errs[ifit]))[0]*=chi2N[ifit];
 465     (*(errs[ifit]))[1]*=chi2N[ifit];
 466     if (ifit>=8) (*(errs[ifit]))[2]*=chi2N[ifit];  // second derivative
 467   }
 468   if (debug){
 469     (*debug)<<"dumpLinear"<<
 470       "alpha="<<fCurrentAlpha<<
 471       "xmean="<<xmean<<
 472       "y0T.="<<params[0]<<
 473       "y1T.="<<params[1]<<
 474       "z0T.="<<params[2]<<
 475       "z1T.="<<params[3]<<
 476       "y0O.="<<params[4]<<
 477       "y1O.="<<params[5]<<
 478       "z0O.="<<params[6]<<
 479       "z1O.="<<params[7]<<
 480       "y0T2.="<<params[8]<<
 481       "y1T2.="<<params[9]<<
 482       "z0T2.="<<params[10]<<
 483       "z1T2.="<<params[11]<<
 484       "y0O2.="<<params[12]<<
 485       "y1O2.="<<params[13]<<
 486       "z0O2.="<<params[14]<<
 487       "z1O2.="<<params[15]<<
 488       "y0TErr.="<<errs[0]<<
 489       "y1TErr.="<<errs[1]<<
 490       "z0TErr.="<<errs[2]<<
 491       "z1TErr.="<<errs[3]<<
 492       "y0OErr.="<<errs[4]<<
 493       "y1OErr.="<<errs[5]<<
 494       "z0OErr.="<<errs[6]<<
 495       "z1OErr.="<<errs[7]<<
 496       "y0T2Err.="<<errs[8]<<
 497       "y1T2Err.="<<errs[9]<<
 498       "z0T2Err.="<<errs[10]<<
 499       "z1T2Err.="<<errs[11]<<
 500       "y0O2Err.="<<errs[12]<<
 501       "y1O2Err.="<<errs[13]<<
 502       "z0O2Err.="<<errs[14]<<
 503       "z1O2Err.="<<errs[15]<<
 504       "chi2.="<<&chi2N<<
 505       "\n";
 506   }
 507   //
 508   //                    delta  alpha   y0     z0    ky   kz
 509   Double_t x[6]={0,0,0,0,0,0};
 510   x[1]=fCurrentAlpha;
 511   x[2]=(*params[0])[0];
 512   x[3]=(*params[2])[0];
 513   x[4]=(*params[0])[1];
 514   x[5]=(*params[2])[1];
 515   //
 516   //Delta y
 517   //
 518   x[0]= (*params[1])[0]-(*params[0])[0];
 519   fLinearTrackDelta[0]->Fill(x);
 520   x[0]/=TMath::Sqrt((*errs[1])[0]*(*errs[1])[0]+(*errs[0])[0]*(*errs[0])[0]);
 521   fLinearTrackPull[0]->Fill(x);
 522   //
 523   //Delta ky
 524   //
 525   x[0]= 100*((*params[1])[1]-(*params[0])[1]);
 526   fLinearTrackDelta[1]->Fill(x);
 527   x[0]/=100*TMath::Sqrt((*errs[1])[1]*(*errs[1])[1]+(*errs[0])[1]*(*errs[0])[1]);
 528   fLinearTrackPull[1]->Fill(x);
 529   //
 530   // ky2_0
 531   //
 532   x[0]= 100*100*((*params[8])[2]);
 533   fLinearTrackDelta[2]->Fill(x);
 534   x[0]/=100*100*TMath::Sqrt((*errs[8])[2]*(*errs[8])[2]);
 535   fLinearTrackPull[2]->Fill(x);
 536   //
 537   // ky2_1
 538   //
 539   x[0]= 100*100*((*params[9])[2]);
 540   fLinearTrackDelta[3]->Fill(x);
 541   x[0]/=100*100*TMath::Sqrt((*errs[9])[2]*(*errs[9])[2]);
 542   fLinearTrackPull[3]->Fill(x);
 543   //
 544   //
 545   //Delta z
 546   //
 547   x[0]= (*params[3])[0]-(*params[2])[0];
 548   fLinearTrackDelta[4]->Fill(x);
 549   x[0]/=TMath::Sqrt((*errs[3])[0]*(*errs[3])[0]+(*errs[2])[0]*(*errs[2])[0]);
 550   fLinearTrackPull[4]->Fill(x);
 551   //
 552   //Delta kz
 553   //
 554   x[0]= 100*((*params[3])[1]-(*params[2])[1]);
 555   fLinearTrackDelta[5]->Fill(x);
 556   x[0]/=100*TMath::Sqrt((*errs[3])[1]*(*errs[3])[1]+(*errs[2])[1]*(*errs[2])[1]);
 557   fLinearTrackPull[5]->Fill(x);
 558   //
 559   // kz2_0
 560   //
 561   x[0]= 100*100*((*params[10])[2]);
 562   fLinearTrackDelta[6]->Fill(x);
 563   x[0]/=100*100*TMath::Sqrt((*errs[10])[2]*(*errs[10])[2]);
 564   fLinearTrackPull[6]->Fill(x);
 565   //
 566   // kz2_1
 567   //
 568   x[0]= 100*100*((*params[11])[2]);
 569   fLinearTrackDelta[7]->Fill(x);
 570   x[0]/=100*100*TMath::Sqrt((*errs[11])[2]*(*errs[11])[2]);
 571   fLinearTrackPull[7]->Fill(x);
 572
 573   //
 574   // rms of track
 575   //
 576   x[0]= chi2N[0];
 577   fLinearTrackDelta[8]->Fill(x);
 578   x[0]= chi2N[1];
 579   fLinearTrackDelta[9]->Fill(x);
 580   x[0]= chi2N[2];
 581   fLinearTrackDelta[10]->Fill(x);
 582   x[0]= chi2N[3];
 583   fLinearTrackDelta[11]->Fill(x);
 584   //
 585
 586
 587   for (Int_t ifit=0; ifit<8;ifit++){
 588     delete fitters[ifit];
 589     delete params[ifit];
 590     delete fitters[ifit+8];
 591     delete params[ifit+8];
 592     delete errs[ifit];
 593     delete errs[ifit+8];
 594   }
 595 }
 596
 597
 598 void AliTPCkalmanFit::Propagate(TTreeSRedirector */*debug*/){
 599   //
 600   // Propagate the Kalman
 601   //
 602 }
 603
 604 void AliTPCkalmanFit::AddCovariance(const char * varName, Double_t sigma){
 605   //
 606   //
 607   //
 608   if (!fCalibCovar) return;
 609   if (!fCalibration) return;
 610   if (!fCalibration->FindObject(varName)) return;
 611   Int_t ncalibs = fCalibration->GetEntries();
 612   TString strVar(varName);
 613   for (Int_t icalib=0;icalib<ncalibs; icalib++){
 614     AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(icalib);
 615     if (strVar.CompareTo(transform->GetName())==0){
 616       (*fCalibCovar)(icalib,icalib)+=sigma*sigma;
 617     }
 618   }
 619 }
 620
 621
 622 void AliTPCkalmanFit::PropagateTime(Int_t time){
 623   //
 624   // Propagate the calibration in time
 625   // - Increase covariance matrix
 626   //
 627   if (!fCalibCovar) return;
 628   Int_t ncalibs = fCalibration->GetEntries();
 629   Double_t deltaT = 0;
 630   if (fLastTimeStamp>0) deltaT = (fLastTimeStamp-time)/(60*60.);  // delta T in hours
 631   fLastTimeStamp = time;
 632   for (Int_t icalib=0;icalib<ncalibs; icalib++){
 633     AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(icalib);
 634     if ((*fCalibCovar)(icalib,icalib)<transform->GetSigmaMax()*transform->GetSigmaMax())
 635       (*fCalibCovar)(icalib,icalib)+=  transform->GetSigma2Time()*TMath::Abs(deltaT);
 636   }
 637 }
 638
 639
 640 void  AliTPCkalmanFit::UpdateLinear(AliTrackPoint &point, TTreeSRedirector *debug){
 641   //
 642   // Update Kalman
 643   //
 644   //
 645   Int_t ncalibs = fCalibration->GetEntries();
 646   Int_t kNmeas  = 2;
 647   Int_t nelem   = ncalibs+4;
 648   TMatrixD &vecXk=*fLinearParam;     // X vector
 649   TMatrixD &covXk=*fLinearCovar;     // X covariance
 650   //
 651   TMatrixD matHk(kNmeas,nelem);     // vector to mesurement
 652   TMatrixD vecYk(kNmeas,1);         // Innovation or measurement residual
 653   TMatrixD vecZk(kNmeas,1);         // Innovation or measurement residual
 654   TMatrixD measR(kNmeas,kNmeas);
 655   TMatrixD matHkT(nelem,kNmeas);    // helper matrix Hk transpose
 656   TMatrixD matSk(kNmeas,kNmeas);    // Innovation (or residual) covariance
 657   TMatrixD matKk(nelem,kNmeas);     // Optimal Kalman gain
 658   TMatrixD covXk2(nelem,nelem);     // helper matrix
 659   TMatrixD covXk3(nelem,nelem);     // helper matrix
 660   TMatrixD mat1(nelem,nelem);
 661   //
 662   // reset matHk
 663   for (Int_t iel=0;iel<nelem;iel++){
 664     for (Int_t ip=0;ip<kNmeas;ip++) {
 665       matHk(ip,iel)=0;
 666     }
 667   }
 668   for (Int_t iel0=0;iel0<nelem;iel0++)
 669     for (Int_t iel1=0;iel1<nelem;iel1++){
 670       if (iel0!=iel1) mat1(iel0,iel1)=0;
 671       if (iel0==iel1) mat1(iel0,iel1)=1;
 672     }
 673   //
 674   // fill matrix Hk
 675   //
 676   Int_t volId = point.GetVolumeID();
 677   Double_t gxyz[3]={point.GetX(), point.GetY(),point.GetZ()};
 678   Double_t rxyz[3]={fCA*gxyz[0]+fSA*gxyz[1],-fSA*gxyz[0]+fCA*gxyz[1] ,point.GetZ()};
 679   //
 680   Double_t dxdydz[3]={0,0,0};
 681   Double_t rdxdydz[3]={0,0,0};
 682
 683   for (Int_t icalib=0;icalib<ncalibs;icalib++){
 684     AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(icalib);
 685     dxdydz[0] = transform->GetDeltaXYZ(0,volId, 1, gxyz[0], gxyz[1], gxyz[2]);
 686     dxdydz[1] = transform->GetDeltaXYZ(1,volId, 1, gxyz[0], gxyz[1], gxyz[2]);
 687     dxdydz[2] = transform->GetDeltaXYZ(2,volId, 1, gxyz[0], gxyz[1], gxyz[2]);
 688     rdxdydz[0] =  fCA*dxdydz[0]+fSA*dxdydz[1];
 689     rdxdydz[1] = -fSA*dxdydz[0]+fCA*dxdydz[1];
 690     rdxdydz[2] =  dxdydz[2];
 691     //
 692     matHk(0,icalib)= rdxdydz[1]-rdxdydz[0]*(*fLinearParam)(ncalibs+1,0);  // shift y + shift x * angleY
 693     matHk(1,icalib)= rdxdydz[2]-rdxdydz[0]*(*fLinearParam)(ncalibs+3,0);  // shift z + shift x * angleZ
 694   }
 695   matHk(0,ncalibs+0)=1;
 696   matHk(0,ncalibs+1)=rxyz[0];
 697   matHk(1,ncalibs+2)=1;
 698   matHk(1,ncalibs+3)=rxyz[0];
 699   //
 700   //
 701   //
 702   vecZk(0,0) =  rxyz[1];
 703   vecZk(1,0) =  rxyz[2];
 704   measR(0,0) = point.GetCov()[1]; measR(0,1)=0;
 705   measR(1,1) = point.GetCov()[2]; measR(1,0)=0;
 706   //
 707   vecYk = vecZk-matHk*vecXk;               // Innovation or measurement residual
 708   matHkT=matHk.T(); matHk.T();
 709   matSk = (matHk*(covXk*matHkT))+measR;    // Innovation (or residual) covariance
 710   matSk.Invert();
 711   matKk = (covXk*matHkT)*matSk;            //  Optimal Kalman gain
 712   //
 713   covXk2= (mat1-(matKk*matHk));
 714   covXk3 =  covXk2*covXk;
 715   //
 716   CheckCovariance(covXk3,100);
 717   vecXk += matKk*vecYk;                    //  updated vector
 718   covXk = covXk3;
 719   if (debug){ // dump debug info
 720     (*debug)<<"updateLinear"<<
 721       //
 722       "point.="<<&point<<
 723       "vecYk.="<<&vecYk<<
 724       "vecZk.="<<&vecZk<<
 725       "measR.="<<&measR<<
 726       "matHk.="<<&matHk<<
 727       "matHkT.="<<&matHkT<<
 728       "matSk.="<<&matSk<<
 729       "matKk.="<<&matKk<<
 730       "covXk2.="<<&covXk2<<
 731       "covXk.="<<&covXk<<
 732       "vecXk.="<<&vecXk<<
 733       "\n";
 734   }
 735 }
 736
 737
 738 AliTrackPointArray * AliTPCkalmanFit::SortPoints(AliTrackPointArray &points){
 739   //
 740   //Creates the array  - points sorted according radius - neccessay for kalman fit
 741   //
 742   //
 743   // 0. choose the frame - rotation angle
 744   //
 745   Int_t npoints = points.GetNPoints();
 746   if (npoints<1) return 0;
 747   Double_t currentAlpha = TMath::ATan2(points.GetY()[npoints-1]-points.GetY()[0], points.GetX()[npoints-1]-points.GetX()[0]);
 748   Double_t ca = TMath::Cos(currentAlpha);
 749   Double_t sa = TMath::Sin(currentAlpha);
 750   //
 751   // 1. sort the points
 752   //
 753   Double_t *rxvector = new Double_t[npoints];
 754   Int_t    *indexes  = new Int_t[npoints];
 755   for (Int_t ipoint=0; ipoint<npoints-1; ipoint++){
 756     rxvector[ipoint]=ca*points.GetX()[ipoint]+sa*points.GetY()[ipoint];
 757   }
 758   TMath::Sort(npoints, rxvector,indexes,kFALSE);
 759   AliTrackPoint point;
 760   AliTrackPointArray *pointsSorted= new AliTrackPointArray(npoints);
 761   for (Int_t ipoint=0; ipoint<npoints; ipoint++){
 762     if (!points.GetPoint(point,indexes[ipoint])) continue;
 763     pointsSorted->AddPoint(ipoint,&point);
 764   }
 765   delete [] rxvector;
 766   delete [] indexes;
 767   return pointsSorted;
 768 }
 769
 770
 771
 772 AliTrackPointArray * AliTPCkalmanFit::MakePointArrayLinear(Double_t alpha, Double_t y0, Double_t z0, Double_t ky, Double_t kz, Double_t err){
 773   //
 774   //
 775   //
 776   AliTrackPointArray array(500);
 777   Float_t cov[10];  // dummy covariance
 778   Int_t npoints=0;
 779   for (Int_t i=0;i<6;i++) cov[i]=0.001;
 780   for (Int_t i=0;i<500;i++){
 781     AliTrackPoint point(0, 0, 0, cov, 0,0,0);
 782     array.AddPoint(npoints, &point);
 783     npoints++;
 784   }
 785   npoints=0;
 786   for (Float_t ir = -245; ir<245; ir++){
 787     //
 788     //
 789     if (TMath::Abs(ir)<80) continue;
 790     Double_t ca = TMath::Cos(alpha);
 791     Double_t sa = TMath::Sin(alpha);
 792     Double_t lx = ir;
 793     Double_t ly = y0+lx*ky+gRandom->Gaus(0,err);
 794     Double_t lz = z0+lx*kz+gRandom->Gaus(0,err);
 795     Double_t gx = lx*ca-ly*sa;
 796     Double_t gy = lx*sa+ly*ca;
 797     Double_t gz = lz;
 798     Double_t galpha = TMath::ATan2(gy,gx);
 799     Int_t isec = TMath::Nint((9*galpha/TMath::Pi()+0.5));
 800     if (isec<0) isec+=18;
 801     if (gz<0) isec+=18;
 802     if (ir>130) isec+=36;
 803     //
 804     AliTrackPoint point(gx, gy, gz, cov, isec,0,0);
 805     array.AddPoint(npoints, &point);
 806     npoints++;
 807   }
 808   AliTrackPointArray *parray = new AliTrackPointArray(npoints);
 809   AliTrackPoint point;
 810   for (Int_t i=0;i<npoints;i++){
 811     array.GetPoint(point,i);
 812     parray->AddPoint(i,&point);
 813   }
 814
 815   return parray;
 816 }
 817
 818 void AliTPCkalmanFit::AddCalibration(AliTPCTransformation * calib){
 819   //
 820   // Add calibration
 821   //
 822   if (!fCalibration) fCalibration = new TObjArray(2000);
 823   fCalibration->AddLast(calib);
 824 }
 825
 826 Int_t AliTPCkalmanFit::GetTransformationIndex(const char * trName){
 827   //
 828   //
 829   //
 830   if (!fCalibration) return -1;
 831   if (!fCalibration->FindObject(trName)) return -1;
 832   //
 833   Int_t ncalibs = fCalibration->GetEntries();
 834   TString strVar(trName);
 835   for (Int_t icalib=0;icalib<ncalibs; icalib++){
 836     AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(icalib);
 837     if (strVar.CompareTo(transform->GetName())==0){
 838       return icalib;
 839     }
 840   }
 841   return -1;
 842 }
 843
 844
 845
 846 void AliTPCkalmanFit::ApplyCalibration(AliTrackPointArray *array, Double_t csign){
 847   //
 848   //
 849   //
 850   if (!fCalibration) return;
 851   Int_t ncalibs = fCalibration->GetEntries();
 852   if (ncalibs==0) return;
 853   Int_t npoints = array->GetNPoints();
 854   for (Int_t ipoint=0; ipoint<npoints; ipoint++){
 855     Int_t volId = array->GetVolumeID()[ipoint];
 856     Double_t xyz[3]={array->GetX()[ipoint], array->GetY()[ipoint],array->GetZ()[ipoint]};
 857     Double_t dxdydz[3]={0,0,0};
 858     for (Int_t icalib=0; icalib<ncalibs; icalib++){
 859       AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(icalib);
 860       dxdydz[0] += transform->GetDeltaXYZ(0,volId, transform->GetParam(), xyz[0], xyz[1], xyz[2]);
 861       dxdydz[1] += transform->GetDeltaXYZ(1,volId, transform->GetParam(), xyz[0], xyz[1], xyz[2]);
 862       dxdydz[2] += transform->GetDeltaXYZ(2,volId, transform->GetParam(), xyz[0], xyz[1], xyz[2]);
 863     }
 864     ((Float_t*)array->GetX())[ipoint]+=csign*dxdydz[0];
 865     ((Float_t*)array->GetY())[ipoint]+=csign*dxdydz[1];
 866     ((Float_t*)array->GetZ())[ipoint]+=csign*dxdydz[2];
 867   }
 868 }
 869
 870 Bool_t AliTPCkalmanFit::DumpCorelation(Double_t threshold,  const char *mask0, const char *mask1){
 871   //
 872   //
 873   //
 874   TMatrixD &mat = *fCalibCovar;
 875   Int_t nrow= mat.GetNrows();
 876   for (Int_t irow=0; irow<nrow; irow++){
 877     AliTPCTransformation * trans0 = GetTransformation(irow);
 878     TString  strName0(trans0->GetName());
 879     if (mask0){
 880       if (!strName0.Contains(mask0)) continue;
 881     }
 882     for (Int_t icol=irow+1; icol<nrow; icol++){
 883       AliTPCTransformation * trans1 = GetTransformation(icol);
 884       TString  strName1(trans1->GetName());
 885       if (mask1){
 886         if (!strName1.Contains(mask1)) continue;
 887       }
 888       //
 889       Double_t diag = TMath::Sqrt(TMath::Abs(mat(irow,irow)*mat(icol,icol)));
 890       if (diag<=0){
 891         printf("Negative covariance\t%d\t%d\t%f\n",irow,icol, mat(irow,icol));
 892         continue;
 893       }
 894       Double_t corr0 = mat(irow,icol)/diag;
 895       if (TMath::Abs(corr0)>threshold){
 896         printf("%d\t%d\t%s\t%s\t%f\t%f\t%f\n", irow,icol, trans0->GetName(), trans1->GetName(),
 897                TMath::Sqrt(mat(irow,irow)), TMath::Sqrt(mat(icol,icol)), corr0);
 898       }
 899     }
 900   }
 901   return (nrow>0);
 902 }
 903
 904 Bool_t AliTPCkalmanFit::DumpCalib(const char *mask, Float_t correlationCut){
 905   //
 906   // Print calibration entries - name, value, error
 907   //
 908   TMatrixD &mat = *fCalibCovar;
 909   Int_t nrow= mat.GetNrows();
 910   TString  strMask(mask);
 911   TVectorD vecCorrSum(nrow);
 912   for (Int_t irow=0; irow<nrow; irow++){
 913     vecCorrSum[irow]=0;
 914     for (Int_t icol=0; icol<nrow; icol++){
 915       if (icol==irow) continue;
 916       Double_t diag = TMath::Sqrt(TMath::Abs(mat(irow,irow)*mat(icol,icol)));
 917       Double_t corr0 = mat(irow,icol)/diag;
 918       vecCorrSum[irow]+=TMath::Abs(corr0);
 919     }
 920     vecCorrSum[irow]*=0.5;
 921   }
 922
 923   for (Int_t irow=0; irow<nrow; irow++){
 924     AliTPCTransformation * trans0 = GetTransformation(irow);
 925     TString  strName(trans0->GetName());
 926     if (mask){
 927       if (!strName.Contains(mask)) continue;
 928     }
 929     if (vecCorrSum[irow]<correlationCut) continue;
 930     printf("%d\t%s\t%f\t%f\t%f\n",
 931            irow,
 932            trans0->GetName(),
 933            (*fCalibParam)(irow,0),
 934            TMath::Sqrt(mat(irow,irow)),  vecCorrSum[irow]);
 935   }
 936   return (nrow>0);
 937 }
 938
 939
 940 Bool_t  AliTPCkalmanFit::CheckCovariance(TMatrixD &mat, Float_t /*maxEl*/){
 941   //
 942   // Check consistency of covariance matrix
 943   // + symetrize coavariance matrix
 944   Bool_t isOK=kTRUE;
 945   Int_t nrow= mat.GetNrows();
 946   for (Int_t irow=0; irow<nrow; irow++){
 947     if (mat(irow,irow)<=0){
 948       printf("Negative covariance\t%d\t%f\n",irow,mat(irow,irow));
 949       isOK=kFALSE;
 950     }
 951 //     if (mat(irow,irow)>maxEl){
 952 //       printf("Too big  covariance\t%d\t%f\n",irow,mat(irow,irow));
 953 //       isOK=kFALSE;
 954 //     }
 955
 956     for (Int_t icol=0; icol<nrow; icol++){
 957       //      if (mat(irow,irow)
 958       Double_t diag = TMath::Sqrt(TMath::Abs(mat(irow,irow)*mat(icol,icol)));
 959       if (diag<=0){
 960         printf("Negative covariance\t%d\t%d\t%f\n",irow,icol, mat(irow,icol));
 961         isOK=kFALSE;
 962         continue;
 963       }
 964       Double_t cov0 = mat(irow,icol)/diag;
 965       Double_t cov1 = mat(icol,irow)/diag;
 966       if (TMath::Abs(cov0)>1 || TMath::Abs(cov1)>1 ){
 967         printf("Covariance Problem %d\t%d\t%f\t%f\n",irow,icol, cov0, cov1);
 968         isOK=kFALSE;
 969       }
 970       if (TMath::Abs(cov0-cov1)>0.0000001){
 971         printf("Asymetry problem %d\t%d\t%f\t%f\n",irow,icol, cov0, cov1);
 972         isOK=kFALSE;
 973       }
 974       //
 975       // symetrize the covariance matrix
 976       Double_t mean = (mat(irow,icol)+ mat(icol,irow))*0.5;
 977       mat(irow,icol)=mean;
 978       mat(icol,irow)=mean;
 979     }
 980   }
 981   return isOK;
 982 }
 983
 984
 985 AliTPCkalmanFit *  AliTPCkalmanFit::Test(Int_t ntracks){
 986   //
 987   // This is test example
 988   //
 989
 990   //
 991   // 1. Setup transformation
 992   //
 993
 994
 995   TVectorD fpar(10);
 996   AliTPCTransformation * transformation=0;
 997   AliTPCkalmanFit * kalmanFit0 =  new AliTPCkalmanFit;
 998   AliTPCkalmanFit * kalmanFit2 =  new AliTPCkalmanFit;
 999   //
1000   // Radial scaling
1001   //
1002   for (Int_t iside=0; iside<=1; iside++)
1003     for (Int_t ipar0=0; ipar0<3; ipar0++)
1004       for (Int_t ipar1=0; ipar1<3; ipar1++){
1005         fpar[0]=ipar0;
1006         fpar[1]=ipar1;
1007         if (ipar0+ipar1==0) continue;
1008         Double_t param = (gRandom->Rndm()-0.5)*0.5;  // generate random parameters
1009         char tname[100];
1010         snprintf(tname,100,"tscalingR%d%dSide%d",ipar0,ipar1,iside);
1011         transformation = new AliTPCTransformation(tname,AliTPCTransformation::BitsSide(iside),"TPCscalingRPol",0,0,1);
1012         transformation->SetParams(0,5*0.25,0,&fpar);
1013         kalmanFit0->AddCalibration(transformation);
1014         transformation = new AliTPCTransformation(tname,AliTPCTransformation::BitsSide(iside),"TPCscalingRPol",0,0,1);
1015         transformation->SetParams(param,0.25,0,&fpar);
1016         kalmanFit2->AddCalibration(transformation);
1017       }
1018
1019
1020   //
1021   // 2. Init transformation
1022   //
1023   kalmanFit2->Init();
1024   kalmanFit0->Init();
1025
1026   //
1027   // 3. Run kalman filter
1028   //
1029   Int_t ncalibs = kalmanFit0->fCalibration->GetEntries();
1030   TVectorD err(ncalibs);
1031   TTreeSRedirector *pcstream = new TTreeSRedirector("kalmanfitTest.root");
1032   for (Int_t i=0;i<ntracks;i++){
1033     if (i%100==0) printf("%d\n",i);
1034      Double_t alpha = gRandom->Rndm()*TMath::TwoPi();
1035      Double_t   y0  = (gRandom->Rndm()-0.5)*180;
1036      Double_t   z0  = (gRandom->Rndm()-0.5)*250*2;
1037      Double_t   ky  = (gRandom->Rndm()-0.5)*1;
1038      Double_t   kz  = (gRandom->Rndm()-0.5)*1;
1039      //generate random TPC track
1040      AliTrackPointArray * array  = AliTPCkalmanFit::MakePointArrayLinear(alpha,y0,z0, ky, kz,0.04);
1041      AliTrackPointArray * arrayB = new AliTrackPointArray(*array);  // backup
1042      kalmanFit2->ApplyCalibration(array,1.);  // misalign ideal track
1043      for (Int_t icalib=0; icalib<ncalibs; icalib++){
1044        err[icalib] = TMath::Sqrt((*kalmanFit0->fCalibCovar)(icalib,icalib));
1045      }
1046      //
1047      (*pcstream)<<"dump0"<<
1048        "alpha="<<alpha<<
1049        "y0="<<y0<<
1050        "z0="<<z0<<
1051        "ky="<<ky<<
1052        "lz="<<kz<<
1053        "p.="<<array<<
1054        "pB.="<<arrayB<<
1055        "cparam.="<<kalmanFit0->fCalibParam<<
1056        "ccovar.="<<kalmanFit0->fCalibCovar<<
1057        "err.="<<&err<<
1058        "gparam.="<<kalmanFit2->fCalibParam<<
1059        "gcovar.="<<kalmanFit2->fCalibCovar<<
1060        "\n";
1061      if (i%20==0) {
1062        kalmanFit0->FitTrackLinear(*array,pcstream); // fit track - dump intermediate results
1063      }else{
1064        kalmanFit0->FitTrackLinear(*array,0);        // fit track + calibration
1065      }
1066      kalmanFit0->DumpTrackLinear(*array,pcstream);    // dump track residuals to the tree + fill histograms
1067   }
1068   pcstream->GetFile()->cd();
1069   kalmanFit0->Write("kalmanFit");
1070   delete pcstream;
1071   return kalmanFit0;
1072 }
1073
1074
1075 // Double_t AliTPCkalmanFit::GetTPCDeltaXYZ(Int_t coord, Int_t volID, Double_t x, Double_t y, Double_t z){
1076 //   //
1077 //   // function for visualization purposes
1078 //   //
1079 //   if (!fCalibration) return 0;
1080 //   Int_t ncalibs = fCalibration->GetEntries();
1081 //   if (ncalibs==0) return 0;
1082 //   Double_t dxdydz[3]={0,0,0};
1083 //   //
1084 //   if (volID<0){
1085 //     Double_t alpha       = TMath::ATan2(y,x);
1086 //     Double_t r           = TMath::Sqrt(y*y+x*x);
1087 //     volID                = TMath::Nint(9*alpha/TMath::Pi()-0.5);
1088 //     if (volID<0) volID+=18;
1089 //     if (z<0) volID+=18;
1090 //     if (r>120) volID+=36;
1091 //   }
1092 //   for (Int_t icalib=0; icalib<ncalibs; icalib++){
1093 //     AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(icalib);
1094 //     Double_t param = (*fCalibParam)(icalib,0);
1095 //     dxdydz[coord] += transform->GetDeltaXYZ(coord,volID, param, x, y,z);
1096 //   }
1097 //   return dxdydz[coord];
1098 // }
1099
1100 // Double_t AliTPCkalmanFit::SGetTPCDeltaXYZ(Int_t coord, Int_t volID, Double_t x, Double_t y, Double_t z){
1101 //   //
1102 //   //
1103 //   //
1104 //   if (AliTPCkalmanFit::fgInstance==0) return 0;
1105 //   return AliTPCkalmanFit::fgInstance->GetTPCDeltaXYZ(coord, volID,x,y,z);
1106 // }
1107
1108
1109
1110
1111
1112 Double_t AliTPCkalmanFit::GetTPCDeltaXYZ(Int_t coord, Int_t volID, Int_t icoordsys, Double_t x, Double_t y, Double_t z){
1113   //
1114   // function for visualization purposes
1115   //
1116   // coord - coordinate for output
1117   //       - 0 -X
1118   //         1 -Y
1119   //         2 -Z
1120   //         3 -R
1121   //         4 -RPhi
1122   //         5 -Z
1123   //
1124   //icoordsys - type of coordinate system for input
1125   //         - 0  - x,y,z
1126   //         - 1  - r,phi,z
1127   //
1128   if (!fCalibration) return 0;
1129   Int_t ncalibs = fCalibration->GetEntries();
1130   if (ncalibs==0) return 0;
1131   Double_t xyz[3]={0,0,0};
1132   Double_t dxdydz[6]={0,0,0,0,0,0};
1133   Double_t alpha=0;
1134   Double_t r=0;
1135   if(icoordsys==0){alpha=TMath::ATan2(y,x); r =TMath::Sqrt(y*y+x*x);}
1136   if(icoordsys==1){alpha=y; r=x;}
1137   Double_t ca    = TMath::Cos(alpha);
1138   Double_t sa    = TMath::Sin(alpha);
1139   if(icoordsys==0){xyz[0]=x; xyz[1]=y; xyz[2]=z;}
1140   if(icoordsys==1){xyz[0]=x*ca; xyz[1]=x*sa; xyz[2]=z;}
1141   //
1142   if (volID<0){
1143     // Double_t alpha       = TMath::ATan2(y,x);
1144     //Double_t r           = TMath::Sqrt(y*y+x*x);
1145     volID                = TMath::Nint(9*alpha/TMath::Pi()-0.5);
1146     if (volID<0) volID+=18;
1147     if (z<0) volID+=18;
1148     if (r>120) volID+=36;
1149   }
1150   for (Int_t icalib=0; icalib<ncalibs; icalib++){
1151     AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(icalib);
1152     Double_t param = (*fCalibParam)(icalib,0);
1153     for (Int_t icoord=0;icoord<6;icoord++){
1154       dxdydz[icoord] += transform->GetDeltaXYZ(icoord,volID, param, xyz[0],xyz[1],xyz[2]);
1155     }
1156   }
1157
1158   return dxdydz[coord];
1159 }
1160
1161
1162 Double_t AliTPCkalmanFit::SGetTPCDeltaXYZ(Int_t coord, Int_t volID, Int_t icoordsys, Double_t x, Double_t y, Double_t z){
1163   //
1164   //
1165   //
1166   if (AliTPCkalmanFit::fgInstance==0) return 0;
1167   return AliTPCkalmanFit::fgInstance->GetTPCDeltaXYZ(coord, volID, icoordsys,x,y,z);
1168 }
1169
1170
1171
1172
1173
1174 Double_t AliTPCkalmanFit::GetTPCtransXYZ(Int_t coord, Int_t volID, Int_t calibID, Int_t icoordsys, Double_t x, Double_t y, Double_t z){
1175
1176   Int_t ncalibs = fCalibration->GetEntries();
1177   if (calibID>=ncalibs) return 0;
1178   //Int_t volID=-1;
1179   //Double_t xyz[3]={x,y,z};
1180   Double_t r=0;
1181   Double_t alpha=0;
1182   if(icoordsys==0){r=TMath::Sqrt(x*x+y*y); alpha = TMath::ATan2(y,x);}
1183   if(icoordsys==1){r=x; alpha = y;}
1184   Double_t ca    = TMath::Cos(alpha);
1185   Double_t sa    = TMath::Sin(alpha);
1186   Double_t xyz[3]={0,0,0};
1187   if(icoordsys==0){xyz[0]=x;xyz[1]=y;xyz[2]=z;}
1188   if(icoordsys==1){xyz[0]=x*ca; xyz[1]=x*sa; xyz[2]=z;}
1189   //xyz[3]=param; xyz[4]=volID;
1190
1191   if (volID<0){
1192     //Double_t alpha       = TMath::ATan2(xyz[1],xyz[0]);
1193     //Double_t r           = TMath::Sqrt(xyz[1]*xyz[1]+xyz[0]*xyz[0]);
1194     volID                = TMath::Nint(9*alpha/TMath::Pi()-0.5);
1195     if (volID<0) volID+=18;
1196     if (xyz[2]<0) volID+=18;
1197     if (r>120) volID+=36;
1198   }
1199   AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(calibID);
1200   //transform->SetInstance(transform);
1201   Double_t param = (*fCalibParam)(calibID,0);
1202   Double_t delta = (Double_t)transform->GetDeltaXYZ(coord,volID, param, xyz[0],xyz[1],xyz[2]);
1203
1204   return delta;
1205 }
1206
1207 Double_t AliTPCkalmanFit::SGetTPCtransXYZ(Int_t coord, Int_t volID, Int_t calibID, Int_t icoordsys, Double_t x, Double_t y, Double_t z){
1208   //
1209   //
1210   //
1211   if (AliTPCkalmanFit::fgInstance==0) return 0;
1212   return AliTPCkalmanFit::fgInstance->GetTPCtransXYZ(coord, volID, calibID,icoordsys,x,y,z);
1213 }
1214
1215
1216 void AliTPCkalmanFit::MakeTreeTrans(TTreeSRedirector *debug, const char *treeName){
1217   //
1218   // Make the Tree before and after current calibration
1219   //
1220   if (!fCalibParam) {
1221     AliError("Kalman Fit not initialized");
1222     return;
1223   }
1224   //
1225   //
1226   //
1227   const Int_t ncalibs = fCalibration->GetEntries();
1228   TMatrixD dxdydz(ncalibs,5);
1229   Double_t * adx    = new Double_t[ncalibs];
1230   Double_t * ady    = new Double_t[ncalibs];
1231   Double_t * adz    = new Double_t[ncalibs];
1232   Double_t * adr    = new Double_t[ncalibs];
1233   Double_t * adrphi = new Double_t[ncalibs];
1234
1235   Double_t x[3];
1236   for (x[0]=-250.;x[0]<=250.;x[0]+=10.){
1237     for (x[1]=-250.;x[1]<=250.;x[1]+=10.){
1238       for (x[2]=-250.;x[2]<=250.;x[2]+=20.) {
1239         Double_t r=TMath::Sqrt(x[0]*x[0]+x[1]*x[1]);
1240         if (r<20) continue;
1241         if (r>260) continue;
1242         //Double_t z = x[2];
1243         Double_t phi=TMath::ATan2(x[1],x[0]);
1244         Double_t ca=TMath::Cos(phi);
1245         Double_t sa=TMath::Sin(phi);
1246         Double_t dx=0;
1247         Double_t dy=0;
1248         Double_t dz=0;
1249         Double_t dr=0;
1250         Double_t rdphi=0;
1251
1252         Int_t volID= TMath::Nint(9*phi/TMath::Pi()-0.5);
1253         if (volID<0) volID+=18;
1254         if (x[2]<0) volID+=18; //C-side
1255         if (r>120) volID+=36; //outer
1256         Double_t volalpha=(volID+0.5)*TMath::Pi()/9;
1257         Double_t cva=TMath::Cos(volalpha);
1258         Double_t sva=TMath::Sin(volalpha);
1259
1260         Double_t lx=x[0]*cva+x[1]*sva;
1261         Double_t ly=-x[0]*sva+x[1]*cva;
1262
1263
1264         for(Int_t icalib=0;icalib<ncalibs;icalib++){
1265           for(Int_t icoord=0;icoord<5;icoord++){
1266             dxdydz(icalib,icoord)= GetTPCtransXYZ(icoord, -1, icalib, 0, x[0], x[1], x[2]);
1267           }
1268         }
1269         dx=GetTPCDeltaXYZ(0, -1, 0, x[0], x[1], x[2]);
1270         dy=GetTPCDeltaXYZ(1, -1, 0, x[0], x[1], x[2]);
1271         dz=GetTPCDeltaXYZ(2, -1, 0, x[0], x[1], x[2]);
1272         dr=GetTPCDeltaXYZ(3, -1, 0, x[0], x[1], x[2]);
1273         rdphi=GetTPCDeltaXYZ(4, -1, 0, x[0], x[1], x[2]);
1274
1275
1276         if(debug){
1277           TTreeStream &cstream=
1278             (*debug)<<treeName<<
1279             "x="<<x[0]<<
1280             "y="<<x[1]<<
1281             "z="<<x[2]<<
1282             "r="<<r<<
1283             "ca="<<ca<<
1284             "sa="<<sa<<
1285             "lx="<<lx<<
1286             "ly="<<ly<<
1287             "sector="<<volID<<
1288             "phi="<<phi<<
1289             "dx="<<dx<<
1290             "dy="<<dy<<
1291             "dz="<<dz<<
1292             "dr="<<dr<<
1293             "rdphi="<<rdphi<<
1294             "dxdydz.="<<&dxdydz;
1295             for (Int_t icalib=0;icalib<ncalibs;icalib++){
1296               AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(icalib);
1297               char tname[1000];
1298               //
1299               snprintf(tname,1000,"dx%s=",transform->GetName());
1300               adx[icalib] =dxdydz(icalib,0);
1301               cstream<<tname<<adx[icalib];
1302               snprintf(tname,1000,"dy%s=",transform->GetName());
1303               ady[icalib] =dxdydz(icalib,1);
1304               cstream<<tname<<ady[icalib];
1305               snprintf(tname,1000,"dz%s=",transform->GetName());
1306               adz[icalib] =dxdydz(icalib,2);
1307               cstream<<tname<<adz[icalib];
1308               //
1309               snprintf(tname,1000,"dr%s=",transform->GetName());
1310               adr[icalib] =dxdydz(icalib,3);
1311               cstream<<tname<<adr[icalib];
1312               snprintf(tname,1000,"rdphi%s=",transform->GetName());
1313               adrphi[icalib] =dxdydz(icalib,4);
1314               cstream<<tname<<adrphi[icalib];
1315             }
1316             cstream<<"\n";
1317         }
1318       }
1319     }
1320     Printf("x0=%f finished",x[0]);
1321   }
1322   delete [] adx;//    = new Double_t[ncalibs];
1323   delete [] ady;//    = new Double_t[ncalibs];
1324   delete [] adz;//    = new Double_t[ncalibs];
1325   delete [] adr;//    = new Double_t[ncalibs];
1326   delete [] adrphi;// = new Double_t[ncalibs];
1327
1328 }