Removing not used parameter and putting the correct flag for the ITS clusters
[u/mrichter/AliRoot.git] / PWG2 / SPECTRA / AliProtonAnalysis.cxx
CommitLineData
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
32ClassImp(AliProtonAnalysis)
33
34//____________________________________________________________________//
35AliProtonAnalysis::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//____________________________________________________________________//
54AliProtonAnalysis::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//____________________________________________________________________//
75AliProtonAnalysis::~AliProtonAnalysis() {
76 //Default destructor
77
78}
79
80//____________________________________________________________________//
81void 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//____________________________________________________________________//
103void 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//____________________________________________________________________//
112TH1D *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//____________________________________________________________________//
124TH1D *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//____________________________________________________________________//
136TH1D *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//____________________________________________________________________//
148TH1D *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//____________________________________________________________________//
160TH1D *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//____________________________________________________________________//
178TH1D *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//____________________________________________________________________//
196TH1D *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//____________________________________________________________________//
220TH1D *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//____________________________________________________________________//
244void 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//____________________________________________________________________//
270Bool_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//____________________________________________________________________//
317Float_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
342Double_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(P*P + fMass*fMass);
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}