1. Modification in order to use it batch mode
[u/mrichter/AliRoot.git] / TPC / CalibMacros / AnalyzeLaserCE.C
index c1a034cedb190b2ea8f67b900501818e27975364..4e88e0be3a84a15bd6e23893f0901c7d0693d8a8 100644 (file)
@@ -3,9 +3,9 @@
 Macro to perform fits of the Laser Central electrode data
 Several fit methods implemented
 
-0. RebuildCE("ce.root"); - rebuild data from the scratch
-   
-
+0. RebuildCE("ce.root","pul.root"); - rebuild data from the scratch
+                                    - the data will be storered in file inname
+                                    
 1. RebuildData() - transform arbitrary layout of the Input data to the internal format
    StoreData();  - The data tree expected in file inname (see variable bellow)
    StoreTree();  - Modify inname and xxside and tcor in order to transform data
@@ -24,6 +24,7 @@ Several fit methods implemented
 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC -I$ALICE_ROOT/STAT");
 gSystem->Load("libSTAT.so");
 .L $ALICE_ROOT/TPC/CalibMacros/AnalyzeLaserCE.C+
+
 // setup aliases
 qaside="CE_Q";
 taside="CE_T";
@@ -67,10 +68,10 @@ Gloabl fit o consist of
 // Control variable - check results
 //
 //
-ffit2~-(timeIn~+tcor~):lx~  - fit value minus input time 
+ffit2~-(timeIn~):lx~  - fit value minus input time 
 
 result cosntruction:
-(timeF2~-ffit2~+fTL~+fInOut~+tcor~):Result~+tcor~
+(timeF2~-ffit2~+fTL~+fInOut~):Result~
 //
 timeF2~-Result~:ffit2~-fTL~-fInOut~
 
@@ -79,11 +80,14 @@ timeF2~-Result~:ffit2~-fTL~-fInOut~
 #include "TString.h"
 #include "TSystem.h"
 #include "TTree.h"
+#include "TCut.h"
 #include "TStatToolkit.h"
 #include "AliTPCCalibViewer.h"
 #include "AliTPCCalibViewerGUI.h"
 #include "AliTPCPreprocessorOnline.h"
 #include "AliTPCCalibCE.h"
+#include "AliTPCCalibPulser.h"
+#include "TStopwatch.h"
 //
 //Define interesting variables - file names
 //
@@ -91,16 +95,17 @@ char * inname = "treeCE.root";  // input file with tree
 //
 // variable name definition in input tree - change it according the input
 //
-TString qaside("CE_A00_Q_05");
-TString taside("CE_A00_T_05");
-TString raside("CE_A00_rms_05");
-TString qcside("CE_C00_Q_05");
-TString tcside("CE_C00_T_05");
-TString rcside("CE_C00_rms_05");
+TString qaside("CE_Q");
+TString taside("CE_T");
+TString raside("CE_RMS");
+TString qcside("CE_Q");
+TString tcside("CE_T");
+TString rcside("CE_RMS");
+
 //
 // correction variable - usually Pulser time
 //
-TString tcor("(sector%36>30)*2");
+TString tcor("-pulCorr");
 
 //
 char * fname  = "treefitCE.root";       // output file with tree
@@ -122,7 +127,7 @@ AliTPCCalPad *calPadOut = 0;            // outlyer CalPad
 //
 // cuts
 //
-const Float_t tThr=0.5;                // max diff in the sector
+const Float_t tThr=0.3;                // max diff in the sector
 const Float_t qThr0=0.5;               // max Q diff in the sector
 const Float_t qThr1=2;                 // max Q diff in the sector
 
@@ -141,7 +146,12 @@ AliTPCCalPad *calPadGXY   = 0;          // global XY missalign (drift velocity g
 AliTPCCalPad *calPadOff   = 0;          // normalization offset fit
 AliTPCCalPad *calPadRes   = 0;          // final calibration  
 
-
+TObjString strFit0="";                     // string fit 0
+TObjString strFit1="";                     // string fit 1
+TObjString strFit2="";                     // string fit 2
+TVectorD vecFit0;                       // parameters fit 0
+TVectorD vecFit1;                       // parameters fit 1
+TVectorD vecFit2;                       // parameters fit 2
 //
 // working variables
 //
