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         fHistRZ(0),
 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 //________________________________________________________________________
 107 AliAnalysisTaskLukeV0::AliAnalysisTaskLukeV0(const char *name) // All data members should be initialised here
 108    :AliAnalysisTaskSE(name),
 109         fOutputList(0),
 110         fTrackCuts(0),
 111         fPIDResponse(0),
 112         fHistCosPA(0),
 113         fHistDCAV0Daughters(0),
 114         fHistDecayL(0),
 115         fHistImpactxyN(0),
 116         fHistImpactzN(0),
 117         fHistImpactxyP(0),
 118         fHistImpactzP(0),
 119         fHistMCLambdacTau(0),
 120         fHistMCLambdaNotcTau(0),
 121         fHistMCLambdaDecayL(0),
 122         fHistMCLambdaNotDecayL(0),
 123         fHistMcNLambdaPrimary(0),
 124         fHistMcNLambda(0),
 125         fHistMcNAntilambda(0),
 126         fHistMcNKshort(0),
 127         fHistMK0(0),
 128         fHistMLa(0),
 129         fHistMLb(0),
 130         fHistNLambda(0),
 131         fHistNV0(0),
 132         fHistPtV0(0),
 133         fHistPVZ(0),
 134         fHistTauLa(0),
 135         fHistV0Z(0),
 136         fHistZ(0),
 137         fHistBetheBlochTPCNeg(0),
 138         fHistBetheBlochTPCPos(0),
 139         fHistImpactxyImpactz(0),
 140         fHistMcPMK0Pt(0),
 141         fHistMcPMLaPt(0),
 142         fHistMcPMLbPt(0),
 143         fHistMcV0MK0Pt(0),
 144         fHistMcV0MLaPt(0),
 145         fHistMcV0MLbPt(0),
 146         fHistMK0Pt(0),
 147         fHistMLaPt(0),
 148         fHistMLbPt(0),
 149         fHistPtArm(0),
 150         fHistRZ(0),
 151         fHistXZ(0),
 152         fHistYZ(0) // The last in the above list should not have a comma after it
 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 //________________________________________________________________________
 162 AliAnalysisTaskLukeV0::~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 //________________________________________________________________________
 173 void 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
 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
 194         fHistCosPA = new        TH1F("fHistCosPA", "Cosine of Pointing Angle of V0s; Cos PA; N(v0s)",202,0.8,1.01);
 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);
 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);
 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);
 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);
 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);
 210         fHistMLa = new  TH1F("fHistMLa","Lambda Mass; M(p#pi^{-}) (GeV/c^{2}); dN/dM (0.125 GeV/c^{2})^{-1}",div,min,max);
 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);
 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);
 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);
 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
 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);
 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);
 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);
 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);
 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);
 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);
 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);
 233         fHistPtArm = new        TH2F("fHistPtArm","Podolanski-Armenteros Plot; #alpha; P_{perp} (GeV/c)",40,-1,1,80,0,0.5);
 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);
 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
 243         fOutputList->Add(fHistCosPA);
 244         fOutputList->Add(fHistDCAV0Daughters);
 245         fOutputList->Add(fHistDecayL);
 246         fOutputList->Add(fHistImpactxyN);
 247         fOutputList->Add(fHistImpactzN);
 248         fOutputList->Add(fHistImpactxyP);
 249         fOutputList->Add(fHistImpactzP);
 250         fOutputList->Add(fHistMCLambdacTau);
 251         fOutputList->Add(fHistMCLambdaNotcTau);
 252         fOutputList->Add(fHistMCLambdaDecayL);
 253         fOutputList->Add(fHistMCLambdaNotDecayL);
 254         fOutputList->Add(fHistMcNLambdaPrimary);
 255         fOutputList->Add(fHistMcNLambda);
 256         fOutputList->Add(fHistMcNAntilambda);
 257         fOutputList->Add(fHistMcNKshort);
 258         fOutputList->Add(fHistMK0);
 259         fOutputList->Add(fHistMLa);
 260         fOutputList->Add(fHistMLb);
 261         fOutputList->Add(fHistNLambda);
 262         fOutputList->Add(fHistNV0);
 263         fOutputList->Add(fHistPtV0);
 264         fOutputList->Add(fHistPVZ);
 265         fOutputList->Add(fHistTauLa);
 266         fOutputList->Add(fHistV0Z);
 267         fOutputList->Add(fHistZ);
 268         fOutputList->Add(fHistBetheBlochTPCNeg);
 269         fOutputList->Add(fHistBetheBlochTPCPos);
 270         fOutputList->Add(fHistImpactxyImpactz);
 271         fOutputList->Add(fHistMcPMK0Pt);
 272         fOutputList->Add(fHistMcPMLaPt);
 273         fOutputList->Add(fHistMcPMLbPt);
 274         fOutputList->Add(fHistMcV0MK0Pt);
 275         fOutputList->Add(fHistMcV0MLaPt);
 276         fOutputList->Add(fHistMcV0MLbPt);
 277         fOutputList->Add(fHistMK0Pt);
 278         fOutputList->Add(fHistMLaPt);
 279         fOutputList->Add(fHistMLbPt);
 280         fOutputList->Add(fHistPtArm);
 281         fOutputList->Add(fHistRZ);
 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 //________________________________________________________________________
 290 void 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);
 300         double isMonteCarlo(true);
 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
 393
 394                         Int_t motherLabel = mcPart->GetMother();
 395                         AliMCParticle *mcMother = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(motherLabel));
 396                         double motherType = -1;
 397                         if(motherLabel >= 0)
 398                         {motherType = mcMother->PdgCode();}
 399
 400                         // this block of code is used to include primary Sigma0 decays as primary lambda/antilambda
 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);
 415                                            if(TMath::Abs(mcPart->Y())<=0.5)
 416                                                 {fHistMcPMLaPt->Fill(mcPart->Pt(),mcPart->M());}
 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);
 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);
 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
 563                 //const double kMLambda = 1.115;
 564                 const double kMProton = 0.938;
 565                 //const double kMPi     = 0.140;
 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]);
 569                 //double pV02 = sqrt(pV0x*pV0x+pV0y*pV0y+pV0z*pV0z);
 570
 571                 //to prevent segfaults when ratios etc taken
 572                 //if(pV02 < 0.01 || pPos2 <0.01 || pNeg2 <0.01)
 573                 //{continue;}
 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
 790                 if(decayLength > 15 )
 791                 {lambdaCandidate = false;}
 792
 793                 // Cuts to look at just lowpt region of lambdas
 794                 //if(v0->Pt() < 0.3 || v0->Pt() > 0.7 || !lambdaCandidate)
 795                 //{continue;}
 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;
 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
 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));
 843
 844                                 if(mcPart2->PdgCode() != kLambda0)
 845                                         {mcLambdaCandidate = false;}
 846                                 if(mcPart2->PdgCode() != kLambda0Bar)
 847                                         {mcAntilambdaCandidate = false;}
 848                                 if(mcPart2->PdgCode() != kK0Short)
 849                                         {mcK0Candidate  =false;}
 850
 851                                 if(mcLambdaCandidate && lambdaCandidate)
 852                                 {
 853                                         fHistMcV0MLaPt->Fill(v0->Pt(),v0->GetEffMass(4,2));
 854                                 }
 855                                 if(mcAntilambdaCandidate && antilambdaCandidate)
 856                                 {
 857                                         fHistMcV0MLbPt->Fill(v0->Pt(),v0->GetEffMass(2,4));
 858                                 }
 859                                 if(mcK0Candidate && kshortCandidate)
 860                                 {
 861                                         fHistMcV0MK0Pt->Fill(v0->Pt(),v0->GetEffMass(2,2));
 862                                 }
 863                                 }
 864
 865                         if(mcLambdaCandidate && lambdaCandidate)
 866                         {
 867                                 fHistMCLambdacTau->Fill(cTauLa);
 868                                 fHistMCLambdaDecayL->Fill(decayLength);
 869                         }
 870                         if(!mcLambdaCandidate && lambdaCandidate)
 871                         {
 872                         fHistMCLambdaNotcTau->Fill(cTauLa);
 873                         fHistMCLambdaNotDecayL->Fill(decayLength);
 874                         }
 875
 876                         }
 877
 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
 893                 fHistDCAV0Daughters->Fill(v0->GetDcaV0Daughters());
 894                 fHistCosPA->Fill(v0->GetV0CosineOfPointingAngle(tPVPosition[0],tPVPosition[1],tPVPosition[2]));
 895                 fHistDecayL->Fill(decayLength);
 896                 fHistTauLa->Fill(cTauLa);
 897
 898                 fHistBetheBlochTPCPos->Fill(TMath::Log10(pPos2),posTrack->GetTPCsignal());
 899                 fHistBetheBlochTPCNeg->Fill(TMath::Log10(pNeg2),negTrack->GetTPCsignal());
 900
 901                 fHistImpactxyN->Fill(nImpactxy);
 902                 fHistImpactzN->Fill(nImpactz);
 903                 fHistImpactxyP->Fill(pImpactxy);
 904                 fHistImpactzP->Fill(pImpactz);
 905
 906                 fHistImpactxyImpactz->Fill(nImpactxy,nImpactz);
 907
 908                 fHistV0Z->Fill(tV0Position[2]);
 909                 fHistZ->Fill(tV0Position[2]-tPVPosition[2]);
 910
 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]);
 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                 {
 938                         fHistMLa->Fill(v0->GetEffMass(dPos,dNeg));
 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));
 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));
 956                         fHistMLbPt->Fill(v0->Pt(),v0->GetEffMass(dPos,dNeg));
 957                 }
 958
 959     }
 960
 961
 962         fHistPVZ->Fill(tPVPosition[2]);
 963         fHistNV0->Fill(nV0s);
 964         fHistNLambda->Fill(nLambda);
 965
 966         // NEW HISTO should be filled before this point, as PostData puts the
 967     // information for this iteration of the UserExec in the container
 968         PostData(1, fOutputList);
 969 }
 970
 971
 972 //________________________________________________________________________
 973 void 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  }