TPC/AliTPCkalmanAlign.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 /*
  17   TPC Kalman filter implementation for TPC sector alignment.
  18   Output of the AliTPCcalibAlign is used as a input for TPC global alignment.
  19   In AliTPCcalibAlign histograms - track parameter matching at sector boundaries are created.
  20   Each sector is aligned with 5 neighborhoud (sectors)
  21     1. Up-down    - 1
  22     2. Left-right - 4
  23
  24   Sector alignment parameters are obtained finding the alignment parameters, minimizing
  25   misalignmnet for all piars fo sectors.
  26
  27   Global minimization-  MakeGlobalAlign
  28
  29
  30   Example usage:
  31   gSystem->Load("libANALYSIS");
  32   gSystem->Load("libTPCcalib");
  33   //
  34   Int_t run=117220;
  35   .x ConfigCalibTrain.C(run)
  36   calibDB = AliTPCcalibDB::Instance()
  37
  38   AliTPCkalmanAlign kalmanAlign("TPC align", "TPC align");  // create the object
  39   kalmanAlign.ReadAlign("CalibObjects.root");               // read the calibration form file
  40   kalmanAlign.MakeGlobalAlign();                            // do kalman alignment
  41   kalmanAlign.DrawDeltaAlign();                             // make QA plot
  42   //
  43
  44
  45 */
  46 #include "TMath.h"
  47 #include "TTreeStream.h"
  48 #include "TGraph.h"
  49 #include "TGraphErrors.h"
  50 #include "TVectorD.h"
  51 #include "TClonesArray.h"
  52 #include "TCut.h"
  53 #include "TChain.h"
  54 #include "AliXRDPROOFtoolkit.h"
  55 #include "TLegend.h"
  56 #include "TCanvas.h"
  57
  58 #include "AliTPCcalibDB.h"
  59 #include "AliTPCCalROC.h"
  60 #include "AliCDBMetaData.h"
  61 #include "AliCDBId.h"
  62 #include "AliCDBManager.h"
  63 #include "AliCDBStorage.h"
  64 #include "AliCDBEntry.h"
  65 #include "AliAlignObjParams.h"
  66 #include "AliTPCROC.h"
  67 #include "AliTracker.h"
  68 #include "TFile.h"
  69 #include "TLinearFitter.h"
  70 #include "AliTPCcalibAlign.h"
  71 #include "TH1.h"
  72 #include "AliTPCCalPad.h"
  73 #include "AliTPCkalmanAlign.h"
  74 #include "TStatToolkit.h"
  75 #include "AliTPCPreprocessorOnline.h"
  76 #include "TPostScript.h"
  77 #include "AliGRPObject.h"
  78
  79 AliTPCkalmanAlign::AliTPCkalmanAlign():
  80   TNamed(),
  81   fCalibAlign(0),     // kalman alignnmnt
  82   fOriginalAlign(0),   // original alignment 0 read for the OCDB
  83   fNewAlign(0),
  84   fPadTime0(0),
  85   fFitCEGlobal(0),
  86   fFitCELocal(0)
  87 {
  88   //
  89   // Default constructor
  90   //
  91   for (Int_t i=0; i<4; i++){
  92     fDelta1D[i]=0;
  93     fCovar1D[i]=0;
  94   }
  95 }
  96
  97 AliTPCkalmanAlign::AliTPCkalmanAlign(const char* name, const char* title):
  98   TNamed(name,title),
  99   fCalibAlign(0),     // kalman alignnmnt
 100   fOriginalAlign(0),   // original alignment 0 read for the OCDB
 101   fNewAlign(0),
 102   fPadTime0(0),
 103   fFitCEGlobal(0),
 104   fFitCELocal(0)
 105 {
 106   //
 107   // Default constructor
 108   //
 109   for (Int_t i=0; i<4; i++){
 110     fDelta1D[i]=0;
 111     fCovar1D[i]=0;
 112   }
 113   fFitCEGlobal = new TObjArray(6);
 114   fFitCELocal  = new TObjArray(6);
 115   for (Int_t ipar=0; ipar<6;ipar++){
 116     fFitCEGlobal->AddAt(new TVectorD(36),ipar);
 117     fFitCELocal->AddAt(new TVectorD(36),ipar);
 118   }
 119 }
 120
 121 void AliTPCkalmanAlign::ReadAlign(const char *fname){
 122   //
 123   // Read old alignment used in the reco
 124   // and the residual histograms
 125   // WE ASSUME that the OCDB path is set in the same way as done in the calibration
 126   //
 127   TFile fcalib(fname);
 128   TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
 129   fCalibAlign=0;
 130   if (array) fCalibAlign=( AliTPCcalibAlign *)array->FindObject("alignTPC");
 131   fCalibAlign = ( AliTPCcalibAlign *)fcalib.Get("alignTPC");
 132   //
 133   // old alignment used
 134   AliCDBEntry * cdbEntry= AliCDBManager::Instance()->Get("TPC/Align/Data");
 135   fOriginalAlign =0;
 136   if (cdbEntry){
 137     fOriginalAlign = (TClonesArray*)cdbEntry->GetObject();
 138   }
 139
 140 }
 141
 142 void AliTPCkalmanAlign::BookAlign1D(TMatrixD &param, TMatrixD &covar, Double_t mean, Double_t sigma){
 143   //
 144   // Book Align 1D parameters and covariance
 145   //
 146   param.ResizeTo(72,1);
 147   covar.ResizeTo(72,72);
 148   for (Int_t i=0;i<72;i++){
 149     param(i,0)=mean;
 150     for (Int_t j=0;j<72;j++) covar(i,j)=0;
 151     covar(i,i)=sigma*sigma;
 152   }
 153 }
 154
 155
 156 void AliTPCkalmanAlign::UpdateAlign1D(Double_t delta, Double_t sigma, Int_t s1, Int_t s2, TMatrixD &vecXk, TMatrixD &covXk){
 157   //
 158   // Update 1D kalman filter
 159   //
 160   const Int_t knMeas=1;
 161   const Int_t knElem=72;
 162   TMatrixD mat1(72,72);            // update covariance matrix
 163   TMatrixD matHk(1,knElem);        // vector to mesurement
 164   TMatrixD vecYk(knMeas,1);        // Innovation or measurement residual
 165   TMatrixD matHkT(knElem,knMeas);  // helper matrix Hk transpose
 166   TMatrixD matSk(knMeas,knMeas);   // Innovation (or residual) covariance
 167   TMatrixD matKk(knElem,knMeas);   // Optimal Kalman gain
 168   TMatrixD covXk2(knElem,knElem);  // helper matrix
 169   TMatrixD covXk3(knElem,knElem);  // helper matrix
 170   TMatrixD vecZk(1,1);
 171   TMatrixD measR(1,1);
 172   vecZk(0,0)=delta;
 173   measR(0,0)=sigma*sigma;
 174   //
 175   // reset matHk
 176   for (Int_t iel=0;iel<knElem;iel++)
 177     for (Int_t ip=0;ip<knMeas;ip++) matHk(ip,iel)=0;
 178   //mat1
 179   for (Int_t iel=0;iel<knElem;iel++) {
 180     for (Int_t jel=0;jel<knElem;jel++) mat1(iel,jel)=0;
 181     mat1(iel,iel)=1;
 182   }
 183   //
 184   matHk(0, s1)=1;
 185   matHk(0, s2)=-1;
 186   vecYk = vecZk-matHk*vecXk;               // Innovation or measurement residual
 187   matHkT=matHk.T(); matHk.T();
 188   matSk = (matHk*(covXk*matHkT))+measR;    // Innovation (or residual) covariance
 189   matSk.Invert();
 190   matKk = (covXk*matHkT)*matSk;            //  Optimal Kalman gain
 191   vecXk += matKk*vecYk;                    //  updated vector
 192   covXk2= (mat1-(matKk*matHk));
 193   covXk3 =  covXk2*covXk;
 194   covXk = covXk3;
 195   Int_t nrows=covXk3.GetNrows();
 196
 197   for (Int_t irow=0; irow<nrows; irow++)
 198     for (Int_t icol=0; icol<nrows; icol++){
 199       // rounding problems - make matrix again symteric
 200       covXk(irow,icol)=(covXk3(irow,icol)+covXk3(icol,irow))*0.5;
 201     }
 202 }
 203
 204 void AliTPCkalmanAlign::UpdateAlign1D(Double_t delta, Double_t sigma, Int_t s1, TMatrixD &vecXk, TMatrixD &covXk){
 205   //
 206   // Update 1D kalman filter - with one measurement
 207   //
 208   const Int_t knMeas=1;
 209   const Int_t knElem=72;
 210   TMatrixD mat1(72,72);            // update covariance matrix
 211   TMatrixD matHk(1,knElem);        // vector to mesurement
 212   TMatrixD vecYk(knMeas,1);        // Innovation or measurement residual
 213   TMatrixD matHkT(knElem,knMeas);  // helper matrix Hk transpose
 214   TMatrixD matSk(knMeas,knMeas);   // Innovation (or residual) covariance
 215   TMatrixD matKk(knElem,knMeas);   // Optimal Kalman gain
 216   TMatrixD covXk2(knElem,knElem);  // helper matrix
 217   TMatrixD covXk3(knElem,knElem);  // helper matrix
 218   TMatrixD vecZk(1,1);
 219   TMatrixD measR(1,1);
 220   vecZk(0,0)=delta;
 221   measR(0,0)=sigma*sigma;
 222   //
 223   // reset matHk
 224   for (Int_t iel=0;iel<knElem;iel++)
 225     for (Int_t ip=0;ip<knMeas;ip++) matHk(ip,iel)=0;
 226   //mat1
 227   for (Int_t iel=0;iel<knElem;iel++) {
 228     for (Int_t jel=0;jel<knElem;jel++) mat1(iel,jel)=0;
 229     mat1(iel,iel)=1;
 230   }
 231   //
 232   matHk(0, s1)=1;
 233   vecYk = vecZk-matHk*vecXk;               // Innovation or measurement residual
 234   matHkT=matHk.T(); matHk.T();
 235   matSk = (matHk*(covXk*matHkT))+measR;    // Innovation (or residual) covariance
 236   matSk.Invert();
 237   matKk = (covXk*matHkT)*matSk;            //  Optimal Kalman gain
 238   vecXk += matKk*vecYk;                    //  updated vector
 239   covXk2= (mat1-(matKk*matHk));
 240   covXk3 =  covXk2*covXk;
 241   covXk = covXk3;
 242   Int_t nrows=covXk3.GetNrows();
 243
 244   for (Int_t irow=0; irow<nrows; irow++)
 245     for (Int_t icol=0; icol<nrows; icol++){
 246       // rounding problems - make matrix again symteric
 247       covXk(irow,icol)=(covXk3(irow,icol)+covXk3(icol,irow))*0.5;
 248     }
 249 }
 250
 251
 252
 253
 254 void AliTPCkalmanAlign::MakeGlobalAlign(){
 255   //
 256   // Combine all pairs of fitters and make global alignemnt
 257   //
 258
 259   AliTPCkalmanAlign &kalmanAlign=*this;
 260   TTreeSRedirector *pcstream = new TTreeSRedirector("AliTPCkalmanAlign.root");
 261   Int_t run = AliCDBManager::Instance()->GetRun();
 262   AliGRPObject * grp = AliTPCcalibDB::Instance()->GetGRP(run);
 263   Float_t bz = AliTracker::GetBz();
 264   UInt_t timeS = grp->GetTimeStart();
 265   UInt_t timeE = grp->GetTimeEnd();
 266   UInt_t  time = (timeS+timeE)/2;
 267
 268   //
 269   // get ce info
 270   //
 271   AliTPCCalPad *padTime0 = AliTPCcalibDB::Instance()->GetPadTime0();
 272   TVectorD paramCE[72];
 273   TMatrixD covarCE[72];
 274   Int_t  statUpDown=0;   // statistic up down
 275   Int_t  statLeftRight=0;   // statistic left-right
 276   Float_t chi2;
 277   for (Int_t isec=0; isec<72; isec++){
 278     AliTPCCalROC * roc0  = padTime0->GetCalROC(isec);
 279     roc0->GlobalFit(0,kFALSE,paramCE[isec],covarCE[isec],chi2,0);
 280     (*pcstream)<<"ceFit"<<
 281       "isec="<<isec<<
 282       "p0.="<<&paramCE[isec]<<
 283       "\n";
 284   }
 285
 286   DumpOldAlignment(pcstream);
 287   const Int_t kMinEntries=400;
 288   TMatrixD vec[5];
 289   TMatrixD cov[5];
 290   kalmanAlign.BookAlign1D(vec[0],cov[0], 0,0.5);
 291   kalmanAlign.BookAlign1D(vec[1],cov[1], 0,0.5);
 292   kalmanAlign.BookAlign1D(vec[2],cov[2], 0,0.01);
 293   kalmanAlign.BookAlign1D(vec[3],cov[3], 0,0.01);
 294   //
 295   TVectorD delta(5);
 296   TVectorD rms(5);
 297   TVectorD stat(5);
 298   TH1 * his=0;
 299   for (Int_t is0=0;is0<72;is0++)
 300     for (Int_t is1=0;is1<72;is1++){
 301       //
 302       //
 303       his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kY,is0,is1);
 304       if (!his) continue;
 305       if (his->GetEntries()<kMinEntries) continue;
 306       delta[0]=his->GetMean();
 307       rms[0]=his->GetRMS();
 308       stat[0]=his->GetEntries();
 309       kalmanAlign.UpdateAlign1D(his->GetMean(),his->GetRMS(),is0,is1, vec[0],cov[0]);
 310       //
 311       his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kZ,is0,is1);
 312       if (!his) continue;
 313       delta[1]=his->GetMean();
 314       rms[1]=his->GetRMS();
 315       stat[1]=his->GetEntries();
 316       kalmanAlign.UpdateAlign1D(his->GetMean(),his->GetRMS(),is0,is1, vec[1],cov[1]);
 317       //
 318       his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kPhi,is0,is1);
 319       if (!his) continue;
 320       delta[2] = his->GetMean();
 321       rms[2]=his->GetRMS();
 322       stat[2]=his->GetEntries();
 323       kalmanAlign.UpdateAlign1D(his->GetMean(),his->GetRMS(),is0,is1, vec[2],cov[2]);
 324       //
 325       his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kTheta,is0,is1);
 326       if (!his) continue;
 327       delta[3] = his->GetMean();
 328       rms[3]=his->GetRMS();
 329       stat[3]=his->GetEntries();
 330       kalmanAlign.UpdateAlign1D(his->GetMean(),his->GetRMS(),is0,is1, vec[3],cov[3]);
 331       if (is1==is0+36) statUpDown+=Int_t(stat[0]);
 332       if (is1==is0+35) statLeftRight+=Int_t(stat[0]);
 333     }
 334
 335   fDelta1D[0] = (TMatrixD*)vec[0].Clone();
 336   fDelta1D[1] = (TMatrixD*)vec[1].Clone();
 337   fDelta1D[2] = (TMatrixD*)vec[2].Clone();
 338   fDelta1D[3] = (TMatrixD*)vec[3].Clone();
 339   //
 340   fCovar1D[0] = (TMatrixD*)cov[0].Clone();
 341   fCovar1D[1] = (TMatrixD*)cov[1].Clone();
 342   fCovar1D[2] = (TMatrixD*)cov[2].Clone();
 343   fCovar1D[3] = (TMatrixD*)cov[3].Clone();
 344   statUpDown/=36;
 345   statLeftRight/=36;
 346   MakeNewAlignment(kTRUE);
 347   FitCE();
 348   for (Int_t is0=0;is0<72;is0++)
 349     for (Int_t is1=0;is1<72;is1++){
 350       Bool_t isPair=kFALSE;
 351       if (TMath::Abs(is0%18-is1%18)<2) isPair=kTRUE;
 352       if (TMath::Abs(is0%18-is1%18)==17) isPair=kTRUE;
 353       if (!isPair) continue;
 354       stat[0]=0; stat[1]=0; stat[2]=0; stat[3]=0;
 355       //
 356       //
 357       his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kY,is0,is1);
 358       if (his){
 359         delta[0]=his->GetMean();
 360         rms[0]=his->GetRMS();
 361         stat[0]=his->GetEntries();
 362       }
 363       //
 364       his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kZ,is0,is1);
 365       if (his) {
 366         delta[1]=his->GetMean();
 367         rms[1]=his->GetRMS();
 368         stat[1]=his->GetEntries();
 369       }
 370       //
 371       his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kPhi,is0,is1);
 372       if (his){
 373         delta[2] = his->GetMean();
 374         rms[2]=his->GetRMS();
 375         stat[2]=his->GetEntries();
 376       }
 377       //
 378       his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kTheta,is0,is1);
 379       if (his){
 380         delta[3] = his->GetMean();
 381         rms[3]=his->GetRMS();
 382         stat[3]=his->GetEntries();
 383       }
 384       TVectorD fceG[8],fceL[6];
 385       if (fFitCEGlobal)
 386         for (Int_t ipar=0; ipar<8;ipar++){
 387           fceG[ipar].ResizeTo(36);
 388           if (ipar<6) fceL[ipar].ResizeTo(36);
 389           if (fFitCEGlobal->At(ipar)){
 390             fceG[ipar]=*((TVectorD*)fFitCEGlobal->At(ipar));
 391             if (ipar<6){
 392               fceL[ipar]=*((TVectorD*)fFitCELocal->At(ipar));
 393             }
 394           }
 395         }
 396
 397       (*pcstream)<<"kalmanAlignDebug"<<
 398         "run="<<run<<            // runs
 399         "time="<<time<<          // time
 400         "timeE="<<timeE<<        // sart of tun time
 401         "timeS="<<timeS<<        // end od run time
 402         "bz="<<bz<<
 403         "is0="<<is0<<
 404         "is1="<<is1<<
 405         "delta.="<<&delta<<      // alignment deltas
 406         "rms.="<<&rms<<          // rms
 407         "stat.="<<&stat<<
 408         "vec0.="<<&vec[0]<<
 409         "vec1.="<<&vec[1]<<
 410         "vec2.="<<&vec[2]<<
 411         "vec3.="<<&vec[3]<<
 412         "pceIn0.="<<&paramCE[is0%36]<<        // default CE parameters
 413         "pceOut0.="<<&paramCE[is0%36+36]<<
 414         "pceIn1.="<<&paramCE[is1%36]<<
 415         "pceOut1.="<<&paramCE[is1%36+36]<<
 416         //                     // current CE parameters form last calibration - not used in Reco
 417         "fceG0.="<<&fceG[0]<<  // global fit of CE  - offset
 418         "fceG1.="<<&fceG[1]<<  // global fit of CE  - Gy gradient
 419         "fceG2.="<<&fceG[2]<<  // global fit of CE  - Gx gradient
 420         "fceG3.="<<&fceG[3]<<  // global fit of CE  - IROC-OROC offset
 421         "fceG4.="<<&fceG[4]<<  // global fit of CE  - commont slope LX
 422         "fceG5.="<<&fceG[5]<<  // global fit of CE  - delta slope LX
 423         "fceG6.="<<&fceG[6]<<  // global fit of CE  - common slope LY
 424         "fceG7.="<<&fceG[7]<<  // global fit of CE  - delta slope LY
 425         //
 426         "fceL0.="<<&fceL[0]<<  // local fit of CE  - offset to mean
 427         "fceL1.="<<&fceL[1]<<  // local fit of CE  - IROC-OROC offset
 428         "fceL2.="<<&fceL[2]<<  // local fit of CE  - common slope LX
 429         "fceL3.="<<&fceL[3]<<  // local fit of CE  - delta slope  LX
 430         "fceL4.="<<&fceL[4]<<  // local fit of CE  - common slope LY
 431         "fceL5.="<<&fceL[5]<<  // local fit of CE  - delta slope  LY
 432         "\n";
 433     }
 434
 435   (*pcstream)<<"runSummary"<<
 436     "run="<<run<<                      // run number
 437     "bz="<<bz<<                        // bz field
 438     "statUpDown="<<statUpDown<<        // stat up-down
 439     "statLeftRight="<<statLeftRight<<  // stat left-right
 440     "\n";
 441
 442   delete pcstream;
 443 }
 444
 445
 446
 447
 448
 449
 450 void AliTPCkalmanAlign::UpdateOCDBTime0( AliTPCCalPad  *pad, Int_t ustartRun, Int_t uendRun,  const char* storagePath ){
 451   //
 452   // Update OCDB
 453   // .x ConfigCalibTrain.C(117116)
 454   // AliTPCcalibDB::Instance()->GetPulserTmean()
 455   // pad->Add(-pad->GetMean())
 456   AliCDBMetaData *metaData= new AliCDBMetaData();
 457   metaData->SetObjectClassName("TObjArray");
 458   metaData->SetResponsible("Marian Ivanov");
 459   metaData->SetBeamPeriod(1);
 460   metaData->SetAliRootVersion("05-25-01"); //root version
 461   metaData->SetComment("Calibration of the z - time misalignment");
 462   AliCDBId* id1=NULL;
 463   id1=new AliCDBId("TPC/Calib/PadTime0", ustartRun, uendRun);
 464   AliCDBStorage* gStorage = AliCDBManager::Instance()->GetStorage(storagePath);
 465   gStorage->Put(pad, (*id1), metaData);
 466 }
 467
 468
 469
 470 void AliTPCkalmanAlign::DrawDeltaAlign(){
 471   //
 472   // Draw the reuslts of the alignment
 473   // Residual misalignment in respect with previous alignment shown
 474   //
 475   //
 476   TFile f("AliTPCkalmanAlign.root","update");
 477   TTree * treeDelta = (TTree*)f.Get("kalmanAlignDebug");
 478   TH1::AddDirectory(0);
 479   //
 480   treeDelta->SetAlias("sector","is0");
 481   treeDelta->SetAlias("dYmeas","10*(delta.fElements[0])");
 482   treeDelta->SetAlias("dZmeas","10*(delta.fElements[1])");
 483   treeDelta->SetAlias("dPhimeas","1000*(delta.fElements[2])");
 484   treeDelta->SetAlias("dThetameas","1000*(delta.fElements[3])");
 485   //
 486   treeDelta->SetAlias("dYfit","10*(vec0.fElements[is0]-vec0.fElements[is1])");
 487   treeDelta->SetAlias("dZfit","10*(vec1.fElements[is0]-vec1.fElements[is1])");
 488   treeDelta->SetAlias("dPhifit","1000*(vec2.fElements[is0]-vec2.fElements[is1])");
 489   treeDelta->SetAlias("dThetafit","1000*(vec3.fElements[is0]-vec3.fElements[is1])");
 490   //
 491   treeDelta->SetMarkerStyle(25);
 492   treeDelta->SetMarkerColor(4);
 493   treeDelta->SetLineColor(4);
 494   const char *type[3]={"up-down","left-right","right-left"};
 495   const char *gropt[3]={"alp","lp","lp"};
 496   const char *varTypeY[3]={"dYmeas-dYfit:sector","dYfit:sector","10*vec0.fElements[is0]:sector"};
 497   const char *varLegendY[3]={"#Delta_{y} fit residual","#Delta_{y} fit value difference","#Delta_{y} sector"};
 498   const char *varTypeZ[3]={"dZmeas-dZfit:sector","dZfit:sector","10*vec1.fElements[is0]:sector"};
 499   const char *varLegendZ[3]={"#Delta_{Z} fit residual","#Delta_{Z} fit value difference","#Delta_{Z} sector"};
 500   const char *varTypeT[3]={"dThetameas-dThetafit:sector","dThetafit:sector","1000*vec3.fElements[is0]:sector"};
 501   const char *varLegendT[3]={"#Delta_{#theta} fit residual","#Delta_{#theta} fit value difference","#Delta_{#theta} sector"};
 502   const char *varTypeP[3]={"dPhimeas-dPhifit:sector","dPhifit:sector","1000*vec2.fElements[is0]:sector"};
 503   const char *varLegendP[3]={"#Delta_{#phi} fit residual","#Delta_{#phi} fit value difference","#Delta_{#phi} sector"};
 504   TLegend *legend = 0;
 505   Int_t grcol[3]={2,1,4};
 506   Int_t entries=0;
 507   TGraph *grDelta[3]={0,0,0};
 508   TCanvas * canvasDy=new TCanvas("canvasDy","canvasDy",1200,800);
 509   TCanvas * canvasDz=new TCanvas("canvasDz","canvasDz",1200,800);
 510   TCanvas * canvasDphi=new TCanvas("canvasDphi","canvasDphi",1200,800);
 511   TCanvas * canvasDtheta=new TCanvas("canvasDtheta","canvasDtheta",1200,800);
 512   canvasDy->Divide(2,2);
 513   canvasDz->Divide(2,2);
 514   canvasDtheta->Divide(2,2);
 515   canvasDphi->Divide(2,2);
 516
 517   //
 518   // Dy
 519   //
 520   canvasDy->cd(1);
 521   treeDelta->Draw("dYmeas:dYfit");
 522   for (Int_t itype=0; itype<3; itype++){
 523     canvasDy->cd(itype+2);
 524     entries=treeDelta->Draw(varTypeY[itype],"is1==is0+36","goff");
 525     grDelta[0]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
 526     entries=treeDelta->Draw(varTypeY[itype],"is1==is0+35","goff");
 527     grDelta[1]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
 528     entries=treeDelta->Draw(varTypeY[itype],"is1==is0+37","goff");
 529     grDelta[2]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
 530     legend = new TLegend(0.5,0.7,0.9,0.9, varLegendY[itype]);
 531     for (Int_t i=0; i<3; i++) {
 532       grDelta[i]->SetMaximum(1.5);
 533       grDelta[i]->SetMinimum(-1);
 534       grDelta[i]->SetTitle(type[i]);
 535       grDelta[i]->SetMarkerColor(grcol[i]);
 536       grDelta[i]->SetLineColor(grcol[i]);
 537       grDelta[i]->SetMarkerStyle(25+i);
 538       grDelta[i]->GetXaxis()->SetTitle("sector");
 539       grDelta[i]->GetYaxis()->SetTitle("#Delta_{y} (mm)");
 540       if (itype==2 && i>0) continue;
 541       grDelta[i]->Draw(gropt[i]);
 542       legend->AddEntry(grDelta[i]);
 543     }
 544     legend->Draw();
 545   }
 546   //
 547   // Dz
 548   //
 549   canvasDz->cd();
 550   canvasDz->cd(1);
 551   treeDelta->Draw("dZmeas:dZfit");
 552   for (Int_t itype=0; itype<3; itype++){
 553     canvasDz->cd(itype+2);
 554     entries=treeDelta->Draw(varTypeZ[itype],"is1==is0+36","goff");
 555     grDelta[0]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
 556     entries=treeDelta->Draw(varTypeZ[itype],"is1==is0+35","goff");
 557     grDelta[1]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
 558     entries=treeDelta->Draw(varTypeZ[itype],"is1==is0+37","goff");
 559     grDelta[2]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
 560     legend = new TLegend(0.5,0.7,0.9,0.9, varLegendZ[itype]);
 561     for (Int_t i=0; i<3; i++) {
 562       grDelta[i]->SetMaximum(1.5);
 563       grDelta[i]->SetMinimum(-1);
 564       grDelta[i]->SetTitle(type[i]);
 565       grDelta[i]->SetMarkerColor(grcol[i]);
 566       grDelta[i]->SetLineColor(grcol[i]);
 567       grDelta[i]->SetMarkerStyle(25+i);
 568       grDelta[i]->GetXaxis()->SetTitle("sector");
 569       grDelta[i]->GetYaxis()->SetTitle("#Delta_{z} (mm)");
 570       if (itype==2 && i>0) continue;
 571       grDelta[i]->Draw(gropt[i]);
 572       legend->AddEntry(grDelta[i]);
 573     }
 574     legend->Draw();
 575   }
 576
 577   //
 578   // Dtheta
 579   //
 580   canvasDtheta->cd(1);
 581   treeDelta->Draw("dThetameas:dThetafit");
 582   for (Int_t itype=0; itype<3; itype++){
 583     canvasDtheta->cd(itype+2);
 584     entries=treeDelta->Draw(varTypeT[itype],"is1==is0+36","goff");
 585     grDelta[0]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
 586     entries=treeDelta->Draw(varTypeT[itype],"is1==is0+35","goff");
 587     grDelta[1]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
 588     entries=treeDelta->Draw(varTypeT[itype],"is1==is0+37","goff");
 589     grDelta[2]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
 590     legend = new TLegend(0.5,0.7,0.9,0.9, varLegendT[itype]);
 591     for (Int_t i=0; i<3; i++) {
 592       grDelta[i]->SetMaximum(4.);
 593       grDelta[i]->SetMinimum(-3.);
 594       grDelta[i]->SetTitle(type[i]);
 595       grDelta[i]->SetMarkerColor(grcol[i]);
 596       grDelta[i]->SetLineColor(grcol[i]);
 597       grDelta[i]->SetMarkerStyle(25+i);
 598       grDelta[i]->GetXaxis()->SetTitle("sector");
 599       grDelta[i]->GetYaxis()->SetTitle("#Delta_{#theta} (mrad)");
 600       if (itype==2 && i>0) continue;
 601       grDelta[i]->Draw(gropt[i]);
 602       legend->AddEntry(grDelta[i]);
 603     }
 604     legend->Draw();
 605   }
 606
 607   //
 608   // Dphi
 609   //
 610   canvasDphi->cd();
 611   canvasDphi->cd(1);
 612   treeDelta->Draw("dPhimeas:dPhifit");
 613   for (Int_t itype=0; itype<3; itype++){
 614     canvasDphi->cd(itype+2);
 615     entries=treeDelta->Draw(varTypeP[itype],"is1==is0+36","goff");
 616     grDelta[0]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
 617     entries=treeDelta->Draw(varTypeP[itype],"is1==is0+35","goff");
 618     grDelta[1]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
 619     entries=treeDelta->Draw(varTypeP[itype],"is1==is0+37","goff");
 620     grDelta[2]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
 621     legend = new TLegend(0.5,0.7,0.9,0.9, varLegendP[itype]);
 622     for (Int_t i=0; i<3; i++) {
 623       grDelta[i]->SetMaximum(4.);
 624       grDelta[i]->SetMinimum(-3.);
 625       grDelta[i]->SetTitle(type[i]);
 626       grDelta[i]->SetMarkerColor(grcol[i]);
 627       grDelta[i]->SetLineColor(grcol[i]);
 628       grDelta[i]->SetMarkerStyle(25+i);
 629       grDelta[i]->GetXaxis()->SetTitle("sector");
 630       grDelta[i]->GetYaxis()->SetTitle("#Delta_{#phi} (mrad)");
 631       if (itype==2 && i>0) continue;
 632       grDelta[i]->Draw(gropt[i]);
 633       legend->AddEntry(grDelta[i]);
 634     }
 635     legend->Draw();
 636   }
 637   //
 638   //
 639   f.cd();
 640   canvasDy->Write();
 641   canvasDz->Write();
 642   canvasDtheta->Write();
 643   canvasDphi->Write();
 644   //
 645   //
 646   //
 647   TPostScript *ps = new TPostScript("alignReport.ps", 112);
 648   ps->NewPage();
 649   canvasDy->Draw();
 650   ps->NewPage();
 651   canvasDz->Draw();
 652   ps->NewPage();
 653   canvasDtheta->Draw();
 654   ps->NewPage();
 655   canvasDphi->Draw();
 656   ps->Close();
 657   delete ps;
 658 }
 659
 660
 661
 662 void AliTPCkalmanAlign::DumpOldAlignment(TTreeSRedirector *pcstream){
 663   // Dump the content of old alignemnt
 664   // Expected that the old alignmnet is loaded
 665   //
 666   if (!fOriginalAlign) return;
 667   //
 668   TVectorD localTrans(3);
 669   TVectorD globalTrans(3);
 670   TVectorD localRot(3);
 671   TVectorD globalRot(3);
 672   AliGeomManager::ELayerID idLayer;
 673   Int_t idModule=0;
 674   //
 675   for (Int_t i=0; i<fOriginalAlign->GetEntries();i++){
 676     AliAlignObjParams *params = (AliAlignObjParams*)fOriginalAlign->At(i);
 677     params->GetVolUID(idLayer,idModule);
 678     params->GetLocalTranslation(localTrans.GetMatrixArray());
 679     params->GetLocalAngles(localRot.GetMatrixArray());
 680     params->GetTranslation(globalTrans.GetMatrixArray());
 681     params->GetAngles(globalRot.GetMatrixArray());
 682     Int_t sector=idModule;
 683     if (idLayer>7) sector+=36;
 684     (*pcstream)<<"oldAlign"<<
 685       //"idLayer="<<idLayer<<
 686       "idModule="<<idModule<<
 687       "sector="<<sector<<
 688       "lT.="<<&localTrans<<
 689       "gT.="<<&localTrans<<
 690       "lR.="<<&localRot<<
 691       "gR.="<<&globalRot<<
 692       "\n";
 693   }
 694 }
 695
 696
 697 void AliTPCkalmanAlign::MakeNewAlignment(Bool_t badd, TTreeSRedirector * pcstream){
 698   //
 699   // make a new Alignment entry
 700   //
 701   if (!fOriginalAlign) return;
 702   //
 703   TVectorD localTrans(3);
 704   TVectorD globalTrans(3);
 705   TVectorD localRot(3);
 706   TVectorD globalRot(3);
 707   //
 708   TVectorD localTransNew(3);   // new entries
 709   TVectorD globalTransNew(3);
 710   TVectorD localRotNew(3);
 711   TVectorD globalRotNew(3);
 712   //
 713   AliGeomManager::ELayerID idLayer;
 714   Int_t idModule=0;
 715   //
 716   fNewAlign = (TClonesArray*)fOriginalAlign->Clone();
 717   for (Int_t i=0; i<fOriginalAlign->GetEntries();i++){
 718     AliAlignObjParams *params = (AliAlignObjParams*)fOriginalAlign->At(i);
 719     //AliAlignObjParams *paramsNew = (AliAlignObjParams*)fNewAlign->At(i);
 720     params->GetVolUID(idLayer,idModule);
 721     Int_t sector=(Int_t)idModule;
 722     if (idLayer>7) sector+=36;
 723     params->GetLocalTranslation(localTrans.GetMatrixArray());
 724     params->GetLocalAngles(localRot.GetMatrixArray());
 725     params->GetTranslation(globalTrans.GetMatrixArray());
 726     params->GetAngles(globalRot.GetMatrixArray());
 727     //
 728     //
 729     //
 730     if (badd){ // addition if
 731       localTransNew=localTrans;
 732       localRotNew=localRot;
 733     }
 734     localTransNew[1]=localTransNew[1]-((*fDelta1D[0])(sector,0));
 735     localRot[0]  =localRot[0]-(*fDelta1D[2])(sector,0);
 736     //
 737     if (pcstream) (*pcstream)<<"alignParams"<<
 738       //"idLayer="<<idLayer<<
 739       "idModule="<<idModule<<
 740       "sector="<<sector<<
 741       "olT.="<<&localTrans<<
 742       "olR.="<<&localRot<<
 743       "ogT.="<<&localTrans<<
 744       "ogR.="<<&globalRot<<
 745       "nlT.="<<&localTransNew<<
 746       "nlR.="<<&localRotNew<<
 747       "ngT.="<<&localTransNew<<
 748       "ngR.="<<&globalRotNew<<
 749       "\n";
 750   }
 751 }
 752
 753
 754
 755 void AliTPCkalmanAlign::DrawAlignmentTrends(){
 756   //
 757   // Draw trends of alingment variables
 758   //
 759   /*
 760     //1.  Create a list of  align data
 761     //
 762     //2. Filter list
 763     AliXRDPROOFtoolkit::FilterList("align.list","AliTPCkalmanAlign.root kalmanAlignDebug",0);
 764
 765   */
 766   AliXRDPROOFtoolkit toolkit;
 767   TChain * chain = toolkit.MakeChainRandom("align.list.Good","kalmanAlignDebug",0,2000);
 768   TChain * chainRef = toolkit.MakeChainRandom("alignRef.list","kalmanAlignDebug",0,2000);
 769   chain->AddFriend(chainRef,"R");
 770   chainRef->AddFriend(chainRef,"T");
 771   //cuts
 772   TCut cutS="stat.fElements[0]>200&&stat.fElements[1]>200&&stat.fElements[3]>200&&stat.fElements[3]>200";   //statistic in the bin
 773   TCut cutST="T.stat.fElements[0]>200&&T.stat.fElements[1]>200&&T.stat.fElements[3]>200&&T.stat.fElements[3]>200";   //statistic in the bin
 774   //  TTree *tree  = chain->CopyTree(cutS);
 775   //TTree *treeR = chainRef->CopyTree(cutST);
 776
 777   TCanvas * canvasDy= new TCanvas("canvasDy","canvasDy");
 778   TH1 *his=0;
 779   TLegend *legend = 0;
 780   //  Int_t grcol[3]={2,1,4};
 781
 782   legend = new TLegend(0.7,0.6,0.9,0.9, "Alignment #Delta_{y}- Up-Down");
 783   for (Int_t isec=0; isec<18; isec+=2){
 784     chain->SetMarkerColor(1+(isec%5));
 785     chain->SetMarkerStyle(isec+20);
 786     chain->Draw("10*(delta.fElements[0]-R.delta.fElements[0]):run",cutS+Form("is1==is0+36&&is0==%d",isec),"profgoff");
 787     his = (TH1*)(chain->GetHistogram()->Clone());
 788     his->SetName(Form("#Delta_{Y} sector %d",isec));
 789     his->SetTitle(Form("#Delta_{Y} sector %d",isec));
 790     his->SetMaximum(1.);
 791     his->SetMinimum(-1.);
 792     his->GetYaxis()->SetTitle("#Delta_{y} (mm)");
 793     his->GetXaxis()->SetTitle("run Number");
 794     if (isec==0) his->Draw("");
 795     if (isec>0) his->Draw("same");
 796     legend->AddEntry(his);
 797   }
 798   legend->Draw();
 799   canvasDy->Draw();
 800 }
 801
 802
 803
 804
 805
 806
 807 void AliTPCkalmanAlign::FitCE(){
 808   //
 809   // fit CE
 810   // 1. Global fit - gy and gx
 811   // 2. Local X fit common
 812   // 3. Sector fit
 813   //
 814   AliTPCPreprocessorOnline * preprocesor = new AliTPCPreprocessorOnline;
 815   //
 816   AliTPCCalPad *padTime0 = AliTPCcalibDB::Instance()->GetPadTime0();
 817   AliTPCCalPad *padNoise = AliTPCcalibDB::Instance()->GetPadNoise();
 818   AliTPCCalPad * ceTmean = AliTPCcalibDB::Instance()->GetCETmean();  // CE information
 819   AliTPCCalPad * ceTrms  = AliTPCcalibDB::Instance()->GetCETrms();
 820   AliTPCCalPad * ceQmean = AliTPCcalibDB::Instance()->GetCEQmean();
 821   AliTPCCalPad * pulserTmean = AliTPCcalibDB::Instance()->GetPulserTmean(); //
 822   AliTPCCalPad * pulserTrms = AliTPCcalibDB::Instance()->GetPulserTrms();
 823   AliTPCCalPad * pulserQmean = AliTPCcalibDB::Instance()->GetPulserQmean();
 824   AliTPCCalPad * dmap0  = AliTPCcalibDB::Instance()->GetDistortionMap(0);   // distortion maps
 825   AliTPCCalPad * dmap1  = AliTPCcalibDB::Instance()->GetDistortionMap(1);
 826   AliTPCCalPad * dmap2  = AliTPCcalibDB::Instance()->GetDistortionMap(2);
 827   pulserTmean->Add(-pulserTmean->GetMean());
 828   //
 829   preprocesor->AddComponent(padTime0->Clone());
 830   preprocesor->AddComponent(padNoise->Clone());
 831   preprocesor->AddComponent(pulserTmean->Clone());
 832   preprocesor->AddComponent(pulserQmean->Clone());
 833   preprocesor->AddComponent(pulserTrms->Clone());
 834   preprocesor->AddComponent(ceTmean->Clone());
 835   preprocesor->AddComponent(ceQmean->Clone());
 836   preprocesor->AddComponent(ceTrms->Clone());
 837   preprocesor->AddComponent(dmap0->Clone());
 838   preprocesor->AddComponent(dmap1->Clone());
 839   preprocesor->AddComponent(dmap2->Clone());
 840   preprocesor->DumpToFile("cetmean.root");
 841
 842   TCut cutNoise="abs(PadNoise.fElements/PadNoise_Median-1)<0.3";
 843   TCut cutPulserT="abs(PulserTrms.fElements/PulserTrms_Median-1)<0.2";
 844   TCut cutPulserQ="abs(PulserQmean.fElements/PulserQmean_Median-1)<0.2";
 845   TCut cutCEQ="CEQmean.fElements>50";
 846   TCut cutCET="abs(CETmean.fElements)<2";
 847   TCut cutAll=cutNoise+cutPulserT+cutPulserQ+cutCEQ+cutCET;
 848   //
 849   //
 850   TFile * f = new TFile("cetmean.root");
 851   TTree * chain = (TTree*) f->Get("calPads");
 852   Int_t entries = chain->Draw("1",cutAll,"goff");
 853   if (entries<200000) return;  // no calibration available - pulser or CE or noise
 854
 855   TStatToolkit toolkit;
 856   Double_t chi2=0;
 857   Int_t    npoints=0;
 858   TVectorD param;
 859   TMatrixD covar;
 860   //
 861   // make a aliases
 862   AliTPCkalmanAlign::MakeAliasCE(chain);
 863   TString  fstringG="";              // global part
 864   //
 865   fstringG+="Gy++";                  // 1 - global y
 866   fstringG+="Gx++";                  // 2 - global x
 867   //
 868   fstringG+="isin++";                // 3 -delta IROC-OROC offset
 869   fstringG+="Lx++";                  // 4 -common slope
 870   fstringG+="Lx*isin++";             // 5 -delta slope
 871   fstringG+="Ly++";                  // 6 -common slope
 872   fstringG+="Ly*isin++";             // 7 -delta slope
 873   TVectorD vecG[2];
 874   TString * strFitG=0;
 875   TString * strFitLX=0;
 876   //
 877   strFitG = TStatToolkit::FitPlane(chain,"deltaT", fstringG.Data(),"sideA"+cutAll, chi2,npoints,vecG[0],covar,-1,0, 10000000, kFALSE);
 878   chain->SetAlias("tfitGA",strFitG->Data());
 879   strFitG->Tokenize("++")->Print();
 880   printf("chi2=%f\n",TMath::Sqrt(chi2/npoints));
 881   //
 882   strFitG = TStatToolkit::FitPlane(chain,"deltaT", fstringG.Data(),"sideC"+cutAll, chi2,npoints,vecG[1],covar,-1,0, 10000000, kFALSE);
 883   chain->SetAlias("tfitGC",strFitG->Data());
 884   strFitG->Tokenize("++")->Print();
 885   printf("chi2=%f\n",TMath::Sqrt(chi2/npoints));
 886   //
 887   AliTPCCalPad *padFitG =AliTPCCalPad::CreateCalPadFit("1++gy/500.++gx/500.++0+++0++0++0++0",vecG[0],vecG[1]);
 888   AliTPCCalPad *padFitLX=AliTPCCalPad::CreateCalPadFit("0++0++0++(sector<36)++(lx-133)/100++(sector<36)*(lx-133)/100.++(ly)/100++(sector<36)*(ly)/100.",vecG[0],vecG[1]);
 889   // swap a side and c side
 890   AliTPCCalPad *padFitGSwap =AliTPCCalPad::CreateCalPadFit("1++gy/500.++gx/500.++0+++0++0++0++0",vecG[1],vecG[0]);
 891   AliTPCCalPad *padFitLXSwap=AliTPCCalPad::CreateCalPadFit("0++0++0++(sector<36)++(lx-133)/100++(sector<36)*(lx-133)/100.++(ly)/100++(sector<36)*(ly)/100.",vecG[1],vecG[0]);
 892   padFitG->SetName("CEG");
 893   padFitLX->SetName("CELX");
 894   padFitGSwap->SetName("CEGS");
 895   padFitLXSwap->SetName("CELXS");
 896   preprocesor->AddComponent(padFitG->Clone());
 897   preprocesor->AddComponent(padFitLX->Clone());
 898   preprocesor->AddComponent(padFitGSwap->Clone());
 899   preprocesor->AddComponent(padFitLXSwap->Clone());
 900   preprocesor->DumpToFile("cetmean.root");   // add it to the file
 901   //
 902   // make local fits
 903   //
 904   f = new TFile("cetmean.root");
 905   chain = (TTree*) f->Get("calPads");
 906   AliTPCkalmanAlign::MakeAliasCE(chain);
 907   TString  fstringL="";              // local fit
 908   //                                 // 0. delta common
 909   fstringL+="isin++";                // 1. delta IROC-OROC offset
 910   fstringL+="Lx++";                  // 2. common slope
 911   fstringL+="Lx*isin++";             // 3. delta slope
 912   fstringL+="Ly++";                  // 4. common slope
 913   fstringL+="Ly*isin++";             // 5. delta slope
 914   TVectorD vecL[36];
 915   TVectorD dummy(6);
 916   AliTPCCalPad *padFitLCE = new AliTPCCalPad("LocalCE","LocalCE");
 917   AliTPCCalPad *padFitTmpCE;
 918   for (Int_t isec=0; isec<36; isec++){
 919     TCut cutSector=Form("(sector%%36)==%d",isec);
 920     strFitLX = TStatToolkit::FitPlane(chain,"deltaT-CEG.fElements-CELX.fElements", fstringL.Data(),cutSector+cutAll+"abs(deltaT-CEG.fElements-CELX.fElements)<0.4", chi2,npoints,vecL[isec],covar,-1,0, 10000000, kFALSE);
 921     printf("sec=%d\tchi2=%f\n",isec,TMath::Sqrt(chi2/npoints));
 922     delete strFitLX;
 923     //
 924     TString fitL=Form("((sector%%36)==%d)++((sector%%36)==%d)*(sector<36)++((sector%%36)==%d)*(lx-133)/100.++((sector%%36)==%d)*(sector<36)*(lx-133)/100.++((sector%%36)==%d)*(ly)/100.++((sector%%36)==%d)*(sector<36)*(ly)/100.",isec,isec,isec,isec,isec,isec);
 925     if (isec<18) padFitTmpCE=AliTPCCalPad::CreateCalPadFit(fitL.Data(),vecL[isec],dummy);
 926     if (isec>=18) padFitTmpCE=AliTPCCalPad::CreateCalPadFit(fitL.Data(),dummy,vecL[isec]);
 927     padFitLCE->Add(padFitTmpCE);
 928   }
 929   //
 930   padFitLCE->SetName("CELocal");
 931   preprocesor->AddComponent(padFitLCE->Clone());
 932   preprocesor->DumpToFile("cetmean.root");   // add it to the file
 933   //
 934   // write data to array
 935   //
 936   fFitCELocal  = new TObjArray(6);
 937   fFitCEGlobal = new TObjArray(8);
 938   for (Int_t ipar=0; ipar<8;ipar++){
 939     //
 940     fFitCEGlobal->AddAt(new TVectorD(36),ipar);
 941     TVectorD &fvecG = *((TVectorD*)fFitCEGlobal->At(ipar));
 942     for (Int_t isec=0; isec<36;isec++){
 943       if (isec<18)  fvecG[isec]=vecG[0][ipar];
 944       if (isec>=18) fvecG[isec]=vecG[1][ipar];
 945     }
 946     if (ipar<6){
 947       fFitCELocal->AddAt(new TVectorD(36),ipar);
 948       TVectorD &fvecL = *((TVectorD*)fFitCELocal->At(ipar));
 949       for (Int_t isec=0; isec<36;isec++){
 950         fvecL[isec]=vecL[isec][ipar];
 951       }
 952     }
 953   }
 954 }
 955
 956 void AliTPCkalmanAlign::MakeAliasCE(TTree * chain){
 957   //
 958   // make a aliases of pad variables
 959   //
 960   chain->SetAlias("side","(-1+(sector%36<18)*2)");
 961   chain->SetAlias("sideA","(sector%36<18)");
 962   chain->SetAlias("sideC","(sector%36>=18)");
 963   chain->SetAlias("isin","(sector<36)");
 964   chain->SetAlias("deltaT","CETmean.fElements-PulserTmean.fElements");
 965   chain->SetAlias("timeP","PulserTmean.fElements");
 966   chain->SetAlias("Gy","(gy.fElements/500.)");
 967   chain->SetAlias("Gx","(gx.fElements/500.)");
 968   chain->SetAlias("Lx","(lx.fElements-133)/100.");   // lx in meters
 969   chain->SetAlias("Ly","(ly.fElements)/100.");
 970   chain->SetAlias("La","(ly.fElements/lx.fElements/0.155)");
 971   chain->SetAlias("deltaT","(CETmean.fElements-PulserTmean.fElements)");
 972 }
 973
 974
 975
 976
 977 void AliTPCkalmanAlign::Update1D(Double_t delta, Double_t sigma, Int_t s1, TMatrixD &vecXk, TMatrixD &covXk){
 978   //
 979   // Update parameters and covariance - with one measurement
 980   //
 981   const Int_t knMeas=1;
 982   Int_t knElem=vecXk.GetNrows();
 983
 984   TMatrixD mat1(knElem,knElem);            // update covariance matrix
 985   TMatrixD matHk(1,knElem);        // vector to mesurement
 986   TMatrixD vecYk(knMeas,1);        // Innovation or measurement residual
 987   TMatrixD matHkT(knElem,knMeas);  // helper matrix Hk transpose
 988   TMatrixD matSk(knMeas,knMeas);   // Innovation (or residual) covariance
 989   TMatrixD matKk(knElem,knMeas);   // Optimal Kalman gain
 990   TMatrixD covXk2(knElem,knElem);  // helper matrix
 991   TMatrixD covXk3(knElem,knElem);  // helper matrix
 992   TMatrixD vecZk(1,1);
 993   TMatrixD measR(1,1);
 994   vecZk(0,0)=delta;
 995   measR(0,0)=sigma*sigma;
 996   //
 997   // reset matHk
 998   for (Int_t iel=0;iel<knElem;iel++)
 999     for (Int_t ip=0;ip<knMeas;ip++) matHk(ip,iel)=0;
1000   //mat1
1001   for (Int_t iel=0;iel<knElem;iel++) {
1002     for (Int_t jel=0;jel<knElem;jel++) mat1(iel,jel)=0;
1003     mat1(iel,iel)=1;
1004   }
1005   //
1006   matHk(0, s1)=1;
1007   vecYk = vecZk-matHk*vecXk;               // Innovation or measurement residual
1008   matHkT=matHk.T(); matHk.T();
1009   matSk = (matHk*(covXk*matHkT))+measR;    // Innovation (or residual) covariance
1010   matSk.Invert();
1011   matKk = (covXk*matHkT)*matSk;            //  Optimal Kalman gain
1012   vecXk += matKk*vecYk;                    //  updated vector
1013   covXk2= (mat1-(matKk*matHk));
1014   covXk3 =  covXk2*covXk;
1015   covXk = covXk3;
1016   Int_t nrows=covXk3.GetNrows();
1017
1018   for (Int_t irow=0; irow<nrows; irow++)
1019     for (Int_t icol=0; icol<nrows; icol++){
1020       // rounding problems - make matrix again symteric
1021       covXk(irow,icol)=(covXk3(irow,icol)+covXk3(icol,irow))*0.5;
1022     }
1023 }
1024
1025
1026 void   AliTPCkalmanAlign::Update1D(TString &input, TString filter, TVectorD &param, TMatrixD & covar, Double_t mean, Double_t sigma){
1027   //
1028   //  Update Parameters
1029   //  input  - variable name description
1030   //  filter - filter string
1031   //  param  - parameter vector
1032   //  covar  - covariance
1033   //  mean   - value to set
1034   //  sigma  - sigma of measurement
1035   TObjArray *array0= input.Tokenize("++");
1036   TObjArray *array1= filter.Tokenize("++");
1037   TMatrixD paramM(param.GetNrows(),1);
1038   for (Int_t i=0; i<array0->GetEntries(); i++){paramM(i,0)=param(i);}
1039
1040   for (Int_t i=0; i<array0->GetEntries(); i++){
1041     Bool_t isOK=kTRUE;
1042     TString str(array0->At(i)->GetName());
1043     for (Int_t j=0; j<array1->GetEntries(); j++){
1044       if (str.Contains(array1->At(j)->GetName())==0) isOK=kFALSE;
1045     }
1046     if (isOK) {
1047       AliTPCkalmanAlign::Update1D(mean, sigma, i, paramM, covar);
1048     }
1049   }
1050   for (Int_t i=0; i<array0->GetEntries(); i++){
1051     param(i)=paramM(i,0);
1052   }
1053 }
1054
1055
1056 TString  AliTPCkalmanAlign::FilterFit(TString &input, TString filter, TVectorD &param, TMatrixD & covar){
1057   //
1058   //
1059   //
1060   TObjArray *array0= input.Tokenize("++");
1061   TObjArray *array1= filter.Tokenize("++");
1062   //TString *presult=new TString("(0");
1063   TString result="(0.0";
1064   for (Int_t i=0; i<array0->GetEntries(); i++){
1065     Bool_t isOK=kTRUE;
1066     TString str(array0->At(i)->GetName());
1067     for (Int_t j=0; j<array1->GetEntries(); j++){
1068       if (str.Contains(array1->At(j)->GetName())==0) isOK=kFALSE;
1069     }
1070     if (isOK) {
1071       result+="+"+str;
1072       result+=Form("*(%f)",param[i+1]);
1073       printf("%f\t%f\t%s\n",param[i+1], TMath::Sqrt(covar(i+1,i+1)),str.Data());
1074     }
1075   }
1076   result+="-0.)";
1077   return result;
1078 }
1079
1080