#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);
 
}