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=117092;
  35   .x ConfigCalibTrain.C(run)
  36
  37   AliTPCkalmanAlign kalmanAlign("TPC align", "TPC align");  // create the object
  38   kalmanAlign.ReadAlign("CalibObjects.root");               // read the calibration form file
  39   kalmanAlign.MakeGlobalAlign();                            // do kalman alignment
  40   kalmanAlign.DrawDeltaAlign();                             // make QA plot
  41   //
  42
  43
  44 */
  45 #include "TMath.h"
  46 #include "TTreeStream.h"
  47 #include "TGraph.h"
  48 #include "TGraphErrors.h"
  49 #include "TVectorD.h"
  50 #include "TClonesArray.h"
  51
  52 #include "AliTPCcalibDB.h"
  53 #include "AliTPCCalROC.h"
  54 #include "AliCDBMetaData.h"
  55 #include "AliCDBId.h"
  56 #include "AliCDBManager.h"
  57 #include "AliCDBStorage.h"
  58 #include "AliCDBEntry.h"
  59 #include "AliAlignObjParams.h"
  60 #include "AliTPCROC.h"
  61 #include "AliTracker.h"
  62 #include "TFile.h"
  63 #include "TLinearFitter.h"
  64 #include "AliTPCcalibAlign.h"
  65 #include "TH1.h"
  66 #include "AliTPCCalPad.h"
  67 #include "AliTPCkalmanAlign.h"
  68 #include "TLegend.h"
  69 #include "TCanvas.h"
  70
  71 AliTPCkalmanAlign::AliTPCkalmanAlign():
  72   TNamed(),
  73   fCalibAlign(0),     // kalman alignnmnt
  74   fOriginalAlign(0),   // original alignment 0 read for the OCDB
  75   fNewAlign(0)
  76 {
  77   //
  78   // Default constructor
  79   //
  80   for (Int_t i=0; i<4; i++){
  81     fDelta1D[i]=0;
  82     fCovar1D[i]=0;
  83   }
  84 }
  85
  86 AliTPCkalmanAlign::AliTPCkalmanAlign(const char* name, const char* title):
  87   TNamed(name,title),
  88   fCalibAlign(0),     // kalman alignnmnt
  89   fOriginalAlign(0),   // original alignment 0 read for the OCDB
  90   fNewAlign(0)     // kalman alignnmnt
  91 {
  92   //
  93   // Default constructor
  94   //
  95   for (Int_t i=0; i<4; i++){
  96     fDelta1D[i]=0;
  97     fCovar1D[i]=0;
  98   }
  99 }
 100
 101 void AliTPCkalmanAlign::ReadAlign(const char *fname){
 102   //
 103   // Read old alignment used in the reco
 104   // and the residual histograms
 105   // WE ASSUME that the OCDB path is set in the same way as done in the calibration
 106   //
 107   TFile fcalib(fname);
 108   TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
 109   fCalibAlign=0;
 110   if (array) fCalibAlign=( AliTPCcalibAlign *)array->FindObject("alignTPC");
 111   fCalibAlign = ( AliTPCcalibAlign *)fcalib.Get("alignTPC");
 112   //
 113   // old alignment used
 114   AliCDBEntry * cdbEntry= AliCDBManager::Instance()->Get("TPC/Align/Data");
 115   fOriginalAlign =0;
 116   if (cdbEntry){
 117     fOriginalAlign = (TClonesArray*)cdbEntry->GetObject();
 118   }
 119
 120 }
 121
 122 void AliTPCkalmanAlign::BookAlign1D(TMatrixD &param, TMatrixD &covar, Double_t mean, Double_t sigma){
 123   //
 124   // Book Align 1D parameters and covariance
 125   //
 126   param.ResizeTo(72,1);
 127   covar.ResizeTo(72,72);
 128   for (Int_t i=0;i<72;i++){
 129     param(i,0)=mean;
 130     for (Int_t j=0;j<72;j++) covar(i,j)=0;
 131     covar(i,i)=sigma*sigma;
 132   }
 133 }
 134
 135
 136 void AliTPCkalmanAlign::UpdateAlign1D(Double_t delta, Double_t sigma, Int_t s1, Int_t s2, TMatrixD &vecXk, TMatrixD &covXk){
 137   //
 138   // Update 1D kalman filter
 139   //
 140   const Int_t knMeas=1;
 141   const Int_t knElem=72;
 142   TMatrixD mat1(72,72);            // update covariance matrix
 143   TMatrixD matHk(1,knElem);        // vector to mesurement
 144   TMatrixD vecYk(knMeas,1);        // Innovation or measurement residual
 145   TMatrixD matHkT(knElem,knMeas);  // helper matrix Hk transpose
 146   TMatrixD matSk(knMeas,knMeas);   // Innovation (or residual) covariance
 147   TMatrixD matKk(knElem,knMeas);   // Optimal Kalman gain
 148   TMatrixD covXk2(knElem,knElem);  // helper matrix
 149   TMatrixD covXk3(knElem,knElem);  // helper matrix
 150   TMatrixD vecZk(1,1);
 151   TMatrixD measR(1,1);
 152   vecZk(0,0)=delta;
 153   measR(0,0)=sigma*sigma;
 154   //
 155   // reset matHk
 156   for (Int_t iel=0;iel<knElem;iel++)
 157     for (Int_t ip=0;ip<knMeas;ip++) matHk(ip,iel)=0;
 158   //mat1
 159   for (Int_t iel=0;iel<knElem;iel++) {
 160     for (Int_t jel=0;jel<knElem;jel++) mat1(iel,jel)=0;
 161     mat1(iel,iel)=1;
 162   }
 163   //
 164   matHk(0, s1)=1;
 165   matHk(0, s2)=-1;
 166   vecYk = vecZk-matHk*vecXk;               // Innovation or measurement residual
 167   matHkT=matHk.T(); matHk.T();
 168   matSk = (matHk*(covXk*matHkT))+measR;    // Innovation (or residual) covariance
 169   matSk.Invert();
 170   matKk = (covXk*matHkT)*matSk;            //  Optimal Kalman gain
 171   vecXk += matKk*vecYk;                    //  updated vector
 172   covXk2= (mat1-(matKk*matHk));
 173   covXk3 =  covXk2*covXk;
 174   covXk = covXk3;
 175 }
 176
 177
 178
 179
 180 void AliTPCkalmanAlign::MakeGlobalAlign(){
 181   //
 182   // Combine all pairs of fitters and make global alignemnt
 183   //
 184   AliTPCkalmanAlign &kalmanAlign=*this;
 185   TTreeSRedirector *pcstream = new TTreeSRedirector("AliTPCkalmanAlign.root");
 186   //
 187   // get ce info
 188   //
 189   AliTPCCalPad *padTime0 = AliTPCcalibDB::Instance()->GetPadTime0();
 190   TVectorD paramCE[72];
 191   TMatrixD covarCE[72];
 192   Float_t chi2;
 193   for (Int_t isec=0; isec<72; isec++){
 194     AliTPCCalROC * roc0  = padTime0->GetCalROC(isec);
 195     roc0->GlobalFit(0,kFALSE,paramCE[isec],covarCE[isec],chi2,0);
 196     (*pcstream)<<"ceFit"<<
 197       "isec="<<isec<<
 198       "p0.="<<&paramCE[isec]<<
 199       "\n";
 200   }
 201
 202   DumpOldAlignment(pcstream);
 203   const Int_t kMinEntries=400;
 204   TMatrixD vec[5];
 205   TMatrixD cov[5];
 206   kalmanAlign.BookAlign1D(vec[0],cov[0], 0,0.5);
 207   kalmanAlign.BookAlign1D(vec[1],cov[1], 0,0.5);
 208   kalmanAlign.BookAlign1D(vec[2],cov[2], 0,0.01);
 209   kalmanAlign.BookAlign1D(vec[3],cov[3], 0,0.01);
 210   //
 211   TVectorD delta(5);
 212   TVectorD rms(5);
 213   TVectorD stat(5);
 214   TH1 * his=0;
 215   for (Int_t is0=0;is0<72;is0++)
 216     for (Int_t is1=0;is1<72;is1++){
 217       //
 218       //
 219       his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kY,is0,is1);
 220       if (!his) continue;
 221       if (his->GetEntries()<kMinEntries) continue;
 222       delta[0]=his->GetMean();
 223       rms[0]=his->GetRMS();
 224       stat[0]=his->GetEntries();
 225       kalmanAlign.UpdateAlign1D(his->GetMean(),his->GetRMS(),is0,is1, vec[0],cov[0]);
 226       //
 227       his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kZ,is0,is1);
 228       if (!his) continue;
 229       delta[1]=his->GetMean();
 230       rms[1]=his->GetRMS();
 231       stat[1]=his->GetEntries();
 232       kalmanAlign.UpdateAlign1D(his->GetMean(),his->GetRMS(),is0,is1, vec[1],cov[1]);
 233       //
 234       his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kPhi,is0,is1);
 235       if (!his) continue;
 236       delta[2] = his->GetMean();
 237       rms[2]=his->GetRMS();
 238       stat[2]=his->GetEntries();
 239       kalmanAlign.UpdateAlign1D(his->GetMean(),his->GetRMS(),is0,is1, vec[2],cov[2]);
 240       //
 241       his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kTheta,is0,is1);
 242       if (!his) continue;
 243       delta[3] = his->GetMean();
 244       rms[3]=his->GetRMS();
 245       stat[3]=his->GetEntries();
 246       kalmanAlign.UpdateAlign1D(his->GetMean(),his->GetRMS(),is0,is1, vec[3],cov[3]);
 247     }
 248
 249   fDelta1D[0] = (TMatrixD*)vec[0].Clone();
 250   fDelta1D[1] = (TMatrixD*)vec[1].Clone();
 251   fDelta1D[2] = (TMatrixD*)vec[2].Clone();
 252   fDelta1D[3] = (TMatrixD*)vec[3].Clone();
 253   //
 254   fCovar1D[0] = (TMatrixD*)cov[0].Clone();
 255   fCovar1D[1] = (TMatrixD*)cov[1].Clone();
 256   fCovar1D[2] = (TMatrixD*)cov[2].Clone();
 257   fCovar1D[3] = (TMatrixD*)cov[3].Clone();
 258
 259   MakeNewAlignment(kTRUE);
 260
 261   for (Int_t is0=0;is0<72;is0++)
 262     for (Int_t is1=0;is1<72;is1++){
 263       Bool_t isPair=kFALSE;
 264       if (TMath::Abs(is0%18-is1%18)<2) isPair=kTRUE;
 265       if (TMath::Abs(is0%18-is1%18)==17) isPair=kTRUE;
 266       if (!isPair) continue;
 267       stat[0]=0; stat[1]=0; stat[2]=0; stat[3]=0;
 268       //
 269       //
 270       his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kY,is0,is1);
 271       if (his){
 272         delta[0]=his->GetMean();
 273         rms[0]=his->GetRMS();
 274         stat[0]=his->GetEntries();
 275       }
 276       //
 277       his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kZ,is0,is1);
 278       if (his) {
 279         delta[1]=his->GetMean();
 280         rms[1]=his->GetRMS();
 281         stat[1]=his->GetEntries();
 282       }
 283       //
 284       his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kPhi,is0,is1);
 285       if (his){
 286         delta[2] = his->GetMean();
 287         rms[2]=his->GetRMS();
 288         stat[2]=his->GetEntries();
 289       }
 290       //
 291       his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kTheta,is0,is1);
 292       if (his){
 293         delta[3] = his->GetMean();
 294         rms[3]=his->GetRMS();
 295         stat[3]=his->GetEntries();
 296       }
 297       Int_t run = AliCDBManager::Instance()->GetRun();
 298       Float_t bz = AliTracker::GetBz();
 299       (*pcstream)<<"kalmanAlignDebug"<<
 300         "run="<<run<<
 301         "bz="<<bz<<
 302         "is0="<<is0<<
 303         "is1="<<is1<<
 304         "delta.="<<&delta<<
 305         "rms.="<<&rms<<
 306         "stat.="<<&stat<<
 307         "vec0.="<<&vec[0]<<
 308         "vec1.="<<&vec[1]<<
 309         "vec2.="<<&vec[2]<<
 310         "vec3.="<<&vec[3]<<
 311         "pceIn0.="<<&paramCE[is0%36]<<
 312         "pceOut0.="<<&paramCE[is0%36+36]<<
 313         "pceIn1.="<<&paramCE[is1%36]<<
 314         "pceOut1.="<<&paramCE[is1%36+36]<<
 315         "\n";
 316     }
 317   delete pcstream;
 318 }
 319
 320
 321
 322
 323
 324
 325 void AliTPCkalmanAlign::UpdateOCDBTime0( AliTPCCalPad  *pad, Int_t ustartRun, Int_t uendRun,  const char* storagePath ){
 326   //
 327   // Update OCDB
 328   // .x ConfigCalibTrain.C(117116)
 329   // AliTPCcalibDB::Instance()->GetPulserTmean()
 330   // pad->Add(-pad->GetMean())
 331   AliCDBMetaData *metaData= new AliCDBMetaData();
 332   metaData->SetObjectClassName("TObjArray");
 333   metaData->SetResponsible("Marian Ivanov");
 334   metaData->SetBeamPeriod(1);
 335   metaData->SetAliRootVersion("05-25-01"); //root version
 336   metaData->SetComment("Calibration of the z - time misalignment");
 337   AliCDBId* id1=NULL;
 338   id1=new AliCDBId("TPC/Calib/PadTime0", ustartRun, uendRun);
 339   AliCDBStorage* gStorage = AliCDBManager::Instance()->GetStorage(storagePath);
 340   gStorage->Put(pad, (*id1), metaData);
 341 }
 342
 343
 344
 345 void AliTPCkalmanAlign::DrawDeltaAlign(){
 346   //
 347   // Draw the reuslts of the alignment
 348   // Residual misalignment in respect with previous alignment shown
 349   //
 350   //
 351   TFile f("AliTPCkalmanAlign.root","update");
 352   TTree * treeDelta = (TTree*)f.Get("kalmanAlignDebug");
 353   TH1::AddDirectory(0);
 354   //
 355   treeDelta->SetAlias("sector","is0");
 356   treeDelta->SetAlias("dYmeas","10*(delta.fElements[0])");
 357   treeDelta->SetAlias("dZmeas","10*(delta.fElements[1])");
 358   treeDelta->SetAlias("dPhimeas","1000*(delta.fElements[2])");
 359   treeDelta->SetAlias("dThetameas","1000*(delta.fElements[3])");
 360   //
 361   treeDelta->SetAlias("dYfit","10*(vec0.fElements[is0]-vec0.fElements[is1])");
 362   treeDelta->SetAlias("dZfit","10*(vec1.fElements[is0]-vec1.fElements[is1])");
 363   treeDelta->SetAlias("dPhifit","1000*(vec2.fElements[is0]-vec2.fElements[is1])");
 364   treeDelta->SetAlias("dThetafit","1000*(vec3.fElements[is0]-vec3.fElements[is1])");
 365   //
 366   treeDelta->SetMarkerStyle(25);
 367   treeDelta->SetMarkerColor(4);
 368   treeDelta->SetLineColor(4);
 369   const char *type[3]={"up-down","left-right","right-left"};
 370   const char *gropt[3]={"alp","lp","lp"};
 371   const char *varTypeY[3]={"dYmeas-dYfit:sector","dYfit:sector","10*vec0.fElements[is0]:sector"};
 372   const char *varLegendY[3]={"#Delta_{y} fit residual","#Delta_{y} fit value difference","#Delta_{y} sector"};
 373   const char *varTypeZ[3]={"dZmeas-dZfit:sector","dZfit:sector","10*vec1.fElements[is0]:sector"};
 374   const char *varLegendZ[3]={"#Delta_{Z} fit residual","#Delta_{Z} fit value difference","#Delta_{Z} sector"};
 375   const char *varTypeT[3]={"dThetameas-dThetafit:sector","dThetafit:sector","1000*vec3.fElements[is0]:sector"};
 376   const char *varLegendT[3]={"#Delta_{#theta} fit residual","#Delta_{#theta} fit value difference","#Delta_{#theta} sector"};
 377   const char *varTypeP[3]={"dPhimeas-dPhifit:sector","dPhifit:sector","1000*vec2.fElements[is0]:sector"};
 378   const char *varLegendP[3]={"#Delta_{#phi} fit residual","#Delta_{#phi} fit value difference","#Delta_{#phi} sector"};
 379   TLegend *legend = 0;
 380   Int_t grcol[3]={2,1,4};
 381   Int_t entries=0;
 382   TGraph *grDelta[3]={0,0,0};
 383   TCanvas * canvasDy=new TCanvas("canvasDy","canvasDy",1200,800);
 384   TCanvas * canvasDz=new TCanvas("canvasDz","canvasDz",1200,800);
 385   TCanvas * canvasDphi=new TCanvas("canvasDphi","canvasDphi",1200,800);
 386   TCanvas * canvasDtheta=new TCanvas("canvasDtheta","canvasDtheta",1200,800);
 387   canvasDy->Divide(2,2);
 388   canvasDz->Divide(2,2);
 389   canvasDtheta->Divide(2,2);
 390   canvasDphi->Divide(2,2);
 391
 392   //
 393   // Dy
 394   //
 395   canvasDy->cd(1);
 396   treeDelta->Draw("dYmeas:dYfit");
 397   for (Int_t itype=0; itype<3; itype++){
 398     canvasDy->cd(itype+2);
 399     entries=treeDelta->Draw(varTypeY[itype],"is1==is0+36","goff");
 400     grDelta[0]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
 401     entries=treeDelta->Draw(varTypeY[itype],"is1==is0+35","goff");
 402     grDelta[1]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
 403     entries=treeDelta->Draw(varTypeY[itype],"is1==is0+37","goff");
 404     grDelta[2]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
 405     legend = new TLegend(0.5,0.7,0.9,0.9, varLegendY[itype]);
 406     for (Int_t i=0; i<3; i++) {
 407       grDelta[i]->SetMaximum(1.5);
 408       grDelta[i]->SetMinimum(-1);
 409       grDelta[i]->SetTitle(type[i]);
 410       grDelta[i]->SetMarkerColor(grcol[i]);
 411       grDelta[i]->SetLineColor(grcol[i]);
 412       grDelta[i]->SetMarkerStyle(25+i);
 413       grDelta[i]->GetXaxis()->SetTitle("sector");
 414       grDelta[i]->GetYaxis()->SetTitle("#Delta_{y} (mm)");
 415       if (itype==2 && i>0) continue;
 416       grDelta[i]->Draw(gropt[i]);
 417       legend->AddEntry(grDelta[i]);
 418     }
 419     legend->Draw();
 420   }
 421   //
 422   // Dz
 423   //
 424   canvasDz->cd();
 425   canvasDz->cd(1);
 426   treeDelta->Draw("dZmeas:dZfit");
 427   for (Int_t itype=0; itype<3; itype++){
 428     canvasDz->cd(itype+2);
 429     entries=treeDelta->Draw(varTypeZ[itype],"is1==is0+36","goff");
 430     grDelta[0]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
 431     entries=treeDelta->Draw(varTypeZ[itype],"is1==is0+35","goff");
 432     grDelta[1]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
 433     entries=treeDelta->Draw(varTypeZ[itype],"is1==is0+37","goff");
 434     grDelta[2]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
 435     legend = new TLegend(0.5,0.7,0.9,0.9, varLegendZ[itype]);
 436     for (Int_t i=0; i<3; i++) {
 437       grDelta[i]->SetMaximum(1.5);
 438       grDelta[i]->SetMinimum(-1);
 439       grDelta[i]->SetTitle(type[i]);
 440       grDelta[i]->SetMarkerColor(grcol[i]);
 441       grDelta[i]->SetLineColor(grcol[i]);
 442       grDelta[i]->SetMarkerStyle(25+i);
 443       grDelta[i]->GetXaxis()->SetTitle("sector");
 444       grDelta[i]->GetYaxis()->SetTitle("#Delta_{z} (mm)");
 445       if (itype==2 && i>0) continue;
 446       grDelta[i]->Draw(gropt[i]);
 447       legend->AddEntry(grDelta[i]);
 448     }
 449     legend->Draw();
 450   }
 451
 452   //
 453   // Dtheta
 454   //
 455   canvasDtheta->cd(1);
 456   treeDelta->Draw("dThetameas:dThetafit");
 457   for (Int_t itype=0; itype<3; itype++){
 458     canvasDtheta->cd(itype+2);
 459     entries=treeDelta->Draw(varTypeT[itype],"is1==is0+36","goff");
 460     grDelta[0]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
 461     entries=treeDelta->Draw(varTypeT[itype],"is1==is0+35","goff");
 462     grDelta[1]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
 463     entries=treeDelta->Draw(varTypeT[itype],"is1==is0+37","goff");
 464     grDelta[2]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
 465     legend = new TLegend(0.5,0.7,0.9,0.9, varLegendT[itype]);
 466     for (Int_t i=0; i<3; i++) {
 467       grDelta[i]->SetMaximum(4.);
 468       grDelta[i]->SetMinimum(-3.);
 469       grDelta[i]->SetTitle(type[i]);
 470       grDelta[i]->SetMarkerColor(grcol[i]);
 471       grDelta[i]->SetLineColor(grcol[i]);
 472       grDelta[i]->SetMarkerStyle(25+i);
 473       grDelta[i]->GetXaxis()->SetTitle("sector");
 474       grDelta[i]->GetYaxis()->SetTitle("#Delta_{#theta} (mrad)");
 475       if (itype==2 && i>0) continue;
 476       grDelta[i]->Draw(gropt[i]);
 477       legend->AddEntry(grDelta[i]);
 478     }
 479     legend->Draw();
 480   }
 481
 482   //
 483   // Dphi
 484   //
 485   canvasDphi->cd();
 486   canvasDphi->cd(1);
 487   treeDelta->Draw("dPhimeas:dPhifit");
 488   for (Int_t itype=0; itype<3; itype++){
 489     canvasDphi->cd(itype+2);
 490     entries=treeDelta->Draw(varTypeP[itype],"is1==is0+36","goff");
 491     grDelta[0]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
 492     entries=treeDelta->Draw(varTypeP[itype],"is1==is0+35","goff");
 493     grDelta[1]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
 494     entries=treeDelta->Draw(varTypeP[itype],"is1==is0+37","goff");
 495     grDelta[2]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
 496     legend = new TLegend(0.5,0.7,0.9,0.9, varLegendP[itype]);
 497     for (Int_t i=0; i<3; i++) {
 498       grDelta[i]->SetMaximum(4.);
 499       grDelta[i]->SetMinimum(-3.);
 500       grDelta[i]->SetTitle(type[i]);
 501       grDelta[i]->SetMarkerColor(grcol[i]);
 502       grDelta[i]->SetLineColor(grcol[i]);
 503       grDelta[i]->SetMarkerStyle(25+i);
 504       grDelta[i]->GetXaxis()->SetTitle("sector");
 505       grDelta[i]->GetYaxis()->SetTitle("#Delta_{#phi} (mrad)");
 506       if (itype==2 && i>0) continue;
 507       grDelta[i]->Draw(gropt[i]);
 508       legend->AddEntry(grDelta[i]);
 509     }
 510     legend->Draw();
 511   }
 512   //
 513   //
 514   f.cd();
 515   canvasDy->Write();
 516   canvasDz->Write();
 517   canvasDtheta->Write();
 518   canvasDphi->Write();
 519   //
 520   canvasDy->SaveAs("alignDy.pdf");
 521   canvasDz->SaveAs("alignDz.pdf");
 522   canvasDtheta->SaveAs("alignDtheta.pdf");
 523   canvasDphi->SaveAs("alignDphi.pdf");
 524 }
 525
 526
 527
 528 void AliTPCkalmanAlign::DumpOldAlignment(TTreeSRedirector *pcstream){
 529   // Dump the content of old alignemnt
 530   // Expected that the old alignmnet is loaded
 531   //
 532   if (!fOriginalAlign) return;
 533   //
 534   TVectorD localTrans(3);
 535   TVectorD globalTrans(3);
 536   TVectorD localRot(3);
 537   TVectorD globalRot(3);
 538   AliGeomManager::ELayerID idLayer;
 539   Int_t idModule=0;
 540   //
 541   for (Int_t i=0; i<fOriginalAlign->GetEntries();i++){
 542     AliAlignObjParams *params = (AliAlignObjParams*)fOriginalAlign->At(i);
 543     params->GetVolUID(idLayer,idModule);
 544     params->GetLocalTranslation(localTrans.GetMatrixArray());
 545     params->GetLocalAngles(localRot.GetMatrixArray());
 546     params->GetTranslation(globalTrans.GetMatrixArray());
 547     params->GetAngles(globalRot.GetMatrixArray());
 548     Int_t sector=idModule;
 549     if (idLayer>7) sector+=36;
 550     (*pcstream)<<"oldAlign"<<
 551       //"idLayer="<<idLayer<<
 552       "idModule="<<idModule<<
 553       "sector="<<sector<<
 554       "lT.="<<&localTrans<<
 555       "gT.="<<&localTrans<<
 556       "lR.="<<&localRot<<
 557       "gR.="<<&globalRot<<
 558       "\n";
 559   }
 560 }
 561
 562
 563 void AliTPCkalmanAlign::MakeNewAlignment(Bool_t badd, TTreeSRedirector * pcstream){
 564   //
 565   // make a new Alignment entry
 566   //
 567   if (!fOriginalAlign) return;
 568   //
 569   TVectorD localTrans(3);
 570   TVectorD globalTrans(3);
 571   TVectorD localRot(3);
 572   TVectorD globalRot(3);
 573   //
 574   TVectorD localTransNew(3);   // new entries
 575   TVectorD globalTransNew(3);
 576   TVectorD localRotNew(3);
 577   TVectorD globalRotNew(3);
 578   //
 579   AliGeomManager::ELayerID idLayer;
 580   Int_t idModule=0;
 581   //
 582   fNewAlign = (TClonesArray*)fOriginalAlign->Clone();
 583   for (Int_t i=0; i<fOriginalAlign->GetEntries();i++){
 584     AliAlignObjParams *params = (AliAlignObjParams*)fOriginalAlign->At(i);
 585     //AliAlignObjParams *paramsNew = (AliAlignObjParams*)fNewAlign->At(i);
 586     params->GetVolUID(idLayer,idModule);
 587     Int_t sector=(Int_t)idModule;
 588     if (idLayer>7) sector+=36;
 589     params->GetLocalTranslation(localTrans.GetMatrixArray());
 590     params->GetLocalAngles(localRot.GetMatrixArray());
 591     params->GetTranslation(globalTrans.GetMatrixArray());
 592     params->GetAngles(globalRot.GetMatrixArray());
 593     //
 594     //
 595     //
 596     if (badd){ // addition if
 597       localTransNew=localTrans;
 598       localRotNew=localRot;
 599     }
 600     localTransNew[1]=localTransNew[1]-((*fDelta1D[0])(sector,0));
 601     localRot[0]  =localRot[0]-(*fDelta1D[2])(sector,0);
 602     //
 603     if (pcstream) (*pcstream)<<"alignParams"<<
 604       //"idLayer="<<idLayer<<
 605       "idModule="<<idModule<<
 606       "sector="<<sector<<
 607       "olT.="<<&localTrans<<
 608       "olR.="<<&localRot<<
 609       "ogT.="<<&localTrans<<
 610       "ogR.="<<&globalRot<<
 611       "nlT.="<<&localTransNew<<
 612       "nlR.="<<&localRotNew<<
 613       "ngT.="<<&localTransNew<<
 614       "ngR.="<<&globalRotNew<<
 615       "\n";
 616   }
 617
 618
 619 }