e2d2636c 1 /**************************************************************************
2 * Author : Nicole Alice Martin (nicole.alice.martin@cern.ch) *
3 * *
4 * Contributors are mentioned in the code where appropriate. *
5 * *
6 * Permission to use, copy, modify and distribute this software and its *
7 * documentation strictly for non-commercial purposes is hereby granted *
8 * without fee, provided that the above copyright notice appears in all *
9 * copies and that both the copyright notice and this permission notice *
10 * appear in the supporting documentation. The authors make no claims *
11 * about the suitability of this software for any purpose. It is *
12 * provided "as is" without express or implied warranty. *
13 **************************************************************************/
14
15 //-----------------------------------------------------------------
16 // AliAnalysisTaskLambdaNAOD class
17 // task for the investigation of (anti-)lambda-n bound state
18 // uses the V0 finders, based on AODs or ESDS
19 //-----------------------------------------------------------------
20
21
22 class TTree;
23 class TParticle;
24 class TVector3;
25
26 class AliESDVertex;
27 class AliESDv0;
28
29 class AliAODVertex;
30 class AliAODv0;
31
32 #include "TChain.h"
33 #include <Riostream.h>
34 #include <iostream>
35 #include <fstream>
36 #include <string>
37 #include "TList.h"
38 #include "TH1.h"
39 #include "TH2.h"
40 #include "TH3.h"
41 #include "THnSparse.h"
42 #include "TCanvas.h"
43 #include "TMath.h"
44 #include "TLegend.h"
45 #include "TFile.h"
46 #include <TRandom3.h>
47 #include "TPDGCode.h"
48 #include <TDatabasePDG.h>
49
50 #include "AliAnalysisManager.h"
51 #include "AliAnalysisTask.h"
52 #include "AliInputEventHandler.h"
53 #include "AliLog.h"
54 #include "AliCentrality.h"
55 #include "AliV0vertexer.h"
56 #include "AliPIDResponse.h"
57 #include "AliMultiplicity.h"
58 #include "AliVertexerTracks.h"
59
60 #include "AliVEvent.h"
61 #include "AliVTrack.h"
62
63 #include "AliESDInputHandler.h"
64 #include "AliESDEvent.h"
65 #include "AliESDtrackCuts.h"
66 #include "AliESDpid.h"
67 #include "AliESDv0.h"
68
69 #include "AliAODInputHandler.h"
70 #include "AliAODEvent.h"
71 #include "AliAODTrack.h"
72 #include "AliAODv0.h"
73
74 #include "AliCFContainer.h"
75 #include "AliKFParticle.h"
76 #include "AliKFVertex.h"
77
78 #include "AliMCEventHandler.h"
79 #include "AliMCEvent.h"
80 #include "AliStack.h"
81
82 #include "AliAnalysisTaskLambdaNAOD.h"
83 using namespace std;
84
85
86 ClassImp(AliAnalysisTaskLambdaNAOD)
87
88
89 //________________________________________________________________________
90 AliAnalysisTaskLambdaNAOD::AliAnalysisTaskLambdaNAOD()
91 : AliAnalysisTaskSE(),
92 fAnalysisType("ESD"),
93 fEventHandler(0),
94 fESDtrackCutsV0(0),
95 fESDpid(0),
96 fPIDResponse(0),
97 fTreeV0(0),
98 fHistNumberOfEvents(0),
99 fHistTriggerStat(0),
100 fHistLambdaNeutronPtGen(0),
101 fHistAntiLambdaNeutronPtGen(0),
102 fHistLambdaNeutronInvaMassGen(0),
103 fHistAntiLambdaNeutronInvaMassGen(0),
104 fHistLambdaNeutronDecayLengthGen(0),
105 fHistAntiLambdaNeutronDecayLengthGen(0),
106 fHistLambdaNeutronPtAso(0),
107 fHistLambdaNeutronPtAsoCuts(0),
108 fHistAntiLambdaNeutronPtAso(0),
109 fHistAntiLambdaNeutronPtAsoCuts(0),
110 fHistLambdaNeutronInvaMassAso(0),
111 fHistAntiLambdaNeutronInvaMassAso(0),
112 fHistLambdaNeutronDecayLengthAso(0),
113 fHistAntiLambdaNeutronDecayLengthAso(0),
114 fof(0),
115 fHistArmenterosPodolanskiDeuteronPion(0),
116 fHistArmenterosPodolanskiAntiDeuteronPion(0),
117 fHistDeDxQA(0),
118 //fHistdEdxMC(0),
119 fNTriggers(5),
120 fMCtrue(0),
121 fOnlyQA(0),
122 fTriggerFired(),
1d281cbe 123 //fV0object(NULL),
0cea21e4 124 fItrk(0),
e2d2636c 125 fOutputContainer(NULL)
126 {
127 // default Constructor
128
129 // Define input and output slots here
130 }
131
132 //________________________________________________________________________
133 AliAnalysisTaskLambdaNAOD::AliAnalysisTaskLambdaNAOD(const char *name)
134 : AliAnalysisTaskSE(name),
135 fAnalysisType("ESD"),
136 fEventHandler(0),
137 fESDtrackCutsV0(0),
138 fESDpid(0),
139 fPIDResponse(0),
140 fTreeV0(0),
141 fHistNumberOfEvents(0),
142 fHistTriggerStat(0),
143 fHistLambdaNeutronPtGen(0),
144 fHistAntiLambdaNeutronPtGen(0),
145 fHistLambdaNeutronInvaMassGen(0),
146 fHistAntiLambdaNeutronInvaMassGen(0),
147 fHistLambdaNeutronDecayLengthGen(0),
148 fHistAntiLambdaNeutronDecayLengthGen(0),
149 fHistLambdaNeutronPtAso(0),
150 fHistLambdaNeutronPtAsoCuts(0),
151 fHistAntiLambdaNeutronPtAso(0),
152 fHistAntiLambdaNeutronPtAsoCuts(0),
153 fHistLambdaNeutronInvaMassAso(0),
154 fHistAntiLambdaNeutronInvaMassAso(0),
155 fHistLambdaNeutronDecayLengthAso(0),
156 fHistAntiLambdaNeutronDecayLengthAso(0),
157 fof(0),
158 fHistArmenterosPodolanskiDeuteronPion(0),
159 fHistArmenterosPodolanskiAntiDeuteronPion(0),
160 fHistDeDxQA(0),
161 //fHistdEdxMC(0),
162 fNTriggers(5),
163 fMCtrue(0),
164 fOnlyQA(0),
165 fTriggerFired(),
1d281cbe 166 //fV0object(NULL),
0cea21e4 167 fItrk(0),
e2d2636c 168 fOutputContainer(NULL)
169 {
170 // Constructor
171
172 // Define input and output slots here
173 // Input slot #0 works with a TChain
174 DefineInput(0, TChain::Class());
175 // Output slot #0 writes into a TH1 container
176 DefineOutput(1, TObjArray::Class());
177 DefineOutput(2, TTree::Class());
178
179 //ESD Track Cuts for v0's
180 fESDtrackCutsV0 = new AliESDtrackCuts("AliESDtrackCutsV0","AliESDtrackCutsV0");
181 fESDtrackCutsV0->SetAcceptKinkDaughters(kFALSE);
182 fESDtrackCutsV0->SetMinNClustersTPC(60);
183 fESDtrackCutsV0->SetMaxChi2PerClusterTPC(5);
184 fESDtrackCutsV0->SetRequireTPCRefit(kTRUE);
185 // fESDtrackCutsV0->SetMinNClustersITS(1); // to be tested if this cut is not too strong
186 fESDtrackCutsV0->SetEtaRange(-0.9,0.9);
187
188
189 fMCtrue = kTRUE;
190 fOnlyQA = kTRUE;
191
192 }
193
194 //____________________________________________________________
195 AliAnalysisTaskLambdaNAOD::~AliAnalysisTaskLambdaNAOD(){
196 //
197 // Destructor
198 //
199 //if (fESD) delete fESD;
200 if (fESDtrackCutsV0) delete fESDtrackCutsV0;
201
202 }
203
204 //____________________________________________________________
205 const Int_t AliAnalysisTaskLambdaNAOD::fgkPdgCode[] = {
206 211, //PionPlus
207 -211, //PionMinus
208 2212, //Proton
209 -2212, //Anti-Proton
210 1000010020, //Deuteron
211 -1000010020, //Anti-Deuteron
212 1000020030, //Helium3
213 -1000020030, //Anti-Helium3
214 3122, //Lambda
215 -3122, //Anti-Lambda
216 1010000020, //LambdaNeutron
217 -1010000020, //Anti-Lambda-Neutron
218 1010010030, //Hypertriton
219 -1010010030 //Anti-Hypertriton
220 };
221
222 //____________________________________________________________
223 const Double_t AliAnalysisTaskLambdaNAOD::fgkMass[] = {
224 0.13957, //Pion
225 0.93827, //Proton
226 1.875612, //Deuteron
227 2.808921, //Triton
228 2.80941 //Helium3 Quelle: Wolfram Alpha
229 };
230
231 //________________________________________________________________________
232 void AliAnalysisTaskLambdaNAOD::UserCreateOutputObjects(){
233
234 // New PID object
235 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
236 AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
237 fPIDResponse = inputHandler->GetPIDResponse();
238
239 // Create histograms for MC
240 //generated histogramms
241 fHistLambdaNeutronPtGen = new TH1F("fHistLambdaNeutronPtGen", "Generated #Lambdan", 201, 0., 10.1);
242 fHistLambdaNeutronPtGen->GetYaxis()->SetTitle("Counts");
243 fHistLambdaNeutronPtGen->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
244
245 fHistAntiLambdaNeutronPtGen = new TH1F("fHistAntiLambdaNeutronPtGen", "Generated #bar{#Lambdan}", 201, 0., 10.1);
246 fHistAntiLambdaNeutronPtGen->GetYaxis()->SetTitle("Counts");
247 fHistAntiLambdaNeutronPtGen->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
248
249 fHistLambdaNeutronInvaMassGen = new TH1F("fHistLambdaNeutronInvaMassGen", "Generated #Lambdan ", 100, 2.0, 2.1);
250 fHistLambdaNeutronInvaMassGen->GetYaxis()->SetTitle("Counts");
251 fHistLambdaNeutronInvaMassGen->GetXaxis()->SetTitle("Invariant mass (d #pi^{-}) (GeV/#it{c}^{2})");
252
253 fHistAntiLambdaNeutronInvaMassGen = new TH1F("fHistAntiLambdaNeutronInvaMassGen", "Generated #bar{#Lambdan}", 100, 2.0, 2.1);
254 fHistAntiLambdaNeutronInvaMassGen->GetYaxis()->SetTitle("Counts");
255 fHistAntiLambdaNeutronInvaMassGen->GetXaxis()->SetTitle("Invariant mass (#bar{d} #pi^{+}) (GeV/#it{c}^{2})");
256
257 fHistLambdaNeutronDecayLengthGen = new TH1F("fHistLambdaNeutronDecayLengthGen", "Generated #Lambdan", 401, 0., 400.1);
258 fHistLambdaNeutronDecayLengthGen->GetYaxis()->SetTitle("Counts");
259 fHistLambdaNeutronDecayLengthGen->GetXaxis()->SetTitle("#it{decay length} (cm)");
260
261 fHistAntiLambdaNeutronDecayLengthGen = new TH1F("fHistAntiLambdaNeutronDecayLengthGen", "Generated #bar{#Lambdan}", 401, 0., 400.1);
262 fHistAntiLambdaNeutronDecayLengthGen->GetYaxis()->SetTitle("Counts");
263 fHistAntiLambdaNeutronDecayLengthGen->GetXaxis()->SetTitle("#it{decay length} (cm)");
264
265 //assoziated (reconstracted) histogramms
266 fHistLambdaNeutronPtAso = new TH1F("fHistLambdaNeutronPtAso", "Associated #Lambdan", 201, 0., 10.1);
267 fHistLambdaNeutronPtAso->GetYaxis()->SetTitle("Counts");
268 fHistLambdaNeutronPtAso->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
269
270 fHistLambdaNeutronPtAsoCuts = new TH1F("fHistLambdaNeutronPtAsoCuts", "Associated #Lambdan Cuts", 201, 0., 10.1);
271 fHistLambdaNeutronPtAsoCuts->GetYaxis()->SetTitle("Counts");
272 fHistLambdaNeutronPtAsoCuts->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
273
274 fHistAntiLambdaNeutronPtAsoCuts = new TH1F("fHistAntiLambdaNeutronPtAsoCuts", " associated #bar{#Lambdan} Cuts", 201, 0., 10.1);
275 fHistAntiLambdaNeutronPtAsoCuts->GetYaxis()->SetTitle("Counts");
276 fHistAntiLambdaNeutronPtAsoCuts->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
277
278 fHistAntiLambdaNeutronPtAso = new TH1F("fHistAntiLambdaNeutronPtAso", " associated #bar{#Lambdan}", 201, 0., 10.1);
279 fHistAntiLambdaNeutronPtAso->GetYaxis()->SetTitle("Counts");
280 fHistAntiLambdaNeutronPtAso->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
281
282 fHistLambdaNeutronInvaMassAso = new TH1F("fHistLambdaNeutronInvaMassAso", "Associated #Lambdan", 100, 2.0, 2.1);
283 fHistLambdaNeutronInvaMassAso->GetYaxis()->SetTitle("Counts");
284 fHistLambdaNeutronInvaMassAso->GetXaxis()->SetTitle("Invariant mass (d #pi^{-}) (GeV/#it{c}^{2})");
285
286 fHistAntiLambdaNeutronInvaMassAso = new TH1F("fHistAntiLambdaNeutronInvaMassAso", " Associated #bar{#Lambdan}", 100, 2.0, 2.1);
287 fHistAntiLambdaNeutronInvaMassAso->GetYaxis()->SetTitle("Counts");
288 fHistAntiLambdaNeutronInvaMassAso->GetXaxis()->SetTitle("Invariant mass (#bar{d} #pi^{+}) (GeV/#it{c}^{2})");
289
290 fHistLambdaNeutronDecayLengthAso = new TH1F("fHistLambdaNeutronDecayLengthAso", "Associated #Lambdan", 401, 0., 400.1);
291 fHistLambdaNeutronDecayLengthAso->GetYaxis()->SetTitle("Counts");
292 fHistLambdaNeutronDecayLengthAso->GetXaxis()->SetTitle("#it{decay length} (cm)");
293
294 fHistAntiLambdaNeutronDecayLengthAso = new TH1F("fHistAntiLambdaNeutronDecayLengthAso", "Associated #bar{#Lambdan}", 401, 0., 400.1);
295 fHistAntiLambdaNeutronDecayLengthAso->GetYaxis()->SetTitle("Counts");
296 fHistAntiLambdaNeutronDecayLengthAso->GetXaxis()->SetTitle("#it{decay length} (cm)");
297
298 fHistLambdaNeutronPtAso = new TH1F("fHistLambdaNeutronPtAso", "Associated #Lambdan", 201, 0., 10.1);
299 fHistLambdaNeutronPtAso->GetYaxis()->SetTitle("Counts");
300 fHistLambdaNeutronPtAso->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
301
302 //Tree
303 //------------ Tree and branch definitions ----------------//
304 fTreeV0 = new TTree("tree", "fTreeV0");
305 fTreeV0->Branch("fItrk", &fItrk, "fItrk/I");
306
1d281cbe 307 //fTreeV0->Branch("fV0object",&fV0object,"fV0object[fItrk]");
e2d2636c 308 fTreeV0->Branch("fV0finder",fV0finder,"fV0finder[fItrk]/I");
309 fTreeV0->Branch("fkMB",fkMB,"fkMB[fItrk]/I");
310 fTreeV0->Branch("fkCentral",fkCentral,"fkCentral[fItrk]/I");
311 fTreeV0->Branch("fkSemiCentral",fkSemiCentral,"fkSemiCentral[fItrk]/I");
1d281cbe 312 //fTreeV0->Branch("fkEMCEJE",fkEMCEJE,"fkEMCEJE[fItrk]/I");
313 //fTreeV0->Branch("fkEMCEGA",fkEMCEGA,"fkEMCEGA[fItrk]/I");
e2d2636c 314
315 fTreeV0->Branch("fPtotN",fPtotN,"fPtotN[fItrk]/D");
316 fTreeV0->Branch("fPtotP",fPtotP,"fPtotP[fItrk]/D");
317 fTreeV0->Branch("fMotherPt",fMotherPt,"fMotherPt[fItrk]/D");
318 fTreeV0->Branch("fdEdxN",fdEdxN,"fdEdxN[fItrk]/D");
319 fTreeV0->Branch("fdEdxP",fdEdxP,"fdEdxP[fItrk]/D");
320 fTreeV0->Branch("fSignN",fSignN,"fSignN[fItrk]/D");
321 fTreeV0->Branch("fSignP",fSignP,"fSignP[fItrk]/D");
322
323 fTreeV0->Branch("fDCAv0",fDCAv0,"fDCAv0[fItrk]/F"); //Dca v0 Daughters
324 fTreeV0->Branch("fCosinePAv0",fCosinePAv0,"fCosinePAv0[fItrk]/F"); //Cosine of Pionting Angle
325 fTreeV0->Branch("fDecayRadiusTree",fDecayRadiusTree,"fDecayRadiusTree[fItrk]/F"); //decay radius
326
327 fTreeV0->Branch("fChargeComboDeuteronPionTree",fChargeComboDeuteronPionTree,"fChargeComboDeuteronPionTree[fItrk]/I");
328 fTreeV0->Branch("fInvaMassDeuteronPionTree",fInvaMassDeuteronPionTree,"fInvaMassDeuteronPionTree[fItrk]/D"); //invariant mass
329
330 fTreeV0->Branch("fIsCorrectlyAssociated",fIsCorrectlyAssociated,"fIsCorrectlyAssociated[fItrk]/O"); //associated hypertriton
331
1d281cbe 332 //fTreeV0->Branch("fAmenterosAlphaTree",fAmenterosAlphaTree,"fAmenterosAlphaTree[fItrk]/D");
333 //fTreeV0->Branch("fAmenterosQtTree",fAmenterosQtTree,"fAmenterosQtTree[fItrk]/D");
e2d2636c 334 fTreeV0->Branch("fRotationTree",fRotationTree,"fRotationTree[fItrk]/I");
335
336 //Armenteros-Podolanski
337 fHistArmenterosPodolanskiDeuteronPion= new TH2F("fHistArmenterosPodolanskiDeuteronPion", "Armenteros-Podolanski d #pi^{-}", 200,-1.0,1.0, 500,0,1);
338 fHistArmenterosPodolanskiDeuteronPion->GetXaxis()->SetTitle("#alpha");
339 fHistArmenterosPodolanskiDeuteronPion->GetYaxis()->SetTitle("q_{t}");
340 fHistArmenterosPodolanskiDeuteronPion->SetMarkerStyle(kFullCircle);
341
342 fHistArmenterosPodolanskiAntiDeuteronPion= new TH2F("fHistArmenterosPodolanskiAntiDeuteronPion", "Armenteros-Podolanski #bar{d} #pi^{+}", 200,-1.0,1.0, 500,0,1);
343 fHistArmenterosPodolanskiAntiDeuteronPion->GetXaxis()->SetTitle("#alpha");
344 fHistArmenterosPodolanskiAntiDeuteronPion->GetYaxis()->SetTitle("q_{t}");
345 fHistArmenterosPodolanskiAntiDeuteronPion->SetMarkerStyle(kFullCircle);
346
347 // control histograms
348 fof = new TH2F("fof", "OnTheFlyStatus ",5,0.5,5.5,2,-0.5,1.5);
349 fof->GetYaxis()->SetBinLabel(1,"offline");
350 fof->GetYaxis()->SetBinLabel(2,"onTheFly");
351 fof->GetXaxis()->SetBinLabel(1,"total");
352 fof->GetXaxis()->SetBinLabel(2,"dcaCut");
353 fof->GetXaxis()->SetBinLabel(3,"cosCut");
354 fof->GetXaxis()->SetBinLabel(4,"nucleonPID");
355 fof->GetXaxis()->SetBinLabel(5,"pionPID");
356 //fof->GetXaxis()->SetBinLabel(6,"decayRadiusCut");
357 //fof->SetMarkerStyle(kFullCircle);
358
359 //histogram to count number of events
360 fHistNumberOfEvents = new TH1F("fHistNumberOfEvents", "Number of events", 11, -0.5, 10.5);
361 fHistNumberOfEvents ->GetXaxis()->SetTitle("Centrality");
362 fHistNumberOfEvents ->GetYaxis()->SetTitle("Entries");
363
364 //trigger statitics histogram
365 fHistTriggerStat = new TH1F("fHistTriggerStat","Trigger statistics", fNTriggers,-0.5,fNTriggers-0.5);
366
367 const Char_t* aTriggerNames[] = { "kMB", "kCentral", "kSemiCentral", "kEMCEJE", "kEMCEGA" };
368 for ( Int_t ii=0; ii < fNTriggers; ii++ )
369 fHistTriggerStat->GetXaxis()->SetBinLabel(ii+1, aTriggerNames[ii]);
370
371 //QA dE/dx
372 fHistDeDxQA = new TH3F("fHistDeDxQA", "QA dE/dx", 400, -6.0, 6.0, 500, 0.0, 2000, 15, -0.5, 14.5);
373 fHistDeDxQA->GetYaxis()->SetTitle("TPC Signal");
374 fHistDeDxQA->GetXaxis()->SetTitle("p (GeV/c)");
375 fHistDeDxQA->GetZaxis()->SetBinLabel(0,"all pos v0 tracks");
376 fHistDeDxQA->GetZaxis()->SetBinLabel(1,"all neg v0 tracks");
377 fHistDeDxQA->GetZaxis()->SetBinLabel(2,"all neg deuteron");
378 fHistDeDxQA->GetZaxis()->SetBinLabel(3,"all pos deuteron");
379 fHistDeDxQA->GetZaxis()->SetBinLabel(6,"all selected tracks");
380 fHistDeDxQA->GetZaxis()->SetBinLabel(7,"neg deuteron for deuteron+pion");
381 fHistDeDxQA->GetZaxis()->SetBinLabel(8,"pos pion for deuteron+pion");
382 fHistDeDxQA->GetZaxis()->SetBinLabel(9,"pos deuteron for deuteron+pion");
383 fHistDeDxQA->GetZaxis()->SetBinLabel(10,"neg pion for deuteron+pion");
384
385 //Define and fill the OutputContainer
386 fOutputContainer = new TObjArray(1);
387 fOutputContainer->SetOwner(kTRUE);
388 fOutputContainer->SetName(GetName()) ;
389 fOutputContainer->Add(fHistArmenterosPodolanskiDeuteronPion);
390 fOutputContainer->Add(fHistArmenterosPodolanskiAntiDeuteronPion);
391 fOutputContainer->Add(fof);
392 fOutputContainer->Add(fHistDeDxQA);
393 fOutputContainer->Add(fHistNumberOfEvents);
394 fOutputContainer->Add(fHistTriggerStat);
395 fOutputContainer->Add(fHistLambdaNeutronPtGen);
396 fOutputContainer->Add(fHistAntiLambdaNeutronPtGen);
397 fOutputContainer->Add(fHistLambdaNeutronInvaMassGen);
398 fOutputContainer->Add(fHistAntiLambdaNeutronInvaMassGen);
399 fOutputContainer->Add(fHistLambdaNeutronDecayLengthGen);
400 fOutputContainer->Add(fHistAntiLambdaNeutronDecayLengthGen);
401 fOutputContainer->Add(fHistLambdaNeutronPtAso);
402 fOutputContainer->Add(fHistLambdaNeutronPtAsoCuts);
403 fOutputContainer->Add(fHistAntiLambdaNeutronPtAso);
404 fOutputContainer->Add(fHistAntiLambdaNeutronPtAsoCuts);
405 fOutputContainer->Add(fHistLambdaNeutronInvaMassAso);
406 fOutputContainer->Add(fHistAntiLambdaNeutronInvaMassAso);
407 fOutputContainer->Add(fHistLambdaNeutronDecayLengthAso);
408 fOutputContainer->Add(fHistAntiLambdaNeutronDecayLengthAso);
409 }
410
411 //________________________________________________________________________
412 void AliAnalysisTaskLambdaNAOD::UserExec(Option_t *){
413
414 // Main loop
415 // Called for each event
416
417 //Data
418 //get Event-Handler for the trigger information
419 fEventHandler= dynamic_cast<AliInputEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
420 if (!fEventHandler) {
421 AliError("Could not get InputHandler");
422 //return -1;
423 return;
424 }
425
426 // Monte Carlo
427 AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
428 if (!eventHandler){
429 //Printf("ERROR: Could not retrieve MC event handler");
430 fMCtrue = kFALSE;
431 }
432
433 AliMCEvent* mcEvent = 0x0;
434 AliStack* stack = 0x0;
435 if (eventHandler) mcEvent = eventHandler->MCEvent();
436 if (!mcEvent){
437 //Printf("ERROR: Could not retrieve MC event");
438 if (fMCtrue) return;
439 }
440
441 if (fMCtrue){
442 stack = mcEvent->Stack();
443 if (!stack) return;
444 }
445
446 //look for the generated particles
447 MCGenerated(stack);
448
1d281cbe 449
e2d2636c 450 if (SetupEvent() < 0) {
451 PostData(1, fOutputContainer);
452 return;
1d281cbe 453 }
e2d2636c 454
455 //DATA
456 AliESDEvent *fESDevent = 0x0;
457 AliAODEvent *fAODevent = 0x0;
458
459
460 if(!fMCtrue){
461 if(!fPIDResponse) {
462 AliError("Cannot get pid response");
463 return;
464 }
465 }
466
467
0cea21e4 468 Int_t centrality = -1;
e2d2636c 469 Double_t vertex[3] = {-100.0, -100.0, -100.0};
470
471 //Initialisation of the event and basic event cuts:
472 //1.) ESDs:
473 if (fAnalysisType == "ESD") {
474
475 fESDevent=dynamic_cast<AliESDEvent*>(InputEvent());
476 if (!fESDevent) {
477 AliWarning("ERROR: fESDevent not available \n");
478 return;
479 }
480
481 //Basic event cuts:
482 //1.) vertex existence
483 const AliESDVertex *vertexESD = fESDevent->GetPrimaryVertexTracks();
484 if (vertexESD->GetNContributors()<1)
485 {
486 // SPD vertex
487 vertexESD = fESDevent->GetPrimaryVertexSPD();
488 if(vertexESD->GetNContributors()<1) {
489 PostData(1, fOutputContainer);
490 return;
491 }
492 }
493
494 vertexESD->GetXYZ(vertex);
495
496 //2. vertex position within 10 cm
497 if (TMath::Abs(vertexESD->GetZv()) > 10) return;
498
499 //centrality selection in PbPb
500 if (fESDevent->GetEventSpecie() == 4) //species == 4 == PbPb
501 { // PbPb
502 AliCentrality *esdCentrality = fESDevent->GetCentrality();
503 centrality = esdCentrality->GetCentralityClass10("V0M"); // centrality percentile determined with V0
504 if(!fMCtrue){
0cea21e4 505 if (centrality < 0. || centrality > 8. ) return; //select only events with centralities between 0 and 80 %
e2d2636c 506 }
507 }
508
0cea21e4 509 //cout << "centrality "<< centrality << endl;
510 // count number of events
511 fHistNumberOfEvents->Fill(centrality);
e2d2636c 512
513 }
514
515 //2.) AODs:
516 else if (fAnalysisType == "AOD") {
517
518 fAODevent = dynamic_cast<AliAODEvent*>( InputEvent() );
519 if (!fAODevent) {
520 AliWarning("ERROR: lAODevent not available \n");
521 return;
522 }
523
524 const AliAODVertex *vertexAOD = fAODevent->GetPrimaryVertex();
525 if (vertexAOD->GetNContributors()<1)
526 {
527 // SPD vertex
528 vertexAOD = fAODevent->GetPrimaryVertexSPD();
529 if(vertexAOD->GetNContributors()<1) {
530 PostData(1, fOutputContainer);
531 return;
532 }
533 }
534 vertexAOD->GetXYZ(vertex);
535
536 //2. vertex position within 10 cm
537 if (TMath::Abs(vertex[2]) > 10) return;
538
539 //centrality selection
540 AliCentrality *aodCentrality = fAODevent->GetCentrality();
541 centrality = aodCentrality->GetCentralityClass10("V0M"); // centrality percentile determined with V0
542 if (centrality < 0 || centrality > 8) return; //select only events with centralities between 0 and 80 %
543
544 // count number of events
545 fHistNumberOfEvents->Fill(centrality);
546
547 } else {
548
549 Printf("Analysis type (ESD or AOD) not specified \n");
550 return;
551
552 }
553
554
555 //v0 loop
556
557
558 fItrk = 0;
559
560 Int_t runNumber = 0;
561 runNumber = (InputEvent())->GetRunNumber();
562
563 Int_t nTrackMultiplicity = -1;
564 nTrackMultiplicity = (InputEvent())->GetNumberOfTracks();
565
566
567 for (Int_t ivertex=0; ivertex<(InputEvent()->GetNumberOfV0s()); ivertex++) { //beginn v0 loop
568
569 AliESDv0 * v0ESD = 0x0;
570 AliAODv0 * v0AOD = 0x0;
571
572 Bool_t v0ChargesCorrect = kTRUE;
573 Bool_t testTrackCuts = kFALSE;
0cea21e4 574 //Bool_t testFilterBit = kFALSE;
e2d2636c 575
576 Int_t of = 7;
577 Int_t onFlyStatus = 5;
578
579 Float_t dcaV0 = -1;
580 Float_t cosPointing= 2;
581 Float_t decayRadius= -1;
582
583 AliVTrack * trackN = 0x0;
584 AliVTrack * trackP = 0x0;
585
586 Double_t ptotN = -1000;
587 Double_t ptotP = -1000;
588
589 Int_t chargeComboDeuteronPion = -1;
590
591 Double_t momentumPion[3]={0,0,0};
592 Double_t momentumPionRot[3]={0,0,0};
593 Double_t momentumDeuteron[3]={0,0,0};
594 Double_t momentumDeuteronRot[3]={0,0,0};
595 Double_t momentumMother[3] = {0,0,0};
596
597 Double_t transversMomentumPion = 0;
598 Double_t transversMomentumDeuteron = 0;
599
600 Double_t transversMomentumMother = 0;
601 Double_t longitudinalMomentumMother = 0;
602
603 Double_t totalEnergyMother = 0;
604 Double_t energyPion = 0;
605 Double_t energyDeuteron = 0;
606
607 Double_t rapidity = 2;
608
609 TVector3 vecPion(0,0,0);
610 TVector3 vecDeuteron(0,0,0);
611 TVector3 vecMother(0,0,0);
612
613 Double_t alpha = 2;
614 Double_t qt = -1;
615
616 Double_t thetaPion = 0;
617 Double_t thetaDeuteron = 0;
618
619 Double_t phi=0;
620 Double_t invaMassDeuteronPion = 0;
621
622 Double_t e12 = 0;
623 Double_t r12 = 0;
624 Double_t d22 = 0;
625 Double_t dr22 = 0;
626
627 //Tree variables
628 fV0finder[fItrk] = -1;
629 fkMB[fItrk] = -1;
630 fkCentral[fItrk] = -1;
631 fkSemiCentral[fItrk] = -1;
1d281cbe 632 //fkEMCEJE[fItrk] = -1;
633 //fkEMCEGA[fItrk] = -1;
e2d2636c 634
635 fPtotN[fItrk] = -1000;
636 fPtotP[fItrk] = -1000;
637 fMotherPt[fItrk] = -1000;
638
639 fdEdxN[fItrk] = -1;
640 fdEdxP[fItrk] = -1;
641 fSignN[fItrk] = 0;
642 fSignP[fItrk] = 0;
643
644 fDCAv0[fItrk] = -1;
645 fCosinePAv0[fItrk] = -2;
646 fDecayRadiusTree[fItrk] = -1;
647
648 fInvaMassDeuteronPionTree[fItrk] = 0;
649 fChargeComboDeuteronPionTree[fItrk] = -1;
650
1d281cbe 651 //fAmenterosAlphaTree[fItrk] = 2;
652 //fAmenterosQtTree[fItrk] = -1;
e2d2636c 653
654 //Get v0 object
655 if(fAnalysisType == "ESD")v0ESD = fESDevent->GetV0(ivertex);
656 if(fAnalysisType == "AOD")v0AOD = fAODevent->GetV0(ivertex);
657
658
659 //distinguish between the two V0 finders: 0 offline V0 finder / 1 online V0 finder
660
661 if(fAnalysisType == "ESD") of = v0ESD->GetOnFlyStatus();
662 if(fAnalysisType == "AOD") of = v0AOD->GetOnFlyStatus();
663
664 if(of) onFlyStatus= 1;
665 if(!of) onFlyStatus= 0;
666
667 if(onFlyStatus==0)fof->Fill(1,0);
668 if(onFlyStatus==1)fof->Fill(1,1);
669
670 //Get dca, cos of pointing angle and decay radius
671
672
673 if(fAnalysisType == "ESD")
674 {
675 dcaV0 = v0ESD->GetDcaV0Daughters();
676 cosPointing = v0ESD->GetV0CosineOfPointingAngle();
677 decayRadius = v0ESD->GetRr();
678 }
679
680 if(fAnalysisType == "AOD")
681 {
682 dcaV0 = v0AOD->DcaV0Daughters();
683 cosPointing = v0AOD->CosPointingAngle(vertex);
684 decayRadius = v0AOD->DecayLengthV0(vertex);
685 }
686
687
688 // select coresponding tracks
689 if(fAnalysisType == "ESD")
690 {
691 trackN = fESDevent->GetTrack(v0ESD->GetIndex(0));
692 trackP = fESDevent->GetTrack(v0ESD->GetIndex(1));
693
694 if (trackN->Charge() > 0) // switch because of bug in V0 interface
695 {
696 trackN = fESDevent->GetTrack(v0ESD->GetIndex(1));
697 trackP = fESDevent->GetTrack(v0ESD->GetIndex(0));
698 v0ChargesCorrect = kFALSE;
699 }
700 //Track-Cuts
701 testTrackCuts = TrackCuts(trackN,testTrackCuts);
702 if(testTrackCuts == kFALSE) continue;
703 testTrackCuts = TrackCuts(trackP,testTrackCuts);
704 if(testTrackCuts == kFALSE) continue;
705 }
706
707 if(fAnalysisType == "AOD")
708 {
709 trackN = dynamic_cast<AliVTrack*>(v0AOD->GetDaughter(0));
710 trackP = dynamic_cast<AliVTrack*>(v0AOD->GetDaughter(1));
711
712 if (trackN->Charge() > 0) // switch because of bug in V0 interface
713 {
714 trackN = dynamic_cast<AliVTrack*>(v0AOD->GetDaughter(1));
715 trackP = dynamic_cast<AliVTrack*>(v0AOD->GetDaughter(0));
716 v0ChargesCorrect = kFALSE;
717 }
718 //Test-Filterbit - only use for track base analysis -- not for V0 candidates -- the AOD V0 candidates sould be a copy of the ESD candidates, BUT NOT for AOD0086 -> here there wqas a wider cos-cutto have more candidates in the AODs
719 //testFilterBit = FilterBit(trackN,testFilterBit);
720 //if(testFilterBit == kFALSE) continue;
721 //testFilterBit = FilterBit(trackP,testFilterBit);
722 //if(testFilterBit == kFALSE) continue;
723 //Track-Cuts
1d281cbe 724 /*testTrackCuts = TrackCuts(trackN,testTrackCuts);
e2d2636c 725 if(testTrackCuts == kFALSE) continue;
726 testTrackCuts = TrackCuts(trackP,testTrackCuts);
1d281cbe 727 if(testTrackCuts == kFALSE) continue;*/
e2d2636c 728 }
729
730 //Get the total momentum for each track ---> at the inner readout of the TPC???? // momenta a always positive
731 if(fAnalysisType == "AOD") {
732 ptotN = trackN->P(); //GetInnerParam()->GetP(); ******************FIXME: InnerParam
733 ptotP = trackP->P(); //GetInnerParam()->GetP(); ******************FIXME: InnerParam
734 }
735
736 if(fAnalysisType == "ESD") {
737 ptotN = MomentumInnerParam(trackN,ptotN);
738 ptotP = MomentumInnerParam(trackP,ptotP);
739 }
740
741 //fill QA dEdx with all V0 candidates
742 if(trackP) fHistDeDxQA->Fill(ptotP*trackP->Charge(), trackP->GetTPCsignal(),0);
743 if(trackN) fHistDeDxQA->Fill(ptotN*trackN->Charge(), trackN->GetTPCsignal(),1);
744
745 if (dcaV0 > 3 || dcaV0 < 1e-20) continue;
746
747 //Check how much the dca cut reduces the statistic (background) for the different VO-finder
748 if(onFlyStatus==0)fof->Fill(2,0);
749 if(onFlyStatus==1)fof->Fill(2,1);
750
751 if (cosPointing < 0.9) continue;
752
753 //Check how much the cos-of-the-pointing-angle-cut reduces the statistic (background) for the different VO-finder
754 if(onFlyStatus==0)fof->Fill(3,0);
755 if(onFlyStatus==1)fof->Fill(3,1);
756
757 //deuteron PID via specific energy loss in the TPC
758 Bool_t isDeuteron[3] = {kFALSE,kFALSE,kFALSE}; //0 = posDeuteron, 1 = negDeuteron, 2 = trackN is deuteron
759
760 DeuteronPID(trackP,trackN,ptotP,ptotN,runNumber,isDeuteron);
761
762 //if(!isDeuteron[0] && !isDeuteron[1]) continue;
763 if((isDeuteron[0]==kFALSE) && (isDeuteron[1]==kFALSE)) continue;
764
765
766 //Check how much the nuclei PID cuts reduce the statistics (background) for the two V0-finders
767 if(onFlyStatus==0)fof->Fill(4,0);
768 if(onFlyStatus==1)fof->Fill(4,1);
769
770 //Fill the QA dEdx with deuterons and helium3 after the nuclei PID cut
771 if(isDeuteron[1]) fHistDeDxQA->Fill(ptotN*trackN->Charge(), trackN->GetTPCsignal(),2);
772 if(isDeuteron[0]) fHistDeDxQA->Fill(ptotP*trackP->Charge(), trackP->GetTPCsignal(),3);
773
774 //deuteron PID via specific energy loss in the TPC
775 Bool_t isPion[2] = {kFALSE,kFALSE}; //0 = posPion, 1 = negPion
776
777 PionPID(trackP,trackN,ptotP,ptotN,runNumber,isPion);
778
779 //if(isDeuteron[0] && !isPion[1]) continue; //pos deuteron and neg Pion
780 //if(isDeuteron[1] && !isPion[0]) continue; //neg deuteron and pos Pion
781
782 //Check how much the pion PID cuts reduce the statistics (background) for the two V0-finders
783 if(onFlyStatus==0)fof->Fill(5,0);
784 if(onFlyStatus==1)fof->Fill(5,1);
785
786
787 //Save the different charge combinations to differentiat between particles and anit-particles:
788 // -/+ = Anti-deuteron + positive Pion
789 // -/- = Anti-deuteron + negative Pion
790 // +/- = deuteron + negative Pion
791 // +/+ = deuteron + positive Pion
792
793 //Like-sign
794 if (trackN->Charge()<0 && trackP->Charge()<0 && isDeuteron[1]==kTRUE && isPion[1]==kTRUE) chargeComboDeuteronPion = 1; // -/-
795 if (trackN->Charge()>0 && trackP->Charge()>0 && isDeuteron[0]==kTRUE && isPion[0]==kTRUE ) chargeComboDeuteronPion = 3; // +/+
796
797 //unlinke-sign
798 if (trackN->Charge()<0 && trackP->Charge()>0 && isDeuteron[1]==kTRUE && isPion[0]==kTRUE ) chargeComboDeuteronPion = 0; // -/+
799 git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGLF/STRANGENESS/Hypernuclei/AliAnalysisTaskLambdaNAOD.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame