]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TPC/AliTPCPreprocessorOffline.cxx
Compatibility with ROOT trunk
[u/mrichter/AliRoot.git] / TPC / AliTPCPreprocessorOffline.cxx
index 428fcd3f9eaf20eb428a2e9d4cbb57d4f58c2b40..6af31a50e1c5df08e3c219fa03fdbe89fbfdcbb2 100644 (file)
@@ -75,6 +75,7 @@
 #include "AliTPCParamSR.h"
 #include "AliTPCcalibTimeGain.h"
 #include "AliTPCcalibGainMult.h"
+#include "AliTPCcalibAlign.h"
 #include "AliSplineFit.h"
 #include "AliTPCComposedCorrection.h"
 #include "AliTPCExBTwist.h"
 #include "TChain.h"
 #include "TCut.h"
 #include "AliTrackerBase.h"
+#include "AliTracker.h"
 #include "AliTPCPreprocessorOffline.h"
+#include "AliTPCCorrectionFit.h"
 
-
+using std::endl;
+using std::cout;
 ClassImp(AliTPCPreprocessorOffline)
 
 AliTPCPreprocessorOffline::AliTPCPreprocessorOffline():
@@ -108,7 +112,10 @@ AliTPCPreprocessorOffline::AliTPCPreprocessorOffline():
   fGainCosmic(0),       // calibration component for cosmic
   fGainMult(0),
   fAlignTree(0),        // alignment tree
-  fSwitchOnValidation(kFALSE) // flag to switch on validation of OCDB parameters
+  fSwitchOnValidation(kFALSE), // flag to switch on validation of OCDB parameters
+  fMinGain(2.0),
+  fMaxGain(3.0),
+  fMaxVdriftCorr(0.03)
 {
   //
   // default constructor
@@ -131,8 +138,8 @@ void AliTPCPreprocessorOffline::GetRunRange(AliTPCcalibTime * const  timeDrift){
   TObjArray *hisArray =timeDrift->GetHistoDrift();
   {for (Int_t i=0; i<hisArray->GetEntriesFast(); i++){
     THnSparse* addHist=(THnSparse*)hisArray->UncheckedAt(i);
-    if (addHist->GetEntries()<fMinEntries) continue;
     if (!addHist) continue;
+    if (addHist->GetEntries()<fMinEntries) continue;
     TH1D* histo    =addHist->Projection(3);
     TH1D* histoTime=addHist->Projection(0);
     printf("%s\t%f\t%d\t%d\n",histo->GetName(), histo->GetEntries(),histo->FindFirstBinAbove(0),histo->FindLastBinAbove(0));
@@ -248,15 +255,21 @@ void AliTPCPreprocessorOffline::UpdateOCDBDrift( Int_t ustartRun, Int_t uendRun,
   gStorage->Put(fVdriftArray, (*id1), metaData);
 }
 
-Bool_t AliTPCPreprocessorOffline::ValidateTimeGain(Double_t minGain, Double_t maxGain)
+Bool_t AliTPCPreprocessorOffline::ValidateTimeGain()
 {
   //
   // Validate time gain corrections 
   //
   Printf("ValidateTimeGain..." );
+  Float_t minGain = fMinGain;
+  Float_t maxGain = fMaxGain;
 
   TGraphErrors *gr = (TGraphErrors*)fGainArray->FindObject("TGRAPHERRORS_MEAN_GAIN_BEAM_ALL");
-  if(!gr) return kFALSE;
+  if (!gr) {
+    gr = (TGraphErrors*)fGainArray->FindObject("TGRAPHERRORS_MEAN_GAIN_COSMIC_ALL");
+    if (!gr) return kFALSE;
+    Printf("Assuming given run is a cosmic run. Using gain calibration from Fermi-plateau muons.");
+  }
   if(gr->GetN()<1) return kFALSE;
 
   // check whether gain in the range
@@ -270,20 +283,28 @@ return kTRUE;
 }
 
 
-Bool_t AliTPCPreprocessorOffline::ValidateTimeDrift(Double_t maxVDriftCorr)
+Bool_t AliTPCPreprocessorOffline::ValidateTimeDrift()
 {
   //
   // Validate time drift velocity corrections 
   //
   Printf("ValidateTimeDrift..." );
 
+  Float_t maxVDriftCorr = fMaxVdriftCorr;
+
   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;
   }
@@ -379,7 +400,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++){
@@ -426,6 +453,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);
@@ -434,10 +465,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);
@@ -471,9 +498,9 @@ void AliTPCPreprocessorOffline::AddAlignmentGraphs(  TObjArray * vdriftArray, Al
   arrayTRD=timeDrift->GetAlignTRDTPC();
   arrayTOF=timeDrift->GetAlignTOFTPC();
 
-  if (arrayITS->GetEntries()>0) mstatITS= AliTPCcalibDButil::MakeStatRelKalman(arrayITS,0.9,50,0.025);
-  if (arrayTOF->GetEntries()>0) mstatTOF= AliTPCcalibDButil::MakeStatRelKalman(arrayTOF,0.9,1000,0.025);
-  if (arrayTRD->GetEntries()>0) mstatTRD= AliTPCcalibDButil::MakeStatRelKalman(arrayTRD,0.9,50,0.025);
+  if (arrayITS->GetEntries()>0) mstatITS= AliTPCcalibDButil::MakeStatRelKalman(arrayITS,0.7,50,fMaxVdriftCorr);
+  if (arrayTOF->GetEntries()>0) mstatTOF= AliTPCcalibDButil::MakeStatRelKalman(arrayTOF,0.7,1000,fMaxVdriftCorr);
+  if (arrayTRD->GetEntries()>0) mstatTRD= AliTPCcalibDButil::MakeStatRelKalman(arrayTRD,0.7,50,fMaxVdriftCorr);
   //
   TObjArray * arrayITSP= AliTPCcalibDButil::SmoothRelKalman(arrayITS,*mstatITS, 0, 5.);
   TObjArray * arrayITSM= AliTPCcalibDButil::SmoothRelKalman(arrayITS,*mstatITS, 1, 5.);
@@ -544,9 +571,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");
@@ -815,6 +842,7 @@ void AliTPCPreprocessorOffline::CalibTimeGain(const Char_t* fileName, Int_t star
   AnalyzeAttachment(startRunNumber,endRunNumber);
   AnalyzePadRegionGain();
   AnalyzeGainMultiplicity();
+  AnalyzeGainChamberByChamber();
   //
   // 3. Make control plots
   //
@@ -1084,7 +1112,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;
@@ -1099,7 +1128,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;
@@ -1118,6 +1148,26 @@ Bool_t AliTPCPreprocessorOffline::AnalyzeGainMultiplicity() {
 
 }
 
+Bool_t AliTPCPreprocessorOffline::AnalyzeGainChamberByChamber(){
+  //
+  // get chamber by chamber gain
+  //
+  TGraphErrors *grShort  = fGainMult->GetGainPerChamber(0);
+  TGraphErrors *grMedium = fGainMult->GetGainPerChamber(1);
+  TGraphErrors *grLong   = fGainMult->GetGainPerChamber(2);
+  if (grShort==0x0 || grMedium==0x0 || grLong==0x0) {
+    delete grShort;
+    delete grMedium;
+    delete grLong;
+    return kFALSE;
+  }
+
+  fGainArray->AddLast(grShort);
+  fGainArray->AddLast(grMedium);
+  fGainArray->AddLast(grLong);
+
+  return kTRUE;
+}
 
 void AliTPCPreprocessorOffline::UpdateOCDBGain(Int_t startRunNumber, Int_t endRunNumber, const Char_t *storagePath){
   //
@@ -1161,9 +1211,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;    
@@ -1190,7 +1242,7 @@ void AliTPCPreprocessorOffline::MakeQAPlot(Float_t  FPtoMIPratio) {
 
 void AliTPCPreprocessorOffline::MakeFitTime(){
   //
-  // mak aligment fit - store results in the file
+  // make aligment fit - store results in the file
   //
   const Int_t kMinEntries=1000;
   MakeChainTime();
@@ -1229,7 +1281,7 @@ void AliTPCPreprocessorOffline::MakeFitTime(){
   fstringFast+="dITS*FAlignRot1++";
   fstringFast+="dITS*FAlignRot2++";
   
-  TCut cutFit="entries>10&&abs(mean)>0.00001";
+  TCut cutFit="entries>10&&abs(mean)>0.00001&&rms>0";
   fAlignTree->SetAlias("err","rms");
 
   TString *strDeltaITS = TStatToolkit::FitPlaneConstrain(fAlignTree,"mean:err", fstringFast.Data(),cutFit, chi2,npoints,param,covar,-1,0, npointsMax, 1);
@@ -1252,21 +1304,40 @@ void AliTPCPreprocessorOffline::MakeFitTime(){
 
 void AliTPCPreprocessorOffline::MakeChainTime(){
   //
-  TFile f("TPCTimeObjects.root");
+  //
+  //
+  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");
+  AliTPCcalibTime  *calibTime = 0;
+  TObjArray * array = (TObjArray*)f.Get("TPCCalib");
+  if (array){
+    calibTime = (AliTPCcalibTime *)array->FindObject("calibTime");
+  } else {
+    calibTime = (AliTPCcalibTime*)f.Get("calibTime");
+  }
   if (!calibTime) return;
+  AliTPCCorrectionFit::CreateAlignMaps(AliTracker::GetBz(), run);
   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);
+    AliTPCCorrection::MakeDistortionMap(his,pcstream, Form("ITS%s",hname[ihis]),run,85.,ihis,3);
+  }
+  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,3);
   }
   ihis=2;
   his = calibTime->GetResHistoTPCITS(ihis);
@@ -1274,7 +1345,7 @@ void AliTPCPreprocessorOffline::MakeChainTime(){
     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);
+    AliTPCCorrection::MakeDistortionMap(his,pcstream, Form("ITS%s",hname[ihis]),run,85.,ihis,3);
   }
   ihis=0;
   his = calibTime->GetResHistoTPCvertex(ihis);
@@ -1282,7 +1353,7 @@ void AliTPCPreprocessorOffline::MakeChainTime(){
     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);
+    AliTPCCorrection::MakeDistortionMap(his,pcstream, Form("Vertex%s",hname[ihis]),run,0.,ihis,3);
   }
   ihis=2;
   his = calibTime->GetResHistoTPCvertex(ihis);
@@ -1290,7 +1361,34 @@ void AliTPCPreprocessorOffline::MakeChainTime(){
     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);
+    AliTPCCorrection::MakeDistortionMap(his,pcstream, Form("Vertex%s",hname[ihis]),run,0.,ihis,3);
+
+  }
+  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,3);
+
+  }
+  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,3);
+
+  }
+  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,3);
 
   }
   delete pcstream;
@@ -1311,6 +1409,46 @@ Double_t AliTPCPreprocessorOffline::EvalAt(Double_t phi, Double_t refX, Double_t
 }
 
 
+Double_t AliTPCPreprocessorOffline::EvalAtPar(Double_t phi0, Double_t snp, Double_t refX, Double_t theta, Int_t corr, Int_t ptype, Int_t nsteps){
+  //
+  // Fit the distortion along the line with the parabolic model
+  // Parameters:
+  //  phi0 - phi at the entrance of the TPC
+  //  snp  - local inclination angle at the entrance of the TPC
+  //  refX - ref X where the distortion is evanluated
+  //  theta
+  //  
+  static TLinearFitter fitter(3,"pol2"); 
+  fitter.ClearPoints();
+  if (nsteps<3) nsteps=3;
+  Double_t deltaX=(245-85)/(nsteps);
+  for (Int_t istep=0; istep<(nsteps+1); istep++){
+    //
+    Double_t localX =85.+deltaX*istep;
+    Double_t localPhi=phi0+deltaX*snp*istep;
+    Double_t sector = 9*localPhi/TMath::Pi();
+    if (sector<0) sector+=18;
+    Double_t y=AliTPCCorrection::GetCorrSector(sector,localX,theta,1,corr);
+    Double_t dlocalX=AliTPCCorrection::GetCorrSector(sector,localX,theta,0,corr);
+    Double_t x[1]={localX-dlocalX};
+    fitter.AddPoint(x,y);
+  }
+  fitter.Eval();
+  Double_t par[3];
+  par[0]=fitter.GetParameter(0);
+  par[1]=fitter.GetParameter(1);
+  par[2]=fitter.GetParameter(2);
+
+  if (ptype==0) return par[0]+par[1]*refX+par[2]*refX*refX;
+  if (ptype==2) return par[1]+2*par[2]*refX;
+  if (ptype==4) return par[2];
+  return 0;
+}
+
+
+
+
+
 
 void AliTPCPreprocessorOffline::MakePrimitivesTime(){
   //
@@ -1455,7 +1593,7 @@ void AliTPCPreprocessorOffline::CreateAlignTime(TString fstring, TVectorD paramC
   rAngles[2]=paramC[index+1]*0.001; rAngles[6] =-paramC[index+1]*0.001;
   matrixGlobal->SetRotation(rAngles);
   //
-  AliTPCCalibGlobalMisalignment *fitAlignTime  =new  AliTPCCalibGlobalMisalignment;
+  AliTPCCalibGlobalMisalignment *fitAlignTime  =0;
   fitAlignTime  =new  AliTPCCalibGlobalMisalignment;
   fitAlignTime->SetName("FitAlignTime");
   fitAlignTime->SetTitle("FitAlignTime");
@@ -1503,6 +1641,3 @@ void AliTPCPreprocessorOffline::CreateAlignTime(TString fstring, TVectorD paramC
   f->Close();
 }
 
-
-
-