[u/mrichter/AliRoot.git] / PWG2 / SPECTRA / AliProtonAnalysis.cxx

/**************************************************************************
 * Author: Panos Christakoglou.                                           *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* $Id$ */

//-----------------------------------------------------------------
//                 AliProtonAnalysis class
//   This is the class to deal with the proton analysis
//   Origin: Panos Christakoglou, UOA-CERN, Panos.Christakoglou@cern.ch
//-----------------------------------------------------------------
#include <Riostream.h>
#include <TFile.h>
#include <TSystem.h>
#include <TF1.h>
#include <TH2F.h>
#include <TH1D.h>

#include "AliProtonAnalysis.h"

#include <AliAODEvent.h>
#include <AliESDEvent.h>
#include <AliLog.h>
#include <AliPID.h>

ClassImp(AliProtonAnalysis)

//____________________________________________________________________//
AliProtonAnalysis::AliProtonAnalysis() : 
  TObject(), 
  fNBinsY(0), fMinY(0), fMaxY(0),
  fNBinsPt(0), fMinPt(0), fMaxPt(0),
  fMinTPCClusters(0), fMinITSClusters(0),
  fMaxChi2PerTPCCluster(0), fMaxChi2PerITSCluster(0),
  fMaxCov11(0), fMaxCov22(0), fMaxCov33(0), fMaxCov44(0), fMaxCov55(0),
  fMaxSigmaToVertex(0),
  fMinTPCClustersFlag(kFALSE), fMinITSClustersFlag(kFALSE),
  fMaxChi2PerTPCClusterFlag(kFALSE), fMaxChi2PerITSClusterFlag(kFALSE),
  fMaxCov11Flag(kFALSE), fMaxCov22Flag(kFALSE), fMaxCov33Flag(kFALSE), fMaxCov44Flag(kFALSE), fMaxCov55Flag(kFALSE),
  fMaxSigmaToVertexFlag(kFALSE),
  fITSRefitFlag(kFALSE), fTPCRefitFlag(kFALSE),
  fFunctionProbabilityFlag(kFALSE), 
  fElectronFunction(0), fMuonFunction(0),
  fPionFunction(0), fKaonFunction(0), fProtonFunction(0),
  fHistYPtProtons(0), fHistYPtAntiProtons(0) {
  //Default constructor
  for(Int_t i = 0; i < 5; i++) fPartFrac[i] = 0.0;
}

