FitRodShift.C

Go to the documentation of this file.

#if !defined(__CINT__) || defined(__MAKECINT__)
#include "TFile.h"
#include "TTreeStream.h"
#include "TMath.h" 
#include "TGraph.h" 
#include "TRandom.h"
#include "TTree.h"
#include "TF1.h"
#include "TH2F.h"
#include "TH3F.h"
#include "TAxis.h"
#include "TPad.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "AliTPCParamSR.h"
#include "TDatabasePDG.h"
#include "AliTPCFCVoltError3D.h"
#include "AliTPCROCVoltError3D.h"
#include "AliTPCComposedCorrection.h"
#include "AliTPCBoundaryVoltError.h"
#include "AliCDBEntry.h"
#include "AliTPCROC.h"
#include <TStatToolkit.h>
#include "TCut.h"
#include "TGraphErrors.h"
#include "AliTPCCalibGlobalMisalignment.h"
#endif

TTreeSRedirector *pcWorkspace= 0;    // Workspace to store the fit parameters and QA histograms
TTree * treeDY= 0;                       //distortion tree - Dy
TTree * treeDZ= 0;                       //distortion tree - Dz
// fit models ...
AliTPCFCVoltError3D *rotOFC  =0;          // fit models
AliTPCFCVoltError3D *rodOFC1 =0;
AliTPCFCVoltError3D *rodOFC2 =0;
AliTPCFCVoltError3D *rotIFC  =0;
AliTPCFCVoltError3D *rodIFC1 =0;
AliTPCFCVoltError3D *rodIFC2 =0;
AliTPCFCVoltError3D *rodAll  =0;        //all lookup initialized
//
AliTPCComposedCorrection *corFit0 = 0;  // composed correction
void FitFunctionQA();
void MakeQA();
void DrawAlignParam();

void PrintFit(TString fitString, TString filter){

  TObjArray *arr = fitString.Tokenize("++");
  Int_t entries=arr->GetEntries();
  for (Int_t i=0; i<entries; i++){
    TString aaa = arr->At(i)->GetName();
    if (aaa.Contains(filter)==1){
      printf("%d\t%s\n",i,aaa.Data());
    }    
  }
}


void LoadModels(){

  Int_t volt = 1;
  TFile f("OCDB-RodShifts.root");
  AliCDBEntry *entry = (AliCDBEntry*) f.Get("AliCDBEntry");
  if (entry) { // load file    
    AliTPCComposedCorrection *cor = (AliTPCComposedCorrection*) entry->GetObject();    
    cor->Print();
    TCollection *iter = cor->GetCorrections();    
    rotOFC = (AliTPCFCVoltError3D*)iter->FindObject("rotOFC");
    rodOFC1 = (AliTPCFCVoltError3D*)iter->FindObject("rodOFC1");
    rodOFC2 = (AliTPCFCVoltError3D*)iter->FindObject("rodOFC2");
    rotIFC = (AliTPCFCVoltError3D*)iter->FindObject("rotIFC");
    rodIFC1 = (AliTPCFCVoltError3D*)iter->FindObject("rodIFC1");
    rodIFC2 = (AliTPCFCVoltError3D*)iter->FindObject("rodIFC2");
    rodAll = (AliTPCFCVoltError3D*)iter->FindObject("rodAll");
    // rotOFC->CreateHistoDZinXY(1,500,500)->Draw("surf2");
  } else {    
    // OFC 
    rotOFC = new AliTPCFCVoltError3D();
    rotOFC->SetOmegaTauT1T2(0,1,1);
    rotOFC->SetRotatedClipVoltA(1,volt);
    rotOFC->SetRotatedClipVoltC(1,volt);
    //
    rodOFC1 = new AliTPCFCVoltError3D();
    rodOFC1->SetOmegaTauT1T2(0,1,1);
    rodOFC1->SetRodVoltShiftA(18,volt);
    rodOFC1->SetRodVoltShiftC(18,volt);
    //
    rodOFC2 = new AliTPCFCVoltError3D();
    rodOFC2->SetOmegaTauT1T2(0,1,1);
    rodOFC2->SetCopperRodShiftA(18,volt);
    rodOFC2->SetCopperRodShiftC(18,volt);    
    // IFC     
    rotIFC = new AliTPCFCVoltError3D();
    rotIFC->SetOmegaTauT1T2(0,1,1);
    rotIFC->SetRotatedClipVoltA(0,volt);
    rotIFC->SetRotatedClipVoltC(0,volt);
    //
    rodIFC1 = new AliTPCFCVoltError3D();
    rodIFC1->SetOmegaTauT1T2(0,1,1);
    rodIFC1->SetRodVoltShiftA(0,volt);
    rodIFC1->SetRodVoltShiftC(0,volt);
    //
    rodIFC2 = new AliTPCFCVoltError3D();
    rodIFC2->SetOmegaTauT1T2(0,1,1);
    rodIFC2->SetCopperRodShiftA(0,volt);
    rodIFC2->SetCopperRodShiftC(0,volt);
    // dummy object with all inits
    rodAll = new AliTPCFCVoltError3D();
    Double_t volt0=0.0000000001;
    rodAll->SetRotatedClipVoltA(1,volt0);
    rodAll->SetRotatedClipVoltC(1,volt0);
    rodAll->SetRodVoltShiftA(18,volt0);
    rodAll->SetRodVoltShiftC(18,volt0);
    rodAll->SetCopperRodShiftA(18,volt0);
    rodAll->SetCopperRodShiftC(18,volt0);    
    rodAll->SetRotatedClipVoltA(0,volt0);
    rodAll->SetRotatedClipVoltC(0,volt0);
    rodAll->SetRodVoltShiftA(0,volt0);
    rodAll->SetRodVoltShiftC(0,volt0);
    rodAll->SetCopperRodShiftA(0,volt0);
    rodAll->SetCopperRodShiftC(0,volt0);
    rodAll->SetOmegaTauT1T2(0,1,1);
    //
    //
    // Initialization of the lookup tables
    //
    printf(" ------- OFC rotated clip:\n"); rotOFC->InitFCVoltError3D();
    printf(" ------- OFC rod & strip:\n");  rodOFC1->InitFCVoltError3D();
    printf(" ------- OFC copper rod:\n");   rodOFC2->InitFCVoltError3D();
    printf(" ------- IFC rotated clip:\n"); rotIFC->InitFCVoltError3D();
    printf(" ------- IFC rod & strip:\n");  rodIFC1->InitFCVoltError3D();
    printf(" ------- IFC copper rod:\n");   rodIFC2->InitFCVoltError3D();
    printf(" ------- Dummy all:\n");        rodAll->InitFCVoltError3D();
    // give names
    rotOFC->SetName("rotOFC");rotOFC->SetTitle("rotOFC");
    rodOFC1->SetName("rodOFC1");rodOFC1->SetTitle("rodOFC1");
    rodOFC2->SetName("rodOFC2");rodOFC2->SetTitle("rodOFC2");
    rotIFC->SetName("rotIFC");rotIFC->SetTitle("rotIFC");
    rodIFC1->SetName("rodIFC1");rodIFC1->SetTitle("rodIFC1");
    rodIFC2->SetName("rodIFC2");rodIFC2->SetTitle("rodIFC2");
    rodAll->SetName("rodAll");rodAll->SetTitle("rodAll");
    //
    // save in file
    AliTPCComposedCorrection *corFit0 = new AliTPCComposedCorrection();
    TObjArray *cs = new TObjArray();
    cs->Add(rotIFC); cs->Add(rotOFC);
    cs->Add(rodIFC1); cs->Add(rodOFC1);
    cs->Add(rodIFC2); cs->Add(rodOFC2);
    cs->Add(rodAll);
    corFit0->SetCorrections(cs);
    corFit0->SetOmegaTauT1T2(0,1.,1.);
    corFit0->Print();    
    corFit0->StoreInOCDB(0,1000000,"testForFit");    
  }
  //
  AliTPCCorrection::AddVisualCorrection(rotOFC,0); 
  AliTPCCorrection::AddVisualCorrection(rodOFC1,1); 
  AliTPCCorrection::AddVisualCorrection(rodOFC2,2); 
  AliTPCCorrection::AddVisualCorrection(rotIFC,3); 
  AliTPCCorrection::AddVisualCorrection(rodIFC1,4); 
  AliTPCCorrection::AddVisualCorrection(rodIFC2,5); 
  AliTPCCorrection::AddVisualCorrection(rodAll,6); 
}


