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