]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG2/SPECTRA/AliProtonAnalysis.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Making GetInputBlock method const.
[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>
3f6d0c08 27#include <TH1I.h>
e4358d7f 28#include <TParticle.h>
734d2c12 29
30#include "AliProtonAnalysis.h"
31
2b748670 32#include <AliExternalTrackParam.h>
ee4ca40d 33#include <AliAODEvent.h>
734d2c12 34#include <AliESDEvent.h>
35#include <AliLog.h>
ee4ca40d 36#include <AliPID.h>
e4358d7f 37#include <AliStack.h>
39f2a708 38#include <AliCFContainer.h>
39#include <AliCFEffGrid.h>
40#include <AliCFEffGrid.h>
e4358d7f 41
734d2c12 42ClassImp(AliProtonAnalysis)
43
44//____________________________________________________________________//
45AliProtonAnalysis::AliProtonAnalysis() :
46 TObject(),
47 fNBinsY(0), fMinY(0), fMaxY(0),
48 fNBinsPt(0), fMinPt(0), fMaxPt(0),
49 fMinTPCClusters(0), fMinITSClusters(0),
50 fMaxChi2PerTPCCluster(0), fMaxChi2PerITSCluster(0),
51 fMaxCov11(0), fMaxCov22(0), fMaxCov33(0), fMaxCov44(0), fMaxCov55(0),
52 fMaxSigmaToVertex(0),
53 fMinTPCClustersFlag(kFALSE), fMinITSClustersFlag(kFALSE),
54 fMaxChi2PerTPCClusterFlag(kFALSE), fMaxChi2PerITSClusterFlag(kFALSE),
55 fMaxCov11Flag(kFALSE), fMaxCov22Flag(kFALSE), fMaxCov33Flag(kFALSE), fMaxCov44Flag(kFALSE), fMaxCov55Flag(kFALSE),
56 fMaxSigmaToVertexFlag(kFALSE),
57 fITSRefitFlag(kFALSE), fTPCRefitFlag(kFALSE),
aafecd8b 58 fFunctionProbabilityFlag(kFALSE),
59 fElectronFunction(0), fMuonFunction(0),
60 fPionFunction(0), fKaonFunction(0), fProtonFunction(0),
39f2a708 61 fUseTPCOnly(kFALSE), fHistEvents(0), fHistYPtProtons(0), fHistYPtAntiProtons(0),
62 fCorrectionList2D(0), fEfficiencyList1D(0), fCorrectionList1D(0) {
734d2c12 63 //Default constructor
64 for(Int_t i = 0; i < 5; i++) fPartFrac[i] = 0.0;
39f2a708 65 fCorrectionList2D = new TList();
66 fEfficiencyList1D = new TList();
67 fCorrectionList1D = new TList();
734d2c12 68}
69
70//____________________________________________________________________//
71AliProtonAnalysis::AliProtonAnalysis(Int_t nbinsY, Float_t fLowY, Float_t fHighY,Int_t nbinsPt, Float_t fLowPt, Float_t fHighPt) :
72 TObject(),
73 fNBinsY(nbinsY), fMinY(fLowY), fMaxY(fHighY),
74 fNBinsPt(nbinsPt), fMinPt(fLowPt), fMaxPt(fHighPt),
aafecd8b 75 fMinTPCClusters(0), fMinITSClusters(0),
76 fMaxChi2PerTPCCluster(0), fMaxChi2PerITSCluster(0),
77 fMaxCov11(0), fMaxCov22(0), fMaxCov33(0), fMaxCov44(0), fMaxCov55(0),
78 fMaxSigmaToVertex(0),
79 fMinTPCClustersFlag(kFALSE), fMinITSClustersFlag(kFALSE),
80 fMaxChi2PerTPCClusterFlag(kFALSE), fMaxChi2PerITSClusterFlag(kFALSE),
81 fMaxCov11Flag(kFALSE), fMaxCov22Flag(kFALSE), fMaxCov33Flag(kFALSE), fMaxCov44Flag(kFALSE), fMaxCov55Flag(kFALSE),
82 fMaxSigmaToVertexFlag(kFALSE),
83 fITSRefitFlag(kFALSE), fTPCRefitFlag(kFALSE),
84 fFunctionProbabilityFlag(kFALSE),
85 fElectronFunction(0), fMuonFunction(0),
86 fPionFunction(0), fKaonFunction(0), fProtonFunction(0),
39f2a708 87 fUseTPCOnly(kFALSE), fHistEvents(0), fHistYPtProtons(0), fHistYPtAntiProtons(0),
88 fCorrectionList2D(0), fEfficiencyList1D(0), fCorrectionList1D(0) {
734d2c12 89 //Default constructor
90
3f6d0c08 91 fHistEvents = new TH1I("fHistEvents","Analyzed events",1,0,1);
92
734d2c12 93 fHistYPtProtons = new TH2F("fHistYPtProtons","y-Pt Protons",fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
94 fHistYPtProtons->SetStats(kTRUE);
95 fHistYPtProtons->GetYaxis()->SetTitle("P_{T} [GeV]");
96 fHistYPtProtons->GetXaxis()->SetTitle("y");
97 fHistYPtProtons->GetXaxis()->SetTitleColor(1);
98
99 fHistYPtAntiProtons = new TH2F("fHistYPtAntiProtons","y-Pt Antiprotons",fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
100 fHistYPtAntiProtons->SetStats(kTRUE);
101 fHistYPtAntiProtons->GetYaxis()->SetTitle("P_{T} [GeV]");
102 fHistYPtAntiProtons->GetXaxis()->SetTitle("y");
103 fHistYPtAntiProtons->GetXaxis()->SetTitleColor(1);
104}
105
106//____________________________________________________________________//
107AliProtonAnalysis::~AliProtonAnalysis() {
108 //Default destructor
39f2a708 109 if(fHistEvents) delete fHistEvents;
110 if(fHistYPtProtons) delete fHistYPtProtons;
111 if(fHistYPtAntiProtons) delete fHistYPtAntiProtons;
112 if(fCorrectionList2D) delete fCorrectionList2D;
113 if(fEfficiencyList1D) delete fEfficiencyList1D;
114 if(fCorrectionList1D) delete fCorrectionList1D;
734d2c12 115}
116
117//____________________________________________________________________//
118void AliProtonAnalysis::InitHistograms(Int_t nbinsY, Float_t fLowY, Float_t fHighY, Int_t nbinsPt, Float_t fLowPt, Float_t fHighPt) {
119 fNBinsY = nbinsY;
120 fMinY = fLowY;
121 fMaxY = fHighY;
122 fNBinsPt = nbinsPt;
123 fMinPt = fLowPt;
124 fMaxPt = fHighPt;
125
3f6d0c08 126 fHistEvents = new TH1I("fHistEvents","Anallyzed events",1,0,1);
127
734d2c12 128 fHistYPtProtons = new TH2F("fHistYPtProtons","y-Pt Protons",fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
129 fHistYPtProtons->SetStats(kTRUE);
130 fHistYPtProtons->GetYaxis()->SetTitle("P_{T} [GeV]");
131 fHistYPtProtons->GetXaxis()->SetTitle("y");
132 fHistYPtProtons->GetXaxis()->SetTitleColor(1);
133
134 fHistYPtAntiProtons = new TH2F("fHistYPtAntiProtons","y-Pt Antiprotons",fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
135 fHistYPtAntiProtons->SetStats(kTRUE);
136 fHistYPtAntiProtons->GetYaxis()->SetTitle("P_{T} [GeV]");
137 fHistYPtAntiProtons->GetXaxis()->SetTitle("y");
138 fHistYPtAntiProtons->GetXaxis()->SetTitleColor(1);
139}
140
141//____________________________________________________________________//
2b748670 142Bool_t AliProtonAnalysis::ReadFromFile(const char* filename) {
143 Bool_t status = kTRUE;
144
734d2c12 145 TFile *file = TFile::Open(filename);
2b748670 146 if(!file) {
147 cout<<"Could not find the input file "<<filename<<endl;
148 status = kFALSE;
149 }
150
e4358d7f 151 TList *list = (TList *)file->Get("outputList1");
2b748670 152 if(list) {
153 cout<<"Retrieving objects from the list "<<list->GetName()<<"..."<<endl;
154 fHistYPtProtons = (TH2F *)list->At(0);
155 fHistYPtAntiProtons = (TH2F *)list->At(1);
3f6d0c08 156 fHistEvents = (TH1I *)list->At(2);
2b748670 157 }
158 else if(!list) {
159 cout<<"Retrieving objects from the file... "<<endl;
160 fHistYPtProtons = (TH2F *)file->Get("fHistYPtProtons");
161 fHistYPtAntiProtons = (TH2F *)file->Get("fHistYPtAntiProtons");
3f6d0c08 162 fHistEvents = (TH1I *)file->Get("fHistEvents");
2b748670 163 }
3f6d0c08 164 if((!fHistYPtProtons)||(!fHistYPtAntiProtons)||(!fHistEvents)) {
2b748670 165 cout<<"Input containers were not found!!!"<<endl;
166 status = kFALSE;
167 }
168 else {
169 fHistYPtProtons->Sumw2();
170 fHistYPtAntiProtons->Sumw2();
171 }
172
173 return status;
734d2c12 174}
175
176//____________________________________________________________________//
177TH1D *AliProtonAnalysis::GetProtonYHistogram() {
3f6d0c08 178 Int_t nAnalyzedEvents = GetNumberOfAnalyzedEvents();
734d2c12 179 TH1D *fYProtons = (TH1D *)fHistYPtProtons->ProjectionX("fYProtons",0,fHistYPtProtons->GetYaxis()->GetNbins(),"e");
180 fYProtons->SetStats(kFALSE);
3f6d0c08 181 fYProtons->GetYaxis()->SetTitle("(1/N_{events})(dN/dy)");
734d2c12 182 fYProtons->SetTitle("dN/dy protons");
183 fYProtons->SetMarkerStyle(kFullCircle);
184 fYProtons->SetMarkerColor(4);
3f6d0c08 185 if(nAnalyzedEvents > 0)
186 fYProtons->Scale(1./nAnalyzedEvents);
734d2c12 187
188 return fYProtons;
189}
190
191//____________________________________________________________________//
192TH1D *AliProtonAnalysis::GetAntiProtonYHistogram() {
3f6d0c08 193 Int_t nAnalyzedEvents = GetNumberOfAnalyzedEvents();
734d2c12 194 TH1D *fYAntiProtons = (TH1D *)fHistYPtAntiProtons->ProjectionX("fYAntiProtons",0,fHistYPtAntiProtons->GetYaxis()->GetNbins(),"e");
195 fYAntiProtons->SetStats(kFALSE);
3f6d0c08 196 fYAntiProtons->GetYaxis()->SetTitle("(1/N_{events})(dN/dy)");
734d2c12 197 fYAntiProtons->SetTitle("dN/dy antiprotons");
198 fYAntiProtons->SetMarkerStyle(kFullCircle);
199 fYAntiProtons->SetMarkerColor(4);
3f6d0c08 200 if(nAnalyzedEvents > 0)
201 fYAntiProtons->Scale(1./nAnalyzedEvents);
734d2c12 202
203 return fYAntiProtons;
204}
205
206//____________________________________________________________________//
207TH1D *AliProtonAnalysis::GetProtonPtHistogram() {
3f6d0c08 208 Int_t nAnalyzedEvents = GetNumberOfAnalyzedEvents();
734d2c12 209 TH1D *fPtProtons = (TH1D *)fHistYPtProtons->ProjectionY("fPtProtons",0,fHistYPtProtons->GetXaxis()->GetNbins(),"e");
210 fPtProtons->SetStats(kFALSE);
3f6d0c08 211 fPtProtons->GetYaxis()->SetTitle("(1/N_{events})(dN/dP_{T})");
734d2c12 212 fPtProtons->SetTitle("dN/dPt protons");
213 fPtProtons->SetMarkerStyle(kFullCircle);
214 fPtProtons->SetMarkerColor(4);
3f6d0c08 215 if(nAnalyzedEvents > 0)
216 fPtProtons->Scale(1./nAnalyzedEvents);
734d2c12 217
218 return fPtProtons;
219}
220
221//____________________________________________________________________//
222TH1D *AliProtonAnalysis::GetAntiProtonPtHistogram() {
3f6d0c08 223 Int_t nAnalyzedEvents = GetNumberOfAnalyzedEvents();
734d2c12 224 TH1D *fPtAntiProtons = (TH1D *)fHistYPtAntiProtons->ProjectionY("fPtAntiProtons",0,fHistYPtProtons->GetXaxis()->GetNbins(),"e");
225 fPtAntiProtons->SetStats(kFALSE);
3f6d0c08 226 fPtAntiProtons->GetYaxis()->SetTitle("(1/N_{events})(dN/dP_{T})");
734d2c12 227 fPtAntiProtons->SetTitle("dN/dPt antiprotons");
228 fPtAntiProtons->SetMarkerStyle(kFullCircle);
229 fPtAntiProtons->SetMarkerColor(4);
3f6d0c08 230 if(nAnalyzedEvents > 0)
231 fPtAntiProtons->Scale(1./nAnalyzedEvents);
734d2c12 232
233 return fPtAntiProtons;
234}
235
236//____________________________________________________________________//
237TH1D *AliProtonAnalysis::GetYRatioHistogram() {
238 TH1D *fYProtons = GetProtonYHistogram();
239 TH1D *fYAntiProtons = GetAntiProtonYHistogram();
240
241 TH1D *hRatioY = new TH1D("hRatioY","",fYProtons->GetNbinsX(),fYProtons->GetXaxis()->GetXmin(),fYProtons->GetXaxis()->GetXmax());
242 hRatioY->Divide(fYAntiProtons,fYProtons,1.0,1.0);
243 hRatioY->SetMarkerStyle(kFullCircle);
244 hRatioY->SetMarkerColor(4);
245 hRatioY->GetYaxis()->SetTitle("#bar{p}/p");
246 hRatioY->GetYaxis()->SetTitleOffset(1.4);
247 hRatioY->GetXaxis()->SetTitle("y");
248 hRatioY->GetXaxis()->SetTitleColor(1);
249 hRatioY->SetStats(kFALSE);
250
251 return hRatioY;
252}
253
254//____________________________________________________________________//
255TH1D *AliProtonAnalysis::GetPtRatioHistogram() {
256 TH1D *fPtProtons = GetProtonPtHistogram();
257 TH1D *fPtAntiProtons = GetAntiProtonPtHistogram();
258
259 TH1D *hRatioPt = new TH1D("hRatioPt","",fPtProtons->GetNbinsX(),fPtProtons->GetXaxis()->GetXmin(),fPtProtons->GetXaxis()->GetXmax());
260 hRatioPt->Divide(fPtAntiProtons,fPtProtons,1.0,1.0);
261 hRatioPt->SetMarkerStyle(kFullCircle);
262 hRatioPt->SetMarkerColor(4);
263 hRatioPt->GetYaxis()->SetTitle("#bar{p}/p");
264 hRatioPt->GetYaxis()->SetTitleOffset(1.4);
265 hRatioPt->GetXaxis()->SetTitle("P_{T} [GeV/c]");
266 hRatioPt->GetXaxis()->SetTitleColor(1);
267 hRatioPt->SetStats(kFALSE);
268
269 return hRatioPt;
270}
271
272//____________________________________________________________________//
273TH1D *AliProtonAnalysis::GetYAsymmetryHistogram() {
274 TH1D *fYProtons = GetProtonYHistogram();
275 TH1D *fYAntiProtons = GetAntiProtonYHistogram();
276
277 TH1D *hsum = new TH1D("hsumY","",fYProtons->GetNbinsX(),fYProtons->GetXaxis()->GetXmin(),fYProtons->GetXaxis()->GetXmax());
278 hsum->Add(fYProtons,fYAntiProtons,1.0,1.0);
279
280 TH1D *hdiff = new TH1D("hdiffY","",fYProtons->GetNbinsX(),fYProtons->GetXaxis()->GetXmin(),fYProtons->GetXaxis()->GetXmax());
281 hdiff->Add(fYProtons,fYAntiProtons,1.0,-1.0);
282
283 TH1D *hAsymmetryY = new TH1D("hAsymmetryY","",fYProtons->GetNbinsX(),fYProtons->GetXaxis()->GetXmin(),fYProtons->GetXaxis()->GetXmax());
284 hAsymmetryY->Divide(hdiff,hsum,2.0,1.);
285 hAsymmetryY->SetMarkerStyle(kFullCircle);
286 hAsymmetryY->SetMarkerColor(4);
287 hAsymmetryY->GetYaxis()->SetTitle("A_{p}");
288 hAsymmetryY->GetYaxis()->SetTitleOffset(1.4);
289 hAsymmetryY->GetXaxis()->SetTitle("y");
290 hAsymmetryY->GetXaxis()->SetTitleColor(1);
291 hAsymmetryY->SetStats(kFALSE);
292
293 return hAsymmetryY;
294}
295
296//____________________________________________________________________//
297TH1D *AliProtonAnalysis::GetPtAsymmetryHistogram() {
298 TH1D *fPtProtons = GetProtonPtHistogram();
299 TH1D *fPtAntiProtons = GetAntiProtonPtHistogram();
300
301 TH1D *hsum = new TH1D("hsumPt","",fPtProtons->GetNbinsX(),fPtProtons->GetXaxis()->GetXmin(),fPtProtons->GetXaxis()->GetXmax());
302 hsum->Add(fPtProtons,fPtAntiProtons,1.0,1.0);
303
304 TH1D *hdiff = new TH1D("hdiffPt","",fPtProtons->GetNbinsX(),fPtProtons->GetXaxis()->GetXmin(),fPtProtons->GetXaxis()->GetXmax());
305 hdiff->Add(fPtProtons,fPtAntiProtons,1.0,-1.0);
306
307 TH1D *hAsymmetryPt = new TH1D("hAsymmetryPt","",fPtProtons->GetNbinsX(),fPtProtons->GetXaxis()->GetXmin(),fPtProtons->GetXaxis()->GetXmax());
308 hAsymmetryPt->Divide(hdiff,hsum,2.0,1.);
309 hAsymmetryPt->SetMarkerStyle(kFullCircle);
310 hAsymmetryPt->SetMarkerColor(4);
311 hAsymmetryPt->GetYaxis()->SetTitle("A_{p}");
312 hAsymmetryPt->GetYaxis()->SetTitleOffset(1.4);
313 hAsymmetryPt->GetXaxis()->SetTitle("P_{T} [GeV/c]");
314 hAsymmetryPt->GetXaxis()->SetTitleColor(1);
315 hAsymmetryPt->SetStats(kFALSE);
316
317 return hAsymmetryPt;
318}
319
aafecd8b 320//____________________________________________________________________//
321Double_t AliProtonAnalysis::GetParticleFraction(Int_t i, Double_t p) {
8b8b0b7a 322 Double_t partFrac=0;
aafecd8b 323 if(fFunctionProbabilityFlag) {
324 if(i == 0) partFrac = fElectronFunction->Eval(p);
325 if(i == 1) partFrac = fMuonFunction->Eval(p);
326 if(i == 2) partFrac = fPionFunction->Eval(p);
327 if(i == 3) partFrac = fKaonFunction->Eval(p);
328 if(i == 4) partFrac = fProtonFunction->Eval(p);
329 }
330 else partFrac = fPartFrac[i];
331
332 return partFrac;
333}
334
734d2c12 335//____________________________________________________________________//
336void AliProtonAnalysis::Analyze(AliESDEvent* fESD) {
e4358d7f 337 //Main analysis part - ESD
3f6d0c08 338 fHistEvents->Fill(0); //number of analyzed events
2b748670 339 Double_t Pt = 0.0, P = 0.0;
734d2c12 340 Int_t nGoodTracks = fESD->GetNumberOfTracks();
341 for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) {
342 AliESDtrack* track = fESD->GetTrack(iTracks);
2b748670 343 Double_t probability[5];
344
345 if(IsAccepted(track)) {
346 if(fUseTPCOnly) {
347 AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
348 if(!tpcTrack) continue;
349 Pt = tpcTrack->Pt();
350 P = tpcTrack->P();
351
352 //pid
353 track->GetTPCpid(probability);
354 Double_t rcc = 0.0;
355 for(Int_t i = 0; i < AliPID::kSPECIES; i++)
356 rcc += probability[i]*GetParticleFraction(i,P);
357 if(rcc == 0.0) continue;
358 Double_t w[5];
359 for(Int_t i = 0; i < AliPID::kSPECIES; i++)
360 w[i] = probability[i]*GetParticleFraction(i,P)/rcc;
361 Long64_t fParticleType = TMath::LocMax(AliPID::kSPECIES,w);
362 if(fParticleType == 4) {
363 if(tpcTrack->Charge() > 0)
364 fHistYPtProtons->Fill(Rapidity(tpcTrack->Px(),tpcTrack->Py(),tpcTrack->Pz()),Pt);
365 else if(tpcTrack->Charge() < 0)
366 fHistYPtAntiProtons->Fill(Rapidity(tpcTrack->Px(),tpcTrack->Py(),tpcTrack->Pz()),Pt);
367 }//proton check
368 }//TPC only tracks
369 else if(!fUseTPCOnly) {
370 Pt = track->Pt();
371 P = track->P();
734d2c12 372
2b748670 373 //pid
374 track->GetESDpid(probability);
375 Double_t rcc = 0.0;
376 for(Int_t i = 0; i < AliPID::kSPECIES; i++)
377 rcc += probability[i]*GetParticleFraction(i,P);
378 if(rcc == 0.0) continue;
379 Double_t w[5];
380 for(Int_t i = 0; i < AliPID::kSPECIES; i++)
381 w[i] = probability[i]*GetParticleFraction(i,P)/rcc;
382 Long64_t fParticleType = TMath::LocMax(AliPID::kSPECIES,w);
383 if(fParticleType == 4) {
384 //cout<<"(Anti)protons found..."<<endl;
385 if(track->Charge() > 0)
386 fHistYPtProtons->Fill(Rapidity(track->Px(),track->Py(),track->Pz()),Pt);
387 else if(track->Charge() < 0)
388 fHistYPtAntiProtons->Fill(Rapidity(track->Px(),track->Py(),track->Pz()),Pt);
389 }//proton check
390 }//combined tracking
734d2c12 391 }//cuts
392 }//track loop
393}
394
ee4ca40d 395//____________________________________________________________________//
396void AliProtonAnalysis::Analyze(AliAODEvent* fAOD) {
e4358d7f 397 //Main analysis part - AOD
3f6d0c08 398 fHistEvents->Fill(0); //number of analyzed events
ee4ca40d 399 Int_t nGoodTracks = fAOD->GetNumberOfTracks();
400 for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) {
401 AliAODTrack* track = fAOD->GetTrack(iTracks);
402 Double_t Pt = track->Pt();
403 Double_t P = track->P();
404
405 //pid
406 Double_t probability[10];
407 track->GetPID(probability);
408 Double_t rcc = 0.0;
409 for(Int_t i = 0; i < AliPID::kSPECIESN; i++) rcc += probability[i]*GetParticleFraction(i,P);
410 if(rcc == 0.0) continue;
738619fd 411 Double_t w[10];
ee4ca40d 412 for(Int_t i = 0; i < AliPID::kSPECIESN; i++) w[i] = probability[i]*GetParticleFraction(i,P)/rcc;
413 Long64_t fParticleType = TMath::LocMax(AliPID::kSPECIESN,w);
414 if(fParticleType == 4) {
415 if(track->Charge() > 0) fHistYPtProtons->Fill(track->Y(fParticleType),Pt);
416 else if(track->Charge() < 0) fHistYPtAntiProtons->Fill(track->Y(fParticleType),Pt);
417 }//proton check
418 }//track loop
419}
420
e4358d7f 421//____________________________________________________________________//
422void AliProtonAnalysis::Analyze(AliStack* stack) {
423 //Main analysis part - MC
3f6d0c08 424 fHistEvents->Fill(0); //number of analyzed events
e4358d7f 425 for(Int_t i = 0; i < stack->GetNprimary(); i++) {
426 TParticle *particle = stack->Particle(i);
427 if(particle->Pt() < 0.1) continue;
428 if(TMath::Abs(particle->Eta()) > 1.0) continue;
429 Int_t pdgcode = particle->GetPdgCode();
430 if(pdgcode == 2212) fHistYPtProtons->Fill(Rapidity(particle->Px(),
431 particle->Py(),
432 particle->Pz()),
433 particle->Pt());
434 if(pdgcode == -2212) fHistYPtAntiProtons->Fill(Rapidity(particle->Px(),
435 particle->Py(),
436 particle->Pz()),
437 particle->Pt());
438 }//particle loop
439}
440
734d2c12 441//____________________________________________________________________//
442Bool_t AliProtonAnalysis::IsAccepted(AliESDtrack* track) {
443 // Checks if the track is excluded from the cuts
2b748670 444 Double_t Pt = 0.0, Px = 0.0, Py = 0.0, Pz = 0.0;
445 if(fUseTPCOnly) {
446 AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
447 if(!tpcTrack) {
448 Pt = 0.0; Px = 0.0; Py = 0.0; Pz = 0.0;
449 }
450 else {
451 Pt = tpcTrack->Pt();
452 Px = tpcTrack->Px();
453 Py = tpcTrack->Py();
454 Pz = tpcTrack->Pz();
455 }
456 }
457 else{
458 Pt = track->Pt();
459 Px = track->Px();
460 Py = track->Py();
461 Pz = track->Pz();
462 }
463
734d2c12 464 Int_t fIdxInt[200];
465 Int_t nClustersITS = track->GetITSclusters(fIdxInt);
466 Int_t nClustersTPC = track->GetTPCclusters(fIdxInt);
467
468 Float_t chi2PerClusterITS = -1;
469 Float_t chi2PerClusterTPC = -1;
470 if (nClustersTPC!=0)
471 chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC);
472
473 Double_t extCov[15];
474 track->GetExternalCovariance(extCov);
475
a5375b97 476 if(fMinITSClustersFlag)
734d2c12 477 if(nClustersITS < fMinITSClusters) return kFALSE;
478 if(fMinTPCClustersFlag)
479 if(nClustersTPC < fMinTPCClusters) return kFALSE;
480 if(fMaxChi2PerTPCClusterFlag)
481 if(chi2PerClusterTPC > fMaxChi2PerTPCCluster) return kFALSE;
482 if(fMaxChi2PerITSClusterFlag)
483 if(chi2PerClusterITS > fMaxChi2PerITSCluster) return kFALSE;
484 if(fMaxCov11Flag)
485 if(extCov[0] > fMaxCov11) return kFALSE;
486 if(fMaxCov22Flag)
487 if(extCov[2] > fMaxCov22) return kFALSE;
488 if(fMaxCov33Flag)
489 if(extCov[5] > fMaxCov33) return kFALSE;
490 if(fMaxCov44Flag)
491 if(extCov[9] > fMaxCov44) return kFALSE;
492 if(fMaxCov55Flag)
493 if(extCov[14] > fMaxCov55) return kFALSE;
494 if(fMaxSigmaToVertexFlag)
495 if(GetSigmaToVertex(track) > fMaxSigmaToVertex) return kFALSE;
496 if(fITSRefitFlag)
497 if ((track->GetStatus() & AliESDtrack::kITSrefit) == 0) return kFALSE;
498 if(fTPCRefitFlag)
499 if ((track->GetStatus() & AliESDtrack::kTPCrefit) == 0) return kFALSE;
500
501 if((Pt < fMinPt) || (Pt > fMaxPt)) return kFALSE;
2b748670 502 if((Rapidity(Px,Py,Pz) < fMinY) || (Rapidity(Px,Py,Pz) > fMaxY)) return kFALSE;
734d2c12 503
504 return kTRUE;
505}
506
507//____________________________________________________________________//
508Float_t AliProtonAnalysis::GetSigmaToVertex(AliESDtrack* esdTrack) {
509 // Calculates the number of sigma to the vertex.
510
511 Float_t b[2];
512 Float_t bRes[2];
513 Float_t bCov[3];
f2eeda10 514 if(fUseTPCOnly)
515 esdTrack->GetImpactParametersTPC(b,bCov);
516 else
517 esdTrack->GetImpactParameters(b,bCov);
518
734d2c12 519 if (bCov[0]<=0 || bCov[2]<=0) {
520 //AliDebug(1, "Estimated b resolution lower or equal zero!");
521 bCov[0]=0; bCov[2]=0;
522 }
523 bRes[0] = TMath::Sqrt(bCov[0]);
524 bRes[1] = TMath::Sqrt(bCov[2]);
525
526 if (bRes[0] == 0 || bRes[1] ==0) return -1;
527
528 Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
529
530 if (TMath::Exp(-d * d / 2) < 1e-10) return 1000;
531
532 d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
533
534 return d;
535}
536
3f6d0c08 537//____________________________________________________________________//
2b748670 538Double_t AliProtonAnalysis::Rapidity(Double_t Px, Double_t Py, Double_t Pz) {
3f6d0c08 539 //returns the rapidity of the proton - to be removed
734d2c12 540 Double_t fMass = 9.38270000000000048e-01;
541
2b748670 542 Double_t P = TMath::Sqrt(TMath::Power(Px,2) +
543 TMath::Power(Py,2) +
544 TMath::Power(Pz,2));
734d2c12 545 Double_t energy = TMath::Sqrt(P*P + fMass*fMass);
546 Double_t y = -999;
2b748670 547 if(energy != Pz)
548 y = 0.5*TMath::Log((energy + Pz)/(energy - Pz));
734d2c12 549
550 return y;
551}
3f6d0c08 552
553//____________________________________________________________________//
554Bool_t AliProtonAnalysis::PrintMean(TH1 *hist, Double_t edge) {
555 //calculates the mean value of the ratio/asymmetry within \pm edge
556 Double_t sum = 0.0;
557 Int_t nentries = 0;
558 //calculate the mean
559 for(Int_t i = 0; i < hist->GetXaxis()->GetNbins(); i++) {
560 Double_t x = hist->GetBinCenter(i+1);
561 Double_t y = hist->GetBinContent(i+1);
562 if(TMath::Abs(x) < edge) {
563 sum += y;
564 nentries += 1;
565 }
566 }
567 Double_t mean = 0.0;
568 if(nentries != 0)
569 mean = sum/nentries;
570
571 //calculate the error
572 for(Int_t i = 0; i < hist->GetXaxis()->GetNbins(); i++) {
573 Double_t x = hist->GetBinCenter(i+1);
574 Double_t y = hist->GetBinContent(i+1);
575 if(TMath::Abs(x) < edge) {
576 sum += TMath::Power((mean - y),2);
577 nentries += 1;
578 }
579 }
580
581 Double_t error = 0.0;
582 if(nentries != 0)
583 error = TMath::Sqrt(sum)/nentries;
584
585 cout<<"========================================="<<endl;
586 cout<<"Input distribution: "<<hist->GetName()<<endl;
587 cout<<"Interval used: -"<<edge<<" -> "<<edge<<endl;
588 cout<<"Mean value :"<<mean<<endl;
589 cout<<"Error: "<<error<<endl;
590 cout<<"========================================="<<endl;
591
592 return 0;
593}
594
595//____________________________________________________________________//
596Bool_t AliProtonAnalysis::PrintYields(TH1 *hist, Double_t edge) {
597 //calculates the (anti)proton yields within the \pm edge
598 Double_t sum = 0.0, sumerror = 0.0;
599 Double_t error = 0.0;
600 for(Int_t i = 0; i < hist->GetXaxis()->GetNbins(); i++) {
601 Double_t x = hist->GetBinCenter(i+1);
602 Double_t y = hist->GetBinContent(i+1);
603 if(TMath::Abs(x) < edge) {
604 sum += y;
605 sumerror += TMath::Power(hist->GetBinError(i+1),2);
606 }
607 }
608
609 error = TMath::Sqrt(sumerror);
610
611 cout<<"========================================="<<endl;
612 cout<<"Input distribution: "<<hist->GetName()<<endl;
613 cout<<"Interval used: -"<<edge<<" -> "<<edge<<endl;
614 cout<<"Yields :"<<sum<<endl;
615 cout<<"Error: "<<error<<endl;
616 cout<<"========================================="<<endl;
617
618 return 0;
619}
620
39f2a708 621//____________________________________________________________________//
622Bool_t AliProtonAnalysis::ReadCorrectionContainer(const char* filename) {
623 // Reads the outout of the correction framework task
624 // Creates the correction maps
625 // Puts the results in the different TList objects
626 Bool_t status = kTRUE;
627
628 TFile *file = TFile::Open(filename);
629 if(!file) {
630 cout<<"Could not find the input CORRFW file "<<filename<<endl;
631 status = kFALSE;
632 }
633
634 AliCFContainer *corrfwContainer = (AliCFContainer*) (file->Get("container"));
635 if(!corrfwContainer) {
636 cout<<"CORRFW container not found!"<<endl;
637 status = kFALSE;
638 }
639
640 Int_t nSteps = corrfwContainer->GetNStep();
641 TH2D *gYPt[4];
642 //currently the GRID is formed by the y-pT parameters
643 //Add Vz as a next step
644 Int_t iRap = 0, iPt = 1;
645 for(Int_t iStep = 0; iStep < nSteps; iStep++) {
646 gYPt[iStep] = corrfwContainer->ShowProjection(iRap,iPt,iStep);
647 //fCorrectionList2D->Add(gYPt[iStep]);
648 }
649
650 //construct the efficiency grid from the data container
651 TString gTitle = 0;
652 AliCFEffGrid *efficiency[3]; //The efficiency array has nStep-1 entries!!!
653 TH1D *gEfficiency[2][3]; //efficiency as a function of pT and of y (raws - [2])
654 TH1D *gCorrection[2][3]; //efficiency as a function of pT and of y (raws - [2])
655
656 //Get the 2D efficiency maps
657 for(Int_t iStep = 1; iStep < nSteps; iStep++) {
658 gTitle = "Efficiency_Step0_Step"; gTitle += iStep;
659 efficiency[iStep] = new AliCFEffGrid(gTitle.Data(),
660 gTitle.Data(),*corrfwContainer);
661 efficiency[iStep]->CalculateEfficiency(iStep,0); //eff= step[i]/step0
662 fCorrectionList2D->Add(efficiency[iStep]);
663 }
664 //Get the projection of the efficiency maps
665 for(Int_t iParameter = 0; iParameter < 2; iParameter++) {
666 for(Int_t iStep = 1; iStep < nSteps; iStep++) {
667 gEfficiency[iParameter][iStep-1] = efficiency[iStep]->Project(iParameter);
668 fEfficiencyList1D->Add(gEfficiency[iParameter][iStep-1]);
669 gTitle = "Correction_Parameter"; gTitle += iParameter+1;
670 gTitle += "_Step0_Step"; gTitle += iStep;
671 gCorrection[iParameter][iStep-1] = new TH1D(gTitle.Data(),
672 gTitle.Data(),
673 gEfficiency[iParameter][iStep-1]->GetNbinsX(),
674 gEfficiency[iParameter][iStep-1]->GetXaxis()->GetXmin(),
675 gEfficiency[iParameter][iStep-1]->GetXaxis()->GetXmax());
676 //initialisation of the correction
677 for(Int_t iBin = 1; iBin <= gEfficiency[iParameter][iStep-1]->GetNbinsX(); iBin++)
678 gCorrection[iParameter][iStep-1]->SetBinContent(iBin,1.0);
679 }//step loop
680 }//parameter loop
681 //Calculate the 1D correction parameters as a function of y and pT
682 for(Int_t iParameter = 0; iParameter < 2; iParameter++) {
683 for(Int_t iStep = 1; iStep < nSteps; iStep++) {
684 gCorrection[iParameter][iStep-1]->Divide(gEfficiency[iParameter][iStep-1]);
685 fCorrectionList1D->Add(gCorrection[iParameter][iStep-1]);
686 }
687 }
688}
689
690
691
692
693
694
695
696
3f6d0c08 697
698