TPC/CalibMacros/MakeAlignCalPad.C

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