@@ -152,6 +162,7 @@ TTree * tree=0;                    // working tree
 void LoadViewer();
 void RebuildData();   // transform the input data to the fit format 
 void MakeFit();       // make fits
+void MakeFitPulser(); // make fit for pulser correction data
 //
 //   internal functions
 //
@@ -162,6 +173,7 @@ void StoreData();     // store current data
 void StoreTree();     // store fit data in the output tree
 
 
+
 void AnalyzeLaser(){
   //
   //
@@ -170,16 +182,46 @@ void AnalyzeLaser(){
   MakeAliases1();
 }
 
+void MakeFitPulser(){
+  TStatToolkit stat;
+  Int_t npoints;
+  Double_t chi2;
+  TVectorD vec0,vec1,vec2;
+  TMatrixD mat;
+  TString fitvar="P_T.fElements";
+  TCut out("abs(P_T.fElements/P_T_Mean.fElements-1.001)<.002");
+  TCut fitadd("P_T.fElements>446");
+  TString fstring="";
+  fstring+="(sector>=0&&sector<9)++";
+  fstring+="(sector>=9&&sector<18)++";
+  fstring+="(sector>=18&&sector<27)++";
+  fstring+="(sector>=27&&sector<36)++";
+  fstring+="(sector>=36&&sector<45)++";
+  fstring+="(sector>=45&&sector<54)++";
+  fstring+="(sector>=54&&sector<63)++";
+  fstring+="(sector>=63&&sector<72)";
+
+  TString *pulOff =stat.FitPlane(tree,fitvar.Data(),fstring.Data(),(out+fitadd).GetTitle(),chi2,npoints,vec0,mat);
+
+  tree->SetAlias("pul0",pulOff->Data());
+  tree->SetAlias("pulCorr","P_T.fElements-pul0");
+  tree->SetAlias("pulOut",out.GetTitle());
+}
 
 void MakeFit(){
   //
   //
   LoadData();
-  LoadViewer();
+  //LoadViewer();
+  //
+  makePad = new AliTPCCalibViewer(fname);
+  tree = makePad->GetTree();
+  MakeAliases1();
+  //
+  //  MakeFitPulser();
   TStatToolkit stat;
   Int_t npoints;
   Double_t chi2;
-  TVectorD vec0,vec1,vec2;
   TMatrixD mat;
   TString fstring="";
   //
@@ -204,46 +246,18 @@ void MakeFit(){
   fstring+="(side*inn*lxr)++";//                        - opposite      //14
   fstring+="(lxr^2)++";       // zr          second     - common        //15
   fstring+="(side*lxr^2)++";  //                        - opposite      //16
+  fstring+="(inn*lxr^2)++";       // zr          second     - common        //15
+  fstring+="(inn*side*lxr^2)++";  //                        - opposite      //16
   //
-  TString *fit0 =stat.FitPlane(tree,"dt",fstring.Data(),"cutF&&cutCE",chi2,npoints,vec0,mat);
+  printf("Fit0\t start\n");
+  TString *fit0 =stat.FitPlane(tree,"dt",fstring.Data(),"cutF&&cutCE",chi2,npoints,vecFit0,mat,0.9);
   tree->SetAlias("f0",fit0->Data());
+  strFit0 = fit0->Data();
+  printf("Fit0\t end\n");
+  printf("Chi2/npoints\t=%f\n",TMath::Sqrt(chi2/npoints));
+  fit0->Tokenize("++")->Print();
   //
-  // Common "deformation" tendencies
-  //
-  fstring+="(sin(atan2(gy.fElements,gx.fElements)))++";
-  fstring+="(cos(atan2(gy.fElements,gx.fElements)))++";
-  //
-  fstring+="(sin(atan2(gy.fElements,gx.fElements)*2))++";
-  fstring+="(cos(atan2(gy.fElements,gx.fElements)*2))++";
-  fstring+="(sin(atan2(gy.fElements,gx.fElements)*3))++";
-  fstring+="(cos(atan2(gy.fElements,gx.fElements)*3))++";
-  //
-  fstring+="(sin(atan2(gy.fElements,gx.fElements)*2))*lxr++";
-  fstring+="(cos(atan2(gy.fElements,gx.fElements)*2))*lxr++";
-  fstring+="(sin(atan2(gy.fElements,gx.fElements)*3))*lxr++";
-  fstring+="(cos(atan2(gy.fElements,gx.fElements)*3))*lxr++";
-  //
-
-  TString *fit1 =stat.FitPlane(tree,"dt",fstring.Data(),"cutF&&cutCE",chi2,npoints,vec1,mat);
-  tree->SetAlias("f1",fit1->Data());
-  //
-  // Central electrode "deformation"
-  //
-  fstring+="(side*sin(atan2(gy.fElements,gx.fElements)))++";
-  fstring+="(side*cos(atan2(gy.fElements,gx.fElements)))++";
-  //
-  fstring+="(side*sin(atan2(gy.fElements,gx.fElements)*2))++";
-  fstring+="(side*cos(atan2(gy.fElements,gx.fElements)*2))++";
-  fstring+="(side*sin(atan2(gy.fElements,gx.fElements)*3))++";
-  fstring+="(side*cos(atan2(gy.fElements,gx.fElements)*3))++";
-  // 
-  fstring+="(side*sin(atan2(gy.fElements,gx.fElements)*2))*lxr++";
-  fstring+="(side*cos(atan2(gy.fElements,gx.fElements)*2))*lxr++";
-  fstring+="(side*sin(atan2(gy.fElements,gx.fElements)*3))*lxr++";
-  fstring+="(side*cos(atan2(gy.fElements,gx.fElements)*3))*lxr++";
-   
-  TString *fit2 =stat.FitPlane(tree,"dt",fstring.Data(),"cutF&&abs(dt-f0)<0.7&&cutCE",chi2,npoints,vec2,mat);
-  tree->SetAlias("f2",fit2->Data());
+  printf("Global tendencies extraction\n");
   //
   // Extract variables
   //
@@ -257,9 +271,9 @@ void MakeFit(){
   TString fitTL("0");  // side offsets
   //
   fitOff+="+";
-  fitOff+=vec2[0];
+  fitOff+=vecFit0[0];
   fitOff+="+side*";
-  fitOff+=vec2[1];
+  fitOff+=vecFit0[1];
   {
   for(Int_t i=0;i<arr->GetEntriesFast();i++){
     if (!arr->At(i)) continue;
@@ -273,30 +287,24 @@ void MakeFit(){
     if (fitstr->Contains("tl^2")){
       fitTL+="+";
       fitTL+=(*fitstr)+"*";
-      fitTL+=vec2[i+1];
+      fitTL+=vecFit0[i+1];
     }
     if (isGXY){
       fitGXY+="+";
       fitGXY+=(*fitstr)+"*";
-      fitGXY+=vec2[i+1];
+      fitGXY+=vecFit0[i+1];
     }
-    //if (isQ){
-    //  //
-    //  fitQ+="+";
-    //  fitQ+=(*fitstr)+"*";
-    //  fitQ+=vec2[i+1];
-    //}
     //
     if (isLX&&!isRot&&!isIn){
       fitLX+="+";
       fitLX+=(*fitstr)+"*";
-      fitLX+=vec2[i+1];
+      fitLX+=vecFit0[i+1];
     }
     //
     if (!isRot&&isIn){
       fitInOut+="+";
       fitInOut+=(*fitstr)+"*";
-      fitInOut+=vec2[i+1];
+      fitInOut+=vecFit0[i+1];
     }
   }
   }
@@ -307,10 +315,57 @@ void MakeFit(){
   tree->SetAlias("fOff",fitOff.Data());
   //tree->SetAlias("fQ",fitQ.Data());
   tree->SetAlias("fTL",fitTL.Data());
+
+  //
   //
+  // Common "deformation" tendencies
+  //
+  fstring+="(sin(atan2(gy.fElements,gx.fElements)))++";
+  fstring+="(cos(atan2(gy.fElements,gx.fElements)))++";
+  //
+  fstring+="(sin(atan2(gy.fElements,gx.fElements)*2))++";
+  fstring+="(cos(atan2(gy.fElements,gx.fElements)*2))++";
+  fstring+="(sin(atan2(gy.fElements,gx.fElements)*3))++";
+  fstring+="(cos(atan2(gy.fElements,gx.fElements)*3))++";
+  //
+  fstring+="(sin(atan2(gy.fElements,gx.fElements)*2))*lxr++";
+  fstring+="(cos(atan2(gy.fElements,gx.fElements)*2))*lxr++";
+  fstring+="(sin(atan2(gy.fElements,gx.fElements)*3))*lxr++";
+  fstring+="(cos(atan2(gy.fElements,gx.fElements)*3))*lxr++";
   //
-  // fits
+
+  TString *fit1 =stat.FitPlane(tree,"dt",fstring.Data(),"cutF&&cutCE",chi2,npoints,vecFit1,mat,0.9);
+  tree->SetAlias("f1",fit1->Data());
+  strFit1 = fit1->Data();
+  printf("Fit1\t end\n");
+  printf("Chi2/npoints\t=%f\n",TMath::Sqrt(chi2/npoints));
+  fit1->Tokenize("++")->Print();
+
+  //
+  // Central electrode "deformation"
+  //
+  fstring+="(side*sin(atan2(gy.fElements,gx.fElements)))++";
+  fstring+="(side*cos(atan2(gy.fElements,gx.fElements)))++";
+  //
+  fstring+="(side*sin(atan2(gy.fElements,gx.fElements)*2))++";
+  fstring+="(side*cos(atan2(gy.fElements,gx.fElements)*2))++";
+  fstring+="(side*sin(atan2(gy.fElements,gx.fElements)*3))++";
+  fstring+="(side*cos(atan2(gy.fElements,gx.fElements)*3))++";
   // 
+  fstring+="(side*sin(atan2(gy.fElements,gx.fElements)*2))*lxr++";
+  fstring+="(side*cos(atan2(gy.fElements,gx.fElements)*2))*lxr++";
+  fstring+="(side*sin(atan2(gy.fElements,gx.fElements)*3))*lxr++";
+  fstring+="(side*cos(atan2(gy.fElements,gx.fElements)*3))*lxr++";
+   
+  TString *fit2 =stat.FitPlane(tree,"dt",fstring.Data(),"cutF&&abs(dt-f0)<0.7&&cutCE",chi2,npoints,vecFit2,mat,0.9);
+  strFit2 = fit2->Data();
+  printf("Fit2\t end\n");
+  printf("Chi2/npoints\t=%f\n",TMath::Sqrt(chi2/npoints));
+  fit2->Tokenize("++")->Print();
+  tree->SetAlias("f2",fit2->Data());
+  //
+  //
+  //
   calPad0      = makePad->GetCalPad("f0","1", "ffit0");
   calPad1      = makePad->GetCalPad("f1","1", "ffit1");
   calPad2      = makePad->GetCalPad("f2","1", "ffit2");
@@ -322,6 +377,8 @@ void MakeFit(){
   calPadOff    = makePad->GetCalPad("fOff","1", "fOff");  
 }
 
+
+
 void LoadViewer(){
   //
   // Load calib Viewer
@@ -342,21 +399,32 @@ void RebuildData(){
   //
   // transform the input data to the fit format 
   //
+  TStopwatch timer;
   makePad = new AliTPCCalibViewer(inname);
   tree = makePad->GetTree();
+  //
+
+  timer.Print();
+  timer.Continue();
+  printf("-1.   MaQkeFitPulser\n");
+  MakeFitPulser();
+  timer.Print();
+  timer.Continue();
   MakeAliases0(); //
   //
   printf("0.   GetCalPads\n");
-
+  timer.Print();
   calPadCor = makePad->GetCalPad("tcor","1", "tcor");
-  calPadOut = makePad->GetCalPad("1","!((cutA||cutC)&&abs(ly.fElements/lx.fElements)<0.155)", "out");
-  calPadIn  = makePad->GetCalPad("dt-tcor","1","timeIn");
-  calPadF1  = calPadIn->GlobalFit("timeF1", calPadOut,kTRUE,0,0.9);
+  calPadOut = makePad->GetCalPad("1","!((cutA||cutC)&&abs(ly.fElements/lx.fElements)<0.16)", "out");
+  calPadIn  = makePad->GetCalPad("dt","1","timeIn");
+  
+  calPadF1  = calPadIn->GlobalFit("timeF1", calPadOut,kFALSE,0,0,0.8);  // use robust fit  
   calPadQIn = makePad->GetCalPad("qa*(sector%36<18)+qc*(sector%36>17)","1","qIn");
   //
   //
   printf("1.   Create outlyer map\n");
-
+  timer.Print(); 
+  timer.Continue();
   //
   // Update outlyer maps
   //
@@ -367,35 +435,42 @@ void RebuildData(){
       if (TMath::Abs(val0-val1)>tThr) calPadOut->GetCalROC(isector)->SetValue(ich,1);
     }
   }
-  printf("1.   Fit sector parabolas\n");
-
-  calPadF1  = calPadIn->GlobalFit("timeF1", calPadOut,kTRUE,0,0.9);
-  calPadF2  = calPadIn->GlobalFit("timeF2", calPadOut,kTRUE,1,0.9);
-  calPadQF1 = calPadQIn->GlobalFit("qF1", calPadOut,kTRUE,0,0.9);
+  printf("2.   Fit sector parabolas\n");
+  timer.Print();
+  timer.Continue();
+  calPadF1  = calPadIn->GlobalFit("timeF1", calPadOut,kFALSE,0,0,0.9);
+  calPadF2  = calPadIn->GlobalFit("timeF2", calPadOut,kFALSE,1,0,0.9);
+  calPadQF1 = calPadQIn->GlobalFit("qF1", calPadOut,kFALSE,1);
   calPadQF2 = calPadQIn->GlobalFit("qF2", calPadOut,kFALSE,1);
 
   //
   // Update outlyer maps
   //
-  printf("2.   Update outlyer map\n");
+  printf("3.   Update outlyer map\n");
   for (Int_t isector=0;isector<72; isector++){
     for (UInt_t ich=0;ich<calPadIn->GetCalROC(isector)->GetNchannels();ich++){
       Float_t val0= calPadQIn->GetCalROC(isector)->GetValue(ich);
       Float_t val1= calPadQF2->GetCalROC(isector)->GetValue(ich);
-      if (val1<0.1)  {
-       calPadOut->GetCalROC(isector)->SetValue(ich,1);
+      if (val1<=0)  {
        continue;
       }
       if (TMath::Abs(val0/val1)<qThr0) calPadOut->GetCalROC(isector)->SetValue(ich,1);
       if (TMath::Abs(val0/val1)>qThr1) calPadOut->GetCalROC(isector)->SetValue(ich,1);
     }
   }
-  printf("3.  Redo fit of the of parabola \n");
+  printf("4.  Redo fit of the of parabola \n");
+  timer.Print();
+  timer.Continue();
+  //
+  //
+  AliTPCCalPad *calPadInCor  = makePad->GetCalPad("dt","1","timeIn");
+  calPadF1  = calPadInCor->GlobalFit("timeF1", calPadOut,kFALSE,0,0.9);
+  calPadF2  = calPadInCor->GlobalFit("timeF2", calPadOut,kFALSE,1,0.9);
+  calPadQF1 = calPadQIn->GlobalFit("qF1", calPadOut,kFALSE,0,0,0.9);
+  calPadQF2 = calPadQIn->GlobalFit("qF2", calPadOut,kFALSE,1,0,0.9);
+  printf("5. End\n");
+  timer.Print();
 
-  calPadF1  = calPadIn->GlobalFit("timeF1", calPadOut,kTRUE,0,0.9);
-  calPadF2  = calPadIn->GlobalFit("timeF2", calPadOut,kTRUE,1,0.9);
-  calPadQF1 = calPadQIn->GlobalFit("qF1", calPadOut,kTRUE,0,0.9);
-  calPadQF2 = calPadQIn->GlobalFit("qF2", calPadOut,kFALSE,1);
 }
 
 void LoadData(){
@@ -503,7 +578,7 @@ void MakeAliases0(){
   tree->SetAlias("dt","(ta)*(sector%36<18)+(tc)*(sector%36>17)+tcor");
   tree->SetAlias("cutA","qa>30&&qa<400&&abs(ta)<2&&ra>0.5&&ra<2");
   tree->SetAlias("cutC","qc>30&&qc<400&&abs(tc)<2&&rc>0.5&&rc<2");
-  tree->SetAlias("cutF","(pad.fElements%4==0)&&(row.fElements%3==0)");
+  tree->SetAlias("cutF","((row.fElements+pad.fElements+sector)%19==0)");  // use just part of the statistic - 5 %
   tree->SetAlias("cutCE","V.out.fElements");
   //
   // fit param aliases
@@ -523,11 +598,12 @@ void MakeAliases1(){
   //
   tree->SetAlias("tcor","tcor.fElements");          // correction variable  
   tree->SetAlias("side","1-(sector%36>17)*2");
-  tree->SetAlias("dt","timeIn.fElements+tcor");
+  tree->SetAlias("dt","timeIn.fElements");
   //
   tree->SetAlias("cutA","out.fElements==1");
   tree->SetAlias("cutC","out.fElements==1");
-  tree->SetAlias("cutF","(pad.fElements%5==0)&&(row.fElements%4==0)");
+  tree->SetAlias("cutF","((row.fElements+pad.fElements+sector.fElements)%19==0)");  // use just part of the statistic - 5 %
+
   tree->SetAlias("cutCE","out.fElements<0.5");
   //
   // fit param aliases
@@ -547,7 +623,6 @@ void MakeRes()
   // make final calibration
   //
   AliTPCCalPad * calPadRes0 =new AliTPCCalPad(*calPadIn);
-  calPadRes0->Add(calPadCor);       // add correction
   calPadRes0->Add(calPad2,-1);      // remove global fit
   calPadRes  = calPadRes0->GlobalFit("Result", calPadOut,kTRUE,1,0.9);
   //
@@ -568,13 +643,12 @@ void MakeRes()
     }
   }
   calPadRes->Add(calPadCor,-1);     // remove back correction (e.g Pulser time 0)
-
 }
 
 
 
 
-void RebuildCE(char *finname){
+void RebuildCE(char *finname, char *pulname){
   //
   // Transformation from the CE to the visualization-analisys output
   //
@@ -588,9 +662,42 @@ void RebuildCE(char *finname){
   padtime->SetName("CE_T");
   padRMS->SetName("CE_RMS");
   padq->SetName("CE_Q");
+
+  TFile f2(pulname);
+  AliTPCCalibPulser *pul = (AliTPCCalibPulser*)f2.Get("AliTPCCalibPulser");
+  AliTPCCalPad *pultime = new AliTPCCalPad((TObjArray*)pul->GetCalPadT0());
+  AliTPCCalPad *pulRMS = new AliTPCCalPad((TObjArray*)pul->GetCalPadRMS());
+  AliTPCCalPad *pulq = new AliTPCCalPad((TObjArray*)pul->GetCalPadQ());
+  pultime->SetName("P_T");
+  pulRMS->SetName("P_RMS");
+  pulq->SetName("P_Q");
+
+
   AliTPCPreprocessorOnline * preprocesor = new AliTPCPreprocessorOnline;
   preprocesor->AddComponent(padtime);
   preprocesor->AddComponent(padq);
   preprocesor->AddComponent(padRMS);
+  preprocesor->AddComponent(pultime);
+  preprocesor->AddComponent(pulq);
+  preprocesor->AddComponent(pulRMS);
   preprocesor->DumpToFile(inname);
 }
+
+
+void AnalyzeLaserCE(){
+  RebuildCE("ce.root","pul.root");
+  StoreData();
+  StoreTree();
+  RebuildData();
+  StoreData();
+  StoreTree();
+  MakeFit();
+  StoreData();
+  StoreTree();
+  MakeRes();
+}
+
+//void AddFile(char *fname, char *aname){
+//  TFile f(fname);
+//  TTree 
+//}