]> 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
In AliMUONTrackerIO:
[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
dd3fa486 19// Origin: Panos Christakoglou | Panos.Christakoglou@cern.ch
734d2c12 20//-----------------------------------------------------------------
21#include <Riostream.h>
22#include <TFile.h>
23#include <TSystem.h>
aafecd8b 24#include <TF1.h>
251e4034 25#include <TH2D.h>
734d2c12 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>
251e4034 40#include <AliCFDataGrid.h>
6667f3a7 41#include <AliESDVertex.h>
e4358d7f 42
734d2c12 43ClassImp(AliProtonAnalysis)
44
45//____________________________________________________________________//
46AliProtonAnalysis::AliProtonAnalysis() :
dd3fa486 47 TObject(), fAnalysisEtaMode(kFALSE),
734d2c12 48 fNBinsY(0), fMinY(0), fMaxY(0),
49 fNBinsPt(0), fMinPt(0), fMaxPt(0),
50 fMinTPCClusters(0), fMinITSClusters(0),
51 fMaxChi2PerTPCCluster(0), fMaxChi2PerITSCluster(0),
52 fMaxCov11(0), fMaxCov22(0), fMaxCov33(0), fMaxCov44(0), fMaxCov55(0),
f24a70b0 53 fMaxSigmaToVertex(0), fMaxSigmaToVertexTPC(0),
96f84c25 54 fMaxDCAXY(0), fMaxDCAXYTPC(0),
55 fMaxDCAZ(0), fMaxDCAZTPC(0),
56 fMaxConstrainChi2(0),
734d2c12 57 fMinTPCClustersFlag(kFALSE), fMinITSClustersFlag(kFALSE),
58 fMaxChi2PerTPCClusterFlag(kFALSE), fMaxChi2PerITSClusterFlag(kFALSE),
24421eb6 59 fMaxCov11Flag(kFALSE), fMaxCov22Flag(kFALSE),
60 fMaxCov33Flag(kFALSE), fMaxCov44Flag(kFALSE), fMaxCov55Flag(kFALSE),
f24a70b0 61 fMaxSigmaToVertexFlag(kFALSE), fMaxSigmaToVertexTPCFlag(kFALSE),
96f84c25 62 fMaxDCAXYFlag(kFALSE), fMaxDCAXYTPCFlag(kFALSE),
63 fMaxDCAZFlag(kFALSE), fMaxDCAZTPCFlag(kFALSE),
64 fMaxConstrainChi2Flag(kFALSE),
734d2c12 65 fITSRefitFlag(kFALSE), fTPCRefitFlag(kFALSE),
f24a70b0 66 fESDpidFlag(kFALSE), fTPCpidFlag(kFALSE),
0008a5a6 67 fPointOnITSLayer1Flag(0), fPointOnITSLayer2Flag(0),
68 fPointOnITSLayer3Flag(0), fPointOnITSLayer4Flag(0),
69 fPointOnITSLayer5Flag(0), fPointOnITSLayer6Flag(0),
aafecd8b 70 fFunctionProbabilityFlag(kFALSE),
71 fElectronFunction(0), fMuonFunction(0),
72 fPionFunction(0), fKaonFunction(0), fProtonFunction(0),
55f9a666 73 fUseTPCOnly(kFALSE), fUseHybridTPC(kFALSE),
251e4034 74 fProtonContainer(0), fAntiProtonContainer(0),
75 fHistEvents(0), fHistYPtProtons(0), fHistYPtAntiProtons(0),
76 fEffGridListProtons(0), fCorrectionListProtons2D(0),
77 fEfficiencyListProtons1D(0), fCorrectionListProtons1D(0),
78 fEffGridListAntiProtons(0), fCorrectionListAntiProtons2D(0),
79 fEfficiencyListAntiProtons1D(0), fCorrectionListAntiProtons1D(0),
80 fCorrectProtons(0), fCorrectAntiProtons(0) {
734d2c12 81 //Default constructor
82 for(Int_t i = 0; i < 5; i++) fPartFrac[i] = 0.0;
83}
84
85//____________________________________________________________________//
86AliProtonAnalysis::AliProtonAnalysis(Int_t nbinsY, Float_t fLowY, Float_t fHighY,Int_t nbinsPt, Float_t fLowPt, Float_t fHighPt) :
dd3fa486 87 TObject(), fAnalysisEtaMode(kFALSE),
734d2c12 88 fNBinsY(nbinsY), fMinY(fLowY), fMaxY(fHighY),
89 fNBinsPt(nbinsPt), fMinPt(fLowPt), fMaxPt(fHighPt),
aafecd8b 90 fMinTPCClusters(0), fMinITSClusters(0),
91 fMaxChi2PerTPCCluster(0), fMaxChi2PerITSCluster(0),
92 fMaxCov11(0), fMaxCov22(0), fMaxCov33(0), fMaxCov44(0), fMaxCov55(0),
f24a70b0 93 fMaxSigmaToVertex(0), fMaxSigmaToVertexTPC(0),
96f84c25 94 fMaxDCAXY(0), fMaxDCAXYTPC(0),
95 fMaxDCAZ(0), fMaxDCAZTPC(0),
96 fMaxConstrainChi2(0),
aafecd8b 97 fMinTPCClustersFlag(kFALSE), fMinITSClustersFlag(kFALSE),
98 fMaxChi2PerTPCClusterFlag(kFALSE), fMaxChi2PerITSClusterFlag(kFALSE),
24421eb6 99 fMaxCov11Flag(kFALSE), fMaxCov22Flag(kFALSE),
100 fMaxCov33Flag(kFALSE), fMaxCov44Flag(kFALSE), fMaxCov55Flag(kFALSE),
f24a70b0 101 fMaxSigmaToVertexFlag(kFALSE), fMaxSigmaToVertexTPCFlag(kFALSE),
96f84c25 102 fMaxDCAXYFlag(kFALSE), fMaxDCAXYTPCFlag(kFALSE),
103 fMaxDCAZFlag(kFALSE), fMaxDCAZTPCFlag(kFALSE),
104 fMaxConstrainChi2Flag(kFALSE),
aafecd8b 105 fITSRefitFlag(kFALSE), fTPCRefitFlag(kFALSE),
f24a70b0 106 fESDpidFlag(kFALSE), fTPCpidFlag(kFALSE),
0008a5a6 107 fPointOnITSLayer1Flag(0), fPointOnITSLayer2Flag(0),
108 fPointOnITSLayer3Flag(0), fPointOnITSLayer4Flag(0),
109 fPointOnITSLayer5Flag(0), fPointOnITSLayer6Flag(0),
aafecd8b 110 fFunctionProbabilityFlag(kFALSE),
111 fElectronFunction(0), fMuonFunction(0),
112 fPionFunction(0), fKaonFunction(0), fProtonFunction(0),
55f9a666 113 fUseTPCOnly(kFALSE), fUseHybridTPC(kFALSE),
251e4034 114 fProtonContainer(0), fAntiProtonContainer(0),
115 fHistEvents(0), fHistYPtProtons(0), fHistYPtAntiProtons(0),
116 fEffGridListProtons(0), fCorrectionListProtons2D(0),
117 fEfficiencyListProtons1D(0), fCorrectionListProtons1D(0),
118 fEffGridListAntiProtons(0), fCorrectionListAntiProtons2D(0),
119 fEfficiencyListAntiProtons1D(0), fCorrectionListAntiProtons1D(0),
120 fCorrectProtons(0), fCorrectAntiProtons(0){
734d2c12 121 //Default constructor
3f6d0c08 122 fHistEvents = new TH1I("fHistEvents","Analyzed events",1,0,1);
123
3e6c06f4 124 fHistYPtProtons = new TH2D("fHistYPtProtons","y-Pt Protons",
125 fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
734d2c12 126 fHistYPtProtons->SetStats(kTRUE);
127 fHistYPtProtons->GetYaxis()->SetTitle("P_{T} [GeV]");
128 fHistYPtProtons->GetXaxis()->SetTitle("y");
129 fHistYPtProtons->GetXaxis()->SetTitleColor(1);
130
3e6c06f4 131 fHistYPtAntiProtons = new TH2D("fHistYPtAntiProtons","y-Pt Antiprotons",
132 fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
734d2c12 133 fHistYPtAntiProtons->SetStats(kTRUE);
134 fHistYPtAntiProtons->GetYaxis()->SetTitle("P_{T} [GeV]");
135 fHistYPtAntiProtons->GetXaxis()->SetTitle("y");
136 fHistYPtAntiProtons->GetXaxis()->SetTitleColor(1);
251e4034 137
138 //setting up the containers
139 Int_t iBin[2];
140 iBin[0] = nbinsY;
141 iBin[1] = nbinsPt;
142 Double_t *binLimY = new Double_t[nbinsY+1];
143 Double_t *binLimPt = new Double_t[nbinsPt+1];
144 //values for bin lower bounds
145 for(Int_t i = 0; i <= nbinsY; i++)
146 binLimY[i]=(Double_t)fLowY + (fHighY - fLowY) /nbinsY*(Double_t)i;
147 for(Int_t i = 0; i <= nbinsPt; i++)
148 binLimPt[i]=(Double_t)fLowPt + (fHighPt - fLowPt) /nbinsPt*(Double_t)i;
149
150 fProtonContainer = new AliCFContainer("containerProtons",
151 "container for protons",
152 1,2,iBin);
153 fProtonContainer->SetBinLimits(0,binLimY); //rapidity
154 fProtonContainer->SetBinLimits(1,binLimPt); //pT
155 fAntiProtonContainer = new AliCFContainer("containerAntiProtons",
156 "container for antiprotons",
157 1,2,iBin);
158 fAntiProtonContainer->SetBinLimits(0,binLimY); //rapidity
159 fAntiProtonContainer->SetBinLimits(1,binLimPt); //pT
734d2c12 160}
161
162//____________________________________________________________________//
163AliProtonAnalysis::~AliProtonAnalysis() {
164 //Default destructor
39f2a708 165 if(fHistEvents) delete fHistEvents;
166 if(fHistYPtProtons) delete fHistYPtProtons;
167 if(fHistYPtAntiProtons) delete fHistYPtAntiProtons;
251e4034 168 if(fProtonContainer) delete fProtonContainer;
169 if(fAntiProtonContainer) delete fAntiProtonContainer;
170
251e4034 171 if(fEffGridListProtons) delete fEffGridListProtons;
cdb3530f 172 if(fCorrectionListProtons2D) delete fCorrectionListProtons2D;
173 if(fEfficiencyListProtons1D) delete fEfficiencyListProtons1D;
174 if(fCorrectionListProtons1D) delete fCorrectionListProtons1D;
251e4034 175 if(fEffGridListAntiProtons) delete fEffGridListAntiProtons;
cdb3530f 176 if(fCorrectionListAntiProtons2D) delete fCorrectionListAntiProtons2D;
177 if(fEfficiencyListAntiProtons1D) delete fEfficiencyListAntiProtons1D;
178 if(fCorrectionListAntiProtons1D) delete fCorrectionListAntiProtons1D;
251e4034 179 if(fCorrectProtons) delete fCorrectProtons;
180 if(fCorrectAntiProtons) delete fCorrectAntiProtons;
734d2c12 181}
182
183//____________________________________________________________________//
251e4034 184void AliProtonAnalysis::InitAnalysisHistograms(Int_t nbinsY,
185 Float_t fLowY, Float_t fHighY,
186 Int_t nbinsPt,
187 Float_t fLowPt, Float_t fHighPt) {
734d2c12 188 fNBinsY = nbinsY;
189 fMinY = fLowY;
190 fMaxY = fHighY;
191 fNBinsPt = nbinsPt;
192 fMinPt = fLowPt;
193 fMaxPt = fHighPt;
194
3f6d0c08 195 fHistEvents = new TH1I("fHistEvents","Anallyzed events",1,0,1);
196
251e4034 197 fHistYPtProtons = new TH2D("fHistYPtProtons","y-Pt Protons",
198 fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
734d2c12 199 fHistYPtProtons->SetStats(kTRUE);
200 fHistYPtProtons->GetYaxis()->SetTitle("P_{T} [GeV]");
201 fHistYPtProtons->GetXaxis()->SetTitle("y");
202 fHistYPtProtons->GetXaxis()->SetTitleColor(1);
203
251e4034 204 fHistYPtAntiProtons = new TH2D("fHistYPtAntiProtons","y-Pt Antiprotons",
205 fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
734d2c12 206 fHistYPtAntiProtons->SetStats(kTRUE);
207 fHistYPtAntiProtons->GetYaxis()->SetTitle("P_{T} [GeV]");
208 fHistYPtAntiProtons->GetXaxis()->SetTitle("y");
209 fHistYPtAntiProtons->GetXaxis()->SetTitleColor(1);
251e4034 210
211 //setting up the containers
212 Int_t iBin[2];
213 iBin[0] = nbinsY;
214 iBin[1] = nbinsPt;
215 Double_t *binLimY = new Double_t[nbinsY+1];
216 Double_t *binLimPt = new Double_t[nbinsPt+1];
217 //values for bin lower bounds
218 for(Int_t i = 0; i <= nbinsY; i++)
219 binLimY[i]=(Double_t)fLowY + (fHighY - fLowY) /nbinsY*(Double_t)i;
220 for(Int_t i = 0; i <= nbinsPt; i++)
221 binLimPt[i]=(Double_t)fLowPt + (fHighPt - fLowPt) /nbinsPt*(Double_t)i;
222
223 fProtonContainer = new AliCFContainer("containerProtons",
224 "container for protons",
225 1,2,iBin);
226 fProtonContainer->SetBinLimits(0,binLimY); //rapidity
227 fProtonContainer->SetBinLimits(1,binLimPt); //pT
228 fAntiProtonContainer = new AliCFContainer("containerAntiProtons",
229 "container for antiprotons",
230 1,2,iBin);
231 fAntiProtonContainer->SetBinLimits(0,binLimY); //rapidity
232 fAntiProtonContainer->SetBinLimits(1,binLimPt); //pT
734d2c12 233}
234
235//____________________________________________________________________//
2b748670 236Bool_t AliProtonAnalysis::ReadFromFile(const char* filename) {
237 Bool_t status = kTRUE;
238
734d2c12 239 TFile *file = TFile::Open(filename);
2b748670 240 if(!file) {
241 cout<<"Could not find the input file "<<filename<<endl;
242 status = kFALSE;
243 }
244
e4358d7f 245 TList *list = (TList *)file->Get("outputList1");
2b748670 246 if(list) {
247 cout<<"Retrieving objects from the list "<<list->GetName()<<"..."<<endl;
251e4034 248 fHistYPtProtons = (TH2D *)list->At(0);
249 fHistYPtAntiProtons = (TH2D *)list->At(1);
3f6d0c08 250 fHistEvents = (TH1I *)list->At(2);
251e4034 251 fProtonContainer = (AliCFContainer *)list->At(3);
252 fAntiProtonContainer = (AliCFContainer *)list->At(4);
2b748670 253 }
254 else if(!list) {
255 cout<<"Retrieving objects from the file... "<<endl;
251e4034 256 fHistYPtProtons = (TH2D *)file->Get("fHistYPtProtons");
257 fHistYPtAntiProtons = (TH2D *)file->Get("fHistYPtAntiProtons");
3f6d0c08 258 fHistEvents = (TH1I *)file->Get("fHistEvents");
251e4034 259 fProtonContainer = (AliCFContainer *)file->Get("containerProtons");
260 fAntiProtonContainer = (AliCFContainer *)file->Get("containerAntiProtons");
2b748670 261 }
251e4034 262 if((!fHistYPtProtons)||(!fHistYPtAntiProtons)||(!fHistEvents)
263 ||(!fProtonContainer)||(!fAntiProtonContainer)) {
2b748670 264 cout<<"Input containers were not found!!!"<<endl;
265 status = kFALSE;
266 }
267 else {
3e23254a 268 //fHistYPtProtons = fProtonContainer->ShowProjection(0,1,0);
269 //fHistYPtAntiProtons = fAntiProtonContainer->ShowProjection(0,1,0);
270 fHistYPtProtons->Sumw2();
271 fHistYPtAntiProtons->Sumw2();
2b748670 272 }
273
274 return status;
734d2c12 275}
276
277//____________________________________________________________________//
278TH1D *AliProtonAnalysis::GetProtonYHistogram() {
3f6d0c08 279 Int_t nAnalyzedEvents = GetNumberOfAnalyzedEvents();
251e4034 280
3e23254a 281 TH1D *fYProtons = (TH1D *)fHistYPtProtons->ProjectionX("fYProtons",0,fHistYPtProtons->GetYaxis()->GetNbins(),"");
282 //TH1D *fYProtons = fProtonContainer->ShowProjection(0,0); //variable-step
251e4034 283
734d2c12 284 fYProtons->SetStats(kFALSE);
3f6d0c08 285 fYProtons->GetYaxis()->SetTitle("(1/N_{events})(dN/dy)");
734d2c12 286 fYProtons->SetTitle("dN/dy protons");
287 fYProtons->SetMarkerStyle(kFullCircle);
288 fYProtons->SetMarkerColor(4);
3f6d0c08 289 if(nAnalyzedEvents > 0)
251e4034 290 fYProtons->Scale(1./nAnalyzedEvents);
291
734d2c12 292 return fYProtons;
293}
294
295//____________________________________________________________________//
296TH1D *AliProtonAnalysis::GetAntiProtonYHistogram() {
3f6d0c08 297 Int_t nAnalyzedEvents = GetNumberOfAnalyzedEvents();
3e23254a 298
299 TH1D *fYAntiProtons = (TH1D *)fHistYPtAntiProtons->ProjectionX("fYAntiProtons",0,fHistYPtAntiProtons->GetYaxis()->GetNbins(),"");
300 //TH1D *fYAntiProtons = fAntiProtonContainer->ShowProjection(0,0);//variable-step
251e4034 301
734d2c12 302 fYAntiProtons->SetStats(kFALSE);
3f6d0c08 303 fYAntiProtons->GetYaxis()->SetTitle("(1/N_{events})(dN/dy)");
734d2c12 304 fYAntiProtons->SetTitle("dN/dy antiprotons");
305 fYAntiProtons->SetMarkerStyle(kFullCircle);
306 fYAntiProtons->SetMarkerColor(4);
3f6d0c08 307 if(nAnalyzedEvents > 0)
308 fYAntiProtons->Scale(1./nAnalyzedEvents);
734d2c12 309
310 return fYAntiProtons;
311}
312
313//____________________________________________________________________//
314TH1D *AliProtonAnalysis::GetProtonPtHistogram() {
3f6d0c08 315 Int_t nAnalyzedEvents = GetNumberOfAnalyzedEvents();
3e23254a 316
317 TH1D *fPtProtons = (TH1D *)fHistYPtProtons->ProjectionY("fPtProtons",0,fHistYPtProtons->GetXaxis()->GetNbins(),"");
318 //TH1D *fPtProtons = fProtonContainer->ShowProjection(1,0); //variable-step
251e4034 319
734d2c12 320 fPtProtons->SetStats(kFALSE);
3f6d0c08 321 fPtProtons->GetYaxis()->SetTitle("(1/N_{events})(dN/dP_{T})");
734d2c12 322 fPtProtons->SetTitle("dN/dPt protons");
323 fPtProtons->SetMarkerStyle(kFullCircle);
324 fPtProtons->SetMarkerColor(4);
3f6d0c08 325 if(nAnalyzedEvents > 0)
326 fPtProtons->Scale(1./nAnalyzedEvents);
734d2c12 327
328 return fPtProtons;
329}
330
331//____________________________________________________________________//
332TH1D *AliProtonAnalysis::GetAntiProtonPtHistogram() {
3f6d0c08 333 Int_t nAnalyzedEvents = GetNumberOfAnalyzedEvents();
3e23254a 334
335 TH1D *fPtAntiProtons = (TH1D *)fHistYPtAntiProtons->ProjectionY("fPtAntiProtons",0,fHistYPtProtons->GetXaxis()->GetNbins(),"");
336 //TH1D *fPtAntiProtons = fAntiProtonContainer->ShowProjection(1,0); //variable-step
251e4034 337
734d2c12 338 fPtAntiProtons->SetStats(kFALSE);
3f6d0c08 339 fPtAntiProtons->GetYaxis()->SetTitle("(1/N_{events})(dN/dP_{T})");
734d2c12 340 fPtAntiProtons->SetTitle("dN/dPt antiprotons");
341 fPtAntiProtons->SetMarkerStyle(kFullCircle);
342 fPtAntiProtons->SetMarkerColor(4);
3f6d0c08 343 if(nAnalyzedEvents > 0)
344 fPtAntiProtons->Scale(1./nAnalyzedEvents);
734d2c12 345
346 return fPtAntiProtons;
347}
348
251e4034 349//____________________________________________________________________//
350TH1D *AliProtonAnalysis::GetProtonCorrectedYHistogram() {
351 Int_t nAnalyzedEvents = GetNumberOfAnalyzedEvents();
352
353 TH1D *fYProtons = fCorrectProtons->Project(0); //0: rapidity
354
355 fYProtons->SetStats(kFALSE);
356 fYProtons->GetYaxis()->SetTitle("(1/N_{events})(dN/dy)");
357 fYProtons->GetXaxis()->SetTitle("y");
358 fYProtons->SetTitle("dN/dy protons");
359 fYProtons->SetMarkerStyle(kFullCircle);
360 fYProtons->SetMarkerColor(4);
361 if(nAnalyzedEvents > 0)
362 fYProtons->Scale(1./nAnalyzedEvents);
363
364 return fYProtons;
365}
366
367//____________________________________________________________________//
368TH1D *AliProtonAnalysis::GetAntiProtonCorrectedYHistogram() {
369 Int_t nAnalyzedEvents = GetNumberOfAnalyzedEvents();
370
371 TH1D *fYAntiProtons = fCorrectAntiProtons->Project(0); //0: rapidity
372
373 fYAntiProtons->SetStats(kFALSE);
374 fYAntiProtons->GetYaxis()->SetTitle("(1/N_{events})(dN/dy)");
375 fYAntiProtons->GetXaxis()->SetTitle("y");
376 fYAntiProtons->SetTitle("dN/dy protons");
377 fYAntiProtons->SetMarkerStyle(kFullCircle);
378 fYAntiProtons->SetMarkerColor(4);
379 if(nAnalyzedEvents > 0)
380 fYAntiProtons->Scale(1./nAnalyzedEvents);
381
382 return fYAntiProtons;
383}
384
385//____________________________________________________________________//
386TH1D *AliProtonAnalysis::GetProtonCorrectedPtHistogram() {
387 Int_t nAnalyzedEvents = GetNumberOfAnalyzedEvents();
388
389 TH1D *fPtProtons = fCorrectProtons->Project(0); //0: rapidity
390
391 fPtProtons->SetStats(kFALSE);
392 fPtProtons->GetYaxis()->SetTitle("(1/N_{events})(dN/dP_{T})");
393 fPtProtons->GetXaxis()->SetTitle("P_{T} [GeV/c]");
394 fPtProtons->SetTitle("dN/dPt protons");
395 fPtProtons->SetMarkerStyle(kFullCircle);
396 fPtProtons->SetMarkerColor(4);
397 if(nAnalyzedEvents > 0)
398 fPtProtons->Scale(1./nAnalyzedEvents);
399
400 return fPtProtons;
401}
402
403//____________________________________________________________________//
404TH1D *AliProtonAnalysis::GetAntiProtonCorrectedPtHistogram() {
405 Int_t nAnalyzedEvents = GetNumberOfAnalyzedEvents();
406
407 TH1D *fPtAntiProtons = fCorrectAntiProtons->Project(0); //0: rapidity
408
409 fPtAntiProtons->SetStats(kFALSE);
410 fPtAntiProtons->GetYaxis()->SetTitle("(1/N_{events})(dN/dP_{T})");
411 fPtAntiProtons->GetXaxis()->SetTitle("P_{T} [GeV/c]");
412 fPtAntiProtons->SetTitle("dN/dPt antiprotons");
413 fPtAntiProtons->SetMarkerStyle(kFullCircle);
414 fPtAntiProtons->SetMarkerColor(4);
415 if(nAnalyzedEvents > 0)
416 fPtAntiProtons->Scale(1./nAnalyzedEvents);
417
418 return fPtAntiProtons;
419}
420
734d2c12 421//____________________________________________________________________//
422TH1D *AliProtonAnalysis::GetYRatioHistogram() {
3e23254a 423 //Returns the rapidity dependence of the ratio (uncorrected)
734d2c12 424 TH1D *fYProtons = GetProtonYHistogram();
425 TH1D *fYAntiProtons = GetAntiProtonYHistogram();
426
427 TH1D *hRatioY = new TH1D("hRatioY","",fYProtons->GetNbinsX(),fYProtons->GetXaxis()->GetXmin(),fYProtons->GetXaxis()->GetXmax());
428 hRatioY->Divide(fYAntiProtons,fYProtons,1.0,1.0);
429 hRatioY->SetMarkerStyle(kFullCircle);
430 hRatioY->SetMarkerColor(4);
431 hRatioY->GetYaxis()->SetTitle("#bar{p}/p");
432 hRatioY->GetYaxis()->SetTitleOffset(1.4);
433 hRatioY->GetXaxis()->SetTitle("y");
434 hRatioY->GetXaxis()->SetTitleColor(1);
435 hRatioY->SetStats(kFALSE);
436
437 return hRatioY;
438}
439
3e23254a 440//____________________________________________________________________//
441TH1D *AliProtonAnalysis::GetYRatioCorrectedHistogram(TH2D *gCorrectionProtons,
442 TH2D *gCorrectionAntiProtons) {
443 //Returns the rapidity dependence of the ratio (corrected)
444 fHistYPtProtons->Multiply(gCorrectionProtons);
445 TH1D *fYProtons = GetProtonYHistogram();
446 fHistYPtAntiProtons->Multiply(gCorrectionAntiProtons);
447 TH1D *fYAntiProtons = GetAntiProtonYHistogram();
448
449 TH1D *hRatioY = new TH1D("hRatioY","",fYProtons->GetNbinsX(),fYProtons->GetXaxis()->GetXmin(),fYProtons->GetXaxis()->GetXmax());
450 hRatioY->Divide(fYAntiProtons,fYProtons,1.0,1.0);
451 hRatioY->SetMarkerStyle(kFullCircle);
452 hRatioY->SetMarkerColor(4);
453 hRatioY->GetYaxis()->SetTitle("#bar{p}/p");
454 hRatioY->GetYaxis()->SetTitleOffset(1.4);
455 hRatioY->GetXaxis()->SetTitle("y");
456 hRatioY->GetXaxis()->SetTitleColor(1);
457 hRatioY->SetStats(kFALSE);
458
459 return hRatioY;
460}
461
734d2c12 462//____________________________________________________________________//
463TH1D *AliProtonAnalysis::GetPtRatioHistogram() {
3e23254a 464 //Returns the pT dependence of the ratio (uncorrected)
734d2c12 465 TH1D *fPtProtons = GetProtonPtHistogram();
466 TH1D *fPtAntiProtons = GetAntiProtonPtHistogram();
467
468 TH1D *hRatioPt = new TH1D("hRatioPt","",fPtProtons->GetNbinsX(),fPtProtons->GetXaxis()->GetXmin(),fPtProtons->GetXaxis()->GetXmax());
469 hRatioPt->Divide(fPtAntiProtons,fPtProtons,1.0,1.0);
470 hRatioPt->SetMarkerStyle(kFullCircle);
471 hRatioPt->SetMarkerColor(4);
472 hRatioPt->GetYaxis()->SetTitle("#bar{p}/p");
473 hRatioPt->GetYaxis()->SetTitleOffset(1.4);
474 hRatioPt->GetXaxis()->SetTitle("P_{T} [GeV/c]");
475 hRatioPt->GetXaxis()->SetTitleColor(1);
476 hRatioPt->SetStats(kFALSE);
477
478 return hRatioPt;
479}
480
3e23254a 481//____________________________________________________________________//
482TH1D *AliProtonAnalysis::GetPtRatioCorrectedHistogram(TH2D *gCorrectionProtons,
483 TH2D *gCorrectionAntiProtons) {
484 //Returns the Pt dependence of the ratio (corrected)
485 fHistYPtProtons->Multiply(gCorrectionProtons);
486 TH1D *fPtProtons = GetProtonPtHistogram();
487 fHistYPtAntiProtons->Multiply(gCorrectionAntiProtons);
488 TH1D *fPtAntiProtons = GetAntiProtonPtHistogram();
489
490 TH1D *hRatioPt = new TH1D("hRatioPt","",fPtProtons->GetNbinsX(),fPtProtons->GetXaxis()->GetXmin(),fPtProtons->GetXaxis()->GetXmax());
491 hRatioPt->Divide(fPtAntiProtons,fPtProtons,1.0,1.0);
492 hRatioPt->SetMarkerStyle(kFullCircle);
493 hRatioPt->SetMarkerColor(4);
494 hRatioPt->GetYaxis()->SetTitle("#bar{p}/p");
495 hRatioPt->GetYaxis()->SetTitleOffset(1.4);
496 hRatioPt->GetXaxis()->SetTitle("y");
497 hRatioPt->GetXaxis()->SetTitleColor(1);
498 hRatioPt->SetStats(kFALSE);
499
500 return hRatioPt;
501}
502
734d2c12 503//____________________________________________________________________//
504TH1D *AliProtonAnalysis::GetYAsymmetryHistogram() {
3e23254a 505 //Returns the rapidity dependence of the asymmetry (uncorrected)
734d2c12 506 TH1D *fYProtons = GetProtonYHistogram();
507 TH1D *fYAntiProtons = GetAntiProtonYHistogram();
508
509 TH1D *hsum = new TH1D("hsumY","",fYProtons->GetNbinsX(),fYProtons->GetXaxis()->GetXmin(),fYProtons->GetXaxis()->GetXmax());
510 hsum->Add(fYProtons,fYAntiProtons,1.0,1.0);
511
512 TH1D *hdiff = new TH1D("hdiffY","",fYProtons->GetNbinsX(),fYProtons->GetXaxis()->GetXmin(),fYProtons->GetXaxis()->GetXmax());
513 hdiff->Add(fYProtons,fYAntiProtons,1.0,-1.0);
514
515 TH1D *hAsymmetryY = new TH1D("hAsymmetryY","",fYProtons->GetNbinsX(),fYProtons->GetXaxis()->GetXmin(),fYProtons->GetXaxis()->GetXmax());
516 hAsymmetryY->Divide(hdiff,hsum,2.0,1.);
517 hAsymmetryY->SetMarkerStyle(kFullCircle);
518 hAsymmetryY->SetMarkerColor(4);
519 hAsymmetryY->GetYaxis()->SetTitle("A_{p}");
520 hAsymmetryY->GetYaxis()->SetTitleOffset(1.4);
521 hAsymmetryY->GetXaxis()->SetTitle("y");
522 hAsymmetryY->GetXaxis()->SetTitleColor(1);
523 hAsymmetryY->SetStats(kFALSE);
524
525 return hAsymmetryY;
526}
527
528//____________________________________________________________________//
529TH1D *AliProtonAnalysis::GetPtAsymmetryHistogram() {
3e23254a 530 //Returns the pT dependence of the asymmetry (uncorrected)
734d2c12 531 TH1D *fPtProtons = GetProtonPtHistogram();
532 TH1D *fPtAntiProtons = GetAntiProtonPtHistogram();
533
534 TH1D *hsum = new TH1D("hsumPt","",fPtProtons->GetNbinsX(),fPtProtons->GetXaxis()->GetXmin(),fPtProtons->GetXaxis()->GetXmax());
535 hsum->Add(fPtProtons,fPtAntiProtons,1.0,1.0);
536
537 TH1D *hdiff = new TH1D("hdiffPt","",fPtProtons->GetNbinsX(),fPtProtons->GetXaxis()->GetXmin(),fPtProtons->GetXaxis()->GetXmax());
538 hdiff->Add(fPtProtons,fPtAntiProtons,1.0,-1.0);
539
540 TH1D *hAsymmetryPt = new TH1D("hAsymmetryPt","",fPtProtons->GetNbinsX(),fPtProtons->GetXaxis()->GetXmin(),fPtProtons->GetXaxis()->GetXmax());
541 hAsymmetryPt->Divide(hdiff,hsum,2.0,1.);
542 hAsymmetryPt->SetMarkerStyle(kFullCircle);
543 hAsymmetryPt->SetMarkerColor(4);
544 hAsymmetryPt->GetYaxis()->SetTitle("A_{p}");
545 hAsymmetryPt->GetYaxis()->SetTitleOffset(1.4);
546 hAsymmetryPt->GetXaxis()->SetTitle("P_{T} [GeV/c]");
547 hAsymmetryPt->GetXaxis()->SetTitleColor(1);
548 hAsymmetryPt->SetStats(kFALSE);
549
550 return hAsymmetryPt;
551}
552
aafecd8b 553//____________________________________________________________________//
554Double_t AliProtonAnalysis::GetParticleFraction(Int_t i, Double_t p) {
8b8b0b7a 555 Double_t partFrac=0;
aafecd8b 556 if(fFunctionProbabilityFlag) {
557 if(i == 0) partFrac = fElectronFunction->Eval(p);
558 if(i == 1) partFrac = fMuonFunction->Eval(p);
559 if(i == 2) partFrac = fPionFunction->Eval(p);
560 if(i == 3) partFrac = fKaonFunction->Eval(p);
561 if(i == 4) partFrac = fProtonFunction->Eval(p);
562 }
563 else partFrac = fPartFrac[i];
564
565 return partFrac;
566}
567
734d2c12 568//____________________________________________________________________//
6667f3a7 569void AliProtonAnalysis::Analyze(AliESDEvent* esd,
570 const AliESDVertex *vertex) {
e4358d7f 571 //Main analysis part - ESD
3f6d0c08 572 fHistEvents->Fill(0); //number of analyzed events
251e4034 573 Double_t containerInput[2] ;
2b748670 574 Double_t Pt = 0.0, P = 0.0;
6667f3a7 575 Int_t nGoodTracks = esd->GetNumberOfTracks();
734d2c12 576 for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) {
6667f3a7 577 AliESDtrack* track = esd->GetTrack(iTracks);
2b748670 578 Double_t probability[5];
df1e8a78 579 AliESDtrack trackTPC;
2b748670 580
32a2ebb6 581 //in case it's a TPC only track relate it to the proper vertex
df1e8a78 582 if((fUseTPCOnly)&&(!fUseHybridTPC)) {
583 Float_t p[2],cov[3];
584 track->GetImpactParametersTPC(p,cov);
585 if (p[0]==0 && p[1]==0)
6667f3a7 586 track->RelateToVertexTPC(((AliESDEvent*)esd)->GetPrimaryVertexTPC(),esd->GetMagneticField(),kVeryBig);
df1e8a78 587 if (!track->FillTPCOnlyTrack(trackTPC)) {
32a2ebb6 588 continue;
df1e8a78 589 }
590 track = &trackTPC ;
591 }
592
6667f3a7 593 if(IsAccepted(esd,vertex,track)) {
2b748670 594 if(fUseTPCOnly) {
595 AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
596 if(!tpcTrack) continue;
597 Pt = tpcTrack->Pt();
598 P = tpcTrack->P();
599
600 //pid
601 track->GetTPCpid(probability);
602 Double_t rcc = 0.0;
603 for(Int_t i = 0; i < AliPID::kSPECIES; i++)
604 rcc += probability[i]*GetParticleFraction(i,P);
605 if(rcc == 0.0) continue;
606 Double_t w[5];
607 for(Int_t i = 0; i < AliPID::kSPECIES; i++)
608 w[i] = probability[i]*GetParticleFraction(i,P)/rcc;
609 Long64_t fParticleType = TMath::LocMax(AliPID::kSPECIES,w);
610 if(fParticleType == 4) {
251e4034 611 if(tpcTrack->Charge() > 0) {
612 fHistYPtProtons->Fill(Rapidity(tpcTrack->Px(),
613 tpcTrack->Py(),
614 tpcTrack->Pz()),
615 Pt);
616 //fill the container
617 containerInput[0] = Rapidity(tpcTrack->Px(),
618 tpcTrack->Py(),
619 tpcTrack->Pz());
620 containerInput[1] = Pt;
621 fProtonContainer->Fill(containerInput,0);
622 }//protons
623 else if(tpcTrack->Charge() < 0) {
624 fHistYPtAntiProtons->Fill(Rapidity(tpcTrack->Px(),
625 tpcTrack->Py(),
626 tpcTrack->Pz()),
627 Pt);
628 //fill the container
629 containerInput[0] = Rapidity(tpcTrack->Px(),
630 tpcTrack->Py(),
631 tpcTrack->Pz());
632 containerInput[1] = Pt;
633 fAntiProtonContainer->Fill(containerInput,0);
634 }//antiprotons
2b748670 635 }//proton check
636 }//TPC only tracks
637 else if(!fUseTPCOnly) {
638 Pt = track->Pt();
639 P = track->P();
734d2c12 640
2b748670 641 //pid
642 track->GetESDpid(probability);
643 Double_t rcc = 0.0;
644 for(Int_t i = 0; i < AliPID::kSPECIES; i++)
645 rcc += probability[i]*GetParticleFraction(i,P);
646 if(rcc == 0.0) continue;
647 Double_t w[5];
648 for(Int_t i = 0; i < AliPID::kSPECIES; i++)
649 w[i] = probability[i]*GetParticleFraction(i,P)/rcc;
650 Long64_t fParticleType = TMath::LocMax(AliPID::kSPECIES,w);
651 if(fParticleType == 4) {
652 //cout<<"(Anti)protons found..."<<endl;
251e4034 653 if(track->Charge() > 0) {
654 fHistYPtProtons->Fill(Rapidity(track->Px(),
655 track->Py(),
656 track->Pz()),
657 Pt);
658 //fill the container
659 containerInput[0] = Rapidity(track->Px(),
660 track->Py(),
661 track->Pz());
662 containerInput[1] = Pt;
663 fProtonContainer->Fill(containerInput,0);
664 }//protons
665 else if(track->Charge() < 0) {
666 fHistYPtAntiProtons->Fill(Rapidity(track->Px(),
667 track->Py(),
668 track->Pz()),
669 Pt);
670 //fill the container
671 containerInput[0] = Rapidity(track->Px(),
672 track->Py(),
673 track->Pz());
674 containerInput[1] = Pt;
675 fAntiProtonContainer->Fill(containerInput,0);
676 }//antiprotons
2b748670 677 }//proton check
678 }//combined tracking
734d2c12 679 }//cuts
680 }//track loop
681}
682
ee4ca40d 683//____________________________________________________________________//
684void AliProtonAnalysis::Analyze(AliAODEvent* fAOD) {
e4358d7f 685 //Main analysis part - AOD
3f6d0c08 686 fHistEvents->Fill(0); //number of analyzed events
ee4ca40d 687 Int_t nGoodTracks = fAOD->GetNumberOfTracks();
688 for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) {
689 AliAODTrack* track = fAOD->GetTrack(iTracks);
690 Double_t Pt = track->Pt();
691 Double_t P = track->P();
692
693 //pid
694 Double_t probability[10];
695 track->GetPID(probability);
696 Double_t rcc = 0.0;
697 for(Int_t i = 0; i < AliPID::kSPECIESN; i++) rcc += probability[i]*GetParticleFraction(i,P);
698 if(rcc == 0.0) continue;
738619fd 699 Double_t w[10];
ee4ca40d 700 for(Int_t i = 0; i < AliPID::kSPECIESN; i++) w[i] = probability[i]*GetParticleFraction(i,P)/rcc;
701 Long64_t fParticleType = TMath::LocMax(AliPID::kSPECIESN,w);
702 if(fParticleType == 4) {
251e4034 703 if(track->Charge() > 0)
704 fHistYPtProtons->Fill(track->Y(fParticleType),Pt);
705 else if(track->Charge() < 0)
706 fHistYPtAntiProtons->Fill(track->Y(fParticleType),Pt);
ee4ca40d 707 }//proton check
708 }//track loop
709}
710
e4358d7f 711//____________________________________________________________________//
da8c4c1c 712void AliProtonAnalysis::Analyze(AliStack* stack,
713 Bool_t iInclusive) {
e4358d7f 714 //Main analysis part - MC
3f6d0c08 715 fHistEvents->Fill(0); //number of analyzed events
3e23254a 716
717 Int_t nParticles = 0;
718 //inclusive protons -
719 if(iInclusive) nParticles = stack->GetNtrack();
720 else nParticles = stack->GetNprimary();
721
722 for(Int_t i = 0; i < nParticles; i++) {
e4358d7f 723 TParticle *particle = stack->Particle(i);
55f9a666 724 if(!particle) continue;
537afcc3 725
da8c4c1c 726 //in case of inclusive protons reject the secondaries from hadronic inter.
727 if(particle->GetUniqueID() == 13) continue;
728
537afcc3 729 if(TMath::Abs(particle->Eta()) > 1.0) continue;
730 if((particle->Pt() > fMaxPt)||(particle->Pt() < fMinPt)) continue;
731 if((Rapidity(particle->Px(),particle->Py(),particle->Pz()) > fMaxY)||(Rapidity(particle->Px(),particle->Py(),particle->Pz()) < fMinY)) continue;
732
e4358d7f 733 Int_t pdgcode = particle->GetPdgCode();
734 if(pdgcode == 2212) fHistYPtProtons->Fill(Rapidity(particle->Px(),
735 particle->Py(),
736 particle->Pz()),
737 particle->Pt());
738 if(pdgcode == -2212) fHistYPtAntiProtons->Fill(Rapidity(particle->Px(),
739 particle->Py(),
740 particle->Pz()),
741 particle->Pt());
251e4034 742 }//particle loop
e4358d7f 743}
744
3e23254a 745//____________________________________________________________________//
746Bool_t AliProtonAnalysis::IsInPhaseSpace(AliESDtrack* track) {
747 // Checks if the track is outside the analyzed y-Pt phase space
748 Double_t Pt = 0.0, Px = 0.0, Py = 0.0, Pz = 0.0;
dd3fa486 749 Double_t eta = 0.0;
750
3e23254a 751 if(fUseTPCOnly) {
752 AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
753 if(!tpcTrack) {
dd3fa486 754 Pt = 0.0; Px = 0.0; Py = 0.0; Pz = 0.0; eta = -10.0;
3e23254a 755 }
756 else {
757 Pt = tpcTrack->Pt();
758 Px = tpcTrack->Px();
759 Py = tpcTrack->Py();
760 Pz = tpcTrack->Pz();
dd3fa486 761 eta = tpcTrack->Eta();
3e23254a 762 }
763 }
764 else {
765 Pt = track->Pt();
766 Px = track->Px();
767 Py = track->Py();
768 Pz = track->Pz();
dd3fa486 769 eta = track->Eta();
3e23254a 770 }
771
772 if((Pt < fMinPt) || (Pt > fMaxPt)) return kFALSE;
dd3fa486 773 if(fAnalysisEtaMode) {
774 if((eta < fMinY) || (eta > fMaxY))
775 return kFALSE;
776 }
777 else {
778 if((Rapidity(Px,Py,Pz) < fMinY) || (Rapidity(Px,Py,Pz) > fMaxY))
779 return kFALSE;
780 }
3e23254a 781
782 return kTRUE;
783}
784
734d2c12 785//____________________________________________________________________//
6667f3a7 786Bool_t AliProtonAnalysis::IsAccepted(AliESDEvent *esd,
787 const AliESDVertex *vertex,
788 AliESDtrack* track) {
734d2c12 789 // Checks if the track is excluded from the cuts
2b748670 790 Double_t Pt = 0.0, Px = 0.0, Py = 0.0, Pz = 0.0;
6667f3a7 791 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0}; //The impact parameters and their covariance.
792
793 if((fUseTPCOnly)&&(!fUseHybridTPC)) {
794 AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
795 if(!tpcTrack) {
796 Pt = 0.0; Px = 0.0; Py = 0.0; Pz = 0.0;
797 dca[0] = -100.; dca[1] = -100.;
798 cov[0] = -100.; cov[1] = -100.; cov[2] = -100.;
799 }
800 else {
801 Pt = tpcTrack->Pt();
802 Px = tpcTrack->Px();
803 Py = tpcTrack->Py();
804 Pz = tpcTrack->Pz();
805 tpcTrack->PropagateToDCA(vertex,
806 esd->GetMagneticField(),
807 100.,dca,cov);
808 }
809 }
810 else if(fUseHybridTPC) {
55f9a666 811 AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
812 if(!tpcTrack) {
813 Pt = 0.0; Px = 0.0; Py = 0.0; Pz = 0.0;
6667f3a7 814 dca[0] = -100.; dca[1] = -100.;
815 cov[0] = -100.; cov[1] = -100.; cov[2] = -100.;
55f9a666 816 }
817 else {
818 Pt = tpcTrack->Pt();
819 Px = tpcTrack->Px();
820 Py = tpcTrack->Py();
821 Pz = tpcTrack->Pz();
6667f3a7 822 tpcTrack->PropagateToDCA(vertex,
823 esd->GetMagneticField(),
824 100.,dca,cov);
55f9a666 825 }
826 }
2b748670 827 else{
828 Pt = track->Pt();
829 Px = track->Px();
830 Py = track->Py();
831 Pz = track->Pz();
6667f3a7 832 track->PropagateToDCA(vertex,
833 esd->GetMagneticField(),
834 100.,dca,cov);
2b748670 835 }
836
734d2c12 837 Int_t fIdxInt[200];
838 Int_t nClustersITS = track->GetITSclusters(fIdxInt);
839 Int_t nClustersTPC = track->GetTPCclusters(fIdxInt);
840
841 Float_t chi2PerClusterITS = -1;
ef1a8dbd 842 if (nClustersITS!=0)
843 chi2PerClusterITS = track->GetITSchi2()/Float_t(nClustersITS);
844 Float_t chi2PerClusterTPC = -1;
845 if (nClustersTPC!=0)
846 chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC);
847
848 Double_t extCov[15];
849 track->GetExternalCovariance(extCov);
850
0008a5a6 851 if(fPointOnITSLayer1Flag)
852 if(!track->HasPointOnITSLayer(0)) return kFALSE;
853 if(fPointOnITSLayer2Flag)
854 if(!track->HasPointOnITSLayer(1)) return kFALSE;
855 if(fPointOnITSLayer3Flag)
856 if(!track->HasPointOnITSLayer(2)) return kFALSE;
857 if(fPointOnITSLayer4Flag)
858 if(!track->HasPointOnITSLayer(3)) return kFALSE;
859 if(fPointOnITSLayer5Flag)
860 if(!track->HasPointOnITSLayer(4)) return kFALSE;
861 if(fPointOnITSLayer6Flag)
862 if(!track->HasPointOnITSLayer(5)) return kFALSE;
ef1a8dbd 863 if(fMinITSClustersFlag)
864 if(nClustersITS < fMinITSClusters) return kFALSE;
f24a70b0 865 if(fMaxChi2PerITSClusterFlag)
866 if(chi2PerClusterITS > fMaxChi2PerITSCluster) return kFALSE;
ef1a8dbd 867 if(fMinTPCClustersFlag)
868 if(nClustersTPC < fMinTPCClusters) return kFALSE;
869 if(fMaxChi2PerTPCClusterFlag)
870 if(chi2PerClusterTPC > fMaxChi2PerTPCCluster) return kFALSE;
ef1a8dbd 871 if(fMaxCov11Flag)
872 if(extCov[0] > fMaxCov11) return kFALSE;
873 if(fMaxCov22Flag)
874 if(extCov[2] > fMaxCov22) return kFALSE;
875 if(fMaxCov33Flag)
876 if(extCov[5] > fMaxCov33) return kFALSE;
877 if(fMaxCov44Flag)
878 if(extCov[9] > fMaxCov44) return kFALSE;
879 if(fMaxCov55Flag)
880 if(extCov[14] > fMaxCov55) return kFALSE;
881 if(fMaxSigmaToVertexFlag)
882 if(GetSigmaToVertex(track) > fMaxSigmaToVertex) return kFALSE;
f24a70b0 883 if(fMaxSigmaToVertexTPCFlag)
884 if(GetSigmaToVertex(track) > fMaxSigmaToVertexTPC) return kFALSE;
3e6c06f4 885 if(fMaxDCAXYFlag)
6667f3a7 886 if(TMath::Abs(dca[0]) > fMaxDCAXY) return kFALSE;
3e6c06f4 887 if(fMaxDCAXYTPCFlag)
6667f3a7 888 if(TMath::Abs(dca[0]) > fMaxDCAXYTPC) return kFALSE;
3e6c06f4 889 if(fMaxDCAZFlag)
6667f3a7 890 if(TMath::Abs(dca[1]) > fMaxDCAZ) return kFALSE;
3e6c06f4 891 if(fMaxDCAZTPCFlag)
6667f3a7 892 if(TMath::Abs(dca[1]) > fMaxDCAZTPC) return kFALSE;
3e6c06f4 893 if(fMaxConstrainChi2Flag) {
894 if(track->GetConstrainedChi2() > 0)
895 if(TMath::Log(track->GetConstrainedChi2()) > fMaxConstrainChi2) return kFALSE;
896 }
ef1a8dbd 897 if(fITSRefitFlag)
898 if ((track->GetStatus() & AliESDtrack::kITSrefit) == 0) return kFALSE;
899 if(fTPCRefitFlag)
900 if ((track->GetStatus() & AliESDtrack::kTPCrefit) == 0) return kFALSE;
f24a70b0 901 if(fESDpidFlag)
902 if ((track->GetStatus() & AliESDtrack::kESDpid) == 0) return kFALSE;
903 if(fTPCpidFlag)
904 if ((track->GetStatus() & AliESDtrack::kTPCpid) == 0) return kFALSE;
ef1a8dbd 905
ef1a8dbd 906 return kTRUE;
907}
908
734d2c12 909//____________________________________________________________________//
910Float_t AliProtonAnalysis::GetSigmaToVertex(AliESDtrack* esdTrack) {
911 // Calculates the number of sigma to the vertex.
912
913 Float_t b[2];
914 Float_t bRes[2];
915 Float_t bCov[3];
55f9a666 916 if((fUseTPCOnly)&&(!fUseHybridTPC))
f2eeda10 917 esdTrack->GetImpactParametersTPC(b,bCov);
918 else
919 esdTrack->GetImpactParameters(b,bCov);
920
734d2c12 921 if (bCov[0]<=0 || bCov[2]<=0) {
922 //AliDebug(1, "Estimated b resolution lower or equal zero!");
923 bCov[0]=0; bCov[2]=0;
924 }
925 bRes[0] = TMath::Sqrt(bCov[0]);
926 bRes[1] = TMath::Sqrt(bCov[2]);
927
928 if (bRes[0] == 0 || bRes[1] ==0) return -1;
929
930 Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
931
932 if (TMath::Exp(-d * d / 2) < 1e-10) return 1000;
933
934 d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
935
936 return d;
937}
938
3f6d0c08 939//____________________________________________________________________//
2b748670 940Double_t AliProtonAnalysis::Rapidity(Double_t Px, Double_t Py, Double_t Pz) {
3f6d0c08 941 //returns the rapidity of the proton - to be removed
734d2c12 942 Double_t fMass = 9.38270000000000048e-01;
943
2b748670 944 Double_t P = TMath::Sqrt(TMath::Power(Px,2) +
945 TMath::Power(Py,2) +
946 TMath::Power(Pz,2));
734d2c12 947 Double_t energy = TMath::Sqrt(P*P + fMass*fMass);
948 Double_t y = -999;
2b748670 949 if(energy != Pz)
950 y = 0.5*TMath::Log((energy + Pz)/(energy - Pz));
734d2c12 951
952 return y;
953}
3f6d0c08 954
955//____________________________________________________________________//
956Bool_t AliProtonAnalysis::PrintMean(TH1 *hist, Double_t edge) {
957 //calculates the mean value of the ratio/asymmetry within \pm edge
958 Double_t sum = 0.0;
959 Int_t nentries = 0;
960 //calculate the mean
961 for(Int_t i = 0; i < hist->GetXaxis()->GetNbins(); i++) {
962 Double_t x = hist->GetBinCenter(i+1);
963 Double_t y = hist->GetBinContent(i+1);
964 if(TMath::Abs(x) < edge) {
965 sum += y;
966 nentries += 1;
967 }
968 }
969 Double_t mean = 0.0;
970 if(nentries != 0)
971 mean = sum/nentries;
972
973 //calculate the error
974 for(Int_t i = 0; i < hist->GetXaxis()->GetNbins(); i++) {
975 Double_t x = hist->GetBinCenter(i+1);
976 Double_t y = hist->GetBinContent(i+1);
977 if(TMath::Abs(x) < edge) {
978 sum += TMath::Power((mean - y),2);
979 nentries += 1;
980 }
981 }
982
983 Double_t error = 0.0;
984 if(nentries != 0)
985 error = TMath::Sqrt(sum)/nentries;
986
987 cout<<"========================================="<<endl;
988 cout<<"Input distribution: "<<hist->GetName()<<endl;
989 cout<<"Interval used: -"<<edge<<" -> "<<edge<<endl;
990 cout<<"Mean value :"<<mean<<endl;
991 cout<<"Error: "<<error<<endl;
992 cout<<"========================================="<<endl;
993
994 return 0;
995}
996
997//____________________________________________________________________//
998Bool_t AliProtonAnalysis::PrintYields(TH1 *hist, Double_t edge) {
999 //calculates the (anti)proton yields within the \pm edge
1000 Double_t sum = 0.0, sumerror = 0.0;
1001 Double_t error = 0.0;
1002 for(Int_t i = 0; i < hist->GetXaxis()->GetNbins(); i++) {
1003 Double_t x = hist->GetBinCenter(i+1);
1004 Double_t y = hist->GetBinContent(i+1);
1005 if(TMath::Abs(x) < edge) {
1006 sum += y;
1007 sumerror += TMath::Power(hist->GetBinError(i+1),2);
1008 }
1009 }
1010
1011 error = TMath::Sqrt(sumerror);
1012
1013 cout<<"========================================="<<endl;
1014 cout<<"Input distribution: "<<hist->GetName()<<endl;
1015 cout<<"Interval used: -"<<edge<<" -> "<<edge<<endl;
1016 cout<<"Yields :"<<sum<<endl;
1017 cout<<"Error: "<<error<<endl;
1018 cout<<"========================================="<<endl;
1019
1020 return 0;
1021}
1022
251e4034 1023//____________________________________________________________________//
1024void AliProtonAnalysis::Correct(Int_t step) {
1025 //Applies the correction maps to the initial containers
1026 fCorrectProtons = new AliCFDataGrid("correctProtons",
1027 "corrected data",
1028 *fProtonContainer);
1029 fCorrectProtons->ApplyEffCorrection(*(AliCFEffGrid *)fEffGridListProtons->At(step));
1030
1031 fCorrectAntiProtons = new AliCFDataGrid("correctAntiProtons",
1032 "corrected data",
1033 *fAntiProtonContainer);
1034 fCorrectAntiProtons->ApplyEffCorrection(*(AliCFEffGrid *)fEffGridListAntiProtons->At(step));
1035}
1036
39f2a708 1037//____________________________________________________________________//
1038Bool_t AliProtonAnalysis::ReadCorrectionContainer(const char* filename) {
1039 // Reads the outout of the correction framework task
1040 // Creates the correction maps
1041 // Puts the results in the different TList objects
1042 Bool_t status = kTRUE;
1043
1044 TFile *file = TFile::Open(filename);
1045 if(!file) {
1046 cout<<"Could not find the input CORRFW file "<<filename<<endl;
1047 status = kFALSE;
1048 }
1049
cdb3530f 1050 //________________________________________//
1051 //Protons
251e4034 1052 fEffGridListProtons = new TList();
cdb3530f 1053 fCorrectionListProtons2D = new TList();
1054 fEfficiencyListProtons1D = new TList();
1055 fCorrectionListProtons1D = new TList();
1056
1057 AliCFContainer *corrfwContainerProtons = (AliCFContainer*) (file->Get("containerProtons"));
1058 if(!corrfwContainerProtons) {
1059 cout<<"CORRFW container for protons not found!"<<endl;
39f2a708 1060 status = kFALSE;
1061 }
1062
cdb3530f 1063 Int_t nSteps = corrfwContainerProtons->GetNStep();
39f2a708 1064 TH2D *gYPt[4];
1065 //currently the GRID is formed by the y-pT parameters
1066 //Add Vz as a next step
1067 Int_t iRap = 0, iPt = 1;
9b168f47 1068 AliCFEffGrid *effProtonsStep0Step1 = new AliCFEffGrid("eff10",
1069 "effProtonsStep0Step1",
1070 *corrfwContainerProtons);
1071 effProtonsStep0Step1->CalculateEfficiency(1,0); //eff= step1/step0
251e4034 1072 fEffGridListProtons->Add(effProtonsStep0Step1);
9b168f47 1073 gYPt[0] = effProtonsStep0Step1->Project(iRap,iPt);
1074 fCorrectionListProtons2D->Add(gYPt[0]);
1075
1076 AliCFEffGrid *effProtonsStep0Step2 = new AliCFEffGrid("eff20",
1077 "effProtonsStep0Step2",
1078 *corrfwContainerProtons);
1079 effProtonsStep0Step2->CalculateEfficiency(2,0); //eff= step2/step0
251e4034 1080 fEffGridListProtons->Add(effProtonsStep0Step2);
9b168f47 1081 gYPt[1] = effProtonsStep0Step2->Project(iRap,iPt);
1082 fCorrectionListProtons2D->Add(gYPt[1]);
1083
1084 AliCFEffGrid *effProtonsStep0Step3 = new AliCFEffGrid("eff30",
1085 "effProtonsStep0Step3",
1086 *corrfwContainerProtons);
1087 effProtonsStep0Step3->CalculateEfficiency(3,0); //eff= step1/step0
251e4034 1088 fEffGridListProtons->Add(effProtonsStep0Step3);
9b168f47 1089 gYPt[2] = effProtonsStep0Step3->Project(iRap,iPt);
1090 fCorrectionListProtons2D->Add(gYPt[2]);
1091
1092 TH1D *gEfficiency[2][3]; //efficiency as a function of pT and of y (raws-[2])
1093 TH1D *gCorrection[2][3]; //efficiency as a function of pT and of y (raws-[2])
39f2a708 1094 TString gTitle = 0;
cdb3530f 1095 //Get the projection of the efficiency maps
9b168f47 1096 for(Int_t iParameter = 0; iParameter < 2; iParameter++) {
1097 gEfficiency[iParameter][0] = effProtonsStep0Step1->Project(iParameter);
1098 gTitle = "ProtonsEfficiency_Parameter"; gTitle += iParameter+1;
1099 gTitle += "_Step0_Step1";
1100 gEfficiency[iParameter][0]->SetName(gTitle.Data());
1101 fEfficiencyListProtons1D->Add(gEfficiency[iParameter][0]);
1102 gTitle = "ProtonsCorrection_Parameter"; gTitle += iParameter+1;
1103 gTitle += "_Step0_Step1";
1104 gCorrection[iParameter][0] = new TH1D(gTitle.Data(),
1105 gTitle.Data(),
1106 gEfficiency[iParameter][0]->GetNbinsX(),
1107 gEfficiency[iParameter][0]->GetXaxis()->GetXmin(),
1108 gEfficiency[iParameter][0]->GetXaxis()->GetXmax());
1109 //initialisation of the correction
1110 for(Int_t iBin = 1; iBin <= gEfficiency[iParameter][0]->GetNbinsX(); iBin++)
1111 gCorrection[iParameter][0]->SetBinContent(iBin,1.0);
1112
1113 gEfficiency[iParameter][1] = effProtonsStep0Step2->Project(iParameter);
1114 gTitle = "ProtonsEfficiency_Parameter"; gTitle += iParameter+1;
1115 gTitle += "_Step0_Step2";
1116 gEfficiency[iParameter][1]->SetName(gTitle.Data());
1117 fEfficiencyListProtons1D->Add(gEfficiency[iParameter][1]);
1118 gTitle = "ProtonsCorrection_Parameter"; gTitle += iParameter+1;
1119 gTitle += "_Step0_Step2";
1120 gCorrection[iParameter][1] = new TH1D(gTitle.Data(),
1121 gTitle.Data(),
1122 gEfficiency[iParameter][1]->GetNbinsX(),
1123 gEfficiency[iParameter][1]->GetXaxis()->GetXmin(),
1124 gEfficiency[iParameter][1]->GetXaxis()->GetXmax());
1125 //initialisation of the correction
1126 for(Int_t iBin = 1; iBin <= gEfficiency[iParameter][1]->GetNbinsX(); iBin++)
1127 gCorrection[iParameter][1]->SetBinContent(iBin,1.0);
1128
1129 gEfficiency[iParameter][2] = effProtonsStep0Step3->Project(iParameter);
1130 gTitle = "ProtonsEfficiency_Parameter"; gTitle += iParameter+1;
1131 gTitle += "_Step0_Step3";
1132 gEfficiency[iParameter][2]->SetName(gTitle.Data());
1133 fEfficiencyListProtons1D->Add(gEfficiency[iParameter][2]);
1134 gTitle = "ProtonsCorrection_Parameter"; gTitle += iParameter+1;
1135 gTitle += "_Step0_Step3";
1136 gCorrection[iParameter][2] = new TH1D(gTitle.Data(),
1137 gTitle.Data(),
1138 gEfficiency[iParameter][2]->GetNbinsX(),
1139 gEfficiency[iParameter][2]->GetXaxis()->GetXmin(),
1140 gEfficiency[iParameter][2]->GetXaxis()->GetXmax());
1141 //initialisation of the correction
1142 for(Int_t iBin = 1; iBin <= gEfficiency[iParameter][2]->GetNbinsX(); iBin++)
1143 gCorrection[iParameter][2]->SetBinContent(iBin,1.0);
cdb3530f 1144 }//parameter loop
1145 //Calculate the 1D correction parameters as a function of y and pT
1146 for(Int_t iParameter = 0; iParameter < 2; iParameter++) {
1147 for(Int_t iStep = 1; iStep < nSteps; iStep++) {
1148 gCorrection[iParameter][iStep-1]->Divide(gEfficiency[iParameter][iStep-1]);
1149 fCorrectionListProtons1D->Add(gCorrection[iParameter][iStep-1]);
1150 }
1151 }
1152
1153 //________________________________________//
1154 //AntiProtons
251e4034 1155 fEffGridListAntiProtons = new TList();
cdb3530f 1156 fCorrectionListAntiProtons2D = new TList();
1157 fEfficiencyListAntiProtons1D = new TList();
1158 fCorrectionListAntiProtons1D = new TList();
1159
1160 AliCFContainer *corrfwContainerAntiProtons = (AliCFContainer*) (file->Get("containerAntiProtons"));
1161 if(!corrfwContainerAntiProtons) {
1162 cout<<"CORRFW container for antiprotons not found!"<<endl;
1163 status = kFALSE;
1164 }
1165
1166 nSteps = corrfwContainerAntiProtons->GetNStep();
1167 //currently the GRID is formed by the y-pT parameters
1168 //Add Vz as a next step
9b168f47 1169 AliCFEffGrid *effAntiProtonsStep0Step1 = new AliCFEffGrid("eff10",
1170 "effAntiProtonsStep0Step1",
1171 *corrfwContainerAntiProtons);
1172 effAntiProtonsStep0Step1->CalculateEfficiency(1,0); //eff= step1/step0
251e4034 1173 fEffGridListAntiProtons->Add(effAntiProtonsStep0Step1);
9b168f47 1174 gYPt[0] = effAntiProtonsStep0Step1->Project(iRap,iPt);
1175 fCorrectionListAntiProtons2D->Add(gYPt[0]);
1176
1177 AliCFEffGrid *effAntiProtonsStep0Step2 = new AliCFEffGrid("eff20",
1178 "effAntiProtonsStep0Step2",
1179 *corrfwContainerAntiProtons);
1180 effAntiProtonsStep0Step2->CalculateEfficiency(2,0); //eff= step2/step0
251e4034 1181 fEffGridListAntiProtons->Add(effAntiProtonsStep0Step2);
9b168f47 1182 gYPt[1] = effAntiProtonsStep0Step2->Project(iRap,iPt);
1183 fCorrectionListAntiProtons2D->Add(gYPt[1]);
1184
1185 AliCFEffGrid *effAntiProtonsStep0Step3 = new AliCFEffGrid("eff30",
1186 "effAntiProtonsStep0Step3",
1187 *corrfwContainerAntiProtons);
1188 effAntiProtonsStep0Step3->CalculateEfficiency(3,0); //eff= step1/step0
251e4034 1189 fEffGridListAntiProtons->Add(effAntiProtonsStep0Step3);
9b168f47 1190 gYPt[2] = effAntiProtonsStep0Step3->Project(iRap,iPt);
1191 fCorrectionListAntiProtons2D->Add(gYPt[2]);
cdb3530f 1192
39f2a708 1193 //Get the projection of the efficiency maps
9b168f47 1194 for(Int_t iParameter = 0; iParameter < 2; iParameter++) {
1195 gEfficiency[iParameter][0] = effAntiProtonsStep0Step1->Project(iParameter);
1196 gTitle = "AntiProtonsEfficiency_Parameter"; gTitle += iParameter+1;
1197 gTitle += "_Step0_Step1";
1198 gEfficiency[iParameter][0]->SetName(gTitle.Data());
1199 fEfficiencyListAntiProtons1D->Add(gEfficiency[iParameter][0]);
1200 gTitle = "AntiProtonsCorrection_Parameter"; gTitle += iParameter+1;
1201 gTitle += "_Step0_Step1";
1202 gCorrection[iParameter][0] = new TH1D(gTitle.Data(),
1203 gTitle.Data(),
1204 gEfficiency[iParameter][0]->GetNbinsX(),
1205 gEfficiency[iParameter][0]->GetXaxis()->GetXmin(),
1206 gEfficiency[iParameter][0]->GetXaxis()->GetXmax());
1207 //initialisation of the correction
1208 for(Int_t iBin = 1; iBin <= gEfficiency[iParameter][0]->GetNbinsX(); iBin++)
1209 gCorrection[iParameter][0]->SetBinContent(iBin,1.0);
1210
1211 gEfficiency[iParameter][1] = effAntiProtonsStep0Step2->Project(iParameter);
1212 gTitle = "AntiProtonsEfficiency_Parameter"; gTitle += iParameter+1;
1213 gTitle += "_Step0_Step2";
1214 gEfficiency[iParameter][1]->SetName(gTitle.Data());
1215 fEfficiencyListAntiProtons1D->Add(gEfficiency[iParameter][1]);
1216 gTitle = "AntiProtonsCorrection_Parameter"; gTitle += iParameter+1;
1217 gTitle += "_Step0_Step2";
1218 gCorrection[iParameter][1] = new TH1D(gTitle.Data(),
1219 gTitle.Data(),
1220 gEfficiency[iParameter][1]->GetNbinsX(),
1221 gEfficiency[iParameter][1]->GetXaxis()->GetXmin(),
1222 gEfficiency[iParameter][1]->GetXaxis()->GetXmax());
1223 //initialisation of the correction
1224 for(Int_t iBin = 1; iBin <= gEfficiency[iParameter][1]->GetNbinsX(); iBin++)
1225 gCorrection[iParameter][1]->SetBinContent(iBin,1.0);
1226
1227 gEfficiency[iParameter][2] = effAntiProtonsStep0Step3->Project(iParameter);
1228 gTitle = "AntiProtonsEfficiency_Parameter"; gTitle += iParameter+1;
1229 gTitle += "_Step0_Step3";
1230 gEfficiency[iParameter][2]->SetName(gTitle.Data());
1231 fEfficiencyListAntiProtons1D->Add(gEfficiency[iParameter][2]);
1232 gTitle = "AntiProtonsCorrection_Parameter"; gTitle += iParameter+1;
1233 gTitle += "_Step0_Step3";
1234 gCorrection[iParameter][2] = new TH1D(gTitle.Data(),
1235 gTitle.Data(),
1236 gEfficiency[iParameter][2]->GetNbinsX(),
1237 gEfficiency[iParameter][2]->GetXaxis()->GetXmin(),
1238 gEfficiency[iParameter][2]->GetXaxis()->GetXmax());
1239 //initialisation of the correction
1240 for(Int_t iBin = 1; iBin <= gEfficiency[iParameter][2]->GetNbinsX(); iBin++)
1241 gCorrection[iParameter][2]->SetBinContent(iBin,1.0);
39f2a708 1242 }//parameter loop
1243 //Calculate the 1D correction parameters as a function of y and pT
1244 for(Int_t iParameter = 0; iParameter < 2; iParameter++) {
1245 for(Int_t iStep = 1; iStep < nSteps; iStep++) {
1246 gCorrection[iParameter][iStep-1]->Divide(gEfficiency[iParameter][iStep-1]);
cdb3530f 1247 fCorrectionListAntiProtons1D->Add(gCorrection[iParameter][iStep-1]);
39f2a708 1248 }
1249 }
ef1a8dbd 1250
1251 return status;
39f2a708 1252}
24421eb6 1253
39f2a708 1254
1255
1256
1257
1258
1259
3f6d0c08 1260
1261