void RegisterAlignFunction(){

  TGeoHMatrix matrixX;
  TGeoHMatrix matrixY;
  TGeoHMatrix matrixZ;
  TGeoRotation rot0;
  TGeoRotation rot1;
  TGeoRotation rot2;  //transformation matrices
  TGeoRotation rot90;  //transformation matrices
  matrixX.SetDx(0.1); matrixY.SetDy(0.1); matrixZ.SetDz(0.1); //1 mm translation
  rot0.SetAngles(0.001*TMath::RadToDeg(),0,0);
  rot1.SetAngles(0,0.001*TMath::RadToDeg(),0);
  rot2.SetAngles(0,0,0.001*TMath::RadToDeg());
  //how to get rot02 ?
  rot90.SetAngles(0,90,0);
  rot2.MultiplyBy(&rot90,kTRUE);
  rot90.SetAngles(0,-90,0);
  rot2.MultiplyBy(&rot90,kFALSE);
  AliTPCCalibGlobalMisalignment *alignRot0  =new  AliTPCCalibGlobalMisalignment;
  alignRot0->SetAlignGlobal(&rot0);
  AliTPCCalibGlobalMisalignment *alignRot1=new  AliTPCCalibGlobalMisalignment;
  alignRot1->SetAlignGlobal(&rot1);
  AliTPCCalibGlobalMisalignment *alignRot2=new  AliTPCCalibGlobalMisalignment;
  alignRot2->SetAlignGlobal(&rot2);
  //
  AliTPCCalibGlobalMisalignment *alignTrans0  =new  AliTPCCalibGlobalMisalignment;
  alignTrans0->SetAlignGlobal(&matrixX);
  AliTPCCalibGlobalMisalignment *alignTrans1=new  AliTPCCalibGlobalMisalignment;
  alignTrans1->SetAlignGlobal(&matrixY);
  AliTPCCalibGlobalMisalignment *alignTrans2=new  AliTPCCalibGlobalMisalignment;
  alignTrans2->SetAlignGlobal(&matrixZ);
  AliTPCCorrection::AddVisualCorrection(alignTrans0  ,10);
  AliTPCCorrection::AddVisualCorrection(alignTrans1  ,11);
  AliTPCCorrection::AddVisualCorrection(alignTrans2  ,12);
  AliTPCCorrection::AddVisualCorrection(alignRot0    ,13);
  AliTPCCorrection::AddVisualCorrection(alignRot1    ,14);
  AliTPCCorrection::AddVisualCorrection(alignRot2    ,15);
  TObjArray arrAlign(6);
  arrAlign.AddAt(alignTrans0->Clone(),0);
  arrAlign.AddAt(alignTrans1->Clone(),1);
  arrAlign.AddAt(alignTrans2->Clone(),2);
  arrAlign.AddAt(alignRot0->Clone(),3);
  arrAlign.AddAt(alignRot1->Clone(),4);
  arrAlign.AddAt(alignRot2->Clone(),5);
  //combAlign.SetCorrections((TObjArray*)arrAlign.Clone());
}


