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 }
 196
 197
 198
 199
 200 void AliTPCkalmanAlign::MakeGlobalAlign(){
 201   //
 202   // Combine all pairs of fitters and make global alignemnt
 203   //
 204
 205   AliTPCkalmanAlign &kalmanAlign=*this;
 206   TTreeSRedirector *pcstream = new TTreeSRedirector("AliTPCkalmanAlign.root");
 207   Int_t run = AliCDBManager::Instance()->GetRun();
 208   AliGRPObject * grp = AliTPCcalibDB::Instance()->GetGRP(run);
 209   Float_t bz = AliTracker::GetBz();
 210   UInt_t timeS = grp->GetTimeStart();
 211   UInt_t timeE = grp->GetTimeEnd();
 212   UInt_t  time = (timeS+timeE)/2;
 213
 214   //
 215   // get ce info
 216   //
 217   AliTPCCalPad *padTime0 = AliTPCcalibDB::Instance()->GetPadTime0();
 218   TVectorD paramCE[72];
 219   TMatrixD covarCE[72];
 220   Int_t  statUpDown=0;   // statistic up down
 221   Int_t  statLeftRight=0;   // statistic left-right
 222   Float_t chi2;
 223   for (Int_t isec=0; isec<72; isec++){
 224     AliTPCCalROC * roc0  = padTime0->GetCalROC(isec);
 225     roc0->GlobalFit(0,kFALSE,paramCE[isec],covarCE[isec],chi2,0);
 226     (*pcstream)<<"ceFit"<<
 227       "isec="<<isec<<
 228       "p0.="<<&paramCE[isec]<<
 229       "\n";
 230   }
 231
 232   DumpOldAlignment(pcstream);
 233   const Int_t kMinEntries=400;
 234   TMatrixD vec[5];
 235   TMatrixD cov[5];
 236   kalmanAlign.BookAlign1D(vec[0],cov[0], 0,0.5);
 237   kalmanAlign.BookAlign1D(vec[1],cov[1], 0,0.5);
 238   kalmanAlign.BookAlign1D(vec[2],cov[2], 0,0.01);
 239   kalmanAlign.BookAlign1D(vec[3],cov[3], 0,0.01);
 240   //
 241   TVectorD delta(5);
 242   TVectorD rms(5);
 243   TVectorD stat(5);
 244   TH1 * his=0;
 245   for (Int_t is0=0;is0<72;is0++)
 246     for (Int_t is1=0;is1<72;is1++){
 247       //
 248       //
 249       his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kY,is0,is1);
 250       if (!his) continue;
 251       if (his->GetEntries()<kMinEntries) continue;
 252       delta[0]=his->GetMean();
 253       rms[0]=his->GetRMS();
 254       stat[0]=his->GetEntries();
 255       kalmanAlign.UpdateAlign1D(his->GetMean(),his->GetRMS(),is0,is1, vec[0],cov[0]);
 256       //
 257       his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kZ,is0,is1);
 258       if (!his) continue;
 259       delta[1]=his->GetMean();
 260       rms[1]=his->GetRMS();
 261       stat[1]=his->GetEntries();
 262       kalmanAlign.UpdateAlign1D(his->GetMean(),his->GetRMS(),is0,is1, vec[1],cov[1]);
 263       //
 264       his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kPhi,is0,is1);
 265       if (!his) continue;
 266       delta[2] = his->GetMean();
 267       rms[2]=his->GetRMS();
 268       stat[2]=his->GetEntries();
 269       kalmanAlign.UpdateAlign1D(his->GetMean(),his->GetRMS(),is0,is1, vec[2],cov[2]);
 270       //
 271       his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kTheta,is0,is1);
 272       if (!his) continue;
 273       delta[3] = his->GetMean();
 274       rms[3]=his->GetRMS();
 275       stat[3]=his->GetEntries();
 276       kalmanAlign.UpdateAlign1D(his->GetMean(),his->GetRMS(),is0,is1, vec[3],cov[3]);
 277       if (is1==is0+36) statUpDown+=Int_t(stat[0]);
 278       if (is1==is0+35) statLeftRight+=Int_t(stat[0]);
 279     }
 280
 281   fDelta1D[0] = (TMatrixD*)vec[0].Clone();
 282   fDelta1D[1] = (TMatrixD*)vec[1].Clone();
 283   fDelta1D[2] = (TMatrixD*)vec[2].Clone();
 284   fDelta1D[3] = (TMatrixD*)vec[3].Clone();
 285   //
 286   fCovar1D[0] = (TMatrixD*)cov[0].Clone();
 287   fCovar1D[1] = (TMatrixD*)cov[1].Clone();
 288   fCovar1D[2] = (TMatrixD*)cov[2].Clone();
 289   fCovar1D[3] = (TMatrixD*)cov[3].Clone();
 290   statUpDown/=36;
 291   statLeftRight/=36;
 292   MakeNewAlignment(kTRUE);
 293   FitCE();
 294   for (Int_t is0=0;is0<72;is0++)
 295     for (Int_t is1=0;is1<72;is1++){
 296       Bool_t isPair=kFALSE;
 297       if (TMath::Abs(is0%18-is1%18)<2) isPair=kTRUE;
 298       if (TMath::Abs(is0%18-is1%18)==17) isPair=kTRUE;
 299       if (!isPair) continue;
 300       stat[0]=0; stat[1]=0; stat[2]=0; stat[3]=0;
 301       //
 302       //
 303       his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kY,is0,is1);
 304       if (his){
 305         delta[0]=his->GetMean();
 306         rms[0]=his->GetRMS();
 307         stat[0]=his->GetEntries();
 308       }
 309       //
 310       his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kZ,is0,is1);
 311       if (his) {
 312         delta[1]=his->GetMean();
 313         rms[1]=his->GetRMS();
 314         stat[1]=his->GetEntries();
 315       }
 316       //
 317       his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kPhi,is0,is1);
 318       if (his){
 319         delta[2] = his->GetMean();
 320         rms[2]=his->GetRMS();
 321         stat[2]=his->GetEntries();
 322       }
 323       //
 324       his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kTheta,is0,is1);
 325       if (his){
 326         delta[3] = his->GetMean();
 327         rms[3]=his->GetRMS();
 328         stat[3]=his->GetEntries();
 329       }
 330       TVectorD fceG[8],fceL[6];
 331       if (fFitCEGlobal)
 332         for (Int_t ipar=0; ipar<8;ipar++){
 333           fceG[ipar].ResizeTo(36);
 334           if (ipar<6) fceL[ipar].ResizeTo(36);
 335           if (fFitCEGlobal->At(ipar)){
 336             fceG[ipar]=*((TVectorD*)fFitCEGlobal->At(ipar));
 337             if (ipar<6){
 338               fceL[ipar]=*((TVectorD*)fFitCELocal->At(ipar));
 339             }
 340           }
 341         }
 342
 343       (*pcstream)<<"kalmanAlignDebug"<<
 344         "run="<<run<<            // runs
 345         "time="<<time<<          // time
 346         "timeE="<<timeE<<        // sart of tun time
 347         "timeS="<<timeS<<        // end od run time
 348         "bz="<<bz<<
 349         "is0="<<is0<<
 350         "is1="<<is1<<
 351         "delta.="<<&delta<<      // alignment deltas
 352         "rms.="<<&rms<<          // rms
 353         "stat.="<<&stat<<
 354         "vec0.="<<&vec[0]<<
 355         "vec1.="<<&vec[1]<<
 356         "vec2.="<<&vec[2]<<
 357         "vec3.="<<&vec[3]<<
 358         "pceIn0.="<<&paramCE[is0%36]<<        // default CE parameters
 359         "pceOut0.="<<&paramCE[is0%36+36]<<
 360         "pceIn1.="<<&paramCE[is1%36]<<
 361         "pceOut1.="<<&paramCE[is1%36+36]<<
 362         //                     // current CE parameters form last calibration - not used in Reco
 363         "fceG0.="<<&fceG[0]<<  // global fit of CE  - offset
 364         "fceG1.="<<&fceG[1]<<  // global fit of CE  - Gy gradient
 365         "fceG2.="<<&fceG[2]<<  // global fit of CE  - Gx gradient
 366         "fceG3.="<<&fceG[3]<<  // global fit of CE  - IROC-OROC offset
 367         "fceG4.="<<&fceG[4]<<  // global fit of CE  - commont slope LX
 368         "fceG5.="<<&fceG[5]<<  // global fit of CE  - delta slope LX
 369         "fceG6.="<<&fceG[6]<<  // global fit of CE  - common slope LY
 370         "fceG7.="<<&fceG[7]<<  // global fit of CE  - delta slope LY
 371         //
 372         "fceL0.="<<&fceL[0]<<  // local fit of CE  - offset to mean
 373         "fceL1.="<<&fceL[1]<<  // local fit of CE  - IROC-OROC offset
 374         "fceL2.="<<&fceL[2]<<  // local fit of CE  - common slope LX
 375         "fceL3.="<<&fceL[3]<<  // local fit of CE  - delta slope  LX
 376         "fceL4.="<<&fceL[4]<<  // local fit of CE  - common slope LY
 377         "fceL5.="<<&fceL[5]<<  // local fit of CE  - delta slope  LY
 378         "\n";
 379     }
 380
 381   (*pcstream)<<"runSummary"<<
 382     "run="<<run<<                      // run number
 383     "bz="<<bz<<                        // bz field
 384     "statUpDown="<<statUpDown<<        // stat up-down
 385     "statLeftRight="<<statLeftRight<<  // stat left-right
 386     "\n";
 387
 388   delete pcstream;
 389 }
 390
 391
 392
 393
 394
 395
 396 void AliTPCkalmanAlign::UpdateOCDBTime0( AliTPCCalPad  *pad, Int_t ustartRun, Int_t uendRun,  const char* storagePath ){
 397   //
 398   // Update OCDB
 399   // .x ConfigCalibTrain.C(117116)
 400   // AliTPCcalibDB::Instance()->GetPulserTmean()
 401   // pad->Add(-pad->GetMean())
 402   AliCDBMetaData *metaData= new AliCDBMetaData();
 403   metaData->SetObjectClassName("TObjArray");
 404   metaData->SetResponsible("Marian Ivanov");
 405   metaData->SetBeamPeriod(1);
 406   metaData->SetAliRootVersion("05-25-01"); //root version
 407   metaData->SetComment("Calibration of the z - time misalignment");
 408   AliCDBId* id1=NULL;
 409   id1=new AliCDBId("TPC/Calib/PadTime0", ustartRun, uendRun);
 410   AliCDBStorage* gStorage = AliCDBManager::Instance()->GetStorage(storagePath);
 411   gStorage->Put(pad, (*id1), metaData);
 412 }
 413
 414
 415
 416 void AliTPCkalmanAlign::DrawDeltaAlign(){
 417   //
 418   // Draw the reuslts of the alignment
 419   // Residual misalignment in respect with previous alignment shown
 420   //
 421   //
 422   TFile f("AliTPCkalmanAlign.root","update");
 423   TTree * treeDelta = (TTree*)f.Get("kalmanAlignDebug");
 424   TH1::AddDirectory(0);
 425   //
 426   treeDelta->SetAlias("sector","is0");
 427   treeDelta->SetAlias("dYmeas","10*(delta.fElements[0])");
 428   treeDelta->SetAlias("dZmeas","10*(delta.fElements[1])");
 429   treeDelta->SetAlias("dPhimeas","1000*(delta.fElements[2])");
 430   treeDelta->SetAlias("dThetameas","1000*(delta.fElements[3])");
 431   //
 432   treeDelta->SetAlias("dYfit","10*(vec0.fElements[is0]-vec0.fElements[is1])");
 433   treeDelta->SetAlias("dZfit","10*(vec1.fElements[is0]-vec1.fElements[is1])");
 434   treeDelta->SetAlias("dPhifit","1000*(vec2.fElements[is0]-vec2.fElements[is1])");
 435   treeDelta->SetAlias("dThetafit","1000*(vec3.fElements[is0]-vec3.fElements[is1])");
 436   //
 437   treeDelta->SetMarkerStyle(25);
 438   treeDelta->SetMarkerColor(4);
 439   treeDelta->SetLineColor(4);
 440   const char *type[3]={"up-down","left-right","right-left"};
 441   const char *gropt[3]={"alp","lp","lp"};
 442   const char *varTypeY[3]={"dYmeas-dYfit:sector","dYfit:sector","10*vec0.fElements[is0]:sector"};
 443   const char *varLegendY[3]={"#Delta_{y} fit residual","#Delta_{y} fit value difference","#Delta_{y} sector"};
 444   const char *varTypeZ[3]={"dZmeas-dZfit:sector","dZfit:sector","10*vec1.fElements[is0]:sector"};
 445   const char *varLegendZ[3]={"#Delta_{Z} fit residual","#Delta_{Z} fit value difference","#Delta_{Z} sector"};
 446   const char *varTypeT[3]={"dThetameas-dThetafit:sector","dThetafit:sector","1000*vec3.fElements[is0]:sector"};
 447   const char *varLegendT[3]={"#Delta_{#theta} fit residual","#Delta_{#theta} fit value difference","#Delta_{#theta} sector"};
 448   const char *varTypeP[3]={"dPhimeas-dPhifit:sector","dPhifit:sector","1000*vec2.fElements[is0]:sector"};
 449   const char *varLegendP[3]={"#Delta_{#phi} fit residual","#Delta_{#phi} fit value difference","#Delta_{#phi} sector"};
 450   TLegend *legend = 0;
 451   Int_t grcol[3]={2,1,4};
 452   Int_t entries=0;
 453   TGraph *grDelta[3]={0,0,0};
 454   TCanvas * canvasDy=new TCanvas("canvasDy","canvasDy",1200,800);
 455   TCanvas * canvasDz=new TCanvas("canvasDz","canvasDz",1200,800);
 456   TCanvas * canvasDphi=new TCanvas("canvasDphi","canvasDphi",1200,800);
 457   TCanvas * canvasDtheta=new TCanvas("canvasDtheta","canvasDtheta",1200,800);
 458   canvasDy->Divide(2,2);
 459   canvasDz->Divide(2,2);
 460   canvasDtheta->Divide(2,2);
 461   canvasDphi->Divide(2,2);
 462
 463   //
 464   // Dy
 465   //
 466   canvasDy->cd(1);
 467   treeDelta->Draw("dYmeas:dYfit");
 468   for (Int_t itype=0; itype<3; itype++){
 469     canvasDy->cd(itype+2);
 470     entries=treeDelta->Draw(varTypeY[itype],"is1==is0+36","goff");
 471     grDelta[0]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
 472     entries=treeDelta->Draw(varTypeY[itype],"is1==is0+35","goff");
 473     grDelta[1]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
 474     entries=treeDelta->Draw(varTypeY[itype],"is1==is0+37","goff");
 475     grDelta[2]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
 476     legend = new TLegend(0.5,0.7,0.9,0.9, varLegendY[itype]);
 477     for (Int_t i=0; i<3; i++) {
 478       grDelta[i]->SetMaximum(1.5);
 479       grDelta[i]->SetMinimum(-1);
 480       grDelta[i]->SetTitle(type[i]);
 481       grDelta[i]->SetMarkerColor(grcol[i]);
 482       grDelta[i]->SetLineColor(grcol[i]);
 483       grDelta[i]->SetMarkerStyle(25+i);
 484       grDelta[i]->GetXaxis()->SetTitle("sector");
 485       grDelta[i]->GetYaxis()->SetTitle("#Delta_{y} (mm)");
 486       if (itype==2 && i>0) continue;
 487       grDelta[i]->Draw(gropt[i]);
 488       legend->AddEntry(grDelta[i]);
 489     }
 490     legend->Draw();
 491   }
 492   //
 493   // Dz
 494   //
 495   canvasDz->cd();
 496   canvasDz->cd(1);
 497   treeDelta->Draw("dZmeas:dZfit");
 498   for (Int_t itype=0; itype<3; itype++){
 499     canvasDz->cd(itype+2);
 500     entries=treeDelta->Draw(varTypeZ[itype],"is1==is0+36","goff");
 501     grDelta[0]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
 502     entries=treeDelta->Draw(varTypeZ[itype],"is1==is0+35","goff");
 503     grDelta[1]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
 504     entries=treeDelta->Draw(varTypeZ[itype],"is1==is0+37","goff");
 505     grDelta[2]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
 506     legend = new TLegend(0.5,0.7,0.9,0.9, varLegendZ[itype]);
 507     for (Int_t i=0; i<3; i++) {
 508       grDelta[i]->SetMaximum(1.5);
 509       grDelta[i]->SetMinimum(-1);
 510       grDelta[i]->SetTitle(type[i]);
 511       grDelta[i]->SetMarkerColor(grcol[i]);
 512       grDelta[i]->SetLineColor(grcol[i]);
 513       grDelta[i]->SetMarkerStyle(25+i);
 514       grDelta[i]->GetXaxis()->SetTitle("sector");
 515       grDelta[i]->GetYaxis()->SetTitle("#Delta_{z} (mm)");
 516       if (itype==2 && i>0) continue;
 517       grDelta[i]->Draw(gropt[i]);
 518       legend->AddEntry(grDelta[i]);
 519     }
 520     legend->Draw();
 521   }
 522
 523   //
 524   // Dtheta
 525   //
 526   canvasDtheta->cd(1);
 527   treeDelta->Draw("dThetameas:dThetafit");
 528   for (Int_t itype=0; itype<3; itype++){
 529     canvasDtheta->cd(itype+2);
 530     entries=treeDelta->Draw(varTypeT[itype],"is1==is0+36","goff");
 531     grDelta[0]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
 532     entries=treeDelta->Draw(varTypeT[itype],"is1==is0+35","goff");
 533     grDelta[1]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
 534     entries=treeDelta->Draw(varTypeT[itype],"is1==is0+37","goff");
 535     grDelta[2]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
 536     legend = new TLegend(0.5,0.7,0.9,0.9, varLegendT[itype]);
 537     for (Int_t i=0; i<3; i++) {
 538       grDelta[i]->SetMaximum(4.);
 539       grDelta[i]->SetMinimum(-3.);
 540       grDelta[i]->SetTitle(type[i]);
 541       grDelta[i]->SetMarkerColor(grcol[i]);
 542       grDelta[i]->SetLineColor(grcol[i]);
 543       grDelta[i]->SetMarkerStyle(25+i);
 544       grDelta[i]->GetXaxis()->SetTitle("sector");
 545       grDelta[i]->GetYaxis()->SetTitle("#Delta_{#theta} (mrad)");
 546       if (itype==2 && i>0) continue;
 547       grDelta[i]->Draw(gropt[i]);
 548       legend->AddEntry(grDelta[i]);
 549     }
 550     legend->Draw();
 551   }
 552
 553   //
 554   // Dphi
 555   //
 556   canvasDphi->cd();
 557   canvasDphi->cd(1);
 558   treeDelta->Draw("dPhimeas:dPhifit");
 559   for (Int_t itype=0; itype<3; itype++){
 560     canvasDphi->cd(itype+2);
 561     entries=treeDelta->Draw(varTypeP[itype],"is1==is0+36","goff");
 562     grDelta[0]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
 563     entries=treeDelta->Draw(varTypeP[itype],"is1==is0+35","goff");
 564     grDelta[1]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
 565     entries=treeDelta->Draw(varTypeP[itype],"is1==is0+37","goff");
 566     grDelta[2]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
 567     legend = new TLegend(0.5,0.7,0.9,0.9, varLegendP[itype]);
 568     for (Int_t i=0; i<3; i++) {
 569       grDelta[i]->SetMaximum(4.);
 570       grDelta[i]->SetMinimum(-3.);
 571       grDelta[i]->SetTitle(type[i]);
 572       grDelta[i]->SetMarkerColor(grcol[i]);
 573       grDelta[i]->SetLineColor(grcol[i]);
 574       grDelta[i]->SetMarkerStyle(25+i);
 575       grDelta[i]->GetXaxis()->SetTitle("sector");
 576       grDelta[i]->GetYaxis()->SetTitle("#Delta_{#phi} (mrad)");
 577       if (itype==2 && i>0) continue;
 578       grDelta[i]->Draw(gropt[i]);
 579       legend->AddEntry(grDelta[i]);
 580     }
 581     legend->Draw();
 582   }
 583   //
 584   //
 585   f.cd();
 586   canvasDy->Write();
 587   canvasDz->Write();
 588   canvasDtheta->Write();
 589   canvasDphi->Write();
 590   //
 591   //
 592   //
 593   TPostScript *ps = new TPostScript("alignReport.ps", 112);
 594   ps->NewPage();
 595   canvasDy->Draw();
 596   ps->NewPage();
 597   canvasDz->Draw();
 598   ps->NewPage();
 599   canvasDtheta->Draw();
 600   ps->NewPage();
 601   canvasDphi->Draw();
 602   ps->Close();
 603   delete ps;
 604 }
 605
 606
 607
 608 void AliTPCkalmanAlign::DumpOldAlignment(TTreeSRedirector *pcstream){
 609   // Dump the content of old alignemnt
 610   // Expected that the old alignmnet is loaded
 611   //
 612   if (!fOriginalAlign) return;
 613   //
 614   TVectorD localTrans(3);
 615   TVectorD globalTrans(3);
 616   TVectorD localRot(3);
 617   TVectorD globalRot(3);
 618   AliGeomManager::ELayerID idLayer;
 619   Int_t idModule=0;
 620   //
 621   for (Int_t i=0; i<fOriginalAlign->GetEntries();i++){
 622     AliAlignObjParams *params = (AliAlignObjParams*)fOriginalAlign->At(i);
 623     params->GetVolUID(idLayer,idModule);
 624     params->GetLocalTranslation(localTrans.GetMatrixArray());
 625     params->GetLocalAngles(localRot.GetMatrixArray());
 626     params->GetTranslation(globalTrans.GetMatrixArray());
 627     params->GetAngles(globalRot.GetMatrixArray());
 628     Int_t sector=idModule;
 629     if (idLayer>7) sector+=36;
 630     (*pcstream)<<"oldAlign"<<
 631       //"idLayer="<<idLayer<<
 632       "idModule="<<idModule<<
 633       "sector="<<sector<<
 634       "lT.="<<&localTrans<<
 635       "gT.="<<&localTrans<<
 636       "lR.="<<&localRot<<
 637       "gR.="<<&globalRot<<
 638       "\n";
 639   }
 640 }
 641
 642
 643 void AliTPCkalmanAlign::MakeNewAlignment(Bool_t badd, TTreeSRedirector * pcstream){
 644   //
 645   // make a new Alignment entry
 646   //
 647   if (!fOriginalAlign) return;
 648   //
 649   TVectorD localTrans(3);
 650   TVectorD globalTrans(3);
 651   TVectorD localRot(3);
 652   TVectorD globalRot(3);
 653   //
 654   TVectorD localTransNew(3);   // new entries
 655   TVectorD globalTransNew(3);
 656   TVectorD localRotNew(3);
 657   TVectorD globalRotNew(3);
 658   //
 659   AliGeomManager::ELayerID idLayer;
 660   Int_t idModule=0;
 661   //
 662   fNewAlign = (TClonesArray*)fOriginalAlign->Clone();
 663   for (Int_t i=0; i<fOriginalAlign->GetEntries();i++){
 664     AliAlignObjParams *params = (AliAlignObjParams*)fOriginalAlign->At(i);
 665     //AliAlignObjParams *paramsNew = (AliAlignObjParams*)fNewAlign->At(i);
 666     params->GetVolUID(idLayer,idModule);
 667     Int_t sector=(Int_t)idModule;
 668     if (idLayer>7) sector+=36;
 669     params->GetLocalTranslation(localTrans.GetMatrixArray());
 670     params->GetLocalAngles(localRot.GetMatrixArray());
 671     params->GetTranslation(globalTrans.GetMatrixArray());
 672     params->GetAngles(globalRot.GetMatrixArray());
 673     //
 674     //
 675     //
 676     if (badd){ // addition if
 677       localTransNew=localTrans;
 678       localRotNew=localRot;
 679     }
 680     localTransNew[1]=localTransNew[1]-((*fDelta1D[0])(sector,0));
 681     localRot[0]  =localRot[0]-(*fDelta1D[2])(sector,0);
 682     //
 683     if (pcstream) (*pcstream)<<"alignParams"<<
 684       //"idLayer="<<idLayer<<
 685       "idModule="<<idModule<<
 686       "sector="<<sector<<
 687       "olT.="<<&localTrans<<
 688       "olR.="<<&localRot<<
 689       "ogT.="<<&localTrans<<
 690       "ogR.="<<&globalRot<<
 691       "nlT.="<<&localTransNew<<
 692       "nlR.="<<&localRotNew<<
 693       "ngT.="<<&localTransNew<<
 694       "ngR.="<<&globalRotNew<<
 695       "\n";
 696   }
 697 }
 698
 699
 700
 701 void AliTPCkalmanAlign::DrawAlignmentTrends(){
 702   //
 703   // Draw trends of alingment variables
 704   //
 705   /*
 706     //1.  Create a list of  align data
 707     //
 708     //2. Filter list
 709     AliXRDPROOFtoolkit::FilterList("align.list","AliTPCkalmanAlign.root kalmanAlignDebug",0);
 710
 711   */
 712   AliXRDPROOFtoolkit toolkit;
 713   TChain * chain = toolkit.MakeChainRandom("align.list.Good","kalmanAlignDebug",0,2000);
 714   TChain * chainRef = toolkit.MakeChainRandom("alignRef.list","kalmanAlignDebug",0,2000);
 715   chain->AddFriend(chainRef,"R");
 716   chainRef->AddFriend(chainRef,"T");
 717   //cuts
 718   TCut cutS="stat.fElements[0]>200&&stat.fElements[1]>200&&stat.fElements[3]>200&&stat.fElements[3]>200";   //statistic in the bin
 719   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
 720   //  TTree *tree  = chain->CopyTree(cutS);
 721   //TTree *treeR = chainRef->CopyTree(cutST);
 722
 723   TCanvas * canvasDy= new TCanvas("canvasDy","canvasDy");
 724   TH1 *his=0;
 725   TLegend *legend = 0;
 726   //  Int_t grcol[3]={2,1,4};
 727
 728   legend = new TLegend(0.7,0.6,0.9,0.9, "Alignment #Delta_{y}- Up-Down");
 729   for (Int_t isec=0; isec<18; isec+=2){
 730     chain->SetMarkerColor(1+(isec%5));
 731     chain->SetMarkerStyle(isec+20);
 732     chain->Draw("10*(delta.fElements[0]-R.delta.fElements[0]):run",cutS+Form("is1==is0+36&&is0==%d",isec),"profgoff");
 733     his = (TH1*)(chain->GetHistogram()->Clone());
 734     his->SetName(Form("#Delta_{Y} sector %d",isec));
 735     his->SetTitle(Form("#Delta_{Y} sector %d",isec));
 736     his->SetMaximum(1.);
 737     his->SetMinimum(-1.);
 738     his->GetYaxis()->SetTitle("#Delta_{y} (mm)");
 739     his->GetXaxis()->SetTitle("run Number");
 740     if (isec==0) his->Draw("");
 741     if (isec>0) his->Draw("same");
 742     legend->AddEntry(his);
 743   }
 744   legend->Draw();
 745   canvasDy->Draw();
 746 }
 747
 748
 749
 750
 751
 752
 753 void AliTPCkalmanAlign::FitCE(){
 754   //
 755   // fit CE
 756   // 1. Global fit - gy and gx
 757   // 2. Local X fit common
 758   // 3. Sector fit
 759   //
 760   AliTPCPreprocessorOnline * preprocesor = new AliTPCPreprocessorOnline;
 761   //
 762   AliTPCCalPad *padTime0 = AliTPCcalibDB::Instance()->GetPadTime0();
 763   AliTPCCalPad *padNoise = AliTPCcalibDB::Instance()->GetPadNoise();
 764   AliTPCCalPad * ceTmean = AliTPCcalibDB::Instance()->GetCETmean();  // CE information
 765   AliTPCCalPad * ceTrms  = AliTPCcalibDB::Instance()->GetCETrms();
 766   AliTPCCalPad * ceQmean = AliTPCcalibDB::Instance()->GetCEQmean();
 767   AliTPCCalPad * pulserTmean = AliTPCcalibDB::Instance()->GetPulserTmean(); //
 768   AliTPCCalPad * pulserTrms = AliTPCcalibDB::Instance()->GetPulserTrms();
 769   AliTPCCalPad * pulserQmean = AliTPCcalibDB::Instance()->GetPulserQmean();
 770   AliTPCCalPad * dmap0  = AliTPCcalibDB::Instance()->GetDistortionMap(0);   // distortion maps
 771   AliTPCCalPad * dmap1  = AliTPCcalibDB::Instance()->GetDistortionMap(1);
 772   AliTPCCalPad * dmap2  = AliTPCcalibDB::Instance()->GetDistortionMap(2);
 773   pulserTmean->Add(-pulserTmean->GetMean());
 774   //
 775   preprocesor->AddComponent(padTime0->Clone());
 776   preprocesor->AddComponent(padNoise->Clone());
 777   preprocesor->AddComponent(pulserTmean->Clone());
 778   preprocesor->AddComponent(pulserQmean->Clone());
 779   preprocesor->AddComponent(pulserTrms->Clone());
 780   preprocesor->AddComponent(ceTmean->Clone());
 781   preprocesor->AddComponent(ceQmean->Clone());
 782   preprocesor->AddComponent(ceTrms->Clone());
 783   preprocesor->AddComponent(dmap0->Clone());
 784   preprocesor->AddComponent(dmap1->Clone());
 785   preprocesor->AddComponent(dmap2->Clone());
 786   preprocesor->DumpToFile("cetmean.root");
 787
 788   TCut cutNoise="abs(PadNoise.fElements/PadNoise_Median-1)<0.3";
 789   TCut cutPulserT="abs(PulserTrms.fElements/PulserTrms_Median-1)<0.2";
 790   TCut cutPulserQ="abs(PulserQmean.fElements/PulserQmean_Median-1)<0.2";
 791   TCut cutCEQ="CEQmean.fElements>50";
 792   TCut cutCET="abs(CETmean.fElements)<2";
 793   TCut cutAll=cutNoise+cutPulserT+cutPulserQ+cutCEQ+cutCET;
 794   //
 795   //
 796   TFile * f = new TFile("cetmean.root");
 797   TTree * chain = (TTree*) f->Get("calPads");
 798   Int_t entries = chain->Draw("1",cutAll,"goff");
 799   if (entries<200000) return;  // no calibration available - pulser or CE or noise
 800
 801   TStatToolkit toolkit;
 802   Double_t chi2=0;
 803   Int_t    npoints=0;
 804   TVectorD param;
 805   TMatrixD covar;
 806   //
 807   // make a aliases
 808   AliTPCkalmanAlign::MakeAliasCE(chain);
 809   TString  fstringG="";              // global part
 810   //
 811   fstringG+="Gy++";                  // 1 - global y
 812   fstringG+="Gx++";                  // 2 - global x
 813   //
 814   fstringG+="isin++";                // 3 -delta IROC-OROC offset
 815   fstringG+="Lx++";                  // 4 -common slope
 816   fstringG+="Lx*isin++";             // 5 -delta slope
 817   fstringG+="Ly++";                  // 6 -common slope
 818   fstringG+="Ly*isin++";             // 7 -delta slope
 819   TVectorD vecG[2];
 820   TString * strFitG=0;
 821   TString * strFitLX=0;
 822   //
 823   strFitG = TStatToolkit::FitPlane(chain,"deltaT", fstringG.Data(),"sideA"+cutAll, chi2,npoints,vecG[0],covar,-1,0, 10000000, kFALSE);
 824   chain->SetAlias("tfitGA",strFitG->Data());
 825   strFitG->Tokenize("++")->Print();
 826   printf("chi2=%f\n",TMath::Sqrt(chi2/npoints));
 827   //
 828   strFitG = TStatToolkit::FitPlane(chain,"deltaT", fstringG.Data(),"sideC"+cutAll, chi2,npoints,vecG[1],covar,-1,0, 10000000, kFALSE);
 829   chain->SetAlias("tfitGC",strFitG->Data());
 830   strFitG->Tokenize("++")->Print();
 831   printf("chi2=%f\n",TMath::Sqrt(chi2/npoints));
 832   //
 833   AliTPCCalPad *padFitG =AliTPCCalPad::CreateCalPadFit("1++gy/500.++gx/500.++0+++0++0++0++0",vecG[0],vecG[1]);
 834   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]);
 835   // swap a side and c side
 836   AliTPCCalPad *padFitGSwap =AliTPCCalPad::CreateCalPadFit("1++gy/500.++gx/500.++0+++0++0++0++0",vecG[1],vecG[0]);
 837   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]);
 838   padFitG->SetName("CEG");
 839   padFitLX->SetName("CELX");
 840   padFitGSwap->SetName("CEGS");
 841   padFitLXSwap->SetName("CELXS");
 842   preprocesor->AddComponent(padFitG->Clone());
 843   preprocesor->AddComponent(padFitLX->Clone());
 844   preprocesor->AddComponent(padFitGSwap->Clone());
 845   preprocesor->AddComponent(padFitLXSwap->Clone());
 846   preprocesor->DumpToFile("cetmean.root");   // add it to the file
 847   //
 848   // make local fits
 849   //
 850   f = new TFile("cetmean.root");
 851   chain = (TTree*) f->Get("calPads");
 852   AliTPCkalmanAlign::MakeAliasCE(chain);
 853   TString  fstringL="";              // local fit
 854   //                                 // 0. delta common
 855   fstringL+="isin++";                // 1. delta IROC-OROC offset
 856   fstringL+="Lx++";                  // 2. common slope
 857   fstringL+="Lx*isin++";             // 3. delta slope
 858   fstringL+="Ly++";                  // 4. common slope
 859   fstringL+="Ly*isin++";             // 5. delta slope
 860   TVectorD vecL[36];
 861   TVectorD dummy(6);
 862   AliTPCCalPad *padFitLCE = new AliTPCCalPad("LocalCE","LocalCE");
 863   AliTPCCalPad *padFitTmpCE;
 864   for (Int_t isec=0; isec<36; isec++){
 865     TCut cutSector=Form("(sector%36)==%d",isec);
 866     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);
 867     printf("sec=%d\tchi2=%f\n",isec,TMath::Sqrt(chi2/npoints));
 868     //
 869     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);
 870     if (isec<18) padFitTmpCE=AliTPCCalPad::CreateCalPadFit(fitL.Data(),vecL[isec],dummy);
 871     if (isec>=18) padFitTmpCE=AliTPCCalPad::CreateCalPadFit(fitL.Data(),dummy,vecL[isec]);
 872     padFitLCE->Add(padFitTmpCE);
 873   }
 874   //
 875   padFitLCE->SetName("CELocal");
 876   preprocesor->AddComponent(padFitLCE->Clone());
 877   preprocesor->DumpToFile("cetmean.root");   // add it to the file
 878   //
 879   // write data to array
 880   //
 881   fFitCELocal  = new TObjArray(6);
 882   fFitCEGlobal = new TObjArray(8);
 883   for (Int_t ipar=0; ipar<8;ipar++){
 884     //
 885     fFitCEGlobal->AddAt(new TVectorD(36),ipar);
 886     TVectorD &fvecG = *((TVectorD*)fFitCEGlobal->At(ipar));
 887     for (Int_t isec=0; isec<36;isec++){
 888       if (isec<18)  fvecG[isec]=vecG[0][ipar];
 889       if (isec>=18) fvecG[isec]=vecG[1][ipar];
 890     }
 891     if (ipar<6){
 892       fFitCELocal->AddAt(new TVectorD(36),ipar);
 893       TVectorD &fvecL = *((TVectorD*)fFitCELocal->At(ipar));
 894       for (Int_t isec=0; isec<36;isec++){
 895         fvecL[isec]=vecL[isec][ipar];
 896       }
 897     }
 898   }
 899 }
 900
 901 void AliTPCkalmanAlign::MakeAliasCE(TTree * chain){
 902   //
 903   // make a aliases of pad variables
 904   //
 905   chain->SetAlias("side","(-1+(sector%36<18)*2)");
 906   chain->SetAlias("sideA","(sector%36<18)");
 907   chain->SetAlias("sideC","(sector%36>=18)");
 908   chain->SetAlias("isin","(sector<36)");
 909   chain->SetAlias("deltaT","CETmean.fElements-PulserTmean.fElements");
 910   chain->SetAlias("timeP","PulserTmean.fElements");
 911   chain->SetAlias("Gy","(gy.fElements/500.)");
 912   chain->SetAlias("Gx","(gx.fElements/500.)");
 913   chain->SetAlias("Lx","(lx.fElements-133)/100.");   // lx in meters
 914   chain->SetAlias("Ly","(ly.fElements)/100.");
 915   chain->SetAlias("La","(ly.fElements/lx.fElements/0.155)");
 916   chain->SetAlias("deltaT","(CETmean.fElements-PulserTmean.fElements)");
 917 }