TPC/CalibMacros/MakeLookup.C

   1 //
   2 // Make lookup of distortion TPC distortion +(ITS and TRD)
   3 // Input: Residual histograms obtained in the AliTPCcalibTime
   4
   5 // Residual histograms:
   6 //   1. TPC-ITS  - entrance of the TPC
   7 //   2. TPC-ITS  - at the vertex
   8 //   3. TPC-TRD  - outer wall of the TPC
   9
  10 // Histogram binning:
  11 //   1. Theta    - fP3
  12 //   2. Phi      - global phi at the entrance (case 1,2) and at the outer wall of TPC (case 3)
  13 //   3. snp(phi) - fP2 - local inclination angle at reference X
  14
  15 // Output value:
  16 //   Mean residuals, rms and number of entries in each bing
  17 //
  18 //
  19
  20
  21 /*
  22   Example usage:
  23   gROOT->Macro("~/rootlogon.C")
  24   gSystem->AddIncludePath("-I$ALICE_ROOT/STAT")
  25   gSystem->AddIncludePath("-I$ALICE_ROOT/TPC")
  26   gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros")
  27   gSystem->Load("libANALYSIS");
  28   gSystem->Load("libTPCcalib");
  29   .L  $ALICE_ROOT/TPC/CalibMacros/MakeLookup.C+
  30   Int_t run=115888
  31   //     .x $ALICE_ROOT/ANALYSIS/CalibMacros/Pass0/ConfigCalibTrain.C(run,"local:///lustre/alice/alien/alice/data/2010/OCDB")
  32
  33   MakeLookup(run,0);
  34
  35   //Local check of procedure
  36   TTreeSRedirector * pcstream  = new TTreeSRedirector("mean.root");
  37   TFile f("CalibObjects.root");
  38   AliTPCcalibTime  *calibTime= f.Get("calibTime");
  39   THnSparse * his = calibTime->GetResHistoTPCITS(0);
  40   his->GetAxis(1)->SetRangeUser(-1,1);
  41   his->GetAxis(2)->SetRange(0,1000000);
  42   his->GetAxis(3)->SetRangeUser(-0.3,0.3);
  43   //
  44
  45
  46
  47 */
  48
  49 #if !defined(__CINT__) || defined(__MAKECINT__)
  50 #include "THnSparse.h"
  51 #include "TLatex.h"
  52 #include "TCanvas.h"
  53 #include "TLegend.h"
  54 #include "TSystem.h"
  55 #include "TFile.h"
  56 #include "TChain.h"
  57 #include "TCut.h"
  58 #include "TH3.h"
  59 #include "TProfile3D.h"
  60 #include "TMath.h"
  61 #include "TVectorD.h"
  62 #include "TMatrixD.h"
  63 #include "TStatToolkit.h"
  64 #include "TTreeStream.h"
  65 #include "AliExternalTrackParam.h"
  66 #include "AliESDfriend.h"
  67 #include "AliTPCcalibTime.h"
  68 #include "TROOT.h"
  69 #include "AliXRDPROOFtoolkit.h"
  70 #include "AliTPCCorrection.h"
  71 #include "AliTPCExBTwist.h"
  72 #include "AliTPCGGVoltError.h"
  73 #include "AliTPCComposedCorrection.h"
  74 #include "AliTPCExBConical.h"
  75 #include "TPostScript.h"
  76 #include "TStyle.h"
  77 #include "AliTrackerBase.h"
  78 #include "AliTPCCalibGlobalMisalignment.h"
  79 #include  "AliTPCExBEffective.h"
  80 #include  "AliTPCExBBShape.h"
  81 #endif
  82
  83 TChain * chain=0;
  84 void MakeLookup(Int_t run, Int_t mode);
  85 //void MakeLookupHisto(THnSparse * his, TTreeSRedirector *pcstream, const char* hname, Int_t run);
  86 void FitLookup(TChain *chainIn, const char *prefix, TVectorD &vecA, TVectorD &vecC, TVectorD& vecStatA, TVectorD &vecStatD, TCut cut, TObjArray *picArray);
  87 void MakeFits(Int_t run);
  88 void MakeFitTree();
  89 void DrawDistortionDy(TCut cutUser="", Double_t ymin=-0.6, Double_t ymax=0.6);
  90 void DrawDistortionMaps(const char *fname="mean.root");
  91 TCanvas *  DrawDistortionMaps(TTree * treedy, TTree * treedsnp, TTree* treed1Pt, const char * name, const char *title);
  92 AliTPCComposedCorrection * MakeComposedCorrection();
  93 void MakeLaserTree();
  94 void AddEffectiveCorrection(AliTPCComposedCorrection* comp);
  95
  96 void MakeLookup(Int_t run, Int_t mode){
  97   //
  98   // make a lookup tree with mean values
  99   // 5. make laser tree
 100   // 4.
 101   gSystem->AddIncludePath("-I$ALICE_ROOT/STAT");
 102   gSystem->Load("libANALYSIS");
 103   gSystem->Load("libTPCcalib");
 104   if (mode==5) {MakeLaserTree();return;};
 105   if (mode==4) {DrawDistortionMaps();return;}
 106   if (mode==1){
 107     DrawDistortionDy("abs(snp)<0.25");
 108     DrawDistortionMaps();
 109     return;
 110   }
 111   if (mode==2) return  MakeFitTree();
 112   TFile f("TPCTimeObjects.root");
 113   AliTPCcalibTime  *calibTime= (AliTPCcalibTime*)f.Get("calibTime");
 114   if (!calibTime){
 115      TFile*  f2 = new TFile("CalibObjects.root");
 116      calibTime= (AliTPCcalibTime*)f2->Get("calibTime");
 117      if (!calibTime){
 118        TFile*  f3 = new TFile("../CalibObjects.root");
 119        calibTime= (AliTPCcalibTime*)f3->Get("calibTime");
 120      }
 121   }
 122   //
 123   TTreeSRedirector * pcstream  = new TTreeSRedirector("mean.root");
 124   THnSparse * his = 0;
 125   const char * hname[5]={"dy","dz","dsnp","dtheta","d1pt"};
 126   if (calibTime){
 127     for (Int_t ihis=0; ihis<5;ihis++){
 128       his = calibTime->GetResHistoTPCITS(ihis);
 129       his->GetAxis(1)->SetRangeUser(-1.1,1.1);
 130       his->GetAxis(2)->SetRange(0,1000000);
 131       his->GetAxis(3)->SetRangeUser(-0.5,0.5);
 132       AliTPCCorrection::MakeDistortionMap(his,pcstream, Form("ITS%s",hname[ihis]),run);
 133       //
 134       his = calibTime->GetResHistoTPCvertex(ihis);
 135       his->GetAxis(1)->SetRangeUser(-1.1,1.1);
 136       his->GetAxis(2)->SetRange(0,1000000);
 137       his->GetAxis(3)->SetRangeUser(-0.5,0.5);
 138       AliTPCCorrection::MakeDistortionMap(his,pcstream, Form("Vertex%s",hname[ihis]),run);
 139       //
 140       his = calibTime->GetResHistoTPCTRD(ihis);
 141       his->GetAxis(1)->SetRangeUser(-1.1,1.1);
 142       his->GetAxis(2)->SetRange(0,1000000);
 143       his->GetAxis(3)->SetRangeUser(-0.5,0.5);
 144       AliTPCCorrection::MakeDistortionMap(his,pcstream, Form("TRD%s",hname[ihis]),run);
 145       his = calibTime->GetResHistoTPCTOF(ihis);
 146       if (his){
 147         his->GetAxis(1)->SetRangeUser(-1.1,1.1);
 148         his->GetAxis(2)->SetRange(0,1000000);
 149         his->GetAxis(3)->SetRangeUser(-0.5,0.5);
 150         AliTPCCorrection::MakeDistortionMap(his,pcstream, Form("TOF%s",hname[ihis]),run);
 151       }
 152     }
 153   }
 154   delete pcstream;
 155   gSystem->Exec("echo `pwd`/mean.root > distort.txt");
 156   MakeFits(run);
 157 }
 158
 159
 160 void MakeFits(Int_t run){
 161   //
 162   // Make the fits of distortion
 163   //   store fit results and QA pictures in the file distortFit.root
 164   //
 165   TCut cut="entries>50&&rms>0";
 166   TTreeSRedirector *pcstream = new TTreeSRedirector("distortFit.root");
 167   AliXRDPROOFtoolkit tool;
 168   TVectorD vecA[100],vecC[100], vecStatA[100],vecStatC[100];
 169   TObjArray * picArray= new TObjArray();
 170   //
 171   TChain *chainITSdy      = tool.MakeChain("distort.txt","ITSdy",0,100000);
 172   TChain *chainITSdsnp     = tool.MakeChain("distort.txt","ITSdsnp",0,100000);
 173   TChain *chainTRDdy      = tool.MakeChain("distort.txt","TRDdy",0,100000);
 174   TChain *chainTRDdsnp     = tool.MakeChain("distort.txt","TRDdsnp",0,100000);
 175   TChain *chainVertexdy   = tool.MakeChain("distort.txt","Vertexdy",0,100000);
 176   TChain *chainVertexdsnp  = tool.MakeChain("distort.txt","Vertexdsnp",0,100000);
 177
 178   (*pcstream)<<"fits"<<"run="<<run;
 179
 180   FitLookup(chainITSdy, "yits-tpc",vecA[0],vecC[0],vecStatA[0],vecStatC[0],cut, picArray);
 181   (*pcstream)<<"fits"<<"itsdyA.="<<&vecA[0]<<"itsdyC.="<<&vecC[0];
 182   (*pcstream)<<"fits"<<"itsdyAS.="<<&vecStatA[0]<<"itsdyCS.="<<&vecStatC[0];
 183   FitLookup(chainITSdsnp, "snpits-tpc",vecA[1],vecC[1],vecStatA[1],vecStatC[1],cut, picArray);
 184   (*pcstream)<<"fits"<<"itsdsnpA.="<<&vecA[1]<<"itsdsnpC.="<<&vecC[1];
 185   (*pcstream)<<"fits"<<"itsdsnpAS.="<<&vecStatA[1]<<"itsdsnpCS.="<<&vecStatC[1];
 186   //
 187   FitLookup(chainTRDdy, "ytrd-tpc",vecA[5],vecC[5],vecStatA[5],vecStatC[5],cut, picArray);
 188   (*pcstream)<<"fits"<<"trddyA.="<<&vecA[5]<<"trddyC="<<&vecC[5];
 189   (*pcstream)<<"fits"<<"trddyAS.="<<&vecStatA[5]<<"trddyCS.="<<&vecStatC[5];
 190   FitLookup(chainTRDdsnp, "snptrd-tpc",vecA[6],vecC[6],vecStatA[6],vecStatC[6],cut, picArray);
 191   (*pcstream)<<"fits"<<"trddsnpA.="<<&vecA[6]<<"trddsnpC.="<<&vecC[6];
 192   (*pcstream)<<"fits"<<"trddsnpAS.="<<&vecStatA[6]<<"trddsnpCS.="<<&vecStatC[6];
 193   //
 194   FitLookup(chainVertexdy, "yvertex-tpc",vecA[10],vecC[10],vecStatA[10],vecStatC[10],cut, picArray);
 195   (*pcstream)<<"fits"<<"vertexdyA.="<<&vecA[10]<<"vertexdyC.="<<&vecC[10];
 196   (*pcstream)<<"fits"<<"vertexdyAS.="<<&vecStatA[10]<<"vertexdyCS.="<<&vecStatC[10];
 197   FitLookup(chainVertexdsnp, "snpvertex-tpc",vecA[11],vecC[11],vecStatA[11],vecStatC[11],cut, picArray);
 198   (*pcstream)<<"fits"<<"vertexdsnpA.="<<&vecA[11]<<"vertexdsnpC.="<<&vecC[11];
 199   (*pcstream)<<"fits"<<"vertexdsnpAS.="<<&vecStatA[11]<<"vertexdsnpCS.="<<&vecStatC[11];
 200   //
 201   (*pcstream)<<"fits"<<"\n";
 202
 203   pcstream->GetFile()->cd();
 204   for (Int_t i=0;i<picArray->GetEntries(); i++){
 205     TObject * obj = picArray->At(i);
 206     if (obj) obj->Write(obj->GetName());
 207   }
 208   //
 209   //
 210   //
 211   chainITSdy->AddFriend(chainITSdy,"ITSdy");
 212   chainITSdy->AddFriend(chainITSdsnp,"ITSdsnp");
 213   chainITSdy->AddFriend(chainTRDdy,"TRDdy");
 214   chainITSdy->AddFriend(chainTRDdsnp,"TRDdsnp");
 215   chainITSdy->AddFriend(chainVertexdy,"Vertexdy");
 216   chainITSdy->AddFriend(chainVertexdsnp,"Vertexdsnp");
 217   TTree * tree = chainITSdy->CloneTree();
 218   tree->Write("distortionTree");
 219   delete pcstream;
 220
 221   //
 222   //
 223 }
 224
 225
 226
 227
 228
 229 void FitLookup(TChain *chainIn, const char *prefix, TVectorD &vecA, TVectorD &vecC, TVectorD& vecStatA, TVectorD &vecStatC,  TCut cut, TObjArray *picArray){
 230   //  TCut cut="entries>100&&rms>0";
 231   vecStatA.ResizeTo(6);
 232   vecStatC.ResizeTo(6);
 233   vecA.ResizeTo(10);
 234   vecC.ResizeTo(10);
 235   Int_t  npointsMax=30000000;
 236   TStatToolkit toolkit;
 237   Double_t chi2A=0;
 238   Double_t chi2C=0;
 239   Int_t    npointsA=0;
 240   Int_t    npointsC=0;
 241   TVectorD paramA;
 242   TVectorD paramC;
 243   TMatrixD covar;
 244   TString *strFitYA;
 245   TString *strFitYC;
 246
 247   TString  fstring="";   // magnetic part
 248   //fstring+="(1)++";             //0    - offset value
 249   fstring+="(cos(phi))++";        //1    - cos part
 250   fstring+="(sin(phi))++";        //2    - sin part
 251   //
 252   fstring+="(theta)++";           //3
 253   fstring+="(theta*cos(phi))++";  //4
 254   fstring+="(theta*sin(phi))++";  //5
 255   //
 256   fstring+="(snp)++";             //6   - delta X(radial)  coeficients - offset
 257   fstring+="(snp*cos(phi))++";    //7   - delta X(radial) coeficient
 258   fstring+="(snp*sin(phi))++";    //8
 259   //
 260   //
 261   strFitYA = TStatToolkit::FitPlane(chainIn,"mean", fstring.Data(),cut+"theta>0", chi2A,npointsA,paramA,covar,-1,0, npointsMax, kFALSE);
 262   strFitYC = TStatToolkit::FitPlane(chainIn,"mean", fstring.Data(),cut+"theta<0", chi2C,npointsC,paramC,covar,-1,0, npointsMax,kFALSE);
 263   strFitYA->Tokenize("++")->Print();
 264   strFitYC->Tokenize("++")->Print();
 265   chainIn->SetAlias("dyA",strFitYA->Data());
 266   chainIn->SetAlias("dyC",strFitYC->Data());
 267   //
 268   TH1 * his=0;
 269   vecA.ResizeTo(paramA.GetNrows());
 270   vecC.ResizeTo(paramC.GetNrows());
 271   vecA=paramA;
 272   vecC=paramC;
 273   // A side stat
 274   vecStatA[0]=npointsA;
 275   vecStatA[1]=TMath::Sqrt(chi2A/npointsA);
 276   chainIn->Draw("mean",cut+"theta>0");
 277   his=chainIn->GetHistogram();
 278   his->SetName(Form("Orig #Delta%sA",prefix));
 279   his->SetTitle(Form("Orig #Delta%sA",prefix));
 280   picArray->AddLast(his->Clone());
 281   vecStatA[2] = his->GetMean();
 282   vecStatA[3] = his->GetRMS();
 283   chainIn->Draw("mean-dyA",cut+"theta>0");
 284   his=chainIn->GetHistogram();
 285   his->SetName(Form("Corr. #Delta%sA",prefix));
 286   his->SetTitle(Form("Corr. #Delta%sA",prefix));
 287   picArray->AddLast(his->Clone());
 288   vecStatA[4] = his->GetMean();
 289   vecStatA[5] = his->GetRMS();
 290   //
 291   // C side stat
 292   vecStatC[0]=npointsC;
 293   vecStatC[1]=TMath::Sqrt(chi2C/npointsC);
 294   chainIn->Draw("mean",cut+"theta<0");
 295   his=chainIn->GetHistogram();
 296   his->SetName(Form("Orig #Delta%sC",prefix));
 297   his->SetTitle(Form("Orig #Delta%sC",prefix));
 298   picArray->AddLast(his->Clone());
 299   vecStatC[2] = his->GetMean();
 300   vecStatC[3] = his->GetRMS();
 301   chainIn->Draw("mean-dyC",cut+"theta<0");
 302   his=chainIn->GetHistogram();
 303   his->SetName(Form("Corr. #Delta%sC",prefix));
 304   his->SetTitle(Form("Corr. #Delta%sC",prefix));
 305   picArray->AddLast(his->Clone());
 306   vecStatC[4] = his->GetMean();
 307   vecStatC[5] = his->GetRMS();
 308
 309   vecStatA.Print();
 310   vecStatC.Print();
 311 }
 312
 313
 314
 315 void DrawDistortionDy(TCut cutUser, Double_t ymin, Double_t ymax){
 316   //
 317   //
 318   //
 319   TFile fplus("meanBplus.root");
 320   TFile fminus("meanBminus.root");
 321   TTree * titsDyPlus=  (TTree*)fplus.Get("ITSdy");
 322   TTree * titsDyMinus= (TTree*)fminus.Get("ITSdy");
 323   TTree * ttrdDyPlus=  (TTree*)fplus.Get("TRDdy");
 324   TTree * ttrdDyMinus= (TTree*)fminus.Get("TRDdy");
 325   //
 326   TTree * titsDsnpPlus=  (TTree*)fplus.Get("ITSdsnp");
 327   TTree * titsDsnpMinus= (TTree*)fminus.Get("ITSdsnp");
 328   TTree * ttrdDsnpPlus=  (TTree*)fplus.Get("TRDdsnp");
 329   TTree * ttrdDsnpMinus= (TTree*)fminus.Get("TRDdsnp");
 330   TTree * titsDthetaPlus=  (TTree*)fplus.Get("ITSdtheta");
 331   TTree * titsDthetaMinus= (TTree*)fminus.Get("ITSdtheta");
 332   TTree * ttrdDthetaPlus=  (TTree*)fplus.Get("TRDdtheta");
 333   TTree * ttrdDthetaMinus= (TTree*)fminus.Get("TRDdtheta");
 334   TCut cut="rms>0&&entries>100&&abs(theta)<1&&abs(snp)<0.3";
 335   cut+=cutUser;
 336   //
 337   titsDyPlus->AddFriend(titsDyMinus,"TM.");
 338   titsDyPlus->AddFriend(titsDsnpMinus,"Tsnp.");
 339   titsDyPlus->AddFriend(titsDthetaMinus,"Ttheta.");
 340   ttrdDyPlus->AddFriend(ttrdDyMinus,"TM.");
 341   titsDsnpPlus->AddFriend(titsDsnpMinus,"TM.");
 342   ttrdDsnpPlus->AddFriend(ttrdDsnpMinus,"TM.");
 343
 344   //
 345   //
 346   //
 347   Int_t entries=0;
 348   TGraph * graphITSYC[3];
 349   TGraph * graphITSYA[3];
 350   TGraph * graphTRDYC[3];
 351   TGraph * graphTRDYA[3];
 352   //
 353   //
 354   entries=titsDyPlus->Draw("mean-TM..mean:phi",cut+"theta<-0.1","goff");
 355   graphITSYC[0] = new TGraph(entries, titsDyPlus->GetV2(),titsDyPlus->GetV1());
 356   titsDyMinus->Draw("mean:phi",cut+"theta<-0.1","goff");
 357   graphITSYC[2] = new TGraph(entries, titsDyMinus->GetV2(),titsDyMinus->GetV1());
 358   titsDyPlus->Draw("mean:phi",cut+"theta<-0.1","goff");
 359   graphITSYC[1] = new TGraph(entries, titsDyPlus->GetV2(),titsDyPlus->GetV1());
 360   //
 361   entries=titsDyPlus->Draw("mean-TM..mean:phi",cut+"theta>0.1","goff");
 362   graphITSYA[0] = new TGraph(entries, titsDyPlus->GetV2(),titsDyPlus->GetV1());
 363   titsDyMinus->Draw("mean:phi",cut+"theta>0.1","goff");
 364   graphITSYA[2] = new TGraph(entries, titsDyMinus->GetV2(),titsDyMinus->GetV1());
 365   titsDyPlus->Draw("mean:phi",cut+"theta>0.1","goff");
 366   graphITSYA[1] = new TGraph(entries, titsDyPlus->GetV2(),titsDyPlus->GetV1());
 367   //
 368   entries=ttrdDyPlus->Draw("mean-TM..mean:phi",cut+"theta<-0.1&&TM..entries>100","goff");
 369   graphTRDYC[0] = new TGraph(entries, ttrdDyPlus->GetV2(),ttrdDyPlus->GetV1());
 370   ttrdDyMinus->Draw("mean:phi",cut+"theta<-0.1","goff");
 371   graphTRDYC[2] = new TGraph(entries, ttrdDyMinus->GetV2(),ttrdDyMinus->GetV1());
 372   ttrdDyPlus->Draw("mean:phi",cut+"theta<-0.1","goff");
 373   graphTRDYC[1] = new TGraph(entries, ttrdDyPlus->GetV2(),ttrdDyPlus->GetV1());
 374   //
 375   entries=ttrdDyPlus->Draw("mean-TM..mean:phi",cut+"theta>0.1&&TM..entries>100","goff");
 376   graphTRDYA[0] = new TGraph(entries, ttrdDyPlus->GetV2(),ttrdDyPlus->GetV1());
 377   ttrdDyMinus->Draw("mean:phi",cut+"theta>0.1","goff");
 378   graphTRDYA[2] = new TGraph(entries, ttrdDyMinus->GetV2(),ttrdDyMinus->GetV1());
 379   ttrdDyPlus->Draw("mean:phi",cut+"theta>0.1","goff");
 380   graphTRDYA[1] = new TGraph(entries, ttrdDyPlus->GetV2(),ttrdDyPlus->GetV1());
 381   //
 382
 383
 384   //
 385   {for (Int_t i=0; i<3; i++){
 386     graphITSYC[i]->SetMaximum(ymax+0.2); graphITSYC[i]->SetMinimum(ymin);
 387     graphITSYC[i]->GetXaxis()->SetTitle("#phi"); graphITSYC[i]->GetYaxis()->SetTitle("#Delta r#phi (cm)");
 388     graphITSYC[i]->SetMarkerColor(i+1); graphITSYC[i]->SetMarkerStyle(20+i);
 389     graphITSYA[i]->SetMaximum(ymax+0.2); graphITSYA[i]->SetMinimum(ymin);
 390     graphITSYA[i]->GetXaxis()->SetTitle("#phi"); graphITSYA[i]->GetYaxis()->SetTitle("#Delta r#phi (cm)");
 391     graphITSYA[i]->SetMarkerColor(i+1); graphITSYA[i]->SetMarkerStyle(20+i);
 392     graphTRDYC[i]->SetMaximum(ymax+0.2); graphTRDYC[i]->SetMinimum(ymin);
 393     graphTRDYC[i]->GetXaxis()->SetTitle("#phi"); graphTRDYC[i]->GetYaxis()->SetTitle("#Delta r#phi (cm)");
 394     graphTRDYC[i]->SetMarkerColor(i+1); graphTRDYC[i]->SetMarkerStyle(20+i);
 395     graphTRDYA[i]->SetMaximum(ymax+0.2); graphTRDYA[i]->SetMinimum(ymin);
 396     graphTRDYA[i]->GetXaxis()->SetTitle("#phi"); graphTRDYA[i]->GetYaxis()->SetTitle("#Delta r#phi (cm)");
 397     graphTRDYA[i]->SetMarkerColor(i+1); graphTRDYA[i]->SetMarkerStyle(20+i);
 398     }
 399   }
 400   TLegend *legend=0;
 401   TString cname="cdeltaY";
 402   TString dirName=gSystem->GetFromPipe("dirname `pwd` | xargs basename ");
 403   TString splus=gSystem->GetFromPipe("ls -al meanBplus.root | gawk '{print \"File: \"$8\" \"$6\" \"$7\" \";}'");
 404   TString sminus=gSystem->GetFromPipe("ls -al meanBminus.root | gawk '{print \"File: \"$8\" \"$6\" \"$7\" \";}'");
 405
 406   TCanvas *cdeltaY = new TCanvas("cdeltaY","cdeltaY",1300,800);
 407   cdeltaY->Divide(2,2);
 408
 409   cdeltaY->cd(1);
 410   graphITSYC[0]->Draw("ap");
 411   graphITSYC[1]->Draw("p");
 412   graphITSYC[2]->Draw("p");
 413   legend = new TLegend(0.6,0.6,1.0,1.0,"ITS-TPC C side #Delta r-#phi");
 414   legend->AddEntry("",dirName.Data());
 415   legend->AddEntry("",splus.Data());
 416   legend->AddEntry("",sminus.Data());
 417   legend->AddEntry(graphITSYC[0],"B_{plus}-B_{minus}");
 418   legend->AddEntry(graphITSYC[1],"B_{plus}");
 419   legend->AddEntry(graphITSYC[2],"B_{minus}");
 420   legend->Draw();
 421
 422   cdeltaY->cd(2);
 423   graphITSYA[0]->Draw("ap");
 424   graphITSYA[1]->Draw("p");
 425   graphITSYA[2]->Draw("p");
 426   legend = new TLegend(0.6,0.6,1.0,1.0,"ITS-TPC A side #Delta r-#phi");
 427   legend->AddEntry(graphITSYA[0],"B_{plus}-B_{minus}");
 428   legend->AddEntry(graphITSYA[1],"B_{plus}");
 429   legend->AddEntry(graphITSYA[2],"B_{minus}");
 430   legend->Draw();
 431
 432   //
 433   cdeltaY->cd(3);
 434   graphTRDYC[0]->Draw("ap");
 435   graphTRDYC[1]->Draw("p");
 436   graphTRDYC[2]->Draw("p");
 437   legend = new TLegend(0.6,0.6,1.0,1.0,"TRD-TPC C side #Delta r-#phi");
 438   legend->AddEntry(graphTRDYC[0],"B_{plus}-B_{minus}");
 439   legend->AddEntry(graphTRDYC[1],"B_{plus}");
 440   legend->AddEntry(graphTRDYC[2],"B_{minus}");
 441   legend->Draw();
 442
 443   cdeltaY->cd(4);
 444   graphTRDYA[0]->Draw("ap");
 445   graphTRDYA[1]->Draw("p");
 446   graphTRDYA[2]->Draw("p");
 447   legend = new TLegend(0.6,0.6,1.0,1.0,"TRD-TPC A side #Delta r-#phi");
 448   legend->AddEntry(graphTRDYA[0],"B_{plus}-B_{minus}");
 449   legend->AddEntry(graphTRDYA[1],"B_{plus}");
 450   legend->AddEntry(graphTRDYA[2],"B_{minus}");
 451   legend->Draw();
 452   cdeltaY->SaveAs(dirName+"_deltaY.pdf");
 453
 454 }
 455
 456
 457
 458
 459
 460
 461 void MakeFitTree(){
 462   //
 463   // 1. Initialize ocdb e.g
 464   //     .x $ALICE_ROOT/ANALYSIS/CalibMacros/Pass0/ConfigCalibTrain.C(114972,)
 465   //     .x $ALICE_ROOT/ANALYSIS/CalibMacros/Pass0/ConfigCalibTrain.C(114972,"local:///lustre/alice/alien/alice/data/2010/OCDB")
 466   //
 467   //
 468
 469   AliTPCComposedCorrection *cc= MakeComposedCorrection();
 470   TObjArray * corr = (TObjArray*)(cc->GetCorrections());
 471   //  corr->AddLast(cc);
 472   const Int_t kskip=23;
 473   //
 474   TFile f("mean.root");
 475   TTree * tree= 0;
 476   tree = (TTree*)f.Get("ITSdy");
 477   printf("ITSdy\n");
 478   AliTPCCorrection::MakeTrackDistortionTree(tree,0,0,corr,kskip);
 479   tree = (TTree*)f.Get("TRDdy");
 480   printf("TRDdy\n");
 481   AliTPCCorrection::MakeTrackDistortionTree(tree,1,0,corr,kskip);
 482   tree = (TTree*)f.Get("Vertexdy");
 483   printf("Vertexdy\n");
 484   AliTPCCorrection::MakeTrackDistortionTree(tree,2,0,corr,kskip);
 485   tree = (TTree*)f.Get("TOFdy");
 486   printf("TOFdy\n");
 487   if (tree) AliTPCCorrection::MakeTrackDistortionTree(tree,3,0,corr,kskip);
 488
 489   //
 490   tree = (TTree*)f.Get("ITSdsnp");
 491   printf("ITSdsnp\n");
 492   AliTPCCorrection::MakeTrackDistortionTree(tree,0,2,corr,kskip);
 493   tree = (TTree*)f.Get("TRDdsnp");
 494   printf("TRDdsnp\n");
 495   AliTPCCorrection::MakeTrackDistortionTree(tree,1,2,corr,kskip);
 496   tree = (TTree*)f.Get("Vertexdsnp");
 497   printf("Vertexdsnp\n");
 498   AliTPCCorrection::MakeTrackDistortionTree(tree,2,2,corr,kskip);
 499   //
 500
 501   tree = (TTree*)f.Get("ITSd1pt");
 502   printf("ITSd1pt\n");
 503   if (tree) AliTPCCorrection::MakeTrackDistortionTree(tree,0,4,corr,kskip);
 504
 505   tree = (TTree*)f.Get("TOFdz");
 506   printf("TOFdz\n");
 507   if (tree) AliTPCCorrection::MakeTrackDistortionTree(tree,3,0,corr,kskip);
 508
 509
 510   tree = (TTree*)f.Get("ITSdz");
 511   printf("ITSdz\n");
 512   AliTPCCorrection::MakeTrackDistortionTree(tree,0,1,corr,kskip);
 513   tree = (TTree*)f.Get("Vertexdz");
 514   printf("Vertexdz\n");
 515   AliTPCCorrection::MakeTrackDistortionTree(tree,2,1,corr,kskip);
 516
 517   //
 518   tree = (TTree*)f.Get("ITSdtheta");
 519   printf("ITSdtheta\n");
 520   AliTPCCorrection::MakeTrackDistortionTree(tree,0,3,corr,kskip);
 521   //
 522   tree = (TTree*)f.Get("Vertexdtheta");
 523   printf("Vertexdtheta\n");
 524   AliTPCCorrection::MakeTrackDistortionTree(tree,2,3,corr,kskip);
 525
 526
 527 }
 528
 529
 530
 531 void MakeGlobalFit(){
 532   AliXRDPROOFtoolkit tool;
 533   TChain *chain      = tool.MakeChain("distortion.txt","fit",0,100000);
 534
 535   TCut cutS="rms>0&&entries>100";
 536   TCut cut=cutS+"(ptype==0||ptype==2||ptype==3||(ptype==4&&dtype==0))";
 537   Int_t  npointsMax=30000000;
 538   TStatToolkit toolkit;
 539   Double_t chi2=0;
 540   Int_t    npoints=0;
 541   TVectorD param;
 542   TMatrixD covar;
 543   chain->SetAlias("side","(1+(theta>0)*2)");
 544   TString  fstring="";     // magnetic part
 545   fstring+="tX1++";        //1    - twist x in mrad
 546   fstring+="tY1++";        //2    - twist y in mrad
 547   fstring+="ExBCon++";     //3    - custom
 548   fstring+="ExBCon*cos(phi)++";     //3    - custom
 549   fstring+="ExBCon*sin(phi)++";     //3    - custom
 550   fstring+="ggA++";        //4    - gating grid
 551   fstring+="ggA*cos(phi)++";        //4    - gating grid
 552   fstring+="ggA*sin(phi)++";        //4    - gating grid
 553   fstring+="ggC++";        //5    - gating grid
 554   fstring+="ggC*cos(phi)++";        //5    - gating grid
 555   fstring+="ggC*sin(phi)++";        //5    - gating grid
 556   fstring+="(ptype==2)++";
 557   fstring+="(ptype==3)++";
 558   fstring+="(ptype==4)++";
 559   fstring+="(ptype==4)*bz++";
 560
 561   TString *strDelta = TStatToolkit::FitPlane(chain,"mean:rms", fstring.Data(),cut+"(Entry$%13==0)", chi2,npoints,param,covar,-1,0, npointsMax, kTRUE);
 562   chain->SetAlias("delta",strDelta->Data());
 563   strDelta->Tokenize("++")->Print();
 564
 565 }
 566
 567
 568 void MakeGlobalFitRelative(Int_t highFrequency=0){
 569   //
 570   // Make a global fit of ExB
 571   // To get rid of the misalignment errors -
 572   //   Use relative change of deltas for 2 different filed settings
 573   //
 574
 575   AliXRDPROOFtoolkit tool;
 576   //  TChain *chain         = tool.MakeChain("distortion.txt","fit",0,100000);
 577   TChain *chainPlus     = tool.MakeChain("distortionPlus.txt","fit",0,100000);
 578   TChain *chainMinus    = tool.MakeChain("distortionMinus.txt","fit",0,100000);
 579   TChain *chain0        = tool.MakeChain("distortionField0.txt","fit",0,100000);
 580
 581   chainPlus->AddFriend(chainMinus,"M");
 582   chainPlus->AddFriend(chain0,"F0");
 583
 584   TCut cut="rms>0&&M.rms>0&&entries>100&&dtype==0&&(ptype==0||ptype==2||ptype==4)";
 585   chainPlus->SetAlias("deltaM","mean-M.mean");
 586   Int_t  npointsMax=30000000;
 587   TStatToolkit toolkit;
 588   Double_t chi2=0;
 589   Int_t    npoints=0;
 590   TVectorD param;
 591   TMatrixD covar;
 592   //
 593   TString  fstring="";     // fit string
 594   fstring+="(ptype==0)*cos(phi)++";   //     - primitive gx movement
 595   fstring+="(ptype==0)*sin(phi)++";   //     - primitive gy movement
 596
 597   fstring+="(tX1-M.tX1)++";           // 1    - twist x in mrad
 598   fstring+="(tY1-M.tY1)++";           // 2    - twist y in mrad
 599   //
 600   fstring+="(ExBConA-M.ExBConA)++";   // 3 conical shape A side
 601   fstring+="(ExBConC-M.ExBConC)++";   // 4 conical shape C side
 602   //
 603   fstring+="(ggA-M.ggA)++";   // 3 conical shape A side
 604   fstring+="(ggC-M.ggC)++";   // 4 conical shape C side
 605
 606   if (highFrequency) {for (Int_t i=1; i<highFrequency; i++){
 607       fstring+=Form("((ggA-M.ggA)*cos(%d*phi))++",i);
 608       fstring+=Form("((ggA-M.ggA)*sin(%d*phi))++",i);
 609       //
 610       fstring+=Form("((ggC-M.ggC)*cos(%d*phi))++",i);
 611       fstring+=Form("((ggC-M.ggC)*sin(%d*phi))++",i);
 612       if (i>1){
 613         fstring+=Form("((ExBConA-M.ExBConA)*cos(%d*phi))++",i);
 614         fstring+=Form("((ExBConA-M.ExBConA)*sin(%d*phi))++",i);
 615         fstring+=Form("((ExBConC-M.ExBConC)*cos(%d*phi))++",i);
 616         fstring+=Form("((ExBConC-M.ExBConC)*sin(%d*phi))++",i);
 617       }
 618     }}
 619
 620   TString *strDelta = TStatToolkit::FitPlane(chainPlus,"deltaM:rms", fstring.Data(),cut+"(Entry$%7==0)", chi2,npoints,param,covar,-1,0, npointsMax, kTRUE);
 621   chainPlus->SetAlias("delta",strDelta->Data());
 622   TObjArray *fitArray = strDelta->Tokenize("++");
 623   fitArray->Print();
 624
 625   //
 626   // Draw results
 627   //
 628   TH1::AddDirectory(0);
 629   TLatex *lfit=new TLatex;
 630
 631
 632   TCanvas * canvasdY= new TCanvas(Form("deltaY%d",highFrequency),Form("deltaY%d",highFrequency),1200,800);
 633   TCanvas * canvasdSnp= new TCanvas(Form("deltaSnp%d",highFrequency),Form("deltaSnp%d",highFrequency),1200,800);
 634   TCanvas * canvasd1pt= new TCanvas(Form("delta1pt%d",highFrequency),Form("delta1pt%d",highFrequency),1200,800);
 635
 636   canvasdY->Divide(3,2);
 637   chainPlus->SetMarkerStyle(25);
 638   chainPlus->SetMarkerSize(0.5);
 639   chainPlus->SetMarkerColor(1);
 640   canvasdY->cd(1);
 641   chainPlus->Draw("deltaM:delta",cut+"theta>0&&abs(snp)<0.2&&ptype==0","");
 642   canvasdY->cd(2);
 643   chainPlus->SetMarkerColor(1);
 644   chainPlus->Draw("deltaM-delta>>deltaCorr(50,-0.4,0.4)",cut+"theta>0&&abs(snp)<0.2&&ptype==0","err");
 645   chainPlus->SetMarkerColor(3);
 646   chainPlus->Draw("deltaM",cut+"theta>0&&abs(snp)<0.2&&ptype==0","errsame");
 647   canvasdY->cd(3);
 648   chainPlus->SetMarkerColor(1);
 649   chainPlus->Draw("deltaM:phi",cut+"theta>0&&abs(snp)<0.2&&ptype==0","");
 650   chainPlus->SetMarkerColor(2);
 651   chainPlus->Draw("deltaM-delta:phi",cut+"theta>0&&abs(snp)<0.2&&ptype==0","same");
 652   canvasdY->cd(6);
 653   chainPlus->SetMarkerStyle(26);
 654   chainPlus->SetMarkerColor(1);
 655   chainPlus->Draw("deltaM:phi",cut+"theta<0&&abs(snp)<0.2&&ptype==0","");
 656   chainPlus->SetMarkerColor(4);
 657   chainPlus->Draw("deltaM-delta:phi",cut+"theta<0&&abs(snp)<0.2&&ptype==0","same");
 658
 659
 660
 661   canvasdSnp->Divide(3,2);
 662   chainPlus->SetMarkerStyle(25);
 663   chainPlus->SetMarkerSize(0.5);
 664   chainPlus->SetMarkerColor(1);
 665   canvasdSnp->cd(1);
 666   chainPlus->Draw("1000*deltaM:1000*delta",cut+"theta>0&&abs(snp)<0.2&&ptype==2","");
 667   canvasdSnp->cd(2);
 668   chainPlus->SetMarkerColor(1);
 669   chainPlus->Draw("1000*(deltaM-delta)>>deltaCorr(50,-4.,4)",cut+"theta>0&&abs(snp)<0.2&&ptype==2","err");
 670   chainPlus->SetMarkerColor(3);
 671   chainPlus->Draw("1000*deltaM",cut+"theta>0&&abs(snp)<0.2&&ptype==2","errsame");
 672   canvasdSnp->cd(3);
 673   chainPlus->SetMarkerColor(1);
 674   chainPlus->Draw("1000*(deltaM):phi",cut+"theta>0&&abs(snp)<0.2&&ptype==2","");
 675   chainPlus->SetMarkerColor(2);
 676   chainPlus->Draw("1000*(deltaM-delta):phi",cut+"theta>0&&abs(snp)<0.2&&ptype==2","same");
 677   canvasdSnp->cd(6);
 678   chainPlus->SetMarkerStyle(26);
 679   chainPlus->SetMarkerColor(1);
 680   chainPlus->Draw("1000*(deltaM):phi",cut+"theta<0&&abs(snp)<0.2&&ptype==2","");
 681   chainPlus->SetMarkerColor(4);
 682   chainPlus->Draw("1000*(deltaM-delta):phi",cut+"theta<0&&abs(snp)<0.2&&ptype==2","same");
 683
 684
 685
 686   canvasd1pt->Divide(3,2);
 687   chainPlus->SetMarkerStyle(25);
 688   chainPlus->SetMarkerSize(0.5);
 689   chainPlus->SetMarkerColor(1);
 690   canvasd1pt->cd(1);
 691   chainPlus->Draw("deltaM:delta",cut+"theta>0&&abs(snp)<0.2&&ptype==4","");
 692   canvasd1pt->cd(2);
 693   chainPlus->SetMarkerColor(1);
 694   chainPlus->Draw("deltaM-delta>>deltaCorr(50,-0.15,0.15)",cut+"theta>0&&abs(snp)<0.2&&ptype==4","err");
 695   chainPlus->GetHistogram()->Fit("gaus");
 696   chainPlus->SetMarkerColor(3);
 697   chainPlus->Draw("deltaM",cut+"theta>0&&abs(snp)<0.2&&ptype==4","errsame");
 698   canvasd1pt->cd(3);
 699   chainPlus->SetMarkerColor(1);
 700   chainPlus->Draw("deltaM:phi",cut+"theta>0&&abs(snp)<0.2&&ptype==4","");
 701   chainPlus->SetMarkerColor(2);
 702   chainPlus->Draw("deltaM-delta:phi",cut+"theta>0&&abs(snp)<0.2&&ptype==4","same");
 703   canvasd1pt->cd(6);
 704   chainPlus->SetMarkerStyle(26);
 705   chainPlus->SetMarkerColor(1);
 706   chainPlus->Draw("deltaM:phi",cut+"theta<0&&abs(snp)<0.2&&ptype==4","");
 707   chainPlus->SetMarkerColor(4);
 708   chainPlus->Draw("deltaM-delta:phi",cut+"theta<0&&abs(snp)<0.2&&ptype==4","same");
 709
 710   {for (Int_t ipad=0; ipad<3;ipad++){
 711     if (ipad==0) canvasdY->cd(5);
 712     if (ipad==1) canvasdSnp->cd(5);
 713     if (ipad==2) canvasd1pt->cd(5);
 714     lfit->SetTextAlign(12); lfit->SetTextSize(0.04);
 715     {for (Int_t i=1; i<fitArray->GetEntries(); i++){
 716         lfit->DrawLatex(0.1,1-i/9.,fitArray->At(i)->GetName());
 717       }}
 718     if (ipad==0) canvasdY->cd(4);
 719     if (ipad==1) canvasdSnp->cd(4);
 720     if (ipad==2) canvasd1pt->cd(4);
 721     lfit->DrawLatex(0.1,0.9,"Global fit - TPC ITS matching in r-#phi");
 722     lfit->DrawLatex(0.1,0.8,"Residual minimization:");
 723     lfit->DrawLatex(0.1,0.7,"r#phi_{0.5T}-r#phi_{-0.5T} (cm)");
 724     lfit->DrawLatex(0.1,0.6,"sin(r#phi)_{0.5T}-sin(r#phi)_{-0.5T} (mrad)");
 725     lfit->DrawLatex(0.1,0.5,"1/pt_{0.5T}-1/pt_{-0.5T} (1/GeV)");
 726     }}
 727
 728
 729   TFile *fpic = new TFile(Form("fitPictures%d.root",highFrequency),"update");
 730   canvasd1pt->Write();
 731   canvasdY->Write();
 732   canvasdSnp->Write();
 733   //
 734   canvasd1pt->SaveAs(Form("fitd1pt%d.pdf",highFrequency));
 735   canvasdY->SaveAs(Form("fitdy%d.pdf",highFrequency));
 736   canvasdSnp->SaveAs(Form("fitdsnp%d.pdf",highFrequency));
 737   delete fpic;
 738 }
 739
 740 void DrawDistortionMaps(const char *fname){
 741   TFile f(fname);
 742   TTree *ITSdy=(TTree*)f.Get("ITSdy");
 743   TTree *ITSdsnp=(TTree*)f.Get("ITSdsnp");
 744   TTree *ITSd1pt=(TTree*)f.Get("ITSd1pt");
 745   TTree *TRDdy=(TTree*)f.Get("TRDdy");
 746   TTree *TRDdsnp=(TTree*)f.Get("TRDdsnp");
 747   TTree *TRDd1pt=(TTree*)f.Get("TRDd1pt");
 748   TTree *Vertexdy=(TTree*)f.Get("Vertexdy");
 749   TTree *Vertexdsnp=(TTree*)f.Get("Vertexdsnp");
 750   TTree *Vertexd1pt=(TTree*)f.Get("Vertexd1pt");
 751   TCanvas *cits = DrawDistortionMaps(ITSdy,ITSdsnp,ITSd1pt,"ITS","ITS");
 752   cits->SaveAs("matchingITS.pdf");
 753   TCanvas *ctrd = DrawDistortionMaps(TRDdy,TRDdsnp,TRDd1pt,"TRD","TRD");
 754   ctrd->SaveAs("matchingTRD.pdf");
 755   TCanvas *cvertex = DrawDistortionMaps(Vertexdy,Vertexdsnp,Vertexd1pt,"Vertex","Vertex");
 756   cvertex->SaveAs("matchingVertex.pdf");
 757   //
 758   //
 759   TPostScript *ps = new TPostScript("matching.ps", 112);
 760   gStyle->SetOptTitle(1);
 761   gStyle->SetOptStat(0);
 762   ps->NewPage();
 763   cits->Draw();
 764   ps->NewPage();
 765   ctrd->Draw();
 766   ps->NewPage();
 767   cvertex->Draw();
 768   ps->Close();
 769   delete ps;
 770 }
 771
 772
 773 TCanvas * DrawDistortionMaps(TTree * treedy, TTree * treedsnp, TTree* treed1pt, const char * name, const char *title){
 774   //
 775   //
 776   //
 777   TH1::AddDirectory(0);
 778   TCanvas *cdist = new TCanvas(name,name,1200,800);
 779   cdist->Divide(3,2);
 780   cdist->Draw("");
 781   TH1 * his=0;
 782
 783   cdist->cd(1);
 784   treedy->Draw("10*mean","rms>0&&abs(snp)<0.25&&entries>100","");
 785   his = (TH1*)treedy->GetHistogram()->Clone();
 786   his->SetName("dY"); his->SetXTitle("#Delta_{r#phi} (cm)");
 787   his->Draw();
 788   //
 789   cdist->cd(2);
 790   treedsnp->Draw("1000*mean","rms>0&&abs(snp)<0.25&&entries>200","");
 791   his = (TH1*)treedsnp->GetHistogram()->Clone();
 792   his->SetName("dsnp"); his->SetXTitle("#Delta_{sin(#phi)}");
 793   his->Draw();
 794   //
 795   cdist->cd(3);
 796   treed1pt->Draw("mean","rms>0&&abs(snp)<0.15&&entries>400","");
 797   his = (TH1*)treed1pt->GetHistogram()->Clone();
 798   his->SetName("d1pt"); his->SetXTitle("#Delta_{1/pt} (1/GeV)");
 799   his->Draw();
 800   //
 801   treedy->SetMarkerSize(1);
 802   treedsnp->SetMarkerSize(1);
 803   treed1pt->SetMarkerSize(1);
 804   {for (Int_t type=0; type<3; type++){
 805       cdist->cd(4+type);
 806       TTree * tree =treedy;
 807       if (type==1) tree=treedsnp;
 808       if (type==2) tree=treed1pt;
 809       Int_t counter=0;
 810       for (Double_t theta=-1; theta<=1; theta+=0.2){
 811         if (TMath::Abs(theta)<0.11) continue;
 812         TCut cut=Form("rms>0&&abs(snp)<0.25&&entries>20&&abs(theta-%f)<0.1",theta);
 813         tree->SetMarkerStyle(20+counter);
 814         Option_t *option=(counter==0)? "": "same";
 815         if (theta>0) tree->SetMarkerColor(2);
 816         if (theta<0) tree->SetMarkerColor(4);
 817         if (type==0) tree->Draw("10*mean:phi>>his(45,-pi,pi)",cut,"profgoff");
 818         if (type==1) tree->Draw("1000*mean:phi>>his(45,-pi,pi)",cut,"profgoff");
 819         if (type==2) tree->Draw("mean:phi>>his(45,-pi,pi)",cut,"profgoff");
 820         his = (TH1*)tree->GetHistogram()->Clone();
 821         his->SetName(Form("%d_%d",counter,type));
 822         his->SetXTitle("#phi");
 823         his->SetMaximum(4);
 824         his->SetMinimum(-4);
 825         if (type==2){
 826           his->SetMaximum(0.06);
 827           his->SetMinimum(-0.06);
 828         }
 829         his->Draw(option);
 830         counter++;
 831       }
 832     }
 833   }
 834   //  cdist->SaveAs(Form("%s.pdf"),name);
 835   return cdist;
 836 }
 837
 838
 839
 840
 841
 842 AliTPCComposedCorrection * MakeComposedCorrection(){
 843
 844   //
 845   Double_t bzField=AliTrackerBase::GetBz();
 846   AliMagF* magF= (AliMagF*)(TGeoGlobalMagField::Instance()->GetField());
 847   Double_t vdrift = 2.6; // [cm/us]   // to be updated: per second (ideally)
 848   Double_t ezField = 400; // [V/cm]   // to be updated: never (hopefully)
 849   Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ;
 850   Double_t T1 = 1.0;
 851   Double_t T2 = 1.0;
 852   //
 853   //
 854
 855
 856   AliTPCExBTwist *twistX1= new AliTPCExBTwist;
 857   twistX1->SetXTwist(0.001);  // 1 mrad twist in x
 858   twistX1->SetName("tX1");
 859   //
 860   AliTPCExBTwist *twistY1= new AliTPCExBTwist;
 861   twistY1->SetYTwist(0.001);   // 1 mrad twist in Y
 862   twistY1->SetName("tY1");
 863   //
 864   AliTPCGGVoltError* ggErrorA = new AliTPCGGVoltError;
 865   ggErrorA->SetDeltaVGGA(40.);  // delta GG - 1 mm
 866   ggErrorA->SetName("ggA");
 867   //
 868   AliTPCGGVoltError* ggErrorC = new AliTPCGGVoltError;
 869   ggErrorC->SetDeltaVGGC(40.);  // delta GG - 1 mm
 870   ggErrorC->SetName("ggC");
 871   // conical free param
 872   //
 873   AliTPCCalibGlobalMisalignment *shiftX=new  AliTPCCalibGlobalMisalignment;
 874   shiftX->SetXShift(0.1);
 875   shiftX->SetName("shiftX");
 876   AliTPCCalibGlobalMisalignment *shiftY=new  AliTPCCalibGlobalMisalignment;
 877   shiftY->SetYShift(0.1);
 878   shiftY->SetName("shiftY");
 879
 880   AliTPCExBBShape *exb11 = new AliTPCExBBShape;
 881   exb11->SetBField(magF);
 882   exb11->SetName("exb11");
 883   AliTPCExBBShape *exbT1 = new AliTPCExBBShape;
 884   exbT1->SetBField(magF);
 885   exbT1->SetName("exbT1");
 886   AliTPCExBBShape *exbT2 = new AliTPCExBBShape;
 887   exbT2->SetBField(magF);
 888   exbT2->SetName("exbT2");
 889
 890   TObjArray * corr = new TObjArray;
 891   corr->AddLast(twistX1);
 892   corr->AddLast(twistY1);
 893   corr->AddLast(ggErrorA);
 894   corr->AddLast(ggErrorC);
 895   corr->AddLast(shiftX);
 896   corr->AddLast(shiftY);
 897   corr->AddLast(exb11);
 898   corr->AddLast(exbT1);
 899   corr->AddLast(exbT2);
 900   //
 901   AliTPCComposedCorrection *cc= new AliTPCComposedCorrection ;
 902   cc->SetCorrections((TObjArray*)(corr));
 903   AddEffectiveCorrection(cc);
 904   cc->SetOmegaTauT1T2(wt,T1,T2);
 905   exbT1->SetOmegaTauT1T2(wt,T1+0.1,T2);
 906   exbT2->SetOmegaTauT1T2(wt,T1,T2+0.1);
 907   cc->Init();
 908   //  cc->SetMode(1);
 909   cc->Print("DA"); // Print used correction classes
 910   cc->SetName("Comp");
 911   return cc;
 912 }
 913
 914
 915
 916 void AddEffectiveCorrection(AliTPCComposedCorrection* comp){
 917   //
 918   //
 919   //
 920   TMatrixD polA(18,4);
 921   TMatrixD polC(18,4);
 922   TMatrixD valA(18,1);
 923   TMatrixD valC(18,1);
 924   Int_t counter=0;
 925   { for (Int_t iside=0; iside<=1; iside++){
 926       {for (Int_t ir=0; ir<3; ir++)
 927           for (Int_t id=0; id<3; id++) {
 928             polA(counter,0)=0;
 929             polA(counter,1)=ir;
 930             polA(counter,2)=id;
 931             polC(counter,0)=0;
 932             polC(counter,1)=ir;
 933             polC(counter,2)=id;
 934             valA(counter,0)=0;
 935             valC(counter,0)=0;
 936             counter++;
 937           }
 938       }
 939     }
 940   }
 941   counter=0;
 942   TObjArray * array = (TObjArray*) comp->GetCorrections();
 943   { for (Int_t iside=0; iside<=1; iside++)
 944       for (Int_t ir=0; ir<3; ir++)
 945         for (Int_t id=0; id<3; id++) {
 946           AliTPCExBEffective *eff=new AliTPCExBEffective;
 947           valA*=0;
 948           valC*=0;
 949           eff->SetName(Form("eff%d_%d_%d",iside,ir,id));
 950           if (iside==0) valA(counter,0)=0.1;
 951           if (iside==1) valC(counter,0)=0.1;
 952           eff->SetPolynoms(&polA,&polC);
 953           eff->SetCoeficients(&valA,&valC);
 954           valA.Print();
 955           valC.Print();
 956           array->AddLast(eff);
 957           counter++;
 958         }
 959   }
 960 }
 961
 962
 963
 964 void MakeLaserTree(){
 965   //
 966   //
 967   //
 968   AliTPCComposedCorrection *cc= MakeComposedCorrection();
 969   TObjArray * corr = (TObjArray*)(cc->GetCorrections());
 970   //  corr->AddLast(cc);
 971   TFile f("laserMean.root");
 972   TTree* tree =(TTree*)f.Get("Mean");
 973   AliTPCCorrection::MakeLaserDistortionTree(tree,corr,0);
 974
 975 }