]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG2/SPECTRA/LambdaK0PbPb/AliAnalysisTaskLukeV0.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Fix Coverity
[u/mrichter/AliRoot.git] / PWG2 / SPECTRA / LambdaK0PbPb / AliAnalysisTaskLukeV0.cxx
CommitLineData
45b90328 1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/* $Id: AliAnalysisTaskLukeV0.cxx 46301 2011-01-06 14:25:27Z agheata $ */
17
18/* AliAnalysisTaskLukeV0.cxx
19 *
20 * Task analysing lambda, antilambda & K0 spectra
21 *
22 * Based on tutorial example from offline pages
23 * Edited by Arvinder Palaha
24 * Adapted by Luke Hanratty
137b565d 25 * ------------------------*
45b90328 26 *
27 */
28
29#include "AliAnalysisTaskLukeV0.h"
30
31#include "Riostream.h"
32#include "TChain.h"
33#include "TTree.h"
34#include "TH1F.h"
35#include "TH2F.h"
36#include "TCanvas.h"
37#include "TList.h"
38#include "TPDGCode.h"
39
40#include "AliAnalysisTaskSE.h"
41#include "AliAnalysisManager.h"
42#include "AliStack.h"
43#include "AliESDtrackCuts.h"
44#include "AliESDEvent.h"
45#include "AliESDv0.h"
46#include "AliESDInputHandler.h"
47#include "AliAODEvent.h"
48#include "AliMCEvent.h"
49#include "AliMCVertex.h"
50#include "AliPID.h"
51#include "AliPIDResponse.h"
52
53ClassImp(AliAnalysisTaskLukeV0)
54
55//________________________________________________________________________
56AliAnalysisTaskLukeV0::AliAnalysisTaskLukeV0() // All data members should be initialised here
57 :AliAnalysisTaskSE(),
58 fOutputList(0),
59 fTrackCuts(0),
60 fPIDResponse(0),
45b90328 61 fHistCosPA(0),
45b90328 62 fHistDCAV0Daughters(0),
63 fHistDecayL(0),
45b90328 64 fHistImpactxyN(0),
65 fHistImpactzN(0),
66 fHistImpactxyP(0),
67 fHistImpactzP(0),
889bc324 68 fHistMCLambdacTau(0),
69 fHistMCLambdaNotcTau(0),
70 fHistMCLambdaDecayL(0),
71 fHistMCLambdaNotDecayL(0),
45b90328 72 fHistMcNLambdaPrimary(0),
73 fHistMcNLambda(0),
74 fHistMcNAntilambda(0),
75 fHistMcNKshort(0),
45b90328 76 fHistMK0(0),
45b90328 77 fHistMLa(0),
45b90328 78 fHistMLb(0),
45b90328 79 fHistNLambda(0),
80 fHistNV0(0),
45b90328 81 fHistPtV0(0),
82 fHistPVZ(0),
83 fHistTauLa(0),
45b90328 84 fHistV0Z(0),
85 fHistZ(0),
45b90328 86 fHistBetheBlochTPCNeg(0),
87 fHistBetheBlochTPCPos(0),
45b90328 88 fHistImpactxyImpactz(0),
45b90328 89 fHistMcPMK0Pt(0),
90 fHistMcPMLaPt(0),
91 fHistMcPMLbPt(0),
45b90328 92 fHistMcV0MK0Pt(0),
93 fHistMcV0MLaPt(0),
94 fHistMcV0MLbPt(0),
45b90328 95 fHistMK0Pt(0),
45b90328 96 fHistMLaPt(0),
45b90328 97 fHistMLbPt(0),
45b90328 98 fHistPtArm(0),
45b90328 99 fHistRZ(0),
45b90328 100 fHistXZ(0),
101 fHistYZ(0) // The last in the above list should not have a comma after it
102{
103 // Dummy constructor ALWAYS needed for I/O.
104}
105
106//________________________________________________________________________
107AliAnalysisTaskLukeV0::AliAnalysisTaskLukeV0(const char *name) // All data members should be initialised here
108 :AliAnalysisTaskSE(name),
889bc324 109 fOutputList(0),
45b90328 110 fTrackCuts(0),
111 fPIDResponse(0),
45b90328 112 fHistCosPA(0),
45b90328 113 fHistDCAV0Daughters(0),
114 fHistDecayL(0),
45b90328 115 fHistImpactxyN(0),
116 fHistImpactzN(0),
117 fHistImpactxyP(0),
118 fHistImpactzP(0),
889bc324 119 fHistMCLambdacTau(0),
120 fHistMCLambdaNotcTau(0),
121 fHistMCLambdaDecayL(0),
122 fHistMCLambdaNotDecayL(0),
45b90328 123 fHistMcNLambdaPrimary(0),
124 fHistMcNLambda(0),
125 fHistMcNAntilambda(0),
126 fHistMcNKshort(0),
45b90328 127 fHistMK0(0),
45b90328 128 fHistMLa(0),
45b90328 129 fHistMLb(0),
45b90328 130 fHistNLambda(0),
131 fHistNV0(0),
45b90328 132 fHistPtV0(0),
133 fHistPVZ(0),
134 fHistTauLa(0),
45b90328 135 fHistV0Z(0),
136 fHistZ(0),
45b90328 137 fHistBetheBlochTPCNeg(0),
138 fHistBetheBlochTPCPos(0),
45b90328 139 fHistImpactxyImpactz(0),
45b90328 140 fHistMcPMK0Pt(0),
141 fHistMcPMLaPt(0),
142 fHistMcPMLbPt(0),
45b90328 143 fHistMcV0MK0Pt(0),
144 fHistMcV0MLaPt(0),
145 fHistMcV0MLbPt(0),
45b90328 146 fHistMK0Pt(0),
45b90328 147 fHistMLaPt(0),
45b90328 148 fHistMLbPt(0),
45b90328 149 fHistPtArm(0),
45b90328 150 fHistRZ(0),
45b90328 151 fHistXZ(0),
889bc324 152 fHistYZ(0) // The last in the above list should not have a comma after it
45b90328 153{
154 // Constructor
155 // Define input and output slots here (never in the dummy constructor)
156 // Input slot #0 works with a TChain - it is connected to the default input container
157 // Output slot #1 writes into a TH1 container
158 DefineOutput(1, TList::Class()); // for output list
159}
160
161//________________________________________________________________________
162AliAnalysisTaskLukeV0::~AliAnalysisTaskLukeV0()
163{
164 // Destructor. Clean-up the output list, but not the histograms that are put inside
165 // (the list is owner and will clean-up these histograms). Protect in PROOF case.
166 if (fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
167 delete fOutputList;
168 }
169 if (fTrackCuts) delete fTrackCuts;
170}
171
172//________________________________________________________________________
173void AliAnalysisTaskLukeV0::UserCreateOutputObjects()
174{
175 // Create histograms
176 // Called once (on the worker node)
177
178 fOutputList = new TList();
179 fOutputList->SetOwner(); // IMPORTANT!
180
181 fTrackCuts = new AliESDtrackCuts();
182
183 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
184 AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
185 fPIDResponse = inputHandler->GetPIDResponse();
186
45b90328 187 // lambda plot parameters
188 int div = 96;
189 float max = 1.2;
190 float min = 1.08;
191
192 // Create remaining histograms
193 // TH1F first
45b90328 194 fHistCosPA = new TH1F("fHistCosPA", "Cosine of Pointing Angle of V0s; Cos PA; N(v0s)",202,0.8,1.01);
45b90328 195 fHistDCAV0Daughters = new TH1F("fHistDCAV0Daughters", "DCA between V0 daughters; DCA (cm); N V0s", 100, 0, 2);
196 fHistDecayL = new TH1F("fHistDecayL", "Distance between V0 and PV; Distance(cm); N(v0s)",200,-0.1,30);
45b90328 197 fHistImpactxyN = new TH1F("fHistImpactxyN", "RSM DCA between negative particle and primary vertex in xy plane; RSM DCA (cm); N(v0s)",100,0,1);
198 fHistImpactzN = new TH1F("fHistImpactzN", "RSM DCA between negative particle and primary vertex in z direction; RSM DCA (cm); N(v0s)",100,0,1);
199 fHistImpactxyP = new TH1F("fHistImpactxyP", "RSM DCA between positive particle and primary vertex in xy plane; RSM DCA (cm); N(v0s)",100,0,1);
200 fHistImpactzP = new TH1F("fHistImpactzP", "RSM DCA between positive particle and primary vertex in z direction; RSM DCA (cm); N(v0s)",100,0,1);
889bc324 201 fHistMCLambdacTau = new TH1F("fHistMCLambdacTau", "Lifetime under Lambda mass hypothesis of MC lambda; Lifetime(s); N(v0s)",200,0,100);
202 fHistMCLambdaNotcTau = new TH1F("fHistMCLambdaNotcTau", "Lifetime under Lambda mass hypothesis of MC lambda background; Lifetime(s); N(v0s)",200,0,100);
203 fHistMCLambdaDecayL = new TH1F("fHistMCLambdaDecayL", "Distance between V0 and PV of MC Lambda; Distance(cm); N(v0s)",200,-0.1,30);
204 fHistMCLambdaNotDecayL = new TH1F("fHistMCLambdaNotDecayL", "Distance between V0 and PV of MC Lambda background; Distance(cm); N(v0s)",200,-0.1,30);
45b90328 205 fHistMcNLambdaPrimary = new TH1F("fHistMcNLambdaPrimary","Number of primary lambdas in MC; NLambdas; i",6,-0.25,2.25);
206 fHistMcNLambda = new TH1F("fHistMcNLambda","Number of lambdas in MC; NLambdas; i",31,-0.5,30);
207 fHistMcNAntilambda = new TH1F("fHistMcNAntilambda","Number of antilambdas in MC; NAntiLambdas; i",31,-0.5,30);
208 fHistMcNKshort = new TH1F("fHistMcNKshort","Number of K0s in MC; NKshort; i",31,-0.5,30);
45b90328 209 fHistMK0 = new TH1F("fHistMK0","K0Short Mass; M(#pi^{+}#pi^{-}) (GeV/c^{2}); dN/dM (0.12 GeV/c^{2})^{-1}",140,0.414,0.582);
45b90328 210 fHistMLa = new TH1F("fHistMLa","Lambda Mass; M(p#pi^{-}) (GeV/c^{2}); dN/dM (0.125 GeV/c^{2})^{-1}",div,min,max);
45b90328 211 fHistMLb = new TH1F("fHistMLb","AntiLambda Mass; M(#bar{p}#pi^{+}) (GeV/c^{2}); dN/dM (0.125 GeV/c^{2})^{-1}",div,min,max);
45b90328 212 fHistNLambda = new TH1F("fHistNLambda", "Number of lambda per event; N(lambda); N(events)",50,-0.5,49.5);
213 fHistNV0 = new TH1F("fHistNV0","V0 frequency distribution; Number of V0 Candidates",1000,0,100000);
45b90328 214 fHistPtV0 = new TH1F("fHistPtV0","V0 P_{T}; P_{T} (GeV/c);dN/dP_{T} (GeV/c)^{-1}",40,0.,4.);
215 fHistPVZ = new TH1F("fHistPVZ","Z primary; Z (cm); Counts",100,-10,10);
216 fHistTauLa = new TH1F("fHistTauLa", "Lifetime under Lambda mass hypothesis; Lifetime(s); N(v0s)",200,0,100);
45b90328 217 fHistV0Z = new TH1F("fHistV0Z","Z decay; Z (cm); Counts",100,-10,10);
218 fHistZ = new TH1F("fHistZ","Z decay - Z primary; Z (cm); Counts",100,-10,10);
219
220 //TH2F follow
45b90328 221 fHistBetheBlochTPCNeg = new TH2F("fHistBetheBlochTPCNeg","-dE/dX against Momentum for negative daughter from TPC; Log10 P (GeV); -dE/dx (keV/cm ?)",1000,-1,1,1000,0,200);
222 fHistBetheBlochTPCPos = new TH2F("fHistBetheBlochTPCPos","-dE/dX against Momentum for positive daughter from TPC; Log10 P (GeV); -dE/dx (keV/cm ?)",1000,-1,1,1000,0,200);
45b90328 223 fHistImpactxyImpactz = new TH2F("fHistImpactxyImpactz", "RSM DCA between negative particle and primary vertex in xy plane; RSM DCA xy (cm); RSM DCA z (cm)",100,0,1,100,0,10);
45b90328 224 fHistMcPMK0Pt = new TH2F("fHistMcPMK0Pt","Monte Carlo primary K0 Mass versus Pt; P_{perp} (GeV/c); K0 Mass (GeV/c^2)",200,0,10,140,0.414,0.582);
225 fHistMcPMLaPt = new TH2F("fHistMcPMLaPt","Monte Carlo primary (& sigma0) Lambda Mass versus Pt; P_{perp} (GeV/c); M(p#pi^{-}) (GeV/c^2)",200,0,10,96,1.08,1.2);
226 fHistMcPMLbPt = new TH2F("fHistMcPMLbPt","Monte Carlo primary (& sigma0) AntiLambda Mass versus Pt; P_{perp} (GeV/c); M(#bar{p}#pi^{+}) (GeV/c^2)",200,0,10,96,1.08,1.2);
45b90328 227 fHistMcV0MK0Pt = new TH2F("fHistMcV0MK0Pt","Monte Carlo V0s passing cuts. K0 Mass versus Pt; P_{perp} (GeV/c); K0 Mass (GeV/c^2)",200,0,10,140,0.414,0.582);
228 fHistMcV0MLaPt = new TH2F("fHistMcV0MLaPt","Monte Carlo V0s passing cuts. Lambda Mass versus Pt; P_{perp} (GeV/c); Lambda Mass (GeV/c^2)",200,0,10,96,1.08,1.2);
229 fHistMcV0MLbPt = new TH2F("fHistMcV0MLbPt","Monte Carlo V0s passing cuts. Antilambda Mass versus Pt; P_{perp} (GeV/c); Antilambda Mass (GeV/c^2)",200,0,10,96,1.08,1.2);
45b90328 230 fHistMK0Pt = new TH2F("fHistMK0Pt","K0 Mass versus Pt; P_{perp} (GeV/c); K0 Mass (GeV/c^2)",200,0,10,140,0.414,0.582);
45b90328 231 fHistMLaPt = new TH2F("fHistMLaPt","Lambda Mass versus Pt; P_{perp} (GeV/c); M(p#pi^{-}) (GeV/c^2)",200,0,10,96,1.08,1.2);
45b90328 232 fHistMLbPt = new TH2F("fHistMLbPt","AntiLambda Mass versus Pt; P_{perp} (GeV/c); M(#bar{p}#pi^{+}) (GeV/c^2)",200,0,10,96,1.08,1.2);
45b90328 233 fHistPtArm = new TH2F("fHistPtArm","Podolanski-Armenteros Plot; #alpha; P_{perp} (GeV/c)",40,-1,1,80,0,0.5);
45b90328 234 fHistPtV0Z = new TH2F("fHistPtV0Z","Pt of V0 vs Z position; Pt (GeV/c); Z (cm)",200,-0.1,1.9,200,-50,50);
235 fHistRZ = new TH2F("fHistRZ","R decay versus Z decay; Z (cm); R (cm)",100,-50,50,120,0,220);
45b90328 236 fHistXZ = new TH2F("fHistXZ","X decay versus Z decay; Z (cm); X (cm)",100,-50,50,200,-200,200);
237 fHistYZ = new TH2F("fHistYZ","Y decay versus Z decay; Z (cm); Y (cm)",100,-50,50,200,-200,200);
238
239
240
241 // All histograms must be added to output list
242
45b90328 243 fOutputList->Add(fHistCosPA);
45b90328 244 fOutputList->Add(fHistDCAV0Daughters);
245 fOutputList->Add(fHistDecayL);
45b90328 246 fOutputList->Add(fHistImpactxyN);
247 fOutputList->Add(fHistImpactzN);
248 fOutputList->Add(fHistImpactxyP);
249 fOutputList->Add(fHistImpactzP);
889bc324 250 fOutputList->Add(fHistMCLambdacTau);
251 fOutputList->Add(fHistMCLambdaNotcTau);
252 fOutputList->Add(fHistMCLambdaDecayL);
253 fOutputList->Add(fHistMCLambdaNotDecayL);
45b90328 254 fOutputList->Add(fHistMcNLambdaPrimary);
255 fOutputList->Add(fHistMcNLambda);
256 fOutputList->Add(fHistMcNAntilambda);
257 fOutputList->Add(fHistMcNKshort);
45b90328 258 fOutputList->Add(fHistMK0);
45b90328 259 fOutputList->Add(fHistMLa);
45b90328 260 fOutputList->Add(fHistMLb);
45b90328 261 fOutputList->Add(fHistNLambda);
262 fOutputList->Add(fHistNV0);
45b90328 263 fOutputList->Add(fHistPtV0);
264 fOutputList->Add(fHistPVZ);
265 fOutputList->Add(fHistTauLa);
45b90328 266 fOutputList->Add(fHistV0Z);
267 fOutputList->Add(fHistZ);
45b90328 268 fOutputList->Add(fHistBetheBlochTPCNeg);
269 fOutputList->Add(fHistBetheBlochTPCPos);
45b90328 270 fOutputList->Add(fHistImpactxyImpactz);
45b90328 271 fOutputList->Add(fHistMcPMK0Pt);
272 fOutputList->Add(fHistMcPMLaPt);
273 fOutputList->Add(fHistMcPMLbPt);
45b90328 274 fOutputList->Add(fHistMcV0MK0Pt);
275 fOutputList->Add(fHistMcV0MLaPt);
276 fOutputList->Add(fHistMcV0MLbPt);
45b90328 277 fOutputList->Add(fHistMK0Pt);
45b90328 278 fOutputList->Add(fHistMLaPt);
45b90328 279 fOutputList->Add(fHistMLbPt);
45b90328 280 fOutputList->Add(fHistPtArm);
45b90328 281 fOutputList->Add(fHistRZ);
45b90328 282 fOutputList->Add(fHistXZ);
283 fOutputList->Add(fHistYZ);
284
285
286 PostData(1, fOutputList); // Post data for ALL output slots >0 here, to get at least an empty histogram
287}
288
289//________________________________________________________________________
290void AliAnalysisTaskLukeV0::UserExec(Option_t *)
291{
292 // Main loop
293 // Called for each event
294
295 // paramaters used for most cuts, to minimise editing
296 double cutCosPa(0.998), cutcTau(2);
297 double cutNImpact(-999), cutDCA(0.4);
298 double cutBetheBloch(3);
299 double cutMinNClustersTPC(70), cutMaxChi2PerClusterTPC(-999);
889bc324 300 double isMonteCarlo(true);
45b90328 301 double cutEta(0.8);
302
303 //Track Cuts set here
304 if(cutMinNClustersTPC != -999)
305 {(fTrackCuts->SetMinNClustersTPC(int(cutMinNClustersTPC)));}
306 if(cutMaxChi2PerClusterTPC != -999)
307 {fTrackCuts->SetMaxChi2PerClusterTPC(cutMaxChi2PerClusterTPC);}
308 fTrackCuts->SetAcceptKinkDaughters(kFALSE);
309 fTrackCuts->SetRequireTPCRefit(kTRUE);
310
311
312 // Create pointer to reconstructed event
313
314 AliVEvent *event = InputEvent();
315 if (!event) { Printf("ERROR: Could not retrieve event"); return; }
316
317 // create pointer to event
318 AliESDEvent* fESD = dynamic_cast<AliESDEvent*>(event);
319 if (!fESD) {
320 AliError("Cannot get the ESD event");
321 return;
322 }
323
324 /*********************************************************************/
325 // MONTE CARLO SECTION
326 // This section loops over all MC tracks
327
328 int nLambdaMC = 0;
329 int nAntilambdaMC = 0;
330 int nKshortMC = 0;
331
332 if(isMonteCarlo)
333 {
334
335 // If the task accesses MC info, this can be done as in the commented block below:
336
337 // Create pointer to reconstructed event
338 AliMCEvent *mcEvent = MCEvent();
339 if (!mcEvent)
340 {
341 Printf("ERROR: Could not retrieve MC event");
342 //return;
343 }
344 else
345 {
346 Printf("MC particles: %d", mcEvent->GetNumberOfTracks());
347 }
348
349 // set up a stack for use in check for primary/stable particles
350 AliStack* mcStack = mcEvent->Stack();
351 if( !mcStack ) { Printf( "Stack not available"); return; }
352
353 AliMCVertex *mcpv = (AliMCVertex *) mcEvent->GetPrimaryVertex();
354 Double_t mcpvPos[3];
355 if (mcpv != 0)
356 {
357 mcpv->GetXYZ(mcpvPos);
358 }
359 else
360 {
361 Printf("ERROR: Could not resolve MC primary vertex");
362 return;
363 }
364
365 //loop over all MC tracks
366 for(Int_t iMCtrack = 0; iMCtrack < mcEvent->GetNumberOfTracks(); iMCtrack++)
367 {
368
369 //booleans to check if track is La, Lb, K0 and primary
370 bool lambdaMC = false;
371 bool antilambdaMC = false;
372 bool kshortMC = false;
373 bool isprimaryMC = false;
374
375 AliMCParticle *mcPart = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(iMCtrack));
376
377 if(mcPart->PdgCode() == kLambda0)
378 {
379 lambdaMC = true;
380 nLambdaMC++;
381 }
382 else if(mcPart->PdgCode() == kK0Short)
383 {
384 kshortMC = true;
385 nKshortMC++;
386 }
387 else if(mcPart->PdgCode() == kLambda0Bar)
388 {
389 antilambdaMC = true;
390 nAntilambdaMC++;
391 }
392
45b90328 393
45b90328 394 Int_t motherLabel = mcPart->GetMother();
45b90328 395 AliMCParticle *mcMother = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(motherLabel));
45b90328 396 double motherType = -1;
397 if(motherLabel >= 0)
398 {motherType = mcMother->PdgCode();}
399
889bc324 400 // this block of code is used to include primary Sigma0 decays as primary lambda/antilambda
45b90328 401 bool sigma0MC = false;
402 if(motherType == 3212 || motherType == -3212)
403 {
404 if(mcEvent->IsPhysicalPrimary(motherLabel))
405 {sigma0MC = true;}
406 }
407
408
409 if(mcEvent->IsPhysicalPrimary(iMCtrack) || sigma0MC)
410 {
411 isprimaryMC = true;
412 if(lambdaMC)
413 {
414 fHistMcNLambdaPrimary->Fill(1);
45b90328 415 if(TMath::Abs(mcPart->Y())<=0.5)
889bc324 416 {fHistMcPMLaPt->Fill(mcPart->Pt(),mcPart->M());}
45b90328 417 }
418 if(antilambdaMC)
419 {
420 if(TMath::Abs(mcPart->Y())<=0.5)
421 {fHistMcPMLbPt->Fill(mcPart->Pt(),mcPart->M());}
422 }
423 if(kshortMC)
424 {
425 if(TMath::Abs(mcPart->Y())<=0.5)
426 {fHistMcPMK0Pt->Fill(mcPart->Pt(),mcPart->M());}
427 }
428 }
429 else
430 {
431 isprimaryMC = false;
432 if(lambdaMC)
433 {
434 fHistMcNLambdaPrimary->Fill(2);
45b90328 435 }
436 }
437
438 }
439
440
441 }
442
443
444 fHistMcNLambda->Fill(nLambdaMC);
445 fHistMcNAntilambda->Fill(nAntilambdaMC);
446 fHistMcNKshort->Fill(nKshortMC);
447
448 //END OF MONTE CARLO SECTION
449 /*********************************************************************/
450
451
452
453
454
455 // Do some fast cuts first
456 // check for good reconstructed vertex
457 if(!(fESD->GetPrimaryVertex()->GetStatus())) return;
458 // if vertex is from spd vertexZ, require more stringent cut
459 if (fESD->GetPrimaryVertex()->IsFromVertexerZ()) {
460 if (fESD->GetPrimaryVertex()->GetDispersion()>0.02 || fESD->GetPrimaryVertex()->GetZRes()>0.25 ) return; // bad vertex from VertexerZ
461 }
462
463 Double_t tV0Position[3];
464 Double_t tPVPosition[3];
465 Double_t radius;
466
467 // physics selection
468 UInt_t maskIsSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
469 if(!maskIsSelected)
470 {
471 //printf("Event failed physics selection\n");
472 return;
473 }
474
475 //if additional initial cuts wanted, can set conditions here
476 bool isCut = (fESD->GetNumberOfTracks()==0);
477 if (isCut)
478 {return;}
479
480 //gets primary vertex for the event
481 const AliESDVertex *kPv = ((AliESDEvent *)fESD)->GetPrimaryVertex();
482 if ( kPv != 0 )
483 {
484 tPVPosition[0] = kPv->GetXv();
485 tPVPosition[1] = kPv->GetYv();
486 tPVPosition[2] = kPv->GetZv();
487 if( tPVPosition[2] == 0. )
488 {
489 //printf("WARNING: Primary vertex a Z = 0, aborting\n");
490 return;
491 }
492 }
493 else
494 {
495 //printf("ERROR: Primary vertex not found\n");
496 return;
497 }
498 if( !kPv->GetStatus())
499 {return;}
500
501
502
503
504 int nLambda(0);
505 int nV0s(0);
45b90328 506
507 // V0 loop - runs over every v0 in the event
508 for (Int_t iV0 = 0; iV0 < fESD->GetNumberOfV0s(); iV0++)
509 {
510
511 AliESDv0 *v0 = fESD->GetV0(iV0);
512 if (!v0)
513 {
514 printf("ERROR: Could not receive v0 %d\n", iV0);
515 continue;
516 }
517
518 bool lambdaCandidate = true;
519 bool antilambdaCandidate = true;
520 bool kshortCandidate = true;
521
522 // keep only events of interest for fHistMLa plots
523
524 if (v0->GetEffMass(4,2) < 1.08 || v0->GetEffMass(4,2) > 1.2 || TMath::Abs(v0->Y(3122))>0.5 )
525 {lambdaCandidate = false;}
526 if (v0->GetEffMass(2,4) < 1.08 || v0->GetEffMass(2,4) > 1.2 || TMath::Abs(v0->Y(-3122))>0.5)
527 {antilambdaCandidate = false;}
528 if (v0->GetEffMass(2,2) < 0.414 || v0->GetEffMass(2,2) > 0.582 || TMath::Abs(v0->Y(310))>0.5)
529 {kshortCandidate = false;}
530 if (v0->GetOnFlyStatus())
531 {continue;}
532
533 if(!isMonteCarlo)
534 {if(lambdaCandidate == false && antilambdaCandidate == false && kshortCandidate == false)
535 {continue;}}
536
537
538 //gets details of the v0
539 v0->GetXYZ(tV0Position[0],tV0Position[1],tV0Position[2]);
540 radius = TMath::Sqrt(tV0Position[0]*tV0Position[0]+tV0Position[1]*tV0Position[1]);
541
542 double decayLength = (sqrt((tV0Position[0]-tPVPosition[0])*(tV0Position[0]-tPVPosition[0])+(tV0Position[1]-tPVPosition[1])*(tV0Position[1]-tPVPosition[1])+(tV0Position[2]-tPVPosition[2])*(tV0Position[2]-tPVPosition[2])));
543 double cTauLa = decayLength*(v0->GetEffMass(4,2))/(v0->P());
544 double cTauLb = decayLength*(v0->GetEffMass(2,4))/(v0->P());
545 double cTauK0 = decayLength*(v0->GetEffMass(2,2))/(v0->P());
546
547 Int_t indexP, indexN;
548 indexP = TMath::Abs(v0->GetPindex());
549 AliESDtrack *posTrack = ((AliESDEvent*)fESD)->GetTrack(indexP);
550 indexN = TMath::Abs(v0->GetNindex());
551 AliESDtrack *negTrack = ((AliESDEvent*)fESD)->GetTrack(indexN);
552
553 if(!posTrack || !negTrack)
554 {continue;}
555
556 double pTrackMomentum[3];
557 double nTrackMomentum[3];
558 double pV0x, pV0y, pV0z;
559 posTrack->GetConstrainedPxPyPz(pTrackMomentum);
560 negTrack->GetConstrainedPxPyPz(nTrackMomentum);
561 v0->GetPxPyPz(pV0x, pV0y, pV0z);
562
889bc324 563 //const double kMLambda = 1.115;
45b90328 564 const double kMProton = 0.938;
889bc324 565 //const double kMPi = 0.140;
45b90328 566
567 double pPos2 = sqrt(pTrackMomentum[0]*pTrackMomentum[0]+pTrackMomentum[1]*pTrackMomentum[1]+pTrackMomentum[2]*pTrackMomentum[2]);
568 double pNeg2 = sqrt(nTrackMomentum[0]*nTrackMomentum[0]+nTrackMomentum[1]*nTrackMomentum[1]+nTrackMomentum[2]*nTrackMomentum[2]);
889bc324 569 //double pV02 = sqrt(pV0x*pV0x+pV0y*pV0y+pV0z*pV0z);
45b90328 570
571 //to prevent segfaults when ratios etc taken
889bc324 572 //if(pV02 < 0.01 || pPos2 <0.01 || pNeg2 <0.01)
573 //{continue;}
45b90328 574
575 Float_t pImpactxy(0), pImpactz(0);
576 Float_t nImpactxy(0), nImpactz(0);
577 posTrack->GetImpactParameters(pImpactxy,pImpactz);
578 negTrack->GetImpactParameters(nImpactxy,nImpactz);
579 nImpactxy = sqrt((nImpactxy*nImpactxy));
580 nImpactz = sqrt((nImpactz *nImpactz ));
581 pImpactxy = sqrt((pImpactxy*pImpactxy));
582 pImpactz = sqrt((pImpactz *pImpactz ));
583
584 /*********************************************************************/
585 // Cuts are implemented here.
586
587 if(!(fTrackCuts->IsSelected(posTrack)) || !(fTrackCuts->IsSelected(negTrack)))
588 {
589 lambdaCandidate = false;
590 antilambdaCandidate = false;
591 kshortCandidate = false;
592 }
593
594 //extra cut to account for difference between p2 & p1 data
595 if(nImpactxy < 0.1 || pImpactxy < 0.1)
596 {
597 lambdaCandidate = false;
598 antilambdaCandidate = false;
599 kshortCandidate = false;
600 }
601
602 //psuedorapidity cut
603 if(cutEta != -999)
604 {
605 if(TMath::Abs(posTrack->Eta()) > cutEta || TMath::Abs(negTrack->Eta()) >cutEta)
606 {
607 lambdaCandidate = false;
608 antilambdaCandidate = false;
609 kshortCandidate = false;
610 }
611 }
612
613 //pointing angle cut
614 if(cutCosPa != -999)
615 {
616 if (v0->GetV0CosineOfPointingAngle(tPVPosition[0],tPVPosition[1],tPVPosition[2]) < cutCosPa)
617 {
618 lambdaCandidate = false;
619 antilambdaCandidate = false;
620 kshortCandidate = false;
621 }
622 }
623
624 //lifetime cut
625 if(cutcTau != -999)
626 {
627 if(cTauLa < cutcTau)
628 {
629 lambdaCandidate = false;
630 }
631 if(cTauLb < cutcTau)
632 {
633 antilambdaCandidate = false;
634 }
635 if(cTauK0 < cutcTau)
636 {
637 kshortCandidate = false;
638 }
639
640 }
641
642 // Impact paramater cut (on neg particle)
643 if(cutNImpact != -999)
644 {
645 if(nImpactxy < cutNImpact || nImpactz < cutNImpact)
646 {
647 lambdaCandidate = false;
648 }
649 if(pImpactxy < cutNImpact || pImpactz < cutNImpact)
650 {
651 antilambdaCandidate = false;
652 }
653 }
654
655
656 // DCA between daughterscut
657 if(cutDCA != -999)
658 {
659 if(v0->GetDcaV0Daughters() > cutDCA)
660 {
661 lambdaCandidate = false;
662 antilambdaCandidate = false;
663 kshortCandidate = false;
664 }
665 }
666
667 // Bethe Bloch cut. Made sightly complicated as options for crude cuts still included. Should probably reduce to just 'official' cuts
668 if(cutBetheBloch != -999)
669 {
670 if(posTrack->GetTPCsignal() <0 || negTrack->GetTPCsignal()<0)
671 {continue;}
672
673 if(lambdaCandidate)
674 {
675 if(cutBetheBloch > 0)
676 {
677 if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(posTrack, AliPID::kProton)) > cutBetheBloch )
678 {lambdaCandidate = false;}
679 if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(negTrack, AliPID::kPion)) > cutBetheBloch )
680 {lambdaCandidate = false;}
681 }
682
683 if(cutBetheBloch == -4)
684 {
685 if(isMonteCarlo)
686 {
687 double beta2 = TMath::Power((pPos2/TMath::Sqrt((pPos2*pPos2+0.9*0.9))),2);
688 double gamma2 = 1.0/(1.0-beta2);
689 if(posTrack->GetTPCsignal() < (2.0/beta2)*(TMath::Log(1e6*beta2*gamma2)-beta2))
690 {lambdaCandidate = false;}
691 }
692 else
693 {
694 double beta2 = TMath::Power((pPos2/TMath::Sqrt((pPos2*pPos2+kMProton*kMProton))),2);
695 double gamma2 = 1.0/(1.0-beta2);
696 if(posTrack->GetTPCsignal() < (2.3/beta2)*(TMath::Log(1e6*beta2*gamma2)-beta2))
697 {lambdaCandidate = false;}
698 }
699 }
700
701 }
702
703 if(antilambdaCandidate)
704 {
705 if(cutBetheBloch > 0)
706 {
707 if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(negTrack, AliPID::kProton)) > cutBetheBloch )
708 {antilambdaCandidate = false;}
709 if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(posTrack, AliPID::kPion)) > cutBetheBloch )
710 {antilambdaCandidate = false;}
711 }
712
713 if(cutBetheBloch == -4)
714 {
715 if(isMonteCarlo)
716 {
717 double beta2 = TMath::Power((pNeg2/TMath::Sqrt((pNeg2*pNeg2+0.9*0.9))),2);
718 double gamma2 = 1.0/(1.0-beta2);
719 if(negTrack->GetTPCsignal() < (2.0/beta2)*(TMath::Log(1e6*beta2*gamma2)-beta2))
720 {antilambdaCandidate = false;}
721 }
722 else
723 {
724 double beta2 = TMath::Power((pNeg2/TMath::Sqrt((pNeg2*pNeg2+0.9*0.9))),2);
725 double gamma2 = 1.0/(1.0-beta2);
726 if(negTrack->GetTPCsignal() < (2.3/beta2)*(TMath::Log(1e6*beta2*gamma2)-beta2))
727 {antilambdaCandidate = false;}
728 }
729 }
730
731 }
732
733 if(kshortCandidate)
734 {
735 if(cutBetheBloch > 0)
736 {
737 if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(negTrack, AliPID::kPion)) > cutBetheBloch )
738 {kshortCandidate = false;}
739 if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(posTrack, AliPID::kPion)) > cutBetheBloch )
740 {kshortCandidate = false;}
741 }
742
743
744 if(cutBetheBloch == -4)
745 {
746 double par0 = 0.20;
747 double par1 = 4.2;
748 double par2 = 1000000;
749
750 if(isMonteCarlo)
751 {
752 double beta2 = TMath::Power((pNeg2/TMath::Sqrt((pNeg2*pNeg2+par0*par0))),2);
753 double gamma2 = 1.0/(1.0-beta2);
754 if(negTrack->GetTPCsignal() > (par1/beta2)*(TMath::Log(par2*beta2*gamma2)-beta2) && TMath::Log10(pNeg2) > -0.6)
755 {kshortCandidate = false;}
756
757 beta2 = TMath::Power((pPos2/TMath::Sqrt((pPos2*pPos2+par0*par0))),2);
758 gamma2 = 1.0/(1.0-beta2);
759 if(posTrack->GetTPCsignal() > (par1/beta2)*(TMath::Log(par2*beta2*gamma2)-beta2) && TMath::Log10(pNeg2) > -0.6)
760 {kshortCandidate = false;}
761 }
762 else
763 {
764 double beta2 = TMath::Power((pNeg2/TMath::Sqrt((pNeg2*pNeg2+par0*par0))),2);
765 double gamma2 = 1.0/(1.0-beta2);
766 if(negTrack->GetTPCsignal() > (par1/beta2)*(TMath::Log(par2*beta2*gamma2)-beta2) && TMath::Log10(pNeg2) > -0.6)
767 {kshortCandidate = false;}
768
769 beta2 = TMath::Power((pPos2/TMath::Sqrt((pPos2*pPos2+par0*par0))),2);
770 gamma2 = 1.0/(1.0-beta2);
771 if(posTrack->GetTPCsignal() > (par1/beta2)*(TMath::Log(par2*beta2*gamma2)-beta2) && TMath::Log10(pPos2) > -0.6)
772 {kshortCandidate = false;}
773 }
774 }
775
776 }
777 }
778
779 // Selection of random cuts which I've been playing with
780 /*if(nImpactxy > 3 || pImpactxy > 2)
781 {
782 lambdaCandidate = false;
783 }*/
784
785 /*if(nImpactxy < 0.4 || pImpactxy < 0.4 || nImpactxy > 2.5 || pImpactxy >2.5)
786 {
787 antilambdaCandidate = false;
788 } */
789
889bc324 790 if(decayLength > 15 )
791 {lambdaCandidate = false;}
45b90328 792
793 // Cuts to look at just lowpt region of lambdas
889bc324 794 //if(v0->Pt() < 0.3 || v0->Pt() > 0.7 || !lambdaCandidate)
795 //{continue;}
45b90328 796
797 // cuts to just look at signal/background region of lambda mass
798 //if(!((v0->GetEffMass(4,2) >= 1.096 && v0->GetEffMass(4,2) < 1.106) || (v0->GetEffMass(4,2) >= 1.126 && v0->GetEffMass(4,2) < 1.136) ))
799 //if(!(v0->GetEffMass(4,2) >= 1.106 && v0->GetEffMass(4,2) < 1.126 ))
800 //{continue;}
801 /*********************************************************************/
802
803 /*********************************************************************/
804 // MONTE CARLO SECTION 2
805 // this section looks at the individual V0s
806
807 bool mcLambdaCandidate = true;
808 bool mcAntilambdaCandidate = true;
809 bool mcK0Candidate = true;
45b90328 810 bool realParticle = true;
811
812 if(isMonteCarlo)
813 {
814
815 AliMCEvent *mcEvent = MCEvent();
816 AliStack* mcStack = mcEvent->Stack();
817 if( !mcStack ) { Printf( "Stack not available"); return; }
818
819 TParticle *negParticle = mcStack->Particle( TMath::Abs(negTrack->GetLabel()));
820 TParticle *posParticle = mcStack->Particle( TMath::Abs(posTrack->GetLabel()));
821
822 Int_t negParticleMotherLabel = negParticle->GetFirstMother();
823 Int_t posParticleMotherLabel = posParticle->GetFirstMother();
824
825 if( negParticleMotherLabel == -1 || posParticleMotherLabel == -1)
826 {
827 realParticle = false;
828 mcLambdaCandidate = false;
829 mcAntilambdaCandidate = false;
830 mcK0Candidate =false;
831 }
832
45b90328 833 if( negParticleMotherLabel != posParticleMotherLabel)
834 {
835 mcLambdaCandidate = false;
836 mcAntilambdaCandidate = false;
837 mcK0Candidate =false;
838 }
839
840 if(realParticle == true)
841 {
842 AliMCParticle *mcPart2 = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(negParticleMotherLabel));
889bc324 843
45b90328 844 if(mcPart2->PdgCode() != kLambda0)
845 {mcLambdaCandidate = false;}
846 if(mcPart2->PdgCode() != kLambda0Bar)
847 {mcAntilambdaCandidate = false;}
848 if(mcPart2->PdgCode() != kK0Short)
849 {mcK0Candidate =false;}
889bc324 850
45b90328 851 if(mcLambdaCandidate && lambdaCandidate)
852 {
45b90328 853 fHistMcV0MLaPt->Fill(v0->Pt(),v0->GetEffMass(4,2));
854 }
855 if(mcAntilambdaCandidate && antilambdaCandidate)
856 {
45b90328 857 fHistMcV0MLbPt->Fill(v0->Pt(),v0->GetEffMass(2,4));
858 }
859 if(mcK0Candidate && kshortCandidate)
860 {
45b90328 861 fHistMcV0MK0Pt->Fill(v0->Pt(),v0->GetEffMass(2,2));
862 }
45b90328 863 }
864
889bc324 865 if(mcLambdaCandidate && lambdaCandidate)
45b90328 866 {
889bc324 867 fHistMCLambdacTau->Fill(cTauLa);
868 fHistMCLambdaDecayL->Fill(decayLength);
45b90328 869 }
889bc324 870 if(!mcLambdaCandidate && lambdaCandidate)
45b90328 871 {
889bc324 872 fHistMCLambdaNotcTau->Fill(cTauLa);
873 fHistMCLambdaNotDecayL->Fill(decayLength);
45b90328 874 }
889bc324 875
876 }
877
45b90328 878
879 // END OF MONTE CARLO SECTION 2
880 /*********************************************************************/
881
882 //remove all non-candidates
883 if(lambdaCandidate == false && antilambdaCandidate == false && kshortCandidate == false)
884 {continue;}
885
886
887 //count number of valid v0s
888 nV0s+=1;
889
890 /*********************************************************************/
891 //This section fills histograms
892
45b90328 893 fHistDCAV0Daughters->Fill(v0->GetDcaV0Daughters());
45b90328 894 fHistCosPA->Fill(v0->GetV0CosineOfPointingAngle(tPVPosition[0],tPVPosition[1],tPVPosition[2]));
895 fHistDecayL->Fill(decayLength);
45b90328 896 fHistTauLa->Fill(cTauLa);
45b90328 897
889bc324 898 fHistBetheBlochTPCPos->Fill(TMath::Log10(pPos2),posTrack->GetTPCsignal());
899 fHistBetheBlochTPCNeg->Fill(TMath::Log10(pNeg2),negTrack->GetTPCsignal());
45b90328 900
45b90328 901 fHistImpactxyN->Fill(nImpactxy);
902 fHistImpactzN->Fill(nImpactz);
903 fHistImpactxyP->Fill(pImpactxy);
904 fHistImpactzP->Fill(pImpactz);
905
45b90328 906 fHistImpactxyImpactz->Fill(nImpactxy,nImpactz);
907
908 fHistV0Z->Fill(tV0Position[2]);
909 fHistZ->Fill(tV0Position[2]-tPVPosition[2]);
889bc324 910
45b90328 911 fHistRZ->Fill(tV0Position[2],radius);
912 fHistPtV0Z->Fill(v0->Pt(),tV0Position[2]);
913
914 fHistPtArm->Fill(v0->AlphaV0(),v0->PtArmV0());
915 fHistXZ->Fill(tV0Position[2],tV0Position[0]);
916 fHistYZ->Fill(tV0Position[2],tV0Position[1]);
45b90328 917 fHistPtV0->Fill(v0->Pt());
918
919 //effective mass histograms
920
921 //sets assumed particle type of pos/neg daughters.
922 // 0 = electron, 1 = Muon, 2 = pion, 3 = kaon, 4 = proton.
923 int dPos = 0;
924 int dNeg = 0;
925
926 // v0->ChangeMassHypothesis(kLambda0);
927 dPos = 4;
928 dNeg = 2;
929 if(v0->GetEffMass(dPos,dNeg) > 1.11 && v0->GetEffMass(dPos,dNeg) < 1.13)
930 {
931 if(!(v0->GetOnFlyStatus()))
932 {
933 nLambda++;
934 }
935 }
936 if(lambdaCandidate)
937 {
889bc324 938 fHistMLa->Fill(v0->GetEffMass(dPos,dNeg));
45b90328 939 fHistMLaPt->Fill(v0->Pt(),v0->GetEffMass(dPos,dNeg));
940 }
941
942 // v0->ChangeMassHypothesis(kK0Short);
943 dPos = 2;
944 dNeg = 2;
945 if(kshortCandidate)
946 {
947 fHistMK0->Fill(v0->GetEffMass(dPos,dNeg));
45b90328 948 fHistMK0Pt->Fill(v0->Pt(),v0->GetEffMass(dPos,dNeg));
949 }
950 // v0->ChangeMassHypothesis(kLambda0Bar);
951 dPos = 2;
952 dNeg = 4;
953 if(antilambdaCandidate)
954 {
955 fHistMLb->Fill(v0->GetEffMass(dPos,dNeg));
45b90328 956 fHistMLbPt->Fill(v0->Pt(),v0->GetEffMass(dPos,dNeg));
957 }
889bc324 958
45b90328 959 }
960
45b90328 961
962 fHistPVZ->Fill(tPVPosition[2]);
963 fHistNV0->Fill(nV0s);
45b90328 964 fHistNLambda->Fill(nLambda);
965
889bc324 966 // NEW HISTO should be filled before this point, as PostData puts the
45b90328 967 // information for this iteration of the UserExec in the container
968 PostData(1, fOutputList);
969}
970
971
972//________________________________________________________________________
973void AliAnalysisTaskLukeV0::Terminate(Option_t *)
974{
975 // Draw result to screen, or perform fitting, normalizations
976 // Called once at the end of the query
977
978 fOutputList = dynamic_cast<TList*> (GetOutputData(1));
979 if(!fOutputList) { Printf("ERROR: could not retrieve TList fOutputList"); return; }
980
981 }