]> git.uio.no Git - u/mrichter/AliRoot.git/commitdiff
new performance classes (Simone Schuchmann)
authorjotwinow <jotwinow@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sat, 7 Nov 2009 04:05:14 +0000 (04:05 +0000)
committerjotwinow <jotwinow@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sat, 7 Nov 2009 04:05:14 +0000 (04:05 +0000)
PWG1/PWG1LinkDef.h
PWG1/TPC/AliPerfAnalyzeInvPt.cxx [new file with mode: 0755]
PWG1/TPC/AliPerfAnalyzeInvPt.h [new file with mode: 0755]
PWG1/TPC/AliPerformancePtCalib.cxx [new file with mode: 0755]
PWG1/TPC/AliPerformancePtCalib.h [new file with mode: 0755]
PWG1/TPC/AliPerformancePtCalibMC.cxx [new file with mode: 0755]
PWG1/TPC/AliPerformancePtCalibMC.h [new file with mode: 0755]

index 0ac305fb63ac284b3ac120ff05a13726b6966779..94649acab6c9d2538dfcf9ede28ca764a8809a00 100644 (file)
@@ -48,6 +48,9 @@
 #pragma link C++ class AliPerformanceTPC+;
 #pragma link C++ class AliPerformanceMC+;
 #pragma link C++ class AliPerformanceMatch+;
+#pragma link C++ class AliPerformancePtCalib+;
+#pragma link C++ class AliPerformancePtCalibMC+;
+#pragma link C++ class AliPerfAnalyzeInvPt+;
 
 #pragma link C++ class AliAlignmentDataFilterITS+;
 #pragma link C++ class AliTrackMatchingTPCITSCosmics+;
