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