TPC/CalibMacros/MakeAlignCalPad.C

   1 /// \file MakeAlignCalPad.C
   2 ///
   3 /// \author marian.ivanov@cern.ch
   4 ///
   5 /// Macro to create  alignment/distortion maps
   6 /// As a input output of AliTPCcalibAlign is used.
   7 ///
   8 /// Algorithm:
   9 ///
  10 /// In the setup without the magnetic field the tracks are fitted using the linear track model.
  11 /// ( in y-rphi coordinate the primary vertex is also used as a constrain)
  12 /// Residuals (deltas0  between the track and clusters in Y and in z direction are filled in the 4 D histogram:
  13 /// Delta: phi(180 bins): localX(53 bins): tan(phi): tan(theta)(10 bins)
  14 ///
  15 /// Distortion map are extracted form the residual histograms as a mean value at each bin.
  16 /// Linear fits are then performed for each pad - delta as function of theta
  17 ///
  18 /// ~~~{.cpp}
  19 /// Delta Ymeas = offsetY+slopeY*tan(theta)
  20 /// Delta Zmeas = offsetZ+slopeZ*tan(theta)
  21 /// ~~~
  22 ///
  23 /// Resulting residuals exported into the OCDB are:
  24 ///
  25 /// ~~~{.cpp}
  26 /// DeltaY  = offsetY
  27 /// DeltaZ  = offsetZ
  28 /// DeltaR  = slopeZ;
  29 /// ~~~
  30 ///
  31 /// Example usage:
  32 ///
  33 /// make calpad+ make report ps file:
  34 ///
  35 /// ~~~
  36 /// aliroot -b -q ~/NimStyle.C ../ConfigCalibTrain.C\(119037\) $ALICE_ROOT/TPC/CalibMacros/MakeAlignCalPad.C\(1\)
  37 /// ~~~
  38 ///
  39 /// make only report ps file:
  40 ///
  41 /// ~~~
  42 /// aliroot -b -q ~/NimStyle.C $ALICE_ROOT/TPC/CalibMacros/MakeAlignCalPad.C\(3\)
  43 /// ~~~
  44 ///
  45 /// Making fit - iterative procedure - see below:
  46 ///
  47 /// ~~~{.cpp}
  48 /// gROOT->Macro("~/rootlogon.C");
  49 /// gSystem->Load("libANALYSIS");
  50 /// gSystem->Load("libSTAT");
  51 /// gSystem->Load("libTPCcalib");
  52 /// gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros -I$ALICE_ROOT/TPC/TPC -I$ALICE_ROOT/STAT");
  53 /// .L $ALICE_ROOT/TPC/CalibMacros/MakeAlignCalPad.C+
  54 /// // load OCDB
  55 ///
  56 /// gROOT->Macro("../ConfigCalibTrain.C(119037)");
  57 /// //2
  58 /// InitTPCalign();
  59 /// MakeFits();      // this is logn proceure 30 minutes
  60 ///
  61 /// //UpdateOCDB(0,AliCDBRunRange::Infinity());
  62 /// //
  63 /// gROOT->Macro("~/NimStyle.C")
  64 /// ~~~
  65
  66 #if !defined(__CINT__) || defined(__MAKECINT__)
  67 #include "TH1D.h"
  68 #include "TH2F.h"
  69 #include "THnSparse.h"
  70 #include "TVectorD.h"
  71 #include "TTreeStream.h"
  72 #include "TFile.h"
  73 #include "AliTPCcalibAlign.h"
  74 #include "AliTPCROC.h"
  75 #include "AliTPCCalROC.h"
  76 #include "AliTPCCalPad.h"
  77 #include "TF1.h"
  78 #include "TH2.h"
  79 #include "TH3.h"
  80 #include "TROOT.h"
  81 #include "TProfile.h"
  82 #include "AliTPCPreprocessorOnline.h"
  83 #include "AliTPCcalibDB.h"
  84 #include "AliTPCkalmanAlign.h"
  85 #include "TPostScript.h"
  86 #include "TLegend.h"
  87 #include "TCut.h"
  88 #include "TCanvas.h"
  89 #include "TStyle.h"
  90 #include "AliTPCExBEffectiveSector.h"
  91 #include "AliTPCComposedCorrection.h"
  92 //
  93 #include "AliCDBMetaData.h"
  94 #include "AliCDBId.h"
  95 #include "AliCDBManager.h"
  96 #include "AliCDBStorage.h"
  97 #include "AliCDBEntry.h"
  98 #include "TStatToolkit.h"
  99 #endif
 100
 101
 102 TObject  *palign=0;
 103 TTree * treeP=0;
 104 TTree * treeM=0;
 105 TTree * tree0=0;
 106 TTree * treePZ=0;
 107 TTree * treeMZ=0;
 108 TTree * tree0Z=0;
 109 AliTPCROC *proc = AliTPCROC::Instance();
 110 AliTPCExBEffectiveSector * geffCorr=0;
 111
 112 Double_t GetCorr(Double_t sector, Double_t localX, Double_t kZ,Int_t type);
 113
 114 void MakeFits();
 115 void InitTPCalign();
 116
 117
 118 void MakeAlignCalPad(Int_t mode){
 119   /// Make AlignCalpad and make report
 120
 121   gSystem->Load("libANALYSIS");
 122   gSystem->Load("libSTAT");
 123   gSystem->Load("libTPCcalib");
 124   printf("Make report mode\t%d\n", mode);
 125   if (mode<1) {
 126     InitTPCalign();
 127     MakeFits();
 128   }
 129   printf("Make report\n");
 130 }
 131
 132
 133 void InitTPCalign(){
 134   /// read the TPC alignment
 135
 136   TFile fcalib("CalibObjects.root");
 137   TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
 138   if (array){
 139     palign = ( AliTPCcalibAlign *)array->FindObject("alignTPC");
 140   }else{
 141     palign = ( AliTPCcalibAlign *)fcalib.Get("alignTPC");
 142   }
 143   if (!palign){
 144      TFile fcalib2("TPCAlignObjects.root");
 145      palign = ( AliTPCcalibAlign *)fcalib2.Get("alignTPC");
 146   }
 147 }
 148
 149
 150 void MakeFits(){
 151   AliTPCcalibAlign * align =  (AliTPCcalibAlign *)palign;
 152   THnSparse * hdY = align->GetClusterDelta(0);
 153   THnSparse * hdZ = align->GetClusterDelta(1);
 154   AliTPCExBEffectiveSector::MakeResidualMap(hdY,"clusterDY.root");
 155   AliTPCExBEffectiveSector::MakeResidualMap(hdZ,"clusterDZ.root");
 156 }
 157
 158 void FitdY(TTree * tree){
 159   ///
 160
 161   AliTPCROC * roc = AliTPCROC::Instance();
 162   Double_t xquadrant = roc->GetPadRowRadii(36,53);
 163   Double_t xIO       = ( roc->GetPadRowRadii(0,roc->GetNRows(0)-1)+roc->GetPadRowRadii(36,0))*0.5;
 164
 165   //
 166   tree->SetAlias("diffS","(sector-int(sector)-0.5)");
 167   tree->SetAlias("diffQ",Form("(localX-%f)",xquadrant));
 168   tree->SetAlias("diffIO",Form("(localX-%f)",xIO));
 169   tree->SetAlias("iroc",Form("(localX<%f)",xIO));   // bool IROC
 170   tree->SetAlias("dqLR0",Form("(localX<%f)*(-1+(diffS<0)*2)",xIO));  //sign LeftRight
 171   tree->SetAlias("dqLR2",Form("(localX>%f)*(-1+(diffS<0)*2)",xquadrant));  //sign LeftRight up
 172   //
 173   tree->SetAlias("diffIFC","abs(localX-80)");
 174   tree->SetAlias("diffOFC","abs(250-localX)");
 175   tree->SetAlias("drift","(1-abs(kZ*localX)/250.)");
 176   //
 177   tree->SetAlias("dIFC2",Form("drift/(1+abs(diffIFC^2/(%f^2)))",10.));
 178   tree->SetAlias("dIFC",Form("drift/(1+abs(diffIFC/(%f)))",5.));
 179   tree->SetAlias("dOFC2",Form("drift/(1+abs(diffOFC^2/(%f^2)))",10.));
 180   tree->SetAlias("dOFC",Form("drift/(1+abs(diffOFC/(%f)))",5.));
 181
 182
 183
 184   TStatToolkit toolkit;
 185   Double_t chi2=0;
 186   Int_t    npoints=0;
 187   TVectorD param;
 188   TMatrixD covar;
 189   TString  fstringG="";
 190   fstringG+="(iroc)++";           //1 iroc shift
 191   fstringG+="(dqLR0)++";
 192   fstringG+="(dqLR2)++";
 193   fstringG+="(diffIO*iroc)++";     // iroc rotation
 194   fstringG+="(diffIO*dqLR0)++";
 195   fstringG+="(diffQ*dqLR2)++";
 196   //
 197   fstringG+="(dIFC+dIFC2)++";            //9
 198   fstringG+="(diffS*(dIFC+dIFC2))++";
 199   fstringG+="(diffS^2*(dIFC+dIFC2))++";
 200   fstringG+="(dIFC-dIFC2)++";            //9
 201   fstringG+="(diffS*(dIFC-dIFC2))++";
 202   fstringG+="(diffS^2*(dIFC-dIFC2))++";
 203   //
 204   //
 205   fstringG+="(dOFC+dOFC2)++";            //9
 206   fstringG+="(diffS*(dOFC+dOFC2))++";
 207   fstringG+="(diffS^2*(dOFC+dOFC2))++";
 208   fstringG+="(dOFC-dOFC2)++";            //9
 209   fstringG+="(diffS*(dOFC-dOFC2))++";
 210   fstringG+="(diffS^2*(dOFC-dOFC2))++";
 211   //
 212   //
 213   //
 214   TTreeSRedirector *pcstream=new TTreeSRedirector("dyFit.root");
 215   {for (Int_t isec=0; isec<18; isec++){
 216       {for (Int_t iside=0; iside<=1; iside++){
 217           TCut cutS=Form("abs(sector-%d.5)<0.5&&kZ*%d>0&&entries>2000",isec,(iside>0)?-1:1);
 218           TString *strFitG = TStatToolkit::FitPlane(tree,"mean", fstringG.Data(),cutS, chi2,npoints,param,covar,-1,0, 10000000, kTRUE);
 219           tree->SetAlias("fitG",strFitG->Data());
 220           tree->Draw("mean-fitG",cutS,"");
 221           //
 222           printf("%d\t%d\t%f\n",iside,isec,TMath::Sqrt(chi2/npoints));
 223           (*pcstream)<<"fitDy"<<
 224             "iside="<<iside<<
 225             "sec="<<isec<<
 226             "chi2="<<chi2<<
 227             "npoints="<<npoints<<
 228             "p.="<<&param<<
 229             "\n";
 230         }
 231       }
 232     }
 233   }
 234   delete pcstream;
 235 }
 236
 237
 238
 239 void DumpDerivative(TH3 * his){
 240   ///
 241
 242   Int_t nbins0=his->GetXaxis()->GetNbins();
 243   Int_t nbins1=his->GetYaxis()->GetNbins();
 244   Int_t nbins2=his->GetZaxis()->GetNbins();
 245   TTreeSRedirector *pcstream = new TTreeSRedirector("aaa.root");
 246   for (Int_t ibin0=1; ibin0<nbins0; ibin0++)
 247     for (Int_t ibin1=1; ibin1<nbins1; ibin1++)
 248       for (Int_t ibin2=1; ibin2<nbins2; ibin2++){
 249         Float_t x= his->GetXaxis()->GetBinCenter(ibin0);
 250         Float_t y= his->GetYaxis()->GetBinCenter(ibin1);
 251         Float_t z= his->GetZaxis()->GetBinCenter(ibin2);
 252         Float_t c= his->GetBinContent(ibin0,ibin1,ibin2);
 253         Float_t c0M=his->GetBinContent(ibin0-1,ibin1,ibin2);
 254         Float_t c0P=his->GetBinContent(ibin0+1,ibin1,ibin2);
 255         Float_t c1M=his->GetBinContent(ibin0,ibin1-1,ibin2);
 256         Float_t c1P=his->GetBinContent(ibin0,ibin1+1,ibin2);
 257         Float_t c2M=his->GetBinContent(ibin0,ibin1,ibin2-1);
 258         Float_t c2P=his->GetBinContent(ibin0,ibin1,ibin2+1);
 259         printf("%f\t%f\t%f\t%f\n",x,y,z,c);
 260         (*pcstream)<<"aaa"<<
 261           "x="<<x<<
 262           "y="<<y<<
 263           "z="<<z<<
 264           "c="<<c<<
 265           "c0M="<<c0M<<
 266           "c0P="<<c0P<<
 267           "c1M="<<c1M<<
 268           "c1P="<<c1P<<
 269           "c2M="<<c2M<<
 270           "c2P="<<c2P<<
 271           "\n";
 272       }
 273   delete pcstream;
 274 }
 275
 276 Double_t GetCorr(Double_t sector, Double_t localX, Double_t kZ, Int_t type){
 277   /// calculate the correction at given position - check the geffCorr
 278
 279   Double_t phi=sector*TMath::Pi()/9.;
 280   Double_t gx = localX*TMath::Cos(phi);
 281   Double_t gy = localX*TMath::Sin(phi);
 282   Double_t gz = localX*kZ;
 283   Int_t nsector=(gz>0) ? 0:18;
 284   //
 285   //
 286   //
 287   Float_t distPoint[3]={gx,gy,gz};
 288   geffCorr->DistortPoint(distPoint, nsector);
 289   Double_t r0=TMath::Sqrt(gx*gx+gy*gy);
 290   Double_t r1=TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
 291   Double_t phi0=TMath::ATan2(gy,gx);
 292   Double_t phi1=TMath::ATan2(distPoint[1],distPoint[0]);
 293   if (type==0) return r1-r0;
 294   if (type==1) return (phi1-phi0)*r0;
 295   return phi1-phi0;
 296 }
 297
 298
 299
 300 void LoadDistortionTrees(){
 301   /// Load distortion tree
 302
 303   TFile *fp = new TFile("clusterDYPlus.root");
 304   TFile *fm = new TFile("clusterDYMinus.root");
 305   TFile *f0 = new TFile("clusterDY0.root");
 306   TFile *fpz= new TFile("clusterDZPlus.root");
 307   TFile *fmz= new TFile("clusterDZMinus.root");
 308   TFile *f0z=new TFile("clusterDZ0.root");
 309   treeP=(TTree*)fp->Get("delta");
 310   treeM=(TTree*)fm->Get("delta");
 311   tree0=(TTree*)f0->Get("delta");
 312   treePZ=(TTree*)fpz->Get("delta");
 313   treeMZ=(TTree*)fmz->Get("delta");
 314   tree0Z=(TTree*)f0z->Get("delta");
 315   tree0->AddFriend(treeP,"P");
 316   tree0->AddFriend(treeM,"M");
 317   tree0->AddFriend(treePZ,"PZ");
 318   tree0->AddFriend(treeMZ,"MZ");
 319   tree0->AddFriend(tree0Z,"Z");
 320   tree0->SetMarkerStyle(25);
 321   tree0->SetMarkerSize(0.4);
 322 }
 323
 324 void UpdateEffSectorOCDB(){
 325   /// Incremeantal update ot the correction maps
 326   /// corrections on top of previous corrections
 327
 328   TFile fp("clusterDYPlus.root");
 329   TFile fm("clusterDYMinus.root");
 330   TFile f0("clusterDY0.root");
 331   TFile f0z("clusterDZ0.root");
 332
 333   TH3F *his3D0=(TH3F*)f0.Get("his3D");
 334   TH3F *his3DP=(TH3F*)fp.Get("his3D");
 335   TH3F *his3DM=(TH3F*)fm.Get("his3D");
 336   TH3F *his3DZ=(TH3F*)f0z.Get("his3D");
 337   TH3F *hisR=0;   //half of difference (Plus-Minus) scaled by c1
 338   hisR=(TH3F*)his3DP->Clone();
 339   hisR->Add(his3DM,-1);
 340   hisR->Scale(0.5);
 341   hisR->Scale(1/0.3); // c1 factor to be added there //is sign correct?
 342   //
 343   //
 344   AliTPCExBEffectiveSector *effSector = new  AliTPCExBEffectiveSector;
 345   effSector->SetName("EffSector");
 346   effSector->SetTitle("EffSector");
 347   effSector->fCorrectionRPhi=(TH3F*)his3D0->Clone();
 348   effSector->fCorrectionR=(TH3F*)hisR->Clone();
 349   effSector->fCorrectionZ=(TH3F*)his3DZ->Clone();
 350   geffCorr=effSector;
 351   AliTPCCorrection::AddVisualCorrection(geffCorr,0);
 352   //
 353   //
 354   //
 355   AliTPCcalibDB * calib = AliTPCcalibDB::Instance();
 356   TObjArray *arr = calib->GetTPCComposedCorrectionArray();
 357   {for (Int_t i=0; i<arr->GetEntries(); i++){
 358     AliTPCComposedCorrection * corr=(AliTPCComposedCorrection*)arr->At(i);
 359     if (!corr) continue;
 360     AliTPCExBEffectiveSector * corrSec=( AliTPCExBEffectiveSector *)corr->GetCorrections()->FindObject("EffSector");
 361     if (i==0) geffCorr=(AliTPCExBEffectiveSector *)corrSec;
 362     if (!corrSec) corr->GetCorrections()->Add(effSector->Clone());
 363     if (corrSec){
 364       if (corrSec->fCorrectionRPhi) corrSec->fCorrectionRPhi->Add(effSector->fCorrectionRPhi);
 365       else corrSec->fCorrectionRPhi=(TH3F*)effSector->fCorrectionRPhi->Clone();
 366       if (corrSec->fCorrectionR) corrSec->fCorrectionR->Add(effSector->fCorrectionR);
 367       else corrSec->fCorrectionR=(TH3F*)effSector->fCorrectionR->Clone();
 368       if (corrSec->fCorrectionZ) corrSec->fCorrectionR->Add(effSector->fCorrectionZ);
 369       else corrSec->fCorrectionZ=(TH3F*)effSector->fCorrectionZ->Clone();
 370     }
 371     if (i==0) AliTPCCorrection::AddVisualCorrection(corrSec,1);
 372     }}
 373 // make OCDB entries
 374   TString userName=gSystem->GetFromPipe("echo $USER");
 375   TString date=gSystem->GetFromPipe("date");
 376   TString ocdbStorage="local:////";
 377   ocdbStorage+=gSystem->GetFromPipe("pwd");
 378   ocdbStorage+="/OCDB";
 379   //
 380   AliCDBMetaData *metaData= new AliCDBMetaData();
 381   metaData->SetObjectClassName("TObjArray");
 382   metaData->SetResponsible("Marian Ivanov");
 383   metaData->SetBeamPeriod(1);
 384   metaData->SetAliRootVersion("05-27-04"); //root version
 385   metaData->SetComment(Form("Correction calibration. User: %s\n Data%s",userName.Data(),date.Data()));
 386   AliCDBId* id1=NULL;
 387   id1=new AliCDBId("TPC/Calib/Correction", 0, AliCDBRunRange::Infinity());
 388   AliCDBStorage* gStorage = 0;
 389
 390   gStorage=AliCDBManager::Instance()->GetStorage((ocdbStorage+"Update").Data());
 391   gStorage->Put(arr, (*id1), metaData);
 392
 393
 394 }
 395
 396
 397 void DrawDiff(){
 398   TFile f0("../mergeField0/mean.root");
 399   TFile fP("../mergePlus/mean.root");
 400   TFile fM("../mergeMinus/mean.root");
 401   //
 402   TTree *itsdy0=(TTree*)f0.Get("ITSdy");
 403   TTree *itsdyP=(TTree*)fP.Get("ITSdy");
 404   TTree *itsdyM=(TTree*)fM.Get("ITSdy");
 405   TTree *vdy0=(TTree*)f0.Get("Vertexdy");
 406   TTree *vdyP=(TTree*)fP.Get("Vertexdy");
 407   TTree *vdyM=(TTree*)fM.Get("Vertexdy");
 408   itsdy0->SetMarkerStyle(25);
 409   itsdy0->SetMarkerSize(0.3);
 410   itsdy0->AddFriend(itsdyP,"P");
 411   itsdy0->AddFriend(itsdyM,"M");
 412   itsdy0->AddFriend(vdy0,"V");
 413   itsdy0->AddFriend(vdyP,"VP");
 414   itsdy0->AddFriend(vdyM,"VM");
 415   itsdy0->SetMarkerStyle(25);
 416   itsdy0->SetMarkerSize(0.4);
 417
 418 }
 419
 420
 421
 422
 423
 424 void MakePlotDeltaZ(){
 425   ///
 426
 427   TCut cut="entries>500&&PZ.entries>500&&MZ.entries>500";
 428   TCanvas *cA = new TCanvas("deltaZA","deltaZA",900,700);
 429   TCanvas *cC = new TCanvas("deltaZC","deltaZC",900,700);
 430   TCanvas *cARef = new TCanvas("deltaZARef","deltaZARef",1000,800);
 431   TCanvas *cCRef = new TCanvas("deltaZCRef","deltaZCRef",1000,800);
 432   //
 433   TH1::AddDirectory(0);
 434   TH1 * his=0;
 435   cA->Divide(4,5);
 436   for (Int_t isec=0; isec<18; isec++){
 437     cA->cd(isec+1);
 438     TCut cutS=Form("abs(sector-0.5-%d)<0.1",isec);
 439     tree0->Draw("Z.mean*10.:localX:abs(kZ)",cutS+"entries>4000&&P.entries>200&&localX<250&&kZ>0&&abs(kZ)<1","colz");
 440     his=(TH1*)tree0->GetHistogram()->Clone();
 441     his->GetXaxis()->SetTitle("r (cm)");
 442     his->GetYaxis()->SetTitle("#Delta_{z} (mm)");
 443     his->GetZaxis()->SetTitle("tan($theta)");
 444     his->Draw("colz");
 445   }
 446   cA->cd(19);
 447
 448   cC->Divide(4,5);
 449   for (Int_t isec=0; isec<18; isec++){
 450     cC->cd(isec+1);
 451     TCut cutS=Form("abs(sector-0.5-%d)<0.1",isec);
 452     tree0->Draw("Z.mean*10.:localX:abs(kZ)",cutS+"entries>4000&&P.entries>200&&localX<250&&kZ<0&&abs(kZ)<1","colz");
 453     his=(TH1*)tree0->GetHistogram()->Clone();
 454     his->GetXaxis()->SetTitle("r (cm)");
 455     his->GetYaxis()->SetTitle("#Delta_{z} (mm)");
 456     his->GetZaxis()->SetTitle("tan($theta)");
 457     his->Draw("colz");
 458
 459   }
 460
 461   cARef->Divide(4,5);
 462   for (Int_t isec=0; isec<18; isec++){
 463     cARef->cd(isec+1);
 464     TCut cutS=Form("abs(sector-0.5-%d)<0.1",isec);
 465     tree0->Draw("(PZ.mean-MZ.mean)*10.:localX:abs(kZ)",cutS+cut+"kZ>0&&abs(kZ)<1","colz");
 466     his=(TH1*)tree0->GetHistogram()->Clone();
 467     his->GetXaxis()->SetTitle("r (cm)");
 468     his->GetYaxis()->SetTitle("#Delta_{z} (mm)");
 469     his->GetZaxis()->SetTitle("tan($theta)");
 470     his->Draw("colz");
 471
 472   }
 473
 474   cCRef->Divide(4,5);
 475   for (Int_t isec=0; isec<18; isec++){
 476     cCRef->cd(isec+1);
 477     TCut cutS=Form("abs(sector-0.5-%d)<0.1",isec);
 478     tree0->Draw("(PZ.mean-MZ.mean)*10.:localX:abs(kZ)",cutS+cut+"kZ<0&&abs(kZ)<1","colz");
 479     his=(TH1*)tree0->GetHistogram()->Clone();
 480     his->GetXaxis()->SetTitle("r (cm)");
 481     his->GetYaxis()->SetTitle("#Delta_{z} (mm)");
 482     his->GetZaxis()->SetTitle("tan($theta)");
 483     his->Draw("colz");
 484   }
 485
 486   TPostScript *ps = new TPostScript("distortionZ.ps", 112);
 487   ps->NewPage();
 488   cA->Draw();
 489   ps->NewPage();
 490   cA->Draw();
 491   ps->NewPage();
 492   cC->Draw();
 493   ps->NewPage();
 494   cARef->Draw();
 495   ps->NewPage();
 496   cCRef->Draw();
 497   ps->Close();
 498 }
 499
 500
 501
 502
 503 void MakeAlign(){
 504   /// make  sector alignment - using Kalman filter method -AliTPCkalmanAlign
 505   /// Combined information is used, mean residuals are minimized:
 506   ///
 507   /// 1. TPC-TPC sector alignment
 508   /// 2. TPC-ITS alignment
 509   /// 3. TPC vertex alignment
 510
 511   TFile fcalib("../mergeField0/TPCAlignObjects.root");
 512   AliTPCcalibAlign * align = ( AliTPCcalibAlign *)fcalib.Get("alignTPC");
 513   TFile f0("../mergeField0/mean.root");
 514
 515   //
 516   TTree *itsdy=(TTree*)f0.Get("ITSdy");
 517   TTree *itsdp=(TTree*)f0.Get("ITSdsnp");
 518   TTree *itsdz=(TTree*)f0.Get("ITSdz");
 519   TTree *itsdt=(TTree*)f0.Get("ITSdtheta");
 520   //
 521   TTree *vdy=(TTree*)f0.Get("Vertexdy");
 522   TTree *vds=(TTree*)f0.Get("Vertexdsnp");
 523   TTree *vdz=(TTree*)f0.Get("Vertexdz");
 524   TTree *vdt=(TTree*)f0.Get("Vertexdtheta");
 525
 526   itsdy->AddFriend(itsdp,"Snp");
 527   itsdy->AddFriend(itsdz,"Z");
 528   itsdy->AddFriend(itsdt,"T");
 529   //
 530   itsdy->AddFriend(vdy,"V");
 531   itsdy->AddFriend(vds,"VSnp");
 532   itsdy->AddFriend(vdz,"VZ");
 533   itsdy->AddFriend(vdt,"VT");
 534   itsdy->SetMarkerStyle(25);
 535   itsdy->SetMarkerSize(0.4);
 536
 537   TCut cutQ="entries>500&&abs(theta)<0.8&&abs(snp)<0.2";
 538   TH1F his1("hdeltaY1","hdeltaY1",100,-0.5,0.5);
 539   TMatrixD vecAlign(72,1);
 540   TMatrixD covAlign(72,72);
 541   AliTPCkalmanAlign::BookAlign1D(vecAlign,covAlign,0,0.0005);
 542   TVectorD vecITSY(72);
 543   TVectorD vecITSS(72);
 544   TVectorD vecVS(72);
 545   TVectorD vecITSTan(72);
 546   TVectorD vecVTan(72);
 547   {for (Int_t isec0=0; isec0<36; isec0++){
 548       Double_t phi0=(isec0%18+0.5)*TMath::Pi()/9.;
 549       if (phi0>TMath::Pi()) phi0-=TMath::TwoPi();
 550       Int_t iside0=(isec0%36<18)? 0:1;
 551       TCut cutSector=Form("abs(%f-phi)<0.14",phi0);
 552       TCut cutSide = (iside0==0)? "theta>0":"theta<0";
 553       itsdy->Draw("mean",cutQ+cutSector+cutSide);
 554       Double_t meanITSY=itsdy->GetHistogram()->GetMean()/83.6;
 555       vecITSY[isec0]=meanITSY;
 556       vecITSY[isec0+36]=meanITSY;
 557       itsdy->Draw("Snp.mean",cutQ+cutSector+cutSide);
 558       Double_t meanITSS=itsdy->GetHistogram()->GetMean();
 559       vecITSS[isec0]=meanITSS;
 560       vecITSS[isec0+36]=meanITSS;
 561       itsdy->Draw("VSnp.mean",cutQ+cutSector+cutSide);
 562       Double_t meanVS=itsdy->GetHistogram()->GetMean();
 563       vecVS[isec0]=meanVS;
 564       vecVS[isec0+36]=meanVS;
 565     }
 566   }
 567   AliTPCROC * roc = AliTPCROC::Instance();
 568   Double_t fXIROC = (roc->GetPadRowRadii(0,0)+roc->GetPadRowRadii(0,roc->GetNRows(0)-1))*0.5;
 569   Double_t fXOROC = (roc->GetPadRowRadii(36,0)+roc->GetPadRowRadii(36,roc->GetNRows(36)-1))*0.5;
 570   Double_t fXmiddle   = ( roc->GetPadRowRadii(0,0)+roc->GetPadRowRadii(36,roc->GetNRows(36)-1))*0.5;
 571   Double_t fXIO       = ( roc->GetPadRowRadii(0,roc->GetNRows(0)-1)+roc->GetPadRowRadii(36,0))*0.5;
 572
 573   TTreeSRedirector *pcstream=new TTreeSRedirector("combAlign.root");
 574
 575   {
 576     for (Int_t iter=0; iter<2; iter++){
 577     for (Int_t isec0=0; isec0<72; isec0++){
 578     for (Int_t isec1=0; isec1<72; isec1++){
 579       TH1 * his = align->GetHisto(AliTPCcalibAlign::kY,isec0,isec1);
 580       TH1 * hisPhi = align->GetHisto(AliTPCcalibAlign::kPhi,isec0,isec1);
 581       if (!his) continue;
 582       if (his->GetEntries()<100) continue;
 583       Double_t xref=fXIO;
 584       if (isec0<36&&isec1<36) xref=fXIROC;
 585       if (isec0>=36&&isec1>=36) xref=fXOROC;
 586       Double_t meanTPC=his->GetMean()/xref;
 587       Double_t meanTPCPhi=hisPhi->GetMean();
 588       Double_t meanITS0=vecITSY[isec0];
 589       Double_t meanITS1=vecITSY[isec1];
 590       Double_t meanITSS0=vecITSS[isec0];
 591       Double_t meanITSS1=vecITSS[isec1];
 592       Double_t meanVS0=vecVS[isec0];
 593       Double_t meanVS1=vecVS[isec1];
 594       AliTPCkalmanAlign::UpdateAlign1D(meanTPC, 0.001, isec0,isec1,  vecAlign,covAlign);
 595       AliTPCkalmanAlign::UpdateAlign1D(meanTPCPhi, 0.001, isec0,isec1,  vecAlign,covAlign);
 596       AliTPCkalmanAlign::UpdateAlign1D(meanITS1-meanITS0, 0.001, isec0,isec1,  vecAlign,covAlign);
 597       AliTPCkalmanAlign::UpdateAlign1D(meanITSS1-meanITSS0, 0.001, isec0,isec1,  vecAlign,covAlign);
 598       AliTPCkalmanAlign::UpdateAlign1D(meanVS1-meanVS0, 0.001, isec0,isec1,  vecAlign,covAlign);
 599       //printf("isec0\t%d\tisec1\t%d\t%f\t%f\t%f\n",isec0,isec1, meanTPC, meanITS0,meanITS1);
 600       Double_t kalman0= vecAlign(isec0,0);
 601       Double_t kalman1= vecAlign(isec1,0);
 602       if (iter>0) (*pcstream)<<"align"<<
 603         "iter="<<iter<<
 604         "xref="<<xref<<
 605         "isec0="<<isec0<<
 606         "isec1="<<isec1<<
 607         "mTPC="<<meanTPC<<
 608         "mTPCPhi="<<meanTPCPhi<<
 609         "mITS0="<<meanITS0<<
 610         "mITS1="<<meanITS1<<
 611         "mITSS0="<<meanITSS0<<
 612         "mITSS1="<<meanITSS1<<
 613         "mVS0="<<meanVS0<<
 614         "mVS1="<<meanVS1<<
 615         "k0="<<kalman0<<
 616         "k1="<<kalman1<<
 617         "\n";
 618     }
 619     }
 620     }
 621   }
 622   pcstream->GetFile()->cd();
 623   vecAlign.Write("alignPhiMean");
 624   covAlign.Write("alingPhiCovar");
 625   delete pcstream;
 626   TFile f("combAlign.root");
 627   TTree * treeA = (TTree*)f.Get("align");
 628   treeA->SetMarkerStyle(25);
 629   treeA->SetMarkerSize(0.5);
 630 }
 631
 632
 633