AOD analysis for protons - analysis train
[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>
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>
32#include <AliLog.h>
ee4ca40d 33#include <AliPID.h>
734d2c12 34
35ClassImp(AliProtonAnalysis)
36
37//____________________________________________________________________//
38AliProtonAnalysis::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//____________________________________________________________________//
60AliProtonAnalysis::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//____________________________________________________________________//
93AliProtonAnalysis::~AliProtonAnalysis() {
94 //Default destructor
95
96}
97
98//____________________________________________________________________//
99void 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//____________________________________________________________________//
121void 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//____________________________________________________________________//
130TH1D *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//____________________________________________________________________//
142TH1D *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//____________________________________________________________________//
154TH1D *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//____________________________________________________________________//
166TH1D *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//____________________________________________________________________//
178TH1D *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//____________________________________________________________________//
196TH1D *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//____________________________________________________________________//
214TH1D *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//____________________________________________________________________//
238TH1D *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//____________________________________________________________________//
262Double_t AliProtonAnalysis::GetParticleFraction(Int_t i, Double_t p) {
263 Double_t partFrac;
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//____________________________________________________________________//
277void 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;
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//____________________________________________________________________//
304void 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;
318 Double_t w[5];
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//____________________________________________________________________//
329Bool_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//____________________________________________________________________//
376Float_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
401Double_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(P*P + fMass*fMass);
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}