void LoadTrees(){

  TFile *fy = new TFile("clusterDY.root");
  TFile *fz = new TFile("clusterDZ.root");
  treeDY= (TTree*)fy->Get("delta");
  treeDZ= (TTree*)fz->Get("delta");
  TCut cutAll = "entries>3000&&abs(kZ)<1&&localX>80&&localX<246&&abs(sector-int(sector)-0.5)<0.41&&abs(localX-134)>2&&rmsG<0.5&&abs(localX-189)>3";

  treeDY->Draw(">>outListY",cutAll,"entryList");
  TEntryList *elistOutY = (TEntryList*)gDirectory->Get("outListY");
  treeDY->SetEntryList(elistOutY);
  treeDZ->Draw(">>outListZ",cutAll,"entryList");
  TEntryList *elistOutZ = (TEntryList*)gDirectory->Get("outListZ");
  treeDZ->SetEntryList(elistOutZ);
  //
  // Make aliases
  //
  treeDY->SetAlias("dsec","(sector-int(sector)-0.5)");
  treeDY->SetAlias("dsec0","(sector-int(sector))");
  treeDY->SetAlias("signy","(-1.+2*(sector-int(sector)>0.5))");
  treeDY->SetAlias("dx","(localX-165.5)");   // distance X to ref. plane
  treeDY->SetAlias("dxm","0.001*(localX-165.5)");
  treeDY->SetAlias("rx","(localX/166.5)");   //ratio distrnace to reference plane
  treeDY->SetAlias("dq1","(((q1==0)*(-rx)+q1*(1-rx))*signy)");
  //
  {for (Int_t isec=0; isec<18; isec++){ // sectors
      treeDY->SetAlias(Form("sec%d",isec), Form("(abs(sector-%3.1lf)<0.5)",isec+0.5));
      treeDZ->SetAlias(Form("sec%d",isec), Form("(abs(sector-%3.1lf)<0.5)",isec+0.5));
    }}
  treeDY->SetAlias("gy","localX*sin(pi*sector/9)");
  treeDY->SetAlias("gx","localX*cos(pi*sector/9)");
  treeDY->SetAlias("side","(-1.+2.*(kZ>0))");  
  treeDY->SetAlias("drphi","mean");
  treeDY->SetMarkerStyle(25);
}

Double_t DeltaLookup(Double_t sector, Double_t localX, Double_t kZ, Double_t xref, Int_t value, Int_t corr){

  Double_t distortion    = AliTPCCorrection::GetCorrSector(sector,localX,kZ,value,corr);
  Double_t distortionRef = AliTPCCorrection::GetCorrSector(sector,xref,kZ,value,corr)*localX/xref;
  return distortion-distortionRef;
}

void MakeAliases(){

  AliTPCROC * roc = AliTPCROC::Instance();
  Double_t xref  = ( roc->GetPadRowRadii(0,0)+roc->GetPadRowRadii(36,roc->GetNRows(36)-1))*0.5;
  treeDY->SetAlias("iroc","(localX<134)");  // IROC
  treeDY->SetAlias("oroc","(localX>134)");  // OROC
  treeDY->SetAlias("q1","abs(localX-161)<28"); // OROC first haf
  treeDY->SetAlias("q2","(localX>189)");       // OROC second half
  //
  treeDY->SetAlias("dIFC","abs(localX-83)");    // distance to IFC
  treeDY->SetAlias("dOFC","abs(localX-250)");   // distance to OFC
  treeDY->SetAlias("errY","(1.+1/(1+(dIFC/2.))+1/(1+dOFC/2.))");

  //
  // rotated clip - IFC --------------------
  treeDY->SetAlias("rotClipI","DeltaLookup(sector,localX,kZ,165.6,1,3+0)");
  // rotated clip - OFC --------------------
  treeDY->SetAlias("rotClipO","DeltaLookup(sector,localX,kZ,165.6,1,0+0)");
  for (Int_t isec=0;isec<18; isec++){    
    //
    // alignment IROC - shift cm - rotation mrad
    treeDY->SetAlias(Form("dyI_%d",isec),Form("(iroc*sec%d)",isec));
    treeDY->SetAlias(Form("drotI_%d",isec),Form("(iroc*sec%d*(localX-%3.1lf)/1000.)",isec,xref));
    // alignment OROC
    treeDY->SetAlias(Form("dyO_%d",isec),Form("(sec%d*oroc-sec%d*localX/165.6)",isec,isec)); // fix local Y shift - do not use it
    treeDY->SetAlias(Form("drotO_%d",isec),Form("(sec%d*oroc*(localX-%3.1lf)/1000.)",isec,xref));
    //
    // quadrant shift OROC 
    treeDY->SetAlias(Form("dqO0_%d",isec),Form("(sec%d*sign(dsec)*(q1-localX/165.6))",isec));
    // quadrant delta rotation OFC inner
    treeDY->SetAlias(Form("dqOR0_%d",isec),Form("(sec%d*sign(dsec)*q1*(localX-165.6)/1000.)",isec));
    //
    treeDY->SetAlias(Form("dqO1_%d",isec),Form("(q2*sec%d*sign(dsec))",isec));  // delta
    treeDY->SetAlias(Form("dqOR1_%d",isec),Form("(q2*sec%d*sign(dsec)*(localX-165.6)/1000.)",isec));
    treeDY->SetAlias(Form("dqO2_%d",isec),Form("(q2*sec%d)",isec)); //common
    treeDY->SetAlias(Form("dqOR2_%d",isec),Form("(q2*sec%d*(localX-165.6)/1000.)",isec));
    //
    //
    // shifted rod -  OFC
    treeDY->SetAlias(Form("rodStripO_%d",isec),Form("DeltaLookup(sector-%d,localX,kZ,165.6,1,1+0)",isec));
    // Shifted cooper rod OFC 
    treeDY->SetAlias(Form("coppRodO_%d",isec),Form("DeltaLookup(sector-%d,localX,kZ,165.6,1,2+0)",isec)); 
    // shifted rod - IFC 
    treeDY->SetAlias(Form("rodStripI_%d",isec),Form("DeltaLookup(sector-%d,localX,kZ,165.6,1,4+0)",isec)); 
    // Shifted cooper rod IFC 
    treeDY->SetAlias(Form("coppRodI_%d",isec),Form("DeltaLookup(sector-%d,localX,kZ,165.6,1,5+0)",isec)); 
  }
  //
  // fitted correction - Using the saved coposed corrections
  //
  treeDY->SetAlias("dAlign","DeltaLookup(sector,localX,kZ,165.6,1,1000+0)");
  treeDY->SetAlias("dQuadrant","DeltaLookup(sector,localX,kZ,165.6,1,1100+0)");
  treeDY->SetAlias("dField","DeltaLookup(sector,localX,kZ,165.6,1,100+0)");
  treeDY->SetAlias("dAll","(dAlign+dQuadrant+dField)");
}