diff --git a/PWG1/TPC/AliPerfAnalyzeInvPt.cxx b/PWG1/TPC/AliPerfAnalyzeInvPt.cxx
new file mode 100755 (executable)
index 0000000..d9559ca
--- /dev/null
@@ -0,0 +1,514 @@
+#include "TNamed.h"
+#include "TList.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TH3F.h"
+#include "TCanvas.h"
+#include "TMath.h"
+#include "TGraphErrors.h"
+#include "TF1.h"
+
+//#include "AliPerformancePtCalibMC.h"
+//#include "AliPerformancePtCalib.h"
+#include "AliPerfAnalyzeInvPt.h"
+
+ClassImp(AliPerfAnalyzeInvPt)
+
+
+// fit functions
+//____________________________________________________________________________________________________________________________________________
+   Double_t AliPerfAnalyzeInvPt::Polynomial(Double_t *x, Double_t *par)
+{
+
+
+   if (x[0] > -par[4] && x[0] < par[4]) {
+      TF1::RejectPoint();
+      return 0;
+   }
+   return  par[2]+par[0]*pow((x[0]-par[1]),2)+par[3]*pow((x[0]-par[1]),4);
+  
+}
+//____________________________________________________________________________________________________________________________________________
+Double_t AliPerfAnalyzeInvPt::PolynomialRejP(Double_t *x, Double_t *par)
+{
+   Double_t pos  = par[5];
+   Double_t neg = - par[5];
+   pos += par[4];
+   neg += par[4];
+  
+   if (x[0] > neg && x[0] < pos) {
+      TF1::RejectPoint();
+      return 0;
+   }
+   return   par[2]+par[0]*pow((x[0]-par[1]),2)+par[3]*pow((x[0]-par[1]),4);
+}
+//____________________________________________________________________________________________________________________________________________
+Double_t AliPerfAnalyzeInvPt::InvGauss(Double_t *x, Double_t *par)
+{
+   if (x[0] > -par[6] && x[0] < par[6]) {
+      TF1::RejectPoint();
+      return 0;
+   }
+
+   return par[3]+par[0]*TMath::Exp(-0.5*(TMath::Power((x[0]-par[1])/par[2], 2.0)))+par[4]*pow((x[0]-par[1]),2)+par[5]*pow((x[0]-par[1]),4) ;
+}
+//____________________________________________________________________________________________________________________________________________
+Double_t AliPerfAnalyzeInvPt::InvGaussRejP(Double_t *x, Double_t *par)
+{
+   Double_t pos  = par[7];//0.12;
+   Double_t neg = - par[7];//0.12;
+   pos += par[6];
+   neg += par[6];
+  
+   if (x[0] > neg && x[0] < pos) {
+      TF1::RejectPoint();
+      return 0;
+   }
+
+
+   return par[3]+par[0]*TMath::Exp(-0.5*(TMath::Power((x[0]-par[1])/par[2], 2.0)))+par[4]*pow((x[0]-par[1]),2)+par[5]*pow((x[0]-par[1]),4) ;
+}
+
+
+//_____________________________________________________________________________________________________________________________________________
+AliPerfAnalyzeInvPt::AliPerfAnalyzeInvPt():
+   TNamed("AliPerfAnalyzeInvPt","AliPerfAnalyzeInvPt"),
+   fNThetaBins(0), 
+   fNPhiBins(0),
+   fRange(0),
+   fExclRange(0),
+   fFitGaus(0) ,
+   fHistH2InvPtTheta(0),
+   fHistH2InvPtPhi(0), 
+   fGrMinPosTheta(0),
+   fGrMinPosPhi(0),
+   fFitMinPos(0),
+   fFitMinPosRejP(0),
+   fFitInvGauss(0),
+   fFitInvGaussRejP(0)
+{
+   // Default constructor
+  
+   fFitGaus = kFALSE;
+   fNThetaBins = 0;
+   fNPhiBins = 0;
+   fRange = 0;
+   fExclRange = 0;
+   fFitGaus = 0;
+   
+   // projection histos
+   TH1D *fHistFitTheta[100];
+   TH1D *fHistFitPhi[100];
+
+   for(Int_t i=0;i<100;i++){
+      
+      fHistFitTheta[i] = NULL;
+      fHistFitPhi[i] = NULL;
+   }
+  
+}
+//_____________________________________________________________________________________________________________________________________________
+AliPerfAnalyzeInvPt::AliPerfAnalyzeInvPt(Char_t* name="AliAnalyzeInvPt",Char_t* title="AliAnalyzeInvPt"): TNamed(name, title),
+   fNThetaBins(0), 
+   fNPhiBins(0),
+   fRange(0),
+   fExclRange(0),
+   fFitGaus(0) ,
+   fHistH2InvPtTheta(0),
+   fHistH2InvPtPhi(0), 
+   fGrMinPosTheta(0),
+   fGrMinPosPhi(0),
+   fFitMinPos(0),
+   fFitMinPosRejP(0),
+   fFitInvGauss(0),
+   fFitInvGaussRejP(0)
+   {
+      //  Double_t fThetaBins[100] = {0};
+      //   Double_t fPhiBins[100] = {0};
+  
+      fFitGaus = kFALSE;
+      fNThetaBins = 0;
+      fNPhiBins =0;
+      fRange = 0;
+      fExclRange = 0;
+      fFitGaus = 0;
+      // projection histos
+      TH1D *fHistFitTheta[100];
+      TH1D *fHistFitPhi[100];
+
+      for(Int_t i=0;i<100;i++){
+    
+        fHistFitTheta[i] = NULL;
+        fHistFitPhi[i] = NULL;
+      }
+      fHistH2InvPtTheta = NULL;
+      fHistH2InvPtPhi = NULL; 
+      fGrMinPosTheta= NULL;
+      fGrMinPosPhi= NULL;
+   }
+
+
+   //______________________________________________________________________________________________________________________________________
+void AliPerfAnalyzeInvPt::InitHistos(Double_t *binsXTheta,Double_t *fitParamTheta,Double_t *errFitParamTheta,Double_t *binsXPhi,Double_t *fitParamPhi,Double_t *errFitParamPhi){// init Histos
+
+
+      
+  
+      fGrMinPosTheta = new TGraphErrors(fNThetaBins,binsXTheta,fitParamTheta,0, errFitParamTheta);  
+      fGrMinPosTheta->SetMarkerStyle(20);
+      fGrMinPosTheta->SetMarkerColor(2);
+      fGrMinPosTheta->SetLineColor(2);
+      fGrMinPosTheta->GetYaxis()->SetTitle("min pos (Gev/c)^{-1}");
+      fGrMinPosTheta->GetXaxis()->SetTitle("#theta bin no.");
+      fGrMinPosTheta->GetYaxis()->SetTitleOffset(1.2);   
+      fGrMinPosTheta->SetTitle("#theta bins ");
+
+      fGrMinPosPhi = new TGraphErrors(fNPhiBins,binsXPhi,fitParamPhi,0,errFitParamPhi);  
+      fGrMinPosPhi->SetMarkerStyle(20);
+      fGrMinPosPhi->SetMarkerColor(4);
+      fGrMinPosPhi->SetLineColor(4);
+      fGrMinPosPhi->GetYaxis()->SetTitle("min pos (Gev/c)^{-1}");
+      fGrMinPosPhi->GetXaxis()->SetTitle("#phi bin no.");
+      fGrMinPosPhi->GetYaxis()->SetTitleOffset(1.2);   
+      fGrMinPosPhi->SetTitle("#phi bins ");
+}
+
+//______________________________________________________________________________________________________________________________________
+
+void AliPerfAnalyzeInvPt::InitFitFcn(){
+      // fit functions
+      fFitMinPos = new TF1("fFitMinPos", Polynomial,-4.0,4.0,5);
+      fFitMinPos->SetLineColor(4);
+      fFitMinPos->SetLineWidth(1);
+      fFitMinPos->SetParameter(0,1.0);
+      fFitMinPos->SetParameter(1,0.0);
+      fFitMinPos->SetParameter(2,0.0);
+
+      fFitMinPosRejP = new TF1("fFitMinPosRejP",PolynomialRejP,-4.0,4.0,6);
+      fFitMinPosRejP->SetLineColor(2);
+      fFitMinPosRejP->SetLineWidth(1);
+      fFitMinPosRejP->SetParameter(0,1.0);
+      fFitMinPosRejP->SetParameter(1,0.0);
+      fFitMinPosRejP->SetParameter(2,0.0);
+  
+      fFitInvGauss = new TF1("fFitInvGauss", InvGauss,-4.0,4.0,6);
+      fFitInvGauss->SetLineColor(4);
+      fFitInvGauss->SetLineWidth(1);
+      fFitInvGauss->SetParameter(2,1.0);
+  
+      fFitInvGaussRejP = new TF1("fFitInvGaussRejP", InvGaussRejP,-4.0,4.0,7);
+      fFitInvGaussRejP->SetLineColor(2);
+      fFitInvGaussRejP->SetLineWidth(1);
+      fFitInvGaussRejP->SetParameter(2,1.0);
+   }
+//______________________________________________________________________________________________________________________________________
+void AliPerfAnalyzeInvPt::StartAnalysis(TH2F *histThetaInvPt, TH2F *histPhiInvPt, TObjArray* aFolderObj){//start Ana
+
+  
+   
+   if(!histThetaInvPt) {
+      Printf("warning: no 1/pt histogram to analyse in theta bins!");
+   }
+   if(!histPhiInvPt) {
+      Printf("warning: no 1/pt histogram to analyse in phit bins!");
+   }
+
+   Double_t thetaBins[9] = {0.77,0.97,1.17,1.37,1.57,1.77,1.97,2.17,2.37};                 // theta bins
+   Double_t phiBins[13] = {0.0,0.5,1.0,1.5,2.0,2.5,3.0,3.5,4.0,4.5,5.0,5.5,6.5};           // phi bins
+  
+   Int_t nThetaBins = 9;
+   Int_t nPhiBins = 13;
+   Double_t range = 1.0;        // fit range
+   Double_t  exclRange  = 0.13; // range of point rejection in fit
+   Bool_t fitGaus=kFALSE;       // fit with Gaussian or with polynomial (kFALSE)
+   
+   if(fNThetaBins == 0){
+      fNThetaBins = nThetaBins;
+      for(Int_t k = 0;k<fNThetaBins;k++){
+        fThetaBins[k] = thetaBins[k];
+      }
+      Printf("warning: for theta bins no user intput! bins are set to default values.");
+   }
+
+   if(fNPhiBins == 0){
+      fNPhiBins = nPhiBins;
+      for(Int_t k = 0;k<fNPhiBins;k++){
+        fPhiBins[k] = phiBins[k];
+      }
+    
+      Printf("warning: for phi bins no user intput! bins are set to default values.");
+    
+   }
+
+   if(fRange==0){
+      fRange = range;
+      fExclRange = exclRange;
+      fFitGaus = fitGaus;
+      Printf("warning: no fit range is set. default fitting conditions are used.");
+   }
+
+      
+   //param arrays
+   Double_t fitParamTheta[100],fitParamPhi[100];
+   Double_t errFitParamTheta[100], errFitParamPhi[100];
+   Double_t binsXTheta[100],binsXPhi[100];
+
+   fHistH2InvPtTheta  = histThetaInvPt;
+   fHistH2InvPtPhi    = histPhiInvPt;
+
+   InitFitFcn();
+        
+   TCanvas *theCan =new TCanvas("theCan","invPt theta bins",1200,900);
+   theCan->Divide((fNThetaBins+2)/2,2);
+   TCanvas *phiCan =new TCanvas("phiCan","invPt phi bins",1200,900);
+   phiCan->Divide((fNPhiBins+2)/2,2);
+   Int_t countPad = 1;
+        
+   // analyse 1/pt in bins of theta 
+   for(Int_t i=0;i<fNThetaBins;i++){
+      TString name = "fit_theta_";
+      name +=i;
+      Int_t firstBin= fHistH2InvPtTheta->GetYaxis()->FindBin(fThetaBins[i]);
+      if(i>0) firstBin +=1;
+      Int_t lastBin = fHistH2InvPtTheta->GetYaxis()->FindBin(fThetaBins[i+1]);
+      if( i == fNThetaBins-1) {
+        firstBin= fHistH2InvPtTheta->GetYaxis()->FindBin(fThetaBins[0]);
+        lastBin = fHistH2InvPtTheta->GetYaxis()->FindBin(fThetaBins[i]);
+      }
+    
+      fHistH2InvPtTheta->SetName(name.Data());
+      fHistFitTheta[i] =  (TH1F*)fHistH2InvPtTheta->ProjectionX("_px",firstBin,lastBin,"e");
+      
+      Char_t titleTheta[50];
+      if(i == fNThetaBins-1) sprintf(titleTheta,"1/pt (GeV/c) integrated over #theta");
+      else  sprintf(titleTheta,"1/pt (GeV/c) for #theta range: %1.3f - %1.3f",fThetaBins[i],fThetaBins[i+1]);
+      
+      fHistFitTheta[i]->SetTitle(titleTheta);
+   
+      Double_t invPtMinPos  = 0;
+      Double_t invPtMinPosErr = 0;
+      Double_t invPtMinPosImpr  = 0;
+      Double_t invPtMinPosErrImpr = 0;
+   
+
+        
+      if(fFitGaus==kFALSE){
+        Printf("making polynomial fit in 1/pt in theta bins");
+        theCan->cd(countPad);
+        MakeFit(fHistFitTheta[i],fFitMinPos, invPtMinPos,invPtMinPosErr, fExclRange,fRange);
+        MakeFitBetter(fHistFitTheta[i],fFitMinPosRejP, invPtMinPosImpr, invPtMinPosErrImpr, invPtMinPos,fExclRange,fRange);
+
+       
+        fHistFitTheta[i]->DrawCopy();
+        fFitMinPos->DrawCopy("L,same");
+        fFitMinPosRejP->DrawCopy("L,same");
+      }
+      else{
+        Printf("making gauss fit in 1/pt in theta bins");
+        theCan->cd(countPad);
+        MakeFitInvGauss(fHistFitTheta[i],fFitInvGauss, invPtMinPos, invPtMinPosErr,fExclRange,fRange);
+        MakeFitInvGaussBetter(fHistFitTheta[i],fFitInvGaussRejP, invPtMinPosImpr, invPtMinPosErrImpr,invPtMinPos,fExclRange,fRange);
+      
+        fHistFitTheta[i]->DrawCopy();
+        fFitInvGauss->DrawCopy("L,same");
+        fFitInvGaussRejP->DrawCopy("L,same");
+      }
+    
+      aFolderObj->Add(fHistFitTheta[i]);
+    
+      fitParamTheta[i] = invPtMinPosImpr;
+      errFitParamTheta[i] = invPtMinPosErrImpr;
+    
+      binsXTheta[i] = i+1.0;
+      countPad++;
+   }
+      
+      
+   countPad = 1;
+
+   
+   // analyse 1/pt in bins of phi 
+  
+   for(Int_t i=0;i<fNPhiBins;i++){
+      TString name = "fit_phi_";
+      name +=i;
+    
+      fHistH2InvPtPhi->SetName(name.Data());
+      Int_t  firstBin = fHistH2InvPtPhi->GetYaxis()->FindBin(fPhiBins[i]);
+      if(i>0) firstBin +=1;
+      Int_t   lastBin =  fHistH2InvPtPhi->GetYaxis()->FindBin(fPhiBins[i+1]);
+      if(i == fNPhiBins-1){
+        firstBin = fHistH2InvPtPhi->GetYaxis()->FindBin(fPhiBins[0]);
+        lastBin =  fHistH2InvPtPhi->GetYaxis()->FindBin(fPhiBins[i]);
+      }
+      fHistFitPhi[i] =  (TH1F*) fHistH2InvPtPhi->ProjectionX("_px",firstBin,lastBin,"e");
+      
+      Char_t titlePhi[50];
+      if(i == fNPhiBins-1) sprintf(titlePhi,"1/pt (GeV/c) integrated over #phi");
+         else  sprintf(titlePhi,"1/pt (GeV/c) for #phi range: %1.3f - %1.3f",fPhiBins[i],fPhiBins[i+1]);
+     
+      fHistFitPhi[i]->SetTitle(titlePhi);
+  
+      Double_t invPtMinPos  = 0;
+      Double_t invPtMinPosErr = 0;
+      Double_t invPtMinPosImpr  = 0;
+      Double_t invPtMinPosErrImpr = 0;
+    
+      if(fFitGaus==kFALSE){
+        Printf("making polynomial fit in 1/pt in phi bins");
+        phiCan->cd(countPad);
+        MakeFit(fHistFitPhi[i],fFitMinPos, invPtMinPos, invPtMinPosErr,fExclRange,fRange);
+        MakeFitBetter(fHistFitPhi[i],fFitMinPosRejP, invPtMinPosImpr, invPtMinPosErrImpr,invPtMinPos,fExclRange,fRange);
+        
+        fHistFitPhi[i]->DrawCopy();
+        fFitMinPos->DrawCopy("L,same");
+        fFitMinPosRejP->DrawCopy("L,same");
+
+      }
+      else {
+        Printf("making gauss fit in 1/pt in phi bins");
+        phiCan->cd(countPad);
+        MakeFitInvGauss(fHistFitPhi[i],fFitInvGauss, invPtMinPos, invPtMinPosErr, exclRange,fRange);
+        MakeFitInvGaussBetter(fHistFitPhi[i],fFitInvGaussRejP, invPtMinPosImpr, invPtMinPosErrImpr, invPtMinPos, fExclRange,fRange);
+        
+        fHistFitPhi[i]->DrawCopy();
+        fFitInvGauss->DrawCopy("L,same");
+        fFitInvGaussRejP->DrawCopy("L,same");
+      }
+    
+      aFolderObj->Add(fHistFitPhi[i]);
+    
+      fitParamPhi[i] = invPtMinPosImpr;
+      errFitParamPhi[i] = invPtMinPosErrImpr;
+    
+      binsXPhi[i] = i+1.0;
+      countPad++;
+   }
+      
+   InitHistos(binsXTheta,fitParamTheta,errFitParamTheta,binsXPhi,fitParamPhi,errFitParamPhi);
+
+   TCanvas *canFitVal = new TCanvas("canFitVal","min pos histos",800,400);
+   canFitVal->Divide(2,1);
+
+   canFitVal->cd(1);
+   fGrMinPosTheta->Draw("ALP");
+   canFitVal->cd(2);
+   fGrMinPosPhi->Draw("ALP");
+
+   Printf("AliPerfAnalyzeInvPt: NOTE: last bin is always fit result  of integral over all angle ranges which have been set by user!");
+   
+   aFolderObj->Add(fGrMinPosTheta);
+   aFolderObj->Add(fGrMinPosPhi);
+   aFolderObj->Add(fHistH2InvPtTheta);
+   aFolderObj->Add(fHistH2InvPtPhi);
+  
+  
+}
+
+
+//____________________________________________________________________________________________________________________________________________
+
+void AliPerfAnalyzeInvPt::MakeFit(TH1 *hproy, TF1 * fitpb, Double_t &mean, Double_t &errMean, Double_t &excl,Double_t &range)
+{
+   fitpb->SetRange(-range,range);
+   fitpb->SetParLimits(1,-0.05,0.05);
+   fitpb->SetParameter(0,1.0);
+   fitpb->FixParameter(4,excl);
+   fitpb->FixParameter(2,0.0);
+    
+   hproy->Fit(fitpb,"RQM");
+   mean = fitpb->GetParameter(1);
+   errMean = fitpb->GetParError(1);
+   if(mean == 0)  errMean = 0;
+}
+//____________________________________________________________________________________________________________________________________________
+void AliPerfAnalyzeInvPt::MakeFitBetter(TH1 *hproy, TF1 * fitpb2, Double_t &mean, Double_t &errMean, Double_t &f, Double_t &excl, Double_t &range)
+{
+   fitpb2->FixParameter(5,excl);
+   fitpb2->FixParameter(4,f);
+   // fitpb2->FixParameter(2,0.0);
+   fitpb2->SetRange(-range+f,range+f);// was 0.25 0.45
+   fitpb2->SetParLimits(1,-0.05,0.05);
+   fitpb2->SetParameter(0,1.0);
+   hproy->Fit(fitpb2,"RQM");
+   mean = fitpb2->GetParameter(1);
+   errMean = fitpb2->GetParError(1);
+   if(mean == 0)   errMean = 0;
+
+}
+//____________________________________________________________________________________________________________________________________________
+void AliPerfAnalyzeInvPt::MakeFitInvGauss(TH1 *hproy, TF1 * fitpb, Double_t &mean, Double_t &errMean, Double_t &excl,Double_t &range)
+{
+   fitpb->FixParameter(6,excl);
+   fitpb->SetRange(-range,range);
+   fitpb->SetParameter(0,-1.0);
+   fitpb->SetParameter(2,1.0);
+   fitpb->SetParameter(3,25000.0);
+   fitpb->SetParameter(4,-1.0);
+   fitpb->SetParLimits(1,-0.02,0.02);
+   hproy->Fit(fitpb,"RQM");
+   mean = fitpb->GetParameter(1);
+   errMean = fitpb->GetParError(1);
+   if(mean == 0)   errMean = 0;
+  
+}
+//____________________________________________________________________________________________________________________________________________
+void AliPerfAnalyzeInvPt::MakeFitInvGaussBetter(TH1 *hproy, TF1 * fitpb2, Double_t &mean, Double_t &errMean, Double_t &f, Double_t &excl, Double_t &range)
+{
+   fitpb2->FixParameter(7,excl);
+   fitpb2->FixParameter(6,f);
+   fitpb2->SetRange(-range+f,range+f);
+   fitpb2->SetParameter(0,-1.0);
+   fitpb2->SetParameter(2,1.0);
+   fitpb2->SetParameter(3,25000.0);
+   fitpb2->SetParameter(4,-1.0);
+   fitpb2->SetParLimits(1,-0.02,0.02);
+   hproy->Fit(fitpb2,"RQM");
+   mean = fitpb2->GetParameter(1);
+   errMean = fitpb2->GetParError(1);
+   if(mean == 0)   errMean = 0;
+   
+}
+
+
+//____________________________________________________________________________________________________________________________________________
+// set variables 
+
+
+void AliPerfAnalyzeInvPt::SetProjBinsPhi(const Double_t *phiBinArray,Int_t nphBins){
+   
+   fNPhiBins = nphBins;
+  
+   for(Int_t k = 0;k<fNPhiBins;k++){
+      fPhiBins[k] = phiBinArray[k];
+   }
+   Printf("number of bins in phi set to %i",fNPhiBins);
+
+}
+//____________________________________________________________________________________________________________________________________________
+void AliPerfAnalyzeInvPt::SetProjBinsTheta(const Double_t *thetaBinArray, Int_t nthBins){
+  
+   fNThetaBins = nthBins;
+   for(Int_t k = 0;k<fNThetaBins;k++){
+      fThetaBins[k] = thetaBinArray[k];
+   }
+   Printf("number of bins in theta set to %i",fNThetaBins);
+
+}
+//____________________________________________________________________________________________________________________________________________
+void AliPerfAnalyzeInvPt::SetMakeFitOption(const Bool_t setGausFit, const Double_t exclusionR,const Double_t fitR ){
+
+  
+   
+   fFitGaus = setGausFit;
+   fExclRange  = exclusionR;
+   fRange = fitR;
+  
+   if(fFitGaus) Printf("set MakeGausFit with fit range %2.3f and exclusion range in 1/pt: %2.3f",fRange,fExclRange);
+   else  Printf("set standard polynomial fit with fit range %2.3f and exclusion range in 1/pt: %2.3f",fRange,fExclRange);
+}
diff --git a/PWG1/TPC/AliPerfAnalyzeInvPt.h b/PWG1/TPC/AliPerfAnalyzeInvPt.h
new file mode 100755 (executable)
index 0000000..afbab55
--- /dev/null
@@ -0,0 +1,73 @@
+#ifndef AliPerfAnalyzeInvPt_h
+#define AliPerfAnalyzeInvPt_h
+
+
+class TNamed;
+class TString;
+class TList;
+class TH1F;
+class TH2F;
+class TGraphErrors;
+
+class AliPerfAnalyzeInvPt : public TNamed {
+
+ public:
+  AliPerfAnalyzeInvPt();
+  AliPerfAnalyzeInvPt(Char_t* name, Char_t* title);
+   virtual ~AliPerfAnalyzeInvPt(){;}
+  void InitHistos(Double_t *binsXTheta,Double_t *fitParamTheta,Double_t *errFitParamTheta,Double_t *binsXPhi,Double_t *fitParamPhi,Double_t *errFitParamPhi);
+   void InitFitFcn();
+  void StartAnalysis(TH2F *histThetaInvPt, TH2F *histPhiInvPt, TObjArray* aFolderObj);
+  void StartAnalysis(TObjArray *aFolderObj);
+  void SetProjBinsPhi(const Double_t *pBins,Int_t sizep);
+  void SetProjBinsTheta(const Double_t *tBins, Int_t sizet);
+  void SetMakeFitOption(const Bool_t setGausFit, const Double_t exclusionR,const Double_t fitR );
+  
+   protected:
+  Double_t fThetaBins[100];
+  Double_t fPhiBins[100];
+  
+  Int_t fNThetaBins;
+  Int_t fNPhiBins ;
+  Double_t fRange;
+  Double_t fExclRange ;
+  Bool_t fFitGaus ;
+
+ // projection histos
+  TH1F *fHistFitTheta[100];
+  TH1F *fHistFitPhi[100];
+
+  TH2F *fHistH2InvPtTheta;
+  TH2F *fHistH2InvPtPhi; 
+  TGraphErrors *fGrMinPosTheta;
+  TGraphErrors *fGrMinPosPhi;
+
+
+ private:
+  TF1 *fFitMinPos;
+  TF1 *fFitMinPosRejP;
+  TF1 *fFitInvGauss;
+  TF1 *fFitInvGaussRejP;
+
+ static Double_t Polynomial(Double_t *x, Double_t *par);
+ static Double_t PolynomialRejP(Double_t *x, Double_t *par);
+ static Double_t InvGauss(Double_t *x, Double_t *par);
+ static Double_t InvGaussRejP(Double_t *x, Double_t *par);
+  void MakeFit(TH1 *dmproy, TF1 * fitpb, Double_t &mean, Double_t &ErrMean,  Double_t &excl,Double_t &range);
+  void MakeFitBetter(TH1 *dmproy, TF1 * fitpb2, Double_t &mean, Double_t &ErrMean, Double_t &f, Double_t &excl, Double_t &range);
+  void MakeFitInvGauss(TH1 *dmproy, TF1 * fitpb2, Double_t &mean, Double_t &ErrMean,Double_t &excl , Double_t &range);
+  void MakeFitInvGaussBetter(TH1 *dmproy, TF1 * fitpb2, Double_t &mean, Double_t &ErrMean, Double_t &f,Double_t &excl,  Double_t &range);
+
+  
+
+  AliPerfAnalyzeInvPt(const AliPerfAnalyzeInvPt&);            // not implemented 
+  AliPerfAnalyzeInvPt& operator=(const AliPerfAnalyzeInvPt&); // not implemented 
+
+  ClassDef(AliPerfAnalyzeInvPt, 1); 
+};
+
+#endif
diff --git a/PWG1/TPC/AliPerformancePtCalib.cxx b/PWG1/TPC/AliPerformancePtCalib.cxx
new file mode 100755 (executable)
index 0000000..dc7d1cd
--- /dev/null
@@ -0,0 +1,593 @@
+#include "TList.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TH3F.h"
+#include "TCanvas.h"
+#include "TMath.h"
+#include "TVector3.h"
+
+
+#include "AliESDVertex.h"
+#include "AliESDEvent.h"
+
+#include "AliESDtrack.h"
+#include "AliESDtrackCuts.h"
+
+
+#include "AliAnalysisTask.h"
+#include "AliAnalysisManager.h"
+#include "AliESDInputHandler.h"
+
+#include "AliPerfAnalyzeInvPt.h"
+#include "AliRecInfoCuts.h"
+#include "AliMCInfoCuts.h"
+
+#include "AliPerformancePtCalib.h"
+
+using namespace std;
+
+ClassImp(AliPerformancePtCalib)
+
+//________________________________________________________________________
+  AliPerformancePtCalib::AliPerformancePtCalib():
+  AliPerformanceObject("AliPerformancePtCalib"),
+  
+ // option parameter for Analyse()
+  fNThetaBins(0), 
+  fNPhiBins(0),
+  fRange(0),
+  fExclRange(0),
+  fFitGaus(0) ,
+  // option for user defined 1/pt shift
+  fShift(0),
+  fDeltaInvP(0),
+  
+  //options for cuts
+  fOptTPC(0),
+  fESDcuts(0),
+  fRefitTPC(0),
+  fRefitITS(0),
+  fDCAcut(0),
+
+  fMinPt(0),
+  fMaxPt(0),
+  fMinNClustersTPC(0),
+  fMaxChi2PerClusterTPC(0),
+  fMaxDCAtoVertexXY(0),
+  fMaxDCAtoVertexZ(0),
+
+  fCutsRC(0),
+  fCutsMC(0),
+
+  fList(0),
+  // histograms
+  fHistInvPtTheta(0),
+  fHistInvPtPhi(0),
+  fHistPtTheta(0),
+  fHistPtPhi(0),
+
+  fHistPtShift0(0),
+  fHistPrimaryVertexPosX(0),
+  fHistPrimaryVertexPosY(0),
+  fHistPrimaryVertexPosZ(0),
+  fHistTrackMultiplicity(0),
+  fHistTrackMultiplicityCuts(0),
+
+  fHistTPCMomentaPosP(0),
+  fHistTPCMomentaNegP(0),
+  fHistTPCMomentaPosPt(0),
+  fHistTPCMomentaNegPt(0),
+  fHistUserPtShift(0),
+
+  //esd track cuts
+  fESDTrackCuts(0),
+  
+  // analysis folder 
+  fAnalysisFolder(0)
+{
+   
+  // Dummy constructor
+
+  fShift = kFALSE;
+   fDeltaInvP = 0.00;
+  //options for cuts
+  fOptTPC =  kTRUE;                      // read TPC tracks yes/no
+  fESDcuts = kTRUE;                      // read ESD track cuts
+  fRefitTPC = kFALSE;                    // require TPC refit
+  fRefitITS = kFALSE;                    // require ITS refit
+  fDCAcut = kTRUE;                       // apply DCA cuts
+
+  fCutsRC = NULL;
+  fCutsMC = NULL;
+
+  fMinPt=0.15;  // GeV/c
+  fMaxPt=1.e10; // GeV/c 
+  fMinNClustersTPC = 50;
+  fMaxChi2PerClusterTPC = 4.0;
+  fMaxDCAtoVertexXY = 2.4; // cm
+  fMaxDCAtoVertexZ  = 3.2; // cm
+  
+ // options for function Analyse()
+  fFitGaus = kFALSE;
+   fNThetaBins = 0;
+   fNPhiBins =0  ;
+   fRange =0;
+   fExclRange =0;
+   fFitGaus =0;
+  
+  Init();
+}
+
+//________________________________________________________________________
+AliPerformancePtCalib::AliPerformancePtCalib(Char_t * name="AliPerformancePtCalib",Char_t* title ="AliPerformancePtCalib") :
+   AliPerformanceObject(name,title),
+   
+   // option parameter for Analyse()
+   fNThetaBins(0), 
+   fNPhiBins(0),
+   fRange(0),
+   fExclRange(0),
+   fFitGaus(0) ,
+  fShift(0),
+  fDeltaInvP(0),
+   
+   //options for cuts
+   fOptTPC(0),
+  fESDcuts(0),
+   fRefitTPC(0),
+  fRefitITS(0),
+  fDCAcut(0),
+
+  fMinPt(0),
+  fMaxPt(0),
+  fMinNClustersTPC(0),
+  fMaxChi2PerClusterTPC(0),
+  fMaxDCAtoVertexXY(0),
+  fMaxDCAtoVertexZ(0),
+
+  fCutsRC(0),
+  fCutsMC(0),
+   fList(0),
+   
+   // histograms
+  fHistInvPtTheta(0),
+  fHistInvPtPhi(0),
+  fHistPtTheta(0),
+  fHistPtPhi(0),
+
+  fHistPtShift0(0),
+  fHistPrimaryVertexPosX(0),
+  fHistPrimaryVertexPosY(0),
+  fHistPrimaryVertexPosZ(0),
+  fHistTrackMultiplicity(0),
+  fHistTrackMultiplicityCuts(0),
+
+  fHistTPCMomentaPosP(0),
+  fHistTPCMomentaNegP(0),
+  fHistTPCMomentaPosPt(0),
+  fHistTPCMomentaNegPt(0),
+   fHistUserPtShift(0),        
+   //esd track cuts                                                                                                                 
+  fESDTrackCuts(0),
+  
+  // analysis folder 
+  fAnalysisFolder(0)
+
+  
+{
+  // Constructor
+  fShift = kFALSE;
+  fDeltaInvP = 0.00;
+  //options for cuts
+  fOptTPC =  kTRUE;                      // read TPC tracks yes/no
+  fESDcuts = kTRUE;                      // read ESD track cuts
+  fRefitTPC = kFALSE;                    // require TPC refit
+  fRefitITS = kFALSE;                    // require ITS refit
+  fDCAcut = kTRUE;                       // apply DCA cuts
+  fCutsRC = NULL;
+  fCutsMC = NULL;
+
+  fMinPt=0.15;  // GeV/c
+  fMaxPt=1.e10; // GeV/c 
+  fMinNClustersTPC = 50;
+  fMaxChi2PerClusterTPC = 4.0;
+  fMaxDCAtoVertexXY = 2.4; // cm
+  fMaxDCAtoVertexZ  = 3.2; // cm
+
+  // options for Analyse()
+ fFitGaus = kFALSE;
+   fNThetaBins = 0;
+   fNPhiBins =0  ;
+   fRange =0;
+   fExclRange =0;
+   fFitGaus =0;
+  
+  Init();
+}
+
+AliPerformancePtCalib::~AliPerformancePtCalib() { 
+//
+// destructor
+//
+if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0; 
+}
+
+//________________________________________________________________________
+void AliPerformancePtCalib::Init() 
+{
+  // Create histograms
+  // Called once
+
+   fList = new TList();
+  // init folder
+  fAnalysisFolder = CreateFolder("folderPt_TPC","Analysis Pt Resolution Folder");
+  fList->Add(fAnalysisFolder);
+  // Primary Vertex:
+  fHistPrimaryVertexPosX       = new TH1F("fHistPrimaryVertexPosX", "Primary Vertex Position X;Primary Vertex Position X (cm);Events",100,-0.5,0.5);
+  fList->Add(fHistPrimaryVertexPosX);
+  fHistPrimaryVertexPosY       = new TH1F("fHistPrimaryVertexPosY", "Primary Vertex Position Y;Primary Vertex Position Y (cm);Events",100,-0.5,0.5);
+  fList->Add(fHistPrimaryVertexPosY);
+  fHistPrimaryVertexPosZ       = new TH1F("fHistPrimaryVertexPosZ", "Primary Vertex Position Z;Primary Vertex Position Z (cm);Events",200,-2.0,2.0);
+  fList->Add(fHistPrimaryVertexPosZ);
+  // Multiplicity:
+  fHistTrackMultiplicity     = new TH1F("fHistTrackMultiplicity", "Multiplicity distribution;Number of tracks;Events", 250, 0, 250);
+  fList->Add(fHistTrackMultiplicity);
+  fHistTrackMultiplicityCuts = new TH1F("fHistTrackMultiplicityCuts", "Multiplicity distribution;Number of tracks after cuts;Events", 250, 0, 250);
+  fList->Add(fHistTrackMultiplicityCuts);
+  // momentum histos
+  //pt shift 0 only needed if shift in 1/pt is applied
+  fHistPtShift0 = new TH1F("fHistPtShift0","1/pt dN/pt vs. pt of ESD track  ",600,0.0,6.0);
+  fList->Add(fHistPtShift0);
+  fHistInvPtTheta = new TH2F("fHistInvPtTheta","#theta vs 1/pt ",900, -4.5, 4.5,300,0.0,3.0);
+  fList->Add(fHistInvPtTheta);
+  fHistInvPtPhi   = new TH2F("fHistInvPtPhi","#phi vs 1/pt",900, -4.5, 4.5,325,0.0,6.5);
+  fList->Add(fHistInvPtPhi);
+  fHistPtTheta    = new TH2F("fHistPtTheta"," #theta vs pt ",300, 0.0, 15.0,300,0.0,3.0);
+  fList->Add(fHistPtTheta);
+  fHistPtPhi      = new TH2F("fHistPtPhi"," #phi vs  pt ",300, 0.0,15.0,325,0.0,6.5);
+  fList->Add(fHistPtPhi);
+            
+   // mom test histos
+  fHistTPCMomentaPosP  =  new TH2F("fHistTPCMomentaPosP","TPC p vs global esd track p pos",300,0.0,15.0,300,0.0,15.0);
+  fList->Add(fHistTPCMomentaPosP);
+  fHistTPCMomentaNegP  =  new TH2F("fHistTPCMomentaNegP","TPC p vs global esd track p neg",300,0.0,15.0,300,0.0,15.0);
+  fList->Add(fHistTPCMomentaNegP);
+  fHistTPCMomentaPosPt =  new TH2F("fHistTPCMomentaPosPt","TPC pt vs global esd track pt pos",300,0.0,15.0,300,0.0,15.0);
+  fList->Add(fHistTPCMomentaPosPt);
+  fHistTPCMomentaNegPt =  new TH2F("fHistTPCMomentaNegPt","TPC pt vs global esd track pt neg",300,0.0,15.0,300,0.0,15.0);
+  fList->Add(fHistTPCMomentaNegPt);
+
+  //user pt shift check
+   fHistUserPtShift = new TH1F("fHistUserPtShift","user defined shift in 1/pt",100,-0.5,1.5);
+   fList->Add(fHistUserPtShift);
+
+  
+  // esd track cuts  
+  fESDTrackCuts = new AliESDtrackCuts("AliESDtrackCuts");
+  
+  //fESDTrackCuts->DefineHistoqgrams(1);
+  fESDTrackCuts->SetRequireSigmaToVertex(kFALSE);
+  fESDTrackCuts->SetRequireTPCRefit(kTRUE);
+  fESDTrackCuts->SetAcceptKinkDaughters(kTRUE);
+  fESDTrackCuts->SetMinNClustersTPC((Int_t)fMinNClustersTPC);
+  fESDTrackCuts->SetMaxChi2PerClusterTPC(fMaxChi2PerClusterTPC);
+  fESDTrackCuts->SetMaxDCAToVertexXY(fMaxDCAtoVertexXY);
+  fESDTrackCuts->SetMaxDCAToVertexZ(fMaxDCAtoVertexZ);
+  fESDTrackCuts->SetDCAToVertex2D(kTRUE);
+  fESDTrackCuts->SetPtRange(fMinPt,fMaxPt); 
+}
+
+//________________________________________________________________________
+void AliPerformancePtCalib::SetPtShift(Double_t shiftVal ) { 
+   if(!(shiftVal==0)) { fShift=kTRUE; fDeltaInvP = shiftVal; } 
+}
+
+//________________________________________________________________________
+void AliPerformancePtCalib::Exec(AliMCEvent*, AliESDEvent* const esdEvent, AliESDfriend*, Bool_t, Bool_t) 
+{
+
+  if (!esdEvent) {
+    Printf("ERROR: Event not available");
+    return;
+  }
+
+  if (!(esdEvent->GetNumberOfTracks())) {
+    Printf(" PtCalib task: There is no track in this event");
+    return;
+  }
+
+  if(fShift) fHistUserPtShift->Fill(fDeltaInvP);
+  
+  fHistTrackMultiplicity->Fill(esdEvent->GetNumberOfTracks());
+
+
+  // read primary vertex info
+  Double_t tPrimaryVtxPosition[3];
+  //  Double_t tPrimaryVtxCov[3];
+  const AliESDVertex *primaryVtx = esdEvent->GetPrimaryVertexTPC();
+  tPrimaryVtxPosition[0] = primaryVtx->GetXv();
+  tPrimaryVtxPosition[1] = primaryVtx->GetYv();
+  tPrimaryVtxPosition[2] = primaryVtx->GetZv();
+  
+  fHistPrimaryVertexPosX->Fill(tPrimaryVtxPosition[0]);
+  fHistPrimaryVertexPosY->Fill(tPrimaryVtxPosition[1]);
+  fHistPrimaryVertexPosZ->Fill(tPrimaryVtxPosition[2]);
+
+
+  //_fill histos for pt spectra and shift of transverse momentum
+  Int_t count=0;
+  for(Int_t j = 0;j<esdEvent->GetNumberOfTracks();j++){// track loop
+    AliESDtrack *ESDTrack = esdEvent->GetTrack(j);
+    if(!ESDTrack) continue;
+    
+    
+    if(fESDcuts == kTRUE){
+        if(!fESDTrackCuts->AcceptTrack(ESDTrack))continue;
+    }
+       
+    //track cuts
+    if(fRefitTPC) if(AddTPCcuts(ESDTrack)) continue;
+    if(fRefitITS) if(AddITScuts(ESDTrack)) continue;
+    if(fDCAcut)   if(AddDCAcuts(ESDTrack)) continue ;
+      
+    // fill histos
+    if(fOptTPC == kTRUE){ //TPC tracks
+      const AliExternalTrackParam *TPCTrack = ESDTrack->GetTPCInnerParam(); 
+      if(!TPCTrack) continue;
+      if(fabs(TPCTrack->Eta())> 0.8) continue;
+      
+      Double_t signedPt = TPCTrack->GetSignedPt();
+      Double_t invPt = 0.0;
+      if(!signedPt==0) {
+       invPt = 1.0/signedPt;
+
+       fHistPtShift0->Fill(fabs(signedPt));
+       
+       if(fShift == kTRUE ){
+         invPt += fDeltaInvP; //shift momentum for tests
+         if(!invPt==0) signedPt = 1.0/invPt;
+         else continue;
+       }
+
+       fHistInvPtTheta->Fill(invPt,TPCTrack->Theta());
+       fHistInvPtPhi->Fill(invPt,TPCTrack->Phi());
+       fHistPtTheta->Fill(fabs(signedPt),TPCTrack->Theta());
+       fHistPtPhi->Fill(fabs(signedPt),TPCTrack->Phi());
+
+
+       Double_t pTPC  = TPCTrack->GetP();
+       Double_t pESD = ESDTrack->GetP();
+       Double_t ptESD  = ESDTrack->GetSignedPt();
+       
+       if(ESDTrack->GetSign()>0){//compare momenta ESD track and TPC track
+         fHistTPCMomentaPosP->Fill(fabs(pESD),fabs(pTPC));
+         fHistTPCMomentaPosPt->Fill(fabs(ptESD),fabs(signedPt));
+       }
+       else{
+         fHistTPCMomentaNegP->Fill(fabs(pESD),fabs(pTPC));
+         fHistTPCMomentaNegPt->Fill(fabs(ptESD),fabs(signedPt));
+       }
+       count++;
+      }
+      else continue;
+    }
+   
+    else{// ESD tracks
+      Double_t invPt = 0.0;
+      Double_t signedPt = ESDTrack->GetSignedPt();
+      if(!signedPt==0){
+       invPt = 1.0/signedPt; 
+
+       fHistPtShift0->Fill(fabs(signedPt));
+         
+       if(fShift == kTRUE ){
+         invPt += fDeltaInvP;//shift momentum for tests
+         if(!invPt==0) signedPt = 1.0/invPt;
+         else continue;
+       }
+       fHistInvPtTheta->Fill(invPt,ESDTrack->Theta());
+       fHistInvPtPhi->Fill(invPt,ESDTrack->Phi());
+       fHistPtTheta->Fill(signedPt,ESDTrack->Theta());
+       fHistPtPhi->Fill(signedPt,ESDTrack->Phi());
+       
+       count++;
+      }
+    }
+  }
+    
+  fHistTrackMultiplicityCuts->Fill(count);
+    
+}    
+
+
+
+//______________________________________________________________________________________________________________________
+Bool_t AliPerformancePtCalib::AddTPCcuts(AliESDtrack *ESDTrack){
+  
+  Bool_t cut = kFALSE;
+  
+  if ((ESDTrack->GetStatus()&AliESDtrack::kTPCrefit)==0) cut=kTRUE; // TPC refit
+  if (ESDTrack->GetTPCNcls()<50) cut=kTRUE; // min. nb. TPC clusters
+  if(cut) return kTRUE;
+  return kFALSE;
+}
+//______________________________________________________________________________________________________________________
+Bool_t AliPerformancePtCalib::AddDCAcuts(AliESDtrack *ESDTrack){
+  
+  Bool_t cut = kFALSE;
+  
+  Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z and impact parameters:
+  ESDTrack->GetImpactParameters(dca,cov);
+  if(TMath::Abs(dca[0])>3. || TMath::Abs(dca[1])>3.) cut=kTRUE;
+  if(ESDTrack->GetKinkIndex(0)>0) cut=kTRUE;
+  if(cut) return kTRUE;
+  return kFALSE;
+}
+
+//______________________________________________________________________________________________________________________
+Bool_t AliPerformancePtCalib::AddITScuts(AliESDtrack *ESDTrack){
+
+  Bool_t cut = kFALSE;
+  
+  if ((ESDTrack->GetStatus()&AliESDtrack::kITSrefit)==0) cut=kTRUE; // ITS refit
+  Int_t clusterITS[200]; 
+  if(ESDTrack->GetITSclusters(clusterITS)<2) cut=kTRUE;  // min. nb. ITS clusters //3
+  
+  if(cut) return kTRUE;
+  return kFALSE;
+}
+
+//______________________________________________________________________________________________________________________
+
+void AliPerformancePtCalib::Analyse()
+{
+  
+  
+  AliPerfAnalyzeInvPt *ana = new  AliPerfAnalyzeInvPt("AliPerfAnalyzeInvPt","AliPerfAnalyzeInvPt");
+  
+     
+  TH1::AddDirectory(kFALSE);
+  ana->SetProjBinsTheta(fThetaBins,fNThetaBins);
+  ana->SetProjBinsPhi(fPhiBins,fNPhiBins);
+  ana->SetMakeFitOption(fFitGaus,fExclRange,fRange);
+  
+  TObjArray *aFolderObj = new TObjArray;
+  ana->StartAnalysis(fHistInvPtTheta,fHistInvPtPhi, aFolderObj);
+  
+  // export objects to analysis folder
+  fAnalysisFolder = ExportToFolder(aFolderObj);
+
+  // delete only TObjArray
+  if(aFolderObj) delete aFolderObj;
+  if(ana) delete ana;
+  
+}
+
+//______________________________________________________________________________________________________________________
+TFolder* AliPerformancePtCalib::ExportToFolder(TObjArray * array) 
+{
+  // recreate folder every time and export objects to new one
+  //
+  AliPerformancePtCalib * comp=this;
+  TFolder *folder = comp->GetAnalysisFolder();
+
+  TString name, title;
+  TFolder *newFolder = 0;
+  Int_t i = 0;
+  Int_t size = array->GetSize();
+
+  if(folder) { 
+    // get name and title from old folder
+    name = folder->GetName();  
+    title = folder->GetTitle();  
+
+    // delete old one
+    delete folder;
+
+    // create new one
+    newFolder = CreateFolder(name.Data(),title.Data());
+    newFolder->SetOwner();
+
+    // add objects to folder
+    while(i < size) {
+      newFolder->Add(array->At(i));
+      i++;
+    }
+  }
+
+  return newFolder;
+}
+
+//______________________________________________________________________________________________________________________
+Long64_t AliPerformancePtCalib::Merge(TCollection* const list) 
+{
+  // Merge list of objects (needed by PROOF)
+
+  if (!list)
+    return 0;
+
+  if (list->IsEmpty())
+    return 1;
+
+  TIterator* iter = list->MakeIterator();
+  TObject* obj = 0;
+
+  // collection of generated histograms
+  Int_t count=0;
+  while((obj = iter->Next()) != 0) 
+    {
+      AliPerformancePtCalib* entry = dynamic_cast<AliPerformancePtCalib*>(obj);
+      if (entry == 0) continue; 
+  
+      fHistInvPtTheta->Add(entry->fHistInvPtTheta);
+      fHistInvPtPhi->Add(entry-> fHistInvPtPhi);
+      fHistPtTheta->Add(entry->fHistPtTheta);
+      fHistPtPhi->Add(entry->fHistPtPhi);
+  
+      fHistPtShift0->Add(entry->fHistPtShift0);
+      fHistPrimaryVertexPosX->Add(entry->fHistPrimaryVertexPosX);
+      fHistPrimaryVertexPosY->Add(entry->fHistPrimaryVertexPosY);
+      fHistPrimaryVertexPosZ->Add(entry->fHistPrimaryVertexPosZ);
+      fHistTrackMultiplicity->Add(entry->fHistTrackMultiplicity);
+      fHistTrackMultiplicityCuts->Add(entry->fHistTrackMultiplicityCuts);
+  
+      fHistTPCMomentaPosP->Add(entry->fHistTPCMomentaPosP);
+      fHistTPCMomentaNegP->Add(entry->fHistTPCMomentaNegP);
+      fHistTPCMomentaPosPt->Add(entry->fHistTPCMomentaPosPt);
+      fHistTPCMomentaNegPt->Add(entry->fHistTPCMomentaNegPt);
+  
+      count++;
+    }
+  
+  return count;
+}
+
+//______________________________________________________________________________________________________________________
+TFolder* AliPerformancePtCalib::CreateFolder(TString name,TString title) { 
+  // create folder for analysed histograms
+  //
+  TFolder *folder = 0;
+  folder = new TFolder(name.Data(),title.Data());
+
+  return folder;
+}
+//______________________________________________________________________________________________________________________
+void AliPerformancePtCalib::SetProjBinsPhi(const Double_t *phiBinArray,Int_t nphBins){
+   
+   fNPhiBins = nphBins;
+  
+   for(Int_t k = 0;k<fNPhiBins;k++){
+      fPhiBins[k] = phiBinArray[k];
+   }
+   Printf("AliPerformancePtCalib: number of bins in phi set to %i",fNPhiBins);
+
+}
+//____________________________________________________________________________________________________________________________________________
+void AliPerformancePtCalib::SetProjBinsTheta(const Double_t *thetaBinArray, Int_t nthBins){
+  
+   fNThetaBins = nthBins;
+   for(Int_t k = 0;k<fNThetaBins;k++){
+      fThetaBins[k] = thetaBinArray[k];
+   }
+   Printf("AliPerformancePtCalib: number of bins in theta set to %i",fNThetaBins);
+
+}
+//____________________________________________________________________________________________________________________________________________
+void AliPerformancePtCalib::SetMakeFitOption(const Bool_t setGausFit, const Double_t exclusionR,const Double_t fitR ){
+
+  
+   fRange = fitR;
+   fFitGaus = setGausFit;
+   fExclRange  = exclusionR;
+  
+   if(fFitGaus) Printf("AliPerformancePtCalib: set MakeGausFit with fit range %2.3f and exclusion range in 1/pt: %2.3f",fRange,fExclRange);
+   else  Printf("AliPerformancePtCalib: set standard polynomial fit with fit range %2.3f and exclusion range in 1/pt: %2.3f",fRange,fExclRange);
+}
diff --git a/PWG1/TPC/AliPerformancePtCalib.h b/PWG1/TPC/AliPerformancePtCalib.h
new file mode 100755 (executable)
index 0000000..6ff08cc
--- /dev/null
@@ -0,0 +1,148 @@
+#ifndef ALIPERFORMANCEPTCALIB_H
+#define ALIPERFORMANCEPTCALIB_H
+
+class TString;
+class TNamed;
+class TCanvas;
+class TH1F;
+class TH2F;
+class TList;
+
+class AliESDVertex;
+class AliESDtrack;
+class AliMCEvent;
+class AliStack;
+class AliTrackReference;
+class AliESDEvent; 
+class AliESDfriend; 
+class AliESDfriendTrack; 
+class AliMCEvent;
+class AliMCParticle;
+class AliMCInfoCuts;
+class AliRecInfoCuts;
+
+#include "AliPerformanceObject.h"
+
+class AliPerformancePtCalib : public AliPerformanceObject {
+ public:
+  AliPerformancePtCalib(); 
+  AliPerformancePtCalib(Char_t* name, Char_t* title);
+  virtual ~AliPerformancePtCalib();
+
+  // Init data members
+  virtual void  Init();
+
+  // Execute analysis
+  virtual void  Exec(AliMCEvent* const mcEvent, AliESDEvent *const esdEvent, AliESDfriend *const esdFriend, const Bool_t bUseMC, const Bool_t bUseESDfriend);
+
+  // Merge output objects (needed by PROOF) 
+  virtual Long64_t Merge(TCollection* const list);
+
+  // Analyse output histograms
+  virtual void Analyse();
+
+  // Get analysis folder
+  virtual TFolder* GetAnalysisFolder() const {return fAnalysisFolder;}
+   
+  Bool_t AddTPCcuts(AliESDtrack *ESDTrack);
+  Bool_t AddITScuts(AliESDtrack *ESDTrack);
+  Bool_t AddDCAcuts(AliESDtrack *ESDTrack);
+  void SetReadTPCTracksOff(){fOptTPC   = kFALSE;}// function  can be used in the macro to set the variables
+  void SetTPCRefit()        {fRefitTPC = kTRUE;} // function  can be used in the macro to set the variables
+  void SetITSRefit()        {fRefitITS = kTRUE;} // function  can be used in the macro to set the variables
+  void SetESDcutsOff()      {fESDcuts  = kFALSE;}// function  can be used in the macro to set the variables
+  void SetDCAcutsOff()      {fDCAcut   = kFALSE;}// function  can be used in the macro to set the variables
+  
+  //void SetPtShift(Double_t shiftVal ) { if(!(shiftVal==0)) {fShift=kTRUE; fDeltaInvP = shiftVal;} }
+  void SetPtShift(Double_t shiftVal);
+
+  void SetESDcutValues(Double_t * esdCutValues){// functions can be used in the macro to set the variables
+    fMinPt                = esdCutValues[0]; 
+    fMaxPt                = esdCutValues[1];
+    fMinNClustersTPC      = esdCutValues[2];
+    fMaxChi2PerClusterTPC = esdCutValues[3];
+    fMaxDCAtoVertexXY     = esdCutValues[4];
+    fMaxDCAtoVertexZ      = esdCutValues[5];
+  }
+   void SetProjBinsPhi(const Double_t *pBins,Int_t sizep);
+   void SetProjBinsTheta(const Double_t *tBins, Int_t sizet);
+   void SetMakeFitOption(const Bool_t setGausFit, const Double_t exclusionR,const Double_t fitR );
+   TList *GetHistoList() {return fList;}
+   
+  // Create folder for analysed histograms
+  TFolder *CreateFolder(TString folder = "folderPtCalib",TString title = "Analysed PtCalib histograms");
+
+  // Export objects to folder
+  TFolder *ExportToFolder(TObjArray * array=0);
+
+  // Selection cuts
+  void SetAliRecInfoCuts(AliRecInfoCuts* const cuts=0) {fCutsRC = cuts;}   
+  void SetAliMCInfoCuts(AliMCInfoCuts* const cuts=0) {fCutsMC = cuts;}  
+   
+  AliRecInfoCuts*  GetAliRecInfoCuts() const {return fCutsRC;}  
+  AliMCInfoCuts*   GetAliMCInfoCuts()  const {return fCutsMC;}
+
+protected:
+   
+  Double_t fThetaBins[100];
+  Double_t fPhiBins[100];
+   
+    Int_t fNThetaBins;
+  Int_t fNPhiBins ;
+  Double_t fRange;
+  Double_t fExclRange ;
+  Bool_t fFitGaus ;
+
+ private:
+
+  Bool_t fShift;
+  Double_t fDeltaInvP;
+  //options for cuts
+  Bool_t fOptTPC;
+  Bool_t fESDcuts;
+  Bool_t fRefitTPC;
+  Bool_t fRefitITS;
+  Bool_t fDCAcut;
+
+  Double_t fMinPt;
+  Double_t fMaxPt;
+  Double_t fMinNClustersTPC;
+  Double_t fMaxChi2PerClusterTPC;
+  Double_t fMaxDCAtoVertexXY;
+  Double_t fMaxDCAtoVertexZ;
+
+  AliRecInfoCuts* fCutsRC;     // selection cuts for reconstructed tracks
+  AliMCInfoCuts*  fCutsMC;     // selection cuts for MC tracks
+
+  TList       *fList;
+  TH2F        *fHistInvPtTheta;
+  TH2F        *fHistInvPtPhi;
+  TH2F        *fHistPtTheta;
+  TH2F        *fHistPtPhi;
+
+  TH1F        *fHistPtShift0;
+  TH1F        *fHistPrimaryVertexPosX;          
+  TH1F        *fHistPrimaryVertexPosY;          
+  TH1F        *fHistPrimaryVertexPosZ;          
+  TH1F        *fHistTrackMultiplicity;          
+  TH1F        *fHistTrackMultiplicityCuts;
+
+  TH2F        *fHistTPCMomentaPosP;
+  TH2F        *fHistTPCMomentaNegP;
+  TH2F        *fHistTPCMomentaPosPt;
+  TH2F        *fHistTPCMomentaNegPt;
+   TH1F        *fHistUserPtShift;
+   
+  AliESDtrackCuts* fESDTrackCuts;
+  
+  // analysis folder 
+  TFolder *fAnalysisFolder; // folder for analysed histograms
+
+  AliPerformancePtCalib(const AliPerformancePtCalib&);            // not implemented 
+  AliPerformancePtCalib& operator=(const AliPerformancePtCalib&); // not implemented 
+
+  ClassDef(AliPerformancePtCalib, 1); 
+};
+
+#endif
diff --git a/PWG1/TPC/AliPerformancePtCalibMC.cxx b/PWG1/TPC/AliPerformancePtCalibMC.cxx
new file mode 100755 (executable)
index 0000000..652b747
--- /dev/null
@@ -0,0 +1,754 @@
+#include "TList.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TH3F.h"
+#include "TCanvas.h"
+#include "TMath.h"
+#include "TVector3.h"
+
+
+#include "AliESDVertex.h"
+#include "AliESDEvent.h"
+
+#include "AliESDtrack.h"
+#include "AliESDtrackCuts.h"
+
+
+#include "AliAnalysisTask.h"
+#include "AliAnalysisManager.h"
+#include "AliESDInputHandler.h"
+#include "AliMCEventHandler.h"
+#include "AliMCEvent.h"
+#include "AliStack.h"
+
+
+#include "AliPerformancePtCalibMC.h"
+#include "AliPerfAnalyzeInvPt.h"
+#include "AliRecInfoCuts.h"
+#include "AliMCInfoCuts.h"
+
+using namespace std;
+
+ClassImp(AliPerformancePtCalibMC)
+
+//________________________________________________________________________
+   AliPerformancePtCalibMC::AliPerformancePtCalibMC() :
+      AliPerformanceObject("AliPerformancePtCalibMC"),
+      // option parameter for Analyse()
+      fNThetaBins(0), 
+      fNPhiBins(0),
+      fRange(0),
+      fExclRange(0),
+      fFitGaus(0) ,
+      fAnaMC(0),
+      // option for user defined 1/pt shift
+      fShift(0),
+      fDeltaInvP (0),
+      //options for cuts
+      fOptTPC(0),
+      fESDcuts(0),
+      fRefitTPC(0),
+      fRefitITS(0),
+      fDCAcut(0),
+      fMinPt(0),
+      fMaxPt(0),
+      fMinNClustersTPC(0),
+      fMaxChi2PerClusterTPC(0),
+      fMaxDCAtoVertexXY(0),
+      fMaxDCAtoVertexZ(0),
+      fCutsRC(0),
+      fCutsMC(0),
+      fList(0),
+      // histograms
+      fHistInvPtTheta(0),
+      fHistInvPtPhi(0),
+      fHistPtTheta(0),
+      fHistPtPhi(0),
+      fHistPtShift0(0),
+      fHistPrimaryVertexPosX(0),
+      fHistPrimaryVertexPosY(0),
+      fHistPrimaryVertexPosZ(0),
+      fHistTrackMultiplicity(0),
+      fHistTrackMultiplicityCuts(0),
+      fHistTPCMomentaPosP(0),
+      fHistTPCMomentaNegP(0),
+      fHistTPCMomentaPosPt(0),
+      fHistTPCMomentaNegPt(0),
+      fHistInvPtThetaMC(0),
+      fHistInvPtPhiMC(0),
+      fHistPtThetaMC(0),
+      fHistPtPhiMC(0),
+      fHistInvPtMCESD(0),
+      fHistInvPtMCTPC(0),
+      fHistPtMCESD(0),
+      fHistPtMCTPC(0),
+      fHistMomresMCESD(0),
+      fHistMomresMCTPC(0),
+      fHistTPCMomentaPosInvPtMC(0),
+      fHistTPCMomentaNegInvPtMC(0),
+      fHistTPCMomentaPosPtMC(0),
+      fHistTPCMomentaNegPtMC(0),
+      fHistESDMomentaPosInvPtMC(0),
+      fHistESDMomentaNegInvPtMC(0),
+      fHistESDMomentaPosPtMC(0), 
+      fHistESDMomentaNegPtMC(0),
+      fHistUserPtShift(0),
+      fESDTrackCuts(0),
+      // analysis folder 
+      fAnalysisFolder(0)
+{
+   // Dummy constructor
+   
+   
+   fShift = kFALSE;
+   fDeltaInvP = 0.00;
+   //options for cuts
+   fOptTPC =  kTRUE;                      // read TPC tracks yes/no
+   fESDcuts = kTRUE;                      // read ESD track cuts
+   fRefitTPC = kFALSE;                    // require TPC refit
+   fRefitITS = kFALSE;                    // require ITS refit
+   fDCAcut = kTRUE;                       // apply DCA cuts
+   
+   fCutsRC = NULL;
+   fCutsMC = NULL;
+   
+   fMinPt=0.15;  // GeV/c
+   fMaxPt=1.e10; // GeV/c 
+   fMinNClustersTPC = 50;
+   fMaxChi2PerClusterTPC = 4.0;
+   fMaxDCAtoVertexXY = 2.4; // cm
+   fMaxDCAtoVertexZ  = 3.0; // cm
+
+   // options for function Analyse()
+   fFitGaus = kFALSE;
+   fNThetaBins = 0;
+   fNPhiBins =0  ;
+   fRange =0;
+   fExclRange =0;
+   fFitGaus =0;
+   fAnaMC = kTRUE;
+   
+   Init();
+} 
+
+//________________________________________________________________________
+AliPerformancePtCalibMC::AliPerformancePtCalibMC(const char *name= "AliPerformancePtCalibMC", const char *title="AliPerformancePtCalibMC") :
+   AliPerformanceObject(name,title),
+   // option parameter for Analyse()
+   fNThetaBins(0), 
+   fNPhiBins(0),
+   fRange(0),
+   fExclRange(0),
+   fFitGaus(0) ,
+   fAnaMC(0),
+   // option for user defined 1/pt shift
+   fShift(0),
+   fDeltaInvP (0),
+   //options for cuts
+   fOptTPC(0),
+   fESDcuts(0),
+   fRefitTPC(0),
+   fRefitITS(0),
+   fDCAcut(0),
+   fMinPt(0),
+   fMaxPt(0),
+   fMinNClustersTPC(0),
+   fMaxChi2PerClusterTPC(0),
+   fMaxDCAtoVertexXY(0),
+   fMaxDCAtoVertexZ(0),
+   fCutsRC(0),
+   fCutsMC(0),
+   fList(0),
+   // histograms
+   fHistInvPtTheta(0),
+   fHistInvPtPhi(0),
+   fHistPtTheta(0),
+   fHistPtPhi(0),
+   fHistPtShift0(0),
+   fHistPrimaryVertexPosX(0),
+   fHistPrimaryVertexPosY(0),
+   fHistPrimaryVertexPosZ(0),
+   fHistTrackMultiplicity(0),
+   fHistTrackMultiplicityCuts(0),
+   fHistTPCMomentaPosP(0),
+   fHistTPCMomentaNegP(0),
+   fHistTPCMomentaPosPt(0),
+   fHistTPCMomentaNegPt(0),
+   fHistInvPtThetaMC(0),
+   fHistInvPtPhiMC(0),
+   fHistPtThetaMC(0),
+   fHistPtPhiMC(0),
+   fHistInvPtMCESD(0),
+   fHistInvPtMCTPC(0),
+   fHistPtMCESD(0),
+   fHistPtMCTPC(0),
+   fHistMomresMCESD(0),
+   fHistMomresMCTPC(0),
+   fHistTPCMomentaPosInvPtMC(0),
+   fHistTPCMomentaNegInvPtMC(0),
+   fHistTPCMomentaPosPtMC(0),
+   fHistTPCMomentaNegPtMC(0),
+   fHistESDMomentaPosInvPtMC(0),
+   fHistESDMomentaNegInvPtMC(0),
+   fHistESDMomentaPosPtMC(0), 
+   fHistESDMomentaNegPtMC(0),
+   fHistUserPtShift(0),
+   
+   fESDTrackCuts(0),
+   // analysis folder 
+   fAnalysisFolder(0)
+   
+{
+   // Constructor
+    
+   fShift = kFALSE;
+   fDeltaInvP = 0.00;
+   //options for cuts
+   fOptTPC =  kTRUE;                      // read TPC tracks yes/no
+   fESDcuts = kTRUE;                      // read ESD track cuts
+   fRefitTPC = kFALSE;                    // require TPC refit
+   fRefitITS = kFALSE;                    // require ITS refit
+   fDCAcut = kTRUE;                       // apply DCA cuts
+   fCutsRC = NULL;
+   fCutsMC = NULL;
+
+   fMinPt=0.15;  // GeV/c
+   fMaxPt=1.e10; // GeV/c 
+   fMinNClustersTPC = 50;
+   fMaxChi2PerClusterTPC = 4.0;
+   fMaxDCAtoVertexXY = 2.4; // cm
+   fMaxDCAtoVertexZ  = 3.0; // cm
+
+   // options for function Analyse()
+   fFitGaus = kFALSE;
+   fNThetaBins = 0;
+   fNPhiBins =0  ;
+   fRange =0;
+   fExclRange =0;
+   fFitGaus =0;
+   fAnaMC = kTRUE;
+
+   Init();
+}
+
+//________________________________________________________________________
+AliPerformancePtCalibMC::~AliPerformancePtCalibMC() { 
+//
+// destructor
+//
+if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0; 
+}
+
+//________________________________________________________________________
+void AliPerformancePtCalibMC::Init() 
+{
+   // Create histograms
+   // Called once
+   
+   fList = new TList();
+   // init folder
+   fAnalysisFolder = CreateFolder("folderPt_TPC","Analysis Pt Resolution Folder");
+   fList->Add(fAnalysisFolder);
+   // Primary Vertex:
+   fHistPrimaryVertexPosX       = new TH1F("fHistPrimaryVertexPosX", "Primary Vertex Position X;Primary Vertex Position X (cm);Events",100,-0.5,0.5);
+   fList->Add(fHistPrimaryVertexPosX);
+   fHistPrimaryVertexPosY       = new TH1F("fHistPrimaryVertexPosY", "Primary Vertex Position Y;Primary Vertex Position Y (cm);Events",100,-0.5,0.5);
+   fList->Add(fHistPrimaryVertexPosY);
+   fHistPrimaryVertexPosZ       = new TH1F("fHistPrimaryVertexPosZ", "Primary Vertex Position Z;Primary Vertex Position Z (cm);Events",200,-2.0,2.0);
+   fList->Add(fHistPrimaryVertexPosZ);
+  
+   // Multiplicity:
+   fHistTrackMultiplicity     = new TH1F("fHistTrackMultiplicity", "Multiplicity distribution;Number of tracks;Events", 250, 0, 250);
+   fList->Add(fHistTrackMultiplicity);
+   fHistTrackMultiplicityCuts = new TH1F("fHistTrackMultiplicityCuts", "Multiplicity distribution;Number of tracks after cuts;Events", 250, 0, 250);
+   fList->Add(fHistTrackMultiplicityCuts);
+  
+   // momentum histos
+   fHistPtShift0   = new TH1F("fHistPtShift0","1/pt dN/pt vs. pt of ESD track  ",600,0.0,6.0);
+   fList->Add(fHistPtShift0);
+   fHistInvPtTheta = new TH2F("fHistInvPtTheta","#theta vs 1/pt ",900, -4.5, 4.5,300,0.0,3.0);
+   fList->Add(fHistInvPtTheta);
+   fHistInvPtPhi   = new TH2F("fHistInvPtPhi","#phi vs 1/pt",900, -4.5, 4.5,325,0.0,6.5);
+   fList->Add(fHistInvPtPhi);
+   fHistPtTheta    = new TH2F("fHistPtTheta"," #theta vs pt ",300, 0.0, 15.0,300,0.0,3.0);
+   fList->Add(fHistPtTheta);
+   fHistPtPhi      = new TH2F("fHistPtPhi"," #phi vs  pt ",300, 0.0,15.0,325,0.0,6.5);
+   fList->Add(fHistPtPhi);
+  
+   // mom test histos
+   fHistTPCMomentaPosP  =  new TH2F("fHistTPCMomentaPosP","TPC p vs global esd track p pos",300,0.0,15.0,300,0.0,15.0);
+   fList->Add(fHistTPCMomentaPosP);
+   fHistTPCMomentaNegP  =  new TH2F("fHistTPCMomentaNegP","TPC p vs global esd track p neg",300,0.0,15.0,300,0.0,15.0);
+   fList->Add(fHistTPCMomentaNegP);
+   fHistTPCMomentaPosPt =  new TH2F("fHistTPCMomentaPosPt","TPC pt vs global esd track pt pos",300,0.0,15.0,300,0.0,15.0);
+   fList->Add(fHistTPCMomentaPosPt);
+   fHistTPCMomentaNegPt =  new TH2F("fHistTPCMomentaNegPt","TPC pt vs global esd track pt neg",300,0.0,15.0,300,0.0,15.0);
+   fList->Add(fHistTPCMomentaNegPt);
+   
+   // mom test histos MC
+   fHistTPCMomentaPosInvPtMC = new TH2F("fHistTPCMomentaPosInvPtMC","TPC-MC of 1/pt vs global ESD-MC of 1/pt pos",500, -10.0, 10.0,500, -10.0,10.0);
+   fList->Add(fHistTPCMomentaPosInvPtMC);
+   fHistTPCMomentaNegInvPtMC = new TH2F("fHistTPCMomentaNegInvPtMC","TPC-MC of 1/pt vs global ESD-MC 1/pt neg",500, -10.0, 10.0,500, -10.0, 10.0);
+   fList->Add(fHistTPCMomentaNegInvPtMC);
+   fHistTPCMomentaPosPtMC    = new TH2F("fHistTPCMomentaPosPtMC","TPC-MC of pt vs global ESD-MC of pt pos",600,-4.0,44.0,600,-4.0,44.0);
+   fList->Add(fHistTPCMomentaPosPtMC);
+   fHistTPCMomentaNegPtMC    = new TH2F("fHistTPCMomentaNegPtMC","TPC-MC of pt vs global ESD-MC of pt neg",600,-4.0,44.0,600,-4.0,44.0);
+   fList->Add(fHistTPCMomentaNegPtMC);
+   fHistESDMomentaPosInvPtMC = new TH1F("fHistESDMomentaPosInvPtMC","ESD-MC of 1/pt ",500, -10.0, 10.0);
+   fList->Add(fHistESDMomentaPosInvPtMC);
+   fHistESDMomentaNegInvPtMC = new TH1F("fHistESDMomentaNegInvPtMC","ESD-MC of 1/pt",500, -10.0, 10.0);
+   fList->Add(fHistESDMomentaNegInvPtMC);
+   fHistESDMomentaPosPtMC    = new TH1F("fHistESDMomentaPosPtMC","ESD-MC of pt ",600,-4.0,44.0);
+   fList->Add(fHistESDMomentaPosPtMC);
+   fHistESDMomentaNegPtMC    = new TH1F("fHistESDMomentaNegPtMC","ESD-MC of pt ",600,-4.0,44.0);
+   fList->Add(fHistESDMomentaNegPtMC);
+
+   // MC info
+   fHistInvPtThetaMC = new TH2F("fHistInvPtThetaMC","theta vs inv pt MC",900, -4.5, 4.5,300,0.0,3.0);
+   fList->Add(fHistInvPtThetaMC);
+   fHistInvPtPhiMC   = new TH2F("fHistInvPtPhiMC","phi vs inv pt MC",900, -4.5, 4.5,325,0.0,6.5);
+   fList->Add(fHistInvPtPhiMC);
+   fHistPtThetaMC    = new TH2F("fHistPtThetaMC","theta vs pt MC",300, 0.0, 15.0,300,0.0,3.0);
+   fList->Add(fHistPtThetaMC);
+   fHistPtPhiMC      = new TH2F("fHistPtPhiMC"," phi vs pt MC",300, 0.0,15.0,325,0.0,6.5);
+   fList->Add(fHistPtPhiMC);
+   //correlation histos MC ESD or TPC
+   fHistInvPtMCESD  = new TH2F("fHistInvPtMCESD","inv pt ESD vs MC",900, 0.0, 9.0,900, 0.0, 9.0);
+   fList->Add(fHistInvPtMCESD);
+   fHistPtMCESD     = new TH2F("fHistPtMCESD"," pt ESD vs MC",300, 0.0, 15.0,300, 0.0, 15.0);
+   fList->Add(fHistPtMCESD);
+   fHistInvPtMCTPC  = new TH2F("fHistInvPtMCTPC","inv pt TPC vs MC",900, 0.0, 9.0,900, 0.0, 9.0);
+   fList->Add(fHistInvPtMCTPC);
+   fHistPtMCTPC     = new TH2F("fHistPtMCTPC"," pt TPC vs MC",300, 0.0, 15.0,300, 0.0, 15.0);
+   fList->Add(fHistPtMCTPC);
+   fHistMomresMCESD = new TH2F("fHistMomresMCESD"," (pt ESD - pt MC)/ptMC vs pt MC",300, 0.0, 15.0,400, -2.0, 2.0);
+   fList->Add(fHistMomresMCESD);
+   fHistMomresMCTPC = new TH2F("fHistMomresMCTPC"," (pt TPC - pt MC)/ptMC vs pt MC",300, 0.0, 15.0,400, -2.0, 2.0);
+   fList->Add(fHistMomresMCTPC);
+
+
+   //user pt shift check
+   fHistUserPtShift = new TH1F("fHistUserPtShift","user defined shift in 1/pt",100,-0.5,1.5);
+   fList->Add(fHistUserPtShift);
+   
+   // esd track cuts  
+   fESDTrackCuts = new AliESDtrackCuts("AliESDtrackCuts");
+  
+ //   //fESDTrackCuts->DefineHistoqgrams(1);
+       fESDTrackCuts->SetRequireSigmaToVertex(kFALSE);
+   fESDTrackCuts->SetRequireTPCRefit(kTRUE);
+   fESDTrackCuts->SetAcceptKinkDaughters(kTRUE);
+   fESDTrackCuts->SetMinNClustersTPC((Int_t)fMinNClustersTPC);
+   fESDTrackCuts->SetMaxChi2PerClusterTPC(fMaxChi2PerClusterTPC);
+   fESDTrackCuts->SetMaxDCAToVertexXY(fMaxDCAtoVertexXY);
+   fESDTrackCuts->SetMaxDCAToVertexZ(fMaxDCAtoVertexZ);
+   fESDTrackCuts->SetDCAToVertex2D(kTRUE);
+   fESDTrackCuts->SetPtRange(fMinPt,fMaxPt); 
+  
+  
+}
+
+//________________________________________________________________________
+void AliPerformancePtCalibMC::SetPtShift(Double_t shiftVal ) { 
+   if(!(shiftVal==0)) { fShift=kTRUE; fDeltaInvP = shiftVal; } 
+}
+
+//________________________________________________________________________
+void AliPerformancePtCalibMC::Exec(AliMCEvent* const mcEvent, AliESDEvent* const esdEvent , AliESDfriend*, Bool_t, Bool_t) 
+{
+  AliStack* stack;
+  if (!esdEvent) {
+    Printf("ERROR: Event not available");
+    return;
+  }
+
+  if (!(esdEvent->GetNumberOfTracks())) {
+    Printf(" PtCalibMC task: There is no track in this event");
+    return;
+  }
+  fHistTrackMultiplicity->Fill(esdEvent->GetNumberOfTracks());
+
+  if (!mcEvent) {
+    Printf("ERROR: Could not retrieve MC event");
+    return;
+  }    
+  stack = mcEvent->Stack();
+  if (!stack) {
+    Printf("ERROR: Could not retrieve stack");
+    return;
+  }
+
+  if(fShift) fHistUserPtShift->Fill(fDeltaInvP);
+  
+  // read primary vertex info
+  Double_t tPrimaryVtxPosition[3];
+  // Double_t tPrimaryVtxCov[3];
+  const AliESDVertex *primaryVtx = esdEvent->GetPrimaryVertexTPC();
+  tPrimaryVtxPosition[0] = primaryVtx->GetXv();
+  tPrimaryVtxPosition[1] = primaryVtx->GetYv();
+  tPrimaryVtxPosition[2] = primaryVtx->GetZv();
+  
+  fHistPrimaryVertexPosX->Fill(tPrimaryVtxPosition[0]);
+  fHistPrimaryVertexPosY->Fill(tPrimaryVtxPosition[1]);
+  fHistPrimaryVertexPosZ->Fill(tPrimaryVtxPosition[2]);
+
+  //fill histos for pt spectra and shift of transverse momentum
+  Int_t count=0;
+  for(Int_t j = 0;j<esdEvent->GetNumberOfTracks();j++){
+    AliESDtrack *ESDTrack = esdEvent->GetTrack(j);
+    if(!ESDTrack) continue;
+
+      //esd track cuts
+    if(fESDcuts == kTRUE){
+             if(!fESDTrackCuts->AcceptTrack(ESDTrack)) continue;
+    }
+    
+    //more track cuts
+    if(fRefitTPC) if(AddTPCcuts(ESDTrack)) continue;
+    if(fRefitITS) if(AddITScuts(ESDTrack)) continue;
+    if(fDCAcut)   if(AddDCAcuts(ESDTrack)) continue ;
+    
+    // get MC info 
+    Int_t esdLabel =ESDTrack->GetLabel();
+    if(esdLabel<0) continue;   
+    TParticle *  partMC = stack->Particle(esdLabel);
+    if (!partMC) continue;
+  
+    // fill correlation histos MC ESD
+    Double_t pESD  = ESDTrack->GetP();
+    Double_t ptESD = ESDTrack->GetSignedPt();
+    
+    if(ptESD==0 || partMC->Pt()==0 ) continue;
+    Double_t mcPt = partMC->Pt();
+    Double_t invPtMC = 1.0/mcPt;
+    Int_t signMC = partMC->GetPdgCode();
+    //MC only
+    if(signMC>0) signMC = 1; 
+    else signMC = -1;
+    //fill MC histos
+    fHistInvPtThetaMC->Fill(signMC*fabs(invPtMC),partMC->Theta());
+    fHistInvPtPhiMC->Fill(signMC*fabs(invPtMC),partMC->Phi());
+    fHistPtThetaMC->Fill(fabs(mcPt),partMC->Theta(),fabs(invPtMC));
+    fHistPtPhiMC->Fill(fabs(mcPt),partMC->Phi(),fabs(invPtMC));
+    //correlation histos MC ESD
+    fHistInvPtMCESD->Fill(fabs(invPtMC),fabs(1.0/ptESD));
+    fHistPtMCESD->Fill(fabs(mcPt),fabs(ptESD));
+
+
+    // fill histos
+    if(fOptTPC == kTRUE){ //TPC tracks
+      const AliExternalTrackParam *TPCTrack = ESDTrack->GetTPCInnerParam(); 
+      if(!TPCTrack) continue;
+      if(fabs(TPCTrack->Eta())> 0.8) continue;
+      
+      Double_t signedPt = TPCTrack->GetSignedPt();
+      Double_t invPt = 0.0;
+      if(!signedPt==0) {
+       invPt = 1.0/signedPt;
+       
+       fHistPtShift0->Fill(fabs(signedPt));
+
+       if(fShift == kTRUE ){
+         invPt += fDeltaInvP; //shift momentum for tests
+         if(!invPt==0) signedPt = 1.0/invPt;
+         else continue;
+       }
+
+       fHistInvPtTheta->Fill(invPt,TPCTrack->Theta());
+       fHistInvPtPhi->Fill(invPt,TPCTrack->Phi());
+       fHistPtTheta->Fill(fabs(signedPt),TPCTrack->Theta());
+       fHistPtPhi->Fill(fabs(signedPt),TPCTrack->Phi());
+
+
+       //correlation histos MC TPC
+       fHistInvPtMCTPC->Fill(fabs(invPtMC),fabs(invPt));
+       fHistPtMCTPC->Fill(fabs(mcPt),fabs(signedPt));
+       
+       //compare to MC info
+       Double_t  ptDiffESD = (fabs(ptESD)-fabs(mcPt))/pow(mcPt,2);
+       Double_t  ptDiffTPC = (fabs(signedPt)-fabs(mcPt))/pow(mcPt,2);
+       Double_t  invPtDiffESD = fabs(1.0/ptESD)-1.0/fabs(mcPt);
+       Double_t  invPtDiffTPC = fabs(invPt)-1.0/fabs(mcPt);
+       Double_t pTPC  = TPCTrack->GetP();
+       
+       if(ESDTrack->GetSign()>0){//compare momenta ESD track and TPC track
+         fHistTPCMomentaPosP->Fill(fabs(pESD),fabs(pTPC));
+         fHistTPCMomentaPosPt->Fill(fabs(ptESD),fabs(signedPt));
+         fHistTPCMomentaPosInvPtMC->Fill(invPtDiffESD,invPtDiffTPC);
+         fHistTPCMomentaPosPtMC->Fill(ptDiffESD,ptDiffTPC);
+       }
+       else{
+         fHistTPCMomentaNegP->Fill(fabs(pESD),fabs(pTPC));
+         fHistTPCMomentaNegPt->Fill(fabs(ptESD),fabs(signedPt));
+         fHistTPCMomentaNegInvPtMC->Fill(invPtDiffESD,invPtDiffTPC);
+         fHistTPCMomentaNegPtMC->Fill(ptDiffESD,ptDiffTPC);
+       }
+       fHistMomresMCESD->Fill((fabs(mcPt)-fabs(ptESD))/fabs(mcPt),fabs(mcPt));
+       fHistMomresMCTPC->Fill((fabs(mcPt)-fabs(signedPt))/fabs(mcPt),fabs(mcPt));
+       count++;
+      }
+      else continue;
+    }
+   
+    else{// ESD tracks
+      Double_t invPt = 0.0;
+      
+      if(!ptESD==0) {
+       invPt = 1.0/ptESD; 
+       fHistPtShift0->Fill(fabs(ptESD));
+       
+       if(fShift == kTRUE ){
+         invPt += fDeltaInvP; //shift momentum for tests
+         if(!invPt==0) ptESD = 1.0/invPt; 
+         else continue;
+       }
+       fHistInvPtTheta->Fill(invPt,ESDTrack->Theta());
+       fHistInvPtPhi->Fill(invPt,ESDTrack->Phi());
+       fHistPtTheta->Fill(ptESD,ESDTrack->Theta());
+       fHistPtPhi->Fill(ptESD,ESDTrack->Phi());
+
+       //differences MC ESD tracks
+       Double_t ptDiffESD = (fabs(ptESD)-fabs(mcPt))/pow(mcPt,2);
+       Double_t invPtdiffESD = fabs(1.0/ptESD)-1.0/fabs(mcPt);
+       if(ESDTrack->GetSign()>0){   
+         fHistESDMomentaPosInvPtMC->Fill(invPtdiffESD);
+         fHistESDMomentaPosPtMC->Fill(ptDiffESD);
+       }
+       else{
+         fHistESDMomentaNegInvPtMC->Fill(invPtdiffESD);
+         fHistESDMomentaNegPtMC->Fill(ptDiffESD);
+       }       
+       
+       fHistMomresMCESD->Fill((fabs(mcPt)-fabs(ptESD))/fabs(mcPt),fabs(mcPt));
+       count++;
+      }
+    }
+  }
+    
+  fHistTrackMultiplicityCuts->Fill(count);
+  
+}    
+
+
+
+//______________________________________________________________________________________________________________________
+Bool_t AliPerformancePtCalibMC::AddTPCcuts(AliESDtrack *ESDTrack){
+  
+  Bool_t cut = kFALSE;
+  
+  if ((ESDTrack->GetStatus()&AliESDtrack::kTPCrefit)==0) cut=kTRUE; // TPC refit
+  if (ESDTrack->GetTPCNcls()<50) cut=kTRUE; // min. nb. TPC clusters
+  if(cut) return kTRUE;
+  return kFALSE;
+}
+//______________________________________________________________________________________________________________________
+Bool_t AliPerformancePtCalibMC::AddDCAcuts(AliESDtrack *ESDTrack){
+  
+  Bool_t cut = kFALSE;
+  
+  Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z and impact parameters:
+  ESDTrack->GetImpactParameters(dca,cov);
+  if(TMath::Abs(dca[0])>3. || TMath::Abs(dca[1])>3.) cut=kTRUE;
+  if(ESDTrack->GetKinkIndex(0)>0) cut=kTRUE;
+  if(cut) return kTRUE;
+  return kFALSE;
+}
+
+//______________________________________________________________________________________________________________________
+Bool_t AliPerformancePtCalibMC::AddITScuts(AliESDtrack *ESDTrack){
+
+  Bool_t cut = kFALSE;
+  
+  if ((ESDTrack->GetStatus()&AliESDtrack::kITSrefit)==0) cut=kTRUE; // ITS refit
+  Int_t clusterITS[200]; 
+  if(ESDTrack->GetITSclusters(clusterITS)<2) cut=kTRUE;  // min. nb. ITS clusters //3
+  
+  if(cut) return kTRUE;
+  return kFALSE;
+}
+
+//______________________________________________________________________________________________________________________
+
+void AliPerformancePtCalibMC::Analyse()
+{
+  
+  
+  AliPerfAnalyzeInvPt *ana = new  AliPerfAnalyzeInvPt("AliPerfAnalyzeInvPt","AliPerfAnalyzeInvPt");
+  
+  TH1::AddDirectory(kFALSE);
+  ana->SetProjBinsTheta(fThetaBins,fNThetaBins);
+  ana->SetProjBinsPhi(fPhiBins,fNPhiBins);
+  ana->SetMakeFitOption(fFitGaus,fExclRange,fRange);
+  
+  TObjArray *aFolderObj = new TObjArray;
+  if(fAnaMC){
+     Printf("AliPerformancePtCalibMC: analysing MC!");
+     ana->StartAnalysis(fHistInvPtThetaMC,fHistInvPtPhiMC, aFolderObj);
+  }
+  else {
+     Printf("AliPerformancePtCalibMC: analysing data!");
+     ana->StartAnalysis(fHistInvPtTheta,fHistInvPtPhi, aFolderObj);
+  }
+  
+  // export objects to analysis folder
+  fAnalysisFolder = ExportToFolder(aFolderObj);
+
+  // delete only TObjArray
+  if(aFolderObj) delete aFolderObj;
+  if(ana) delete ana;
+  
+}
+
+//______________________________________________________________________________________________________________________
+TFolder* AliPerformancePtCalibMC::ExportToFolder(TObjArray * array) 
+{
+  // recreate folder avery time and export objects to new one
+  //
+  AliPerformancePtCalibMC * comp=this;
+  TFolder *folder = comp->GetAnalysisFolder();
+
+  TString name, title;
+  TFolder *newFolder = 0;
+  Int_t i = 0;
+  Int_t size = array->GetSize();
+
+  if(folder) { 
+    // get name and title from old folder
+    name = folder->GetName();  
+    title = folder->GetTitle();  
+
+    // delete old one
+    delete folder;
+
+    // create new one
+    newFolder = CreateFolder(name.Data(),title.Data());
+    newFolder->SetOwner();
+
+    // add objects to folder
+    while(i < size) {
+      newFolder->Add(array->At(i));
+      i++;
+    }
+  }
+
+  return newFolder;
+}
+
+//______________________________________________________________________________________________________________________
+Long64_t AliPerformancePtCalibMC::Merge(TCollection* const list) 
+{
+  // Merge list of objects (needed by PROOF)
+
+  if (!list)
+    return 0;
+
+  if (list->IsEmpty())
+    return 1;
+
+  TIterator* iter = list->MakeIterator();
+  TObject* obj = 0;
+
+  // collection of generated histograms
+  Int_t count=0;
+  while((obj = iter->Next()) != 0) 
+    {
+      AliPerformancePtCalibMC* entry = dynamic_cast<AliPerformancePtCalibMC*>(obj);
+      if (entry == 0) continue; 
+  
+      fHistInvPtTheta->Add(entry->fHistInvPtTheta);
+      fHistInvPtPhi->Add(entry-> fHistInvPtPhi);
+      fHistPtTheta->Add(entry->fHistPtTheta);
+      fHistPtPhi->Add(entry->fHistPtPhi);
+
+      fHistInvPtThetaMC->Add(entry->fHistInvPtThetaMC);
+      fHistInvPtPhiMC->Add(entry->fHistInvPtPhiMC);
+      fHistPtThetaMC->Add(entry->fHistPtThetaMC);
+      fHistPtPhiMC->Add(entry->fHistPtPhiMC);
+      fHistInvPtMCESD->Add(entry->fHistInvPtMCESD);
+      fHistPtMCESD->Add(entry->fHistPtMCESD);
+      fHistInvPtMCTPC->Add(entry->fHistInvPtMCTPC);
+      fHistPtMCTPC->Add(entry->fHistPtMCTPC);
+      fHistMomresMCESD->Add(entry->fHistMomresMCESD);
+      fHistMomresMCTPC->Add(entry->fHistMomresMCTPC);
+
+      fHistPtShift0->Add(entry->fHistPtShift0);
+      fHistPrimaryVertexPosX->Add(entry->fHistPrimaryVertexPosX);
+      fHistPrimaryVertexPosY->Add(entry->fHistPrimaryVertexPosY);
+      fHistPrimaryVertexPosZ->Add(entry->fHistPrimaryVertexPosZ);
+      fHistTrackMultiplicity->Add(entry->fHistTrackMultiplicity);
+      fHistTrackMultiplicityCuts->Add(entry->fHistTrackMultiplicityCuts);
+      
+      fHistTPCMomentaPosP->Add(entry->fHistTPCMomentaPosP);
+      fHistTPCMomentaNegP->Add(entry->fHistTPCMomentaNegP);
+      fHistTPCMomentaPosPt->Add(entry->fHistTPCMomentaPosPt);
+      fHistTPCMomentaNegPt->Add(entry->fHistTPCMomentaNegPt);
+      fHistTPCMomentaPosInvPtMC->Add(entry->fHistTPCMomentaPosInvPtMC);
+      fHistTPCMomentaNegInvPtMC->Add(entry->fHistTPCMomentaNegInvPtMC);
+      fHistTPCMomentaPosPtMC->Add(entry->fHistTPCMomentaPosPtMC);
+      fHistTPCMomentaNegPtMC->Add(entry->fHistTPCMomentaNegPtMC);
+      fHistESDMomentaPosInvPtMC->Add(entry->fHistESDMomentaPosInvPtMC);
+      fHistESDMomentaNegInvPtMC->Add(entry->fHistESDMomentaNegInvPtMC);
+      fHistESDMomentaPosPtMC->Add(entry->fHistESDMomentaPosPtMC);
+      fHistESDMomentaNegPtMC->Add(entry->fHistESDMomentaNegPtMC);
+      count++;
+    }
+  
+  return count;
+}
+
+//______________________________________________________________________________________________________________________
+TFolder* AliPerformancePtCalibMC::CreateFolder(TString name,TString title) { 
+  // create folder for analysed histograms
+  //
+  TFolder *folder = 0;
+  folder = new TFolder(name.Data(),title.Data());
+
+  return folder;
+}
+
+
+// set variables for Analyse()
+
+//______________________________________________________________________________________________________________________
+void AliPerformancePtCalibMC::SetProjBinsPhi(const Double_t *phiBinArray,Int_t nphBins){
+   
+   fNPhiBins = nphBins;
+  
+   for(Int_t k = 0;k<fNPhiBins;k++){
+      fPhiBins[k] = phiBinArray[k];
+   }
+   Printf("AliPerformancePtCalibMC: number of bins in phi set to %i",fNPhiBins);
+
+}
+//____________________________________________________________________________________________________________________________________________
+void AliPerformancePtCalibMC::SetProjBinsTheta(const Double_t *thetaBinArray, Int_t nthBins){
+  
+   fNThetaBins = nthBins;
+   for(Int_t k = 0;k<fNThetaBins;k++){
+      fThetaBins[k] = thetaBinArray[k];
+   }
+   Printf("AliPerformancePtCalibMC: number of bins in theta set to %i",fNThetaBins);
+
+}
+//____________________________________________________________________________________________________________________________________________
+void AliPerformancePtCalibMC::SetMakeFitOption(const Bool_t setGausFit, const Double_t exclusionR,const Double_t fitR ){
+
+  
+   
+   fFitGaus = setGausFit;
+   fExclRange  = exclusionR;
+   fRange = fitR;
+  
+   if(fFitGaus) Printf("AliPerformancePtCalibMC:set MakeGausFit with fit range %2.3f and exclusion range in 1/pt: %2.3f",fRange,fExclRange);
+   else  Printf("AliPerformancePtCalibMC: set standard polynomial fit with fit range %2.3f and exclusion range in 1/pt: %2.3f",fRange,fExclRange);
+}
diff --git a/PWG1/TPC/AliPerformancePtCalibMC.h b/PWG1/TPC/AliPerformancePtCalibMC.h
new file mode 100755 (executable)
index 0000000..b627df4
--- /dev/null
@@ -0,0 +1,175 @@
+#ifndef ALIPERFORMANCEPTCALIBMC_H
+#define ALIPERFORMANCEPTCALIBMC_H
+
+class TString;
+class TNamed;
+class TCanvas;
+class TH1F;
+class TH2F;
+class TList;
+
+class AliESDVertex;
+class AliESDtrack;
+class AliMCEvent;
+class AliStack;
+class AliTrackReference;
+class AliESDEvent; 
+class AliESDfriend; 
+class AliESDfriendTrack; 
+class AliMCEvent;
+class AliMCParticle;
+class AliMCInfoCuts;
+class AliRecInfoCuts;
+class AliESDtrackCuts;
+
+#include "AliPerformanceObject.h"
+
+class AliPerformancePtCalibMC : public AliPerformanceObject {
+ public:
+  AliPerformancePtCalibMC();
+  AliPerformancePtCalibMC(const char *name, const char *title);
+   virtual ~AliPerformancePtCalibMC() ;
+
+  // Init data members
+  virtual void  Init();
+
+  // Execute analysis
+  virtual void  Exec(AliMCEvent* const mcEvent, AliESDEvent *const esdEvent, AliESDfriend *const esdFriend, const Bool_t bUseMC, const Bool_t bUseESDfriend);
+
+  // Merge output objects (needed by PROOF) 
+  virtual Long64_t Merge(TCollection* const list);
+
+  // Analyse output histograms
+  virtual void Analyse();
+
+  // Get analysis folder
+  virtual TFolder* GetAnalysisFolder() const {return fAnalysisFolder;}
+
+   Bool_t AddTPCcuts(AliESDtrack *ESDTrack);
+   Bool_t AddITScuts(AliESDtrack *ESDTrack);
+   Bool_t AddDCAcuts(AliESDtrack *ESDTrack);
+
+   void SetReadTPCTracksOff(){fOptTPC   = kFALSE;}// function  can be used in the macro to set the variables
+   void SetTPCRefitOn()      {fRefitTPC = kTRUE;} // function  can be used in the macro to set the variables
+   void SetITSRefitOn()      {fRefitITS = kTRUE;} // function  can be used in the macro to set the variables
+   void SetESDcutsOff()      {fESDcuts  = kFALSE;}// function  can be used in the macro to set the variables
+   void SetDCAcutsOff()      {fDCAcut   = kFALSE;}// function  can be used in the macro to set the variables
+
+   //void SetPtShift(Double_t shiftVal ){if(!(shiftVal==0)) {fShift=kTRUE; fDeltaInvP = shiftVal;} };
+   void SetPtShift(Double_t shiftVal); //shift in 1/pt
+
+   void SetESDcutValues(Double_t * esdCutValues){// functions can be used in the macro to set the variables
+      fMinPt                = esdCutValues[0]; 
+      fMaxPt                = esdCutValues[1];
+      fMinNClustersTPC      = esdCutValues[2];
+      fMaxChi2PerClusterTPC = esdCutValues[3];
+      fMaxDCAtoVertexXY     = esdCutValues[4];
+      fMaxDCAtoVertexZ      = esdCutValues[5];
+   }
+   
+   void SetProjBinsPhi(const Double_t *pBins,Int_t sizep);
+   void SetProjBinsTheta(const Double_t *tBins, Int_t sizet);
+   void SetMakeFitOption(const Bool_t setGausFit, const Double_t exclusionR,const Double_t fitR );
+   void SetAnaMCOff() {fAnaMC = kFALSE;}
+   TList *GetHistoList() {return fList;}
+   
+  // Create folder for analysed histograms
+  TFolder *CreateFolder(TString folder = "folderPtCalib",TString title = "Analysed PtCalib histograms");
+
+  // Export objects to folder
+  TFolder *ExportToFolder(TObjArray * array=0);
+
+  // Selection cuts
+  void SetAliRecInfoCuts(AliRecInfoCuts* const cuts=0) {fCutsRC = cuts;}   
+  void SetAliMCInfoCuts(AliMCInfoCuts* const cuts=0) {fCutsMC = cuts;}
+  
+   
+  AliRecInfoCuts*  GetAliRecInfoCuts() const {return fCutsRC;}  
+  AliMCInfoCuts*   GetAliMCInfoCuts()  const {return fCutsMC;}
+
+protected:
+   
+  Double_t fThetaBins[100];
+  Double_t fPhiBins[100];
+   
+   Int_t fNThetaBins;
+   Int_t fNPhiBins ;
+   Double_t fRange;
+   Double_t fExclRange ;
+   Bool_t fFitGaus ;
+   Bool_t  fAnaMC;
+   
+    
+ private:
+  Bool_t fShift;
+  Double_t fDeltaInvP;
+  //options for cuts
+  Bool_t fOptTPC;
+  Bool_t fESDcuts;
+  Bool_t fRefitTPC;
+  Bool_t fRefitITS;
+  Bool_t fDCAcut;
+  
+  Double_t fMinPt;
+  Double_t fMaxPt;
+  Double_t fMinNClustersTPC;
+  Double_t fMaxChi2PerClusterTPC;
+  Double_t fMaxDCAtoVertexXY;
+  Double_t fMaxDCAtoVertexZ;
+
+  AliRecInfoCuts* fCutsRC;     // selection cuts for reconstructed tracks
+  AliMCInfoCuts*  fCutsMC;     // selection cuts for MC tracks
+
+   
+  TList       *fList;
+  TH2F        *fHistInvPtTheta;
+  TH2F        *fHistInvPtPhi;
+  TH2F        *fHistPtTheta;
+  TH2F        *fHistPtPhi;
+
+  TH1F        *fHistPtShift0;
+  TH1F        *fHistPrimaryVertexPosX;          
+  TH1F        *fHistPrimaryVertexPosY;          
+  TH1F        *fHistPrimaryVertexPosZ;          
+  TH1F        *fHistTrackMultiplicity;          
+  TH1F        *fHistTrackMultiplicityCuts;
+
+  TH2F        *fHistTPCMomentaPosP;
+  TH2F        *fHistTPCMomentaNegP;
+  TH2F        *fHistTPCMomentaPosPt;
+  TH2F        *fHistTPCMomentaNegPt;
+   
+  TH2F        *fHistInvPtThetaMC;
+  TH2F        *fHistInvPtPhiMC;
+  TH2F        *fHistPtThetaMC;
+  TH2F        *fHistPtPhiMC;   
+  TH2F        *fHistInvPtMCESD;
+  TH2F        *fHistInvPtMCTPC;
+  TH2F        *fHistPtMCESD;
+  TH2F        *fHistPtMCTPC;
+  TH2F        *fHistMomresMCESD;   
+  TH2F        *fHistMomresMCTPC;   
+  TH2F        *fHistTPCMomentaPosInvPtMC;
+  TH2F        *fHistTPCMomentaNegInvPtMC;
+  TH2F        *fHistTPCMomentaPosPtMC;
+  TH2F        *fHistTPCMomentaNegPtMC;
+  
+  TH1F        *fHistESDMomentaPosInvPtMC;
+  TH1F        *fHistESDMomentaNegInvPtMC;
+  TH1F        *fHistESDMomentaPosPtMC;
+   TH1F        *fHistESDMomentaNegPtMC;
+   TH1F        *fHistUserPtShift;
+   
+  AliESDtrackCuts* fESDTrackCuts;
+  
+  // analysis folder 
+  TFolder *fAnalysisFolder; // folder for analysed histograms
+
+  AliPerformancePtCalibMC(const AliPerformancePtCalibMC&);            // not implemented 
+  AliPerformancePtCalibMC& operator=(const AliPerformancePtCalibMC&); // not implemented 
+
+  ClassDef(AliPerformancePtCalibMC, 1); 
+};
+
+#endif