[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 <TH2F.h>
#include <TH1D.h>

#include "AliProtonAnalysis.h"

#include <AliESDEvent.h>
#include <AliLog.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),
  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),
  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;
}

//____________________________________________________________________//
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();
	
	  //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]*fPartFrac[i];
      if(rcc == 0.0) continue;
      Double_t w[5];
      for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = probability[i]*fPartFrac[i]/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 
}

//____________________________________________________________________//
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>
	24	#include <TH2F.h>
	25	#include <TH1D.h>
	26
	27	#include "AliProtonAnalysis.h"
	28
	29	#include <AliESDEvent.h>
	30	#include <AliLog.h>
	31
	32	ClassImp(AliProtonAnalysis)
	33
	34	//____________________________________________________________________//
	35	AliProtonAnalysis::AliProtonAnalysis() :
	36	TObject(),
	37	fNBinsY(0), fMinY(0), fMaxY(0),
	38	fNBinsPt(0), fMinPt(0), fMaxPt(0),
	39	fMinTPCClusters(0), fMinITSClusters(0),
	40	fMaxChi2PerTPCCluster(0), fMaxChi2PerITSCluster(0),
	41	fMaxCov11(0), fMaxCov22(0), fMaxCov33(0), fMaxCov44(0), fMaxCov55(0),
	42	fMaxSigmaToVertex(0),
	43	fMinTPCClustersFlag(kFALSE), fMinITSClustersFlag(kFALSE),
	44	fMaxChi2PerTPCClusterFlag(kFALSE), fMaxChi2PerITSClusterFlag(kFALSE),
	45	fMaxCov11Flag(kFALSE), fMaxCov22Flag(kFALSE), fMaxCov33Flag(kFALSE), fMaxCov44Flag(kFALSE), fMaxCov55Flag(kFALSE),
	46	fMaxSigmaToVertexFlag(kFALSE),
	47	fITSRefitFlag(kFALSE), fTPCRefitFlag(kFALSE),
	48	fHistYPtProtons(0), fHistYPtAntiProtons(0) {
	49	//Default constructor
	50	for(Int_t i = 0; i < 5; i++) fPartFrac[i] = 0.0;
	51	}
	52
	53	//____________________________________________________________________//
	54	AliProtonAnalysis::AliProtonAnalysis(Int_t nbinsY, Float_t fLowY, Float_t fHighY,Int_t nbinsPt, Float_t fLowPt, Float_t fHighPt) :
	55	TObject(),
	56	fNBinsY(nbinsY), fMinY(fLowY), fMaxY(fHighY),
	57	fNBinsPt(nbinsPt), fMinPt(fLowPt), fMaxPt(fHighPt),
	58	fHistYPtProtons(0), fHistYPtAntiProtons(0) {
	59	//Default constructor
	60
	61	fHistYPtProtons = new TH2F("fHistYPtProtons","y-Pt Protons",fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
	62	fHistYPtProtons->SetStats(kTRUE);
	63	fHistYPtProtons->GetYaxis()->SetTitle("P_{T} [GeV]");
	64	fHistYPtProtons->GetXaxis()->SetTitle("y");
65	fHistYPtProtons->GetXaxis()->SetTitleColor(1);
66
67	fHistYPtAntiProtons = new TH2F("fHistYPtAntiProtons","y-Pt Antiprotons",fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
68	fHistYPtAntiProtons->SetStats(kTRUE);
69	fHistYPtAntiProtons->GetYaxis()->SetTitle("P_{T} [GeV]");
70	fHistYPtAntiProtons->GetXaxis()->SetTitle("y");
71	fHistYPtAntiProtons->GetXaxis()->SetTitleColor(1);
72	}
73
74	//____________________________________________________________________//
75	AliProtonAnalysis::~AliProtonAnalysis() {
76	//Default destructor
77
78	}
79
80	//____________________________________________________________________//
81	void AliProtonAnalysis::InitHistograms(Int_t nbinsY, Float_t fLowY, Float_t fHighY, Int_t nbinsPt, Float_t fLowPt, Float_t fHighPt) {
82	fNBinsY = nbinsY;
83	fMinY = fLowY;
84	fMaxY = fHighY;
85	fNBinsPt = nbinsPt;
86	fMinPt = fLowPt;
87	fMaxPt = fHighPt;
88
89	fHistYPtProtons = new TH2F("fHistYPtProtons","y-Pt Protons",fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
90	fHistYPtProtons->SetStats(kTRUE);
91	fHistYPtProtons->GetYaxis()->SetTitle("P_{T} [GeV]");
92	fHistYPtProtons->GetXaxis()->SetTitle("y");
93	fHistYPtProtons->GetXaxis()->SetTitleColor(1);
94
95	fHistYPtAntiProtons = new TH2F("fHistYPtAntiProtons","y-Pt Antiprotons",fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
96	fHistYPtAntiProtons->SetStats(kTRUE);
97	fHistYPtAntiProtons->GetYaxis()->SetTitle("P_{T} [GeV]");
98	fHistYPtAntiProtons->GetXaxis()->SetTitle("y");
99	fHistYPtAntiProtons->GetXaxis()->SetTitleColor(1);
100	}
101
102	//____________________________________________________________________//
103	void AliProtonAnalysis::ReadFromFile(const char* filename) {
104	TFile *file = TFile::Open(filename);
105	fHistYPtProtons = (TH2F *)file->Get("fHistYPtProtons");
106	fHistYPtAntiProtons = (TH2F *)file->Get("fHistYPtAntiProtons");
107	fHistYPtProtons->Sumw2();
108	fHistYPtAntiProtons->Sumw2();
109	}
110
111	//____________________________________________________________________//
112	TH1D *AliProtonAnalysis::GetProtonYHistogram() {
113	TH1D fYProtons = (TH1D )fHistYPtProtons->ProjectionX("fYProtons",0,fHistYPtProtons->GetYaxis()->GetNbins(),"e");
114	fYProtons->SetStats(kFALSE);
115	fYProtons->GetYaxis()->SetTitle("dN/dy");
116	fYProtons->SetTitle("dN/dy protons");
117	fYProtons->SetMarkerStyle(kFullCircle);
118	fYProtons->SetMarkerColor(4);
119
120	return fYProtons;
121	}
122
123	//____________________________________________________________________//
124	TH1D *AliProtonAnalysis::GetAntiProtonYHistogram() {
125	TH1D fYAntiProtons = (TH1D )fHistYPtAntiProtons->ProjectionX("fYAntiProtons",0,fHistYPtAntiProtons->GetYaxis()->GetNbins(),"e");
126	fYAntiProtons->SetStats(kFALSE);
127	fYAntiProtons->GetYaxis()->SetTitle("dN/dy");
128	fYAntiProtons->SetTitle("dN/dy antiprotons");
129	fYAntiProtons->SetMarkerStyle(kFullCircle);
130	fYAntiProtons->SetMarkerColor(4);
131
132	return fYAntiProtons;
133	}
134
135	//____________________________________________________________________//
136	TH1D *AliProtonAnalysis::GetProtonPtHistogram() {
137	TH1D fPtProtons = (TH1D )fHistYPtProtons->ProjectionY("fPtProtons",0,fHistYPtProtons->GetXaxis()->GetNbins(),"e");
138	fPtProtons->SetStats(kFALSE);
139	fPtProtons->GetYaxis()->SetTitle("dN/dP_{T}");
140	fPtProtons->SetTitle("dN/dPt protons");
141	fPtProtons->SetMarkerStyle(kFullCircle);
142	fPtProtons->SetMarkerColor(4);
143
144	return fPtProtons;
145	}
146
147	//____________________________________________________________________//
148	TH1D *AliProtonAnalysis::GetAntiProtonPtHistogram() {
149	TH1D fPtAntiProtons = (TH1D )fHistYPtAntiProtons->ProjectionY("fPtAntiProtons",0,fHistYPtProtons->GetXaxis()->GetNbins(),"e");
150	fPtAntiProtons->SetStats(kFALSE);
151	fPtAntiProtons->GetYaxis()->SetTitle("dN/dP_{T}");
152	fPtAntiProtons->SetTitle("dN/dPt antiprotons");
153	fPtAntiProtons->SetMarkerStyle(kFullCircle);
154	fPtAntiProtons->SetMarkerColor(4);
155
156	return fPtAntiProtons;
157	}
158
159	//____________________________________________________________________//
160	TH1D *AliProtonAnalysis::GetYRatioHistogram() {
161	TH1D *fYProtons = GetProtonYHistogram();
162	TH1D *fYAntiProtons = GetAntiProtonYHistogram();
163
164	TH1D *hRatioY = new TH1D("hRatioY","",fYProtons->GetNbinsX(),fYProtons->GetXaxis()->GetXmin(),fYProtons->GetXaxis()->GetXmax());
165	hRatioY->Divide(fYAntiProtons,fYProtons,1.0,1.0);
166	hRatioY->SetMarkerStyle(kFullCircle);
167	hRatioY->SetMarkerColor(4);
168	hRatioY->GetYaxis()->SetTitle("#bar{p}/p");
169	hRatioY->GetYaxis()->SetTitleOffset(1.4);
170	hRatioY->GetXaxis()->SetTitle("y");
171	hRatioY->GetXaxis()->SetTitleColor(1);
172	hRatioY->SetStats(kFALSE);
173
174	return hRatioY;
175	}
176
177	//____________________________________________________________________//
178	TH1D *AliProtonAnalysis::GetPtRatioHistogram() {
179	TH1D *fPtProtons = GetProtonPtHistogram();
180	TH1D *fPtAntiProtons = GetAntiProtonPtHistogram();
181
182	TH1D *hRatioPt = new TH1D("hRatioPt","",fPtProtons->GetNbinsX(),fPtProtons->GetXaxis()->GetXmin(),fPtProtons->GetXaxis()->GetXmax());
183	hRatioPt->Divide(fPtAntiProtons,fPtProtons,1.0,1.0);
184	hRatioPt->SetMarkerStyle(kFullCircle);
185	hRatioPt->SetMarkerColor(4);
186	hRatioPt->GetYaxis()->SetTitle("#bar{p}/p");
187	hRatioPt->GetYaxis()->SetTitleOffset(1.4);
188	hRatioPt->GetXaxis()->SetTitle("P_{T} [GeV/c]");
189	hRatioPt->GetXaxis()->SetTitleColor(1);
190	hRatioPt->SetStats(kFALSE);
191
192	return hRatioPt;
193	}
194
195	//____________________________________________________________________//
196	TH1D *AliProtonAnalysis::GetYAsymmetryHistogram() {
197	TH1D *fYProtons = GetProtonYHistogram();
198	TH1D *fYAntiProtons = GetAntiProtonYHistogram();
199
200	TH1D *hsum = new TH1D("hsumY","",fYProtons->GetNbinsX(),fYProtons->GetXaxis()->GetXmin(),fYProtons->GetXaxis()->GetXmax());
201	hsum->Add(fYProtons,fYAntiProtons,1.0,1.0);
202
203	TH1D *hdiff = new TH1D("hdiffY","",fYProtons->GetNbinsX(),fYProtons->GetXaxis()->GetXmin(),fYProtons->GetXaxis()->GetXmax());
204	hdiff->Add(fYProtons,fYAntiProtons,1.0,-1.0);
205
206	TH1D *hAsymmetryY = new TH1D("hAsymmetryY","",fYProtons->GetNbinsX(),fYProtons->GetXaxis()->GetXmin(),fYProtons->GetXaxis()->GetXmax());
207	hAsymmetryY->Divide(hdiff,hsum,2.0,1.);
208	hAsymmetryY->SetMarkerStyle(kFullCircle);
209	hAsymmetryY->SetMarkerColor(4);
210	hAsymmetryY->GetYaxis()->SetTitle("A_{p}");
211	hAsymmetryY->GetYaxis()->SetTitleOffset(1.4);
212	hAsymmetryY->GetXaxis()->SetTitle("y");
213	hAsymmetryY->GetXaxis()->SetTitleColor(1);
214	hAsymmetryY->SetStats(kFALSE);
215
216	return hAsymmetryY;
217	}
218
219	//____________________________________________________________________//
220	TH1D *AliProtonAnalysis::GetPtAsymmetryHistogram() {
221	TH1D *fPtProtons = GetProtonPtHistogram();
222	TH1D *fPtAntiProtons = GetAntiProtonPtHistogram();
223
224	TH1D *hsum = new TH1D("hsumPt","",fPtProtons->GetNbinsX(),fPtProtons->GetXaxis()->GetXmin(),fPtProtons->GetXaxis()->GetXmax());
225	hsum->Add(fPtProtons,fPtAntiProtons,1.0,1.0);
226
227	TH1D *hdiff = new TH1D("hdiffPt","",fPtProtons->GetNbinsX(),fPtProtons->GetXaxis()->GetXmin(),fPtProtons->GetXaxis()->GetXmax());
228	hdiff->Add(fPtProtons,fPtAntiProtons,1.0,-1.0);
229
230	TH1D *hAsymmetryPt = new TH1D("hAsymmetryPt","",fPtProtons->GetNbinsX(),fPtProtons->GetXaxis()->GetXmin(),fPtProtons->GetXaxis()->GetXmax());
231	hAsymmetryPt->Divide(hdiff,hsum,2.0,1.);
232	hAsymmetryPt->SetMarkerStyle(kFullCircle);
233	hAsymmetryPt->SetMarkerColor(4);
234	hAsymmetryPt->GetYaxis()->SetTitle("A_{p}");
235	hAsymmetryPt->GetYaxis()->SetTitleOffset(1.4);
236	hAsymmetryPt->GetXaxis()->SetTitle("P_{T} [GeV/c]");
237	hAsymmetryPt->GetXaxis()->SetTitleColor(1);
238	hAsymmetryPt->SetStats(kFALSE);
239
240	return hAsymmetryPt;
241	}
242
243	//____________________________________________________________________//
244	void AliProtonAnalysis::Analyze(AliESDEvent* fESD) {
245	//Main analysis part
246	Int_t nGoodTracks = fESD->GetNumberOfTracks();
247	for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) {
248	AliESDtrack* track = fESD->GetTrack(iTracks);
249	if(IsAccepted(track)) {
250	Double_t Pt = track->Pt();
251
252	//pid
253	Double_t probability[5];
254	track->GetESDpid(probability);
255	Double_t rcc = 0.0;
256	for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += probability[i]*fPartFrac[i];
257	if(rcc == 0.0) continue;
258	Double_t w[5];
259	for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = probability[i]*fPartFrac[i]/rcc;
260	Long64_t fParticleType = TMath::LocMax(AliPID::kSPECIES,w);
261	if(fParticleType == 4) {
262	if(track->Charge() > 0) fHistYPtProtons->Fill(Rapidity(track),Pt);
263	else if(track->Charge() < 0) fHistYPtAntiProtons->Fill(Rapidity(track),Pt);
264	}//proton check
265	}//cuts
266	}//track loop
267	}
268
269	//____________________________________________________________________//
270	Bool_t AliProtonAnalysis::IsAccepted(AliESDtrack* track) {
271	// Checks if the track is excluded from the cuts
272	Int_t fIdxInt[200];
273	Int_t nClustersITS = track->GetITSclusters(fIdxInt);
274	Int_t nClustersTPC = track->GetTPCclusters(fIdxInt);
275
276	Float_t chi2PerClusterITS = -1;
277	Float_t chi2PerClusterTPC = -1;
278	if (nClustersTPC!=0)
279	chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC);
280
281	Double_t extCov[15];
282	track->GetExternalCovariance(extCov);
283
284	Double_t Pt = track->Pt();
734d2c12	285
a5375b97	286	if(fMinITSClustersFlag)
734d2c12	287	if(nClustersITS < fMinITSClusters) return kFALSE;
	288	if(fMinTPCClustersFlag)
	289	if(nClustersTPC < fMinTPCClusters) return kFALSE;
	290	if(fMaxChi2PerTPCClusterFlag)
	291	if(chi2PerClusterTPC > fMaxChi2PerTPCCluster) return kFALSE;
	292	if(fMaxChi2PerITSClusterFlag)
	293	if(chi2PerClusterITS > fMaxChi2PerITSCluster) return kFALSE;
	294	if(fMaxCov11Flag)
	295	if(extCov[0] > fMaxCov11) return kFALSE;
	296	if(fMaxCov22Flag)
	297	if(extCov[2] > fMaxCov22) return kFALSE;
	298	if(fMaxCov33Flag)
	299	if(extCov[5] > fMaxCov33) return kFALSE;
	300	if(fMaxCov44Flag)
	301	if(extCov[9] > fMaxCov44) return kFALSE;
	302	if(fMaxCov55Flag)
	303	if(extCov[14] > fMaxCov55) return kFALSE;
	304	if(fMaxSigmaToVertexFlag)
	305	if(GetSigmaToVertex(track) > fMaxSigmaToVertex) return kFALSE;
	306	if(fITSRefitFlag)
	307	if ((track->GetStatus() & AliESDtrack::kITSrefit) == 0) return kFALSE;
	308	if(fTPCRefitFlag)
	309	if ((track->GetStatus() & AliESDtrack::kTPCrefit) == 0) return kFALSE;
	310
	311	if((Pt < fMinPt) \|\| (Pt > fMaxPt)) return kFALSE;
	312
	313	return kTRUE;
	314	}
	315
	316	//____________________________________________________________________//
	317	Float_t AliProtonAnalysis::GetSigmaToVertex(AliESDtrack* esdTrack) {
	318	// Calculates the number of sigma to the vertex.
	319
	320	Float_t b[2];
	321	Float_t bRes[2];
	322	Float_t bCov[3];
	323	esdTrack->GetImpactParameters(b,bCov);
	324	if (bCov[0]<=0 \|\| bCov[2]<=0) {
	325	//AliDebug(1, "Estimated b resolution lower or equal zero!");
	326	bCov[0]=0; bCov[2]=0;
	327	}
	328	bRes[0] = TMath::Sqrt(bCov[0]);
	329	bRes[1] = TMath::Sqrt(bCov[2]);
	330
	331	if (bRes[0] == 0 \|\| bRes[1] ==0) return -1;
	332
	333	Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
	334
	335	if (TMath::Exp(-d * d / 2) < 1e-10) return 1000;
	336
	337	d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
	338
	339	return d;
	340	}
	341
	342	Double_t AliProtonAnalysis::Rapidity(AliESDtrack *track) {
	343	Double_t fMass = 9.38270000000000048e-01;
	344
	345	Double_t P = TMath::Sqrt(TMath::Power(track->Px(),2) +
	346	TMath::Power(track->Py(),2) +
	347	TMath::Power(track->Pz(),2));
	348	Double_t energy = TMath::Sqrt(PP + fMassfMass);
	349	Double_t y = -999;
	350	if(energy != track->Pz())
351	y = 0.5*TMath::Log((energy + track->Pz())/(energy - track->Pz()));
352
353	return y;
354	}