TPC/CalibMacros/MakeGlobalFit.C

   1  /*
   2   gROOT->Macro("~/rootlogon.C")
   3   gSystem->AddIncludePath("-I$ALICE_ROOT/STAT")
   4   gSystem->AddIncludePath("-I$ALICE_ROOT/TPC")
   5   gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros")
   6   gSystem->Load("libANALYSIS");
   7   gSystem->Load("libTPCcalib");
   8
   9   .L $ALICE_ROOT/TPC/CalibMacros/MakeGlobalFit.C+
  10   MakeGlobalFit();
  11
  12   rm matching.ps
  13   aliroot -b -q $ALICE_ROOT/TPC/CalibMacros/MakeGlobalFit.C | tee globalFit.log
  14
  15 */
  16
  17 #if !defined(__CINT__) || defined(__MAKECINT__)
  18 #include "THnSparse.h"
  19 #include "TLatex.h"
  20 #include "TCanvas.h"
  21 #include "TLegend.h"
  22 #include "TSystem.h"
  23 #include "TFile.h"
  24 #include "TChain.h"
  25 #include "TCut.h"
  26 #include "TH3.h"
  27 #include "TH2F.h"
  28 #include "TProfile3D.h"
  29 #include "TMath.h"
  30 #include "TVectorD.h"
  31 #include "TMatrixD.h"
  32 #include "TStatToolkit.h"
  33 #include "TTreeStream.h"
  34 #include "AliExternalTrackParam.h"
  35 #include "AliESDfriend.h"
  36 #include "AliTPCcalibTime.h"
  37 #include "TROOT.h"
  38 #include "AliXRDPROOFtoolkit.h"
  39 #include "AliTPCCorrection.h"
  40 #include "AliTPCExBTwist.h"
  41 #include "AliTPCGGVoltError.h"
  42 #include "AliTPCComposedCorrection.h"
  43 #include "AliTPCExBConical.h"
  44 #include "TPostScript.h"
  45 #include "TStyle.h"
  46 #include "AliTrackerBase.h"
  47 #include "AliTPCCalibGlobalMisalignment.h"
  48 #include "AliTPCExBEffective.h"
  49 #include "TEntryList.h"
  50 #include "AliCDBMetaData.h"
  51 #include "AliCDBId.h"
  52 #include "AliCDBManager.h"
  53 #include "AliCDBStorage.h"
  54 #include "AliTPCcalibDB.h"
  55 #endif
  56
  57 const char* chDetName[5]={"ITS","TRD", "Vertex", "TOF","Laser"};
  58 const char* chParName[5]={"rphi","z", "snp", "tan","1/pt"};
  59 const char* chSideName[2]={"A side","C side"};
  60 Bool_t enableDet[5]={1,1,1,0,1};      // detector  for fit  ITS=0, TRD=1, Vertex=2, Laser=4
  61 Bool_t enableParam[5]={1,0,1,0,1};    //
  62 Bool_t enableSign=kFALSE;
  63 Bool_t useEff0=kFALSE;
  64 Bool_t useEffD=kFALSE;
  65 Bool_t useEffR=kFALSE;
  66 //
  67 TChain *chain    = 0;
  68 TChain *chainRef = 0;
  69 Bool_t printMatrix=kFALSE;
  70 TString *fitString=0;
  71 //
  72 // cuts
  73 //
  74 TCut cut="1";
  75 TCut cutD="1";
  76
  77 void MakeChain();
  78 void MakeAliases();
  79 void MakeCuts();
  80 void MakeFit(TCut cutCustom);
  81 void MakeOCDBEntry(Int_t refRun);
  82 TCanvas* DrawFitITS(const char *name);
  83 TCanvas* DrawFitVertex(const char *name);
  84 TCanvas*  DrawFitLaser(const char *cname);
  85 TCanvas* DrawCorrdY();
  86 TCanvas* DrawCorrdSnp();
  87 TCanvas * DrawFitdY(const char *name);
  88 TCanvas * DrawFitdSnp(const char *name);
  89 void PrintMatch();
  90 TCanvas * MakeComposedCorrection(const char * name);
  91
  92 void MakeAliases(){
  93   chain->SetAlias("isITS","(dtype==0||dtype==2)");
  94   chain->SetAlias("isTRD","(dtype==1)");
  95   chain->SetAlias("isLaser","(dtype==4)");
  96   chain->SetAlias("r","sqrt(gx*gx+gy*gy)");
  97   chain->SetAlias("r250","(sqrt(gx*gx+gy*gy)/250.)");
  98   chain->SetAlias("weight","((dtype==4)*rms*10+rms)");  // downscale laser
  99   chain->SetAlias("side","(1+(theta>0)*2)");
 100   chain->SetAlias("mdelta","(mean-R.mean-isLaser*((dRrec-R.dRrec)*tan(asin(snp))))");
 101   chain->SetAlias("drift","(1-abs(gz/250))");
 102   chain->SetAlias("r0","(r-160)/80");
 103   chain->SetAlias("delta","(0*gx)");
 104 }
 105
 106
 107 void MakeGlobalFit(){
 108   //
 109   //
 110   //
 111   gROOT->Macro("~/rootlogon.C");
 112   //gROOT->Macro("NimStyle.C");
 113   gSystem->AddIncludePath("-I$ALICE_ROOT/STAT");
 114   gSystem->AddIncludePath("-I$ALICE_ROOT/TPC");
 115   gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
 116   gSystem->Load("libANALYSIS");
 117   gSystem->Load("libTPCcalib");
 118
 119   MakeChain();
 120   gStyle->SetOptTitle(1);
 121   gStyle->SetOptStat(0);
 122   TPostScript *ps = new TPostScript("exbFit.eps", 112);
 123   ps->NewPage();
 124   TCanvas *c=0;
 125   useEff0=kFALSE; useEffD=kFALSE; useEffR=kFALSE;
 126   ps->NewPage();
 127   MakeFit("1");
 128   //
 129   ps->NewPage();
 130   c=DrawFitdY("dY-Physical");
 131   c->Update();
 132   ps->NewPage();
 133   c->Update();
 134
 135   ps->NewPage();
 136   c=DrawFitdSnp("dSnp-Physical");
 137   c->Update();
 138
 139   ps->NewPage();
 140   c=DrawFitITS("ITS Physical");
 141   c->Update();
 142
 143   ps->NewPage();
 144   c=DrawFitVertex("Vertex Physical");
 145   c->Update();
 146
 147   ps->NewPage();
 148   c=DrawFitLaser("Laser Physical");
 149   c->Update();
 150
 151   ps->NewPage();
 152   c=MakeComposedCorrection("Correction physical");
 153   c->Update();
 154
 155   //
 156   //
 157   //
 158   //
 159   useEff0=kTRUE; useEffD=kTRUE; useEffR=kTRUE; enableSign=kFALSE;
 160   ps->NewPage();
 161   MakeFit("1");
 162   ps->NewPage();
 163   c=DrawFitdY("dY-Physical+Effective ");
 164   c->Update();
 165
 166   ps->NewPage();
 167   c=DrawFitdSnp("dSnp-Physical+Effective ");
 168   c->Update();
 169
 170   ps->NewPage();
 171   c=DrawFitITS("ITS Physical+Effective");
 172   c->Update();
 173
 174   ps->NewPage();
 175   c=DrawFitVertex("Vertex Physical+Effective");
 176   c->Update();
 177
 178   ps->NewPage();
 179   c=DrawFitLaser("Laser Physical +Effective ");
 180   c->Update();
 181
 182   ps->NewPage();
 183   c=MakeComposedCorrection("Correction physical+Effective");
 184   c->Update();
 185   //
 186   useEff0=kTRUE; useEffD=kTRUE; useEffR=kTRUE; enableSign=kTRUE;
 187   ps->NewPage();
 188   MakeFit("1");
 189   //
 190   ps->NewPage();
 191   c=DrawFitdY("dY-Physical+Effective Sign");
 192   c->Update();
 193
 194   ps->NewPage();
 195   c=DrawFitdSnp("dSnp-Physical+Effective Sign");
 196   c->Update();
 197
 198   ps->NewPage();
 199   c=DrawFitITS("ITS Physical+Effective Sign");
 200   c->Update();
 201
 202   ps->NewPage();
 203   c=DrawFitVertex("Vertex Physical+Effective Sign");
 204   c->Update();
 205
 206   ps->NewPage();
 207   c=DrawFitLaser("Laser Physical +Effective Sign");
 208   c->Update();
 209
 210   ps->NewPage();
 211   c=MakeComposedCorrection("Correction physical+Effective Sign");
 212   c->Update();
 213
 214   ps->Close();
 215   delete ps;
 216
 217   //
 218 }
 219
 220 void MakeChain(){
 221   //
 222   //
 223   TH1::AddDirectory(0);
 224   TFile * f0 =0;      // file 0 field
 225   TFile * fp =0;      // file plus
 226   TFile * fm =0;      // file minus
 227   TTree * tree0=0;
 228   TTree * treeP=0;
 229   TTree * treeM=0;
 230   //
 231   chain    = new TChain("fit","fit");
 232   chainRef = new TChain("fit","fit");
 233   //
 234   //
 235   //
 236   for (Int_t idet=0; idet<5; idet++){
 237     for (Int_t ipar=0; ipar<5; ipar++){
 238       f0= TFile::Open(Form("../mergeField0/distortion%d_%d.root",idet,ipar));
 239       fp= TFile::Open(Form("../mergePlus/distortion%d_%d.root",idet,ipar));
 240       fm= TFile::Open(Form("../mergeMinus/distortion%d_%d.root",idet,ipar));
 241       tree0 = (f0) ? (TTree*)f0->Get("fit"):0;
 242       treeP = (fp) ? (TTree*)fp->Get("fit"):0;
 243       treeM = (fm) ? (TTree*)fm->Get("fit"):0;
 244       //
 245       if ( ipar==0 || ipar==2){
 246         if (tree0 && treeP){
 247           chain->Add(Form("../mergePlus/distortion%d_%d.root",idet,ipar));
 248           chainRef->Add(Form("../mergeField0/distortion%d_%d.root",idet,ipar));
 249         }
 250         if (tree0 && treeM){
 251           chain->Add(Form("../mergeMinus/distortion%d_%d.root",idet,ipar));
 252           chainRef->Add(Form("../mergeField0/distortion%d_%d.root",idet,ipar));
 253         }
 254       }
 255       if ( ipar==1 || ipar==3 || ipar==4){
 256         if (treeP && treeM){
 257           chain->Add(Form("../mergePlus/distortion%d_%d.root",idet,ipar));
 258           chainRef->Add(Form("../mergeMinus/distortion%d_%d.root",idet,ipar));
 259         }
 260       }
 261     }
 262   }
 263   chain->AddFriend(chainRef,"R");
 264   MakeAliases();
 265   MakeCuts();
 266 }
 267
 268
 269 void MakeCuts(){
 270   //
 271   //
 272   //
 273   TCut cutS="((rms>0&&R.rms>0&&entries>0&&R.entries>0))";         // statistic cuts
 274   TCut cutType="((dtype==R.dtype)&&(ptype==R.ptype))";            // corresponding types
 275   TCut cutOut="(ptype==0)*abs(mdelta)<(0.3+rms)||(ptype==0&&abs(mdelta*85)<(0.3+rms*85))";            // corresponding types
 276   //
 277
 278   chain->Draw(">>outList",cutS+cutType+cutOut+"abs(snp)<0.5","entryList");
 279   TEntryList *elistOut = (TEntryList*)gDirectory->Get("outList");
 280   chain->SetEntryList(elistOut);
 281
 282
 283 }
 284
 285
 286 TMatrixD * MakeCorrelation(TMatrixD &matrix){
 287   //
 288   //
 289   //
 290   Int_t nrows = matrix.GetNrows();
 291   TMatrixD * mat = new TMatrixD(nrows,nrows);
 292   for (Int_t irow=0; irow<nrows; irow++)
 293     for (Int_t icol=0; icol<nrows; icol++){
 294       (*mat)(irow,icol)= matrix(irow,icol)/TMath::Sqrt(matrix(irow,irow)*matrix(icol,icol));
 295     }
 296   return mat;
 297 }
 298
 299 void MakeDetCut(){
 300   cutD=Form("((dtype==%d)||(dtype==%d)||(dtype==%d)||(dtype==%d)||(dtype==%d))",enableDet[0]?0:0,enableDet[1]?1:0,enableDet[2]?2:0,enableDet[3]?3:0,enableDet[4]?4:0);
 301   cutD+=Form("((ptype==%d)||(ptype==%d)||(ptype==%d)||(ptype==%d)||(ptype==%d))",enableParam[0]?0:0,enableParam[1]?1:0,enableParam[2]?2:0,enableParam[3]?3:0,enableParam[4]?4:0);
 302 }
 303
 304 void MakeFit(TCut cutCustom){
 305   MakeDetCut();
 306   Int_t  npointsMax=30000000;
 307   TStatToolkit toolkit;
 308   Double_t chi2=0;
 309   Int_t    npoints=0;
 310   TVectorD param;
 311   TMatrixD covar;
 312   //
 313   TString  fstring="";                       // magnetic part
 314   fstring+="(tX1-R.tX1)++";                  // twist X
 315   fstring+="(tY1-R.tY1)++";                   // twist Y
 316   fstring+="(sign(bz)*(tX1-R.tX1))++";       // twist X
 317   fstring+="(sign(bz)*(tY1-R.tY1))++";       // twist Y
 318
 319   {if (enableDet[0]){
 320     fstring+="(isITS*shiftX)++";             // shift X - ITS
 321     fstring+="(isITS*shiftY)++";             // shift Y - ITS
 322     fstring+="(isITS*sign(bz)*shiftX)++";    // shift X - ITS
 323     fstring+="(isITS*sign(bz)*shiftY)++";    // shift Y - ITS
 324     }}
 325   {if (enableDet[1]){
 326     fstring+="(isTRD*shiftX)++";             // shift X - TRD
 327     fstring+="(isTRD*shiftY)++";             // shift Y - TRD
 328     fstring+="(isTRD*sign(bz)*shiftX)++";    // shift X - TRD
 329     fstring+="(isTRD*sign(bz)*shiftY)++";    // shift Y - TRD
 330     }}
 331   //
 332   if (enableDet[4]){
 333   }
 334   TString fstringCustom="";
 335   if (useEff0){
 336     fstring+="(exbT1-exb11-(R.exbT1-R.exb11))++";                // T1 adjustment
 337     fstring+="(exbT2-exb11-(R.exbT2-R.exb11))++";                // T2 adjustment
 338     fstringCustom+="(eff0_0_0-R.eff0_0_0)++";                  // effective correction constant part
 339     fstringCustom+="(eff1_0_0-R.eff1_0_0)++";                  //
 340   }
 341   if (useEffD){
 342     fstringCustom+="(eff0_0_1-R.eff0_0_1)++";                  // effective correction drift part
 343     fstringCustom+="(eff1_0_1-R.eff1_0_1)++";                  //
 344     fstringCustom+="(eff0_0_2-R.eff0_0_2)++";                  // effective correction drift part
 345     fstringCustom+="(eff1_0_2-R.eff1_0_2)++";                  //
 346   }
 347   if (useEffR){
 348     fstringCustom+="(eff0_1_0-R.eff0_1_0)++";                  // effective correction radial part
 349     fstringCustom+="(eff1_1_0-R.eff1_1_0)++";                  //
 350     fstringCustom+="(eff0_1_1-R.eff0_1_1)++";                  // effective correction radial part
 351     fstringCustom+="(eff1_1_1-R.eff1_1_1)++";                  //
 352     fstringCustom+="(eff0_1_2-R.eff0_1_2)++";                  // effective correction radial part
 353     fstringCustom+="(eff1_1_2-R.eff1_1_2)++";                  //
 354     fstringCustom+="(eff0_2_0-R.eff0_2_0)++";                  // effective correction radial part
 355     fstringCustom+="(eff1_2_0-R.eff1_2_0)++";                  //
 356     fstringCustom+="(eff0_2_1-R.eff0_2_1)++";                  // effective correction radial part
 357     fstringCustom+="(eff1_2_1-R.eff1_2_1)++";                  //
 358     fstringCustom+="(eff0_2_2-R.eff0_2_2)++";                  // effective correction radial part
 359     fstringCustom+="(eff1_2_2-R.eff1_2_2)++";                  //
 360   }
 361   if (enableSign){
 362     if (useEff0){
 363       fstringCustom+="sign(bz)*(eff0_0_0-R.eff0_0_0)++";                  // effective correction constant part
 364       fstringCustom+="sign(bz)*(eff1_0_0-R.eff1_0_0)++";                  //
 365     }
 366     if (useEffD){
 367       fstringCustom+="sign(bz)*(eff0_0_1-R.eff0_0_1)++";                  // effective correction drift part
 368       fstringCustom+="sign(bz)*(eff1_0_1-R.eff1_0_1)++";                  //
 369       fstringCustom+="sign(bz)*(eff0_0_2-R.eff0_0_2)++";                  // effective correction drift part
 370       fstringCustom+="sign(bz)*(eff1_0_2-R.eff1_0_2)++";                  //
 371     }
 372     if (useEffR){
 373       fstringCustom+="sign(bz)*(eff0_1_0-R.eff0_1_0)++";                  // effective correction radial part
 374       fstringCustom+="sign(bz)*(eff1_1_0-R.eff1_1_0)++";                  //
 375       fstringCustom+="sign(bz)*(eff0_1_1-R.eff0_1_1)++";                  // effective correction radial part
 376       fstringCustom+="sign(bz)*(eff1_1_1-R.eff1_1_1)++";                  //
 377       fstringCustom+="sign(bz)*(eff0_1_2-R.eff0_1_2)++";                  // effective correction radial part
 378       fstringCustom+="sign(bz)*(eff1_1_2-R.eff1_1_2)++";                  //
 379       fstringCustom+="sign(bz)*(eff0_2_0-R.eff0_2_0)++";                  // effective correction radial part
 380       fstringCustom+="sign(bz)*(eff1_2_0-R.eff1_2_0)++";                  //
 381       fstringCustom+="sign(bz)*(eff0_2_1-R.eff0_2_1)++";                  // effective correction radial part
 382       fstringCustom+="sign(bz)*(eff1_2_1-R.eff1_2_1)++";                  //
 383       fstringCustom+="sign(bz)*(eff0_2_2-R.eff0_2_2)++";                  // effective correction radial part
 384       fstringCustom+="sign(bz)*(eff1_2_2-R.eff1_2_2)++";                  //
 385     }
 386   }
 387   //
 388   fstring+=fstringCustom;
 389   //
 390   TString *strDelta = TStatToolkit::FitPlane(chain,"mdelta:weight", fstring.Data(),cut+cutD+cutCustom, chi2,npoints,param,covar,-1,0, npointsMax, kTRUE);
 391   if (printMatrix) MakeCorrelation(covar)->Print();
 392   chain->SetAlias("delta",strDelta->Data());
 393   strDelta->Tokenize("++")->Print();
 394   fitString = strDelta;
 395   PrintMatch();
 396 }
 397
 398
 399
 400
 401
 402 void PrintMatch(){
 403   //
 404   // Print detector matching info
 405   //
 406   for (Int_t ipar=0; ipar<5; ipar++){
 407     for (Int_t idet=0; idet<5; idet++){
 408       Double_t mean0,rms0,mean1,rms1;
 409       Int_t entries = chain->Draw("delta-mdelta>>rhis",cut+Form("dtype==%d&&ptype==%d",idet,ipar),"goff");
 410       if (entries==0) continue;
 411       TH1 * his = (TH1*)(chain->GetHistogram()->Clone());
 412       mean1=his->GetMean();
 413       rms1 =his->GetRMS();
 414       delete his;
 415       //
 416       entries = chain->Draw("mdelta>>rhis",cut+Form("dtype==%d&&ptype==%d",idet,ipar),"goff");
 417       if (entries==0) continue;
 418       his = (TH1*)(chain->GetHistogram()->Clone());
 419       //
 420       mean0=his->GetMean();
 421       rms0 =his->GetRMS();
 422
 423       printf("ptype==%s\tdtype==%s\tMean=%f -> %f\tRMS=%f -> %f\n",chParName[ipar],chDetName[idet], mean0,mean1,rms0,rms1);
 424       delete his;
 425     }
 426   }
 427 }
 428
 429
 430
 431
 432 TCanvas* DrawFitITS(const char *name){
 433   //
 434   //
 435   //
 436   TLegend *legend=0;
 437   TCanvas *canvas = new TCanvas(name,name,800,800);
 438   canvas->Divide(1,2);
 439   Int_t entries=0;
 440   const char * grname[10]={"A side (B=0.5T)","A side (B=-0.5T)", "C side (B=0.5T)", "C side (B=-0.5T)",0};
 441
 442   TGraph *graphsdyITS[10];
 443   TGraph *graphsdyITSc[10];
 444   entries=chain->Draw("mdelta:phi",cut+"ptype==0&&dtype==0&&bz>0&&theta>0","goff");
 445   graphsdyITS[0]=new TGraph(entries,chain->GetV2(), chain->GetV1());
 446   entries=chain->Draw("mdelta:phi",cut+cut+"ptype==0&&dtype==0&&bz<0&&theta>0","goff");
 447   graphsdyITS[1]=new TGraph(entries,chain->GetV2(), chain->GetV1());
 448   entries=chain->Draw("mdelta:phi",cut+cut+"ptype==0&&dtype==0&&bz>0&&theta<0","goff");
 449   graphsdyITS[2]=new TGraph(entries,chain->GetV2(), chain->GetV1());
 450   entries=chain->Draw("mdelta:phi",cut+cut+"ptype==0&&dtype==0&&bz<0&&theta<0","goff");
 451   graphsdyITS[3]=new TGraph(entries,chain->GetV2(), chain->GetV1());
 452   //
 453   entries=chain->Draw("mdelta-delta:phi",cut+cut+"ptype==0&&dtype==0&&bz>0&&theta>0","goff");
 454   graphsdyITSc[0]=new TGraph(entries,chain->GetV2(), chain->GetV1());
 455   entries=chain->Draw("mdelta-delta:phi",cut+cut+"ptype==0&&dtype==0&&bz<0&&theta>0","goff");
 456   graphsdyITSc[1]=new TGraph(entries,chain->GetV2(), chain->GetV1());
 457   entries=chain->Draw("mdelta-delta:phi",cut+cut+"ptype==0&&dtype==0&&bz>0&&theta<0","goff");
 458   graphsdyITSc[2]=new TGraph(entries,chain->GetV2(), chain->GetV1());
 459   entries=chain->Draw("mdelta-delta:phi",cut+cut+"ptype==0&&dtype==0&&bz<0&&theta<0","goff");
 460   graphsdyITSc[3]=new TGraph(entries,chain->GetV2(), chain->GetV1());
 461
 462   canvas->cd(1);
 463   legend = new TLegend(0.6,0.6,1.0,1.0,"#Delta R#phi TPC-ITS");
 464   for (Int_t i=0; i<4; i++){
 465     graphsdyITS[i]->SetMaximum(0.5);
 466     graphsdyITS[i]->SetMinimum(-0.5);
 467     graphsdyITS[i]->SetMarkerColor(i+1);
 468     graphsdyITS[i]->SetMarkerStyle(25+i);
 469     //    graphsdyITS[i]->SetName(grname[i]);
 470     graphsdyITS[i]->GetXaxis()->SetTitle("#phi");
 471     graphsdyITS[i]->GetYaxis()->SetTitle("#Delta R#Phi (cm)");
 472     if (i==0) graphsdyITS[i]->Draw("ap");
 473     graphsdyITS[i]->Draw("p");
 474     legend->AddEntry(graphsdyITS[i],grname[i]);
 475   }
 476   legend->Draw();
 477   canvas->cd(2);
 478   legend = new TLegend(0.6,0.6,1.0,1.0,"#Delta R#phi TPC-ITS corrected");
 479   for (Int_t i=0; i<4; i++){
 480     graphsdyITSc[i]->SetMaximum(0.5);
 481     graphsdyITSc[i]->SetMinimum(-0.5);
 482     graphsdyITSc[i]->SetMarkerColor(i+1);
 483     graphsdyITSc[i]->SetMarkerStyle(25+i);
 484     //    graphsdyITS[i]->SetName(grname[i]);
 485     graphsdyITSc[i]->GetXaxis()->SetTitle("#phi");
 486     graphsdyITSc[i]->GetYaxis()->SetTitle("#Delta R#Phi (cm)");
 487     if (i==0) graphsdyITS[i]->Draw("ap");
 488     graphsdyITSc[i]->Draw("p");
 489     legend->AddEntry(graphsdyITSc[i],grname[i]);
 490   }
 491   legend->Draw();
 492   return canvas;
 493 }
 494
 495
 496 TCanvas*  DrawFitLaser(const char *cname){
 497   //
 498   //
 499   //
 500   TH1::AddDirectory(0);
 501   TCut cutLaser=cut+"isLaser&&bz<0";
 502   TCanvas *canvas= new TCanvas(cname, cname,800,800);
 503   canvas->Divide(2,2);
 504   canvas->Draw();
 505   TLegend*legend=0;
 506   //
 507   TH1 *his[16]={0};
 508   {for (Int_t icorr=0; icorr<2; icorr++){
 509     for (Int_t iside=0; iside<2; iside++){
 510       canvas->cd(iside+1);
 511       for (Int_t ib=0; ib<4; ib++){
 512         TCut cutB="1";
 513         if (iside==0) cutB=Form("(id==%d&&gz>0)",ib);
 514         if (iside==1) cutB=Form("(id==%d&&gz<0)",ib);
 515         //cutB.Print();
 516         if (icorr==0) chain->Draw("10*mdelta:r",cutLaser+cutB,"prof");
 517         if (icorr==1) chain->Draw("10*(mdelta-delta):r",cutLaser+cutB,"prof");
 518         his[icorr*8+iside*4+ib] = (TH1*)(chain->GetHistogram()->Clone());
 519         his[icorr*8+iside*4+ib]->SetName(Form("B%d%d%d",icorr,iside,ib));
 520         his[icorr*8+iside*4+ib]->SetTitle(Form("Bundle %d",ib));
 521         his[icorr*8+iside*4+ib]->SetMarkerColor(ib+1);
 522         his[icorr*8+iside*4+ib]->SetMarkerStyle(ib+25);
 523         his[icorr*8+iside*4+ib]->SetMarkerSize(0.4);
 524         his[icorr*8+iside*4+ib]->SetMaximum(3);
 525         his[icorr*8+iside*4+ib]->SetMinimum(-3);
 526         his[icorr*8+iside*4+ib]->GetXaxis()->SetTitle("r (cm)");
 527         his[icorr*8+iside*4+ib]->GetYaxis()->SetTitle("#Delta r#phi (mm)");
 528       }
 529     }
 530     }}
 531     //
 532   for (Int_t icorr=0; icorr<2; icorr++){
 533     for (Int_t iside=0; iside<2; iside++){
 534       canvas->cd(icorr*2+iside+1);
 535       legend = new TLegend(0.6,0.6,1.0,1.0,Form("#Delta R#phi Laser-%s",chSideName[iside]));
 536       for (Int_t ib=0; ib<4; ib++){
 537         if (ib==0) his[icorr*2+iside*4+ib]->Draw();
 538         his[icorr*8+iside*4+ib]->Draw("same");
 539         legend->AddEntry(his[icorr*8+iside*4+ib]);
 540       }
 541       legend->Draw();
 542     }
 543   }
 544   return canvas;
 545 }
 546
 547
 548
 549
 550
 551
 552
 553
 554
 555
 556
 557 TCanvas* DrawFitVertex(const char *name){
 558   //
 559   //
 560   //
 561   TLegend *legend=0;
 562   TCanvas *canvas = new TCanvas(name,name,800,800);
 563   canvas->Divide(1,2);
 564   Int_t entries=0;
 565   const char * grname[10]={"A side (B=0.5T)","A side (B=-0.5T)", "C side (B=0.5T)", "C side (B=-0.5T)",0};
 566
 567   TGraph *graphsdyVertex[10];
 568   TGraph *graphsdyVertexc[10];
 569   entries=chain->Draw("mdelta:phi",cut+cut+"ptype==0&&dtype==2&&bz>0&&theta>0","goff");
 570   graphsdyVertex[0]=new TGraph(entries,chain->GetV2(), chain->GetV1());
 571   entries=chain->Draw("mdelta:phi",cut+cut+"ptype==0&&dtype==2&&bz<0&&theta>0","goff");
 572   graphsdyVertex[1]=new TGraph(entries,chain->GetV2(), chain->GetV1());
 573   entries=chain->Draw("mdelta:phi",cut+cut+"ptype==0&&dtype==2&&bz>0&&theta<0","goff");
 574   graphsdyVertex[2]=new TGraph(entries,chain->GetV2(), chain->GetV1());
 575   entries=chain->Draw("mdelta:phi",cut+cut+"ptype==0&&dtype==2&&bz<0&&theta<0","goff");
 576   graphsdyVertex[3]=new TGraph(entries,chain->GetV2(), chain->GetV1());
 577   //
 578   entries=chain->Draw("mdelta-delta:phi",cut+cut+"ptype==0&&dtype==2&&bz>0&&theta>0","goff");
 579   graphsdyVertexc[0]=new TGraph(entries,chain->GetV2(), chain->GetV1());
 580   entries=chain->Draw("mdelta-delta:phi",cut+cut+"ptype==0&&dtype==2&&bz<0&&theta>0","goff");
 581   graphsdyVertexc[1]=new TGraph(entries,chain->GetV2(), chain->GetV1());
 582   entries=chain->Draw("mdelta-delta:phi",cut+cut+"ptype==0&&dtype==2&&bz>0&&theta<0","goff");
 583   graphsdyVertexc[2]=new TGraph(entries,chain->GetV2(), chain->GetV1());
 584   entries=chain->Draw("mdelta-delta:phi",cut+cut+"ptype==0&&dtype==2&&bz<0&&theta<0","goff");
 585   graphsdyVertexc[3]=new TGraph(entries,chain->GetV2(), chain->GetV1());
 586
 587   canvas->cd(1);
 588   legend = new TLegend(0.6,0.6,1.0,1.0,"#Delta R#phi TPC-Vertex");
 589   for (Int_t i=0; i<4; i++){
 590     graphsdyVertex[i]->SetMaximum(1);
 591     graphsdyVertex[i]->SetMinimum(-1);
 592     graphsdyVertex[i]->SetMarkerColor(i+1);
 593     graphsdyVertex[i]->SetMarkerStyle(25+i);
 594     //    graphsdyVertex[i]->SetName(grname[i]);
 595     graphsdyVertex[i]->GetXaxis()->SetTitle("#phi");
 596     graphsdyVertex[i]->GetYaxis()->SetTitle("#Delta R#Phi (cm)");
 597     if (i==0) graphsdyVertex[i]->Draw("ap");
 598     graphsdyVertex[i]->Draw("p");
 599     legend->AddEntry(graphsdyVertex[i],grname[i]);
 600   }
 601   legend->Draw();
 602   canvas->cd(2);
 603   legend = new TLegend(0.6,0.6,1.0,1.0,"#Delta R#phi TPC-Vertex corrected");
 604   for (Int_t i=0; i<4; i++){
 605     graphsdyVertexc[i]->SetMaximum(1);
 606     graphsdyVertexc[i]->SetMinimum(-1);
 607     graphsdyVertexc[i]->SetMarkerColor(i+1);
 608     graphsdyVertexc[i]->SetMarkerStyle(25+i);
 609     //    graphsdyVertex[i]->SetName(grname[i]);
 610     graphsdyVertexc[i]->GetXaxis()->SetTitle("#phi");
 611     graphsdyVertexc[i]->GetYaxis()->SetTitle("#Delta R#Phi (cm)");
 612     if (i==0) graphsdyVertex[i]->Draw("ap");
 613     graphsdyVertexc[i]->Draw("p");
 614     legend->AddEntry(graphsdyVertexc[i],grname[i]);
 615   }
 616   legend->Draw();
 617   return canvas;
 618 }
 619
 620
 621
 622
 623
 624
 625 TCanvas* DrawCorrdY(){
 626
 627   TLegend *legend=0;
 628   TCanvas *canvas = new TCanvas("Corrections","Corrections",800,800);
 629   canvas->Divide(1,3);
 630
 631   TGraph *graphsITS[10];
 632   TGraph *graphsTRD[10];
 633   TGraph *graphsVertex[10];
 634   Int_t entries=0;
 635   const char * grname[10]={"X twist (1 mrad)","Y twist (1 mrad)", "X shift (1 mm)", "Y shift (1 mm)", "drPhi (1 mm)"};
 636   //
 637   entries=chain->Draw("tX1:phi",cut+cut+"ptype==0&&dtype==0&&bz>0","goff");
 638   graphsITS[0]=new TGraph(entries,chain->GetV2(), chain->GetV1());
 639   entries=chain->Draw("tY1:phi",cut+cut+"ptype==0&&dtype==0&&bz>0","goff");
 640   graphsITS[1]=new TGraph(entries,chain->GetV2(), chain->GetV1());
 641   entries=chain->Draw("shiftX:phi",cut+cut+"ptype==0&&dtype==0&&bz>0","goff");
 642   graphsITS[2]=new TGraph(entries,chain->GetV2(), chain->GetV1());
 643   entries=chain->Draw("shiftY:phi",cut+cut+"ptype==0&&dtype==0&&bz>0","goff");
 644   graphsITS[3]=new TGraph(entries,chain->GetV2(), chain->GetV1());
 645   entries=chain->Draw("drPhiA+drPhiC:phi",cut+cut+"ptype==0&&dtype==0&&bz>0","goff");
 646   graphsITS[4]=new TGraph(entries,chain->GetV2(), chain->GetV1());
 647   //
 648   entries=chain->Draw("tX1:phi",cut+cut+"ptype==0&&dtype==1&&bz>0","goff");
 649   graphsTRD[0]=new TGraph(entries,chain->GetV2(), chain->GetV1());
 650   entries=chain->Draw("tY1:phi",cut+cut+"ptype==0&&dtype==1&&bz>0","goff");
 651   graphsTRD[1]=new TGraph(entries,chain->GetV2(), chain->GetV1());
 652   entries=chain->Draw("shiftX:phi",cut+cut+"ptype==0&&dtype==1&&bz>0","goff");
 653   graphsTRD[2]=new TGraph(entries,chain->GetV2(), chain->GetV1());
 654   entries=chain->Draw("shiftY:phi",cut+cut+"ptype==0&&dtype==1&&bz>0","goff");
 655   graphsTRD[3]=new TGraph(entries,chain->GetV2(), chain->GetV1());
 656   entries=chain->Draw("drPhiA+drPhiC:phi",cut+cut+"ptype==0&&dtype==1&&bz>0","goff");
 657   graphsTRD[4]=new TGraph(entries,chain->GetV2(), chain->GetV1());
 658   //
 659   entries=chain->Draw("tX1:phi",cut+cut+"ptype==0&&dtype==2&&bz>0","goff");
 660   graphsVertex[0]=new TGraph(entries,chain->GetV2(), chain->GetV1());
 661   entries=chain->Draw("tY1:phi",cut+cut+"ptype==0&&dtype==2&&bz>0","goff");
 662   graphsVertex[1]=new TGraph(entries,chain->GetV2(), chain->GetV1());
 663   entries=chain->Draw("shiftX:phi",cut+cut+"ptype==0&&dtype==2&&bz>0","goff");
 664   graphsVertex[2]=new TGraph(entries,chain->GetV2(), chain->GetV1());
 665   entries=chain->Draw("shiftY:phi",cut+cut+"ptype==0&&dtype==2&&bz>0","goff");
 666   graphsVertex[3]=new TGraph(entries,chain->GetV2(), chain->GetV1());
 667   entries=chain->Draw("drPhiA+drPhiC:phi",cut+cut+"ptype==0&&dtype==2&&bz>0","goff");
 668   graphsVertex[4]=new TGraph(entries,chain->GetV2(), chain->GetV1());
 669
 670   canvas->cd(1);
 671   //
 672   legend = new TLegend(0.6,0.6,1.0,1.0,"#Delta R#Phi TPC-Vertex");
 673   for (Int_t i=0; i<5; i++){
 674     graphsVertex[i]->SetMaximum(0.2);
 675     graphsVertex[i]->SetMinimum(-0.2);
 676     graphsVertex[i]->SetMarkerColor(i+1);
 677     graphsVertex[i]->SetMarkerStyle(25+i);
 678     graphsVertex[i]->SetName(grname[i]);
 679     graphsVertex[i]->GetXaxis()->SetTitle("#phi");
 680     graphsVertex[i]->GetYaxis()->SetTitle("#DeltaR#Phi (cm)");
 681     if (i==0) graphsVertex[i]->Draw("ap");
 682     graphsVertex[i]->Draw("p");
 683     legend->AddEntry(graphsITS[i],grname[i]);
 684   }
 685   legend->Draw();
 686
 687   canvas->cd(2);
 688   legend = new TLegend(0.6,0.6,1.0,1.0,"#Delta R#Phi ITS-TPC");
 689   for (Int_t i=0; i<5; i++){
 690     graphsITS[i]->SetMaximum(0.2);
 691     graphsITS[i]->SetMinimum(-0.2);
 692     graphsITS[i]->SetMarkerColor(i+1);
 693     graphsITS[i]->SetMarkerStyle(25+i);
 694     graphsITS[i]->SetName(grname[i]);
 695     graphsITS[i]->GetXaxis()->SetTitle("#phi");
 696     graphsITS[i]->GetYaxis()->SetTitle("#Delta R#Phi (cm)");
 697     if (i==0) graphsITS[i]->Draw("ap");
 698     graphsITS[i]->Draw("p");
 699     legend->AddEntry(graphsITS[i],grname[i]);
 700   }
 701   legend->Draw();
 702   //
 703   canvas->cd(3);
 704   legend = new TLegend(0.6,0.6,1.0,1.0,"#Delta R#Phi TPC-TRD");
 705   for (Int_t i=0; i<5; i++){
 706     graphsTRD[i]->SetMaximum(0.2);
 707     graphsTRD[i]->SetMinimum(-0.2);
 708     graphsTRD[i]->SetMarkerColor(i+1);
 709     graphsTRD[i]->SetMarkerStyle(25+i);
 710     graphsTRD[i]->SetName(grname[i]);
 711     graphsTRD[i]->GetXaxis()->SetTitle("#phi");
 712     graphsTRD[i]->GetYaxis()->SetTitle("#DeltaR#Phi (cm)");
 713     if (i==0) graphsTRD[i]->Draw("ap");
 714     graphsTRD[i]->Draw("p");
 715     legend->AddEntry(graphsITS[i],grname[i]);
 716   }
 717   legend->Draw();
 718   return canvas;
 719 }
 720
 721
 722 TCanvas * DrawCorrdSnp(){
 723
 724   TLegend *legend=0;
 725   TCanvas *canvas = new TCanvas("Corrections dSnp","Corrections dSnp",800,800);
 726   canvas->Divide(1,3);
 727
 728   TGraph *graphsITS[10];
 729   TGraph *graphsTRD[10];
 730   TGraph *graphsVertex[10];
 731   Int_t entries=0;
 732   const char * grname[10]={"X twist (1 mrad)","Y twist (1 mrad)", "X shift (1 mm)", "Y shift (1 mm)",0};
 733   //
 734   entries=chain->Draw("1000*tX1:phi",cut+cut+"ptype==2&&dtype==0&&bz>0","goff");
 735   graphsITS[0]=new TGraph(entries,chain->GetV2(), chain->GetV1());
 736   entries=chain->Draw("1000*tY1:phi",cut+cut+"ptype==2&&dtype==0&&bz>0","goff");
 737   graphsITS[1]=new TGraph(entries,chain->GetV2(), chain->GetV1());
 738   entries=chain->Draw("1000*shiftX:phi",cut+cut+"ptype==2&&dtype==0&&bz>0","goff");
 739   graphsITS[2]=new TGraph(entries,chain->GetV2(), chain->GetV1());
 740   entries=chain->Draw("1000*shiftY:phi",cut+cut+"ptype==2&&dtype==0&&bz>0","goff");
 741   graphsITS[3]=new TGraph(entries,chain->GetV2(), chain->GetV1());
 742   //
 743   entries=chain->Draw("1000*tX1:phi",cut+cut+"ptype==2&&dtype==1&&bz>0","goff");
 744   graphsTRD[0]=new TGraph(entries,chain->GetV2(), chain->GetV1());
 745   entries=chain->Draw("1000*tY1:phi",cut+cut+"ptype==2&&dtype==1&&bz>0","goff");
 746   graphsTRD[1]=new TGraph(entries,chain->GetV2(), chain->GetV1());
 747   entries=chain->Draw("1000*shiftX:phi",cut+cut+"ptype==2&&dtype==1&&bz>0","goff");
 748   graphsTRD[2]=new TGraph(entries,chain->GetV2(), chain->GetV1());
 749   entries=chain->Draw("1000*shiftY:phi",cut+cut+"ptype==2&&dtype==1&&bz>0","goff");
 750   graphsTRD[3]=new TGraph(entries,chain->GetV2(), chain->GetV1());
 751   //
 752   entries=chain->Draw("1000*tX1:phi",cut+cut+"ptype==2&&dtype==2&&bz>0","goff");
 753   graphsVertex[0]=new TGraph(entries,chain->GetV2(), chain->GetV1());
 754   entries=chain->Draw("1000*tY1:phi",cut+cut+"ptype==2&&dtype==2&&bz>0","goff");
 755   graphsVertex[1]=new TGraph(entries,chain->GetV2(), chain->GetV1());
 756   entries=chain->Draw("1000*shiftX:phi",cut+cut+"ptype==2&&dtype==2&&bz>0","goff");
 757   graphsVertex[2]=new TGraph(entries,chain->GetV2(), chain->GetV1());
 758   entries=chain->Draw("1000*shiftY:phi",cut+cut+"ptype==2&&dtype==2&&bz>0","goff");
 759   graphsVertex[3]=new TGraph(entries,chain->GetV2(), chain->GetV1());
 760
 761   canvas->cd(1);
 762   //
 763   legend = new TLegend(0.6,0.6,1.0,1.0,"#Delta sin(#phi) TPC-Vertex");
 764   for (Int_t i=0; i<4; i++){
 765     graphsVertex[i]->SetMaximum(2);
 766     graphsVertex[i]->SetMinimum(2);
 767     graphsVertex[i]->SetMarkerColor(i+1);
 768     graphsVertex[i]->SetMarkerStyle(25+i);
 769     graphsVertex[i]->SetName(grname[i]);
 770     graphsVertex[i]->GetXaxis()->SetTitle("#phi");
 771     graphsVertex[i]->GetYaxis()->SetTitle("#Deltasin(#phi) (mrad)");
 772     if (i==0) graphsVertex[i]->Draw("ap");
 773     graphsVertex[i]->Draw("p");
 774     legend->AddEntry(graphsITS[i],grname[i]);
 775   }
 776   legend->Draw();
 777
 778   canvas->cd(2);
 779   legend = new TLegend(0.6,0.6,1.0,1.0,"#Delta sin(#phi) ITS-TPC");
 780   for (Int_t i=0; i<4; i++){
 781     graphsITS[i]->SetMaximum(1);
 782     graphsITS[i]->SetMinimum(-1);
 783     graphsITS[i]->SetMarkerColor(i+1);
 784     graphsITS[i]->SetMarkerStyle(25+i);
 785     graphsITS[i]->SetName(grname[i]);
 786     graphsITS[i]->GetXaxis()->SetTitle("#phi");
 787     graphsITS[i]->GetYaxis()->SetTitle("#Delta sin(#phi) (mrad)");
 788     if (i==0) graphsITS[i]->Draw("ap");
 789     graphsITS[i]->Draw("p");
 790     legend->AddEntry(graphsITS[i],grname[i]);
 791   }
 792   legend->Draw();
 793   //
 794   canvas->cd(3);
 795   legend = new TLegend(0.6,0.6,1.0,1.0,"#Delta sin(#phi) TPC-TRD");
 796   for (Int_t i=0; i<4; i++){
 797     graphsTRD[i]->SetMaximum(1);
 798     graphsTRD[i]->SetMinimum(-1);
 799     graphsTRD[i]->SetMarkerColor(i+1);
 800     graphsTRD[i]->SetMarkerStyle(25+i);
 801     graphsTRD[i]->SetName(grname[i]);
 802     graphsTRD[i]->GetXaxis()->SetTitle("#phi");
 803     graphsTRD[i]->GetYaxis()->SetTitle("#Deltasin(#phi) (mrad)");
 804     if (i==0) graphsTRD[i]->Draw("ap");
 805     graphsTRD[i]->Draw("p");
 806     legend->AddEntry(graphsITS[i],grname[i]);
 807   }
 808   legend->Draw();
 809   return canvas;
 810  }
 811
 812
 813 TCanvas * DrawFitdY(const char *name){
 814   //
 815   //
 816   //
 817   TH1::AddDirectory(0);
 818   TCanvas *canvas = new TCanvas(name,name,800,800);
 819   canvas->Divide(3,5);
 820   for (Int_t idet=0; idet<5; idet++){
 821     chain->SetMarkerColor(4);
 822     chain->SetMarkerStyle(25);
 823     chain->SetMarkerSize(0.3);
 824     chain->SetLineColor(2);
 825     //
 826     canvas->cd(idet*3+1);
 827     chain->Draw("mdelta:delta",cut+Form("ptype==0&&dtype==%d",idet),"");
 828     //
 829     canvas->cd(idet*3+2);
 830     chain->SetLineColor(2);
 831     chain->Draw("mdelta-delta",cut+Form("ptype==0&&dtype==%d",idet),"");
 832     chain->SetLineColor(1);
 833     chain->Draw("mdelta",cut+Form("ptype==0&&dtype==%d",idet),"same");
 834     //
 835     canvas->cd(idet*3+3);
 836     chain->SetMarkerColor(1);
 837     chain->Draw("mdelta:phi",cut+Form("ptype==0&&dtype==%d&&theta>0&&bz>0",idet),"");
 838     chain->SetMarkerColor(2);
 839     chain->Draw("mdelta-delta:phi",cut+Form("ptype==0&&dtype==%d&&theta>0&&bz>0",idet),"same");
 840
 841   }
 842   return canvas;
 843 }
 844
 845 TCanvas * DrawFitdSnp(const char *name){
 846   //
 847   //
 848   //
 849   TH1::AddDirectory(0);
 850   TCanvas *canvas = new TCanvas(name,name,800,800);
 851   canvas->Divide(3,5);
 852   {for (Int_t idet=0; idet<5; idet++){
 853     chain->SetMarkerColor(4);
 854     chain->SetMarkerStyle(25);
 855     chain->SetMarkerSize(0.3);
 856     chain->SetLineColor(2);
 857     //
 858     canvas->cd(idet*3+1);
 859     chain->Draw("1000*mdelta:1000*delta",cut+Form("ptype==2&&dtype==%d",idet),"");
 860     //
 861     canvas->cd(idet*3+2);
 862     chain->SetLineColor(2);
 863     chain->Draw("1000*(mdelta-delta)",cut+Form("ptype==2&&dtype==%d",idet),"");
 864     chain->SetLineColor(1);
 865     chain->Draw("1000*mdelta",cut+Form("ptype==2&&dtype==%d",idet),"same");
 866     //
 867     canvas->cd(idet*3+3);
 868     chain->SetMarkerColor(1);
 869     chain->Draw("1000*(mdelta):phi",cut+Form("ptype==2&&dtype==%d&&theta>0&&bz>0",idet),"");
 870     chain->SetMarkerColor(2);
 871     chain->Draw("1000*(mdelta-delta):phi",cut+Form("ptype==2&&dtype==%d&&theta>0&&bz>0",idet),"same");
 872     }}
 873   return canvas;
 874 }
 875
 876
 877
 878
 879
 880
 881 TCanvas * MakeComposedCorrection(const char *name){
 882
 883   TString  fit = chain->GetAlias("delta");
 884   TObjArray * array = fit.Tokenize("++");
 885   Int_t nfun=array->GetEntries();
 886   Double_t wt = 0.3 ;
 887   Double_t T1 = 1.0;
 888   Double_t T2 = 1.0;
 889   //
 890   // sign independent correction
 891   //
 892   AliTPCExBEffective            *eff      = new  AliTPCExBEffective;
 893   AliTPCCalibGlobalMisalignment *shiftITS = new  AliTPCCalibGlobalMisalignment;
 894   AliTPCExBTwist                *twist    = new  AliTPCExBTwist;
 895   //
 896   // sign dependent correction
 897   //
 898   AliTPCExBEffective            *effS      = new  AliTPCExBEffective;
 899   AliTPCCalibGlobalMisalignment *shiftITSS = new  AliTPCCalibGlobalMisalignment;
 900   AliTPCExBTwist                *twistS    = new  AliTPCExBTwist;
 901
 902   TMatrixD polA(100,4);
 903   TMatrixD polC(100,4);
 904   TMatrixD valA(100,1);
 905   TMatrixD valC(100,1);
 906   TMatrixD valAS(100,1);
 907   TMatrixD valCS(100,1);
 908   Int_t counterA=0;
 909   Int_t counterC=0;
 910   Int_t counterAS=0;
 911   Int_t counterCS=0;
 912
 913   {
 914     for (Int_t i=1; i<nfun;i++){
 915       TObjString fitName=array->At(i)->GetName();
 916       TObjString fitVal= &(fitName.String()[1+fitName.String().Last('(')]);
 917       Double_t value= fitVal.String().Atof();
 918       if (fitName.String().Contains("sign")) continue;
 919       if (fitName.String().Contains("tX1")){
 920         twist->SetXTwist(value*0.001);
 921       }
 922       if (fitName.String().Contains("tY1")){
 923         twist->SetYTwist(value*0.001);
 924       }
 925
 926       if (fitName.String().Contains("shiftX")&&fitName.String().Contains("isITS")){
 927         shiftITS->SetXShift(value*0.1);
 928       }
 929       if (fitName.String().Contains("shiftY")&&fitName.String().Contains("isITS")){
 930         shiftITS->SetYShift(value*0.1);
 931       }
 932
 933       if (fitName.String().Contains("eff")){
 934         Int_t index=fitName.String().First("_")-1;
 935         Int_t side=atoi(&(fitName.String()[index]));
 936         Int_t px  =atoi(&(fitName.String()[index+2]));
 937         Int_t pd  =atoi(&(fitName.String()[index+4]));
 938         Int_t pp  =0;
 939         //printf("%s\t%d\t%d\t%d\t%d\t%f\n",fitName.GetName(),side,px,pd,pp, value);
 940         if (side==0){
 941           polA(counterA,0)=0;
 942           polA(counterA,1)=px;
 943           polA(counterA,2)=pd;
 944           polA(counterA,3)=pp;
 945           valA(counterA,0)=value*0.1;
 946           counterA++;
 947         }
 948         if (side==1){
 949           polC(counterC,0)=0;
 950           polC(counterC,1)=px;
 951           polC(counterC,2)=pd;
 952           polC(counterC,3)=pp;
 953           valC(counterC,0)=value*0.1;
 954           counterC++;
 955         }
 956       }
 957     }
 958   }
 959   polA.ResizeTo(counterA,4);
 960   polC.ResizeTo(counterC,4);
 961   valA.ResizeTo(counterA,1);
 962   valC.ResizeTo(counterC,1);
 963   eff->SetPolynoms(&polA,&polC);
 964   eff->SetCoeficients(&valA,&valC);
 965   eff->SetOmegaTauT1T2(wt,T1,T2);
 966   shiftITS->SetOmegaTauT1T2(wt,T1,T2);
 967   twist->SetOmegaTauT1T2(wt,T1,T2);
 968   //
 969   //
 970   //
 971   {
 972     counterAS=0;
 973     counterCS=0;
 974     for (Int_t i=1; i<nfun;i++){
 975       TObjString fitName=array->At(i)->GetName();
 976       TObjString fitVal= &(fitName.String()[1+fitName.String().Last('(')]);
 977       if (!fitName.String().Contains("sign")) continue;
 978       Double_t value= fitVal.String().Atof();
 979       if (fitName.String().Contains("tX1")){
 980         twistS->SetXTwist(value*0.001);
 981       }
 982       if (fitName.String().Contains("tY1")){
 983         twistS->SetYTwist(value*0.001);
 984       }
 985
 986       if (fitName.String().Contains("shiftX")&&fitName.String().Contains("isITS")){
 987         shiftITSS->SetXShift(value*0.1);
 988       }
 989       if (fitName.String().Contains("shiftY")&&fitName.String().Contains("isITS")){
 990         shiftITSS->SetYShift(value*0.1);
 991       }
 992
 993       if (fitName.String().Contains("eff")){
 994         Int_t index=fitName.String().First("_")-1;
 995         Int_t side=atoi(&(fitName.String()[index]));
 996         Int_t px  =atoi(&(fitName.String()[index+2]));
 997         Int_t pd  =atoi(&(fitName.String()[index+4]));
 998         Int_t pp  =0;
 999         //printf("%s\t%d\t%d\t%d\t%d\t%f\n",fitName.GetName(),side,px,pd,pp, value);
1000         if (side==0){
1001           polA(counterAS,0)=0;
1002           polA(counterAS,1)=px;
1003           polA(counterAS,2)=pd;
1004           polA(counterAS,3)=pp;
1005           valAS(counterAS,0)=value*0.1;
1006           counterAS++;
1007         }
1008         if (side==1){
1009           polC(counterCS,0)=0;
1010           polC(counterCS,1)=px;
1011           polC(counterCS,2)=pd;
1012           polC(counterCS,3)=pp;
1013           valCS(counterCS,0)=value*0.1;
1014           counterCS++;
1015         }
1016       }
1017     }
1018   }
1019   polA.ResizeTo(counterA,4);
1020   polC.ResizeTo(counterC,4);
1021   valA.ResizeTo(counterA,1);
1022   valC.ResizeTo(counterC,1);
1023   effS->SetPolynoms(&polA,&polC);
1024   effS->SetCoeficients(&valAS,&valCS);
1025   effS->SetOmegaTauT1T2(wt,T1,T2);
1026   shiftITSS->SetOmegaTauT1T2(wt,T1,T2);
1027   twistS->SetOmegaTauT1T2(wt,T1,T2);
1028   //
1029   // Make combined correction
1030   //
1031
1032   TObjArray * corr0 = new TObjArray;
1033   TObjArray * corrP = new TObjArray;
1034   TObjArray * corrM = new TObjArray;
1035   AliTPCExBEffective            *eff0      = new  AliTPCExBEffective;
1036   AliTPCCalibGlobalMisalignment *shiftITSP = new  AliTPCCalibGlobalMisalignment;
1037   AliTPCExBTwist                *twistP    = new  AliTPCExBTwist;
1038   AliTPCExBEffective            *effP      = new  AliTPCExBEffective;
1039   AliTPCCalibGlobalMisalignment *shiftITSM = new  AliTPCCalibGlobalMisalignment;
1040   AliTPCExBTwist                *twistM    = new  AliTPCExBTwist;
1041   AliTPCExBEffective            *effM      = new  AliTPCExBEffective;
1042   //
1043   shiftITSP->SetXShift(shiftITS->GetXShift()+shiftITSS->GetXShift());    // shift due to the B field
1044   shiftITSP->SetYShift(shiftITS->GetXShift()+shiftITSS->GetYShift());
1045   shiftITSM->SetXShift(shiftITS->GetXShift()-shiftITSS->GetXShift());
1046   shiftITSM->SetYShift(shiftITS->GetXShift()-shiftITSS->GetYShift());
1047   //
1048   twistP->SetXTwist(twist->GetXTwist()+twistS->GetXTwist());     // twist between field  - both used
1049   twistP->SetYTwist(twist->GetYTwist()+twistS->GetYTwist());     //
1050   twistM->SetXTwist(twist->GetXTwist()-twistS->GetXTwist());
1051   twistM->SetYTwist(twist->GetYTwist()-twistS->GetYTwist());
1052   //
1053   effP->SetPolynoms(&polA,&polC);                                // effective correction
1054   effP->SetCoeficients(&valA,&valC);
1055   effM->SetPolynoms(&polA,&polC);
1056   effM->SetCoeficients(&valA,&valC);
1057   //
1058   eff0->SetPolynoms((TMatrixD*)&polA,(TMatrixD*)&polC);
1059   eff0->SetCoeficients((TMatrixD*)&valA,(TMatrixD*)&valC);
1060   //
1061   //
1062   //
1063
1064   //
1065   corrP->AddLast(shiftITSP);
1066   corrP->AddLast(twistP);
1067   corrP->AddLast(effP);
1068   corrM->AddLast(shiftITSM);
1069   corrM->AddLast(twistM);
1070   corrM->AddLast(effM);
1071   corr0->AddLast(eff0);
1072   //
1073
1074   AliTPCComposedCorrection *comp0= new AliTPCComposedCorrection ;
1075   AliTPCComposedCorrection *compP= new AliTPCComposedCorrection ;
1076   AliTPCComposedCorrection *compM= new AliTPCComposedCorrection ;
1077   comp0->SetCorrections((TObjArray*)(corr0));
1078   compP->SetCorrections((TObjArray*)(corrP));
1079   compM->SetCorrections((TObjArray*)(corrM));
1080   //
1081   comp0->SetOmegaTauT1T2(0*wt,T1,T2);
1082   compP->SetOmegaTauT1T2(wt,T1,T2);
1083   compM->SetOmegaTauT1T2(-wt,T1,T2);
1084
1085   TFile f("correctionsExB.root","update");
1086   f.mkdir(name);
1087   f.cd(name);
1088   printf("\nDump correction B=+0.5T:\n");
1089   compP->Print("da");
1090   printf("\nDump correction B=-0.5T:\n");
1091   compM->Print("da");
1092   //
1093   comp0->Write("ExB-Field0");
1094   compP->Write("ExB-Bplus");
1095   compM->Write("ExB-Bminus");
1096   //
1097   //
1098   TCanvas *c = new TCanvas(name,name,800,800);
1099   c->Divide(4,2);
1100   c->cd(1);
1101   shiftITS->CreateHistoDRPhiinXY()->Draw("surf2");
1102   c->cd(2);
1103   twist->CreateHistoDRPhiinXY()->Draw("surf2");
1104   c->cd(3);
1105   eff->CreateHistoDRPhiinZR()->Draw("surf2");
1106   c->cd(5);
1107   shiftITSS->CreateHistoDRPhiinXY()->Draw("surf2");
1108   c->cd(6);
1109   twistS->CreateHistoDRPhiinXY()->Draw("surf2");
1110   c->cd(7);
1111   effS->CreateHistoDRPhiinZR()->Draw("surf2");
1112   //
1113   c->cd(4);
1114   compP->CreateHistoDRPhiinZR()->Draw("surf2");
1115   c->cd(8);
1116   compM->CreateHistoDRPhiinZR()->Draw("surf2");
1117
1118
1119   return c;
1120 }
1121
1122
1123 void MakeOCDBEntry(Int_t refRun){
1124   //
1125   // make a Correction OCDB entry
1126   // take the fit values writen in config file
1127   //
1128   //
1129   // 1. Read previous value used in calibration
1130   //    OCDB has to be initialized before
1131
1132   gROOT->Macro(Form("ConfigCalibTrain.C(%d)",refRun));  // configuring calib db
1133   gROOT->LoadMacro("AddTaskTPCCalib.C");
1134   gROOT->ProcessLine(Form("ConfigOCDB(%d);",refRun));
1135   AliTPCCorrection *corr  = AliTPCcalibDB::Instance()->GetTPCComposedCorrection();
1136   //
1137   TFile f("correctionsExB.root");
1138   AliTPCComposedCorrection * corrP0=(AliTPCComposedCorrection *)f.Get("/Correction physical/ExB-Field0");
1139   AliTPCComposedCorrection * corrPP=(AliTPCComposedCorrection *)f.Get("/Correction physical/ExB-Bplus");
1140   AliTPCComposedCorrection * corrPM=(AliTPCComposedCorrection *)f.Get("/Correction physical/ExB-Bminus");
1141   //
1142   AliTPCComposedCorrection * corrPE0=(AliTPCComposedCorrection *)f.Get("/Correction physical+Effective/ExB-Field0");
1143   AliTPCComposedCorrection * corrPEP=(AliTPCComposedCorrection *)f.Get("/Correction physical+Effective/ExB-Bplus");
1144   AliTPCComposedCorrection * corrPEM=(AliTPCComposedCorrection *)f.Get("/Correction physical+Effective/ExB-Bminus");
1145   //
1146   corrP0->SetName("Field0 correction");
1147   corrPP->SetName("FieldP correction");
1148   corrPM->SetName("FieldM correction");
1149   corrPE0->SetName("Field0 correction");
1150   corrPEP->SetName("FieldP correction");
1151   corrPEM->SetName("FieldM correction");
1152
1153   TObjArray corrPhysical;
1154   TObjArray corrPhysicalEffective;
1155   //
1156   // add the base correction
1157   corrP0->GetCorrections()->Add(corr);
1158   corrPP->GetCorrections()->Add(corr);
1159   corrPM->GetCorrections()->Add(corr);
1160   corrPE0->GetCorrections()->Add(corr);
1161   corrPEP->GetCorrections()->Add(corr);
1162   corrPEM->GetCorrections()->Add(corr);
1163   //
1164   corrPhysical.AddLast(corrP0);  corrPhysical.AddLast(corrPP);  corrPhysical.AddLast(corrPM);
1165   corrPhysicalEffective.AddLast(corrPE0);  corrPhysicalEffective.AddLast(corrPEP);  corrPhysicalEffective.AddLast(corrPEM);
1166   //
1167   // make OCDB entries
1168   TString userName=gSystem->GetFromPipe("echo $USER");
1169   TString date=gSystem->GetFromPipe("date");
1170   TString ocdbStorage="local:////";
1171   ocdbStorage+=gSystem->GetFromPipe("pwd")+"/OCDB";
1172   //
1173   AliCDBMetaData *metaData= new AliCDBMetaData();
1174   metaData->SetObjectClassName("TObjArray");
1175   metaData->SetResponsible("Marian Ivanov");
1176   metaData->SetBeamPeriod(1);
1177   metaData->SetAliRootVersion("05-25-01"); //root version
1178   metaData->SetComment(Form("Correction calibration. User: %s\n Data%s",userName.Data(),date.Data()));
1179   AliCDBId* id1=NULL;
1180   id1=new AliCDBId("TPC/Calib/Correction", 0, AliCDBRunRange::Infinity());
1181   AliCDBStorage* gStorage = 0;
1182
1183   gStorage=AliCDBManager::Instance()->GetStorage((ocdbStorage+"Physical").Data());
1184   gStorage->Put(&corrPhysical, (*id1), metaData);
1185   gStorage = AliCDBManager::Instance()->GetStorage((ocdbStorage+"PhysicalEffective").Data());
1186   gStorage->Put(&corrPhysicalEffective, (*id1), metaData);
1187 }
1188