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