]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TPC/AliTPCPreprocessorOffline.cxx
Coverity fix.
[u/mrichter/AliRoot.git] / TPC / AliTPCPreprocessorOffline.cxx
index f9b8863e37ac3cb6e34911a4ed4036429b6ef18b..0e4b5335144a53f5ef45e189c8a68daa3877fd5d 100644 (file)
 #include "AliTPCcalibTimeGain.h"
 #include "AliTPCcalibGainMult.h"
 #include "AliSplineFit.h"
+#include "AliTPCComposedCorrection.h"
+#include "AliTPCExBTwist.h"
+#include "AliTPCCalibGlobalMisalignment.h"
+#include "TStatToolkit.h"
+#include "TChain.h"
+#include "TCut.h"
+#include "AliTrackerBase.h"
 #include "AliTPCPreprocessorOffline.h"
 
 
@@ -100,6 +107,7 @@ AliTPCPreprocessorOffline::AliTPCPreprocessorOffline():
   fGainMIP(0),          // calibration component for MIP
   fGainCosmic(0),       // calibration component for cosmic
   fGainMult(0),
+  fAlignTree(0),        // alignment tree
   fSwitchOnValidation(kFALSE) // flag to switch on validation of OCDB parameters
 {
   //
@@ -207,7 +215,15 @@ void AliTPCPreprocessorOffline::CalibTimeVdrift(const Char_t* file, Int_t ustart
     Printf("TPC time drift OCDB parameters out of range!");
     return;
   }
-
+  //
+  //4.b make alignment
+  //
+  MakeFitTime();
+  TFile * ftime= TFile::Open("fitITSVertex.root");
+  if (ftime){
+    TObject * alignmentTime=ftime->Get("FitCorrectionTime");
+    if (alignmentTime) fVdriftArray->AddLast(alignmentTime);
+  }
   //
   //
   // 5. update of OCDB
@@ -262,12 +278,18 @@ Bool_t AliTPCPreprocessorOffline::ValidateTimeDrift(Double_t maxVDriftCorr)
   Printf("ValidateTimeDrift..." );
 
   TGraphErrors* gr = (TGraphErrors*)fVdriftArray->FindObject("ALIGN_ITSB_TPC_DRIFTVD");
+  Printf("ALIGN_ITSB_TPC_DRIFTVD graph = %p",gr);
+
   if(!gr) return kFALSE;
-  if(gr->GetN()<1) return kFALSE;
+  if(gr->GetN()<1)  { 
+    Printf("ALIGN_ITSB_TPC_DRIFTVD number of points = %d",gr->GetN());
+    return kFALSE;
+  }
 
   // check whether drift velocity corrections in the range
   for(Int_t iPoint = 0; iPoint<gr->GetN(); iPoint++) 
   {
+    Printf("Y value from the graph: %f",TMath::Abs(gr->GetY()[iPoint]));
     if(TMath::Abs(gr->GetY()[iPoint]) > maxVDriftCorr)  
       return kFALSE;
   }
@@ -363,7 +385,13 @@ TGraphErrors* AliTPCPreprocessorOffline::FilterGraphMedianAbs(TGraphErrors * gra
   Double_t *erry=new Double_t[npoints0];
   //
   //
-  if (npoints0<kMinPoints) return 0;
+  if (npoints0<kMinPoints) {
+    delete []outx;
+    delete []outy;
+    delete []errx;
+    delete []erry;
+    return 0;
+  }
   for (Int_t iter=0; iter<3; iter++){
     npoints=0;
     for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
@@ -410,6 +438,10 @@ void AliTPCPreprocessorOffline::AddHistoGraphs(  TObjArray * vdriftArray, AliTPC
       THnSparse* newHist=hist->Projection(4,dim);
       newHist->SetName(name);
       TGraphErrors* graph=AliTPCcalibBase::FitSlices(newHist,2,0,400,100,0.05,0.95, kTRUE);
+      if (!graph) {
+       printf("Graph =%s filtered out\n", name.Data());
+       continue;
+      }
       printf("name=%s graph=%i, N=%i\n", name.Data(), graph==0, graph->GetN());
       Int_t pos=name.Index("_");
       name=name(pos,name.Capacity()-pos);
@@ -418,10 +450,6 @@ void AliTPCPreprocessorOffline::AddHistoGraphs(  TObjArray * vdriftArray, AliTPC
       graphName.ToUpper();
       //
       graph = FilterGraphDrift(graph, kErrSigmaCut, kMedianCutAbs);
-      if (!graph) {
-       printf("Graph =%s filtered out\n", name.Data());
-       continue;
-      }
       //
       if (graph){
         graph->SetMarkerStyle(i%8+20);
@@ -528,9 +556,9 @@ void AliTPCPreprocessorOffline::AddLaserGraphs(  TObjArray * vdriftArray, AliTPC
   // add graphs for laser
   //
   const Double_t delayL0L1 = 0.071;  //this is hack for 1/2 weeks
-  THnSparse *hisN=0;
+  //THnSparse *hisN=0;
   TGraphErrors *grLaser[6]={0,0,0,0,0,0};
-  hisN = timeDrift->GetHistVdriftLaserA(0);
+  //hisN = timeDrift->GetHistVdriftLaserA(0);
   if (timeDrift->GetHistVdriftLaserA(0)){
     grLaser[0]=MakeGraphFilter0(timeDrift->GetHistVdriftLaserA(0),0,2,5,delayL0L1);
     grLaser[0]->SetName("GRAPH_MEAN_DELAY_LASER_ALL_A");
@@ -833,6 +861,10 @@ void AliTPCPreprocessorOffline::ReadGainGlobal(const Char_t* fileName){
     fGainCosmic = ( AliTPCcalibTimeGain *)fcalib.Get("calibTimeGainCosmic");
     fGainMult   = ( AliTPCcalibGainMult *)fcalib.Get("calibGainMult");
   }
+  if (!fGainMult){
+    TFile fcalibMult("TPCMultObjects.root");
+    fGainMult   = ( AliTPCcalibGainMult *)fcalibMult.Get("calibGainMult");
+  }
   TH1 * hisT=0;
   Int_t firstBinA =0, lastBinA=0;
 
@@ -1042,8 +1074,18 @@ Bool_t AliTPCPreprocessorOffline::AnalyzeGainMultiplicity() {
   TH1D * meanMax = (TH1D*)arrMax.At(1);
   TH1D * meanTot = (TH1D*)arrTot.At(1);
   Float_t meanMult = histMultMax->GetMean();
-  meanMax->Scale(1./meanMax->GetBinContent(meanMax->FindBin(meanMult)));
-  meanTot->Scale(1./meanTot->GetBinContent(meanTot->FindBin(meanMult)));
+  if(meanMax->GetBinContent(meanMax->FindBin(meanMult))) {
+    meanMax->Scale(1./meanMax->GetBinContent(meanMax->FindBin(meanMult)));
+  }
+  else {
+   return kFALSE;
+  }
+  if(meanTot->GetBinContent(meanTot->FindBin(meanMult))) {
+    meanTot->Scale(1./meanTot->GetBinContent(meanTot->FindBin(meanMult)));
+  }
+  else {
+   return kFALSE;
+  }
   Float_t xMultMax[50];
   Float_t yMultMax[50];
   Float_t yMultErrMax[50];
@@ -1054,7 +1096,8 @@ Bool_t AliTPCPreprocessorOffline::AnalyzeGainMultiplicity() {
   Int_t nCountMax = 0;
   for(Int_t iBin = 1; iBin < meanMax->GetXaxis()->GetNbins(); iBin++) {
     Float_t yValMax = meanMax->GetBinContent(iBin);
-    if (yValMax < 0.01) continue;
+    if (yValMax < 0.7) continue;
+    if (yValMax > 1.3) continue;
     if (meanMax->GetBinError(iBin)/yValMax > 0.01) continue;
     xMultMax[nCountMax] = meanMax->GetXaxis()->GetBinCenter(iBin);
     yMultMax[nCountMax] = yValMax;
@@ -1069,7 +1112,8 @@ Bool_t AliTPCPreprocessorOffline::AnalyzeGainMultiplicity() {
   Int_t nCountTot = 0;
   for(Int_t iBin = 1; iBin < meanTot->GetXaxis()->GetNbins(); iBin++) {
     Float_t yValTot = meanTot->GetBinContent(iBin);
-    if (yValTot < 0.1) continue;
+    if (yValTot < 0.7) continue;
+    if (yValTot > 1.3) continue;
     if (meanTot->GetBinError(iBin)/yValTot > 0.1) continue;
     xMultTot[nCountTot] = meanTot->GetXaxis()->GetBinCenter(iBin);
     yMultTot[nCountTot] = yValTot;
@@ -1131,9 +1175,11 @@ void AliTPCPreprocessorOffline::MakeQAPlot(Float_t  FPtoMIPratio) {
        fGraphCosmic->GetY()[i] *= FPtoMIPratio;        
       }
     }
-    fGraphCosmic->Draw("lp");
-    grfFitCosmic->SetLineColor(2);
-    grfFitCosmic->Draw("lu");
+    fGraphCosmic->Draw("lp"); 
+    if (grfFitCosmic) {
+      grfFitCosmic->SetLineColor(2);
+      grfFitCosmic->Draw("lu");
+    }
     fGainArray->AddLast(gainHistoCosmic);
     //fGainArray->AddLast(canvasCosmic->Clone());
     delete canvasCosmic;    
@@ -1158,5 +1204,357 @@ void AliTPCPreprocessorOffline::MakeQAPlot(Float_t  FPtoMIPratio) {
   }  
 }
 
+void AliTPCPreprocessorOffline::MakeFitTime(){
+  //
+  // mak aligment fit - store results in the file
+  //
+  const Int_t kMinEntries=1000;
+  MakeChainTime();
+  MakePrimitivesTime();
+  if (!fAlignTree) return;
+  if (fAlignTree->GetEntries()<kMinEntries) return;
+  fAlignTree->SetAlias("ptype","type");
+  fAlignTree->SetAlias("hasITS","(1+0)");
+  fAlignTree->SetAlias("dITS","1-2*(refX<40)");
+  fAlignTree->SetAlias("isITS","refX>10");
+  fAlignTree->SetAlias("isVertex","refX<10");
+  // 
+  Int_t  npointsMax=30000000;
+  TStatToolkit toolkit;
+  Double_t chi2=0;
+  Int_t    npoints=0;
+  TVectorD param;
+  TMatrixD covar;
+
+  TString fstringFast="";
+  fstringFast+="FExBTwistX++";
+  fstringFast+="FExBTwistY++";
+  fstringFast+="FAlignRot0D++";
+  fstringFast+="FAlignTrans0D++";
+  fstringFast+="FAlignTrans1D++";
+  //
+  fstringFast+="hasITS*FAlignTrans0++";
+  fstringFast+="hasITS*FAlignTrans1++";
+  fstringFast+="hasITS*FAlignRot0++";
+  fstringFast+="hasITS*FAlignRot1++";
+  fstringFast+="hasITS*FAlignRot2++";
+  //
+  fstringFast+="dITS*FAlignTrans0++";
+  fstringFast+="dITS*FAlignTrans1++";
+  fstringFast+="dITS*FAlignRot0++";
+  fstringFast+="dITS*FAlignRot1++";
+  fstringFast+="dITS*FAlignRot2++";
+  
+  TCut cutFit="entries>10&&abs(mean)>0.00001";
+  fAlignTree->SetAlias("err","rms");
+
+  TString *strDeltaITS = TStatToolkit::FitPlaneConstrain(fAlignTree,"mean:err", fstringFast.Data(),cutFit, chi2,npoints,param,covar,-1,0, npointsMax, 1);
+  strDeltaITS->Tokenize("++")->Print();
+  fAlignTree->SetAlias("fitYFast",strDeltaITS->Data());
+  // 
+  TVectorD paramC= param;
+  TMatrixD covarC= covar;
+  TStatToolkit::Constrain1D(fstringFast,"Trans0D",paramC,covarC,0, 0.1);
+  TStatToolkit::Constrain1D(fstringFast,"Trans1D",paramC,covarC,0, 0.1);
+  TStatToolkit::Constrain1D(fstringFast,"TwistX",paramC,covarC,0, 0.1);
+  TStatToolkit::Constrain1D(fstringFast,"TwistY",paramC,covarC,0, 0.1);
+  TString strFitConst=TStatToolkit::MakeFitString(fstringFast, paramC,covar);
+  fAlignTree->SetAlias("fitYFastC",strFitConst.Data());
+  CreateAlignTime(fstringFast,paramC);
+
+
+}
+
+
+void AliTPCPreprocessorOffline::MakeChainTime(){
+  //
+  TFile f("CalibObjects.root");
+  //  const char *cdtype[7]={"ITS","TRD","Vertex","TOF","TPC","TPC0","TPC1"};
+  //const char *cptype[5]={"dy","dz","dsnp","dtheta","d1pt"}; 
+  const char * hname[5]={"dy","dz","dsnp","dtheta","d1pt"};
+  Int_t run=0;
+  AliTPCcalibTime  *calibTime= (AliTPCcalibTime*) f.Get("calibTime");
+  if (!calibTime) return;
+  TTreeSRedirector *pcstream = new TTreeSRedirector("meanITSVertex.root");
+  //
+  Int_t ihis=0;
+  THnSparse *his = calibTime->GetResHistoTPCITS(ihis);
+  if (his){
+    his->GetAxis(1)->SetRangeUser(-1.1,1.1);
+    his->GetAxis(2)->SetRange(0,1000000);
+    his->GetAxis(3)->SetRangeUser(-0.35,0.35);
+    AliTPCCorrection::MakeDistortionMap(his,pcstream, Form("ITS%s",hname[ihis]),run,85.,ihis,5);
+  }
+  ihis=1;
+  his = calibTime->GetResHistoTPCITS(ihis);
+  if (his){
+    his->GetAxis(1)->SetRangeUser(-1.1,1.1);
+    his->GetAxis(2)->SetRange(0,1000000);
+    his->GetAxis(3)->SetRangeUser(-0.35,0.35);
+    AliTPCCorrection::MakeDistortionMap(his,pcstream, Form("ITS%s",hname[ihis]),run,85.,ihis,5);
+  }
+  ihis=2;
+  his = calibTime->GetResHistoTPCITS(ihis);
+  if (his){
+    his->GetAxis(1)->SetRangeUser(-1.1,1.1);
+    his->GetAxis(2)->SetRange(0,1000000);
+    his->GetAxis(3)->SetRangeUser(-0.35,0.35);
+    AliTPCCorrection::MakeDistortionMap(his,pcstream, Form("ITS%s",hname[ihis]),run,85.,ihis,5);
+  }
+  ihis=0;
+  his = calibTime->GetResHistoTPCvertex(ihis);
+  if (his){
+    his->GetAxis(1)->SetRangeUser(-1.1,1.1);
+    his->GetAxis(2)->SetRange(0,1000000);
+    his->GetAxis(3)->SetRangeUser(-0.35,0.35);
+    AliTPCCorrection::MakeDistortionMap(his,pcstream, Form("Vertex%s",hname[ihis]),run,0.,ihis,5);
+  }
+  ihis=2;
+  his = calibTime->GetResHistoTPCvertex(ihis);
+  if (his){
+    his->GetAxis(1)->SetRangeUser(-1.1,1.1);
+    his->GetAxis(2)->SetRange(0,1000000);
+    his->GetAxis(3)->SetRangeUser(-0.35,0.35);
+    AliTPCCorrection::MakeDistortionMap(his,pcstream, Form("Vertex%s",hname[ihis]),run,0.,ihis,5);
+
+  }
+  ihis=1;
+  his = calibTime->GetResHistoTPCvertex(ihis);
+  if (his){
+    his->GetAxis(1)->SetRangeUser(-1.1,1.1);
+    his->GetAxis(2)->SetRange(0,1000000);
+    his->GetAxis(3)->SetRangeUser(-0.35,0.35);
+    AliTPCCorrection::MakeDistortionMap(his,pcstream, Form("Vertex%s",hname[ihis]),run,0.,ihis,5);
+
+  }
+  ihis=0;
+  his = calibTime->GetResHistoTPCTOF(ihis);
+  if (his){
+    his->GetAxis(1)->SetRangeUser(-1.1,1.1);
+    his->GetAxis(2)->SetRange(0,1000000);
+    his->GetAxis(3)->SetRangeUser(-0.35,0.35);
+    AliTPCCorrection::MakeDistortionMap(his,pcstream, Form("TOF%s",hname[ihis]),run,0.,ihis,10);
+
+  }
+  ihis=0;
+  his = calibTime->GetResHistoTPCTRD(ihis);
+  if (his){
+    his->GetAxis(1)->SetRangeUser(-1.1,1.1);
+    his->GetAxis(2)->SetRange(0,1000000);
+    his->GetAxis(3)->SetRangeUser(-0.35,0.35);
+    AliTPCCorrection::MakeDistortionMap(his,pcstream, Form("TRD%s",hname[ihis]),run,0.,ihis,10);
+
+  }
+  delete pcstream;
+}
+
+
+Double_t AliTPCPreprocessorOffline::EvalAt(Double_t phi, Double_t refX, Double_t theta, Int_t corr, Int_t ptype){
+  //
+  //
+  //
+  Double_t sector = 9*phi/TMath::Pi();
+  if (sector<0) sector+=18;
+  Double_t y85=AliTPCCorrection::GetCorrSector(sector,85,theta,1,corr);
+  Double_t y245=AliTPCCorrection::GetCorrSector(sector,245,theta,1,corr);
+  if (ptype==0) return y85+(y245-y85)*(refX-85.)/(245.-85.);
+  if (ptype==2) return (y245-y85)/(245.-85.);
+  return 0;
+}
+
+
+
+void AliTPCPreprocessorOffline::MakePrimitivesTime(){
+  //
+  // Create primitive transformation to fit
+  //
+  fAlignTree=new TChain("fit","fit");
+  fAlignTree->AddFile("meanITSVertex.root",10000000,"ITSdy");
+  fAlignTree->AddFile("meanITSVertex.root",10000000,"ITSdsnp");
+  fAlignTree->AddFile("meanITSVertex.root",10000000,"Vertexdy");
+  fAlignTree->AddFile("meanITSVertex.root",10000000,"Vertexdsnp");
+  // 
+  AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters();
+  Double_t bzField=AliTrackerBase::GetBz(); 
+  Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us]   // From dataBase: to be updated: per second (ideally)
+  Double_t ezField = 400; // [V/cm]   // to be updated: never (hopefully)
+  Double_t wtP = -10.0 * (bzField) * vdrift /  ezField ; 
+  AliTPCExBTwist *fitExBTwistX= new  AliTPCExBTwist;
+  AliTPCExBTwist *fitExBTwistY= new  AliTPCExBTwist;
+  AliTPCCalibGlobalMisalignment *trans0   =new  AliTPCCalibGlobalMisalignment;
+  AliTPCCalibGlobalMisalignment *trans1   =new  AliTPCCalibGlobalMisalignment;
+  AliTPCCalibGlobalMisalignment *trans0D  =new  AliTPCCalibGlobalMisalignment;
+  AliTPCCalibGlobalMisalignment *trans1D  =new  AliTPCCalibGlobalMisalignment;
+  AliTPCCalibGlobalMisalignment *rot0     =new  AliTPCCalibGlobalMisalignment;
+  AliTPCCalibGlobalMisalignment *rot1     =new  AliTPCCalibGlobalMisalignment;
+  AliTPCCalibGlobalMisalignment *rot2     =new  AliTPCCalibGlobalMisalignment;
+  AliTPCCalibGlobalMisalignment *rot3     =new  AliTPCCalibGlobalMisalignment;
+  //
+  //
+  fitExBTwistX->SetXTwist(0.001);
+  fitExBTwistX->SetOmegaTauT1T2(wtP,1,1);  
+  //
+  fitExBTwistY->SetYTwist(0.001);
+  fitExBTwistY->SetOmegaTauT1T2(wtP,1,1);  
+  //
+  TGeoHMatrix *matrixRot = new TGeoHMatrix; 
+  TGeoHMatrix *matrixX = new TGeoHMatrix; 
+  TGeoHMatrix *matrixY = new TGeoHMatrix; 
+  matrixX->SetDx(0.1);
+  matrixY->SetDy(0.1);
+  Double_t rotAngles0[9]={0};
+  Double_t rotAngles1[9]={0};
+  Double_t rotAngles2[9]={0};
+  //
+  Double_t rotAngles3[9]={0};
+
+  rotAngles0[0]=1; rotAngles0[4]=1; rotAngles0[8]=1;
+  rotAngles1[0]=1; rotAngles1[4]=1; rotAngles1[8]=1;
+  rotAngles2[0]=1; rotAngles2[4]=1; rotAngles2[8]=1;
+  rotAngles3[0]=1; rotAngles3[4]=1; rotAngles3[8]=1;
+
+  rotAngles0[1]=-0.001;rotAngles0[3]=0.001;
+  rotAngles1[5]=-0.001;rotAngles1[7]=0.001;
+  rotAngles2[2]=0.001;rotAngles2[6]=-0.001;
+  rotAngles3[1]=0.001;rotAngles3[3]=-0.001;
+  matrixRot->SetRotation(rotAngles0);
+  rot0->SetAlignGlobal(matrixRot);
+  matrixRot->SetRotation(rotAngles1);
+  rot1->SetAlignGlobal(matrixRot);
+  matrixRot->SetRotation(rotAngles2);
+  rot2->SetAlignGlobal(matrixRot); 
+  matrixRot->SetRotation(rotAngles3);
+  rot3->SetAlignGlobalDelta(matrixRot); 
+  //
+  trans0->SetAlignGlobal(matrixX);
+  trans1->SetAlignGlobal(matrixY);
+  trans0D->SetAlignGlobalDelta(matrixX);
+  trans1D->SetAlignGlobalDelta(matrixY);
+  fitExBTwistX->Init();
+  fitExBTwistY->Init();
+  //
+  fitExBTwistX->AddVisualCorrection((AliTPCExBTwist*)(fitExBTwistX->Clone()),100);
+  fitExBTwistY->AddVisualCorrection((AliTPCExBTwist*)(fitExBTwistY->Clone()),101);
+  //
+  fitExBTwistX->AddVisualCorrection((AliTPCExBTwist*)(rot0->Clone()),102);
+  fitExBTwistY->AddVisualCorrection((AliTPCExBTwist*)(rot1->Clone()),103);
+  fitExBTwistX->AddVisualCorrection((AliTPCExBTwist*)(rot2->Clone()),104);
+  fitExBTwistX->AddVisualCorrection((AliTPCExBTwist*)(rot3->Clone()),105);
+
+  fitExBTwistX->AddVisualCorrection((AliTPCExBTwist*)(trans0->Clone()),106);
+  fitExBTwistX->AddVisualCorrection((AliTPCExBTwist*)(trans1->Clone()),107);
+  fitExBTwistX->AddVisualCorrection((AliTPCExBTwist*)(trans0D->Clone()),108);
+  fitExBTwistX->AddVisualCorrection((AliTPCExBTwist*)(trans1D->Clone()),109);
+  //
+  fAlignTree->SetAlias("FExBTwistX", "AliTPCPreprocessorOffline::EvalAt(phi,refX,theta,100,ptype)+0");
+  fAlignTree->SetAlias("FExBTwistY","AliTPCPreprocessorOffline::EvalAt(phi,refX,theta,101,ptype)+0");
+  fAlignTree->SetAlias("FAlignRot0","AliTPCPreprocessorOffline::EvalAt(phi,refX,theta,102,ptype)+0");
+  fAlignTree->SetAlias("FAlignRot0D","AliTPCPreprocessorOffline::EvalAt(phi,refX,theta,105,ptype)+0");
+  fAlignTree->SetAlias("FAlignRot1","AliTPCPreprocessorOffline::EvalAt(phi,refX,theta,103,ptype)+0");
+  fAlignTree->SetAlias("FAlignRot2","AliTPCPreprocessorOffline::EvalAt(phi,refX,theta,104,ptype)+0");
+  fAlignTree->SetAlias("FAlignTrans0","AliTPCPreprocessorOffline::EvalAt(phi,refX,theta,106,ptype)+0");
+  fAlignTree->SetAlias("FAlignTrans1","AliTPCPreprocessorOffline::EvalAt(phi,refX,theta,107,ptype)+0");
+  fAlignTree->SetAlias("FAlignTrans0D","AliTPCPreprocessorOffline::EvalAt(phi,refX,theta,108,ptype)+0");
+  fAlignTree->SetAlias("FAlignTrans1D","AliTPCPreprocessorOffline::EvalAt(phi,refX,theta,109,ptype)+0");
+  //
+  // test fast function
+  //
+//   fAlignTree->Draw("FExBTwistX:ExBTwistX","isITS&&ptype==0&&abs(snp)<0.05","");
+//   fAlignTree->Draw("FExBTwistY:ExBTwistY","isITS&&ptype==0&&abs(snp)<0.05","");
+//   fAlignTree->Draw("FAlignRot0:alignRot0","isITS&&ptype==0&&abs(snp)<0.05","");
+//   fAlignTree->Draw("FAlignRot1:alignRot1","isITS&&ptype==0&&abs(snp)<0.05","");
+//   fAlignTree->Draw("FAlignRot2:alignRot2","isITS&&ptype==0&&abs(snp)<0.05","");
+//   //
+//   fAlignTree->Draw("FAlignTrans0:alignTrans0","isITS&&ptype==0&&abs(snp)<0.05","");
+//   fAlignTree->Draw("FAlignTrans1:alignTrans1","isITS&&ptype==0&&abs(snp)<0.05","");
+
+} 
+
+
+void AliTPCPreprocessorOffline::CreateAlignTime(TString fstring, TVectorD paramC){
+  //
+  //
+  //
+  //
+  TGeoHMatrix *matrixDelta     = new TGeoHMatrix; 
+  TGeoHMatrix *matrixGlobal    = new TGeoHMatrix; 
+  Double_t rAngles[9];
+  Int_t index=0;
+  //
+  index=TStatToolkit::GetFitIndex(fstring,"FAlignTrans0D");
+  if (index>=0) matrixDelta->SetDx(paramC[index+1]*0.1);
+  index=TStatToolkit::GetFitIndex(fstring,"FAlignTrans1D");
+  if (index>=0) matrixDelta->SetDy(paramC[index+1]*0.1);
+  rAngles[0]=1; rAngles[4]=1; rAngles[8]=1;
+  index=TStatToolkit::GetFitIndex(fstring,"FAlignRot0D");
+  rAngles[1]=-paramC[index+1]*0.001; rAngles[3]=paramC[index+1]*0.001;
+  rAngles[5]=0; rAngles[7] =0;
+  rAngles[2]=0; rAngles[6] =0;
+  matrixDelta->SetRotation(rAngles);
+  //
+  //
+  //
+  index=TStatToolkit::GetFitIndex(fstring,"hasITS*FAlignTrans0");
+  if (index>=0) matrixGlobal->SetDx(paramC[index+1]*0.1);
+  index=TStatToolkit::GetFitIndex(fstring,"hasITS*FAlignTrans1");
+  if (index>=0) matrixGlobal->SetDy(paramC[index+1]*0.1);
+  rAngles[0]=1; rAngles[4]=1; rAngles[8]=1;
+  index=TStatToolkit::GetFitIndex(fstring,"hasITS*FAlignRot0");
+  rAngles[1]=-paramC[index+1]*0.001; rAngles[3]=paramC[index+1]*0.001;
+  index=TStatToolkit::GetFitIndex(fstring,"hasITS*FAlignRot1");  
+  rAngles[5]=-paramC[index+1]*0.001; rAngles[7]=paramC[index+1]*0.001;
+  index=TStatToolkit::GetFitIndex(fstring,"hasITS*FAlignRot2");  
+  rAngles[2]=paramC[index+1]*0.001; rAngles[6] =-paramC[index+1]*0.001;
+  matrixGlobal->SetRotation(rAngles);
+  //
+  AliTPCCalibGlobalMisalignment *fitAlignTime  =0;
+  fitAlignTime  =new  AliTPCCalibGlobalMisalignment;
+  fitAlignTime->SetName("FitAlignTime");
+  fitAlignTime->SetTitle("FitAlignTime");
+  fitAlignTime->SetAlignGlobalDelta(matrixDelta);
+  fitAlignTime->SetAlignGlobal(matrixGlobal);
+  //
+  AliTPCExBTwist * fitExBTwist= new  AliTPCExBTwist;
+  Int_t indexX=TStatToolkit::GetFitIndex(fstring,"ExBTwistX");
+  Int_t indexY=TStatToolkit::GetFitIndex(fstring,"ExBTwistY");  
+  fitExBTwist->SetXTwist(0.001*paramC[indexX+1]);  // 1 mrad twist in x
+  fitExBTwist->SetYTwist(0.001*paramC[indexY+1]);  // 1 mrad twist in x
+  fitExBTwist->SetName("FitExBTwistTime");
+  fitExBTwist->SetTitle("FitExBTwistTime"); 
+  AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters();
+  Double_t bzField=AliTrackerBase::GetBz();
+  Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us]   // From dataBase: to be updated: per second (ideally)
+
+  Double_t ezField = 400; // [V/cm]   // to be updated: never (hopefully)
+  Double_t wt = -10.0 * (bzField) * vdrift /  ezField ; 
+  //
+  fitExBTwist->SetOmegaTauT1T2(wt,1,1);  
+  fitExBTwist->Init();  
+
+  AliTPCComposedCorrection *corrTime =  new AliTPCComposedCorrection;
+  TObjArray *arr = new TObjArray;
+  corrTime->SetCorrections(arr);
+  
+  corrTime->GetCorrections()->Add(fitExBTwist);
+  corrTime->GetCorrections()->Add(fitAlignTime);
+  corrTime->SetName("FitCorrectionTime");
+  corrTime->SetTitle("FitCorrectionTime");
+
+  fitExBTwist->AddVisualCorrection((AliTPCExBTwist*)(fitExBTwist->Clone()),1001);
+  fitAlignTime->AddVisualCorrection((AliTPCExBTwist*)(fitAlignTime->Clone()),1002);
+  fitAlignTime->AddVisualCorrection((AliTPCExBTwist*)(corrTime->Clone()),1003);
+  
+  
+  fAlignTree->SetAlias("ExBTwistTime","AliTPCPreprocessorOffline::EvalAt(phi,refX,theta,1001,ptype)+0");
+  fAlignTree->SetAlias("AlignTime","AliTPCPreprocessorOffline::EvalAt(phi,refX,theta,1002,ptype)+0");
+  fAlignTree->SetAlias("FitCorrectionTime","AliTPCPreprocessorOffline::EvalAt(phi,refX,theta,1003,ptype)+0");
+
+
+  TFile *f = new TFile("fitITSVertex.root","update");
+  corrTime->Write("FitCorrectionTime");
+  f->Close();
+}
+
+