4e88e0be3a84a15bd6e23893f0901c7d0693d8a8
[u/mrichter/AliRoot.git] / TPC / CalibMacros / AnalyzeLaserCE.C
1 /*
2
3 Macro to perform fits of the Laser Central electrode data
4 Several fit methods implemented
5
6 0. RebuildCE("ce.root","pul.root"); - rebuild data from the scratch
7                                     - the data will be storered in file inname
8                                     
9 1. RebuildData() - transform arbitrary layout of the Input data to the internal format
10    StoreData();  - The data tree expected in file inname (see variable bellow)
11    StoreTree();  - Modify inname and xxside and tcor in order to transform data
12
13
14 2. MakeFit();    - Make a fit of the data - already in internal format    
15    StoreData();  - Store
16    StoreTree();
17
18 3. MakeRes();    - Make the final calibration  + combination of different components
19
20 4. LoadViewer(); - Browse the fit parameters
21  
22
23 .x ~/rootlogon.C
24 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC -I$ALICE_ROOT/STAT");
25 gSystem->Load("libSTAT.so");
26 .L $ALICE_ROOT/TPC/CalibMacros/AnalyzeLaserCE.C+
27
28 // setup aliases
29 qaside="CE_Q";
30 taside="CE_T";
31 raside="CE_RMS";
32 qcside="CE_Q";
33 tcside="CE_T";
34 rcside="CE_RMS";
35
36
37
38
39 Calibration viewer variables:
40
41 Result  -  resulting correction
42 out     -  outlyers not used for fit
43 tcor    -  offset specified by user before fitting
44 timeF1  -  sector local fit - plane
45 timeF2  -  sector local fit - parabola
46 timeIn  -  input times
47 qIn     -  input charge
48 out     -  outlyers not used for fit
49 tcor    -  offset specified by user before fitting
50 timeF1  -  sector time local fit - plane
51 timeF2  -  sector time local fit - parabola
52 qF1     -  sector q local fit    - plane
53 qF2     -  sector q local fit    - parabola
54 // fitted values
55 //
56 ffit0   - base fit
57 ffit1   - adding common shifts    - alpha dependendent
58 ffit2   - adding opposite shifts  - alpha dependent
59 //
60 fGXY    -  global fit parameter - XY
61 fInOut  -  global fit parameter - inner-outer sector matching
62 fLX     -  global LX  dependence
63 //
64 Gloabl fit o consist of
65 -fGXY~-fLX~-fTL~-fOff~:ffit0~
66
67 //
68 // Control variable - check results
69 //
70 //
71 ffit2~-(timeIn~):lx~  - fit value minus input time 
72
73 result cosntruction:
74 (timeF2~-ffit2~+fTL~+fInOut~):Result~
75 //
76 timeF2~-Result~:ffit2~-fTL~-fInOut~
77
78
79 */
80 #include "TString.h"
81 #include "TSystem.h"
82 #include "TTree.h"
83 #include "TCut.h"
84 #include "TStatToolkit.h"
85 #include "AliTPCCalibViewer.h"
86 #include "AliTPCCalibViewerGUI.h"
87 #include "AliTPCPreprocessorOnline.h"
88 #include "AliTPCCalibCE.h"
89 #include "AliTPCCalibPulser.h"
90 #include "TStopwatch.h"
91 //
92 //Define interesting variables - file names
93 //
94 char * inname = "treeCE.root";  // input file with tree
95 //
96 // variable name definition in input tree - change it according the input
97 //
98 TString qaside("CE_Q");
99 TString taside("CE_T");
100 TString raside("CE_RMS");
101 TString qcside("CE_Q");
102 TString tcside("CE_T");
103 TString rcside("CE_RMS");
104
105 //
106 // correction variable - usually Pulser time
107 //
108 TString tcor("-pulCorr");
109
110 //
111 char * fname  = "treefitCE.root";       // output file with tree
112 char * oname  = "fitCE.root";           // output file with CalPads fit 
113
114 //
115 //
116 // Input CalPads
117 //
118 AliTPCCalPad *calPadIn  = 0;            // original time pad
119 AliTPCCalPad *calPadF1  = 0;            // original time pad - fit plane
120 AliTPCCalPad *calPadF2  = 0;            // original time pad - fit parabola
121 AliTPCCalPad *calPadQIn = 0;            // original Q pad
122 AliTPCCalPad *calPadQF1 = 0;            // original Q pad
123 AliTPCCalPad *calPadQF2 = 0;            // original Q pad
124
125 AliTPCCalPad *calPadCor = 0;            // base correction CalPad
126 AliTPCCalPad *calPadOut = 0;            // outlyer CalPad
127 //
128 // cuts
129 //
130 const Float_t tThr=0.3;                // max diff in the sector
131 const Float_t qThr0=0.5;               // max Q diff in the sector
132 const Float_t qThr1=2;                 // max Q diff in the sector
133
134 //
135 //
136 // fit Cal Pads
137 AliTPCCalPad *calPad0   = 0;            // global fit 0 - base
138 AliTPCCalPad *calPad1   = 0;            // global fit 1 - common behavior rotation -A-C
139 AliTPCCalPad *calPad2   = 0;            // gloabl fit 2 - CE missalign rotation     A-C
140 //
141 AliTPCCalPad *calPadInOut = 0;          // misaalign in-out
142 AliTPCCalPad *calPadLX    = 0;          // local x missalign
143 AliTPCCalPad *calPadTL   = 0;           // tan 
144 AliTPCCalPad *calPadQ     = 0;          // time (q)  correction
145 AliTPCCalPad *calPadGXY   = 0;          // global XY missalign (drift velocity grad) 
146 AliTPCCalPad *calPadOff   = 0;          // normalization offset fit
147 AliTPCCalPad *calPadRes   = 0;          // final calibration  
148
149 TObjString strFit0="";                     // string fit 0
150 TObjString strFit1="";                     // string fit 1
151 TObjString strFit2="";                     // string fit 2
152 TVectorD vecFit0;                       // parameters fit 0
153 TVectorD vecFit1;                       // parameters fit 1
154 TVectorD vecFit2;                       // parameters fit 2
155 //
156 // working variables
157 //
158 AliTPCCalibViewerGUI * viewer=0;   //viewerGUI
159 AliTPCCalibViewer    *makePad=0;   //viewer
160 TTree * tree=0;                    // working tree
161
162 void LoadViewer();
163 void RebuildData();   // transform the input data to the fit format 
164 void MakeFit();       // make fits
165 void MakeFitPulser(); // make fit for pulser correction data
166 //
167 //   internal functions
168 //
169 void MakeAliases0();  // Make Aliases 0  - for user tree
170 void MakeAliases1();  // Make Aliases 1  - for default tree   
171 void LoadData();      // Load data
172 void StoreData();     // store current data
173 void StoreTree();     // store fit data in the output tree
174
175
176
177 void AnalyzeLaser(){
178   //
179   //
180   //
181   LoadViewer();
182   MakeAliases1();
183 }
184
185 void MakeFitPulser(){
186   TStatToolkit stat;
187   Int_t npoints;
188   Double_t chi2;
189   TVectorD vec0,vec1,vec2;
190   TMatrixD mat;
191   TString fitvar="P_T.fElements";
192   TCut out("abs(P_T.fElements/P_T_Mean.fElements-1.001)<.002");
193   TCut fitadd("P_T.fElements>446");
194   TString fstring="";
195   fstring+="(sector>=0&&sector<9)++";
196   fstring+="(sector>=9&&sector<18)++";
197   fstring+="(sector>=18&&sector<27)++";
198   fstring+="(sector>=27&&sector<36)++";
199   fstring+="(sector>=36&&sector<45)++";
200   fstring+="(sector>=45&&sector<54)++";
201   fstring+="(sector>=54&&sector<63)++";
202   fstring+="(sector>=63&&sector<72)";
203
204   TString *pulOff =stat.FitPlane(tree,fitvar.Data(),fstring.Data(),(out+fitadd).GetTitle(),chi2,npoints,vec0,mat);
205
206   tree->SetAlias("pul0",pulOff->Data());
207   tree->SetAlias("pulCorr","P_T.fElements-pul0");
208   tree->SetAlias("pulOut",out.GetTitle());
209 }
210
211 void MakeFit(){
212   //
213   //
214   LoadData();
215   //LoadViewer();
216   //
217   makePad = new AliTPCCalibViewer(fname);
218   tree = makePad->GetTree();
219   MakeAliases1();
220   //
221   //  MakeFitPulser();
222   TStatToolkit stat;
223   Int_t npoints;
224   Double_t chi2;
225   TMatrixD mat;
226   TString fstring="";
227   //
228   //Basic  correction
229   //
230   fstring+="side++";        // offset on 2 different sides              //1
231   //fstring+="(1/qp)++";      // Q -threshold effect correction           //2
232   //fstring+="(qp)++";        // Q -threshold effect correction           //3
233   fstring+="(inn)++";       //  inner outer misalign   - common         //4 
234   fstring+="(side*inn)++";  //                         - opposite       //5
235   //
236   fstring+="(gyr)++";       // drift velocity gradient - common         //6
237   fstring+="(side*gyr)++";  //                         - opposite       //7
238   fstring+="(gxr)++";       //  global x tilt          - common         //8
239   fstring+="(side*gxr)++";  //                         - opposite       //9
240   //
241   fstring+="tl^2++";        // local phi correction                     //10
242   //
243   fstring+="(lxr)++";       // zr            angle      - common        //11
244   fstring+="(side*lxr)++";  //                          - opposite      //12
245   fstring+="(inn*lxr)++";   // inner outer angle        - common        //13             
246   fstring+="(side*inn*lxr)++";//                        - opposite      //14
247   fstring+="(lxr^2)++";       // zr          second     - common        //15
248   fstring+="(side*lxr^2)++";  //                        - opposite      //16
249   fstring+="(inn*lxr^2)++";       // zr          second     - common        //15
250   fstring+="(inn*side*lxr^2)++";  //                        - opposite      //16
251   //
252   printf("Fit0\t start\n");
253   TString *fit0 =stat.FitPlane(tree,"dt",fstring.Data(),"cutF&&cutCE",chi2,npoints,vecFit0,mat,0.9);
254   tree->SetAlias("f0",fit0->Data());
255   strFit0 = fit0->Data();
256   printf("Fit0\t end\n");
257   printf("Chi2/npoints\t=%f\n",TMath::Sqrt(chi2/npoints));
258   fit0->Tokenize("++")->Print();
259   //
260   printf("Global tendencies extraction\n");
261   //
262   // Extract variables
263   //
264   TString tmpstr = fstring;
265   TObjArray *arr = tmpstr.Tokenize("++");
266   TString fitQ("0");       // q correction 
267   TString fitLX("0");      // lx correction 
268   TString fitInOut("0");   // inner-outer - match
269   TString fitGXY("0");      // global xy fit
270   TString fitOff("0");  // side offsets
271   TString fitTL("0");  // side offsets
272   //
273   fitOff+="+";
274   fitOff+=vecFit0[0];
275   fitOff+="+side*";
276   fitOff+=vecFit0[1];
277   {
278   for(Int_t i=0;i<arr->GetEntriesFast();i++){
279     if (!arr->At(i)) continue;
280     TString *fitstr = new TString(arr->At(i)->GetName());
281     //
282     //Bool_t isQ      = fitstr->Contains("qp)");
283     Bool_t isRot    = fitstr->Contains("sin(")+fitstr->Contains("cos(");
284     Bool_t isLX     = fitstr->Contains("lxr");
285     Bool_t isIn     = fitstr->Contains("inn");
286     Bool_t isGXY    = fitstr->Contains("gxr")+fitstr->Contains("gyr");
287     if (fitstr->Contains("tl^2")){
288       fitTL+="+";
289       fitTL+=(*fitstr)+"*";
290       fitTL+=vecFit0[i+1];
291     }
292     if (isGXY){
293       fitGXY+="+";
294       fitGXY+=(*fitstr)+"*";
295       fitGXY+=vecFit0[i+1];
296     }
297     //
298     if (isLX&&!isRot&&!isIn){
299       fitLX+="+";
300       fitLX+=(*fitstr)+"*";
301       fitLX+=vecFit0[i+1];
302     }
303     //
304     if (!isRot&&isIn){
305       fitInOut+="+";
306       fitInOut+=(*fitstr)+"*";
307       fitInOut+=vecFit0[i+1];
308     }
309   }
310   }
311   //
312   tree->SetAlias("fInOut",fitInOut.Data());
313   tree->SetAlias("fLX",fitLX.Data());
314   tree->SetAlias("fGXY",fitGXY.Data());
315   tree->SetAlias("fOff",fitOff.Data());
316   //tree->SetAlias("fQ",fitQ.Data());
317   tree->SetAlias("fTL",fitTL.Data());
318
319   //
320   //
321   // Common "deformation" tendencies
322   //
323   fstring+="(sin(atan2(gy.fElements,gx.fElements)))++";
324   fstring+="(cos(atan2(gy.fElements,gx.fElements)))++";
325   //
326   fstring+="(sin(atan2(gy.fElements,gx.fElements)*2))++";
327   fstring+="(cos(atan2(gy.fElements,gx.fElements)*2))++";
328   fstring+="(sin(atan2(gy.fElements,gx.fElements)*3))++";
329   fstring+="(cos(atan2(gy.fElements,gx.fElements)*3))++";
330   //
331   fstring+="(sin(atan2(gy.fElements,gx.fElements)*2))*lxr++";
332   fstring+="(cos(atan2(gy.fElements,gx.fElements)*2))*lxr++";
333   fstring+="(sin(atan2(gy.fElements,gx.fElements)*3))*lxr++";
334   fstring+="(cos(atan2(gy.fElements,gx.fElements)*3))*lxr++";
335   //
336
337   TString *fit1 =stat.FitPlane(tree,"dt",fstring.Data(),"cutF&&cutCE",chi2,npoints,vecFit1,mat,0.9);
338   tree->SetAlias("f1",fit1->Data());
339   strFit1 = fit1->Data();
340   printf("Fit1\t end\n");
341   printf("Chi2/npoints\t=%f\n",TMath::Sqrt(chi2/npoints));
342   fit1->Tokenize("++")->Print();
343
344   //
345   // Central electrode "deformation"
346   //
347   fstring+="(side*sin(atan2(gy.fElements,gx.fElements)))++";
348   fstring+="(side*cos(atan2(gy.fElements,gx.fElements)))++";
349   //
350   fstring+="(side*sin(atan2(gy.fElements,gx.fElements)*2))++";
351   fstring+="(side*cos(atan2(gy.fElements,gx.fElements)*2))++";
352   fstring+="(side*sin(atan2(gy.fElements,gx.fElements)*3))++";
353   fstring+="(side*cos(atan2(gy.fElements,gx.fElements)*3))++";
354   // 
355   fstring+="(side*sin(atan2(gy.fElements,gx.fElements)*2))*lxr++";
356   fstring+="(side*cos(atan2(gy.fElements,gx.fElements)*2))*lxr++";
357   fstring+="(side*sin(atan2(gy.fElements,gx.fElements)*3))*lxr++";
358   fstring+="(side*cos(atan2(gy.fElements,gx.fElements)*3))*lxr++";
359    
360   TString *fit2 =stat.FitPlane(tree,"dt",fstring.Data(),"cutF&&abs(dt-f0)<0.7&&cutCE",chi2,npoints,vecFit2,mat,0.9);
361   strFit2 = fit2->Data();
362   printf("Fit2\t end\n");
363   printf("Chi2/npoints\t=%f\n",TMath::Sqrt(chi2/npoints));
364   fit2->Tokenize("++")->Print();
365   tree->SetAlias("f2",fit2->Data());
366   //
367   //
368   //
369   calPad0      = makePad->GetCalPad("f0","1", "ffit0");
370   calPad1      = makePad->GetCalPad("f1","1", "ffit1");
371   calPad2      = makePad->GetCalPad("f2","1", "ffit2");
372   calPadInOut  = makePad->GetCalPad("fInOut","1", "fInOut");
373   calPadLX     = makePad->GetCalPad("fLX","1", "fLX");
374   calPadTL     = makePad->GetCalPad("fTL","1", "fTL");
375   //calPadQ      = makePad->GetCalPad("fQ","1", "fQ");
376   calPadGXY    = makePad->GetCalPad("fGXY","1", "fGXY");
377   calPadOff    = makePad->GetCalPad("fOff","1", "fOff");  
378 }
379
380
381
382 void LoadViewer(){
383   //
384   // Load calib Viewer
385   //
386   TObjArray * array = AliTPCCalibViewerGUI::ShowGUI(fname);
387   viewer = (AliTPCCalibViewerGUI*)array->At(0);
388   makePad = viewer->GetViewer();
389   tree = viewer->GetViewer()->GetTree();
390   MakeAliases1();  
391 }
392
393
394
395
396
397
398 void RebuildData(){
399   //
400   // transform the input data to the fit format 
401   //
402   TStopwatch timer;
403   makePad = new AliTPCCalibViewer(inname);
404   tree = makePad->GetTree();
405   //
406
407   timer.Print();
408   timer.Continue();
409   printf("-1.   MaQkeFitPulser\n");
410   MakeFitPulser();
411   timer.Print();
412   timer.Continue();
413   MakeAliases0(); //
414   //
415   printf("0.   GetCalPads\n");
416   timer.Print();
417   calPadCor = makePad->GetCalPad("tcor","1", "tcor");
418   calPadOut = makePad->GetCalPad("1","!((cutA||cutC)&&abs(ly.fElements/lx.fElements)<0.16)", "out");
419   calPadIn  = makePad->GetCalPad("dt","1","timeIn");
420   
421   calPadF1  = calPadIn->GlobalFit("timeF1", calPadOut,kFALSE,0,0,0.8);  // use robust fit  
422   calPadQIn = makePad->GetCalPad("qa*(sector%36<18)+qc*(sector%36>17)","1","qIn");
423   //
424   //
425   printf("1.   Create outlyer map\n");
426   timer.Print(); 
427   timer.Continue();
428   //
429   // Update outlyer maps
430   //
431   for (Int_t isector=0;isector<72; isector++){
432     for (UInt_t ich=0;ich<calPadIn->GetCalROC(isector)->GetNchannels();ich++){
433       Float_t val0= calPadIn->GetCalROC(isector)->GetValue(ich);
434       Float_t val1= calPadF1->GetCalROC(isector)->GetValue(ich);
435       if (TMath::Abs(val0-val1)>tThr) calPadOut->GetCalROC(isector)->SetValue(ich,1);
436     }
437   }
438   printf("2.   Fit sector parabolas\n");
439   timer.Print();
440   timer.Continue();
441   calPadF1  = calPadIn->GlobalFit("timeF1", calPadOut,kFALSE,0,0,0.9);
442   calPadF2  = calPadIn->GlobalFit("timeF2", calPadOut,kFALSE,1,0,0.9);
443   calPadQF1 = calPadQIn->GlobalFit("qF1", calPadOut,kFALSE,1);
444   calPadQF2 = calPadQIn->GlobalFit("qF2", calPadOut,kFALSE,1);
445
446   //
447   // Update outlyer maps
448   //
449   printf("3.   Update outlyer map\n");
450   for (Int_t isector=0;isector<72; isector++){
451     for (UInt_t ich=0;ich<calPadIn->GetCalROC(isector)->GetNchannels();ich++){
452       Float_t val0= calPadQIn->GetCalROC(isector)->GetValue(ich);
453       Float_t val1= calPadQF2->GetCalROC(isector)->GetValue(ich);
454       if (val1<=0)  {
455         continue;
456       }
457       if (TMath::Abs(val0/val1)<qThr0) calPadOut->GetCalROC(isector)->SetValue(ich,1);
458       if (TMath::Abs(val0/val1)>qThr1) calPadOut->GetCalROC(isector)->SetValue(ich,1);
459     }
460   }
461   printf("4.  Redo fit of the of parabola \n");
462   timer.Print();
463   timer.Continue();
464   //
465   //
466   AliTPCCalPad *calPadInCor  = makePad->GetCalPad("dt","1","timeIn");
467   calPadF1  = calPadInCor->GlobalFit("timeF1", calPadOut,kFALSE,0,0.9);
468   calPadF2  = calPadInCor->GlobalFit("timeF2", calPadOut,kFALSE,1,0.9);
469   calPadQF1 = calPadQIn->GlobalFit("qF1", calPadOut,kFALSE,0,0,0.9);
470   calPadQF2 = calPadQIn->GlobalFit("qF2", calPadOut,kFALSE,1,0,0.9);
471   printf("5. End\n");
472   timer.Print();
473
474 }
475
476 void LoadData(){
477   //
478   // Get Data
479   //
480   TFile f(oname);
481   calPadIn  = (AliTPCCalPad*)f.Get("timeIn");  // original time pad
482   calPadF1  = (AliTPCCalPad*)f.Get("timeF1");  // original time pad - fit plane
483   calPadF2  = (AliTPCCalPad*)f.Get("timeF2");  // original time pad - fit parabola
484   //
485   calPadQIn  = (AliTPCCalPad*)f.Get("qIn");  // original time pad
486   calPadQF1  = (AliTPCCalPad*)f.Get("qF1");  // original time pad - fit plane
487   calPadQF2  = (AliTPCCalPad*)f.Get("qF2");  // original time pad - fit parabola
488   //
489   calPadCor = (AliTPCCalPad*)f.Get("tcor");    // base correction CalPad
490   calPadOut = (AliTPCCalPad*)f.Get("out");     // outlyer CalPad  
491   //
492   calPad0   = (AliTPCCalPad*)f.Get("ffit0");   // global fit 0 - base
493   calPad1   = (AliTPCCalPad*)f.Get("ffit1");   // global fit 1 - common behavior rotation -A-C
494   calPad2   = (AliTPCCalPad*)f.Get("ffit2");   // gloabl fit 2 - CE missalign rotation     A-C
495   calPadInOut = (AliTPCCalPad*)f.Get("fInOut");// misaalign in-out
496   calPadLX    = (AliTPCCalPad*)f.Get("fLX");   // local x missalign
497   calPadTL    = (AliTPCCalPad*)f.Get("fTL");   // local y/x missalign
498   calPadQ     = (AliTPCCalPad*)f.Get("fQ");    // time (q)  correction
499   calPadGXY   = (AliTPCCalPad*)f.Get("fGXY");  // global XY missalign (drift velocity grad) 
500   calPadOff   = (AliTPCCalPad*)f.Get("fOff");  // normalization offset fit
501   calPadRes   = (AliTPCCalPad*)f.Get("Result");  //resulting calibration
502 }
503
504 void StoreData(){
505   //
506   // Store data
507   // 
508   TFile * fstore = new TFile(oname,"recreate");
509   if (calPadIn) calPadIn->Write("timeIn");   // original time pad
510   if (calPadF1) calPadF1->Write("timeF1");   // original time pad - fit plane
511   if (calPadF2) calPadF2->Write("timeF2");   // original time pad - fit parabola
512   //
513   if (calPadQIn) calPadQIn->Write("qIn");   // original time pad
514   if (calPadQF1) calPadQF1->Write("qF1");   // original time pad - fit plane
515   if (calPadQF2) calPadQF2->Write("qF2");   // original time pad - fit parabola
516   //
517   if (calPadCor) calPadCor->Write("tcor");   // base correction CalPad
518   if (calPadOut) calPadOut->Write("out");    // outlyer CalPad  
519   //
520   if (calPad0)   calPad0->Write("ffit0");    // global fit 0 - base
521   if (calPad1)   calPad1->Write("ffit1");    // global fit 1 - common behavior rotation -A-C
522   if (calPad2)   calPad2->Write("ffit2");    // gloabl fit 2 - CE missalign rotation     A-C
523   if (calPadInOut)calPadInOut->Write("fInOut");   // misaalign in-out
524   if (calPadLX)  calPadLX->Write("fLX");     // local x missalign
525   if (calPadTL)  calPadTL->Write("fTL");     // local y/x missalign
526   if (calPadQ)   calPadQ->Write("fQ");       // time (q)  correction
527   if (calPadGXY) calPadGXY->Write("fGXY");   // global XY missalign (drift velocity grad) 
528   if (calPadOff) calPadOff->Write("fOff");   // normalization offset fit
529   if (calPadRes) calPadRes->Write("Result");   //resulting calibration
530   fstore->Close();
531   delete fstore;
532 }
533
534 void StoreTree(){
535   //
536   //
537   //
538   AliTPCPreprocessorOnline * preprocesor = new AliTPCPreprocessorOnline;
539   //
540   if (calPadIn) preprocesor->AddComponent(calPadIn->Clone());
541   if (calPadF1) preprocesor->AddComponent(calPadF1->Clone());   
542   if (calPadF2) preprocesor->AddComponent(calPadF2->Clone());   
543   //
544   if (calPadQIn) preprocesor->AddComponent(calPadQIn->Clone());
545   if (calPadQF1) preprocesor->AddComponent(calPadQF1->Clone());   
546   if (calPadQF2) preprocesor->AddComponent(calPadQF2->Clone());   
547   //
548   if (calPadCor) preprocesor->AddComponent(calPadCor->Clone());   
549   if (calPadOut) preprocesor->AddComponent(calPadOut->Clone());  
550   //
551   if (calPad0)   preprocesor->AddComponent(calPad0->Clone());
552   if (calPad1)   preprocesor->AddComponent(calPad1->Clone());
553   if (calPad2)   preprocesor->AddComponent(calPad2->Clone());
554   if (calPadInOut)preprocesor->AddComponent(calPadInOut->Clone());
555   if (calPadLX)  preprocesor->AddComponent(calPadLX->Clone());
556   if (calPadTL)  preprocesor->AddComponent(calPadTL->Clone());
557   if (calPadQ)   preprocesor->AddComponent(calPadQ->Clone());
558   if (calPadGXY) preprocesor->AddComponent(calPadGXY->Clone());
559   if (calPadOff) preprocesor->AddComponent(calPadOff->Clone());
560   if (calPadRes) preprocesor->AddComponent(calPadRes->Clone());
561   preprocesor->DumpToFile(fname);
562   delete preprocesor;
563 }
564
565
566 void MakeAliases0(){
567   //
568   // Define variables and selection of outliers - for user defined tree
569   //
570   tree->SetAlias("tcor",tcor.Data());          // correction variable
571   tree->SetAlias("ta",taside+".fElements");
572   tree->SetAlias("tc",tcside+".fElements");
573   tree->SetAlias("qa",qaside+".fElements");
574   tree->SetAlias("qc",qcside+".fElements");
575   tree->SetAlias("ra",raside+".fElements");
576   tree->SetAlias("rc",rcside+".fElements");
577   tree->SetAlias("side","1-(sector%36>17)*2");
578   tree->SetAlias("dt","(ta)*(sector%36<18)+(tc)*(sector%36>17)+tcor");
579   tree->SetAlias("cutA","qa>30&&qa<400&&abs(ta)<2&&ra>0.5&&ra<2");
580   tree->SetAlias("cutC","qc>30&&qc<400&&abs(tc)<2&&rc>0.5&&rc<2");
581   tree->SetAlias("cutF","((row.fElements+pad.fElements+sector)%19==0)");  // use just part of the statistic - 5 %
582   tree->SetAlias("cutCE","V.out.fElements");
583   //
584   // fit param aliases
585   //
586   tree->SetAlias("inn","sector<36");
587   tree->SetAlias("gxr","(gx.fElements/250.)"); //
588   tree->SetAlias("gyr","(gy.fElements/250.)"); //
589   tree->SetAlias("lxr","(lx.fElements-133.41)/250.");
590   tree->SetAlias("qp","((sector%36<18)*sqrt(qa)/10.+(sector%36>17)*sqrt(qc)/10.)"); //
591   tree->SetAlias("tl","(ly.fElements/lx.fElements)/0.17");  
592 }
593
594
595 void MakeAliases1(){
596   //
597   // Define variables and selection of outliers -for default usage
598   //
599   tree->SetAlias("tcor","tcor.fElements");          // correction variable  
600   tree->SetAlias("side","1-(sector%36>17)*2");
601   tree->SetAlias("dt","timeIn.fElements");
602   //
603   tree->SetAlias("cutA","out.fElements==1");
604   tree->SetAlias("cutC","out.fElements==1");
605   tree->SetAlias("cutF","((row.fElements+pad.fElements+sector.fElements)%19==0)");  // use just part of the statistic - 5 %
606
607   tree->SetAlias("cutCE","out.fElements<0.5");
608   //
609   // fit param aliases
610   //
611   tree->SetAlias("inn","sector<36");
612   tree->SetAlias("gxr","(gx.fElements/250.)"); //
613   tree->SetAlias("gyr","(gy.fElements/250.)"); //
614   tree->SetAlias("lxr","(lx.fElements-133.41)/250.");
615   tree->SetAlias("qp","(sqrt(qIn.fElements)/10.+(out.fElements>0.5))"); //
616   tree->SetAlias("tl","(ly.fElements/lx.fElements)/0.17");  
617 }
618
619
620 void MakeRes()
621 {
622   //
623   // make final calibration
624   //
625   AliTPCCalPad * calPadRes0 =new AliTPCCalPad(*calPadIn);
626   calPadRes0->Add(calPad2,-1);      // remove global fit
627   calPadRes  = calPadRes0->GlobalFit("Result", calPadOut,kTRUE,1,0.9);
628   //
629   //
630   {
631     Float_t tlmedian =  calPadTL->GetMedian();
632     for (Int_t isector=0;isector<72; isector++){
633       for (UInt_t ich=0;ich<calPadIn->GetCalROC(isector)->GetNchannels();ich++){
634         //
635         //
636         Float_t val0 = calPadRes->GetCalROC(isector)->GetValue(ich);
637         if (TMath::Abs(val0)>0.5) calPadRes->GetCalROC(isector)->SetValue(ich,0);
638         Float_t tl = calPadTL->GetCalROC(isector)->GetValue(ich);
639         Float_t inOut = calPadInOut->GetCalROC(isector)->GetValue(ich);
640         calPadRes->GetCalROC(isector)->SetValue(ich,calPadRes->GetCalROC(isector)->GetValue(ich)+tl+inOut);
641         //
642       }
643     }
644   }
645   calPadRes->Add(calPadCor,-1);     // remove back correction (e.g Pulser time 0)
646 }
647
648
649
650
651 void RebuildCE(char *finname, char *pulname){
652   //
653   // Transformation from the CE to the visualization-analisys output
654   //
655   // finname = CE_Vscan_Run_61684-50_170.root;
656   TFile f(finname);
657   AliTPCCalibCE * ce  = (AliTPCCalibCE*)f.Get("AliTPCCalibCE");
658   //
659   AliTPCCalPad *padtime = new AliTPCCalPad((TObjArray*)ce->GetCalPadT0());
660   AliTPCCalPad *padRMS = new AliTPCCalPad((TObjArray*)ce->GetCalPadRMS());
661   AliTPCCalPad *padq = new AliTPCCalPad((TObjArray*)ce->GetCalPadQ());
662   padtime->SetName("CE_T");
663   padRMS->SetName("CE_RMS");
664   padq->SetName("CE_Q");
665
666   TFile f2(pulname);
667   AliTPCCalibPulser *pul = (AliTPCCalibPulser*)f2.Get("AliTPCCalibPulser");
668   AliTPCCalPad *pultime = new AliTPCCalPad((TObjArray*)pul->GetCalPadT0());
669   AliTPCCalPad *pulRMS = new AliTPCCalPad((TObjArray*)pul->GetCalPadRMS());
670   AliTPCCalPad *pulq = new AliTPCCalPad((TObjArray*)pul->GetCalPadQ());
671   pultime->SetName("P_T");
672   pulRMS->SetName("P_RMS");
673   pulq->SetName("P_Q");
674
675
676   AliTPCPreprocessorOnline * preprocesor = new AliTPCPreprocessorOnline;
677   preprocesor->AddComponent(padtime);
678   preprocesor->AddComponent(padq);
679   preprocesor->AddComponent(padRMS);
680   preprocesor->AddComponent(pultime);
681   preprocesor->AddComponent(pulq);
682   preprocesor->AddComponent(pulRMS);
683   preprocesor->DumpToFile(inname);
684 }
685
686
687 void AnalyzeLaserCE(){
688   RebuildCE("ce.root","pul.root");
689   StoreData();
690   StoreTree();
691   RebuildData();
692   StoreData();
693   StoreTree();
694   MakeFit();
695   StoreData();
696   StoreTree();
697   MakeRes();
698 }
699
700 //void AddFile(char *fname, char *aname){
701 //  TFile f(fname);
702 //  TTree 
703 //}