10 #include "AliESDtrackCuts.h"
11 #include "AliESDpid.h"
12 #include "AliESDVertex.h"
15 #include "AliAnalysisTask.h"
16 #include "AliAnalysisManager.h"
17 #include "AliCentrality.h"
18 #include "AliESDEvent.h"
19 #include "AliESDInputHandler.h"
21 #include "AliMCEventHandler.h"
22 #include "AliMCEvent.h"
25 #include "AliAnalysisTaskAntiHe4.h"
30 ///////////////////////////////////////////////////////////
32 // Analysis for the observation of anti-alpha particles. //
34 ///////////////////////////////////////////////////////////
37 ClassImp(AliAnalysisTaskAntiHe4)
39 //________________________________________________________________________
40 AliAnalysisTaskAntiHe4::AliAnalysisTaskAntiHe4()
41 : AliAnalysisTaskSE(),
44 fHistCentralityClass10(0),
45 fHistCentralityPercentile(0),
47 fHistTriggerStatAfterEventSelection(0),
60 fBBParametersLightParticles(),
61 fBBParametersNuclei(),
65 fESDtrackCutsSharp(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),
78 fHistHelium4PtAsoPrim(0),
79 fHistHelium4PtAsoSec(0),
80 fHistAntiHelium4PtAso(0),
86 // default Constructor
89 // Define input and output slots here
92 //________________________________________________________________________
93 AliAnalysisTaskAntiHe4::AliAnalysisTaskAntiHe4(const char *name)
94 : AliAnalysisTaskSE(name),
97 fHistCentralityClass10(0),
98 fHistCentralityPercentile(0),
100 fHistTriggerStatAfterEventSelection(0),
105 fGraphAlphaSignal(0),
113 fBBParametersLightParticles(),
114 fBBParametersNuclei(),
118 fESDtrackCutsSharp(0),
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),
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());
148 // cuts for candidates
150 fESDtrackCuts = new AliESDtrackCuts("AliESDtrackCuts","AliESDtrackCuts");
152 fESDtrackCuts->SetAcceptKinkDaughters(kFALSE);
153 fESDtrackCuts->SetMinNClustersTPC(70);
154 fESDtrackCuts->SetMaxChi2PerClusterTPC(6);
155 fESDtrackCuts->SetMaxDCAToVertexXY(3);
156 fESDtrackCuts->SetMaxDCAToVertexZ(2);
157 fESDtrackCuts->SetRequireTPCRefit(kTRUE);
158 //fESDtrackCuts->SetRequireITSRefit(kTRUE);
159 fESDtrackCuts->SetMinNClustersITS(1);
160 fESDtrackCuts->SetEtaRange(-1.0,1.0);
162 // cuts for final plots
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);
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);
197 //________________________________________________________________________
198 void AliAnalysisTaskAntiHe4::UserCreateOutputObjects()
203 //histogram to count number of events
205 fHistCentralityClass10 = new TH1F("fHistCentralityClass10", "centrality with class10", 11, -0.5, 10.5);
206 fHistCentralityClass10->GetXaxis()->SetTitle("Centrality");
207 fHistCentralityClass10->GetYaxis()->SetTitle("Entries");
209 fHistCentralityPercentile = new TH1F("fHistCentralityPercentile", "centrality with percentile", 101, -0.1, 100.1);
210 fHistCentralityPercentile->GetXaxis()->SetTitle("Centrality");
211 fHistCentralityPercentile->GetYaxis()->SetTitle("Entries");
213 //trigger statitics histogram
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]);
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]);
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");
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);
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)");
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);
246 fHistTOF3D = new TH3F("fHistTOF3D","mass determination from TOF",400,0.25,4.5,2,-0.5,1.5,2,-1,1);
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)");
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)");
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})");
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})");
265 fGraphAlphaSignal = new TGraph(20);
268 // (0.) dca, (1.) sign, (2.) particle Type, (3.) p_tot
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};
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};
280 fAntiAlpha = new THnF("fAntiAlpha", "AntiAlpha; (0.) dca, (1.) sign, (2.) particle Type, (3.) p_tot", 4, binsAntiAlpha, xminAntiAlpha, xmaxAntiAlpha);
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]);
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)");
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)");
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)");
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)");
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)");
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)");
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)");
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)");
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)");
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)");
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)");
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)");
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)");
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);
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);
379 //------------ Tree and branch definitions ----------------//
380 fTree = new TTree("tree", "alpha tree");
381 //------------ Event variables ------------//
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");
387 //-------------------------------------------//
388 //----------- Track variables --------------//
389 fTree->Branch("fEta",fEta,"fEta[fItrk]/D");
390 fTree->Branch("fKinkIndex",fKinkIndex,"fKinkIndex[fItrk]/I");
391 fTree->Branch("fCentrality",fCentrality,"fCentrality[fItrk]/F");
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");
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");
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");
415 fTree->Branch("fTRDin",fTRDin,"fTRDin[fItrk]/O");
417 fTree->Branch("fDCAXY",fDCAXY,"fDCAXY[fItrk]/F");
418 fTree->Branch("fDCAZ",fDCAZ,"fDCAZ[fItrk]/F");
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,"fMass[fItrk]/F");
426 fTree->Branch("fTime",fTime,"fTime[fItrk]/F");
427 fTree->Branch("fLength",fLength,"fLength[fItrk]/F");
428 fTree->Branch("fSigmaQP",fSigmaQP,"fSigmaQP[fItrk]/D");
430 fTree->Branch("fAssociated",fAssociated,"fAssociated[fItrk]/O");
434 //________________________________________________________________________
435 void AliAnalysisTaskAntiHe4::UserExec(Option_t *)
438 // Called for each event
440 //get Event-Handler for the trigger information
441 fEventHandler= dynamic_cast<AliInputEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
442 if (!fEventHandler) {
443 AliError("Could not get InputHandler");
449 AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
451 //Printf("ERROR: Could not retrieve MC event handler");
455 AliMCEvent* mcEvent = 0x0;
456 AliStack* stack = 0x0;
457 if (eventHandler) mcEvent = eventHandler->MCEvent();
459 //Printf("ERROR: Could not retrieve MC event");
464 stack = mcEvent->Stack();
468 //look for the generated particles
471 fESD=dynamic_cast<AliESDEvent*>(InputEvent());
473 //Printf("ERROR: fESD not available");
477 if (SetupEvent() < 0) {
478 PostData(1, fOutputContainer);
482 const AliESDVertex *vertex = fESD->GetPrimaryVertexTracks();
483 if(vertex->GetNContributors()<1) {
485 vertex = fESD->GetPrimaryVertexSPD();
486 if(vertex->GetNContributors()<1) {
487 PostData(1, fOutputContainer);
492 // check if event is selected by physics selection class
494 Bool_t isSelected = kFALSE;
495 isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
496 if (!isSelected || TMath::Abs(vertex->GetZv()) > 10) {
497 PostData(1, fOutputContainer);
502 //centrality selection in PbPb
504 Int_t centralityClass10 = -1;
505 Float_t centralityPercentile = -1;
507 if (fESD->GetEventSpecie() == 4) { //species == 4 == PbPb
509 AliCentrality *esdCentrality = fESD->GetCentrality();
510 centralityClass10 = esdCentrality->GetCentralityClass10("V0M"); // centrality percentile determined with V0
511 centralityPercentile = esdCentrality->GetCentralityPercentile("V0M"); // centrality percentile determined with V0
513 if (centralityPercentile < 0. || centralityPercentile > 80. ) return; //select only events with centralities between 0 and 80 %
517 if (!fTriggerFired[0] && !fTriggerFired[1] && !fTriggerFired[2]) return; // select only events which pass kMB, kCentral, kSemiCentral
519 fHistCentralityClass10->Fill(centralityClass10);
520 fHistCentralityPercentile->Fill(centralityPercentile);
522 if(fTriggerFired[0] == kTRUE) fHistTriggerStatAfterEventSelection->Fill(0);
523 if(fTriggerFired[1] == kTRUE) fHistTriggerStatAfterEventSelection->Fill(1);
524 if(fTriggerFired[2] == kTRUE) fHistTriggerStatAfterEventSelection->Fill(2);
525 if(fTriggerFired[3] == kTRUE) fHistTriggerStatAfterEventSelection->Fill(3);
526 if(fTriggerFired[4] == kTRUE) fHistTriggerStatAfterEventSelection->Fill(4);
529 if (!fESDpid) fESDpid = ((AliESDInputHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetESDpid();
531 fESDpid = new AliESDpid(); // HACK FOR MC PBPB --> PLEASE REMOVE AS SOON AS POSSIBLE
532 fESDpid->GetTPCResponse().SetBetheBlochParameters(1.28778e+00/50., 3.13539e+01, TMath::Exp(-3.16327e+01), 1.87901e+00, 6.41583e+00);
535 Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z for th // for Anti-Alpha
536 fEvnt = fESD->GetEventNumberInFile();
537 sscanf(fInputHandler->GetTree()->GetCurrentFile()->GetName(),"%s", fName);
541 runNumber = fESD->GetRunNumber();
543 Bool_t fillTree = kFALSE;
544 // Track loop to fill the spectram
545 for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
547 fEventNumber[fItrk] = -1;
550 fCentrality[fItrk] = -1;
551 fTPCNsignal[fItrk] = -1;
552 fTPCnCluster[fItrk] = -1;
553 fChi2PerClusterTPC[fItrk] = -1;
554 fTPCRefit[fItrk] = kFALSE;
555 fTPCsignal0[fItrk] = -1;
556 fTPCsignal1[fItrk] = -1;
557 fTPCsignal2[fItrk] = -1;
558 fTPCsignal3[fItrk] = -1;
559 fTPCSharedClusters[fItrk] = -1;
560 fTPCNclsIter1[fItrk] = -1;
562 fITSsignal[fItrk] = -1;
563 fITSnCluster[fItrk] = -1;
564 fChi2PerClusterITS[fItrk] = -1;
565 fITSRefit[fItrk] = kFALSE;
567 fTOFRefit[fItrk] = kFALSE;
568 fTOFtime[fItrk] = kFALSE;
569 fTOFout[fItrk] = kFALSE;
570 fTOFsignalDz[fItrk] = -1;
571 fTOFsignalDx[fItrk] = -1;
573 fTRDin[fItrk] = kFALSE;
578 fTrkPtot[fItrk] = -1;
579 fTPCPtot[fItrk] = -1;
580 fTrackPt[fItrk] = -1;
586 fSigmaQP[fItrk] = -1;
588 fAssociated[fItrk] = kFALSE;
590 AliESDtrack* track = dynamic_cast<AliESDtrack*>(fESD->GetTrack(iTracks));
591 if (!fESDtrackCuts->AcceptTrack(track)) continue;
593 Double_t nClustersTPCPID=track->GetTPCsignalN();
594 if(nClustersTPCPID < 60) continue;
596 if (!track->GetInnerParam()) continue;
598 Double32_t signal[4] = {0,0,0,0};
601 if (track->GetTPCdEdxInfo()) {
602 track->GetTPCdEdxInfo()->GetTPCSignalRegionInfo(signal,ncl,nrows);
605 Double_t ptot = track->GetInnerParam()->GetP();
606 Double_t ptotInc = track->GetP(); // total momentum of the incoming particle
607 Double_t sign = track->GetSign();
609 Double_t eta = TMath::Abs(track->Eta());
610 TBits shared = track->GetTPCSharedMap();
612 track->GetImpactParameters(dca, cov);
613 Float_t dcaXY = TMath::Sqrt(dca[0]*dca[0]);
614 Float_t dcaXYsign = dca[0];
615 Float_t dcaZ = TMath::Sqrt(dca[1]*dca[1]);
618 track->GetExternalCovariance(cov1);//->GetSigmaQoverP();
620 Double_t tpcSignal = track->GetTPCsignal();
622 //PID via specific energy loss in the TPC
623 //define the arrays for the Bethe-Bloch-Parameters
626 SetBBParameters(runNumber);
628 //define expected signals for the various species
629 Double_t expSignalProton = AliExternalTrackParam::BetheBlochAleph(ptot/0.93827,fBBParametersLightParticles[0],fBBParametersLightParticles[1],fBBParametersLightParticles[2],fBBParametersLightParticles[3],fBBParametersLightParticles[4]);
630 Double_t expSignalDeuteron = AliExternalTrackParam::BetheBlochAleph(ptot/1.875612,fBBParametersLightParticles[0],fBBParametersLightParticles[1],fBBParametersLightParticles[2],fBBParametersLightParticles[3],fBBParametersLightParticles[4]);
631 Double_t expSignalTriton = AliExternalTrackParam::BetheBlochAleph(ptot/2.808921,fBBParametersLightParticles[0],fBBParametersLightParticles[1],fBBParametersLightParticles[2],fBBParametersLightParticles[3],fBBParametersLightParticles[4]);
633 Double_t expSignalHelium3 = 4*AliExternalTrackParam::BetheBlochAleph(2*ptot/2.80941,fBBParametersNuclei[0],fBBParametersNuclei[1],fBBParametersNuclei[2],fBBParametersNuclei[3],fBBParametersNuclei[4]);
634 Double_t expSignalHelium4 = 4*AliExternalTrackParam::BetheBlochAleph(2*ptot/3.728401,fBBParametersNuclei[0],fBBParametersNuclei[1],fBBParametersNuclei[2],fBBParametersNuclei[3],fBBParametersNuclei[4]);
641 Bool_t hasTOFout = kFALSE;
642 Bool_t hasTOFtime = kFALSE;
643 Float_t length = track->GetIntegratedLength();
644 UInt_t status = track->GetStatus();
645 Bool_t isAssociated = kFALSE;
648 time = track->GetTOFsignal() - fESDpid->GetTOFResponse().GetStartTime(track->P());
650 beta = length / (2.99792457999999984e-02 * time);
651 gamma = 1/TMath::Sqrt(1 - beta*beta);
652 mass = ptot/TMath::Sqrt(gamma*gamma - 1); // using inner TPC mom. as approx.
656 // fill tree and print candidates (for short list)
658 Float_t cut = 4*AliExternalTrackParam::BetheBlochAleph(2*ptot/(0.938*3),1.1931,
663 if (fMCtrue) cut = 4*AliExternalTrackParam::BetheBlochAleph(2*ptot/(0.938*3),0.9931,
668 Bool_t IsDeuteron = kFALSE;
669 Bool_t IsTriton = kFALSE;
671 Double_t DeuteronSigma = TMath::Abs(tpcSignal - expSignalDeuteron)/expSignalDeuteron;
672 Double_t TritonSigma = TMath::Abs(tpcSignal - expSignalTriton)/expSignalTriton;
675 if(DeuteronSigma < 0.3 && runNumber < 166500) IsDeuteron = kTRUE;
676 if(TritonSigma < 0.3 && runNumber < 166500) IsTriton = kTRUE;
678 //cout << "isDeuteron: " << IsDeuteron << endl;
680 //if(tpcSignal > cut || IsDeuteron == kTRUE || IsTriton == kTRUE){
681 if(tpcSignal > cut ){
683 //cout << "isDeuteron: " << IsDeuteron << endl;
684 if (eta < 1.0 && tpcSignal < 1000 && dcaZ < 15 && dcaXY < 15 && ptot < 20 && ptot > 1.0 && tpcSignal > 120){// && track->GetTPCsignalN() > 60) { // && ptot > 1.0 &6 tpcSignal > 120
686 cout << "AntiAlphaEvent" << " "
687 << AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()->GetTree()->GetCurrentFile()->GetName() << " "
688 << "event number in file: " << fESD->GetEventNumberInFile()
689 << " Index " << iTracks
691 << " sig: " << tpcSignal <<endl;
696 sscanf(fInputHandler->GetTree()->GetCurrentFile()->GetName(),"%s", fFileName);
697 fEventNumber[fItrk] = fESD->GetEventNumberInFile();
700 fKinkIndex[fItrk] = track->GetKinkIndex(0);
701 fCentrality[fItrk] = centralityPercentile;
703 fTPCNsignal[fItrk] = track->GetTPCsignalN();
704 fTPCnCluster[fItrk] = track->GetTPCNcls();
705 fChi2PerClusterTPC[fItrk] = track->GetTPCchi2()/fTPCnCluster[fItrk];
706 if(status&AliESDtrack::kTPCrefit)
707 fTPCRefit[fItrk] = kTRUE;
708 else fTPCRefit[fItrk] = kFALSE;
709 fTPCsignal0[fItrk] = signal[0];
710 fTPCsignal1[fItrk] = signal[1];
711 fTPCsignal2[fItrk] = signal[2];
712 fTPCsignal3[fItrk] = signal[3];
713 fTPCSharedClusters[fItrk] = shared.CountBits();
714 fTPCNclsIter1[fItrk] = track->GetTPCNclsIter1();
716 fITSsignal[fItrk] = track->GetITSsignal();
717 fITSnCluster[fItrk] = track->GetNcls(0);
718 fChi2PerClusterITS[fItrk] = track->GetITSchi2()/fITSnCluster[fItrk];
719 if(status&AliESDtrack::kITSrefit)
720 fITSRefit[fItrk] = kTRUE;
721 else fITSRefit[fItrk] = kFALSE;
724 if(status&AliESDtrack::kITSrefit)
725 fITSRefit[fItrk] = kTRUE;
726 else fITSRefit[fItrk] = kFALSE;
727 hasTOFout = status&AliESDtrack::kTOFout;
728 hasTOFtime = status&AliESDtrack::kTIME;
729 fTOFtime[fItrk] = hasTOFtime;
730 fTOFout[fItrk] = hasTOFout;
731 fTOFsignalDz[fItrk] = track->GetTOFsignalDz();
732 fTOFsignalDx[fItrk] = track->GetTOFsignalDx();
734 fTRDin[fItrk] = status&AliESDtrack::kTRDin;
736 fDCAZ[fItrk] = dcaXY;
737 fDCAXY[fItrk] = dcaZ;
739 fTrkPtot[fItrk] = track->P();
740 fTPCPtot[fItrk] = ptot;
741 fTrackPt[fItrk] = track->Pt();
742 fDeDx[fItrk] = tpcSignal;
746 fLength[fItrk] = length;
747 fSigmaQP[fItrk] = cov1[15];
749 if (fMCtrue){ //associated
751 Int_t label = track->GetLabel();
752 TParticle *tparticle = stack->Particle(TMath::Abs(label));
754 Bool_t isPrimary = stack->IsPhysicalPrimary(TMath::Abs(label));
755 Bool_t isSecondary = stack->IsSecondaryFromMaterial(TMath::Abs(label));
757 Long_t pdgCode = tparticle->GetPdgCode();
758 Double_t pT =(track->Pt())*2;
760 if(pdgCode == 1000020040)
762 fHistHelium4PtAso->Fill(pT);
763 if(isPrimary) fHistHelium4PtAsoPrim->Fill(pT);
764 if(isSecondary) fHistHelium4PtAsoSec->Fill(pT);
765 isAssociated = kTRUE;
768 if(pdgCode == -1000020040)
770 fHistAntiHelium4PtAso->Fill(pT);
771 isAssociated = kTRUE;
776 fAssociated[fItrk] = isAssociated;
782 // do pid fill histogram for raw ratios
784 // (0.) dca, (1.) sign, (2.) particle Type, (3.) p_tot
786 if (ptot > 0.2 && TMath::Abs(tpcSignal - expSignalProton)/expSignalProton < 0.2) id = 0;
787 if (ptot > 0.3 && TMath::Abs(tpcSignal - expSignalDeuteron)/expSignalDeuteron < 0.2) id = 1;
788 if (ptot > 0.7 && TMath::Abs(tpcSignal - expSignalTriton)/expSignalTriton < 0.2) id = 2;
789 if (ptot > 0.5 && (tpcSignal - expSignalHelium3)/expSignalHelium3 > -0.1 && (tpcSignal - expSignalHelium3)/expSignalHelium3 < 0.2) id = 3;
791 Double_t vecAntiAlpha[4] = {dcaXYsign, sign, static_cast<Double_t>(id), ptotInc};
792 if (id != -1 && tpcSignal > 120) fAntiAlpha->Fill(vecAntiAlpha);
794 // fill final histograms
796 if (!fESDtrackCutsSharp->AcceptTrack(track) || shared.CountBits() > 1 || track->GetTPCsignalN() < 80 || track->GetKinkIndex(0) != 0 || track->GetTPCNclsIter1() < 80) continue;
798 fHistDEDx->Fill(ptot,track->GetTPCsignal(), sign);
799 if (TMath::Abs(tpcSignal - expSignalHelium4)/expSignalHelium4 < 0.2 && time < 99998) fHistTOF3D->Fill(mass,1,sign);
801 // remove overlapping tracks with special dE/dx cut
803 //if (signal[1] < 6) continue;
804 //if (signal[0]/signal[1] > 1.6 || signal[2]/signal[1] > 1.6 || signal[3]/signal[1] > 1.3) continue;
807 fHistDeDx->Fill(ptot, track->GetTPCsignal());
808 if (track->GetTPCsignalN() > 100 &&
809 TMath::Abs(track->Eta()) < 1.0 &&
810 signal[3]/signal[1] > 0.6 &&
811 signal[0]/signal[1] > 0.5 &&
812 signal[3]/signal[1] < 1.2 &&
813 track->GetTPCNclsIter1() > 70 &&
814 track->GetTPCnclsS() < 10) {
815 fHistDeDxSharp->Fill(ptot, track->GetTPCsignal());
819 // dE/dx for specific regions
821 for(Int_t iSig = 0; iSig < 4; iSig++) {
822 if (signal[1] > 6) fHistDeDxRegion->Fill(ptot,signal[iSig]/signal[1],iSig);
827 if((track->GetTPCsignal()-expSignalHelium4)/expSignalHelium4 > -0.15 && (track->GetTPCsignal()-expSignalHelium4)/expSignalHelium4 < 0.15) {
829 hasTOFout = status&AliESDtrack::kTOFout;
830 hasTOFtime = status&AliESDtrack::kTIME;
831 Bool_t hasTOF = kFALSE;
833 if (hasTOFout && hasTOFtime) hasTOF = kTRUE;
835 if (length < 350.) hasTOF = kFALSE;
837 Float_t time0 = fESDpid->GetTOFResponse().GetStartTime(track->P());//fESDpid->GetTOFResponse().GetTimeZero();
839 if (length > 350. && hasTOF == kTRUE && ptot < 4) {
840 time = track->GetTOFsignal() - time0;
842 beta = length / (2.99792457999999984e-02 * time);
844 gamma = 1/TMath::Sqrt(1 - beta*beta);
845 mass = ptot/TMath::Sqrt(gamma*gamma - 1); // using inner TPC mom. as approx.
846 if (TMath::Sqrt(track->GetTOFsignalDz()*track->GetTOFsignalDz() + track->GetTOFsignalDx()*track->GetTOFsignalDx()) < 5.) {
847 fHistAlpha->Fill(mass*mass);
848 if (mass*mass > 3. && mass*mass < 4.) {
849 fHistAlphaSignal->Fill(mass*mass);
850 fGraphAlphaSignal->SetPoint(fNCounter, ptot, track->GetTPCsignal());
853 fHistTOFnuclei->Fill(ptot,beta);
860 }//end loop over tracks
862 if (fillTree) fTree->Fill();
865 PostData(1, fOutputContainer);
869 //________________________________________________________________________
870 void AliAnalysisTaskAntiHe4::Terminate(const Option_t *)
872 // Draw result to the screen
873 // Called once at the end of the query
875 //get output data and darw 'fHistPt'
876 if (!GetOutputData(0)) return;
881 //________________________________________________________________________
882 void AliAnalysisTaskAntiHe4::BinLogAxis(const THn *h, Int_t axisNumber) {
884 // Method for the correct logarithmic binning of histograms
886 TAxis *axis = h->GetAxis(axisNumber);
887 int bins = axis->GetNbins();
889 Double_t from = axis->GetXmin();
890 Double_t to = axis->GetXmax();
891 Double_t *newBins = new Double_t[bins + 1];
894 Double_t factor = pow(to/from, 1./bins);
896 for (int i = 1; i <= bins; i++) {
897 newBins[i] = factor * newBins[i-1];
899 axis->Set(bins, newBins);
904 //________________________________________________________________________
905 void AliAnalysisTaskAntiHe4::BinLogAxis(const TH3 *h, Int_t axisNumber) {
907 // Method for the correct logarithmic binning of histograms
909 TAxis *axis = h->GetXaxis();
910 if (axisNumber == 1) axis = h->GetYaxis();
911 if (axisNumber == 2) axis = h->GetZaxis();
912 int bins = axis->GetNbins();
914 Double_t from = axis->GetXmin();
915 Double_t to = axis->GetXmax();
916 Double_t *newBins = new Double_t[bins + 1];
919 Double_t factor = pow(to/from, 1./bins);
921 for (int i = 1; i <= bins; i++) {
922 newBins[i] = factor * newBins[i-1];
924 axis->Set(bins, newBins);
928 //________________________________________________________________________
929 void AliAnalysisTaskAntiHe4::BinLogAxis(const TH1 *h) {
931 // Method for the correct logarithmic binning of histograms
933 TAxis *axis = h->GetXaxis();
934 int bins = axis->GetNbins();
936 Double_t from = axis->GetXmin();
937 Double_t to = axis->GetXmax();
938 Double_t *newBins = new Double_t[bins + 1];
941 Double_t factor = pow(to/from, 1./bins);
943 for (int i = 1; i <= bins; i++) {
944 newBins[i] = factor * newBins[i-1];
946 axis->Set(bins, newBins);
950 //_____________________________________________________________________________
951 Int_t AliAnalysisTaskAntiHe4::Initialize() {
960 //________________________________________________________________________
961 Int_t AliAnalysisTaskAntiHe4::SetupEvent() {
962 // Setup Reading of event
965 // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
967 // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
969 //Bool_t isTriggered = IsTriggered();
973 //_____________________________________________________________________________
974 void AliAnalysisTaskAntiHe4::ResetEvent() {
977 for (Int_t ii = 0; ii < fNTriggers; ++ii)
978 fTriggerFired[ii] = kFALSE;
982 //________________________________________________________________________
983 Bool_t AliAnalysisTaskAntiHe4::IsTriggered() {
984 // Check if Event is triggered and fill Trigger Histogram
986 if ((fEventHandler->IsEventSelected() & AliVEvent::kMB)) fTriggerFired[0] = kTRUE;
987 if ((fEventHandler->IsEventSelected() & AliVEvent::kCentral)) fTriggerFired[1] = kTRUE;
988 if ((fEventHandler->IsEventSelected() & AliVEvent::kSemiCentral)) fTriggerFired[2] = kTRUE;
989 if ((fEventHandler->IsEventSelected() & AliVEvent::kEMCEJE)) fTriggerFired[3] = kTRUE;
990 if ((fEventHandler->IsEventSelected() & AliVEvent::kEMCEGA)) fTriggerFired[4] = kTRUE;
992 Bool_t isTriggered = kFALSE;
994 for (Int_t ii=0; ii<fNTriggers; ++ii) {
995 if(fTriggerFired[ii]) {
997 fHistTriggerStat->Fill(ii);
1003 //________________________________________________________________________
1004 void AliAnalysisTaskAntiHe4::SetBBParameters(Int_t runNumber){
1006 if(runNumber < 166500){ //LHC10h
1008 fBBParametersLightParticles[0] = 1.45802;
1009 fBBParametersLightParticles[1] = 27.4992;
1010 fBBParametersLightParticles[2] = 4.00313e-15;
1011 fBBParametersLightParticles[3] = 2.48485;
1012 fBBParametersLightParticles[4] = 8.31768;
1014 fBBParametersNuclei[0] = 1.69155;
1015 fBBParametersNuclei[1] = 27.4992;
1016 fBBParametersNuclei[2] = 4.00313e-15;
1017 fBBParametersNuclei[3] = 2.48485;
1018 fBBParametersNuclei[4] = 8.31768;
1022 if(runNumber > 166500){ //LHC11h
1024 fBBParametersLightParticles[0] = 4.69637;//1.11243;
1025 fBBParametersLightParticles[1] = 7.51827;//26.1144;
1026 fBBParametersLightParticles[2] = 0.0183746;//4.00313e-15;
1027 fBBParametersLightParticles[3] = 2.60;//2.72969;
1028 fBBParametersLightParticles[4] = 2.7;//9.15038;
1030 fBBParametersNuclei[0] = 1.4906;
1031 fBBParametersNuclei[1] = 27.9758;
1032 fBBParametersNuclei[2] = 4.00313e-15;
1033 fBBParametersNuclei[3] = 2.50804;
1034 fBBParametersNuclei[4] = 8.31768;
1038 //______________________________________________________________________________
1039 void AliAnalysisTaskAntiHe4::MCGenerated(AliStack* stack)
1042 // Monte Carlo for genenerated particles
1043 if (fMCtrue) //MC loop
1048 for(stackN = 0; stackN < stack->GetNtrack(); stackN++) //loop over stack
1051 Bool_t isPrimary = stack->IsPhysicalPrimary(stackN);
1052 Bool_t isSecondary = stack->IsSecondaryFromMaterial(stackN);
1054 const TParticle *tparticle = stack->Particle(stackN);
1055 Long_t pdgCode = tparticle->GetPdgCode();
1057 Double_t pTGen = tparticle->Pt();
1058 Double_t eta = tparticle->Eta();
1059 //check which particle it is
1062 if(pdgCode == 1000020040)
1064 fHistHelium4PtGen->Fill(pTGen);
1065 if(isPrimary) fHistHelium4PtGenPrim->Fill(pTGen);
1066 if(isSecondary) fHistHelium4PtGenSec->Fill(pTGen);
1067 if(TMath::Abs(eta) < 0.8)fHistHelium4PtGenEta->Fill(pTGen);
1068 if(isPrimary && TMath::Abs(eta) < 0.8)fHistHelium4PtGenPrimEta->Fill(pTGen);
1072 if(pdgCode == -1000020040)
1074 fHistAntiHelium4PtGen->Fill(pTGen);
1075 if(isPrimary) fHistAntiHelium4PtGenPrim->Fill(pTGen);
1076 if(isSecondary) fHistAntiHelium4PtGenSec->Fill(pTGen);
1077 if(TMath::Abs(eta) < 0.8)fHistAntiHelium4PtGenEta->Fill(pTGen);
1081 }//end loop over stack