AliTPCCalibGlobalMisalignment *MakeAlignCorrection(TVectorD paramA, TVectorD paramC,  TMatrixD covar, Double_t chi2){

  AliTPCCalibGlobalMisalignment *alignLocal  =new  AliTPCCalibGlobalMisalignment;  
  Int_t offset=3;
  AliTPCROC * roc = AliTPCROC::Instance();
  Double_t xref  = ( roc->GetPadRowRadii(0,0)+roc->GetPadRowRadii(36,roc->GetNRows(36)-1))*0.5;
  // 
  pcWorkspace->GetFile()->cd();
  TObjArray * array = new TObjArray(72);
  
  for (Int_t side=-1; side<2; side+=2){
    TVectorD &param = (side==1) ? paramA : paramC;
    TGeoHMatrix matrixGX;
    TGeoHMatrix matrixGY;
    matrixGX.SetDx(0.1*param[1]); 
    matrixGY.SetDy(0.1*param[2]); 
    //  
    for (Int_t isec=0; isec<18; isec++){ 
      TGeoHMatrix matrixOROC;
      TGeoHMatrix matrixIROC; 
      TGeoRotation rotIROC;
      TGeoRotation rotOROC;
      Double_t phi= (Double_t(isec)+0.5)*TMath::Pi()/9.;
      //
      Double_t drotIROC = -param[offset+isec+18]*0.001;
      Double_t drotOROC = -param[offset+isec+36]*0.001;
      Double_t dlyIROC  = param[offset+isec];
      Double_t dlyIROC0 = param[offset+isec]+drotIROC*xref;     // shift at x=0
      Double_t dlyOROC  = 0;
      Double_t dlyOROC0 = 0+drotOROC*xref;  // shift at x=0
      //
      Double_t dgxIROC = TMath::Cos(phi)*0+TMath::Sin(phi)*dlyIROC0;
      Double_t dgyIROC = TMath::Sin(phi)*0-TMath::Cos(phi)*dlyIROC0;
      Double_t dgxOROC = TMath::Cos(phi)*0+TMath::Sin(phi)*dlyOROC0;
      Double_t dgyOROC = TMath::Sin(phi)*0-TMath::Cos(phi)*dlyOROC0;    
      Double_t errYIROC=TMath::Sqrt(covar(offset+isec, offset+isec)*chi2);
      Double_t errPhiIROC=TMath::Sqrt(covar(offset+isec+18, offset+isec+18)*chi2);
      Double_t errPhiOROC=TMath::Sqrt(covar(offset+isec+36, offset+isec+36)*chi2);
      matrixIROC.SetDx(dgxIROC);
      matrixIROC.SetDy(dgyIROC);
      matrixOROC.SetDx(dgxOROC);
      matrixOROC.SetDy(dgyOROC);
      rotIROC.SetAngles(drotIROC*TMath::RadToDeg(),0,0);
      matrixIROC.Multiply(&matrixGX);
      matrixIROC.Multiply(&matrixGY);
      matrixIROC.Multiply(&rotIROC);
      rotOROC.SetAngles(drotOROC*TMath::RadToDeg(),0,0);
      matrixOROC.Multiply(&matrixGX);
      matrixOROC.Multiply(&matrixGY);
      matrixOROC.Multiply(&rotOROC);
      if (side>0){
  array->AddAt(matrixIROC.Clone(),isec);
  array->AddAt(matrixOROC.Clone(),isec+36);
      }
      if (side<0){
  array->AddAt(matrixIROC.Clone(),isec+18);
  array->AddAt(matrixOROC.Clone(),isec+36+18);
      }
      Double_t rms=TMath::Sqrt(chi2); //chi2 normalized 1/npoints before
      (*pcWorkspace)<<"align"<<
  "isec="<<isec<<           // sector
  "side="<<side<<           // side
  "phi="<<phi<<             // phi
  "rms="<<rms<<             // rms in cm
  "cov.="<<&covar<<
  // errors
  "errYIROC="<<errYIROC<<      //error in local y
  "errPhiIROC="<<errPhiIROC<<  //error in phi IROC
  "errPhiOROC="<<errPhiOROC<<  //error in phi OROC
  //
  "dlyIROC="<<dlyIROC<<        //delta Local Y at refX=165 -IROC
  "dlyIROC0="<<dlyIROC0<<      //delta Local Y at refX=0   -IROC
  "dgxIROC="<<dgxIROC<<        //
  "dgyIROC="<<dgyIROC<<
  "drotIROC="<<drotIROC<<
  //
  "dlyOROC="<<dlyOROC<<       // assuming 0 at 
  "dlyOROC0="<<dlyOROC0<<
  "dgxOROC="<<dgxOROC<<
  "dgyOROC="<<dgyOROC<<
  "drotOROC="<<drotOROC<<
  "\n";
    } 
  }
  alignLocal->SetAlignSectors(array);
  AliTPCCorrection::AddVisualCorrection(alignLocal,1000);
  /*
    Check:
    alignLocal->CreateHistoDRPhiinXY(10,100,100)->Draw("colz");  // OK
    treeDY->Draw("DeltaLookup(sector,localX,kZ,165.6,1,1001):localX:int(sector)","kZ>0","colz");
    //
    (( (*pcWorkspace)<<"align"))->GetTree()->Draw("dlyOROC-dlyOROC:isec","","*")
    treeDY->Draw("tfitA:DeltaLookup(sector,localX,kZ,165.6,1,1001):sector","kZ>0","colz");

  */
  return alignLocal;
}

