]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGLF/STRANGENESS/Hypernuclei/AliAnalysisTaskAntiHe4.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
added snippets to save histos in header
[u/mrichter/AliRoot.git] / PWGLF / STRANGENESS / Hypernuclei / AliAnalysisTaskAntiHe4.cxx
CommitLineData
06ea1204 1#include <Riostream.h>
110805a7 2#include "TChain.h"
3#include "TTree.h"
4#include "TF1.h"
5#include "TH1F.h"
6#include "TH2F.h"
7#include "TH3F.h"
8#include "THn.h"
9#include "TCanvas.h"
10#include "AliESDtrackCuts.h"
11#include "AliESDpid.h"
12#include "AliESDVertex.h"
13#include "TFile.h"
14
15#include "AliAnalysisTask.h"
16#include "AliAnalysisManager.h"
17#include "AliCentrality.h"
18#include "AliESDEvent.h"
19#include "AliESDInputHandler.h"
20
d1f5b077 21#include "AliMCEventHandler.h"
22#include "AliMCEvent.h"
23#include "AliStack.h"
24
110805a7 25#include "AliAnalysisTaskAntiHe4.h"
26
f30d3e2c 27using std::cout;
28using std::endl;
110805a7 29
30///////////////////////////////////////////////////////////
31// //
32// Analysis for the observation of anti-alpha particles. //
33// //
34///////////////////////////////////////////////////////////
35
36
37ClassImp(AliAnalysisTaskAntiHe4)
38
39//________________________________________________________________________
40AliAnalysisTaskAntiHe4::AliAnalysisTaskAntiHe4()
41 : AliAnalysisTaskSE(),
42 fEventHandler(0),
43 fESD(0),
44 fHistCentralityClass10(0),
45 fHistCentralityPercentile(0),
46 fHistTriggerStat(0),
47 fHistTriggerStatAfterEventSelection(0),
48 fHistDEDx(0),
49 fHistTOF3D(0),
50 fHistAlpha(0),
51 fHistAlphaSignal(0),
52 fGraphAlphaSignal(0),
53 fNCounter(0),
54 fHistDeDx(0),
55 fHistDeDxRegion(0),
56 fHistDeDxSharp(0),
57 fHistTOF2D(0),
58 fHistTOFnuclei(0x0),
59 fNTriggers(5),
60 fBBParametersLightParticles(),
61 fBBParametersNuclei(),
62 fMCtrue(0),
63 fTriggerFired(),
64 fESDtrackCuts(0),
65 fESDtrackCutsSharp(0),
66 fESDpid(0),
67 fAntiAlpha(0),
d1f5b077 68 fHistHelium4PtGen(0),
69 fHistHelium4PtGenPrim(0),
70 fHistHelium4PtGenSec(0),
71 fHistHelium4PtGenEta(0),
72 fHistHelium4PtGenPrimEta(0),
73 fHistAntiHelium4PtGen(0),
74 fHistAntiHelium4PtGenPrim(0),
75 fHistAntiHelium4PtGenSec(0),
76 fHistAntiHelium4PtGenEta(0),
77 fHistHelium4PtAso(0),
78 fHistHelium4PtAsoPrim(0),
79 fHistHelium4PtAsoSec(0),
80 fHistAntiHelium4PtAso(0),
110805a7 81 fTree(0),
439dd8c8 82 fOutputContainer(0),
83 fEvnt(0),
84 fItrk(0)
110805a7 85{
86 // default Constructor
87
88
89 // Define input and output slots here
90}
91
92//________________________________________________________________________
93AliAnalysisTaskAntiHe4::AliAnalysisTaskAntiHe4(const char *name)
94 : AliAnalysisTaskSE(name),
95 fEventHandler(0),
96 fESD(0),
97 fHistCentralityClass10(0),
98 fHistCentralityPercentile(0),
99 fHistTriggerStat(0),
100 fHistTriggerStatAfterEventSelection(0),
101 fHistDEDx(0),
102 fHistTOF3D(0),
103 fHistAlpha(0),
104 fHistAlphaSignal(0),
105 fGraphAlphaSignal(0),
106 fNCounter(0),
107 fHistDeDx(0),
108 fHistDeDxRegion(0),
109 fHistDeDxSharp(0),
110 fHistTOF2D(0),
111 fHistTOFnuclei(0x0),
112 fNTriggers(5),
113 fBBParametersLightParticles(),
114 fBBParametersNuclei(),
115 fMCtrue(0),
116 fTriggerFired(),
117 fESDtrackCuts(0),
118 fESDtrackCutsSharp(0),
119 fESDpid(0),
120 fAntiAlpha(0),
d1f5b077 121 fHistHelium4PtGen(0),
122 fHistHelium4PtGenPrim(0),
123 fHistHelium4PtGenSec(0),
124 fHistHelium4PtGenEta(0),
125 fHistHelium4PtGenPrimEta(0),
126 fHistAntiHelium4PtGen(0),
127 fHistAntiHelium4PtGenPrim(0),
128 fHistAntiHelium4PtGenSec(0),
129 fHistAntiHelium4PtGenEta(0),
130 fHistHelium4PtAso(0),
131 fHistHelium4PtAsoPrim(0),
132 fHistHelium4PtAsoSec(0),
133 fHistAntiHelium4PtAso(0),
110805a7 134 fTree(0),
439dd8c8 135 fOutputContainer(0),
136 fEvnt(0),
137 fItrk(0)
110805a7 138{
139 // Constructor
140
141 // Define input and output slots here
142 // Input slot #0 works with a TChain
143 DefineInput(0, TChain::Class());
144 // Output slot #0 writes into a TH1 container
145 DefineOutput(1, TObjArray::Class());
146 DefineOutput(2, TTree::Class());
147 //
148 // cuts for candidates
149 //
150 fESDtrackCuts = new AliESDtrackCuts("AliESDtrackCuts","AliESDtrackCuts");
151 //
152 fESDtrackCuts->SetAcceptKinkDaughters(kFALSE);
153 fESDtrackCuts->SetMinNClustersTPC(70);
154 fESDtrackCuts->SetMaxChi2PerClusterTPC(5);
155 fESDtrackCuts->SetMaxDCAToVertexXY(3);
156 fESDtrackCuts->SetMaxDCAToVertexZ(2);
157 fESDtrackCuts->SetRequireTPCRefit(kTRUE);
158 //fESDtrackCuts->SetRequireITSRefit(kTRUE);
159 fESDtrackCuts->SetMinNClustersITS(2);
439dd8c8 160 fESDtrackCuts->SetEtaRange(-1.0,1.0);
110805a7 161 //
162 // cuts for final plots
163 //
164 fESDtrackCutsSharp = new AliESDtrackCuts("AliESDtrackCuts","AliESDtrackCuts");
165 fESDtrackCutsSharp->SetAcceptKinkDaughters(kFALSE);
166 fESDtrackCutsSharp->SetMinNClustersTPC(80);
167 fESDtrackCutsSharp->SetMaxChi2PerClusterITS(10);// TO BE INVESTIGATED !!!!!!!!!!!!!!
168 fESDtrackCutsSharp->SetMaxChi2PerClusterTPC(5);
169 fESDtrackCutsSharp->SetRequireTPCRefit(kTRUE);
170 fESDtrackCutsSharp->SetRequireITSRefit(kTRUE);
171 fESDtrackCutsSharp->SetMinNClustersITS(2);
172 fESDtrackCutsSharp->SetMaxDCAToVertexXY(0.1);
173 fESDtrackCutsSharp->SetMaxDCAToVertexZ(0.5);
174 fESDtrackCutsSharp->SetEtaRange(-0.8,0.8);
175
176 //ESD Track cuts from TestFilterRawTask
177 /* fESDtrackCuts = new AliESDtrackCuts("AliESDtrackCuts","AliESDtrackCuts");
178 fESDtrackCuts->SetMinNClustersTPC(80);
179 fESDtrackCuts->SetMinNClustersITS(2);
180 fESDtrackCuts->SetAcceptKinkDaughters(kFALSE);
181 fESDtrackCuts->SetMaxChi2PerClusterTPC(5);
182 fESDtrackCuts->SetRequireTPCRefit(kTRUE);
183 fESDtrackCuts->SetRequireITSRefit(kTRUE);
184 fESDtrackCuts->SetEtaRange(-1,1);
185 fESDtrackCuts->SetMaxDCAToVertexXY(1);
186 fESDtrackCuts->SetMaxDCAToVertexZ(2);
187 //test strickter cuts
188 //fESDtrackCuts->SetMaxDCAToVertexXY(0.1);
189 //fESDtrackCuts->SetMaxDCAToVertexZ(0.5);
190 //fESDtrackCuts->SetEtaRange(-0.8,0.8);
191 */
192
193 fMCtrue = kTRUE;
194
195}
196
197//________________________________________________________________________
198void AliAnalysisTaskAntiHe4::UserCreateOutputObjects()
199{
200 // Create histograms
201 // Called once
202 //
203 //histogram to count number of events
204 //
205 fHistCentralityClass10 = new TH1F("fHistCentralityClass10", "centrality with class10", 11, -0.5, 10.5);
206 fHistCentralityClass10->GetXaxis()->SetTitle("Centrality");
207 fHistCentralityClass10->GetYaxis()->SetTitle("Entries");
208 //
209 fHistCentralityPercentile = new TH1F("fHistCentralityPercentile", "centrality with percentile", 101, -0.1, 100.1);
210 fHistCentralityPercentile->GetXaxis()->SetTitle("Centrality");
211 fHistCentralityPercentile->GetYaxis()->SetTitle("Entries");
212 //
213 //trigger statitics histogram
214 //
215 fHistTriggerStat = new TH1F("fHistTriggerStat","Trigger statistics", fNTriggers,-0.5,fNTriggers-0.5);
216 const Char_t* aTriggerNames[] = { "kMB", "kCentral", "kSemiCentral", "kEMCEJE", "kEMCEGA" };
217 for ( Int_t ii=0; ii < fNTriggers; ii++ )
218 fHistTriggerStat->GetXaxis()->SetBinLabel(ii+1, aTriggerNames[ii]);
219
220 fHistTriggerStatAfterEventSelection = new TH1F("fHistTriggerStatAfterEventSelection","Trigger statistics after event selection", fNTriggers,-0.5,fNTriggers-0.5);
221 for ( Int_t ii=0; ii < fNTriggers; ii++ )
222 fHistTriggerStatAfterEventSelection->GetXaxis()->SetBinLabel(ii+1, aTriggerNames[ii]);
223
224 //dE/dx performance
225 fHistDEDx= new TH3F("fHistDEDx","DEDx",1000,0.01,100,1000,1,2000,2,-2,2);
226 BinLogAxis(fHistDEDx, 0);
227 BinLogAxis(fHistDEDx, 1);
228 fHistDEDx->GetXaxis()->SetTitle("p_{tot}/sign");
229 fHistDEDx->GetYaxis()->SetTitle("TPC signal");
230
231 fHistDeDx = new TH2F("fHistDeDx", "dE/dx", 1000, 0.1, 6.0, 1000, 0.0, 1000);
232 fHistDeDx->GetYaxis()->SetTitle("TPC dE/dx signal (a.u.)");
233 fHistDeDx->GetXaxis()->SetTitle("#frac{p}{z} (GeV/c)");
234 BinLogAxis(fHistDeDx);
235
236 fHistDeDxRegion = new TH3F("fHistDeDxRegion", "dE/dx", 400, 0., 6.0, 300, 0., 3., 4, -0.5, 3.5);
237 fHistDeDxRegion->GetYaxis()->SetTitle("TPC dE/dx signal (a.u.)");
238 fHistDeDxRegion->GetXaxis()->SetTitle("#frac{p}{z} (GeV/c)");
239
240 fHistDeDxSharp = new TH2F("fHistDeDxSharp", "dE/dx", 1000, 0.1, 6.0, 1000, 0.0, 1000);
241 fHistDeDxSharp->GetYaxis()->SetTitle("TPC dE/dx signal (a.u.)");
242 fHistDeDxSharp->GetXaxis()->SetTitle("#frac{p}{z} (GeV/c)");
243 BinLogAxis(fHistDeDxSharp);
244
245 //TOF performance
246 fHistTOF3D = new TH3F("fHistTOF3D","mass determination from TOF",400,0.25,4.5,2,-0.5,1.5,2,-1,1);
247
248 fHistTOF2D = new TH2F("fHistTOF2D", "TOF2D", 500, 0.0, 10., 2250, 0.2, 1.1);
249 fHistTOF2D->GetYaxis()->SetTitle("#beta");
250 fHistTOF2D->GetXaxis()->SetTitle("p (GeV/c)");
251
252 fHistTOFnuclei = new TH2F("fHistTOFnuclei", "TOF", 500, 0.0, 10., 2250, 0.7, 1.);
253 fHistTOFnuclei->GetYaxis()->SetTitle("#beta");
254 fHistTOFnuclei->GetXaxis()->SetTitle("#frac{p}{z} (GeV/c)");
255
256 //alphas
257 fHistAlpha = new TH1F("fHistAlpha", "Anti-Alpha", 22, 1.12, 4.31);
258 fHistAlpha->GetYaxis()->SetTitle("Counts");
259 fHistAlpha->GetXaxis()->SetTitle("#frac{m^{2}}{z^{2}} (GeV^{2}/c^{4})");
260
261 fHistAlphaSignal = new TH1F("fHistAlphaSignal", "Anti-Alpha", 22, 1.12, 4.31);
262 fHistAlphaSignal->GetYaxis()->SetTitle("Counts");
263 fHistAlphaSignal->GetXaxis()->SetTitle("#frac{m^{2}}{z^{2}} (GeV^{2}/c^{4})");
264
265 fGraphAlphaSignal = new TGraph(20);
266 fNCounter = 0;
267 //
268 // (0.) dca, (1.) sign, (2.) particle Type, (3.) p_tot
269 //
270 TString axisNameMult[4]={"dca","sign","particleType","ptot"};
271 TString axisTitleMult[4]={"dca","sign","particleType","ptot"};
272 const Int_t kDcaBins = 76;
273 Double_t binsDca[kDcaBins+1] = {-3,-2.85,-2.7,-2.55,-2.4,-2.25,-2.1,-1.95,-1.8,-1.65,-1.5,-1.35,-1.2,-1.05,-0.9,-0.75,-0.6,-0.45,-0.3,-0.285,-0.27,-0.255,-0.24,-0.225,-0.21,-0.195,-0.18,-0.165,-0.15,-0.135,-0.12,-0.105,-0.09,-0.075,-0.06,-0.045,-0.03,-0.015,0,0.015,0.03,0.045,0.06,0.075,0.09,0.105,0.12,0.135,0.15,0.165,0.18,0.195,0.21,0.225,0.24,0.255,0.27,0.285,0.3,0.45,0.6,0.75,0.9,1.05,1.2,1.35,1.5,1.65,1.8,1.95,2.1,2.25,2.4,2.55,2.7,2.85,3};
274 //
275 // (0.) dca, (1.) sign, (2.) particle Type, (3.) p_tot
276 Int_t binsAntiAlpha[4] = {77, 2, 4, 100};
277 Double_t xminAntiAlpha[4] = {-3, -2, -0.5, 0};
278 Double_t xmaxAntiAlpha[4] = { 3, 2, 3.5, 6};
279 //
280 fAntiAlpha = new THnF("fAntiAlpha", "AntiAlpha; (0.) dca, (1.) sign, (2.) particle Type, (3.) p_tot", 4, binsAntiAlpha, xminAntiAlpha, xmaxAntiAlpha);
281 //
282 fAntiAlpha->GetAxis(0)->Set(kDcaBins, binsDca);
283 for (Int_t iaxis=0; iaxis<4;iaxis++){
284 fAntiAlpha->GetAxis(iaxis)->SetName(axisNameMult[iaxis]);
285 fAntiAlpha->GetAxis(iaxis)->SetTitle(axisTitleMult[iaxis]);
286 }
287 //
d1f5b077 288 //
289 //Create histograms for MC
290 //generated histogramms
291 fHistHelium4PtGen = new TH1F("fHistHelium4PtGen", "Generated ^{4}He", 200, 0.0, 10.0);
292 fHistHelium4PtGen->GetYaxis()->SetTitle("Counts");
293 fHistHelium4PtGen->GetXaxis()->SetTitle("p_{T} (GeV/c)");
294
295 fHistHelium4PtGenPrim = new TH1F("fHistHelium4PtGenPrim", "Primary generated ^{4}He", 200, 0.0, 10.0);
296 fHistHelium4PtGenPrim->GetYaxis()->SetTitle("Counts");
297 fHistHelium4PtGenPrim->GetXaxis()->SetTitle("p_{T} (GeV/c)");
298
299 fHistHelium4PtGenSec = new TH1F("fHistHelium4PtGenSec", "Secondary generated ^{4}He", 200, 0.0, 10.0);
300 fHistHelium4PtGenSec->GetYaxis()->SetTitle("Counts");
301 fHistHelium4PtGenSec->GetXaxis()->SetTitle("p_{T} (GeV/c)");
302
303 fHistHelium4PtGenEta = new TH1F("fHistHelium4PtGenEta", "Generated ^{4}He in #eta < 0.8", 200, 0.0, 10.0);
304 fHistHelium4PtGenEta->GetYaxis()->SetTitle("Counts");
305 fHistHelium4PtGenEta->GetXaxis()->SetTitle("p_{T} (GeV/c)");
306
307 fHistHelium4PtGenPrimEta = new TH1F("fHistHelium4PtGenPrimEta", "Primary generated ^{4}He in #eta < 0.8", 200, 0.0, 10.0);
308 fHistHelium4PtGenPrimEta->GetYaxis()->SetTitle("Counts");
309 fHistHelium4PtGenPrimEta->GetXaxis()->SetTitle("p_{T} (GeV/c)");
310
311 fHistAntiHelium4PtGen = new TH1F("fHistAntiHelium4PtGen", "Generated #bar{^{4}He}", 200, 0.0, 10.0);
312 fHistAntiHelium4PtGen->GetYaxis()->SetTitle("Counts");
313 fHistAntiHelium4PtGen->GetXaxis()->SetTitle("p_{T} (GeV/c)");
314
315 fHistAntiHelium4PtGenPrim = new TH1F("fHistAntiHelium4PtGenPrim", "Primary generated #bar{^{4}He}", 200, 0.0, 10.0);
316 fHistAntiHelium4PtGenPrim->GetYaxis()->SetTitle("Counts");
317 fHistAntiHelium4PtGenPrim->GetXaxis()->SetTitle("p_{T} (GeV/c)");
318
319 fHistAntiHelium4PtGenSec = new TH1F("fHistAntiHelium4PtGenSec", "Secondary generated #bar{^{4}He}", 200, 0.0, 10.0);
320 fHistAntiHelium4PtGenSec->GetYaxis()->SetTitle("Counts");
321 fHistAntiHelium4PtGenSec->GetXaxis()->SetTitle("p_{T} (GeV/c)");
322
323 fHistAntiHelium4PtGenEta = new TH1F("fHistAntiHelium4PtGenEta", "Generated #bar{^{4}He} in #eta < 0.8", 200, 0.0, 10.0);
324 fHistAntiHelium4PtGenEta->GetYaxis()->SetTitle("Counts");
325 fHistAntiHelium4PtGenEta->GetXaxis()->SetTitle("p_{T} (GeV/c)");
326
327 //associated histograms
328 fHistHelium4PtAso = new TH1F("fHistHelium4PtAso", "Associated ^{4}He", 200, 0.0, 10.0);
329 fHistHelium4PtAso->GetYaxis()->SetTitle("Counts");
330 fHistHelium4PtAso->GetXaxis()->SetTitle("p_{T} (GeV/c)");
331
332 fHistHelium4PtAsoPrim = new TH1F("fHistHelium4PtAsoPrim", "Associated prim ^{4}He", 200, 0.0, 10.0);
333 fHistHelium4PtAsoPrim->GetYaxis()->SetTitle("Counts");
334 fHistHelium4PtAsoPrim->GetXaxis()->SetTitle("p_{T} (GeV/c)");
335
336 fHistHelium4PtAsoSec = new TH1F("fHistHelium4PtAsoSec", "Associated sec ^{4}He", 200, 0.0, 10.0);
337 fHistHelium4PtAsoSec->GetYaxis()->SetTitle("Counts");
338 fHistHelium4PtAsoSec->GetXaxis()->SetTitle("p_{T} (GeV/c)");
339
340 fHistAntiHelium4PtAso = new TH1F("fHistAntiHelium4PtAso", "Associated #bar{^{4}He}", 200, 0.0, 10.0);
341 fHistAntiHelium4PtAso->GetYaxis()->SetTitle("Counts");
342 fHistAntiHelium4PtAso->GetXaxis()->SetTitle("p_{T} (GeV/c)");
343
110805a7 344 //
345 //
346 fOutputContainer = new TObjArray(1);
347 fOutputContainer->SetOwner(kTRUE);
348 fOutputContainer->SetName(GetName());
349 fOutputContainer->Add(fHistCentralityClass10);
350 fOutputContainer->Add(fHistCentralityPercentile);
351 fOutputContainer->Add(fHistTriggerStat);
352 fOutputContainer->Add(fHistTriggerStatAfterEventSelection);
353 fOutputContainer->Add(fHistDEDx);
354 fOutputContainer->Add(fHistDeDx);
355 fOutputContainer->Add(fHistDeDxRegion);
356 fOutputContainer->Add(fHistDeDxSharp);
357 fOutputContainer->Add(fAntiAlpha);
358 fOutputContainer->Add(fHistAlpha);
359 fOutputContainer->Add(fHistAlphaSignal);
360 fOutputContainer->Add(fGraphAlphaSignal);
361 fOutputContainer->Add(fHistTOF3D);
362 fOutputContainer->Add(fHistTOF2D);
363 fOutputContainer->Add(fHistTOFnuclei);
d1f5b077 364 fOutputContainer->Add(fHistHelium4PtGen);
365 fOutputContainer->Add(fHistHelium4PtGenPrim);
366 fOutputContainer->Add(fHistHelium4PtGenSec);
367 fOutputContainer->Add(fHistHelium4PtGenEta);
368 fOutputContainer->Add(fHistHelium4PtGenPrimEta);
369 fOutputContainer->Add(fHistAntiHelium4PtGen);
370 fOutputContainer->Add(fHistAntiHelium4PtGenPrim);
371 fOutputContainer->Add(fHistAntiHelium4PtGenSec);
372 fOutputContainer->Add(fHistAntiHelium4PtGenEta);
373 fOutputContainer->Add(fHistHelium4PtAso);
374 fOutputContainer->Add(fHistHelium4PtAsoPrim);
375 fOutputContainer->Add(fHistHelium4PtAsoSec);
376 fOutputContainer->Add(fHistAntiHelium4PtAso);
110805a7 377
378
379 //------------ Tree and branch definitions ----------------//
439dd8c8 380 fTree = new TTree("tree", "alpha tree");
110805a7 381 //------------ Event variables ------------//
439dd8c8 382 fTree->Branch("fName",fName,"fName/C");
383 fTree->Branch("fEvnt",&fEvnt, "fEvnt/I");
384 fTree->Branch("fFileName",fFileName,"fFileName/C");
385 fTree->Branch("fEventNumber",fEventNumber,"fEventNumber/I");
386 fTree->Branch("fItrk",&fItrk, "fItrk/I");
110805a7 387 //-------------------------------------------//
388 //----------- Track variables --------------//
439dd8c8 389 fTree->Branch("fEta",fEta,"fEta[fItrk]/D");
390 fTree->Branch("fKinkIndex",fKinkIndex,"fKinkIndex[fItrk]/I");
4c074a34 391 fTree->Branch("fCentrality",fCentrality,"fCentrality[fItrk]/F");
439dd8c8 392 //
393 fTree->Branch("fTPCnCluster",fTPCnCluster,"fTPCnCluster[fItrk]/s");
394 fTree->Branch("fTPCNsignal",fTPCNsignal,"fTPCNsignal[fItrk]/s");
395 fTree->Branch("fChi2PerClusterTPC",fChi2PerClusterTPC,"fChi2PerClusterTPC[fItrk]/D");
396 fTree->Branch("fTPCRefit",fTPCRefit,"fTPCRefit[fItrk]/O");
397 fTree->Branch("fTPCsignal0",fTPCsignal0,"fTPCsignal0[fItrk]/D");
398 fTree->Branch("fTPCsignal1",fTPCsignal1,"fTPCsignal1[fItrk]/D");
399 fTree->Branch("fTPCsignal2",fTPCsignal2,"fTPCsignal2[fItrk]/D");
400 fTree->Branch("fTPCsignal3",fTPCsignal3,"fTPCsignal3[fItrk]/D");
401 fTree->Branch("fTPCSharedClusters",fTPCSharedClusters,"fTPCSharedClusters[fItrk]/I");
402 fTree->Branch("fTPCNclsIter1",fTPCNclsIter1,"fTPCNclsIter1[fItrk]/s");
403 //
404 fTree->Branch("fITSsignal",fITSsignal,"fITSsignal[fItrk]/D");
405 fTree->Branch("fITSnCluster",fITSnCluster,"fITSnCluster[fItrk]/I");
406 fTree->Branch("fChi2PerClusterITS",fChi2PerClusterITS,"fChi2PerClusterITS[fItrk]/D");
407 fTree->Branch("fITSRefit",fITSRefit,"fITSRefit[fItrk]/O");
d1f5b077 408 //
439dd8c8 409 fTree->Branch("fTOFtime",fTOFtime,"fTOFtime[fItrk]/O");
410 fTree->Branch("fTOFRefit",fTOFRefit,"fTOFRefit[fItrk]/O");
411 fTree->Branch("fTOFout",fTOFout,"fTOFout[fItrk]/O");
412 fTree->Branch("fTOFsignalDz",fTOFsignalDz,"fTOFsignalDz[fItrk]/D");
413 fTree->Branch("fTOFsignalDx",fTOFsignalDx,"fTOFsignalDx[fItrk]/D");
414 //
06ea1204 415 fTree->Branch("fTRDin",fTRDin,"fTRDin[fItrk]/O");
416 //
439dd8c8 417 fTree->Branch("fDCAXY",fDCAXY,"fDCAXY[fItrk]/F");
418 fTree->Branch("fDCAZ",fDCAZ,"fDCAZ[fItrk]/F");
419 //
420 fTree->Branch("fTrkPtot",fTrkPtot,"fTrkPtot[fItrk]/D");
421 fTree->Branch("fTPCPtot",fTPCPtot,"fTPCPtot[fItrk]/D");
422 fTree->Branch("fTrackPt",fTrackPt,"fTrackPt[fItrk]/D");
423 fTree->Branch("fDeDx",fDeDx,"fDeDx[fItrk]/D");
424 fTree->Branch("fSign",fSign,"fSign[fItrk]/D");
425 fTree->Branch("fMass",fMass,"Mass[fItrk]/F");
426 //
427 fTree->Branch("fAssociated",fAssociated,"fAssociated[fItrk]/O");
110805a7 428
429}
430
431//________________________________________________________________________
432void AliAnalysisTaskAntiHe4::UserExec(Option_t *)
433{
434 // Main loop
435 // Called for each event
436
437 //get Event-Handler for the trigger information
438 fEventHandler= dynamic_cast<AliInputEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
439 if (!fEventHandler) {
440 AliError("Could not get InputHandler");
441 //return -1;
442 return;
443 }
d1f5b077 444
445 // Monte Carlo
446 AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
447 if (!eventHandler){
448 //Printf("ERROR: Could not retrieve MC event handler");
449 fMCtrue = kFALSE;
450 }
451
452 AliMCEvent* mcEvent = 0x0;
453 AliStack* stack = 0x0;
454 if (eventHandler) mcEvent = eventHandler->MCEvent();
455 if (!mcEvent){
456 //Printf("ERROR: Could not retrieve MC event");
457 if (fMCtrue) return;
458 }
459
460 if (fMCtrue){
461 stack = mcEvent->Stack();
462 if (!stack) return;
463 }
464
465 //look for the generated particles
466 MCGenerated(stack);
110805a7 467
468 fESD=dynamic_cast<AliESDEvent*>(InputEvent());
469 if (!fESD) {
470 //Printf("ERROR: fESD not available");
471 return;
472 }
473
474 if (SetupEvent() < 0) {
475 PostData(1, fOutputContainer);
476 return;
477 }
478
479 const AliESDVertex *vertex = fESD->GetPrimaryVertexTracks();
480 if(vertex->GetNContributors()<1) {
481 // SPD vertex
482 vertex = fESD->GetPrimaryVertexSPD();
483 if(vertex->GetNContributors()<1) {
484 PostData(1, fOutputContainer);
485 return;
486 }
487
488 }
489 // check if event is selected by physics selection class
490 //
491 Bool_t isSelected = kFALSE;
492 isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
493 if (!isSelected || TMath::Abs(vertex->GetZv()) > 10) {
494 PostData(1, fOutputContainer);
495 return;
496 }
497
498 //
499 //centrality selection in PbPb
500 //
501 Int_t centralityClass10 = -1;
4c074a34 502 Float_t centralityPercentile = -1;
110805a7 503 //
504 if (fESD->GetEventSpecie() == 4) { //species == 4 == PbPb
505 // PbPb
506 AliCentrality *esdCentrality = fESD->GetCentrality();
507 centralityClass10 = esdCentrality->GetCentralityClass10("V0M"); // centrality percentile determined with V0
508 centralityPercentile = esdCentrality->GetCentralityPercentile("V0M"); // centrality percentile determined with V0
509 if(!fMCtrue){
510 if (centralityPercentile < 0. || centralityPercentile > 80. ) return; //select only events with centralities between 0 and 80 %
511 }
512 }
513 //
06ea1204 514 if (!fTriggerFired[0] && !fTriggerFired[1] && !fTriggerFired[2]) return; // select only events which pass kMB, kCentral, kSemiCentral
515 //
110805a7 516 fHistCentralityClass10->Fill(centralityClass10);
517 fHistCentralityPercentile->Fill(centralityPercentile);
518 //
519 if(fTriggerFired[0] == kTRUE) fHistTriggerStatAfterEventSelection->Fill(0);
520 if(fTriggerFired[1] == kTRUE) fHistTriggerStatAfterEventSelection->Fill(1);
521 if(fTriggerFired[2] == kTRUE) fHistTriggerStatAfterEventSelection->Fill(2);
522 if(fTriggerFired[3] == kTRUE) fHistTriggerStatAfterEventSelection->Fill(3);
523 if(fTriggerFired[4] == kTRUE) fHistTriggerStatAfterEventSelection->Fill(4);
524 //
525 //
526 if (!fESDpid) fESDpid = ((AliESDInputHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetESDpid();
527 if (!fESDpid) {
528 fESDpid = new AliESDpid(); // HACK FOR MC PBPB --> PLEASE REMOVE AS SOON AS POSSIBLE
529 fESDpid->GetTPCResponse().SetBetheBlochParameters(1.28778e+00/50., 3.13539e+01, TMath::Exp(-3.16327e+01), 1.87901e+00, 6.41583e+00);
530 }
531 //
532 Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z for th // for Anti-Alpha
439dd8c8 533 fEvnt = fESD->GetEventNumberInFile();
534 sscanf(fInputHandler->GetTree()->GetCurrentFile()->GetName(),"%s", fName);
535 fItrk = 0;
110805a7 536 //
537 Int_t runNumber = 0;
538 runNumber = fESD->GetRunNumber();
539 //
540 Bool_t fillTree = kFALSE;
541 // Track loop to fill the spectram
542 for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
439dd8c8 543
544 fEventNumber[fItrk] = -1;
545
546 fEta[fItrk] = -2;
4c074a34 547 fCentrality[fItrk] = -1;
439dd8c8 548 fTPCNsignal[fItrk] = -1;
549 fTPCnCluster[fItrk] = -1;
550 fChi2PerClusterTPC[fItrk] = -1;
551 fTPCRefit[fItrk] = kFALSE;
552 fTPCsignal0[fItrk] = -1;
553 fTPCsignal1[fItrk] = -1;
554 fTPCsignal2[fItrk] = -1;
555 fTPCsignal3[fItrk] = -1;
556 fTPCSharedClusters[fItrk] = -1;
557 fTPCNclsIter1[fItrk] = -1;
558
559 fITSsignal[fItrk] = -1;
560 fITSnCluster[fItrk] = -1;
561 fChi2PerClusterITS[fItrk] = -1;
562 fITSRefit[fItrk] = kFALSE;
563
564 fTOFRefit[fItrk] = kFALSE;
565 fTOFtime[fItrk] = kFALSE;
566 fTOFout[fItrk] = kFALSE;
567 fTOFsignalDz[fItrk] = -1;
568 fTOFsignalDx[fItrk] = -1;
569
06ea1204 570 fTRDin[fItrk] = kFALSE;
571
439dd8c8 572 fDCAZ[fItrk] = -1;
573 fDCAXY[fItrk] = -1;
574
575 fTrkPtot[fItrk] = -1;
576 fTPCPtot[fItrk] = -1;
577 fTrackPt[fItrk] = -1;
578 fDeDx[fItrk] = -1;
579 fSign[fItrk] = -2;
580 fMass[fItrk] = -1;
581
582 fAssociated[fItrk] = kFALSE;
583
110805a7 584 AliESDtrack* track = dynamic_cast<AliESDtrack*>(fESD->GetTrack(iTracks));
585 if (!fESDtrackCuts->AcceptTrack(track)) continue;
586 //
587 Double_t nClustersTPCPID=track->GetTPCsignalN();
588 if(nClustersTPCPID < 60) continue;
589 //
590 if (!track->GetInnerParam()) continue;
591 //
592 Double32_t signal[4] = {0,0,0,0};
593 Char_t ncl[3];
594 Char_t nrows[3];
595 if (track->GetTPCdEdxInfo()) {
596 track->GetTPCdEdxInfo()->GetTPCSignalRegionInfo(signal,ncl,nrows);
597 }
598 //
599 Double_t ptot = track->GetInnerParam()->GetP();
600 Double_t ptotInc = track->GetP(); // total momentum of the incoming particle
601 Double_t sign = track->GetSign();
602 //
603 Double_t eta = TMath::Abs(track->Eta());
604 TBits shared = track->GetTPCSharedMap();
605
606 track->GetImpactParameters(dca, cov);
607 Float_t dcaXY = TMath::Sqrt(dca[0]*dca[0]);
608 Float_t dcaXYsign = dca[0];
609 Float_t dcaZ = TMath::Sqrt(dca[1]*dca[1]);
610 //
611 //
612 Double_t tpcSignal = track->GetTPCsignal();
613 //
614 //PID via specific energy loss in the TPC
615 //define the arrays for the Bethe-Bloch-Parameters
616 //
617
618 SetBBParameters(runNumber);
619
620 //define expected signals for the various species
621 Double_t expSignalProton = AliExternalTrackParam::BetheBlochAleph(ptot/0.93827,fBBParametersLightParticles[0],fBBParametersLightParticles[1],fBBParametersLightParticles[2],fBBParametersLightParticles[3],fBBParametersLightParticles[4]);
622 Double_t expSignalDeuteron = AliExternalTrackParam::BetheBlochAleph(ptot/1.875612,fBBParametersLightParticles[0],fBBParametersLightParticles[1],fBBParametersLightParticles[2],fBBParametersLightParticles[3],fBBParametersLightParticles[4]);
623 Double_t expSignalTriton = AliExternalTrackParam::BetheBlochAleph(ptot/2.808921,fBBParametersLightParticles[0],fBBParametersLightParticles[1],fBBParametersLightParticles[2],fBBParametersLightParticles[3],fBBParametersLightParticles[4]);
624
625 Double_t expSignalHelium3 = 4*AliExternalTrackParam::BetheBlochAleph(2*ptot/2.80941,fBBParametersNuclei[0],fBBParametersNuclei[1],fBBParametersNuclei[2],fBBParametersNuclei[3],fBBParametersNuclei[4]);
626 Double_t expSignalHelium4 = 4*AliExternalTrackParam::BetheBlochAleph(2*ptot/3.728401,fBBParametersNuclei[0],fBBParametersNuclei[1],fBBParametersNuclei[2],fBBParametersNuclei[3],fBBParametersNuclei[4]);
627
628 //
629 Float_t mass = 0;
630 Float_t time = -1;
631 Float_t beta = 0;
632 Float_t gamma = -1;
633 Bool_t hasTOFout = kFALSE;
634 Bool_t hasTOFtime = kFALSE;
635 Float_t length = track->GetIntegratedLength();
636 UInt_t status = track->GetStatus();
d1f5b077 637 Bool_t isAssociated = kFALSE;
110805a7 638
639 if (length > 350) {
640 time = track->GetTOFsignal() - fESDpid->GetTOFResponse().GetStartTime(track->P());
641 if (time > 0) {
642 beta = length / (2.99792457999999984e-02 * time);
643 gamma = 1/TMath::Sqrt(1 - beta*beta);
644 mass = ptot/TMath::Sqrt(gamma*gamma - 1); // using inner TPC mom. as approx.
645 }
646 }
647 //
648 // fill tree and print candidates (for short list)
649 //
650 Float_t cut = 4*AliExternalTrackParam::BetheBlochAleph(2*ptot/(0.938*3),1.1931,
651 31.9806,
652 5.04114e-11,
653 2.13096,
654 2.38541);
d1f5b077 655 if (fMCtrue) cut = 4*AliExternalTrackParam::BetheBlochAleph(2*ptot/(0.938*3),0.9931,
656 31.9806,
657 5.04114e-11,
658 2.13096,
659 2.38541);
439dd8c8 660 if (eta < 1.0 && tpcSignal > 120 && tpcSignal > cut && tpcSignal < 1000 && track->GetTPCsignalN() > 60 && dcaZ < 15 && dcaXY < 15 && ptot > 1.0 && ptot < 20) {
110805a7 661 //
d1f5b077 662 cout << "AntiAlphaEvent" << " "
110805a7 663 << AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()->GetTree()->GetCurrentFile()->GetName() << " "
664 << "event number in file: " << fESD->GetEventNumberInFile()
665 << " Index " << iTracks
666 << " ptot: " << ptot
667 << " sig: " << tpcSignal <<endl;
668 //
669 fillTree = kTRUE;
670 //
439dd8c8 671
110805a7 672 sscanf(fInputHandler->GetTree()->GetCurrentFile()->GetName(),"%s", fFileName);
439dd8c8 673 fEventNumber[fItrk] = fESD->GetEventNumberInFile();
674
675 fEta[fItrk] = eta;
676 fKinkIndex[fItrk] = track->GetKinkIndex(0);
4c074a34 677 fCentrality[fItrk] = centralityPercentile;
439dd8c8 678
679 fTPCNsignal[fItrk] = track->GetTPCsignalN();
680 fTPCnCluster[fItrk] = track->GetTPCNcls();
681 fChi2PerClusterTPC[fItrk] = track->GetTPCchi2()/fTPCnCluster[fItrk];
682 if(status&AliESDtrack::kTPCrefit)
683 fTPCRefit[fItrk] = kTRUE;
684 else fTPCRefit[fItrk] = kFALSE;
685 fTPCsignal0[fItrk] = signal[0];
686 fTPCsignal1[fItrk] = signal[1];
687 fTPCsignal2[fItrk] = signal[2];
688 fTPCsignal3[fItrk] = signal[3];
689 fTPCSharedClusters[fItrk] = shared.CountBits();
690 fTPCNclsIter1[fItrk] = track->GetTPCNclsIter1();
691
692 fITSsignal[fItrk] = track->GetITSsignal();
693 fITSnCluster[fItrk] = track->GetNcls(0);
694 fChi2PerClusterITS[fItrk] = track->GetITSchi2()/fITSnCluster[fItrk];
695 if(status&AliESDtrack::kITSrefit)
696 fITSRefit[fItrk] = kTRUE;
697 else fITSRefit[fItrk] = kFALSE;
698
699
700 if(status&AliESDtrack::kITSrefit)
701 fITSRefit[fItrk] = kTRUE;
702 else fITSRefit[fItrk] = kFALSE;
110805a7 703 hasTOFout = status&AliESDtrack::kTOFout;
704 hasTOFtime = status&AliESDtrack::kTIME;
439dd8c8 705 fTOFtime[fItrk] = hasTOFtime;
706 fTOFout[fItrk] = hasTOFout;
707 fTOFsignalDz[fItrk] = track->GetTOFsignalDz();
708 fTOFsignalDx[fItrk] = track->GetTOFsignalDx();
709
06ea1204 710 fTRDin[fItrk] = status&AliESDtrack::kTRDin;
711
439dd8c8 712 fDCAZ[fItrk] = dcaXY;
713 fDCAXY[fItrk] = dcaZ;
714
715 fTrkPtot[fItrk] = track->P();
716 fTPCPtot[fItrk] = ptot;
717 fTrackPt[fItrk] = track->Pt();
718 fDeDx[fItrk] = tpcSignal;
719 fSign[fItrk] = sign;
720 fMass[fItrk] = mass;
d1f5b077 721
722 if (fMCtrue){ //associated
723
724 Int_t label = track->GetLabel();
725 TParticle *tparticle = stack->Particle(TMath::Abs(label));
726
727 Bool_t isPrimary = stack->IsPhysicalPrimary(TMath::Abs(label));
728 Bool_t isSecondary = stack->IsSecondaryFromMaterial(TMath::Abs(label));
729
730 Long_t pdgCode = tparticle->GetPdgCode();
731 Double_t pT =(track->Pt())*2;
732
733 if(pdgCode == 1000020040)
734 {
735 fHistHelium4PtAso->Fill(pT);
736 if(isPrimary) fHistHelium4PtAsoPrim->Fill(pT);
737 if(isSecondary) fHistHelium4PtAsoSec->Fill(pT);
738 isAssociated = kTRUE;
739 }
740
741 if(pdgCode == -1000020040)
742 {
743 fHistAntiHelium4PtAso->Fill(pT);
744 isAssociated = kTRUE;
745 }
746
747 }
748
439dd8c8 749 fAssociated[fItrk] = isAssociated;
750
751 fItrk++;
110805a7 752 }
753 //
754 // do pid fill histogram for raw ratios
755 //
756 // (0.) dca, (1.) sign, (2.) particle Type, (3.) p_tot
757 Int_t id = -1;
758 if (ptot > 0.2 && TMath::Abs(tpcSignal - expSignalProton)/expSignalProton < 0.2) id = 0;
759 if (ptot > 0.3 && TMath::Abs(tpcSignal - expSignalDeuteron)/expSignalDeuteron < 0.2) id = 1;
760 if (ptot > 0.7 && TMath::Abs(tpcSignal - expSignalTriton)/expSignalTriton < 0.2) id = 2;
761 if (ptot > 0.5 && (tpcSignal - expSignalHelium3)/expSignalHelium3 > -0.1 && (tpcSignal - expSignalHelium3)/expSignalHelium3 < 0.2) id = 3;
762 //
2942f542 763 Double_t vecAntiAlpha[4] = {dcaXYsign, sign, static_cast<Double_t>(id), ptotInc};
110805a7 764 if (id != -1 && tpcSignal > 120) fAntiAlpha->Fill(vecAntiAlpha);
765 //
766 // fill final histograms
767 //
768 if (!fESDtrackCutsSharp->AcceptTrack(track) || shared.CountBits() > 1 || track->GetTPCsignalN() < 80 || track->GetKinkIndex(0) != 0 || track->GetTPCNclsIter1() < 80) continue;
769 //
770 fHistDEDx->Fill(ptot,track->GetTPCsignal(), sign);
771 if (TMath::Abs(tpcSignal - expSignalHelium4)/expSignalHelium4 < 0.2 && time < 99998) fHistTOF3D->Fill(mass,1,sign);
772 //
773 // remove overlapping tracks with special dE/dx cut
774 //
775 //if (signal[1] < 6) continue;
776 //if (signal[0]/signal[1] > 1.6 || signal[2]/signal[1] > 1.6 || signal[3]/signal[1] > 1.3) continue;
777 //
778 if(sign<0) {
779 fHistDeDx->Fill(ptot, track->GetTPCsignal());
780 if (track->GetTPCsignalN() > 100 &&
439dd8c8 781 TMath::Abs(track->Eta()) < 1.0 &&
110805a7 782 signal[3]/signal[1] > 0.6 &&
783 signal[0]/signal[1] > 0.5 &&
784 signal[3]/signal[1] < 1.2 &&
785 track->GetTPCNclsIter1() > 70 &&
786 track->GetTPCnclsS() < 10) {
787 fHistDeDxSharp->Fill(ptot, track->GetTPCsignal());
788 }
789 }
790 //
791 // dE/dx for specific regions
792 //
793 for(Int_t iSig = 0; iSig < 4; iSig++) {
794 if (signal[1] > 6) fHistDeDxRegion->Fill(ptot,signal[iSig]/signal[1],iSig);
795 }
796 //
797 // alpha TOF plot
798 //
799 if((track->GetTPCsignal()-expSignalHelium4)/expSignalHelium4 > -0.15 && (track->GetTPCsignal()-expSignalHelium4)/expSignalHelium4 < 0.15) {
800 //TOF
801 hasTOFout = status&AliESDtrack::kTOFout;
802 hasTOFtime = status&AliESDtrack::kTIME;
803 Bool_t hasTOF = kFALSE;
804
805 if (hasTOFout && hasTOFtime) hasTOF = kTRUE;
806 //
807 if (length < 350.) hasTOF = kFALSE;
808
809 Float_t time0 = fESDpid->GetTOFResponse().GetStartTime(track->P());//fESDpid->GetTOFResponse().GetTimeZero();
810 //
811 if (length > 350. && hasTOF == kTRUE && ptot < 4) {
812 time = track->GetTOFsignal() - time0;
813 if (time > 0) {
814 beta = length / (2.99792457999999984e-02 * time);
815 if (beta < 0.975) {
816 gamma = 1/TMath::Sqrt(1 - beta*beta);
817 mass = ptot/TMath::Sqrt(gamma*gamma - 1); // using inner TPC mom. as approx.
818 if (TMath::Sqrt(track->GetTOFsignalDz()*track->GetTOFsignalDz() + track->GetTOFsignalDx()*track->GetTOFsignalDx()) < 5.) {
819 fHistAlpha->Fill(mass*mass);
820 if (mass*mass > 3. && mass*mass < 4.) {
821 fHistAlphaSignal->Fill(mass*mass);
822 fGraphAlphaSignal->SetPoint(fNCounter, ptot, track->GetTPCsignal());
823 fNCounter++;
824 }
825 fHistTOFnuclei->Fill(ptot,beta);
826 }
827 }
828 }
829 }
830 }
831
832 }//end loop over tracks
d1f5b077 833
110805a7 834 if (fillTree) fTree->Fill();
835
836 // Post output data.
837 PostData(1, fOutputContainer);
838 PostData(2, fTree);
839}
840
841//________________________________________________________________________
842void AliAnalysisTaskAntiHe4::Terminate(const Option_t *)
843{
844 // Draw result to the screen
845 // Called once at the end of the query
846
847 //get output data and darw 'fHistPt'
848 if (!GetOutputData(0)) return;
849
850}
851
852
853//________________________________________________________________________
854void AliAnalysisTaskAntiHe4::BinLogAxis(const THn *h, Int_t axisNumber) {
855 //
856 // Method for the correct logarithmic binning of histograms
857 //
858 TAxis *axis = h->GetAxis(axisNumber);
859 int bins = axis->GetNbins();
860
861 Double_t from = axis->GetXmin();
862 Double_t to = axis->GetXmax();
863 Double_t *newBins = new Double_t[bins + 1];
864
865 newBins[0] = from;
866 Double_t factor = pow(to/from, 1./bins);
867
868 for (int i = 1; i <= bins; i++) {
869 newBins[i] = factor * newBins[i-1];
870 }
871 axis->Set(bins, newBins);
872 delete [] newBins;
873
874}
875
876//________________________________________________________________________
877void AliAnalysisTaskAntiHe4::BinLogAxis(const TH3 *h, Int_t axisNumber) {
878 //
879 // Method for the correct logarithmic binning of histograms
880 //
881 TAxis *axis = h->GetXaxis();
882 if (axisNumber == 1) axis = h->GetYaxis();
883 if (axisNumber == 2) axis = h->GetZaxis();
884 int bins = axis->GetNbins();
885
886 Double_t from = axis->GetXmin();
887 Double_t to = axis->GetXmax();
888 Double_t *newBins = new Double_t[bins + 1];
889
890 newBins[0] = from;
891 Double_t factor = pow(to/from, 1./bins);
892
893 for (int i = 1; i <= bins; i++) {
894 newBins[i] = factor * newBins[i-1];
895 }
896 axis->Set(bins, newBins);
897 delete [] newBins;
898
899}
900//________________________________________________________________________
901void AliAnalysisTaskAntiHe4::BinLogAxis(const TH1 *h) {
902 //
903 // Method for the correct logarithmic binning of histograms
904 //
905 TAxis *axis = h->GetXaxis();
906 int bins = axis->GetNbins();
907
908 Double_t from = axis->GetXmin();
909 Double_t to = axis->GetXmax();
910 Double_t *newBins = new Double_t[bins + 1];
911
912 newBins[0] = from;
913 Double_t factor = pow(to/from, 1./bins);
914
915 for (int i = 1; i <= bins; i++) {
916 newBins[i] = factor * newBins[i-1];
917 }
918 axis->Set(bins, newBins);
919 delete [] newBins;
920
921}
922//_____________________________________________________________________________
923Int_t AliAnalysisTaskAntiHe4::Initialize() {
924
925
926 // -- Reset Event
927 // ----------------
928 ResetEvent();
929
930 return 0;
931}
932//________________________________________________________________________
933Int_t AliAnalysisTaskAntiHe4::SetupEvent() {
934 // Setup Reading of event
935
936 ResetEvent();
937 // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
938 // -- Get Trigger
939 // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
940
941 //Bool_t isTriggered = IsTriggered();
942 IsTriggered();
943 return 0;
944}
945//_____________________________________________________________________________
946void AliAnalysisTaskAntiHe4::ResetEvent() {
947 // Reset event
948 // -- Reset QA
949 for (Int_t ii = 0; ii < fNTriggers; ++ii)
950 fTriggerFired[ii] = kFALSE;
951
952 return;
953}
954//________________________________________________________________________
955Bool_t AliAnalysisTaskAntiHe4::IsTriggered() {
956 // Check if Event is triggered and fill Trigger Histogram
957
958 if ((fEventHandler->IsEventSelected() & AliVEvent::kMB)) fTriggerFired[0] = kTRUE;
959 if ((fEventHandler->IsEventSelected() & AliVEvent::kCentral)) fTriggerFired[1] = kTRUE;
960 if ((fEventHandler->IsEventSelected() & AliVEvent::kSemiCentral)) fTriggerFired[2] = kTRUE;
961 if ((fEventHandler->IsEventSelected() & AliVEvent::kEMCEJE)) fTriggerFired[3] = kTRUE;
962 if ((fEventHandler->IsEventSelected() & AliVEvent::kEMCEGA)) fTriggerFired[4] = kTRUE;
963
964 Bool_t isTriggered = kFALSE;
965
966 for (Int_t ii=0; ii<fNTriggers; ++ii) {
967 if(fTriggerFired[ii]) {
968 isTriggered = kTRUE;
969 fHistTriggerStat->Fill(ii);
970 }
971 }
972
973 return isTriggered;
974}
975//________________________________________________________________________
976void AliAnalysisTaskAntiHe4::SetBBParameters(Int_t runNumber){
977
978 if(runNumber < 166500){ //LHC10h
979
980 fBBParametersLightParticles[0] = 1.45802;
981 fBBParametersLightParticles[1] = 27.4992;
982 fBBParametersLightParticles[2] = 4.00313e-15;
983 fBBParametersLightParticles[3] = 2.48485;
984 fBBParametersLightParticles[4] = 8.31768;
985
986 fBBParametersNuclei[0] = 1.69155;
987 fBBParametersNuclei[1] = 27.4992;
988 fBBParametersNuclei[2] = 4.00313e-15;
989 fBBParametersNuclei[3] = 2.48485;
990 fBBParametersNuclei[4] = 8.31768;
991
992 }
993
994 if(runNumber > 166500){ //LHC11h
995
996 fBBParametersLightParticles[0] = 4.69637;//1.11243;
997 fBBParametersLightParticles[1] = 7.51827;//26.1144;
998 fBBParametersLightParticles[2] = 0.0183746;//4.00313e-15;
999 fBBParametersLightParticles[3] = 2.60;//2.72969;
1000 fBBParametersLightParticles[4] = 2.7;//9.15038;
1001
1002 fBBParametersNuclei[0] = 1.4906;
1003 fBBParametersNuclei[1] = 27.9758;
1004 fBBParametersNuclei[2] = 4.00313e-15;
1005 fBBParametersNuclei[3] = 2.50804;
1006 fBBParametersNuclei[4] = 8.31768;
1007
1008 }
1009}
d1f5b077 1010//______________________________________________________________________________
1011void AliAnalysisTaskAntiHe4::MCGenerated(AliStack* stack)
1012{
1013
1014 // Monte Carlo for genenerated particles
1015 if (fMCtrue) //MC loop
1016 {
1017
1018 Int_t stackN = 0;
1019
1020 for(stackN = 0; stackN < stack->GetNtrack(); stackN++) //loop over stack
1021 {
1022
1023 Bool_t isPrimary = stack->IsPhysicalPrimary(stackN);
1024 Bool_t isSecondary = stack->IsSecondaryFromMaterial(stackN);
1025
1026 const TParticle *tparticle = stack->Particle(stackN);
1027 Long_t pdgCode = tparticle->GetPdgCode();
1028
1029 Double_t pTGen = tparticle->Pt();
1030 Double_t eta = tparticle->Eta();
1031 //check which particle it is
1032
1033 //Alpha
1034 if(pdgCode == 1000020040)
1035 {
1036 fHistHelium4PtGen->Fill(pTGen);
1037 if(isPrimary) fHistHelium4PtGenPrim->Fill(pTGen);
1038 if(isSecondary) fHistHelium4PtGenSec->Fill(pTGen);
439dd8c8 1039 if(TMath::Abs(eta) < 1.0)fHistHelium4PtGenEta->Fill(pTGen);
1040 if(isPrimary && TMath::Abs(eta) < 1.0)fHistHelium4PtGenPrimEta->Fill(pTGen);
d1f5b077 1041 }
1042
1043 //Anti-Alpha
1044 if(pdgCode == -1000020040)
1045 {
1046 fHistAntiHelium4PtGen->Fill(pTGen);
1047 if(isPrimary) fHistAntiHelium4PtGenPrim->Fill(pTGen);
1048 if(isSecondary) fHistAntiHelium4PtGenSec->Fill(pTGen);
439dd8c8 1049 if(TMath::Abs(eta) < 1.0)fHistAntiHelium4PtGenEta->Fill(pTGen);
d1f5b077 1050 }
1051
1052
1053 }//end loop over stack
1054
1055
1056 }//end MC
1057}