Coverity
[u/mrichter/AliRoot.git] / PWGLF / STRANGENESS / LambdaK0PbPb / AliAnalysisTaskLukeV0.cxx
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
25  * ------------------------*
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
53 ClassImp(AliAnalysisTaskLukeV0)
54
55 //________________________________________________________________________
56 AliAnalysisTaskLukeV0::AliAnalysisTaskLukeV0() // All data members should be initialised here
57    :AliAnalysisTaskSE(),
58     fOutputList(0),
59     fTrackCuts(0),
60         fPIDResponse(0),
61         fHistCosPA(0), 
62         fHistDCAV0Daughters(0), 
63         fHistDecayL(0), 
64         fHistImpactxyN(0), 
65         fHistImpactzN(0), 
66         fHistImpactxyP(0), 
67         fHistImpactzP(0), 
68         fHistMCLambdacTau(0),
69         fHistMCLambdaNotcTau(0),
70         fHistMCLambdaDecayL(0),
71         fHistMCLambdaNotDecayL(0),
72         fHistMcNLambdaPrimary(0),
73         fHistMcNLambda(0),
74         fHistMcNAntilambda(0),
75         fHistMcNKshort(0),
76         fHistMK0(0), 
77         fHistMLa(0), 
78         fHistMLb(0), 
79         fHistNLambda(0), 
80         fHistNV0(0), 
81         fHistPtV0(0), 
82         fHistPVZ(0), 
83         fHistTauLa(0), 
84         fHistV0Z(0), 
85         fHistZ(0),
86         fHistBetheBlochTPCNeg(0), 
87         fHistBetheBlochTPCPos(0), 
88         fHistImpactxyImpactz(0), 
89         fHistMcPMK0Pt(0),
90         fHistMcPMLaPt(0),
91         fHistMcPMLbPt(0),
92         fHistMcV0MK0Pt(0),
93         fHistMcV0MLaPt(0),
94         fHistMcV0MLbPt(0),
95         fHistMK0Pt(0), 
96         fHistMLaPt(0), 
97         fHistMLbPt(0), 
98         fHistPtArm(0), 
99         fHistPtV0Z(0),
100         fHistRZ(0), 
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 //________________________________________________________________________
108 AliAnalysisTaskLukeV0::AliAnalysisTaskLukeV0(const char *name) // All data members should be initialised here
109    :AliAnalysisTaskSE(name),
110         fOutputList(0),
111         fTrackCuts(0),
112         fPIDResponse(0),
113         fHistCosPA(0), 
114         fHistDCAV0Daughters(0), 
115         fHistDecayL(0), 
116         fHistImpactxyN(0), 
117         fHistImpactzN(0), 
118         fHistImpactxyP(0), 
119         fHistImpactzP(0), 
120         fHistMCLambdacTau(0),
121         fHistMCLambdaNotcTau(0),
122         fHistMCLambdaDecayL(0),
123         fHistMCLambdaNotDecayL(0),
124         fHistMcNLambdaPrimary(0),
125         fHistMcNLambda(0),
126         fHistMcNAntilambda(0),
127         fHistMcNKshort(0),
128         fHistMK0(0), 
129         fHistMLa(0), 
130         fHistMLb(0), 
131         fHistNLambda(0), 
132         fHistNV0(0), 
133         fHistPtV0(0), 
134         fHistPVZ(0), 
135         fHistTauLa(0), 
136         fHistV0Z(0), 
137         fHistZ(0),
138         fHistBetheBlochTPCNeg(0), 
139         fHistBetheBlochTPCPos(0), 
140         fHistImpactxyImpactz(0), 
141         fHistMcPMK0Pt(0),
142         fHistMcPMLaPt(0),
143         fHistMcPMLbPt(0),
144         fHistMcV0MK0Pt(0),
145         fHistMcV0MLaPt(0),
146         fHistMcV0MLbPt(0),
147         fHistMK0Pt(0), 
148         fHistMLaPt(0), 
149         fHistMLbPt(0), 
150         fHistPtArm(0),
151         fHistPtV0Z(0),
152         fHistRZ(0), 
153         fHistXZ(0), 
154         fHistYZ(0) // The last in the above list should not have a comma after it
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 //________________________________________________________________________
164 AliAnalysisTaskLukeV0::~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 //________________________________________________________________________
175 void 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
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
196         fHistCosPA = new        TH1F("fHistCosPA", "Cosine of Pointing Angle of V0s; Cos PA; N(v0s)",202,0.8,1.01);
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);
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);
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);
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);
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);
212         fHistMLa = new  TH1F("fHistMLa","Lambda Mass; M(p#pi^{-}) (GeV/c^{2}); dN/dM (0.125 GeV/c^{2})^{-1}",div,min,max);
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);
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);
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);
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
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);
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);
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);
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);
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);
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);
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);
235         fHistPtArm = new        TH2F("fHistPtArm","Podolanski-Armenteros Plot; #alpha; P_{perp} (GeV/c)",40,-1,1,80,0,0.5);
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);
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         
245         fOutputList->Add(fHistCosPA); 
246         fOutputList->Add(fHistDCAV0Daughters); 
247         fOutputList->Add(fHistDecayL); 
248         fOutputList->Add(fHistImpactxyN); 
249         fOutputList->Add(fHistImpactzN); 
250         fOutputList->Add(fHistImpactxyP); 
251         fOutputList->Add(fHistImpactzP); 
252         fOutputList->Add(fHistMCLambdacTau);
253         fOutputList->Add(fHistMCLambdaNotcTau);
254         fOutputList->Add(fHistMCLambdaDecayL);
255         fOutputList->Add(fHistMCLambdaNotDecayL);
256         fOutputList->Add(fHistMcNLambdaPrimary);
257         fOutputList->Add(fHistMcNLambda);
258         fOutputList->Add(fHistMcNAntilambda);
259         fOutputList->Add(fHistMcNKshort);
260         fOutputList->Add(fHistMK0); 
261         fOutputList->Add(fHistMLa); 
262         fOutputList->Add(fHistMLb); 
263         fOutputList->Add(fHistNLambda); 
264         fOutputList->Add(fHistNV0); 
265         fOutputList->Add(fHistPtV0); 
266         fOutputList->Add(fHistPVZ); 
267         fOutputList->Add(fHistTauLa); 
268         fOutputList->Add(fHistV0Z); 
269         fOutputList->Add(fHistZ);
270         fOutputList->Add(fHistBetheBlochTPCNeg); 
271         fOutputList->Add(fHistBetheBlochTPCPos); 
272         fOutputList->Add(fHistImpactxyImpactz); 
273         fOutputList->Add(fHistMcPMK0Pt);
274         fOutputList->Add(fHistMcPMLaPt);
275         fOutputList->Add(fHistMcPMLbPt);
276         fOutputList->Add(fHistMcV0MK0Pt);
277         fOutputList->Add(fHistMcV0MLaPt);
278         fOutputList->Add(fHistMcV0MLbPt);
279         fOutputList->Add(fHistMK0Pt); 
280         fOutputList->Add(fHistMLaPt); 
281         fOutputList->Add(fHistMLbPt); 
282         fOutputList->Add(fHistPtArm); 
283         fOutputList->Add(fHistRZ); 
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 //________________________________________________________________________
292 void 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);
302         double isMonteCarlo(true);
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"); 
344                                 return; 
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));
378                         if (!mcPart) continue;
379
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                         
396                         
397                         Int_t motherLabel = mcPart->GetMother();
398                         AliMCParticle *mcMother = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(motherLabel));
399                         double motherType = -1;
400                         if(motherLabel >= 0 && mcMother)
401                         {motherType = mcMother->PdgCode();}
402                         
403                         // this block of code is used to include primary Sigma0 decays as primary lambda/antilambda
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);
418                                            if(TMath::Abs(mcPart->Y())<=0.5)
419                                                 {fHistMcPMLaPt->Fill(mcPart->Pt(),mcPart->M());}
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);
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     {
487                 tPVPosition[0] = kPv->GetXv();
488                 tPVPosition[1] = kPv->GetYv();
489                 tPVPosition[2] = kPv->GetZv();
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);
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
566                 //const double kMLambda = 1.115;
567                 const double kMProton = 0.938;
568                 //const double kMPi     = 0.140;
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]);
572                 //double pV02 = sqrt(pV0x*pV0x+pV0y*pV0y+pV0z*pV0z);
573                 
574                 //to prevent segfaults when ratios etc taken
575                 //if(pV02 < 0.01 || pPos2 <0.01 || pNeg2 <0.01)
576                 //{continue;}
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                 
793                 if(decayLength > 15 )
794                 {lambdaCandidate = false;}
795                 
796                 // Cuts to look at just lowpt region of lambdas
797                 //if(v0->Pt() < 0.3 || v0->Pt() > 0.7 || !lambdaCandidate)
798                 //{continue;}
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;
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
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));
846                                 if (mcPart2) {
847
848                                 if(mcPart2->PdgCode() != kLambda0)
849                                         {mcLambdaCandidate = false;}
850                                 if(mcPart2->PdgCode() != kLambda0Bar)
851                                         {mcAntilambdaCandidate = false;}
852                                 if(mcPart2->PdgCode() != kK0Short)
853                                         {mcK0Candidate  =false;}
854
855                                 if(mcLambdaCandidate && lambdaCandidate)
856                                 {
857                                         fHistMcV0MLaPt->Fill(v0->Pt(),v0->GetEffMass(4,2)); 
858                                 }
859                                 if(mcAntilambdaCandidate && antilambdaCandidate)
860                                 {
861                                         fHistMcV0MLbPt->Fill(v0->Pt(),v0->GetEffMass(2,4));
862                                 }
863                                 if(mcK0Candidate && kshortCandidate)
864                                 {
865                                         fHistMcV0MK0Pt->Fill(v0->Pt(),v0->GetEffMass(2,2));
866                                 }
867                                 }
868                         }
869
870                         if(mcLambdaCandidate && lambdaCandidate)
871                         {
872                                 fHistMCLambdacTau->Fill(cTauLa);
873                                 fHistMCLambdaDecayL->Fill(decayLength);
874                         }
875                         if(!mcLambdaCandidate && lambdaCandidate)
876                         {
877                         fHistMCLambdaNotcTau->Fill(cTauLa);
878                         fHistMCLambdaNotDecayL->Fill(decayLength);
879                         }
880                                         
881                         }
882                 
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                 
898                 fHistDCAV0Daughters->Fill(v0->GetDcaV0Daughters());
899                 fHistCosPA->Fill(v0->GetV0CosineOfPointingAngle(tPVPosition[0],tPVPosition[1],tPVPosition[2]));
900                 fHistDecayL->Fill(decayLength);
901                 fHistTauLa->Fill(cTauLa);
902                 
903                 fHistBetheBlochTPCPos->Fill(TMath::Log10(pPos2),posTrack->GetTPCsignal());
904                 fHistBetheBlochTPCNeg->Fill(TMath::Log10(pNeg2),negTrack->GetTPCsignal());
905
906                 fHistImpactxyN->Fill(nImpactxy);
907                 fHistImpactzN->Fill(nImpactz);
908                 fHistImpactxyP->Fill(pImpactxy);
909                 fHistImpactzP->Fill(pImpactz);
910                 
911                 fHistImpactxyImpactz->Fill(nImpactxy,nImpactz);
912                 
913                 fHistV0Z->Fill(tV0Position[2]);
914                 fHistZ->Fill(tV0Position[2]-tPVPosition[2]);
915         
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]);
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                 {
943                         fHistMLa->Fill(v0->GetEffMass(dPos,dNeg));  
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));
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));
961                         fHistMLbPt->Fill(v0->Pt(),v0->GetEffMass(dPos,dNeg));
962                 }
963                                 
964     }
965         
966         
967         fHistPVZ->Fill(tPVPosition[2]);
968         fHistNV0->Fill(nV0s);
969         fHistNLambda->Fill(nLambda);    
970         
971         // NEW HISTO should be filled before this point, as PostData puts the
972     // information for this iteration of the UserExec in the container
973         PostData(1, fOutputList);
974 }
975
976
977 //________________________________________________________________________
978 void 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  }