AliTPCCalibGlobalMisalignment *MakeQuadrantCorrection(TVectorD paramA, TVectorD paramC, TMatrixD covar, Double_t chi2){

  AliTPCCalibGlobalMisalignment *alignLocalQuadrant  =new  AliTPCCalibGlobalMisalignment;  
  //
  Int_t offset=3+3*18;
  TVectorD quadrantQ0(36);   //dly+=sign*(*fQuadrantQ0)[isec];  // shift in cm
  TVectorD quadrantRQ0(36);  //dly+=sign*(*fQuadrantRQ0)[isec]*(pos[0]-xref);      
  TVectorD quadrantQ1(36);   //dly+=sign*(*fQuadrantQ1)[isec];  // shift in cm
  TVectorD quadrantRQ1(36);  //dly+=sign*(*fQuadrantRQ1)[isec]*(pos[0]-xref);      
  TVectorD quadrantQ2(36);   //dly+=(*fQuadrantQ2)[isec];  // shift in cm
  TVectorD quadrantRQ2(36);  //dly+=(*fQuadrantRQ2)[isec]*(pos[0]-xref);      

  for (Int_t side=-1; side<2; side+=2){
    Int_t offsetSec= (side==1) ? 0:18;
    TVectorD &param = (side==1) ? paramA : paramC;
    for (Int_t isec=0; isec<18; isec++){ 
      Double_t  q0=-param[offset+isec+0*18];
      Double_t rq0=-param[offset+isec+1*18]*0.001;
      Double_t  q1=-param[offset+isec+2*18];
      Double_t rq1=-param[offset+isec+3*18]*0.001;
      Double_t  q2=-param[offset+isec+4*18];
      Double_t rq2=-param[offset+isec+5*18]*0.001;
      //
      Double_t  sq0=TMath::Sqrt(covar(offset+isec+0*18,offset+isec+0*18)*chi2);
      Double_t srq0=TMath::Sqrt(covar(offset+isec+1*18,offset+isec+1*18)*chi2)*0.001;
      Double_t  sq1=TMath::Sqrt(covar(offset+isec+2*18,offset+isec+2*18)*chi2);
      Double_t srq1=TMath::Sqrt(covar(offset+isec+3*18,offset+isec+3*18)*chi2)*0.001;
      Double_t  sq2=TMath::Sqrt(covar(offset+isec+4*18,offset+isec+4*18)*chi2);
      Double_t srq2=TMath::Sqrt(covar(offset+isec+5*18,offset+isec+5*18)*chi2)*0.001;
      //
      quadrantQ0[offsetSec+isec]=q0;
      quadrantRQ0[offsetSec+isec]=rq0;
      quadrantQ1[offsetSec+isec]=q1;
      quadrantRQ1[offsetSec+isec]=rq1;
      quadrantQ2[offsetSec+isec]=q2;
      quadrantRQ2[offsetSec+isec]=rq2;
      Double_t rms=TMath::Sqrt(chi2); //chi2 normalized 1/npoints before
      (*pcWorkspace)<<"quadrant"<<
  "isec="<<isec<<           // sector
  "side="<<side<<           // side
  "rms="<<rms<<             // rms in cm
  "cov.="<<&covar<<
  //
  "q0="<<q0<<               // quadrant alignment
  "rq0="<<rq0<<
  "q1="<<q1<<
  "rq1="<<rq1<<
  "q2="<<q2<<
  "rq2="<<rq2<<
  "sq0="<<sq0<<             //rms of quadrant parameters 
  "srq0="<<srq0<<
  "sq1="<<sq1<<
  "srq1="<<srq1<<
  "sq2="<<sq2<<
  "srq2="<<srq2<<
  "\n";    
    }
  }
  alignLocalQuadrant->SetQuadranAlign(&quadrantQ0, &quadrantRQ0, &quadrantQ1, &quadrantRQ1, &quadrantQ2, &quadrantRQ2);
  AliTPCCorrection::AddVisualCorrection(alignLocalQuadrant,1100);  
  /*
    (( (*pcWorkspace)<<"quadrant"))->GetTree()->Draw("q0:isec","","*")

    alignLocalQuadrant->CreateHistoDRPhiinXY(10,100,100)->Draw("colz");  //OK
    
    treeDY->Draw("tfitA-DeltaLookup(sector,localX,kZ,165.6,1,1000):DeltaLookup(sector,localX,kZ,165.6,1,1100):int(sector)","kZ>0&&oroc","colz");
  */
  return alignLocalQuadrant;
}

