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