TPC/CalibMacros/FitAlignCombined.C

   1 /*
   2   marian.ivanov@cern.ch
   3   Macro to create  alignment/distortion maps
   4   As a input output of AliTPCcalibAlign and AliTPCcalibTime is used.
   5   distortion lookup tables are used.
   6
   7   Input file mean.root with distortion tree expected to be in directory:
   8   ../mergeField0/mean.root
   9
  10
  11   The ouput file fitAlignCombined.root contains:
  12   1. Resulting (residual) AliTPCCalibMisalignment
  13   2. QA fit plots
  14
  15   Functions documented inside:
  16
  17   RegisterAlignFunction();
  18   MakeAlignFunctionGlobal();
  19   MakeAlignFunctionSector();
  20
  21 */
  22
  23
  24 /*
  25   Example usage:
  26   //
  27   .x ~/NimStyle.C
  28   gROOT->Macro("~/rootlogon.C");
  29   gSystem->Load("libANALYSIS");
  30   gSystem->Load("libSTAT");
  31   gSystem->Load("libTPCcalib");
  32   gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros -I$ALICE_ROOT/TPC/TPC -I$ALICE_ROOT/STAT");
  33   .L $ALICE_ROOT/TPC/CalibMacros/FitAlignCombined.C+
  34   .x ConfigCalibTrain.C(119047)
  35   //
  36   //
  37   FitAlignCombinedCorr();
  38
  39
  40 */
  41
  42
  43
  44 #if !defined(__CINT__) || defined(__MAKECINT__)
  45 #include "TH1D.h"
  46 #include "TH2F.h"
  47 #include "THnSparse.h"
  48 #include "TVectorD.h"
  49 #include "TTreeStream.h"
  50 #include "TFile.h"
  51 #include "TChain.h"
  52 #include "AliTPCcalibAlign.h"
  53 #include "AliTPCROC.h"
  54 #include "AliTPCCalROC.h"
  55 #include "AliTPCCalPad.h"
  56 #include "TF1.h"
  57 #include "TH2.h"
  58 #include "TH3.h"
  59 #include "TROOT.h"
  60 #include "TProfile.h"
  61 #include "AliTPCPreprocessorOnline.h"
  62 #include "AliTPCcalibDB.h"
  63 #include "AliTPCkalmanAlign.h"
  64 #include "TPostScript.h"
  65 #include "TLegend.h"
  66 #include "TCut.h"
  67 #include "TCanvas.h"
  68 #include "TStyle.h"
  69 #include "AliLog.h"
  70 #include "AliTPCExBEffectiveSector.h"
  71 #include "AliTPCComposedCorrection.h"
  72 #include "AliTPCCalibGlobalMisalignment.h"
  73 //
  74 #include "AliCDBMetaData.h"
  75 #include "AliCDBId.h"
  76 #include "AliCDBManager.h"
  77 #include "AliCDBStorage.h"
  78 #include "AliCDBEntry.h"
  79 #include "TStatToolkit.h"
  80 #include "AliAlignObjParams.h"
  81 #include "AliTPCParam.h"
  82 #endif
  83
  84
  85 void DrawFitQA();
  86
  87 AliTPCROC *proc = AliTPCROC::Instance();
  88 AliTPCCalibGlobalMisalignment *combAlignOCDBOld=0;  // old misalignment from OCDB - used for reconstruction
  89 AliTPCCalibGlobalMisalignment *combAlignOCDBNew=0;  // new misalignment
  90 AliTPCCalibGlobalMisalignment *combAlignGlobal=0;   // delta misalignment - global part
  91 AliTPCCalibGlobalMisalignment *combAlignLocal=0;    // delta misalignment - sector/local part
  92 //
  93
  94 TTreeSRedirector *pcstream= 0;    // Workspace to store the fit parameters and QA histograms
  95 //
  96 TChain * chain=0;
  97 TChain * chainZ=0;                // trees with z measurement
  98 TTree * tree=0;
  99 TTree *itsdy=0;
 100 TTree *itsdyP=0;
 101 TTree *itsdyM=0;
 102 TTree *itsdp=0;
 103 TTree *itsdphiP=0;
 104 TTree *itsdphiM=0;
 105
 106 TTree *itsdz=0;
 107 TTree *itsdt=0;
 108 TTree *tofdy=0;
 109 TTree *trddy=0;
 110 //
 111 TTree *vdy=0;
 112 TTree *vdyP=0;
 113 TTree *vdyM=0;
 114 TTree *vds=0;
 115 TTree *vdz=0;
 116 TTree *vdt=0;
 117
 118 TCut cutS="entries>1000&&abs(snp)<0.2&&abs(theta)<1.";
 119
 120
 121 void RegisterAlignFunction(){
 122   //
 123   // Register primitive alignment components.
 124   // Linear conbination of primitev forulas used for fit
 125   // The nominal delta 1 mm in shift and 1 mrad in rotation
 126   // Primitive formulas registeren in AliTPCCoreection::AddvisualCorrection
 127   // 0 - deltaX
 128   // 1 - deltaY
 129   // 2 - deltaZ
 130   // 3 - rot0 (phi)
 131   // 4 - rot1 (theta)
 132   // 5 - rot2
 133   TGeoHMatrix matrixX;
 134   TGeoHMatrix matrixY;
 135   TGeoHMatrix matrixZ;
 136   TGeoRotation rot0;
 137   TGeoRotation rot1;
 138   TGeoRotation rot2;  //transformation matrices
 139   TGeoRotation rot90;  //transformation matrices
 140   matrixX.SetDx(0.1); matrixY.SetDy(0.1); matrixZ.SetDz(0.1); //1 mm translation
 141   rot0.SetAngles(0.001*TMath::RadToDeg(),0,0);
 142   rot1.SetAngles(0,0.001*TMath::RadToDeg(),0);
 143   rot2.SetAngles(0,0,0.001*TMath::RadToDeg());
 144   //how to get rot02 ?
 145   rot90.SetAngles(0,90,0);
 146   rot2.MultiplyBy(&rot90,kTRUE);
 147   rot90.SetAngles(0,-90,0);
 148   rot2.MultiplyBy(&rot90,kFALSE);
 149   AliTPCCalibGlobalMisalignment *alignRot0  =new  AliTPCCalibGlobalMisalignment;
 150   alignRot0->SetAlignGlobal(&rot0);
 151   AliTPCCalibGlobalMisalignment *alignRot1=new  AliTPCCalibGlobalMisalignment;
 152   alignRot1->SetAlignGlobal(&rot1);
 153   AliTPCCalibGlobalMisalignment *alignRot2=new  AliTPCCalibGlobalMisalignment;
 154   alignRot2->SetAlignGlobal(&rot2);
 155   //
 156   AliTPCCalibGlobalMisalignment *alignTrans0  =new  AliTPCCalibGlobalMisalignment;
 157   alignTrans0->SetAlignGlobal(&matrixX);
 158   AliTPCCalibGlobalMisalignment *alignTrans1=new  AliTPCCalibGlobalMisalignment;
 159   alignTrans1->SetAlignGlobal(&matrixY);
 160   AliTPCCalibGlobalMisalignment *alignTrans2=new  AliTPCCalibGlobalMisalignment;
 161   alignTrans2->SetAlignGlobal(&matrixZ);
 162   AliTPCCorrection::AddVisualCorrection(alignTrans0  ,0);
 163   AliTPCCorrection::AddVisualCorrection(alignTrans1  ,1);
 164   AliTPCCorrection::AddVisualCorrection(alignTrans2  ,2);
 165   AliTPCCorrection::AddVisualCorrection(alignRot0    ,3);
 166   AliTPCCorrection::AddVisualCorrection(alignRot1    ,4);
 167   AliTPCCorrection::AddVisualCorrection(alignRot2    ,5);
 168   TObjArray arrAlign(6);
 169   arrAlign.AddAt(alignTrans0->Clone(),0);
 170   arrAlign.AddAt(alignTrans1->Clone(),1);
 171   arrAlign.AddAt(alignTrans2->Clone(),2);
 172   arrAlign.AddAt(alignRot0->Clone(),3);
 173   arrAlign.AddAt(alignRot1->Clone(),4);
 174   arrAlign.AddAt(alignRot2->Clone(),5);
 175   //combAlign.SetCorrections((TObjArray*)arrAlign.Clone());
 176 }
 177
 178 AliTPCCalibGlobalMisalignment * MakeAlignFunctionGlobal(TVectorD paramYGlobal){
 179   //
 180   // Take a fit parameters and make a combined correction
 181   // 1. Take the common part
 182   // 3. Make combined AliTPCCalibGlobalMisalignment - register it
 183   // 4. Compare the aliases with fit values - IT is OK
 184   //
 185   AliTPCCalibGlobalMisalignment *alignGlobal  =new  AliTPCCalibGlobalMisalignment;
 186   TGeoHMatrix matGlobal; // global parameters
 187   TGeoHMatrix matDelta;  // delta A side - C side
 188   //
 189   TGeoHMatrix matrixX;
 190   TGeoHMatrix matrixY;
 191   TGeoHMatrix matrixZ;
 192   TGeoRotation rot0;
 193   TGeoRotation rot1;
 194   TGeoRotation rot2;  //transformation matrices
 195   TGeoRotation rot90;  //transformation matrices
 196   //
 197   //
 198   matrixX.SetDx(0.1*paramYGlobal[1]);
 199   matrixY.SetDy(0.1*paramYGlobal[2]);
 200   rot0.SetAngles(0.001*TMath::RadToDeg()*paramYGlobal[3],0,0);
 201   rot1.SetAngles(0,0.001*TMath::RadToDeg()*paramYGlobal[4],0);
 202   rot2.SetAngles(0,0,0.001*TMath::RadToDeg()*paramYGlobal[5]);
 203   rot90.SetAngles(0,90,0);
 204   rot2.MultiplyBy(&rot90,kTRUE);
 205   rot90.SetAngles(0,-90,0);
 206   rot2.MultiplyBy(&rot90,kFALSE);
 207   matGlobal.Multiply(&matrixX);
 208   matGlobal.Multiply(&matrixY);
 209   matGlobal.Multiply(&rot0);
 210   matGlobal.Multiply(&rot1);
 211   matGlobal.Multiply(&rot2);
 212   alignGlobal->SetAlignGlobal((TGeoMatrix*)matGlobal.Clone());
 213   AliTPCCorrection::AddVisualCorrection(alignGlobal    ,100);
 214   //
 215   /*
 216     chain->Draw("mean-deltaG:int(sector)",cutS+"type==0&&refX==0&&theta>0","profsame");
 217     chain->Draw("mean-deltaG:int(sector+18)",cutS+"type==0&&refX==0&&theta<0","profsame");
 218
 219     chain->Draw("deltaG:fitYGlobal",cutS+"type==0");
 220     chain->Draw("deltaG:fitYGlobal",cutS+"type==2");
 221
 222   */
 223
 224   return alignGlobal;
 225 }
 226
 227
 228 AliTPCCalibGlobalMisalignment * MakeAlignFunctionSector(TVectorD paramYLocal){
 229   //
 230   // Take a fit parameters and make a combined correction:
 231   // Only delta local Y and delta phi are fitted - not sensityvity for other parameters
 232   // Algorithm:
 233   //   1. Loop over sectors
 234   //   2. Make combined AliTPCCalibGlobalMisalignment - register it
 235   //   3. Compare the aliases with fit values - IT is OK
 236   //
 237   AliTPCCalibGlobalMisalignment *alignLocal  =new  AliTPCCalibGlobalMisalignment;
 238   TGeoHMatrix matrixX;
 239   TGeoHMatrix matrixY;
 240   TGeoHMatrix matrixZ;
 241   TGeoRotation rot0;
 242   TGeoRotation rot1;
 243   TGeoRotation rot2;  //transformation matrices
 244   TGeoRotation rot90;  //transformation matrices
 245   TObjArray * array = new TObjArray(72);
 246   TVectorD vecLY(72);
 247   TVectorD vecGY(72);
 248   TVectorD vecGX(72);
 249   TVectorD vecPhi(72);
 250   TVectorD vecSec(72);
 251   //
 252   //
 253   {for (Int_t isec=0; isec<72; isec++){
 254       Int_t sector=isec%18;
 255       Int_t offset=sector*4+1;
 256       if ((isec%36)>=18) offset+=2;
 257       Double_t phi= (Double_t(sector)+0.5)*TMath::Pi()/9.;
 258       //
 259       Double_t dly = paramYLocal[offset];
 260       Double_t dphi= paramYLocal[offset+1];
 261       Double_t dgx =   TMath::Cos(phi)*0+TMath::Sin(phi)*dly;
 262       Double_t dgy =   TMath::Sin(phi)*0-TMath::Cos(phi)*dly;
 263       vecSec[isec]=isec;
 264       vecLY[isec]=dly;
 265       vecGX[isec]=dgx;
 266       vecGY[isec]=dgy;
 267       vecPhi[isec]=dphi;
 268       //
 269       matrixX.SetDx(dgx);
 270       matrixY.SetDy(dgy);
 271       rot0.SetAngles(0.001*TMath::RadToDeg()*dphi,0,0);
 272       TGeoHMatrix matrixSector; // global parameters
 273       matrixSector.Multiply(&matrixX);
 274       matrixSector.Multiply(&matrixY);
 275       matrixSector.Multiply(&rot0);
 276       array->AddAt(matrixSector.Clone(),isec);
 277     }}
 278   alignLocal->SetAlignSectors(array);
 279   AliTPCCorrection::AddVisualCorrection(alignLocal,101);
 280   //
 281   //
 282   TGraph * graphLY = new TGraph(72,vecSec.GetMatrixArray(), vecLY.GetMatrixArray());
 283   TGraph * graphGX = new TGraph(72,vecSec.GetMatrixArray(), vecGX.GetMatrixArray());
 284   TGraph * graphGY = new TGraph(72,vecSec.GetMatrixArray(), vecGY.GetMatrixArray());
 285   TGraph * graphPhi = new TGraph(72,vecSec.GetMatrixArray(), vecPhi.GetMatrixArray());
 286   graphLY->SetMarkerStyle(25); graphGX->SetMarkerStyle(25);  graphGY->SetMarkerStyle(25);   graphPhi->SetMarkerStyle(25);
 287   graphLY->SetName("DeltaLY"); graphLY->SetTitle("#Delta_{ly} (mm)");
 288   graphGX->SetName("DeltaGX"); graphGX->SetTitle("#Delta_{gx} (mm)");
 289   graphGY->SetName("DeltaGY"); graphGY->SetTitle("#Delta_{gy} (mm)");
 290   graphPhi->SetName("DeltaPhi"); graphPhi->SetTitle("#Delta_{phi} (mrad)");
 291   //
 292   graphLY->Write("grDeltaLY");
 293   graphGX->Write("grDeltaGX");
 294   graphGY->Write("grDeltaGY");
 295   graphPhi->Write("grDeltaPhi");
 296   // Check:
 297   /*
 298     graphLY.Draw("alp")
 299     chain->Draw("mean-deltaG:int(sector)",cutS+"type==0&&refX==0&&theta>0","profsame");
 300     chain->Draw("mean-deltaG:int(sector+18)",cutS+"type==0&&refX==0&&theta<0","profsame");
 301
 302     graphPhi.Draw("alp")
 303     chain->Draw("1000*(mean-deltaG):int(sector)",cutS+"type==2&&refX==0&&theta>0","profsame");
 304     chain->Draw("1000*(mean-deltaG):int(sector+18)",cutS+"type==2&&refX==0&&theta<0","profsame");
 305
 306     //
 307     chain->Draw("deltaL:fitYLocal",cutS+"type==2&&refX==0","");  // phi shift - OK
 308     chain->Draw("deltaL:fitYLocal",cutS+"type==0&&refX==0","");  // OK
 309     chain->Draw("deltaL:fitYLocal",cutS+"type==0",""); // OK
 310
 311   */
 312
 313   return alignLocal;
 314 }
 315
 316
 317
 318
 319
 320 void LoadTrees(){
 321   //
 322   // make  sector alignment - using Kalman filter method -AliTPCkalmanAlign
 323   // Combined information is used, mean residuals are minimized:
 324   //
 325   // 1. TPC-ITS alignment
 326   // 2. TPC vertex alignment
 327   //
 328   TFile *f0 = new TFile("../mergeField0/mean.root");
 329   TFile *fP= new TFile("../mergePlus/mean.root");
 330   TFile *fM= new TFile("../mergeMinus/mean.root");
 331   //
 332   itsdy=(TTree*)f0->Get("ITSdy");
 333   itsdyP=(TTree*)fP->Get("ITSdy");
 334   itsdyM=(TTree*)fM->Get("ITSdy");
 335   itsdp=(TTree*)f0->Get("ITSdsnp");
 336   itsdphiP=(TTree*)fP->Get("ITSdsnp");
 337   itsdphiM=(TTree*)fM->Get("ITSdsnp");
 338
 339   itsdz=(TTree*)f0->Get("ITSdz");
 340   itsdt=(TTree*)f0->Get("ITSdtheta");
 341   tofdy=(TTree*)f0->Get("TOFdy");
 342   trddy=(TTree*)f0->Get("TRDdy");
 343   //
 344   vdy=(TTree*)f0->Get("Vertexdy");
 345   vdyP=(TTree*)fP->Get("Vertexdy");
 346   vdyM=(TTree*)fM->Get("Vertexdy");
 347   vds=(TTree*)f0->Get("Vertexdsnp");
 348   vdz=(TTree*)f0->Get("Vertexdz");
 349   vdt=(TTree*)f0->Get("Vertexdtheta");
 350   chain    = new TChain("mean","mean");
 351   chainZ    = new TChain("mean","mean");
 352   chain->AddFile("../mergeField0/mean.root",10000000,"ITSdy");
 353   chain->AddFile("../mergeField0/mean.root",10000000,"ITSdsnp");
 354   chain->AddFile("../mergeField0/mean.root",10000000,"Vertexdy");
 355   chain->AddFile("../mergeField0/mean.root",10000000,"Vertexdsnp");
 356   chain->AddFile("../mergeField0/mean.root",10000000,"TOFdy");
 357   chain->AddFile("../mergeField0/mean.root",10000000,"TRDdy");
 358   //
 359   chainZ->AddFile("../mergeField0/mean.root",10000000,"Vertexdz");
 360   chainZ->AddFile("../mergeField0/mean.root",10000000,"Vertexdtheta");
 361   chainZ->AddFile("../mergeField0/mean.root",10000000,"ITSdz");
 362   chainZ->AddFile("../mergeField0/mean.root",10000000,"ITSdtheta");
 363
 364   //
 365   itsdy->AddFriend(itsdp,"Phi");
 366   itsdy->AddFriend(itsdphiP,"PhiP");
 367   itsdy->AddFriend(itsdphiM,"PhiM");
 368   itsdy->AddFriend(itsdz,"Z");
 369   itsdy->AddFriend(itsdt,"T");
 370   itsdy->AddFriend(itsdyP,"YP");
 371   itsdy->AddFriend(itsdyM,"YM");
 372   //
 373   itsdy->AddFriend(vdy,"V");
 374   itsdy->AddFriend(vdyP,"VP");
 375   itsdy->AddFriend(vdyP,"VM");
 376   itsdy->AddFriend(tofdy,"TOF");
 377   itsdy->AddFriend(trddy,"TRD");
 378   itsdy->AddFriend(vds,"VPhi");
 379   itsdy->AddFriend(vdz,"VZ");
 380   itsdy->AddFriend(vdt,"VT");
 381   itsdy->AddFriend(tofdy,"TOF.");
 382   itsdy->SetMarkerStyle(25);
 383   itsdy->SetMarkerSize(0.4);
 384   tree=itsdy;
 385   tree->SetAlias("side","(-1+2*(theta>0))");
 386   chain->SetAlias("side","(-1+2*(theta>0))");
 387   chain->SetMarkerStyle(25);
 388   chain->SetMarkerSize(0.5);
 389 }
 390
 391
 392 void FitAlignCombinedCorr(){
 393   //
 394   // Fit Global X and globalY shift at vertex and at ITS
 395   //
 396   RegisterAlignFunction();
 397   LoadTrees();
 398   combAlignOCDBOld = AliTPCCalibGlobalMisalignment::CreateOCDBAlign();
 399   //
 400
 401   pcstream= new TTreeSRedirector("fitAlignCombined.root");
 402
 403   TStatToolkit toolkit;
 404   Double_t chi2=0;
 405   Int_t    npoints=0;
 406   // Fit vectors
 407   // global fits
 408   TVectorD vecSec(37);
 409   TVectorD paramYGlobal;
 410   TVectorD paramYLocal;
 411   TMatrixD covar;
 412   // make Aliases
 413   {for (Int_t isec=0; isec<18; isec++){ // sectors
 414       tree->SetAlias(Form("sec%d",isec), Form("(abs(sector-%f)<0.5)",isec+0.5));
 415       chain->SetAlias(Form("sec%d",isec), Form("(abs(sector-%f)<0.5)",isec+0.5));
 416     }}
 417   for (Int_t isec=-1;isec<36; isec++){
 418     vecSec[isec+1]=isec;
 419   }
 420   chain->SetAlias("err","((type==0)*0.1+(type==2)*0.001)");
 421   // delta phi
 422   chain->SetAlias("dpgx","(0+0)");
 423   chain->SetAlias("dpgy","(0+0)");
 424   chain->SetAlias("dpgr0","((AliTPCCorrection::GetCorrSector(sector,245,theta,1,3)-AliTPCCorrection::GetCorrSector(sector,85,theta,1,3))/160)");
 425   chain->SetAlias("dpgr1","((AliTPCCorrection::GetCorrSector(sector,245,theta,1,4)-AliTPCCorrection::GetCorrSector(sector,85,theta,1,4))/160)");
 426   chain->SetAlias("dpgr2","((AliTPCCorrection::GetCorrSector(sector,245,theta,1,5)-AliTPCCorrection::GetCorrSector(sector,85,theta,1,5))/160)");
 427   chain->SetAlias("deltaPG","((AliTPCCorrection::GetCorrSector(sector,245,theta,1,100)-AliTPCCorrection::GetCorrSector(sector,85,theta,1,100))/160)");
 428   chain->SetAlias("deltaPL","((AliTPCCorrection::GetCorrSector(sector,245,theta,1,101)-AliTPCCorrection::GetCorrSector(sector,85,theta,1,101))/160)");
 429   // delta y at 85 cm
 430   chain->SetAlias("dygxT","(AliTPCCorrection::GetCorrSector(sector,85,theta,1,0)+0)");
 431   chain->SetAlias("dygyT","(AliTPCCorrection::GetCorrSector(sector,85,theta,1,1)+0)");
 432   chain->SetAlias("dyr0T","(AliTPCCorrection::GetCorrSector(sector,85,theta,1,3)+0)");
 433   chain->SetAlias("dyr1T","(AliTPCCorrection::GetCorrSector(sector,85,theta,1,4)+0)");
 434   chain->SetAlias("dyr2T","(AliTPCCorrection::GetCorrSector(sector,85,theta,1,5)+0)");
 435   chain->SetAlias("deltaYGT","(AliTPCCorrection::GetCorrSector(sector,85,theta,1,100)+0)");
 436   chain->SetAlias("deltaYLT","(AliTPCCorrection::GetCorrSector(sector,85,theta,1,101)+0)");
 437   //
 438   // delta y at reference X
 439   chain->SetAlias("dygx","(dygxT+0)");   // due global x shift
 440   chain->SetAlias("dygy","(dygyT+0)");   // due global y shift
 441   chain->SetAlias("dyr0","(dyr0T+dpgr0*(refX-85.))");  // due rotation 0
 442   chain->SetAlias("dyr1","(dyr1T+dpgr1*(refX-85.))");  // due rotation 1
 443   chain->SetAlias("dyr2","(dyr2T+dpgr2*(refX-85.))");  // due rotation 2
 444   chain->SetAlias("deltaYG","(deltaYGT+deltaPG*(refX-85.))");  // due global distortion
 445   chain->SetAlias("deltaYL","(deltaYLT+deltaPL*(refX-85.))");  // due local distortion
 446
 447   chain->SetAlias("deltaG","(type==0)*deltaYG+(type==2)*deltaPG");  //alias to global function
 448   chain->SetAlias("deltaL","(type==0)*deltaYL+(type==2)*deltaPL");  //alias to local function
 449   //
 450   TString  fstringGlobal="";
 451   // global part
 452   fstringGlobal+="((type==0)*dygx+0)++";
 453   fstringGlobal+="((type==0)*dygy+0)++";
 454   fstringGlobal+="((type==0)*dyr0+(type==2)*dpgr0)++";
 455   fstringGlobal+="((type==0)*dyr1+(type==2)*dpgr1)++";
 456   fstringGlobal+="((type==0)*dyr2+(type==2)*dpgr2)++";
 457   //
 458   // Make global fits
 459   //
 460   TString *strFitYGlobal = TStatToolkit::FitPlane(chain,"mean:err", fstringGlobal.Data(),cutS+"", chi2,npoints,paramYGlobal,covar,-1,0, 10000000, kTRUE);
 461   combAlignGlobal =MakeAlignFunctionGlobal(paramYGlobal);
 462   printf("chi2=%f\n",TMath::Sqrt(chi2/npoints));
 463   chain->SetAlias("fitYGlobal",strFitYGlobal->Data());
 464   paramYGlobal.Print();
 465   paramYGlobal.Write("paramYGlobal");
 466   //
 467   //
 468   //
 469   //
 470   {for (Int_t isec=0; isec<18; isec++){
 471       //A side
 472       chain->SetAlias(Form("glyASec%d",isec),Form("((type==0)*sec%d*(theta>0))",isec)); // ly shift
 473       chain->SetAlias(Form("gxASec%d",isec),Form("((type==0)*dygx*sec%d*(theta>0))",isec)); // gx shift
 474       chain->SetAlias(Form("gyASec%d",isec),Form("((type==0)*dygy*sec%d*(theta>0))",isec)); // gy shift
 475       chain->SetAlias(Form("gpASec%d",isec),Form("(((type==0)*dyr0+(type==2)*dpgr0)*sec%d)*(theta>0)",isec)); // phi rotation
 476       //C side
 477       chain->SetAlias(Form("glyCSec%d",isec),Form("((type==0)*sec%d*(theta<0))",isec)); // ly shift
 478       chain->SetAlias(Form("gxCSec%d",isec),Form("((type==0)*dygx*sec%d*(theta<0))",isec));
 479       chain->SetAlias(Form("gyCSec%d",isec),Form("((type==0)*dygy*sec%d*(theta<0))",isec));
 480       chain->SetAlias(Form("gpCSec%d",isec),Form("(((type==0)*dyr0+(type==2)*dpgr0)*sec%d)*(theta<0)",isec)); // phi rotation
 481     }}
 482
 483   TString  fstringLocal="";
 484   {for (Int_t isec=0; isec<18; isec++){
 485       //fstringLocal+=Form("((type==0)*dygx*sec%d*(theta>0))++",isec);
 486       //      fstringLocal+=Form("(gxASec%d)++",isec);
 487       // fstringLocal+=Form("(gyASec%d)++",isec);
 488       fstringLocal+=Form("(glyASec%d)++",isec);
 489       fstringLocal+=Form("(gpASec%d)++",isec);
 490       fstringLocal+=Form("(glyCSec%d)++",isec);
 491       //fstringLocal+=Form("(gxCSec%d)++",isec);
 492       //fstringLocal+=Form("(gyCSec%d)++",isec);
 493       fstringLocal+=Form("(gpCSec%d)++",isec);
 494     }}
 495
 496   TString *strFitYLocal = TStatToolkit::FitPlane(chain,"(mean-deltaG):err", fstringLocal.Data(),cutS+"", chi2,npoints,paramYLocal,covar,-1,0, 10000000, kTRUE);
 497   printf("chi2=%f\n",TMath::Sqrt(chi2/npoints));
 498   chain->SetAlias("fitYLocal",strFitYLocal->Data());
 499   combAlignLocal =MakeAlignFunctionSector(paramYLocal);
 500   //paramYLocal.Print();
 501   paramYLocal.Write("paramYLocal");
 502   //
 503   // scaling  - why do we need it?
 504   TString  fstringScale="";
 505   fstringScale+="((type==1)*refX+(type==2))*dsec++";
 506   fstringScale+="((type==1)*refX+(type==2))*dsec*cos(phi)++";
 507   fstringScale+="((type==1)*refX+(type==2))*dsec*sin(phi)++";
 508   fstringScale+="((type==1)*refX+(type==2))*dsec*side++";
 509   fstringScale+="((type==1)*refX+(type==2))*dsec*side*cos(phi)++";
 510   fstringScale+="((type==1)*refX+(type==2))*dsec*side*sin(phi)++";
 511   //
 512   //
 513   //
 514   pcstream->GetFile()->cd();
 515   DrawFitQA();
 516   //
 517   AliTPCCalibGlobalMisalignment * combAlignOCDBNew0 =( AliTPCCalibGlobalMisalignment *)combAlignOCDBOld->Clone();
 518   combAlignOCDBNew0->AddAlign(*combAlignGlobal);
 519   combAlignOCDBNew0->AddAlign(*combAlignLocal);
 520   combAlignOCDBNew= AliTPCCalibGlobalMisalignment::CreateMeanAlign(combAlignOCDBNew0);
 521   //
 522   combAlignGlobal->SetName("fcombAlignGlobal");
 523   combAlignLocal->SetName("fcombAlignLocal");
 524   combAlignOCDBOld->SetName("fcombAlignOCDBOld");
 525   combAlignOCDBNew->SetName("fcombAlignOCDBNew");
 526   combAlignOCDBNew0->SetName("fcombAlignOCDBNew0");
 527   //
 528   combAlignGlobal->Write("fcombAlignGlobal");
 529   combAlignLocal->Write("fcombAlignLocal");
 530   combAlignOCDBOld->Write("fcombAlignOCDBOld");
 531   combAlignOCDBNew->Write("fcombAlignOCDBNew");
 532   combAlignOCDBNew0->Write("fcombAlignOCDBNew0");
 533   combAlignOCDBNew->Print();
 534   delete pcstream;
 535 }
 536
 537 void DrawFitQA(){
 538   //
 539   //
 540   //
 541  //
 542   // MakeQA plot 1D
 543   //
 544   TCanvas c;
 545   c.SetLeftMargin(0.15);
 546   chain->Draw("1000*(mean-deltaG)>>his(100,-1.5,1.5)",cutS+"type==2&&refX==0","");
 547   chain->GetHistogram()->SetName("TPC_VertexDphiGlobal1D");
 548   chain->GetHistogram()->SetTitle("TPC-Vertex #Delta_{#phi} - (Global Fit)");
 549   chain->GetHistogram()->GetXaxis()->SetTitle("#Delta_{#phi} (mrad)");
 550   chain->GetHistogram()->Fit("gaus","qnr");
 551   chain->GetHistogram()->Write();
 552   //
 553   chain->Draw("1000*(mean-deltaG)>>his(100,-1.5,1.5)",cutS+"type==2&&abs(refX-85)<2","");
 554   chain->GetHistogram()->SetName("TPC_ITSDphiGlobal1D");
 555   chain->GetHistogram()->SetTitle("TPC-ITS #Delta_{#phi} - (Global Fit)");
 556   chain->GetHistogram()->GetXaxis()->SetTitle("#Delta_{#phi} (mrad)");
 557   chain->GetHistogram()->Fit("gaus","qnr");
 558   chain->GetHistogram()->Write();
 559   //
 560   chain->Draw("10*(mean-deltaG)>>his(100,-1,1)",cutS+"type==0&&abs(refX-85)<2","");
 561   chain->GetHistogram()->SetName("TPC_ITSDRPhiGlobal1D");
 562   chain->GetHistogram()->SetTitle("TPC-ITS #Delta_{r#phi} - (Global Fit)");
 563   chain->GetHistogram()->GetXaxis()->SetTitle("#Delta_{r#phi} (mm)");
 564   chain->GetHistogram()->Fit("gaus","qnr");
 565   chain->GetHistogram()->Write();
 566   //
 567   chain->Draw("10*(mean-deltaG)>>his(100,-1,1)",cutS+"type==0&&abs(refX-0)<2","");
 568   chain->GetHistogram()->SetName("TPC-VertexDRPhiGlobal1D");
 569   chain->GetHistogram()->SetTitle("TPC-Vertex #Delta_{r#phi} - (Global Fit)");
 570   chain->GetHistogram()->GetXaxis()->SetTitle("#Delta_{r#phi} (mm)");
 571   chain->GetHistogram()->Fit("gaus","qnr");
 572   chain->GetHistogram()->Write();
 573   //
 574   // Make QA plot 3D
 575   //
 576   chain->Draw("1000*(mean-deltaG):sector:abs(theta)",cutS+"theta>0&&type==2&&refX==0","colz");
 577   chain->GetHistogram()->SetName("TPC_VertexDphi3DGlobal3D");
 578   chain->GetHistogram()->SetTitle("TPC-Vertex #Delta_{#phi} - (Global Fit)");
 579   chain->GetHistogram()->GetYaxis()->SetTitle("#Delta_{#phi} (mrad)");
 580   chain->GetHistogram()->GetXaxis()->SetTitle("sector");
 581   chain->GetHistogram()->GetZaxis()->SetTitle("tan(#Theta)");
 582   chain->GetHistogram()->Draw("colz");
 583   chain->GetHistogram()->Write();
 584   //
 585   chain->Draw("1000*(mean-deltaG):sector:abs(theta)",cutS+"theta>0&&type==2&&abs(refX-85)<2","colz");
 586   chain->GetHistogram()->SetName("TPC_ITSDphi3DGlobal3D");
 587   chain->GetHistogram()->SetTitle("TPC-ITS #Delta_{#phi} - (Global Fit)");
 588   chain->GetHistogram()->GetYaxis()->SetTitle("#Delta_{#phi} (mrad)");
 589   chain->GetHistogram()->GetXaxis()->SetTitle("sector");
 590   chain->GetHistogram()->GetZaxis()->SetTitle("tan(#Theta)");
 591   chain->GetHistogram()->Draw("colz");
 592   chain->GetHistogram()->Write();
 593   //
 594   chain->Draw("10*(mean-deltaG):sector:abs(theta)",cutS+"theta>0&&type==0&&abs(refX-85)<2","colz");
 595   chain->GetHistogram()->SetName("TPC_ITSDRphi3DGlobal3D");
 596   chain->GetHistogram()->SetTitle("TPC-ITS #Delta_{r#phi} - (Global Fit)");
 597   chain->GetHistogram()->GetYaxis()->SetTitle("#Delta_{r#phi} (mm)");
 598   chain->GetHistogram()->GetXaxis()->SetTitle("sector");
 599   chain->GetHistogram()->GetZaxis()->SetTitle("tan(#Theta)");
 600   chain->GetHistogram()->Draw("colz");
 601   chain->GetHistogram()->Write();
 602   //
 603   chain->Draw("10*(mean-deltaG):sector:abs(theta)",cutS+"theta>0&&type==0&&abs(refX-0)<2","colz");
 604   chain->GetHistogram()->SetName("TPC_VertexDRphi3DGlobal3D");
 605   chain->GetHistogram()->SetTitle("TPC-Vertex #Delta_{r#phi} - (Global Fit)");
 606   chain->GetHistogram()->GetYaxis()->SetTitle("#Delta_{r#phi} (mm)");
 607   chain->GetHistogram()->GetXaxis()->SetTitle("sector");
 608   chain->GetHistogram()->GetZaxis()->SetTitle("tan(#Theta)");
 609   chain->GetHistogram()->Draw("colz");
 610   chain->GetHistogram()->Write();
 611   //
 612   //
 613   //
 614   chain->Draw("1000*(mean-deltaG-deltaL):sector:abs(theta)",cutS+"theta>0&&type==2&&refX==0","colz");
 615   chain->GetHistogram()->SetName("TPC_VertexDphi3DLocal3D");
 616   chain->GetHistogram()->SetTitle("TPC-Vertex #Delta_{#phi} - (Local Fit)");
 617   chain->GetHistogram()->GetYaxis()->SetTitle("#Delta_{#phi} (mrad)");
 618   chain->GetHistogram()->GetXaxis()->SetTitle("sector");
 619   chain->GetHistogram()->GetZaxis()->SetTitle("tan(#Theta)");
 620   chain->GetHistogram()->Draw("colz");
 621   chain->GetHistogram()->Write();
 622   //
 623   chain->Draw("1000*(mean-deltaG-deltaL):sector:abs(theta)",cutS+"theta>0&&type==2&&abs(refX-85)<2","colz");
 624   chain->GetHistogram()->SetName("TPC_ITSDphi3DLocal3D");
 625   chain->GetHistogram()->SetTitle("TPC-ITS #Delta_{#phi} - (Local Fit)");
 626   chain->GetHistogram()->GetYaxis()->SetTitle("#Delta_{#phi} (mrad)");
 627   chain->GetHistogram()->GetXaxis()->SetTitle("sector");
 628   chain->GetHistogram()->GetZaxis()->SetTitle("tan(#Theta)");
 629   chain->GetHistogram()->Draw("colz");
 630   chain->GetHistogram()->Write();
 631   //
 632   chain->Draw("10*(mean-deltaG-deltaL):sector:abs(theta)",cutS+"theta>0&&type==0&&abs(refX-85)<2","colz");
 633   chain->GetHistogram()->SetName("TPC_ITSDRphi3DLocal3D");
 634   chain->GetHistogram()->SetTitle("TPC-ITS #Delta_{r#phi} - (Local Fit)");
 635   chain->GetHistogram()->GetYaxis()->SetTitle("#Delta_{r#phi} (mm)");
 636   chain->GetHistogram()->GetXaxis()->SetTitle("sector");
 637   chain->GetHistogram()->GetZaxis()->SetTitle("tan(#Theta)");
 638   chain->GetHistogram()->Draw("colz");
 639   chain->GetHistogram()->Write();
 640   //
 641   chain->Draw("10*(mean-deltaG-deltaL):sector:abs(theta)",cutS+"theta>0&&type==0&&abs(refX-0)<2","colz");
 642   chain->GetHistogram()->SetName("TPC_VertexDRphiLocal3D");
 643   chain->GetHistogram()->SetTitle("TPC-Vertex #Delta_{r#phi} - (Local Fit)");
 644   chain->GetHistogram()->GetYaxis()->SetTitle("#Delta_{r#phi} (mm)");
 645   chain->GetHistogram()->GetXaxis()->SetTitle("sector");
 646   chain->GetHistogram()->GetZaxis()->SetTitle("tan(#Theta)");
 647   chain->GetHistogram()->Draw("colz");
 648   chain->GetHistogram()->Write();
 649 }
 650
 651
 652
 653
 654 void FitAlignCombined0(){
 655   //
 656   // Fit Global X and globalY shift at vertex and at ITS
 657   //
 658
 659   TTreeSRedirector *pcstream= new TTreeSRedirector("fitAlignCombined.root");
 660
 661   TStatToolkit toolkit;
 662   Double_t chi2=0;
 663   Int_t    npoints=0;
 664   // Fit vectors
 665   // global fits
 666   TVectorD vecSec(37);
 667   TVectorD paramYVGlobal;
 668   TVectorD paramYITSGlobal;
 669   TVectorD paramPhiVGlobal;
 670   TVectorD paramPhiITSGlobal;
 671   // local fits
 672   TVectorD paramYVLocal;
 673   TVectorD paramPhiVLocal;
 674   TVectorD paramYITSLocal;
 675   TVectorD paramPhiITSLocal;
 676   TMatrixD covar;
 677   // make Aliases
 678   {for (Int_t isec=0; isec<18; isec++){ // sectors
 679       tree->SetAlias(Form("sec%d",isec), Form("(abs(sector-%f)<0.5)",isec+0.5));
 680     }}
 681   for (Int_t isec=-1;isec<36; isec++){
 682     vecSec[isec+1]=isec;
 683   }
 684   //
 685   TString  fstringGlobal="";
 686   fstringGlobal+="side++";
 687   fstringGlobal+="theta++";
 688   fstringGlobal+="cos(phi)*(theta>0)++";  // GX - A side
 689   fstringGlobal+="sin(phi)*(theta>0)++";  // GY - A side
 690   fstringGlobal+="cos(phi)*(theta<0)++";  // GX - C side
 691   fstringGlobal+="sin(phi)*(theta<0)++";  // GY - C side
 692   fstringGlobal+="cos(phi)*(theta)++";    // theta - get rid of rotation
 693   fstringGlobal+="sin(phi)*(theta)++";    // theta - get rid of rotation
 694   //
 695   // Local Fit
 696   //
 697   TString  fstringLocal="";
 698   {for (Int_t sector=0; sector<18; sector++){
 699       fstringLocal+=Form("((sec%d)*theta>0)++",sector);
 700       fstringLocal+=Form("((sec%d)*theta<0)++",sector);
 701     }}
 702   //
 703   // Make global fits
 704   //
 705   TString *strFitYVGlobal = TStatToolkit::FitPlane(tree,"V.mean", fstringGlobal.Data(),cutS+"", chi2,npoints,paramYVGlobal,covar,-1,0, 10000000, kFALSE);
 706   printf("chi2=%f\n",TMath::Sqrt(chi2/npoints));
 707   tree->SetAlias("fitYVGlobal",strFitYVGlobal->Data());
 708   tree->Draw("V.mean:fitYVGlobal",cutS);
 709   //
 710   TString *strFitPhiVGlobal = TStatToolkit::FitPlane(tree,"VPhi.mean", fstringGlobal.Data(),cutS+"", chi2,npoints,paramPhiVGlobal,covar,-1,0, 10000000, kFALSE);
 711   printf("chi2=%f\n",TMath::Sqrt(chi2/npoints));
 712   tree->SetAlias("fitPhiVGlobal",strFitPhiVGlobal->Data());
 713   tree->Draw("VPhi.mean:fitPhiVGlobal",cutS);
 714   //
 715   TString *strFitYITSGlobal = TStatToolkit::FitPlane(tree,"mean", fstringGlobal.Data(),cutS+"", chi2,npoints,paramYITSGlobal,covar,-1,0, 10000000, kFALSE);
 716   printf("chi2=%f\n",TMath::Sqrt(chi2/npoints));
 717   tree->SetAlias("fitYITSGlobal",strFitYITSGlobal->Data());
 718   tree->Draw("mean:fitYITSGlobal",cutS);
 719   TString *strFitPhiITSGlobal = TStatToolkit::FitPlane(tree,"Phi.mean", fstringGlobal.Data(),cutS+"", chi2,npoints,paramPhiITSGlobal,covar,-1,0, 10000000, kFALSE);
 720   printf("chi2=%f\n",TMath::Sqrt(chi2/npoints));
 721   tree->SetAlias("fitPhiITSGlobal",strFitPhiITSGlobal->Data());
 722   tree->Draw("Phi.mean:fitPhiITSGlobal",cutS);
 723   //
 724   // Residuals to the global fit
 725   //
 726   tree->SetAlias("dyVG", "V.mean-fitYVGlobal");
 727   tree->SetAlias("dphiVG","(VPhi.mean-fitPhiVGlobal)");
 728   tree->SetAlias("dyITSG","mean-fitYITSGlobal");
 729   tree->SetAlias("dphiITSG","(Phi.mean-fitPhiITSGlobal)");
 730
 731   TString *strFitYVLocal = TStatToolkit::FitPlane(tree,"dyVG", fstringLocal.Data(),cutS+"", chi2,npoints,paramYVLocal,covar,-1,0, 10000000, kFALSE);
 732   printf("chi2=%f\n",TMath::Sqrt(chi2/npoints));
 733   tree->SetAlias("fitYVLocal",strFitYVLocal->Data());
 734   tree->Draw("dyVG-fitYVLocal:sector:abs(theta)",cutS+"theta>0","colz");
 735   //
 736   TString *strFitPhiVLocal = TStatToolkit::FitPlane(tree,"dphiVG", fstringLocal.Data(),cutS+"", chi2,npoints,paramPhiVLocal,covar,-1,0, 10000000, kFALSE);
 737   printf("chi2=%f\n",TMath::Sqrt(chi2/npoints));
 738   tree->SetAlias("fitPhiVLocal",strFitPhiVLocal->Data());
 739   tree->Draw("dphiVG-fitPhiVLocal:sector:abs(theta)",cutS+"theta>0","colz");
 740   //
 741   //
 742   //
 743   TString *strFitYITSLocal = TStatToolkit::FitPlane(tree,"dyITSG", fstringLocal.Data(),cutS+"", chi2,npoints,paramYITSLocal,covar,-1,0, 10000000, kFALSE);
 744   printf("chi2=%f\n",TMath::Sqrt(chi2/npoints));
 745   tree->SetAlias("fitYITSLocal",strFitYITSLocal->Data());
 746   tree->Draw("dyITSG-fitYITSLocal:sector:abs(theta)",cutS+"theta>0","colz");
 747   //
 748   TString *strFitPhiITSLocal = TStatToolkit::FitPlane(tree,"dphiITSG", fstringLocal.Data(),cutS+"", chi2,npoints,paramPhiITSLocal,covar,-1,0, 10000000, kFALSE);
 749   printf("chi2=%f\n",TMath::Sqrt(chi2/npoints));
 750   tree->SetAlias("fitPhiITSLocal",strFitPhiITSLocal->Data());
 751   tree->Draw("dphiITSG-fitPhiITSLocal:sector:abs(theta)",cutS+"theta>0","colz");
 752
 753   //
 754   TH1D *hisA = 0;
 755   TH1D *hisC = 0;
 756   TVectorD vSec(36);
 757   TVectorD vecDPhi(36);
 758   TVectorD vecDLY(36);
 759   TVectorD vecDGX(36);
 760   TVectorD vecDGY(36);
 761
 762   //
 763   tree->SetAlias("phiMean","(fitPhiVLocal+fitPhiITSLocal+fitPhiVGlobal+fitPhiITSGlobal)*0.5");
 764   tree->SetAlias("yMeanV","((fitYITSLocal+fitYITSGlobal-85*phiMean)+(fitYVLocal+fitYVGlobal))*0.5");
 765   tree->SetAlias("yMeanITS","((fitYITSLocal+fitYITSGlobal)+(fitYVLocal+fitYVGlobal+85*phiMean))*0.5");
 766
 767
 768   tree->Draw("phiMean:int(sector)>>hhhA(18,0.18)",cutS+"abs(theta)<0.15&&theta>0","prof");
 769   hisA = (TH1D*) tree->GetHistogram()->Clone();
 770   tree->Draw("phiMean:int(sector)>>hhhC(18,0.18)",cutS+"abs(theta)<0.15&&theta<0","prof");
 771   hisC = (TH1D*) tree->GetHistogram()->Clone();
 772
 773   for (Int_t isec=0; isec<36; isec++){
 774     vSec[isec]=isec;
 775     if (isec<18)  vecDPhi[isec]=hisA->GetBinContent(isec%18+1);
 776     if (isec>=18) vecDPhi[isec]=hisC->GetBinContent(isec%18+1);
 777   }
 778   tree->Draw("yMeanV:int(sector)>>hhhA(18,0.18)",cutS+"abs(theta)<0.15&&theta>0","prof");
 779   hisA = (TH1D*) tree->GetHistogram()->Clone();
 780   tree->Draw("yMeanV:int(sector)>>hhhC(18,0.18)",cutS+"abs(theta)<0.15&&theta<0","prof");
 781   hisC = (TH1D*) tree->GetHistogram()->Clone();
 782   //
 783   for (Int_t isec=0; isec<36; isec++){
 784     if (isec<18)  vecDLY[isec]=hisA->GetBinContent(isec%18+1);
 785     if (isec>=18) vecDLY[isec]=hisC->GetBinContent(isec%18+1);
 786     vecDGX[isec]=-vecDLY[isec]*TMath::Sin(TMath::Pi()*(isec+0.5)/9.);
 787     vecDGY[isec]= vecDLY[isec]*TMath::Cos(TMath::Pi()*(isec+0.5)/9.);
 788   }
 789
 790
 791
 792
 793   //
 794   // Store results
 795   //
 796   {(*pcstream)<<"fitLocal"<<
 797       //global fits
 798       "sec.="<<&vSec<<
 799       "pYVGlobal.="<<&paramYVGlobal<<
 800       "pYITSGlobal.="<<&paramYITSGlobal<<
 801       "pPhiVGlobal.="<<&paramPhiVGlobal<<
 802       "pPhiITSGlobal.="<<&paramPhiITSGlobal<<
 803       // local fits
 804       "pYVLocal.="<<&paramYVLocal<<
 805       "pPhiVLocal.="<<&paramPhiVLocal<<
 806       "pYITSLocal.="<<&paramYITSLocal<<
 807       "pPhiITSLocal.="<<&paramPhiITSLocal<<
 808       //
 809       // combined
 810       "vDPhi.="<<&vecDPhi<<
 811       "vDLY.="<<&vecDLY<<
 812       "vDGX.="<<&vecDGX<<
 813       "vDGY.="<<&vecDGY<<
 814       "\n";
 815   }
 816   paramYVGlobal.Write("paramYVGlobal");
 817   paramYITSGlobal.Write("paramYITSGlobal");
 818   paramPhiVGlobal.Write("paramPhiVGlobal");
 819   paramPhiITSGlobal.Write("paramPhiITSGlobal");
 820   // local fits
 821   paramYVLocal.Write("paramYVLocal");
 822   paramPhiVLocal.Write("paramPhiVLocal");
 823   paramYITSLocal.Write("paramYITSLocal");
 824   paramPhiITSLocal.Write("paramPhiITSLocal");
 825   vecDPhi.Write("vecDPhi");
 826   vecDLY.Write("vecDLY");
 827   vecDGX.Write("vecDGX");
 828   vecDGY.Write("vecDGY");
 829
 830
 831
 832   tree->Draw("dyVG:sector:abs(theta)",cutS+"theta>0","colz");
 833   tree->GetHistogram()->Write("fitYVGlobal");
 834   tree->Draw("dphiVG:sector:abs(theta)",cutS+"theta>0","colz");
 835   tree->GetHistogram()->Write("fitPhiVGlobal");
 836   tree->Draw("dyITSG:sector:abs(theta)",cutS+"theta>0","colz");
 837   tree->GetHistogram()->Write("fitYITSGlobal");
 838   tree->Draw("dphiITSG:sector:abs(theta)",cutS+"theta>0","colz");
 839   tree->GetHistogram()->Write("fitPhiITSGlobal");
 840   //
 841   tree->Draw("dyVG-fitYVLocal:sector:abs(theta)",cutS+"theta>0","colz");
 842   tree->GetHistogram()->Write("fitYVLocal");
 843   tree->Draw("dphiVG-fitPhiVLocal:sector:abs(theta)",cutS+"theta>0","colz");
 844   tree->GetHistogram()->Write("fitPhiVLocal");
 845   tree->Draw("dyITSG-fitYITSLocal:sector:abs(theta)",cutS+"theta>0","colz");
 846   tree->GetHistogram()->Write("fitYITSLocal");
 847   tree->Draw("dphiITSG-fitPhiITSLocal:sector:abs(theta)",cutS+"theta>0","colz");
 848   tree->GetHistogram()->Write("fitPhiITSLocal");
 849   delete pcstream;
 850 }
 851
 852 void FitAlignCombined(){
 853   //
 854   //
 855   // make  sector alignment - using Kalman filter method -AliTPCkalmanAlign
 856   // Combined information is used, mean residuals are minimized:
 857   //
 858   // 1. TPC-TPC sector alignment
 859   // 2. TPC-ITS alignment
 860   // 3. TPC vertex alignment
 861   TFile fcalib("../mergeField0/TPCAlignObjects.root");
 862   AliTPCcalibAlign * align = ( AliTPCcalibAlign *)fcalib.Get("alignTPC");
 863
 864   TCut cutQ="entries>1000&&abs(theta)<0.5&&abs(snp)<0.2&&YP.entries>50&&YM.entries>50";
 865   TH1F his1("hdeltaY1","hdeltaY1",100,-0.5,0.5);
 866   TMatrixD vecAlign(72,1);
 867   TMatrixD covAlign(72,72);
 868   TMatrixD vecAlignY(72,1);
 869   TMatrixD covAlignY(72,72);
 870   TMatrixD vecAlignTheta(72,1);
 871   TMatrixD covAlignTheta(72,72);
 872   TMatrixD vecAlignZ(72,1);
 873   TMatrixD covAlignZ(72,72);
 874   AliTPCkalmanAlign::BookAlign1D(vecAlign,covAlign,0,0.002);
 875   AliTPCkalmanAlign::BookAlign1D(vecAlignY,covAlignY,0,0.1);
 876   AliTPCkalmanAlign::BookAlign1D(vecAlignTheta,covAlignTheta,0,0.002);
 877   AliTPCkalmanAlign::BookAlign1D(vecAlignZ,covAlignZ,0,0.1);
 878   TVectorD vecITSY(72);
 879   TVectorD vecITSYPM(72);
 880   TVectorD vecITSPhi(72);
 881   TVectorD vecITSPhiPM(72);
 882   TVectorD vecVY(72);
 883   TVectorD vecVS(72);
 884   TVectorD vecITSTheta(72);
 885   TVectorD vecVTheta(72);
 886   {for (Int_t isec0=0; isec0<36; isec0++){
 887       Double_t phi0=(isec0%18+0.5)*TMath::Pi()/9.;
 888       if (phi0>TMath::Pi()) phi0-=TMath::TwoPi();
 889       Int_t iside0=(isec0%36<18)? 0:1;
 890       TCut cutSector=Form("abs(%f-phi)<0.14",phi0);
 891       TCut cutSide = (iside0==0)? "theta>0":"theta<0";
 892       //
 893       itsdy->Draw("mean",cutQ+cutSector+cutSide);
 894       Double_t meanITSY=itsdy->GetHistogram()->GetMean();
 895       vecITSY[isec0]=meanITSY;
 896       vecITSY[isec0+36]=meanITSY;
 897       //
 898       itsdy->Draw("(YP.mean+YM.mean)*0.5",cutQ+cutSector+cutSide);
 899       Double_t meanITSYPM=itsdy->GetHistogram()->GetMean();
 900       vecITSYPM[isec0]=meanITSYPM;
 901       vecITSYPM[isec0+36]=meanITSYPM;
 902       //
 903       itsdy->Draw("Phi.mean",cutQ+cutSector+cutSide);
 904       Double_t meanITSPhi=itsdy->GetHistogram()->GetMean();
 905       vecITSPhi[isec0]=meanITSPhi;
 906       vecITSPhi[isec0+36]=meanITSPhi;
 907       //
 908       itsdy->Draw("(PhiP.mean+PhiM.mean)*0.5",cutQ+cutSector+cutSide);
 909       Double_t meanITSPhiPM=itsdy->GetHistogram()->GetMean();
 910       vecITSPhiPM[isec0]=meanITSPhiPM;
 911       vecITSPhiPM[isec0+36]=meanITSPhiPM;
 912       //
 913       //
 914       itsdy->Draw("VPhi.mean",cutQ+cutSector+cutSide);
 915       Double_t meanVS=itsdy->GetHistogram()->GetMean();
 916       vecVS[isec0]=meanVS;
 917       vecVS[isec0+36]=meanVS;
 918       //
 919       itsdy->Draw("V.mean",cutQ+cutSector+cutSide);
 920       Double_t meanVY=itsdy->GetHistogram()->GetMean();
 921       vecVY[isec0]=meanVY;
 922       vecVY[isec0+36]=meanVY;
 923       //
 924       itsdy->Draw("T.mean",cutQ+cutSector+cutSide);
 925       Double_t meanITST=itsdy->GetHistogram()->GetMean();
 926       vecITSTheta[isec0]=meanITST;
 927       vecITSTheta[isec0+36]=meanITST;
 928       //
 929       itsdy->Draw("VT.mean",cutQ+cutSector+cutSide);
 930       Double_t meanVT=itsdy->GetHistogram()->GetMean();
 931       vecVTheta[isec0]=meanVT;
 932       vecVTheta[isec0+36]=meanVT;
 933     }
 934   }
 935   AliTPCROC * roc = AliTPCROC::Instance();
 936   Double_t fXIROC = (roc->GetPadRowRadii(0,0)+roc->GetPadRowRadii(0,roc->GetNRows(0)-1))*0.5;
 937   Double_t fXOROC = (roc->GetPadRowRadii(36,0)+roc->GetPadRowRadii(36,roc->GetNRows(36)-1))*0.5;
 938   Double_t fXIO       = ( roc->GetPadRowRadii(0,roc->GetNRows(0)-1)+roc->GetPadRowRadii(36,0))*0.5;
 939
 940   TTreeSRedirector *pcstream=new TTreeSRedirector("combAlign.root");
 941
 942   {
 943     for (Int_t iter=0; iter<5; iter++){
 944     for (Int_t isec0=0; isec0<72; isec0++){
 945     for (Int_t isec1=0; isec1<72; isec1++){
 946       //if (isec0%36!=isec0%36) continue;
 947       if (iter==0 && isec0%36!=isec1%36) continue;
 948       if (iter==1 && isec0%18!=isec1%18) continue;
 949       TH1 * his = align->GetHisto(AliTPCcalibAlign::kY,isec0,isec1);
 950       TH1 * hisPhi = align->GetHisto(AliTPCcalibAlign::kPhi,isec0,isec1);
 951       TH1 * hisTheta = align->GetHisto(AliTPCcalibAlign::kTheta,isec0,isec1);
 952       TH1 * hisZ = align->GetHisto(AliTPCcalibAlign::kZ,isec0,isec1);
 953       Int_t side0=(isec0%36)<18 ? 1:-1;
 954       Int_t side1=(isec1%36)<18 ? 1:-1;
 955       Double_t weightTPC= 0.002;
 956       if (isec0%18==isec1%18) weightTPC=0.0005;
 957       if (isec0%36==isec1%36) weightTPC=0.00005;
 958       if (side0!=side1) weightTPC*=2.;
 959       Double_t weightTPCT= 0.001;
 960       if (isec0%18==isec1%18) weightTPC=0.0005;
 961       if (isec0%36==isec1%36) weightTPC=0.00005;
 962       if (TMath::Abs(isec0%18-isec1%18)==1) weightTPC = 0.0005;
 963       if (TMath::Abs(isec0%36-isec1%36)==1) weightTPC = 0.0001;
 964
 965       //
 966       Double_t meanITS0=vecITSY[isec0];
 967       Double_t meanITSPM0=vecITSYPM[isec0];
 968       Double_t meanITSPhi0=vecITSPhi[isec0];
 969       Double_t meanITSPhiPM0=vecITSPhiPM[isec0];
 970       Double_t meanVY0=vecVY[isec0];
 971       Double_t meanVPhi0=vecVS[isec0];
 972       Double_t meanITSTheta0=vecITSTheta[isec0];
 973       Double_t meanVTheta0=vecVTheta[isec0];
 974       //
 975       Double_t meanITS1=vecITSY[isec1];
 976       Double_t meanITSPM1=vecITSYPM[isec1];
 977       Double_t meanITSPhi1=vecITSPhi[isec1];
 978       Double_t meanITSPhiPM1=vecITSPhiPM[isec1];
 979       Double_t meanVY1=vecVY[isec1];
 980       Double_t meanVPhi1=vecVS[isec1];
 981       Double_t meanITSTheta1=vecITSTheta[isec1];
 982       Double_t meanVTheta1=vecVTheta[isec1];
 983       //
 984       Double_t meanPhi0 = (meanITSPhi0+meanVPhi0+2*meanITS0/83.)/4.;
 985       Double_t meanPhi1 = (meanITSPhi1+meanVPhi1+2*meanITS1/83.)/4.;
 986       //
 987       //
 988       if (iter>2 &&isec0==isec1){
 989         AliTPCkalmanAlign::UpdateAlign1D(-meanITS0/83, 0.001, isec0,  vecAlign,covAlign);
 990         AliTPCkalmanAlign::UpdateAlign1D(-meanITSPhi0, 0.001, isec0,  vecAlign,covAlign);
 991         AliTPCkalmanAlign::UpdateAlign1D(-meanVPhi0, 0.001, isec0,  vecAlign,covAlign);
 992         //
 993         AliTPCkalmanAlign::UpdateAlign1D(-meanPhi0, 0.001, isec0,  vecAlign,covAlign);
 994         AliTPCkalmanAlign::UpdateAlign1D(-meanVTheta0, 0.001, isec0,  vecAlignTheta,covAlignTheta);
 995         AliTPCkalmanAlign::UpdateAlign1D(-meanITSTheta0, 0.001, isec0,  vecAlignTheta,covAlignTheta);
 996       }
 997       if (iter>2){
 998         AliTPCkalmanAlign::UpdateAlign1D(meanPhi1-meanPhi0, 0.001, isec0,isec1,  vecAlign,covAlign);
 999         AliTPCkalmanAlign::UpdateAlign1D(meanITSPhiPM1-meanITSPhiPM0, 0.001, isec0,isec1,  vecAlign,covAlign);
1000         AliTPCkalmanAlign::UpdateAlign1D((meanITSPM1-meanITSPM0)/83., 0.001, isec0,isec1,  vecAlign,covAlign);
1001         AliTPCkalmanAlign::UpdateAlign1D(meanVTheta1-meanVTheta0, 0.001, isec0,isec1,  vecAlignTheta,covAlignTheta);
1002         AliTPCkalmanAlign::UpdateAlign1D(meanITSTheta1-meanITSTheta0, 0.001, isec0,isec1,  vecAlignTheta,covAlignTheta);
1003       }
1004       //
1005       if (!his) continue;
1006       if (!hisPhi) continue;
1007       if (!hisTheta) continue;
1008       if (his->GetEntries()<100) continue;
1009       Double_t xref=fXIO;
1010       if (isec0<36&&isec1<36) xref=fXIROC;
1011       if (isec0>=36&&isec1>=36) xref=fXOROC;
1012       Double_t meanTPC=his->GetMean();
1013       Double_t meanTPCPhi=hisPhi->GetMean();
1014       Double_t meanTPCTheta=hisTheta->GetMean();
1015       Double_t meanTPCZ=hisZ->GetMean();
1016       //
1017       //
1018       Double_t kalman0= vecAlign(isec0,0);
1019       Double_t kalman1= vecAlign(isec1,0);
1020       Double_t kalmanY0= vecAlignY(isec0,0);
1021       Double_t kalmanY1= vecAlignY(isec1,0);
1022       Double_t kalmanTheta0= vecAlignTheta(isec0,0);
1023       Double_t kalmanTheta1= vecAlignTheta(isec1,0);
1024       Double_t kalmanZ0= vecAlignZ(isec0,0);
1025       Double_t kalmanZ1= vecAlignZ(isec1,0);
1026       //
1027       //
1028       AliTPCkalmanAlign::UpdateAlign1D((meanTPC)/xref, weightTPC, isec0,isec1,  vecAlign,covAlign);
1029       AliTPCkalmanAlign::UpdateAlign1D(meanTPCPhi, weightTPC*10,isec0,isec1,  vecAlign,covAlign);
1030       //
1031       AliTPCkalmanAlign::UpdateAlign1D(meanTPCTheta, weightTPCT, isec0,isec1,  vecAlignTheta,covAlignTheta);
1032       if (side0==side1) AliTPCkalmanAlign::UpdateAlign1D(meanTPCZ, weightTPCT*100., isec0,isec1,  vecAlignZ,covAlignZ);
1033       //printf("isec0\t%d\tisec1\t%d\t%f\t%f\t%f\n",isec0,isec1, meanTPC, meanITS0,meanITS1);
1034
1035       if (iter>=0) (*pcstream)<<"align"<<
1036         "iter="<<iter<<
1037         "xref="<<xref<<                 // reference X
1038         "isec0="<<isec0<<               // sector number
1039         "isec1="<<isec1<<
1040         "side0="<<side0<<
1041         "side1="<<side1<<
1042         //TPC
1043         "mTPC="<<meanTPC<<              // delta Y  / xref
1044         "mTPCPhi="<<meanTPCPhi<<        // delta Phi
1045         "mTPCZ="<<meanTPCZ<<        // delta Z
1046         "mTPCTheta="<<meanTPCTheta<<    // delta Theta
1047         //ITS
1048         "mITS0="<<meanITS0<<
1049         "mITS1="<<meanITS1<<
1050         "mITSPhi0="<<meanITSPhi0<<
1051         "mITSPhi1="<<meanITSPhi1<<
1052         //
1053         "mITSPM0="<<meanITSPM0<<
1054         "mITSPM1="<<meanITSPM1<<
1055         "mITSPhiPM0="<<meanITSPhiPM0<<
1056         "mITSPhiPM1="<<meanITSPhiPM1<<
1057         //
1058         "mITSTheta0="<<meanITSTheta0<<
1059         "mITSTheta1="<<meanITSTheta1<<
1060         //Vertex
1061         "mVY0="<<meanVY0<<
1062         "mVY1="<<meanVY1<<
1063         "mVPhi0="<<meanVPhi0<<
1064         "mVPhi1="<<meanVPhi1<<
1065         "mVTheta0="<<meanVTheta0<<
1066         "mVTheta1="<<meanVTheta1<<
1067         // Vertex+ITS mean
1068         "mPhi0="<<meanPhi0<<
1069         "mPhi1="<<meanPhi1<<
1070         //Kalman
1071         "kY0="<<kalmanY0<<
1072         "kY1="<<kalmanY1<<
1073         "kPhi0="<<kalman0<<
1074         "kPhi1="<<kalman1<<
1075         "kTheta0="<<kalmanTheta0<<
1076         "kTheta1="<<kalmanTheta1<<
1077         "kZ0="<<kalmanZ0<<
1078         "kZ1="<<kalmanZ1<<
1079         "\n";
1080     }
1081     }
1082     }
1083   }
1084   pcstream->GetFile()->cd();
1085   vecAlign.Write("alignPhiMean");
1086   covAlign.Write("alingPhiCovar");
1087   vecAlignTheta.Write("alignThetaMean");
1088   covAlignTheta.Write("alingThetaCovar");
1089   vecAlignZ.Write("alignZMean");
1090   covAlignZ.Write("alingZCovar");
1091   delete pcstream;
1092   TFile f("combAlign.root");
1093   TTree * treeA = (TTree*)f.Get("align");
1094   treeA->SetMarkerStyle(25);
1095   treeA->SetMarkerSize(0.5);
1096 }
1097
1098
1099
1100
1101
1102 void UpdateOCDBAlign(){
1103   //
1104   // Store resulting OCDB entry
1105   // 0. Setup OCDB to get necccessary old entries - not done here
1106   // 1. Get old OCDB entry
1107   // 2. Get delta alignment
1108   // 3. Add delta alignment
1109   // 4. Store new alignment in
1110
1111   AliCDBEntry * entry = AliCDBManager::Instance()->Get("TPC/Align/Data");
1112   TClonesArray * array = (TClonesArray*)entry->GetObject();
1113   Int_t entries = array->GetEntries();
1114   TClonesArray * arrayNew=(TClonesArray*)array->Clone();
1115   TFile f("combAlign.root");
1116   TMatrixD *matPhi=(TMatrixD*)f.Get("alignPhiMean");
1117   //
1118   //
1119   { for (Int_t i=0;i<entries; i++){
1120       //
1121       //
1122       AliAlignObjParams *alignP = (AliAlignObjParams*)array->UncheckedAt(i);
1123       AliAlignObjParams *alignPNew = (AliAlignObjParams*)arrayNew->UncheckedAt(i);
1124       Int_t imod;
1125       AliGeomManager::ELayerID ilayer;
1126       alignP->GetVolUID(ilayer, imod);
1127       if (ilayer==AliGeomManager::kInvalidLayer) continue;
1128       Int_t sector=imod;
1129       if (ilayer==AliGeomManager::kTPC2) sector+=36;
1130       Double_t transOld[3], rotOld[3];
1131       alignP->GetTranslation(transOld);  // in cm
1132       alignP->GetAngles(rotOld);         // in degrees
1133       printf("%d\t%d\t%d\t",ilayer, imod,sector);
1134       printf("%f mrad \t%f mrad\n",1000*rotOld[2]*TMath::DegToRad(), 1000*((*matPhi)(sector,0)));
1135       alignPNew->SetRotation(rotOld[0],rotOld[1], rotOld[2]+(*matPhi)(sector,0)*TMath::RadToDeg());
1136       alignPNew->Print("kokot");
1137       alignP->Print("kokot");
1138     }
1139   }
1140
1141
1142
1143   TString ocdbStorage="local://"+gSystem->GetFromPipe("pwd")+"/OCDBUpdate";
1144   TString userName=gSystem->GetFromPipe("echo $USER");
1145   TString date=gSystem->GetFromPipe("date");
1146
1147   AliCDBMetaData *metaData= new AliCDBMetaData();
1148   metaData->SetObjectClassName("TObjArray");
1149   metaData->SetResponsible("Marian Ivanov");
1150   metaData->SetBeamPeriod(1);
1151   metaData->SetAliRootVersion("05-26-02"); //root version
1152   metaData->SetComment(Form("Correction calibration. User: %s\n Data%s",userName.Data(),date.Data()));
1153
1154   AliCDBId* id1=NULL;
1155   id1=new AliCDBId("TPC/Align/Data", 0, AliCDBRunRange::Infinity());
1156   AliCDBStorage* gStorage = AliCDBManager::Instance()->GetStorage(ocdbStorage);
1157   gStorage->Put(arrayNew, (*id1), metaData);
1158
1159   TFile falign("falign.root","recreate");
1160   arrayNew->Write("new");
1161   array->Write("old");
1162   falign.Close();
1163 }
1164
1165
1166
1167 void UpdateOCDBAlign0(){
1168   //
1169   // Store resulting OCDB entry
1170   // 0. Setup OCDB to get necccessary old entries - not done here
1171   // 1. Get old OCDB entry
1172   // 2. Get delta alignment
1173   // 3. Add delta alignment
1174   // 4. Store new alignment in
1175
1176   AliCDBEntry * entry = AliCDBManager::Instance()->Get("TPC/Align/Data");
1177   TClonesArray * array = (TClonesArray*)entry->GetObject();
1178   Int_t entries = array->GetEntries();
1179   TClonesArray * arrayNew=(TClonesArray*)array->Clone();
1180   TFile f("fitAlignCombined.root");
1181   TVectorD *vecDPhi = (TVectorD*) f.Get("vecDPhi");
1182   TVectorD *vecDGX = (TVectorD*) f.Get("vecDGX");
1183   TVectorD *vecDGY = (TVectorD*) f.Get("vecDGY");
1184
1185   //
1186   //
1187   { for (Int_t i=0;i<entries; i++){
1188       //
1189       //
1190       AliAlignObjParams *alignP = (AliAlignObjParams*)array->UncheckedAt(i);
1191       AliAlignObjParams *alignPNew = (AliAlignObjParams*)arrayNew->UncheckedAt(i);
1192       Int_t imod;
1193       AliGeomManager::ELayerID ilayer;
1194       alignP->GetVolUID(ilayer, imod);
1195       if (ilayer==AliGeomManager::kInvalidLayer) continue;
1196       Int_t sector=imod;
1197       if (ilayer==AliGeomManager::kTPC2) sector+=36;
1198       Double_t transOld[3], rotOld[3];
1199       alignP->GetTranslation(transOld);  // in cm
1200       alignP->GetAngles(rotOld);         // in degrees
1201       printf("%d\t%d\t%d\t",ilayer, imod,sector);
1202       printf("%f mrad \t%f mrad\n",1000*rotOld[2]*TMath::DegToRad(), 1000*((*vecDPhi)[sector%36]));
1203       alignPNew->SetRotation(rotOld[0],rotOld[1], rotOld[2]-(*vecDPhi)[sector%36]*TMath::RadToDeg());
1204       alignPNew->SetTranslation(transOld[0]-(*vecDGX)[sector%36],transOld[1]-(*vecDGY)[sector%36], transOld[2]);
1205       alignPNew->Print("kokot");
1206       alignP->Print("kokot");
1207     }
1208   }
1209
1210
1211
1212   TString ocdbStorage="local://"+gSystem->GetFromPipe("pwd")+"/OCDBUpdate";
1213   TString userName=gSystem->GetFromPipe("echo $USER");
1214   TString date=gSystem->GetFromPipe("date");
1215
1216   AliCDBMetaData *metaData= new AliCDBMetaData();
1217   metaData->SetObjectClassName("TObjArray");
1218   metaData->SetResponsible("Marian Ivanov");
1219   metaData->SetBeamPeriod(1);
1220   metaData->SetAliRootVersion("05-26-02"); //root version
1221   metaData->SetComment(Form("Correction calibration. User: %s\n Data%s",userName.Data(),date.Data()));
1222
1223   AliCDBId* id1=NULL;
1224   id1=new AliCDBId("TPC/Align/Data", 0, AliCDBRunRange::Infinity());
1225   AliCDBStorage* gStorage = AliCDBManager::Instance()->GetStorage(ocdbStorage);
1226   gStorage->Put(arrayNew, (*id1), metaData);
1227
1228   TFile falign("falign.root","recreate");
1229   arrayNew->Write("new");
1230   array->Write("old");
1231   falign.Close();
1232 }