//____________________________________________________________________//
AliProtonAnalysis::AliProtonAnalysis(Int_t nbinsY, Float_t fLowY, Float_t fHighY,Int_t nbinsPt, Float_t fLowPt, Float_t fHighPt) : 
  TObject(),
  fNBinsY(nbinsY), fMinY(fLowY), fMaxY(fHighY),
  fNBinsPt(nbinsPt), fMinPt(fLowPt), fMaxPt(fHighPt),
  fMinTPCClusters(0), fMinITSClusters(0),
  fMaxChi2PerTPCCluster(0), fMaxChi2PerITSCluster(0),
  fMaxCov11(0), fMaxCov22(0), fMaxCov33(0), fMaxCov44(0), fMaxCov55(0),
  fMaxSigmaToVertex(0),
  fMinTPCClustersFlag(kFALSE), fMinITSClustersFlag(kFALSE),
  fMaxChi2PerTPCClusterFlag(kFALSE), fMaxChi2PerITSClusterFlag(kFALSE),
  fMaxCov11Flag(kFALSE), fMaxCov22Flag(kFALSE), fMaxCov33Flag(kFALSE), fMaxCov44Flag(kFALSE), fMaxCov55Flag(kFALSE),
  fMaxSigmaToVertexFlag(kFALSE),
  fITSRefitFlag(kFALSE), fTPCRefitFlag(kFALSE),
  fFunctionProbabilityFlag(kFALSE), 
  fElectronFunction(0), fMuonFunction(0),
  fPionFunction(0), fKaonFunction(0), fProtonFunction(0),
  fHistYPtProtons(0), fHistYPtAntiProtons(0) {
  //Default constructor

  fHistYPtProtons = new TH2F("fHistYPtProtons","y-Pt Protons",fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
  fHistYPtProtons->SetStats(kTRUE);
  fHistYPtProtons->GetYaxis()->SetTitle("P_{T} [GeV]");
  fHistYPtProtons->GetXaxis()->SetTitle("y");
  fHistYPtProtons->GetXaxis()->SetTitleColor(1);

  fHistYPtAntiProtons = new TH2F("fHistYPtAntiProtons","y-Pt Antiprotons",fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
  fHistYPtAntiProtons->SetStats(kTRUE);
  fHistYPtAntiProtons->GetYaxis()->SetTitle("P_{T} [GeV]");
  fHistYPtAntiProtons->GetXaxis()->SetTitle("y");
  fHistYPtAntiProtons->GetXaxis()->SetTitleColor(1);
} 

//____________________________________________________________________//
AliProtonAnalysis::~AliProtonAnalysis() {
  //Default destructor
  
}

//____________________________________________________________________//
void AliProtonAnalysis::InitHistograms(Int_t nbinsY, Float_t fLowY, Float_t fHighY, Int_t nbinsPt, Float_t fLowPt, Float_t fHighPt) {
  fNBinsY = nbinsY;
  fMinY = fLowY;
  fMaxY = fHighY;
  fNBinsPt = nbinsPt;
  fMinPt = fLowPt;
  fMaxPt = fHighPt;

  fHistYPtProtons = new TH2F("fHistYPtProtons","y-Pt Protons",fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
  fHistYPtProtons->SetStats(kTRUE);
  fHistYPtProtons->GetYaxis()->SetTitle("P_{T} [GeV]");
  fHistYPtProtons->GetXaxis()->SetTitle("y");
  fHistYPtProtons->GetXaxis()->SetTitleColor(1);

  fHistYPtAntiProtons = new TH2F("fHistYPtAntiProtons","y-Pt Antiprotons",fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
  fHistYPtAntiProtons->SetStats(kTRUE);
  fHistYPtAntiProtons->GetYaxis()->SetTitle("P_{T} [GeV]");
  fHistYPtAntiProtons->GetXaxis()->SetTitle("y");
  fHistYPtAntiProtons->GetXaxis()->SetTitleColor(1);
}

//____________________________________________________________________//
void AliProtonAnalysis::ReadFromFile(const char* filename) {
  TFile *file = TFile::Open(filename);
  fHistYPtProtons = (TH2F *)file->Get("fHistYPtProtons");
  fHistYPtAntiProtons = (TH2F *)file->Get("fHistYPtAntiProtons");
  fHistYPtProtons->Sumw2();
  fHistYPtAntiProtons->Sumw2();
}

//____________________________________________________________________//
TH1D *AliProtonAnalysis::GetProtonYHistogram() {
  TH1D *fYProtons = (TH1D *)fHistYPtProtons->ProjectionX("fYProtons",0,fHistYPtProtons->GetYaxis()->GetNbins(),"e"); 
  fYProtons->SetStats(kFALSE);
  fYProtons->GetYaxis()->SetTitle("dN/dy");
  fYProtons->SetTitle("dN/dy protons");
  fYProtons->SetMarkerStyle(kFullCircle);
  fYProtons->SetMarkerColor(4);

  return fYProtons;
}

//____________________________________________________________________//
TH1D *AliProtonAnalysis::GetAntiProtonYHistogram() {
  TH1D *fYAntiProtons = (TH1D *)fHistYPtAntiProtons->ProjectionX("fYAntiProtons",0,fHistYPtAntiProtons->GetYaxis()->GetNbins(),"e"); 
  fYAntiProtons->SetStats(kFALSE);
  fYAntiProtons->GetYaxis()->SetTitle("dN/dy");
  fYAntiProtons->SetTitle("dN/dy antiprotons");
  fYAntiProtons->SetMarkerStyle(kFullCircle);
  fYAntiProtons->SetMarkerColor(4);

  return fYAntiProtons;
}

//____________________________________________________________________//
TH1D *AliProtonAnalysis::GetProtonPtHistogram() {
  TH1D *fPtProtons = (TH1D *)fHistYPtProtons->ProjectionY("fPtProtons",0,fHistYPtProtons->GetXaxis()->GetNbins(),"e"); 
  fPtProtons->SetStats(kFALSE);
  fPtProtons->GetYaxis()->SetTitle("dN/dP_{T}");
  fPtProtons->SetTitle("dN/dPt protons");
  fPtProtons->SetMarkerStyle(kFullCircle);
  fPtProtons->SetMarkerColor(4);

  return fPtProtons;
}

//____________________________________________________________________//
TH1D *AliProtonAnalysis::GetAntiProtonPtHistogram() {
  TH1D *fPtAntiProtons = (TH1D *)fHistYPtAntiProtons->ProjectionY("fPtAntiProtons",0,fHistYPtProtons->GetXaxis()->GetNbins(),"e"); 
  fPtAntiProtons->SetStats(kFALSE);
  fPtAntiProtons->GetYaxis()->SetTitle("dN/dP_{T}");
  fPtAntiProtons->SetTitle("dN/dPt antiprotons");
  fPtAntiProtons->SetMarkerStyle(kFullCircle);
  fPtAntiProtons->SetMarkerColor(4);

  return fPtAntiProtons;
}

//____________________________________________________________________//
TH1D *AliProtonAnalysis::GetYRatioHistogram() {
  TH1D *fYProtons = GetProtonYHistogram();
  TH1D *fYAntiProtons = GetAntiProtonYHistogram();
  
  TH1D *hRatioY = new TH1D("hRatioY","",fYProtons->GetNbinsX(),fYProtons->GetXaxis()->GetXmin(),fYProtons->GetXaxis()->GetXmax());
  hRatioY->Divide(fYAntiProtons,fYProtons,1.0,1.0);
  hRatioY->SetMarkerStyle(kFullCircle);
  hRatioY->SetMarkerColor(4);
  hRatioY->GetYaxis()->SetTitle("#bar{p}/p");
  hRatioY->GetYaxis()->SetTitleOffset(1.4);
  hRatioY->GetXaxis()->SetTitle("y");
  hRatioY->GetXaxis()->SetTitleColor(1);
  hRatioY->SetStats(kFALSE);

  return hRatioY;
}

//____________________________________________________________________//
TH1D *AliProtonAnalysis::GetPtRatioHistogram() {
  TH1D *fPtProtons = GetProtonPtHistogram();
  TH1D *fPtAntiProtons = GetAntiProtonPtHistogram();
  
  TH1D *hRatioPt = new TH1D("hRatioPt","",fPtProtons->GetNbinsX(),fPtProtons->GetXaxis()->GetXmin(),fPtProtons->GetXaxis()->GetXmax());
  hRatioPt->Divide(fPtAntiProtons,fPtProtons,1.0,1.0);
  hRatioPt->SetMarkerStyle(kFullCircle);
  hRatioPt->SetMarkerColor(4);
  hRatioPt->GetYaxis()->SetTitle("#bar{p}/p");
  hRatioPt->GetYaxis()->SetTitleOffset(1.4);
  hRatioPt->GetXaxis()->SetTitle("P_{T} [GeV/c]");
  hRatioPt->GetXaxis()->SetTitleColor(1);
  hRatioPt->SetStats(kFALSE);

  return hRatioPt;
}

//____________________________________________________________________//
TH1D *AliProtonAnalysis::GetYAsymmetryHistogram() {
  TH1D *fYProtons = GetProtonYHistogram();
  TH1D *fYAntiProtons = GetAntiProtonYHistogram();
  
  TH1D *hsum = new TH1D("hsumY","",fYProtons->GetNbinsX(),fYProtons->GetXaxis()->GetXmin(),fYProtons->GetXaxis()->GetXmax());
  hsum->Add(fYProtons,fYAntiProtons,1.0,1.0);

  TH1D *hdiff = new TH1D("hdiffY","",fYProtons->GetNbinsX(),fYProtons->GetXaxis()->GetXmin(),fYProtons->GetXaxis()->GetXmax());
  hdiff->Add(fYProtons,fYAntiProtons,1.0,-1.0);

  TH1D *hAsymmetryY = new TH1D("hAsymmetryY","",fYProtons->GetNbinsX(),fYProtons->GetXaxis()->GetXmin(),fYProtons->GetXaxis()->GetXmax());
  hAsymmetryY->Divide(hdiff,hsum,2.0,1.);
  hAsymmetryY->SetMarkerStyle(kFullCircle);
  hAsymmetryY->SetMarkerColor(4);
  hAsymmetryY->GetYaxis()->SetTitle("A_{p}");
  hAsymmetryY->GetYaxis()->SetTitleOffset(1.4);
  hAsymmetryY->GetXaxis()->SetTitle("y");
  hAsymmetryY->GetXaxis()->SetTitleColor(1);
  hAsymmetryY->SetStats(kFALSE);

  return hAsymmetryY;
}

//____________________________________________________________________//
TH1D *AliProtonAnalysis::GetPtAsymmetryHistogram() {
  TH1D *fPtProtons = GetProtonPtHistogram();
  TH1D *fPtAntiProtons = GetAntiProtonPtHistogram();
  
  TH1D *hsum = new TH1D("hsumPt","",fPtProtons->GetNbinsX(),fPtProtons->GetXaxis()->GetXmin(),fPtProtons->GetXaxis()->GetXmax());
  hsum->Add(fPtProtons,fPtAntiProtons,1.0,1.0);

  TH1D *hdiff = new TH1D("hdiffPt","",fPtProtons->GetNbinsX(),fPtProtons->GetXaxis()->GetXmin(),fPtProtons->GetXaxis()->GetXmax());
  hdiff->Add(fPtProtons,fPtAntiProtons,1.0,-1.0);

  TH1D *hAsymmetryPt = new TH1D("hAsymmetryPt","",fPtProtons->GetNbinsX(),fPtProtons->GetXaxis()->GetXmin(),fPtProtons->GetXaxis()->GetXmax());
  hAsymmetryPt->Divide(hdiff,hsum,2.0,1.);
  hAsymmetryPt->SetMarkerStyle(kFullCircle);
  hAsymmetryPt->SetMarkerColor(4);
  hAsymmetryPt->GetYaxis()->SetTitle("A_{p}");
  hAsymmetryPt->GetYaxis()->SetTitleOffset(1.4);
  hAsymmetryPt->GetXaxis()->SetTitle("P_{T} [GeV/c]");
  hAsymmetryPt->GetXaxis()->SetTitleColor(1);
  hAsymmetryPt->SetStats(kFALSE);

  return hAsymmetryPt;
}

//____________________________________________________________________//
Double_t AliProtonAnalysis::GetParticleFraction(Int_t i, Double_t p) {
  Double_t partFrac=0;
  if(fFunctionProbabilityFlag) {
    if(i == 0) partFrac = fElectronFunction->Eval(p);
    if(i == 1) partFrac = fMuonFunction->Eval(p);
    if(i == 2) partFrac = fPionFunction->Eval(p);
    if(i == 3) partFrac = fKaonFunction->Eval(p);
    if(i == 4) partFrac = fProtonFunction->Eval(p);
  }
  else partFrac = fPartFrac[i];

  return partFrac;
}

//____________________________________________________________________//
void AliProtonAnalysis::Analyze(AliESDEvent* fESD) {
  //Main analysis part
  Int_t nGoodTracks = fESD->GetNumberOfTracks();
  for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) {
    AliESDtrack* track = fESD->GetTrack(iTracks);
    if(IsAccepted(track)) {
      Double_t Pt = track->Pt();
      Double_t P = track->P();
	
      //pid
      Double_t probability[5];
	  track->GetESDpid(probability);
      Double_t rcc = 0.0;
      for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += probability[i]*GetParticleFraction(i,P);
      if(rcc == 0.0) continue;
      Double_t w[5];
      for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = probability[i]*GetParticleFraction(i,P)/rcc;
      Long64_t fParticleType = TMath::LocMax(AliPID::kSPECIES,w);
      if(fParticleType == 4) {
        if(track->Charge() > 0) fHistYPtProtons->Fill(Rapidity(track),Pt);
        else if(track->Charge() < 0) fHistYPtAntiProtons->Fill(Rapidity(track),Pt);
      }//proton check
    }//cuts
  }//track loop 
}

//____________________________________________________________________//
void AliProtonAnalysis::Analyze(AliAODEvent* fAOD) {
  //Main analysis part
  Int_t nGoodTracks = fAOD->GetNumberOfTracks();
  for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) {
    AliAODTrack* track = fAOD->GetTrack(iTracks);
    Double_t Pt = track->Pt();
    Double_t P = track->P();
    
    //pid
    Double_t probability[10];
    track->GetPID(probability);
    Double_t rcc = 0.0;
    for(Int_t i = 0; i < AliPID::kSPECIESN; i++) rcc += probability[i]*GetParticleFraction(i,P);
    if(rcc == 0.0) continue;
    Double_t w[10];
    for(Int_t i = 0; i < AliPID::kSPECIESN; i++) w[i] = probability[i]*GetParticleFraction(i,P)/rcc;
    Long64_t fParticleType = TMath::LocMax(AliPID::kSPECIESN,w);
    if(fParticleType == 4) {
      if(track->Charge() > 0) fHistYPtProtons->Fill(track->Y(fParticleType),Pt);
      else if(track->Charge() < 0) fHistYPtAntiProtons->Fill(track->Y(fParticleType),Pt);
    }//proton check
  }//track loop 
}

//____________________________________________________________________//
Bool_t AliProtonAnalysis::IsAccepted(AliESDtrack* track) {
  // Checks if the track is excluded from the cuts
  Int_t  fIdxInt[200];
  Int_t nClustersITS = track->GetITSclusters(fIdxInt);
  Int_t nClustersTPC = track->GetTPCclusters(fIdxInt);

  Float_t chi2PerClusterITS = -1;
  Float_t chi2PerClusterTPC = -1;
  if (nClustersTPC!=0)
    chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC);

  Double_t extCov[15];
  track->GetExternalCovariance(extCov);

  Double_t Pt = track->Pt();

  if(fMinITSClustersFlag)
    if(nClustersITS < fMinITSClusters) return kFALSE;
  if(fMinTPCClustersFlag)
    if(nClustersTPC < fMinTPCClusters) return kFALSE;
  if(fMaxChi2PerTPCClusterFlag)
    if(chi2PerClusterTPC > fMaxChi2PerTPCCluster) return kFALSE; 
  if(fMaxChi2PerITSClusterFlag)
    if(chi2PerClusterITS > fMaxChi2PerITSCluster) return kFALSE; 
  if(fMaxCov11Flag)
    if(extCov[0] > fMaxCov11) return kFALSE;
  if(fMaxCov22Flag)
    if(extCov[2] > fMaxCov22) return kFALSE;
  if(fMaxCov33Flag)
    if(extCov[5] > fMaxCov33) return kFALSE;
  if(fMaxCov44Flag)
    if(extCov[9] > fMaxCov44) return kFALSE;
  if(fMaxCov55Flag)
    if(extCov[14] > fMaxCov55) return kFALSE;
  if(fMaxSigmaToVertexFlag)
    if(GetSigmaToVertex(track) > fMaxSigmaToVertex) return kFALSE;
  if(fITSRefitFlag)
    if ((track->GetStatus() & AliESDtrack::kITSrefit) == 0) return kFALSE;
  if(fTPCRefitFlag)
    if ((track->GetStatus() & AliESDtrack::kTPCrefit) == 0) return kFALSE;

  if((Pt < fMinPt) || (Pt > fMaxPt)) return kFALSE;

  return kTRUE;
}

//____________________________________________________________________//
Float_t AliProtonAnalysis::GetSigmaToVertex(AliESDtrack* esdTrack) {
  // Calculates the number of sigma to the vertex.
  
  Float_t b[2];
  Float_t bRes[2];
  Float_t bCov[3];
  esdTrack->GetImpactParameters(b,bCov);
  if (bCov[0]<=0 || bCov[2]<=0) {
    //AliDebug(1, "Estimated b resolution lower or equal zero!");
    bCov[0]=0; bCov[2]=0;
  }
  bRes[0] = TMath::Sqrt(bCov[0]);
  bRes[1] = TMath::Sqrt(bCov[2]);
  
  if (bRes[0] == 0 || bRes[1] ==0) return -1;
  
  Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
  
  if (TMath::Exp(-d * d / 2) < 1e-10) return 1000;
  
  d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
  
  return d;
}

Double_t AliProtonAnalysis::Rapidity(AliESDtrack *track) {
  Double_t fMass = 9.38270000000000048e-01;
  
  Double_t P = TMath::Sqrt(TMath::Power(track->Px(),2) + 
                           TMath::Power(track->Py(),2) + 
						   TMath::Power(track->Pz(),2));
  Double_t energy = TMath::Sqrt(P*P + fMass*fMass);
  Double_t y = -999;
  if(energy != track->Pz()) 
    y = 0.5*TMath::Log((energy + track->Pz())/(energy - track->Pz()));

  return y;
}
Commit	Line	Data
734d2c12	1	/**************************************************************************
	2	* Author: Panos Christakoglou. *
	3	* Contributors are mentioned in the code where appropriate. *
	4	* *
	5	* Permission to use, copy, modify and distribute this software and its *
	6	* documentation strictly for non-commercial purposes is hereby granted *
	7	* without fee, provided that the above copyright notice appears in all *
	8	* copies and that both the copyright notice and this permission notice *
	9	* appear in the supporting documentation. The authors make no claims *
	10	* about the suitability of this software for any purpose. It is *
	11	* provided "as is" without express or implied warranty. *
	12	**************************************************************************/
	13
	14	/* $Id$ */
	15
	16	//-----------------------------------------------------------------
	17	// AliProtonAnalysis class
	18	// This is the class to deal with the proton analysis
	19	// Origin: Panos Christakoglou, UOA-CERN, Panos.Christakoglou@cern.ch
	20	//-----------------------------------------------------------------
	21	#include <Riostream.h>
	22	#include <TFile.h>
	23	#include <TSystem.h>
aafecd8b	24	#include <TF1.h>
734d2c12	25	#include <TH2F.h>
	26	#include <TH1D.h>
	27
	28	#include "AliProtonAnalysis.h"
	29
ee4ca40d	30	#include <AliAODEvent.h>
734d2c12	31	#include <AliESDEvent.h>
734d2c12	32	#include <AliLog.h>
ee4ca40d	33	#include <AliPID.h>
734d2c12	34
	35	ClassImp(AliProtonAnalysis)
	36
	37	//____________________________________________________________________//
	38	AliProtonAnalysis::AliProtonAnalysis() :
	39	TObject(),
	40	fNBinsY(0), fMinY(0), fMaxY(0),
	41	fNBinsPt(0), fMinPt(0), fMaxPt(0),
	42	fMinTPCClusters(0), fMinITSClusters(0),
	43	fMaxChi2PerTPCCluster(0), fMaxChi2PerITSCluster(0),
	44	fMaxCov11(0), fMaxCov22(0), fMaxCov33(0), fMaxCov44(0), fMaxCov55(0),
	45	fMaxSigmaToVertex(0),
	46	fMinTPCClustersFlag(kFALSE), fMinITSClustersFlag(kFALSE),
	47	fMaxChi2PerTPCClusterFlag(kFALSE), fMaxChi2PerITSClusterFlag(kFALSE),
	48	fMaxCov11Flag(kFALSE), fMaxCov22Flag(kFALSE), fMaxCov33Flag(kFALSE), fMaxCov44Flag(kFALSE), fMaxCov55Flag(kFALSE),
	49	fMaxSigmaToVertexFlag(kFALSE),
	50	fITSRefitFlag(kFALSE), fTPCRefitFlag(kFALSE),
aafecd8b	51	fFunctionProbabilityFlag(kFALSE),
	52	fElectronFunction(0), fMuonFunction(0),
	53	fPionFunction(0), fKaonFunction(0), fProtonFunction(0),
734d2c12	54	fHistYPtProtons(0), fHistYPtAntiProtons(0) {
	55	//Default constructor
	56	for(Int_t i = 0; i < 5; i++) fPartFrac[i] = 0.0;
	57	}
	58
	59	//____________________________________________________________________//
	60	AliProtonAnalysis::AliProtonAnalysis(Int_t nbinsY, Float_t fLowY, Float_t fHighY,Int_t nbinsPt, Float_t fLowPt, Float_t fHighPt) :
	61	TObject(),
	62	fNBinsY(nbinsY), fMinY(fLowY), fMaxY(fHighY),
	63	fNBinsPt(nbinsPt), fMinPt(fLowPt), fMaxPt(fHighPt),
aafecd8b	64	fMinTPCClusters(0), fMinITSClusters(0),
	65	fMaxChi2PerTPCCluster(0), fMaxChi2PerITSCluster(0),
	66	fMaxCov11(0), fMaxCov22(0), fMaxCov33(0), fMaxCov44(0), fMaxCov55(0),
	67	fMaxSigmaToVertex(0),
	68	fMinTPCClustersFlag(kFALSE), fMinITSClustersFlag(kFALSE),
	69	fMaxChi2PerTPCClusterFlag(kFALSE), fMaxChi2PerITSClusterFlag(kFALSE),
	70	fMaxCov11Flag(kFALSE), fMaxCov22Flag(kFALSE), fMaxCov33Flag(kFALSE), fMaxCov44Flag(kFALSE), fMaxCov55Flag(kFALSE),
	71	fMaxSigmaToVertexFlag(kFALSE),
	72	fITSRefitFlag(kFALSE), fTPCRefitFlag(kFALSE),
	73	fFunctionProbabilityFlag(kFALSE),
	74	fElectronFunction(0), fMuonFunction(0),
	75	fPionFunction(0), fKaonFunction(0), fProtonFunction(0),
734d2c12	76	fHistYPtProtons(0), fHistYPtAntiProtons(0) {
	77	//Default constructor
	78
	79	fHistYPtProtons = new TH2F("fHistYPtProtons","y-Pt Protons",fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
	80	fHistYPtProtons->SetStats(kTRUE);
	81	fHistYPtProtons->GetYaxis()->SetTitle("P_{T} [GeV]");
	82	fHistYPtProtons->GetXaxis()->SetTitle("y");
	83	fHistYPtProtons->GetXaxis()->SetTitleColor(1);
	84
	85	fHistYPtAntiProtons = new TH2F("fHistYPtAntiProtons","y-Pt Antiprotons",fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
	86	fHistYPtAntiProtons->SetStats(kTRUE);
	87	fHistYPtAntiProtons->GetYaxis()->SetTitle("P_{T} [GeV]");
	88	fHistYPtAntiProtons->GetXaxis()->SetTitle("y");
	89	fHistYPtAntiProtons->GetXaxis()->SetTitleColor(1);
	90	}
	91
	92	//____________________________________________________________________//
	93	AliProtonAnalysis::~AliProtonAnalysis() {
	94	//Default destructor
	95
	96	}
	97
	98	//____________________________________________________________________//
	99	void AliProtonAnalysis::InitHistograms(Int_t nbinsY, Float_t fLowY, Float_t fHighY, Int_t nbinsPt, Float_t fLowPt, Float_t fHighPt) {
	100	fNBinsY = nbinsY;
	101	fMinY = fLowY;
	102	fMaxY = fHighY;
	103	fNBinsPt = nbinsPt;
	104	fMinPt = fLowPt;
	105	fMaxPt = fHighPt;
	106
	107	fHistYPtProtons = new TH2F("fHistYPtProtons","y-Pt Protons",fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
	108	fHistYPtProtons->SetStats(kTRUE);
	109	fHistYPtProtons->GetYaxis()->SetTitle("P_{T} [GeV]");
	110	fHistYPtProtons->GetXaxis()->SetTitle("y");
	111	fHistYPtProtons->GetXaxis()->SetTitleColor(1);
	112
	113	fHistYPtAntiProtons = new TH2F("fHistYPtAntiProtons","y-Pt Antiprotons",fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
	114	fHistYPtAntiProtons->SetStats(kTRUE);
	115	fHistYPtAntiProtons->GetYaxis()->SetTitle("P_{T} [GeV]");
	116	fHistYPtAntiProtons->GetXaxis()->SetTitle("y");
	117	fHistYPtAntiProtons->GetXaxis()->SetTitleColor(1);
	118	}
	119
	120	//____________________________________________________________________//
	121	void AliProtonAnalysis::ReadFromFile(const char* filename) {
	122	TFile *file = TFile::Open(filename);
	123	fHistYPtProtons = (TH2F *)file->Get("fHistYPtProtons");
	124	fHistYPtAntiProtons = (TH2F *)file->Get("fHistYPtAntiProtons");
	125	fHistYPtProtons->Sumw2();
	126	fHistYPtAntiProtons->Sumw2();
	127	}
	128
	129	//____________________________________________________________________//
	130	TH1D *AliProtonAnalysis::GetProtonYHistogram() {
	131	TH1D fYProtons = (TH1D )fHistYPtProtons->ProjectionX("fYProtons",0,fHistYPtProtons->GetYaxis()->GetNbins(),"e");
	132	fYProtons->SetStats(kFALSE);
	133	fYProtons->GetYaxis()->SetTitle("dN/dy");
	134	fYProtons->SetTitle("dN/dy protons");
	135	fYProtons->SetMarkerStyle(kFullCircle);
	136	fYProtons->SetMarkerColor(4);
	137
	138	return fYProtons;
	139	}
140
141	//____________________________________________________________________//
142	TH1D *AliProtonAnalysis::GetAntiProtonYHistogram() {
143	TH1D fYAntiProtons = (TH1D )fHistYPtAntiProtons->ProjectionX("fYAntiProtons",0,fHistYPtAntiProtons->GetYaxis()->GetNbins(),"e");
144	fYAntiProtons->SetStats(kFALSE);
145	fYAntiProtons->GetYaxis()->SetTitle("dN/dy");
146	fYAntiProtons->SetTitle("dN/dy antiprotons");
147	fYAntiProtons->SetMarkerStyle(kFullCircle);
148	fYAntiProtons->SetMarkerColor(4);
149
150	return fYAntiProtons;
151	}
152
153	//____________________________________________________________________//
154	TH1D *AliProtonAnalysis::GetProtonPtHistogram() {
155	TH1D fPtProtons = (TH1D )fHistYPtProtons->ProjectionY("fPtProtons",0,fHistYPtProtons->GetXaxis()->GetNbins(),"e");
156	fPtProtons->SetStats(kFALSE);
157	fPtProtons->GetYaxis()->SetTitle("dN/dP_{T}");
158	fPtProtons->SetTitle("dN/dPt protons");
159	fPtProtons->SetMarkerStyle(kFullCircle);
160	fPtProtons->SetMarkerColor(4);
161
162	return fPtProtons;
163	}
164
165	//____________________________________________________________________//
166	TH1D *AliProtonAnalysis::GetAntiProtonPtHistogram() {
167	TH1D fPtAntiProtons = (TH1D )fHistYPtAntiProtons->ProjectionY("fPtAntiProtons",0,fHistYPtProtons->GetXaxis()->GetNbins(),"e");
168	fPtAntiProtons->SetStats(kFALSE);
169	fPtAntiProtons->GetYaxis()->SetTitle("dN/dP_{T}");
170	fPtAntiProtons->SetTitle("dN/dPt antiprotons");
171	fPtAntiProtons->SetMarkerStyle(kFullCircle);
172	fPtAntiProtons->SetMarkerColor(4);
173
174	return fPtAntiProtons;
175	}
176
177	//____________________________________________________________________//
178	TH1D *AliProtonAnalysis::GetYRatioHistogram() {
179	TH1D *fYProtons = GetProtonYHistogram();
180	TH1D *fYAntiProtons = GetAntiProtonYHistogram();
181
182	TH1D *hRatioY = new TH1D("hRatioY","",fYProtons->GetNbinsX(),fYProtons->GetXaxis()->GetXmin(),fYProtons->GetXaxis()->GetXmax());
183	hRatioY->Divide(fYAntiProtons,fYProtons,1.0,1.0);
184	hRatioY->SetMarkerStyle(kFullCircle);
185	hRatioY->SetMarkerColor(4);
186	hRatioY->GetYaxis()->SetTitle("#bar{p}/p");
187	hRatioY->GetYaxis()->SetTitleOffset(1.4);
188	hRatioY->GetXaxis()->SetTitle("y");
189	hRatioY->GetXaxis()->SetTitleColor(1);
190	hRatioY->SetStats(kFALSE);
191
192	return hRatioY;
193	}
194
195	//____________________________________________________________________//
196	TH1D *AliProtonAnalysis::GetPtRatioHistogram() {
197	TH1D *fPtProtons = GetProtonPtHistogram();
198	TH1D *fPtAntiProtons = GetAntiProtonPtHistogram();
199
200	TH1D *hRatioPt = new TH1D("hRatioPt","",fPtProtons->GetNbinsX(),fPtProtons->GetXaxis()->GetXmin(),fPtProtons->GetXaxis()->GetXmax());
201	hRatioPt->Divide(fPtAntiProtons,fPtProtons,1.0,1.0);
202	hRatioPt->SetMarkerStyle(kFullCircle);
203	hRatioPt->SetMarkerColor(4);
204	hRatioPt->GetYaxis()->SetTitle("#bar{p}/p");
205	hRatioPt->GetYaxis()->SetTitleOffset(1.4);
206	hRatioPt->GetXaxis()->SetTitle("P_{T} [GeV/c]");
207	hRatioPt->GetXaxis()->SetTitleColor(1);
208	hRatioPt->SetStats(kFALSE);
209
210	return hRatioPt;
211	}
212
213	//____________________________________________________________________//
214	TH1D *AliProtonAnalysis::GetYAsymmetryHistogram() {
215	TH1D *fYProtons = GetProtonYHistogram();
216	TH1D *fYAntiProtons = GetAntiProtonYHistogram();
217
218	TH1D *hsum = new TH1D("hsumY","",fYProtons->GetNbinsX(),fYProtons->GetXaxis()->GetXmin(),fYProtons->GetXaxis()->GetXmax());
219	hsum->Add(fYProtons,fYAntiProtons,1.0,1.0);
220
221	TH1D *hdiff = new TH1D("hdiffY","",fYProtons->GetNbinsX(),fYProtons->GetXaxis()->GetXmin(),fYProtons->GetXaxis()->GetXmax());
222	hdiff->Add(fYProtons,fYAntiProtons,1.0,-1.0);
223
224	TH1D *hAsymmetryY = new TH1D("hAsymmetryY","",fYProtons->GetNbinsX(),fYProtons->GetXaxis()->GetXmin(),fYProtons->GetXaxis()->GetXmax());
225	hAsymmetryY->Divide(hdiff,hsum,2.0,1.);
226	hAsymmetryY->SetMarkerStyle(kFullCircle);
227	hAsymmetryY->SetMarkerColor(4);
228	hAsymmetryY->GetYaxis()->SetTitle("A_{p}");
229	hAsymmetryY->GetYaxis()->SetTitleOffset(1.4);
230	hAsymmetryY->GetXaxis()->SetTitle("y");
231	hAsymmetryY->GetXaxis()->SetTitleColor(1);
232	hAsymmetryY->SetStats(kFALSE);
233
234	return hAsymmetryY;
235	}
236
237	//____________________________________________________________________//
238	TH1D *AliProtonAnalysis::GetPtAsymmetryHistogram() {
239	TH1D *fPtProtons = GetProtonPtHistogram();
240	TH1D *fPtAntiProtons = GetAntiProtonPtHistogram();
241
242	TH1D *hsum = new TH1D("hsumPt","",fPtProtons->GetNbinsX(),fPtProtons->GetXaxis()->GetXmin(),fPtProtons->GetXaxis()->GetXmax());
243	hsum->Add(fPtProtons,fPtAntiProtons,1.0,1.0);
244
245	TH1D *hdiff = new TH1D("hdiffPt","",fPtProtons->GetNbinsX(),fPtProtons->GetXaxis()->GetXmin(),fPtProtons->GetXaxis()->GetXmax());
246	hdiff->Add(fPtProtons,fPtAntiProtons,1.0,-1.0);
247
248	TH1D *hAsymmetryPt = new TH1D("hAsymmetryPt","",fPtProtons->GetNbinsX(),fPtProtons->GetXaxis()->GetXmin(),fPtProtons->GetXaxis()->GetXmax());
249	hAsymmetryPt->Divide(hdiff,hsum,2.0,1.);
250	hAsymmetryPt->SetMarkerStyle(kFullCircle);
251	hAsymmetryPt->SetMarkerColor(4);
252	hAsymmetryPt->GetYaxis()->SetTitle("A_{p}");
253	hAsymmetryPt->GetYaxis()->SetTitleOffset(1.4);
254	hAsymmetryPt->GetXaxis()->SetTitle("P_{T} [GeV/c]");
255	hAsymmetryPt->GetXaxis()->SetTitleColor(1);
256	hAsymmetryPt->SetStats(kFALSE);
257
258	return hAsymmetryPt;
259	}
260
aafecd8b	261	//____________________________________________________________________//
aafecd8b	262	Double_t AliProtonAnalysis::GetParticleFraction(Int_t i, Double_t p) {
8b8b0b7a	263	Double_t partFrac=0;
aafecd8b	264	if(fFunctionProbabilityFlag) {
	265	if(i == 0) partFrac = fElectronFunction->Eval(p);
	266	if(i == 1) partFrac = fMuonFunction->Eval(p);
	267	if(i == 2) partFrac = fPionFunction->Eval(p);
	268	if(i == 3) partFrac = fKaonFunction->Eval(p);
	269	if(i == 4) partFrac = fProtonFunction->Eval(p);
	270	}
	271	else partFrac = fPartFrac[i];
	272
	273	return partFrac;
	274	}
	275
734d2c12	276	//____________________________________________________________________//
	277	void AliProtonAnalysis::Analyze(AliESDEvent* fESD) {
	278	//Main analysis part
	279	Int_t nGoodTracks = fESD->GetNumberOfTracks();
	280	for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) {
	281	AliESDtrack* track = fESD->GetTrack(iTracks);
	282	if(IsAccepted(track)) {
	283	Double_t Pt = track->Pt();
aafecd8b	284	Double_t P = track->P();
734d2c12	285
aafecd8b	286	//pid
734d2c12	287	Double_t probability[5];
	288	track->GetESDpid(probability);
	289	Double_t rcc = 0.0;
aafecd8b	290	for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += probability[i]*GetParticleFraction(i,P);
734d2c12	291	if(rcc == 0.0) continue;
734d2c12	292	Double_t w[5];
aafecd8b	293	for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = probability[i]*GetParticleFraction(i,P)/rcc;
734d2c12	294	Long64_t fParticleType = TMath::LocMax(AliPID::kSPECIES,w);
	295	if(fParticleType == 4) {
	296	if(track->Charge() > 0) fHistYPtProtons->Fill(Rapidity(track),Pt);
	297	else if(track->Charge() < 0) fHistYPtAntiProtons->Fill(Rapidity(track),Pt);
	298	}//proton check
	299	}//cuts
	300	}//track loop
	301	}
	302
ee4ca40d	303	//____________________________________________________________________//
	304	void AliProtonAnalysis::Analyze(AliAODEvent* fAOD) {
	305	//Main analysis part
	306	Int_t nGoodTracks = fAOD->GetNumberOfTracks();
	307	for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) {
	308	AliAODTrack* track = fAOD->GetTrack(iTracks);
	309	Double_t Pt = track->Pt();
	310	Double_t P = track->P();
	311
	312	//pid
	313	Double_t probability[10];
	314	track->GetPID(probability);
	315	Double_t rcc = 0.0;
	316	for(Int_t i = 0; i < AliPID::kSPECIESN; i++) rcc += probability[i]*GetParticleFraction(i,P);
	317	if(rcc == 0.0) continue;
738619fd	318	Double_t w[10];
ee4ca40d	319	for(Int_t i = 0; i < AliPID::kSPECIESN; i++) w[i] = probability[i]*GetParticleFraction(i,P)/rcc;
	320	Long64_t fParticleType = TMath::LocMax(AliPID::kSPECIESN,w);
	321	if(fParticleType == 4) {
	322	if(track->Charge() > 0) fHistYPtProtons->Fill(track->Y(fParticleType),Pt);
	323	else if(track->Charge() < 0) fHistYPtAntiProtons->Fill(track->Y(fParticleType),Pt);
	324	}//proton check
	325	}//track loop
	326	}
	327
734d2c12	328	//____________________________________________________________________//
	329	Bool_t AliProtonAnalysis::IsAccepted(AliESDtrack* track) {
	330	// Checks if the track is excluded from the cuts
	331	Int_t fIdxInt[200];
	332	Int_t nClustersITS = track->GetITSclusters(fIdxInt);
	333	Int_t nClustersTPC = track->GetTPCclusters(fIdxInt);
	334
	335	Float_t chi2PerClusterITS = -1;
	336	Float_t chi2PerClusterTPC = -1;
	337	if (nClustersTPC!=0)
	338	chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC);
	339
	340	Double_t extCov[15];
	341	track->GetExternalCovariance(extCov);
	342
	343	Double_t Pt = track->Pt();
734d2c12	344
a5375b97	345	if(fMinITSClustersFlag)
734d2c12	346	if(nClustersITS < fMinITSClusters) return kFALSE;
	347	if(fMinTPCClustersFlag)
	348	if(nClustersTPC < fMinTPCClusters) return kFALSE;
	349	if(fMaxChi2PerTPCClusterFlag)
	350	if(chi2PerClusterTPC > fMaxChi2PerTPCCluster) return kFALSE;
	351	if(fMaxChi2PerITSClusterFlag)
	352	if(chi2PerClusterITS > fMaxChi2PerITSCluster) return kFALSE;
	353	if(fMaxCov11Flag)
	354	if(extCov[0] > fMaxCov11) return kFALSE;
	355	if(fMaxCov22Flag)
	356	if(extCov[2] > fMaxCov22) return kFALSE;
	357	if(fMaxCov33Flag)
	358	if(extCov[5] > fMaxCov33) return kFALSE;
	359	if(fMaxCov44Flag)
	360	if(extCov[9] > fMaxCov44) return kFALSE;
	361	if(fMaxCov55Flag)
	362	if(extCov[14] > fMaxCov55) return kFALSE;
	363	if(fMaxSigmaToVertexFlag)
	364	if(GetSigmaToVertex(track) > fMaxSigmaToVertex) return kFALSE;
	365	if(fITSRefitFlag)
	366	if ((track->GetStatus() & AliESDtrack::kITSrefit) == 0) return kFALSE;
	367	if(fTPCRefitFlag)
	368	if ((track->GetStatus() & AliESDtrack::kTPCrefit) == 0) return kFALSE;
	369
	370	if((Pt < fMinPt) \|\| (Pt > fMaxPt)) return kFALSE;
	371
	372	return kTRUE;
	373	}
	374
	375	//____________________________________________________________________//
	376	Float_t AliProtonAnalysis::GetSigmaToVertex(AliESDtrack* esdTrack) {
	377	// Calculates the number of sigma to the vertex.
	378
	379	Float_t b[2];
	380	Float_t bRes[2];
	381	Float_t bCov[3];
	382	esdTrack->GetImpactParameters(b,bCov);
	383	if (bCov[0]<=0 \|\| bCov[2]<=0) {
	384	//AliDebug(1, "Estimated b resolution lower or equal zero!");
	385	bCov[0]=0; bCov[2]=0;
	386	}
	387	bRes[0] = TMath::Sqrt(bCov[0]);
	388	bRes[1] = TMath::Sqrt(bCov[2]);
	389
	390	if (bRes[0] == 0 \|\| bRes[1] ==0) return -1;
	391
	392	Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
	393
	394	if (TMath::Exp(-d * d / 2) < 1e-10) return 1000;
	395
	396	d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
	397
	398	return d;
	399	}
	400
	401	Double_t AliProtonAnalysis::Rapidity(AliESDtrack *track) {
	402	Double_t fMass = 9.38270000000000048e-01;
	403
	404	Double_t P = TMath::Sqrt(TMath::Power(track->Px(),2) +
	405	TMath::Power(track->Py(),2) +
	406	TMath::Power(track->Pz(),2));
	407	Double_t energy = TMath::Sqrt(PP + fMassfMass);
	408	Double_t y = -999;
	409	if(energy != track->Pz())
410	y = 0.5*TMath::Log((energy + track->Pz())/(energy - track->Pz()));
411
412	return y;
413	}