AliTPCFCVoltError3D* MakeEfieldCorrection(TVectorD paramA, TVectorD paramC, TMatrixD covar, Double_t chi2){

  Int_t offset=3+9*18;
  //
  AliTPCFCVoltError3D* corrField = new AliTPCFCVoltError3D;
  //
  Double_t rotIROCA=paramA[offset+0];
  Double_t rotIROCC=paramC[offset+0];
  Double_t rotOROCA=paramA[offset+1];
  Double_t rotOROCC=paramC[offset+1];
  //
  corrField->SetRotatedClipVoltA(0,paramA[offset+0]);
  corrField->SetRotatedClipVoltC(0,paramC[offset+0]);
  corrField->SetRotatedClipVoltA(1,paramA[offset+1]);
  corrField->SetRotatedClipVoltC(1,paramC[offset+1]);
  {for (Int_t isec=0; isec<18; isec++){
      corrField->SetRodVoltShiftA(isec,paramA[offset+2+36+isec]);     // rod shift IFC
      corrField->SetRodVoltShiftA(18+isec,paramA[offset+2+isec]);     // rod shift OFC
      corrField->SetRodVoltShiftC(isec,paramC[offset+2+36+isec]);    // rod shift IFC
      corrField->SetRodVoltShiftC(18+isec,paramC[offset+2+isec]);    // rod shift OFC
      Double_t rodIROCA=paramA[offset+2+36+isec];
      Double_t rodIROCC=paramC[offset+2+36+isec];
      Double_t rodOROCA=paramA[offset+2+isec];
      Double_t rodOROCC=paramC[offset+2+isec];
      //
      Double_t srodIROCA=TMath::Sqrt(covar(offset+2+36+isec,offset+2+36+isec)*chi2);
      Double_t srodOROCA=TMath::Sqrt(covar(offset+2+isec,offset+2+isec)*chi2);
      Double_t phi= (Double_t(isec)+0.5)*TMath::Pi()/9.;
      Double_t rms=TMath::Sqrt(chi2); //chi2 normalized 1/npoints before
      (*pcWorkspace)<<"field"<<
  "isec="<<isec<<           // sector
  //"side="<<side<<           // side
  "phi="<<phi<<             // phi
  "rms="<<rms<<             // rms in cm
  "cov.="<<&covar<<
  //
  "rotIROCA="<<rotIROCA<<  
  "rotIROCC="<<rotIROCC<<
  "rotOROCA="<<rotOROCA<<
  "rotOROCC="<<rotOROCC<<
  //
  "rodIROCA="<<rodIROCA<<
  "rodIROCC="<<rodIROCC<<
  "rodOROCA="<<rodOROCA<<
  "rodOROCC="<<rodOROCC<<
  //
  "srodIROCA="<<srodIROCA<<
  "srodOROCA="<<srodOROCA<<
  "srodIROCC="<<srodIROCA<<
  "srodOROCC="<<srodOROCA<<
  "\n";
    }    
  }
  corrField->SetOmegaTauT1T2(0,1.,1.);
  corrField->InitFCVoltError3D();
  AliTPCCorrection::AddVisualCorrection(corrField,100);  
  return corrField;
}


