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