TPC/CalibMacros/AnalyzeLaserCE.C

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