void FitRodShift(Bool_t flagIFCcopper = kTRUE) {

  LoadTrees();
  LoadModels();
  RegisterAlignFunction();
  MakeAliases();

  pcWorkspace= new TTreeSRedirector("fitAlignLookup.root");
  TFormula::SetMaxima(10000);
  //
  // 1. fit Global
  //
  Double_t chi2G=0;   Int_t npointsG=0;  TVectorD paramG;  TMatrixD covarG;
  TString fstringGlobal="";  
  fstringGlobal+="DeltaLookup(sector,localX,kZ,165.6,1,10-0)++"; // deltaGX
  fstringGlobal+="DeltaLookup(sector,localX,kZ,165.6,1,11-0)++"; // deltaGY
  fstringGlobal+="DeltaLookup(sector,localX,kZ,165.6,1,10-0)*side++"; // deltaGX - sides
  fstringGlobal+="DeltaLookup(sector,localX,kZ,165.6,1,11-0)*side++"; // deltaGY -sides
  //
  TString *strFitGlobal = TStatToolkit::FitPlane(treeDY,"drphi", fstringGlobal.Data(),"1", chi2G,npointsG,paramG,covarG,-1,0, 10000000, kTRUE);
  treeDY->SetAlias("fitYGlobal",strFitGlobal->Data());
  strFitGlobal->Tokenize("++")->Print();
  printf("chi2=%f\n",TMath::Sqrt(chi2G/npointsG));
  // testplots - if necessary - e.g.
  // treeDY->Draw("AliTPCCorrection::GetCorrSector(sector,localX,kZ,1,2):sector:localX","Cut","colz")

  // FIT PREPARATION +++++++++++++++++++++++++++++++++
  // define the fit function

  TString  fstring="";         

  //
  // Alignment
  //
  {
    fstring+="DeltaLookup(sector,localX,kZ,165.6,1,10-0)++"; // deltaGX
    fstring+="DeltaLookup(sector,localX,kZ,165.6,1,11-0)++"; // deltaGY
    // alignment IROC
    for (Int_t i=0;i<18;i++){
      fstring+=Form("dyI_%d++",i);  // alignment - linear shift (in cm)
    }
    for (Int_t i=0;i<18;i++){
      fstring+=Form("drotI_%d++",i);  // alignment - rotation (in mrad)
    }
    // alignment OROC
    for (Int_t i=0;i<18;i++){ // shift of OROC is not allowed
      //fstring+=Form("dyO_%d++",i);  // alignment - linear shift (in cm)
    }
    for (Int_t i=0;i<18;i++){
      fstring+=Form("drotO_%d++",i);  // alignment - rotation (in mrad)
    }
    //
    for (Int_t i=0;i<18;i++){
      fstring+=Form("dqO0_%d++",i);  // alignment - quadrant shift OROC in
    }
    for (Int_t i=0;i<18;i++){
      fstring+=Form("dqOR0_%d++",i);  // alignment - quadrant rotation OROC in
    }
    for (Int_t i=0;i<18;i++){
      fstring+=Form("dqO1_%d++",i);  // alignment - quadrant shift OROC out
    }
    for (Int_t i=0;i<18;i++){
      fstring+=Form("dqOR1_%d++",i);  // alignment - quadrant rotation OROC out
    }
    for (Int_t i=0;i<18;i++){
      fstring+=Form("dqO2_%d++",i);  // alignment - quadrant shift OROC out
    }
    for (Int_t i=0;i<18;i++){
      fstring+=Form("dqOR2_%d++",i);  // alignment - quadrant rotation OROC out
    }
  }
  //
  // Field distortion
  //
  {
    // rotated clip - IFC --------------------
    treeDY->SetAlias("rotClipI","DeltaLookup(sector,localX,kZ,165.6,1,3+0)");
    fstring+="rotClipI++";
    // rotated clip - OFC --------------------
    treeDY->SetAlias("rotClipO","DeltaLookup(sector,localX,kZ,165.6,1,0+0)");
    fstring+="rotClipO++";
    
    // shifted rod -  OFC
    for (Int_t i=0;i<18;i++){
      fstring+=Form("rodStripO_%d++",i);
    }    
    // Shifted cooper rod OFC     
    for (Int_t i=0;i<18;i++){
      fstring+=Form("coppRodO_%d++",i);
    }    
    // shifted rod - IFC 
    for (Int_t i=0;i<18;i++){
      fstring+=Form("rodStripI_%d++",i);
    }    
    // Shifted cooper rod IFC 
    if (flagIFCcopper) for (Int_t i=0;i<18;i++){
  fstring+=Form("coppRodI_%d++",i);
    }  
  }
  //fstring+=fstringGlobal;
 
  // FIT ++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  Double_t chi2A=0;   Int_t npointsA=0;  TVectorD paramA;  TMatrixD covarA;
  Double_t chi2C=0;   Int_t npointsC=0;  TVectorD paramC;  TMatrixD covarC;
  
  
  printf("Fitting A side\n");
  TCut cutAllA = "Cut&&kZ>0";
  TString *strFitGA = TStatToolkit::FitPlane(treeDY,"drphi:errY", fstring.Data(),"kZ>0", chi2A,npointsA,paramA,covarA,-1,0, 10000000, kTRUE);
  treeDY->SetAlias("tfitA",strFitGA->Data());
  // strFitGA->Tokenize("++")->Print();
  printf("chi2=%f\n",TMath::Sqrt(chi2A/npointsA));
  printf("Sigma Y:\t%f (mm)\n",10.*TMath::Sqrt(covarA(3,3)*chi2A/npointsA));
  printf("IROC Sigma Angle:\t%f (mrad)\n",TMath::Sqrt(covarA(3+18,3+18)*chi2A/npointsA));
  printf("OROC Sigma Angle:\t%f (mrad)\n",TMath::Sqrt(covarA(3+36,3+36)*chi2A/npointsA));

  printf("Fitting C side\n");
  TCut cutAllC = "Cut&&kZ<0";
  TString *strFitGC = TStatToolkit::FitPlane(treeDY,"drphi:errY", fstring.Data(),"kZ<0", chi2C,npointsC,paramC,covarC,-1,0, 10000000, kTRUE);
  treeDY->SetAlias("tfitC",strFitGC->Data());
  //  strFitGC->Tokenize("++")->Print();
  printf("chi2=%f\n",TMath::Sqrt(chi2C/npointsC));
  printf("Sigma Y:\t%f (mm)\n",10.*TMath::Sqrt(covarC(3,3)*chi2C/npointsC));
  printf("IROC Sigma Angle:\t%f (mrad)\n",TMath::Sqrt(covarC(3+18,3+18)*chi2C/npointsC));
  printf("OROC Sigma Angle:\t%f (mrad)\n",TMath::Sqrt(covarC(3+36,3+36)*chi2C/npointsC));

  //
  // make combined correction
  //
  TH2 *hdist =0;
  //
  AliTPCCalibGlobalMisalignment * alignLocal = MakeAlignCorrection(paramA, paramC, covarA, chi2A/npointsA);
  alignLocal->SetName("alignLocal");
  alignLocal->Write("alignLocal");
  hdist = alignLocal->CreateHistoDRPhiinXY(10,250,250);
  hdist->SetName("AlignLocalAside"); hdist->SetTitle("Alignment map (A side)");
  hdist->Write("AlignLocalAside");
  hdist = alignLocal->CreateHistoDRPhiinXY(-10,250,250);
  hdist->SetName("AlignLocalCside"); hdist->SetTitle("Alignment map (C side)");
  hdist->Write("AlignLocalCside");
  //
  AliTPCCalibGlobalMisalignment *alignQuadrant = MakeQuadrantCorrection(paramA, paramC,  covarA, chi2A/npointsA);
  alignQuadrant->SetName("alignQuadrant");
  alignQuadrant->Write("alignQuadrant");
  hdist = alignQuadrant->CreateHistoDRPhiinXY(10,250,250);
  hdist->SetName("AlignQuadrantAside"); hdist->SetTitle("Quadrant Alignment map (A side)");
  hdist->Write("AlignQuadrantAside");
  hdist = alignQuadrant->CreateHistoDRPhiinXY(-10,250,250);
  hdist->SetName("AlignQuadrantCside"); hdist->SetTitle("Quadrant Alignment map (C side)");
  hdist->Write("AlignQuadrantCside");
  //
  AliTPCFCVoltError3D* corrField = MakeEfieldCorrection(paramA, paramC, covarA, chi2A/npointsA);
  corrField->SetName("corrField");
  corrField->Write("corrField");

  hdist = corrField->CreateHistoDRPhiinXY(10,250,250);
  hdist->SetName("AlignEfieldAside"); hdist->SetTitle("Efield Alignment map (A side)");
  hdist->Write("AlignEfieldAside");
  hdist = corrField->CreateHistoDRPhiinXY(-10,250,250);
  hdist->SetName("AlignEfieldCside"); hdist->SetTitle("Efield Alignment map (C side)");
  hdist->Write("AlignEfieldCside");


  //
  delete pcWorkspace;
  MakeQA();
  return;  
}


void MakeQA(){
  LoadTrees();
  LoadModels();
  RegisterAlignFunction();
  MakeAliases();
  TFile f("fitAlignLookup.root","update");  
  AliTPCCalibGlobalMisalignment*     alignLocal =  (AliTPCCalibGlobalMisalignment*)f.Get("alignLocal");
  AliTPCCalibGlobalMisalignment*  alignQuadrant =  (AliTPCCalibGlobalMisalignment*)f.Get("alignQuadrant");
  AliTPCFCVoltError3D      *corrField= (AliTPCFCVoltError3D*)f.Get("corrField");
  AliTPCCorrection::AddVisualCorrection(alignLocal,1000);
  AliTPCCorrection::AddVisualCorrection(alignQuadrant,1100);
  AliTPCCorrection::AddVisualCorrection(corrField,100);
  FitFunctionQA();
  DrawAlignParam();
  f.Close();
}

void FitFunctionQA(){

  TH1 *his=0;
  TCanvas *canvasDist= new TCanvas("FitQA","fitQA",1200,800);
  canvasDist->Divide(2,2);
  //
  canvasDist->cd(1)->SetLogy(kFALSE); canvasDist->cd(1)->SetRightMargin(0.15);
  treeDY->Draw("10*mean:sector:localX","kZ<0&&localX<134","colz"); 
  his= (TH1*)treeDY->GetHistogram()->Clone();   his->SetName("DeltaRPhi1");
  his->GetXaxis()->SetTitle("Sector"); 
  his->GetZaxis()->SetTitle("R (cm)"); 
  his->GetYaxis()->SetTitle("#Delta_{r#phi} (mm)");
  his->Draw("cont2"); 
  his->Draw("colz");
  //
  canvasDist->cd(2)->SetLogy(kFALSE); canvasDist->cd(2)->SetRightMargin(0.15);
  treeDY->Draw("10*mean:sector:localX","kZ<0&&localX>160","colz"); 
  his= (TH1*)treeDY->GetHistogram()->Clone();   his->SetName("DeltaRPhi2");
  his->GetXaxis()->SetTitle("Sector"); 
  his->GetZaxis()->SetTitle("R (cm)"); 
  his->GetYaxis()->SetTitle("#Delta_{r#phi} (mm)");
  his->Draw("cont2"); 
  his->Draw("colz");
  //
  canvasDist->cd(3)->SetLogy(kFALSE); canvasDist->cd(3)->SetRightMargin(0.15);
  treeDY->Draw("10*(mean-dAll):sector:localX","kZ<0&&localX<134","colz"); 
  his= (TH1*)treeDY->GetHistogram()->Clone();   his->SetName("DeltaRPhi3"); his->SetTitle("Delta #RPhi");
  his->GetXaxis()->SetTitle("Sector"); 
  his->GetZaxis()->SetTitle("R (cm)"); 
  his->GetYaxis()->SetTitle("#Delta_{r#phifit}-#Delta_{r#phi} (mm)"); 
  his->Draw("cont2"); 
  his->Draw("colz");
  //
  canvasDist->cd(4)->SetLogy(kFALSE); canvasDist->cd(4)->SetRightMargin(0.15);
  treeDY->SetMarkerColor(1);
  treeDY->Draw("10*mean:10*dAll:localX","","colz"); 
  his= (TH1*)treeDY->GetHistogram()->Clone();   his->SetName("DeltaRPhi4");
  his->GetXaxis()->SetTitle("Fit value #Delta_{r#phi} (mm)"); 
  his->GetYaxis()->SetTitle("#Delta_{r#phi} (mm)"); 
  his->GetZaxis()->SetTitle("R (cm)"); 
  his->Draw("cont2");
  his->Draw("colz");
  //
  canvasDist->Write("FitQA");
}

void DrawAlignParam(){

  TFile f("fitAlignLookup.root");
  TTree * treeAlign=(TTree*)f.Get("align");
  TTree * treeQuadrant=(TTree*)f.Get("quadrant");
  TCanvas *canvasAlign=new TCanvas("align","align",1000,800); 
  TH1 *his=0;
  canvasAlign->Divide(2,2);
  treeAlign->SetMarkerStyle(25);
  treeQuadrant->SetMarkerStyle(25);
  gStyle->SetOptStat(kTRUE);
  //
  canvasAlign->cd(1); canvasAlign->cd(1)->SetRightMargin(0.15);
  treeAlign->Draw("10*dlyIROC*side:isec:1+side","","colz");
  his= (TH1*)treeAlign->GetHistogram()->Clone();   his->SetName("dlyIROC");his->SetTitle("IROC Alignment #Delta_{r#phi}");
  his->GetXaxis()->SetTitle("sector"); 
  his->GetYaxis()->SetTitle("#Delta_{r#phi} (mm)"); 
  his->GetZaxis()->SetTitle("side");
  his->Draw("cont2"); 
  his->Draw("colz");
  //
  canvasAlign->cd(2); canvasAlign->cd(2)->SetRightMargin(0.15);
  treeAlign->Draw("1000*drotIROC*side:isec:1+side","","colz");
  his= (TH1*)treeAlign->GetHistogram()->Clone();   his->SetName("drotIROC");his->SetTitle("IROC Angular Alignment #Delta_{r#phi}");
  his->GetXaxis()->SetTitle("sector"); 
  his->GetYaxis()->SetTitle("#Delta_{r#phi} (mrad)"); 
  his->GetZaxis()->SetTitle("side");
  his->Draw("cont2");
  his->Draw("colz");
  //
  canvasAlign->cd(4);canvasAlign->cd(4)->SetRightMargin(0.15);
  treeAlign->Draw("1000*drotOROC*side:isec:1+side","","colz");
  his= (TH1*)treeAlign->GetHistogram()->Clone();   his->SetName("drotOROC");his->SetTitle("OROC Angular Alignment #Delta_{r#phi}");
  his->GetXaxis()->SetTitle("sector"); 
  his->GetYaxis()->SetTitle("#Delta_{r#phi} (mrad)"); 
  his->GetZaxis()->SetTitle("side");
  his->Draw("cont2");
  his->Draw("colz");
  //
  canvasAlign->cd(3);canvasAlign->cd(3)->SetRightMargin(0.15);
  treeQuadrant->Draw("10*q2*side:isec:1+side","","colz");
  his= (TH1*)treeQuadrant->GetHistogram()->Clone();   his->SetName("drphiOROC");his->SetTitle("OROC Alignment  Outer Quadrant #Delta_{r#phi}");
  his->GetXaxis()->SetTitle("sector"); 
  his->GetYaxis()->SetTitle("#Delta_{r#phi} (mm)"); 
  his->GetZaxis()->SetTitle("side");
  his->Draw("cont2");
  his->Draw("colz");


}