PWG2/SPECTRA/AliAnalysisTaskESDStrangeMC.cxx

   1
   2 /**************************************************************************
   3  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   4  *                                                                        *
   5  * Author: The ALICE Off-line Project.                                    *
   6  * Contributors are mentioned in the code where appropriate.              *
   7  *                                                                        *
   8  * Permission to use, copy, modify and distribute this software and its   *
   9  * documentation strictly for non-commercial purposes is hereby granted   *
  10  * without fee, provided that the above copyright notice appears in all   *
  11  * copies and that both the copyright notice and this permission notice   *
  12  * appear in the supporting documentation. The authors make no claims     *
  13  * about the suitability of this software for any purpose. It is          *
  14  * provided "as is" without express or implied warranty.                  *
  15  **************************************************************************/
  16
  17 // macro to study V0s, with Monte Carlo information access
  18 // loops over ESD files, and creates AliAODv0
  19 // Author: H.Ricaud, Helene.Ricaud@IReS.in2p3.fr
  20
  21
  22 #include "Riostream.h"
  23 #include "TChain.h"
  24 #include "TTree.h"
  25 #include "TH1F.h"
  26 #include "TH2F.h"
  27 #include "TList.h"
  28 #include "TMath.h"
  29 #include "TCanvas.h"
  30
  31 #include "AliAnalysisTask.h"
  32 #include "AliAnalysisManager.h"
  33
  34 #include "AliESDVertex.h"
  35 #include "AliESDEvent.h"
  36 #include "AliESDInputHandler.h"
  37 #include "AliESDtrack.h"
  38 #include "AliESDv0.h"
  39 #include "AliAODv0.h"
  40 #include "AliAODTrack.h"
  41
  42 #include "AliMCEventHandler.h"
  43 #include "AliMCEvent.h"
  44 #include "AliStack.h"
  45
  46 #include "AliLog.h"
  47
  48 #include "AliAnalysisTaskESDStrangeMC.h"
  49
  50
  51 ClassImp(AliAnalysisTaskESDStrangeMC)
  52
  53 //________________________________________________________________________
  54 AliAnalysisTaskESDStrangeMC::AliAnalysisTaskESDStrangeMC(const char *name)
  55   : AliAnalysisTask(name, ""), fESD(0), fListHist(),
  56     fHistPtMC(0),
  57     fHistMCMultiplicity(0),
  58     fHistMCPtVsYK0s(0),
  59     fHistMCPtVsYLambda(0),
  60     fHistMCPtVsYAntiLambda(0),
  61     fHistTrackPerEvent(0),
  62     fHistMCDaughterTrack(0),
  63     fHistPrimaryVertexX(0),
  64     fHistPrimaryVertexY(0),
  65     fHistPrimaryVertexZ(0),
  66     fHistDcaPosToPrimVertex(0),
  67     fHistDcaNegToPrimVertex(0),
  68     fHistDcaPosToPrimVertexZoom(0),
  69     fHistDcaNegToPrimVertexZoom(0),
  70     fHistRadiusV0(0),
  71     fHistDecayLengthV0(0),
  72     fHistDcaV0Daughters(0),
  73     fHistChi2(0),
  74     fHistCosPointAngle(0),
  75     fHistCosPointAngleZoom(0),
  76     fHistPtVsYK0s(0),
  77     fHistPtVsYK0sMI(0),
  78     fHistPtVsYLambda(0),
  79     fHistPtVsYLambdaMI(0),
  80     fHistPtVsYAntiLambda(0),
  81     fHistPtVsYAntiLambdaMI(0),
  82     fHistMassK0(0),
  83     fHistMassK0MI(0),
  84     fHistMassLambda(0),
  85     fHistMassLambdaMI(0),
  86     fHistMassAntiLambda(0),
  87     fHistMassAntiLambdaMI(0),
  88     fHistMassVsRadiusK0(0),
  89     fHistMassVsRadiusK0MI(0),
  90     fHistMassVsRadiusLambda(0),
  91     fHistMassVsRadiusLambdaMI(0),
  92     fHistMassVsRadiusAntiLambda(0),
  93     fHistMassVsRadiusAntiLambdaMI(0),
  94     fHistArmenterosPodolanski(0),
  95     fHistArmenterosPodolanskiMI(0),
  96     fHistAsMcPtK0(0),
  97     fHistAsMcPtK0MI(0),
  98     fHistAsMcPtLambda(0),
  99     fHistAsMcPtLambdaMI(0),
 100     fHistAsMcPtAntiLambda(0),
 101     fHistAsMcPtAntiLambdaMI(0),
 102     fHistAsMcPtZoomK0(0),
 103     fHistAsMcPtZoomK0MI(0),
 104     fHistAsMcPtZoomLambda(0),
 105     fHistAsMcPtZoomLambdaMI(0),
 106     fHistPidMcMassK0(0),
 107     fHistPidMcMassK0MI(0),
 108     fHistPidMcMassLambda(0),
 109     fHistPidMcMassLambdaMI(0),
 110     fHistPidMcMassAntiLambda(0),
 111     fHistPidMcMassAntiLambdaMI(0),
 112     fHistAsMcMassK0(0),
 113     fHistAsMcMassK0MI(0),
 114     fHistAsMcMassLambda(0),
 115     fHistAsMcMassLambdaMI(0),
 116     fHistAsMcMassAntiLambda(0),
 117     fHistAsMcMassAntiLambdaMI(0),
 118     fHistAsMcMassVsRadiusK0(0),
 119     fHistAsMcMassVsRadiusK0MI(0),
 120     fHistAsMcMassVsRadiusLambda(0),
 121     fHistAsMcMassVsRadiusLambdaMI(0),
 122     fHistAsMcMassVsRadiusAntiLambda(0),
 123     fHistAsMcMassVsRadiusAntiLambdaMI(0),
 124     fHistAsMcResxK0(0),
 125     fHistAsMcResyK0(0),
 126     fHistAsMcReszK0(0),
 127     fHistAsMcResrVsRadiusK0(0),
 128     fHistAsMcReszVsRadiusK0(0),
 129     fHistAsMcResxK0MI(0),
 130     fHistAsMcResyK0MI(0),
 131     fHistAsMcReszK0MI(0),
 132     fHistAsMcResrVsRadiusK0MI(0),
 133     fHistAsMcReszVsRadiusK0MI(0),
 134     fHistAsMcResxLambda(0),
 135     fHistAsMcResyLambda(0),
 136     fHistAsMcReszLambda(0),
 137     fHistAsMcResrVsRadiusLambda(0),
 138     fHistAsMcReszVsRadiusLambda(0),
 139     fHistAsMcResxLambdaMI(0),
 140     fHistAsMcResyLambdaMI(0),
 141     fHistAsMcReszLambdaMI(0),
 142     fHistAsMcResrVsRadiusLambdaMI(0),
 143     fHistAsMcReszVsRadiusLambdaMI(0),
 144     fHistAsMcResxAntiLambda(0),
 145     fHistAsMcResyAntiLambda(0),
 146     fHistAsMcReszAntiLambda(0),
 147     fHistAsMcResrVsRadiusAntiLambda(0),
 148     fHistAsMcReszVsRadiusAntiLambda(0),
 149     fHistAsMcResxAntiLambdaMI(0),
 150     fHistAsMcResyAntiLambdaMI(0),
 151     fHistAsMcReszAntiLambdaMI(0),
 152     fHistAsMcResrVsRadiusAntiLambdaMI(0),
 153     fHistAsMcReszVsRadiusAntiLambdaMI(0)
 154
 155 {
 156   // Constructor
 157
 158   // Define input and output slots here
 159   // Input slot #0 works with a TChain
 160   DefineInput(0, TChain::Class());
 161   // Output slot #0 writes into a TH1 container
 162   //DefineOutput(0, TH1F::Class());
 163   // Output slot #0 writes into a TList container
 164   DefineOutput(0, TList::Class());
 165 }
 166
 167 //________________________________________________________________________
 168 void AliAnalysisTaskESDStrangeMC::ConnectInputData(Option_t *)
 169 {
 170   // Connect ESD or AOD here
 171   // Called once
 172
 173   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
 174   if (!tree) {
 175     Printf("ERROR: Could not read chain from input slot 0");
 176   } else {
 177     // Disable all branches, we want to process only MC
 178     tree->SetBranchStatus("*", kFALSE);
 179     tree->SetBranchStatus("fTracks.*", kTRUE);
 180     tree->SetBranchStatus("fV0s.*", kTRUE);
 181
 182     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 183
 184     if (!esdH) {
 185       Printf("ERROR: Could not get ESDInputHandler");
 186     } else
 187       fESD = esdH->GetEvent();
 188   }
 189 }
 190
 191 //________________________________________________________________________
 192 void AliAnalysisTaskESDStrangeMC::CreateOutputObjects()
 193 {
 194   // Create histograms
 195   // Called once
 196
 197   fListHist = new TList();
 198
 199
 200   //***************
 201   // MC histograms
 202   //***************
 203   fHistPtMC = new TH1F("h1PtMC", "P_{T} distribution", 15, 0.1, 3.1);
 204   fHistPtMC->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 205   fHistPtMC->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 206   fHistPtMC->SetMarkerStyle(kFullCircle);
 207   fListHist->Add(fHistPtMC);
 208
 209   // Multiplicity
 210   fHistMCMultiplicity           = new TH1F("h1MCMultiplicity", "MC Multiplicity;Ntracks;Count", 201, -0.5, 200.5);
 211   fListHist->Add(fHistMCMultiplicity);
 212
 213   // Pt and rapidity distribution:
 214   fHistMCPtVsYK0s               = new TH2F("h2MCPtVsYK0s", "K^{0} candidates;p_{t} (GeV/c);rapidity",150,0,15,20,-10,10);
 215   fListHist->Add(fHistMCPtVsYK0s);
 216
 217   fHistMCPtVsYLambda            = new TH2F("h2MCPtVsYLambda", "#Lambda^{0} candidates;p_{t} (GeV/c);rapidity",150,0,15,20,-10,10);
 218   fListHist->Add(fHistMCPtVsYLambda);
 219
 220   fHistMCPtVsYAntiLambda        = new TH2F("h2MCPtVsYAntiLambda", "#bar{#Lambda}^{0} candidates;p_{t} (GeV/c);rapidity",150,0,15,20,-10,10);
 221   fListHist->Add(fHistMCPtVsYAntiLambda);
 222
 223
 224   //***********************************
 225   // Reconstructed particles histograms
 226   //***********************************
 227
 228   // multiplicity
 229   fHistTrackPerEvent           = new TH1F("h1TrackPerEvent", "Tracks per event;Number of Tracks;Number of Events",50,0,50);
 230   fListHist->Add(fHistTrackPerEvent);
 231
 232   fHistMCDaughterTrack         = new TH1F("h1MCDaughterTrack","Distribution of mc id for daughters;id tags;Counts",15,0,15);
 233   fListHist->Add(fHistMCDaughterTrack);
 234
 235   // Primary Vertex:
 236   fHistPrimaryVertexX          = new TH1F("h1PrimaryVertexX", "Primary Vertex Position X;Primary Vertex Position X (cm);Events",40,-1,1);
 237   fListHist->Add(fHistPrimaryVertexX);
 238
 239   fHistPrimaryVertexY          = new TH1F("h1PrimaryVertexY", "Primary Vertex Position Y;Primary Vertex Position Y (cm);Events",40,-1,1);
 240   fListHist->Add(fHistPrimaryVertexY);
 241
 242   fHistPrimaryVertexZ          = new TH1F("h1PrimaryVertexZ", "Primary Vertex Position Z;Primary Vertex Position Z (cm);Events",60,-5,5);
 243   fListHist->Add(fHistPrimaryVertexZ);
 244
 245   // Cut checks:
 246   fHistDcaPosToPrimVertex      = new TH2F("h2DcaPosToPrimVertex", "Positive V0 daughter;dca(cm);Status",500,0,5,2,-0.5,1.5);
 247   fListHist->Add(fHistDcaPosToPrimVertex);
 248
 249   fHistDcaNegToPrimVertex      = new TH2F("h2DcaNegToPrimVertex", "Negative V0 daughter;dca(cm);Status",500,0,5,2,-0.5,1.5);
 250   fListHist->Add(fHistDcaNegToPrimVertex);
 251
 252   fHistDcaPosToPrimVertexZoom  = new TH2F("h2DcaPosToPrimVertexZoom", "Positive V0 daughter;dca(cm);Status",100,0,0.1,2,-0.5,1.5);
 253   fListHist->Add(fHistDcaPosToPrimVertexZoom);
 254
 255   fHistDcaNegToPrimVertexZoom  = new TH2F("h2DcaNegToPrimVertexZoom", "Negative V0 daughter;dca(cm);Status",100,0,0.1,2,-0.5,1.5);
 256   fListHist->Add(fHistDcaNegToPrimVertexZoom);
 257
 258   fHistRadiusV0                = new TH2F("h2RadiusV0", "Radius;Radius(cm);Status",1000,0,100,2,-0.5,1.5);
 259   fListHist->Add(fHistRadiusV0);
 260
 261   fHistDecayLengthV0           = new TH2F("h2DecayLengthV0", "V0s decay Length;decay length(cm);Status", 200, 0, 100,2,-0.5,1.5);
 262   fListHist->Add(fHistDecayLengthV0);
 263
 264   fHistDcaV0Daughters          = new TH2F("h2DcaV0Daughters", "DCA between daughters;dca(cm);Status", 160, 0, 4,2,-0.5,1.5);
 265   fListHist->Add(fHistDcaV0Daughters);
 266
 267   fHistChi2                    = new TH2F("h2Chi2", "V0s chi2;chi2;Status", 33, 0, 33,2,-0.5,1.5);
 268   fListHist->Add(fHistChi2);
 269
 270   fHistCosPointAngle           = new TH2F("h2CosPointAngle", "Cosine of V0's pointing angle", 100,0,1,2,-0.5,1.5);
 271   fListHist->Add(fHistCosPointAngle);
 272
 273   fHistCosPointAngleZoom       = new TH2F("h2CosPointAngleZoom", "Cosine of V0's pointing angle", 100,0.9,1,2,-0.5,1.5);
 274   fListHist->Add(fHistCosPointAngleZoom);
 275
 276   // Pt and rapidity distribution:
 277   fHistPtVsYK0s                = new TH2F("h2PtVsYK0s", "K^{0} candidates;p_{t} (GeV/c);rapidity",150,0,15,30,-1.5,1.5);
 278   fListHist->Add(fHistPtVsYK0s);
 279   fHistPtVsYK0sMI              = new TH2F("h2PtVsYK0sMI", "K^{0} MI candidates;p_{t} (GeV/c);rapidity",150,0,15,30,-1.5,1.5);
 280   fListHist->Add(fHistPtVsYK0sMI);
 281
 282   fHistPtVsYLambda             = new TH2F("h2PtVsYLambda", "#Lambda^{0} candidates;p_{t} (GeV/c);rapidity",150,0,15,30,-1.5,1.5);
 283   fListHist->Add(fHistPtVsYLambda);
 284   fHistPtVsYLambdaMI           = new TH2F("h2PtVsYLambdaMI", "#Lambda^{0} MI candidates;p_{t} (GeV/c);rapidity",150,0,15,30,-1.5,1.5);
 285   fListHist->Add(fHistPtVsYLambdaMI);
 286
 287   fHistPtVsYAntiLambda         = new TH2F("h2PtVsYAntiLambda", "#bar{#Lambda}^{0} candidates;p_{t} (GeV/c);rapidity",150,0,15,30,-1.5,1.5);
 288   fListHist->Add(fHistPtVsYAntiLambda);
 289   fHistPtVsYAntiLambdaMI       = new TH2F("h2PtVsYAntiLambdaMI", "#bar{#Lambda}^{0} MI candidates;p_{t} (GeV/c);rapidity",150,0,15,30,-1.5,1.5);
 290   fListHist->Add(fHistPtVsYAntiLambdaMI);
 291
 292   // Mass:
 293   fHistMassK0                   = new TH1F("h1MassK0", "K^{0} candidates;M(#pi^{+}#pi^{-}) (GeV/c^{2});Counts", 100, 0.4, 0.6);
 294   fListHist->Add(fHistMassK0);
 295   fHistMassK0MI                 = new TH1F("h1MassK0MI", "K^{0} candidates;M(#pi^{+}#pi^{-}) (GeV/c^{2});Counts", 100, 0.4, 0.6);
 296   fListHist->Add(fHistMassK0MI);
 297
 298   fHistMassLambda               = new TH1F("h1MassLambda", "#Lambda^{0} candidates;M(p#pi^{-}) (GeV/c^{2});Counts", 75, 1.05, 1.2);
 299   fListHist->Add(fHistMassLambda);
 300   fHistMassLambdaMI             = new TH1F("h1MassLambdaMI", "#Lambda^{0} candidates;M(p#pi^{-}) (GeV/c^{2});Counts", 75, 1.05, 1.2);
 301   fListHist->Add(fHistMassLambdaMI);
 302
 303   fHistMassAntiLambda           = new TH1F("h1MassAntiLambda", "#bar{#Lambda}^{0} candidates;M(#bar{p}#pi^{+}) (GeV/c^{2});Counts", 75, 1.05, 1.2);
 304   fListHist->Add(fHistMassAntiLambda);
 305   fHistMassAntiLambdaMI         = new TH1F("h1MassAntiLambdaMI", "#bar{#Lambda}^{0} candidates;M(#bar{p}#pi^{+}) (GeV/c^{2});Counts", 75, 1.05, 1.2);
 306   fListHist->Add(fHistMassAntiLambdaMI);
 307
 308   // invariant mass vs radius
 309   const Double_t radius[10] = {0.0,2.5,2.9,3.9,7.6,15.0,23.9,37.8,42.8,100.0};
 310   Int_t NbinRadius        = 9;
 311   Int_t NbinInvMassLambda = 300;
 312
 313   fHistMassVsRadiusK0           = new TH2F("h2MassVsRadiusK0", "K^{0} reconstructed;radius (cm);M(#pi^{+}#pi^{-}) (GeV/c^{2})",NbinRadius,radius, 200, 0.4, 0.6);
 314   fListHist->Add(fHistMassVsRadiusK0);
 315
 316   fHistMassVsRadiusK0MI         = new TH2F("h2MassVsRadiusK0MI", "K^{0} MI reconstructed;radius (cm);M(#pi^{+}#pi^{-}) (GeV/c^{2})",NbinRadius,radius, 200, 0.4, 0.6);
 317   fListHist->Add(fHistMassVsRadiusK0MI);
 318
 319   fHistMassVsRadiusLambda       = new TH2F("h2MassVsRadiusLambda", "#Lambda reconstructed;radius (cm);M(p#pi^{-}) (GeV/c^{2})",NbinRadius,radius, 140, 1.06, 1.2);
 320   fListHist->Add(fHistMassVsRadiusLambda);
 321
 322   fHistMassVsRadiusLambdaMI     = new TH2F("h2MassVsRadiusLambdaMI", "#Lambda MI reconstructed;radius (cm);M(p#pi^{-}) (GeV/c^{2})",NbinRadius,radius, 140, 1.06, 1.2);
 323   fListHist->Add(fHistMassVsRadiusLambdaMI);
 324
 325   fHistMassVsRadiusAntiLambda   = new TH2F("h2MassVsRadiusAntiLambda", "#bar{#Lambda} reconstructed;radius (cm);M(#bar{p}#pi^{+}) (GeV/c^{2})",NbinRadius,radius, 140, 1.06, 1.2);
 326   fListHist->Add(fHistMassVsRadiusAntiLambda);
 327
 328   fHistMassVsRadiusAntiLambdaMI = new TH2F("h2MassVsRadiusAntiLambdaMI", "#bar{#Lambda} reconstructed;radius (cm);M(#bar{p}#pi^{+}) (GeV/c^{2})",NbinRadius,radius, 140, 1.06, 1.2);
 329   fListHist->Add(fHistMassVsRadiusAntiLambdaMI);
 330
 331
 332   fHistArmenterosPodolanski     = new TH2F("h2ArmenterosPodolanski","Armenteros-Podolanski phase space;#alpha;p_{t} arm",100,-1.0,1.0,50,0,0.5);
 333   fHistArmenterosPodolanskiMI   = new TH2F("h2ArmenterosPodolanskiMI","Armenteros-Podolanski phase space;#alpha;p_{t} arm",100,-1.0,1.0,50,0,0.5);
 334
 335
 336   //********************************
 337   // Associated particles histograms
 338   //********************************
 339   //Pt distribution
 340   fHistAsMcPtK0                = new TH1F("h1AsMcPtK0", "K^{0} associated;p_{t} (GeV/c);Counts", 150, 0, 15);
 341   fListHist->Add(fHistAsMcPtK0);
 342   fHistAsMcPtK0MI              = new TH1F("h1AsMcPtK0MI", "K^{0} associated;p_{t} (GeV/c);Counts", 150, 0, 15);
 343   fListHist->Add(fHistAsMcPtK0MI);
 344
 345   fHistAsMcPtLambda            = new TH1F("h1AsMcPtLambda", "#Lambda^{0} associated;p_{t} (GeV/c);Counts", 150, 0, 15);
 346   fListHist->Add(fHistAsMcPtLambda);
 347   fHistAsMcPtLambdaMI          = new TH1F("h1AsMcPtLambdaMI", "#Lambda^{0} associated;p_{t} (GeV/c);Counts", 150, 0, 15);
 348   fListHist->Add(fHistAsMcPtLambdaMI);
 349
 350   fHistAsMcPtAntiLambda        = new TH1F("h1AsMcPtAntiLambda", "#bar{#Lambda}^{0} associated;p_{t} (GeV/c);Counts", 150, 0, 15);
 351   fListHist->Add(fHistAsMcPtAntiLambda);
 352   fHistAsMcPtAntiLambdaMI      = new TH1F("h1AsMcPtAntiLambdaMI", "#bar{#Lambda}^{0} associated;p_{t} (GeV/c);Counts", 150, 0, 15);
 353   fListHist->Add(fHistAsMcPtAntiLambdaMI);
 354
 355   fHistAsMcPtZoomK0            = new TH1F("h1AsMcPtZoomK0", "K^{0} candidates in -1 <y<1;p_{t} (GeV/c);Counts",20,0,1);
 356   fListHist->Add(fHistAsMcPtZoomK0);
 357   fHistAsMcPtZoomK0MI         = new TH1F("h1AsMcPtZoomK0MI", "K^{0} MI candidates in -1 <y<1;p_{t} (GeV/c);Counts",20,0,1);
 358   fListHist->Add(fHistAsMcPtZoomK0MI);
 359
 360   fHistAsMcPtZoomLambda        = new TH1F("h1AsMcPtZoomLambda", "#Lambda^{0} candidates in -1 <y<1;p_{t} (GeV/c);Counts",20,0,1);
 361   fListHist->Add(fHistAsMcPtZoomLambda);
 362   fHistAsMcPtZoomLambdaMI      = new TH1F("h1AsMcPtZoomLambdaMI", "#Lambda^{0} MI candidates in -1 <y<1;p_{t} (GeV/c);Counts",20,0,1);
 363   fListHist->Add(fHistAsMcPtZoomLambdaMI);
 364
 365
 366   // Mass
 367   fHistPidMcMassK0             = new TH1F("h1PidMcMassK0", "K^{0} MC PId checked;M(#pi^{+}#pi^{-}) (GeV/c^{2});Counts", 100, 0.4, 0.6);
 368   fListHist->Add(fHistPidMcMassK0);
 369   fHistPidMcMassK0MI           = new TH1F("h1PidMcMassK0MI", "K^{0} MC PId checked;M(#pi^{+}#pi^{-}) (GeV/c^{2});Counts", 100, 0.4, 0.6);
 370   fListHist->Add(fHistPidMcMassK0MI);
 371
 372   fHistPidMcMassLambda         = new TH1F("h1PidMcMassLambda", "#Lambda^{0} MC PId checked;M(p#pi^{-}) (GeV/c^{2});Counts", 75, 1.05, 1.2);
 373   fListHist->Add(fHistPidMcMassLambda);
 374   fHistPidMcMassLambdaMI       = new TH1F("h1PidMcMassLambdaMI", "#Lambda^{0} MC PId checked;M(p#pi^{-}) (GeV/c^{2});Counts", 75, 1.05, 1.2);
 375   fListHist->Add(fHistPidMcMassLambdaMI);
 376
 377   fHistPidMcMassAntiLambda     = new TH1F("h1PidMcMassAntiLambda", "#bar{#Lambda}^{0} MC PId checked;M(#bar{p}#pi^{+}) (GeV/c^{2});Counts", 75, 1.05, 1.2);
 378   fListHist->Add(fHistPidMcMassAntiLambda);
 379   fHistPidMcMassAntiLambdaMI   = new TH1F("h1PidMcMassAntiLambdaMI", "#bar{#Lambda}^{0} MC PId checked;M(#bar{p}#pi^{+}) (GeV/c^{2});Counts", 75, 1.05, 1.2);
 380   fListHist->Add(fHistPidMcMassAntiLambdaMI);
 381
 382   fHistAsMcMassK0              = new TH1F("h1AsMcMassK0", "K^{0} associated;M(#pi^{+}#pi^{-}) (GeV/c^{2});Counts", 100, 0.4, 0.6);
 383   fListHist->Add(fHistAsMcMassK0);
 384   fHistAsMcMassK0MI            = new TH1F("h1AsMcMassK0MI", "K^{0} associated;M(#pi^{+}#pi^{-}) (GeV/c^{2});Counts", 100, 0.4, 0.6);
 385   fListHist->Add(fHistAsMcMassK0MI);
 386
 387   fHistAsMcMassLambda          = new TH1F("h1AsMcMassLambda", "#Lambda^{0} associated;M(p#pi^{-}) (GeV/c^{2});Counts", 75, 1.05, 1.2);
 388   fListHist->Add(fHistAsMcMassLambda);
 389   fHistAsMcMassLambdaMI        = new TH1F("h1AsMcMassLambdaMI", "#Lambda^{0} associated;M(p#pi^{-}) (GeV/c^{2});Counts", 75, 1.05, 1.2);
 390   fListHist->Add(fHistAsMcMassLambdaMI);
 391
 392   fHistAsMcMassAntiLambda      = new TH1F("h1AsMcMassAntiLambda", "#bar{#Lambda}^{0} associated;M(#bar{p}#pi^{+}) (GeV/c^{2});Counts", 75, 1.05, 1.2);
 393   fListHist->Add(fHistAsMcMassAntiLambda);
 394   fHistAsMcMassAntiLambdaMI    = new TH1F("h1AsMcMassAntiLambdaMI", "#bar{#Lambda}^{0} associated;M(#bar{p}#pi^{+}) (GeV/c^{2});Counts", 75, 1.05, 1.2);
 395   fListHist->Add(fHistAsMcMassAntiLambdaMI);
 396
 397
 398   // invariant mass vs radius
 399   fHistAsMcMassVsRadiusK0           = new TH2F("h2AsMcMassVsRadiusK0", "K^{0} associated;radius (cm);M(#pi^{+}#pi^{-}) (GeV/c^{2})",NbinRadius,radius, 500, 0.47, 0.52);
 400   fListHist->Add(fHistAsMcMassVsRadiusK0);
 401
 402   fHistAsMcMassVsRadiusK0MI         = new TH2F("h2AsMcMassVsRadiusK0MI", "K^{0} MI associated;radius (cm);M(#pi^{+}#pi^{-}) (GeV/c^{2})",NbinRadius,radius, 500, 0.47, 0.52);
 403   fListHist->Add(fHistAsMcMassVsRadiusK0MI);
 404
 405   fHistAsMcMassVsRadiusLambda       = new TH2F("h2AsMcMassVsRadiusLambda", "#Lambda associated;radius (cm);M(p#pi^{-}) (GeV/c^{2})",NbinRadius,radius, NbinInvMassLambda, 1.10, 1.13);
 406   fListHist->Add(fHistAsMcMassVsRadiusLambda);
 407
 408   fHistAsMcMassVsRadiusLambdaMI     = new TH2F("h2AsMcMassVsRadiusLambdaMI", "#Lambda MI associated;radius (cm);M(p#pi^{-}) (GeV/c^{2})",NbinRadius,radius, NbinInvMassLambda, 1.10, 1.13);
 409   fListHist->Add(fHistAsMcMassVsRadiusLambdaMI);
 410
 411   fHistAsMcMassVsRadiusAntiLambda   = new TH2F("h2AsMcMassVsRadiusAntiLambda", "#bar{#Lambda} associated;radius (cm);M(#bar{p}#pi^{+}) (GeV/c^{2})",NbinRadius,radius,NbinInvMassLambda , 1.10, 1.13);
 412   fListHist->Add(fHistAsMcMassVsRadiusAntiLambda);
 413
 414   fHistAsMcMassVsRadiusAntiLambdaMI = new TH2F("h2AsMcMassVsRadiusAntiLambdaMI", "#bar{#Lambda} MI associated;radius (cm);M(#bar{p}#pi^{+}) (GeV/c^{2})",NbinRadius,radius,NbinInvMassLambda , 1.10, 1.13);
 415   fListHist->Add(fHistAsMcMassVsRadiusAntiLambdaMI);
 416
 417
 418   // Resolution
 419   fHistAsMcResxK0                   = new TH1F("h1AsMcResxK0", "K^{0} associated;#Delta x (cm);Counts", 50, -0.25, 0.25);
 420   fListHist->Add(fHistAsMcResxK0);
 421   fHistAsMcResyK0                   = new TH1F("h1AsMcResyK0", "K^{0} associated;#Delta y (cm);Counts", 50, -0.25, 0.25);
 422   fListHist->Add(fHistAsMcResyK0);
 423   fHistAsMcReszK0                   = new TH1F("h1AsMcReszK0", "K^{0} associated;#Delta z (cm);Counts", 50, -0.25, 0.25);
 424   fListHist->Add(fHistAsMcReszK0);
 425   fHistAsMcResrVsRadiusK0           = new TH2F("h2AsMcResrVsRadiusK0", "K^{0} associated;Radius (cm);#Delta r (cm)",8,radius, 50, -0.25, 0.25);
 426   fListHist->Add(fHistAsMcResrVsRadiusK0);
 427   fHistAsMcReszVsRadiusK0           = new TH2F("h2AsMcReszVsRadiusK0", "K^{0} associated;Radius (cm);#Delta z (cm)",8,radius, 50, -0.25, 0.25);
 428   fListHist->Add(fHistAsMcReszVsRadiusK0);
 429
 430   fHistAsMcResxK0MI                 = new TH1F("h1AsMcResxK0MI", "K^{0} MI associated;#Delta x (cm);Counts", 50, -0.25, 0.25);
 431   fListHist->Add(fHistAsMcResxK0MI);
 432   fHistAsMcResyK0MI                 = new TH1F("h1AsMcResyK0MI", "K^{0} MI associated;#Delta y (cm);Counts", 50, -0.25, 0.25);
 433   fListHist->Add(fHistAsMcResyK0MI);
 434   fHistAsMcReszK0MI                 = new TH1F("h1AsMcReszK0MI", "K^{0} MI associated;#Delta z (cm);Counts", 50, -0.25, 0.25);
 435   fListHist->Add(fHistAsMcReszK0MI);
 436   fHistAsMcResrVsRadiusK0MI         = new TH2F("h2AsMcResrVsRadiusK0MI", "K^{0} MI associated;Radius (cm);#Delta r (cm)",8,radius, 50, -0.25, 0.25);
 437   fListHist->Add(fHistAsMcResrVsRadiusK0MI);
 438   fHistAsMcReszVsRadiusK0MI         = new TH2F("h2AsMcReszVsRadiusK0MI", "K^{0} MI associated;Radius (cm);#Delta z (cm)",8,radius, 50, -0.25, 0.25);
 439   fListHist->Add(fHistAsMcReszVsRadiusK0MI);
 440
 441   fHistAsMcResxLambda             = new TH1F("h1AsMcResxLambda", "#Lambda^{0} associated;#Delta x (cm);Counts", 50, -0.25, 0.25);
 442   fListHist->Add(fHistAsMcResxLambda);
 443   fHistAsMcResyLambda             = new TH1F("h1AsMcResyLambda", "#Lambda^{0} associated;#Delta y (cm);Counts", 50, -0.25, 0.25);
 444   fListHist->Add(fHistAsMcResyLambda);
 445   fHistAsMcReszLambda             = new TH1F("h1AsMcReszLambda", "#Lambda^{0} associated;#Delta z (cm);Counts", 50, -0.25, 0.25);
 446   fListHist->Add(fHistAsMcReszLambda);
 447   fHistAsMcResrVsRadiusLambda     = new TH2F("h2AsMcResrVsRadiusLambda", "#Lambda^{0} associated;Radius (cm);#Delta r (cm)",8,radius, 50, -0.25, 0.25);
 448   fListHist->Add(fHistAsMcResrVsRadiusLambda);
 449   fHistAsMcReszVsRadiusLambda     = new TH2F("h2AsMcReszVsRadiusLambda", "#Lambda^{0} associated;Radius (cm);#Delta z (cm)",8,radius, 50, -0.25, 0.25);
 450   fListHist->Add(fHistAsMcReszVsRadiusLambda);
 451
 452   fHistAsMcResxLambdaMI           = new TH1F("h1AsMcResxLambdaMI", "#Lambda^{0} MI associated;#Delta x (cm);Counts", 50, -0.25, 0.25);
 453   fListHist->Add(fHistAsMcResxLambdaMI);
 454   fHistAsMcResyLambdaMI           = new TH1F("h1AsMcResyLambdaMI", "#Lambda^{0} MI associated;#Delta y (cm);Counts", 50, -0.25, 0.25);
 455   fListHist->Add(fHistAsMcResyLambdaMI);
 456   fHistAsMcReszLambdaMI           = new TH1F("h1AsMcReszLambdaMI", "#Lambda^{0} MI associated;#Delta z (cm);Counts", 50, -0.25, 0.25);
 457   fListHist->Add(fHistAsMcReszLambdaMI);
 458   fHistAsMcResrVsRadiusLambdaMI   = new TH2F("h2AsMcResrVsRadiusLambdaMI", "#Lambda^{0} MI associated;Radius (cm);#Delta r (cm)",8,radius, 50, -0.25, 0.25);
 459   fListHist->Add(fHistAsMcResrVsRadiusLambdaMI);
 460   fHistAsMcReszVsRadiusLambdaMI   = new TH2F("h2AsMcReszVsRadiusLambdaMI", "#Lambda^{0} MI associated;Radius (cm);#Delta z (cm)",8,radius, 50, -0.25, 0.25);
 461   fListHist->Add(fHistAsMcReszVsRadiusLambdaMI);
 462
 463   fHistAsMcResxAntiLambda         = new TH1F("h1AsMcResxAntiLambda", "#bar{#Lambda}^{0} associated;#Delta x (cm);Counts", 50, -0.25, 0.25);
 464   fListHist->Add(fHistAsMcResxAntiLambda);
 465   fHistAsMcResyAntiLambda         = new TH1F("h1AsMcResyAntiLambda", "#bar{#Lambda}^{0} associated;#Delta y (cm);Counts", 50, -0.25, 0.25);
 466   fListHist->Add(fHistAsMcResyAntiLambda);
 467   fHistAsMcReszAntiLambda         = new TH1F("h1AsMcReszAntiLambda", "#bar{#Lambda}^{0} associated;#Delta z (cm);Counts", 50, -0.25, 0.25);
 468   fListHist->Add(fHistAsMcReszAntiLambda);
 469   fHistAsMcResrVsRadiusAntiLambda = new TH2F("h2AsMcResrVsRadiusAntiLambda", "#bar{#Lambda}^{0} associated;Radius (cm);#Delta r (cm)",8,radius, 50, -0.25, 0.25);
 470   fListHist->Add(fHistAsMcResrVsRadiusAntiLambda);
 471   fHistAsMcReszVsRadiusAntiLambda = new TH2F("h2AsMcReszVsRadiusAntiLambda", "#bar{#Lambda}^{0} associated;Radius (cm);#Delta z (cm)",8,radius, 50, -0.25, 0.25);
 472   fListHist->Add(fHistAsMcReszVsRadiusAntiLambda);
 473
 474   fHistAsMcResxAntiLambdaMI         = new TH1F("h1AsMcResxAntiLambdaMI", "#bar{#Lambda}^{0} MI associated;#Delta x (cm);Counts", 50, -0.25, 0.25);
 475   fListHist->Add(fHistAsMcResxAntiLambdaMI);
 476   fHistAsMcResyAntiLambdaMI         = new TH1F("h1AsMcResyAntiLambdaMI", "#bar{#Lambda}^{0} MI associated;#Delta y (cm);Counts", 50, -0.25, 0.25);
 477   fListHist->Add(fHistAsMcResyAntiLambdaMI);
 478   fHistAsMcReszAntiLambdaMI         = new TH1F("h1AsMcReszAntiLambdaMI", "#bar{#Lambda}^{0} MI associated;#Delta z (cm);Counts", 50, -0.25, 0.25);
 479   fListHist->Add(fHistAsMcReszAntiLambdaMI);
 480   fHistAsMcResrVsRadiusAntiLambdaMI = new TH2F("h2AsMcResrVsRadiusAntiLambdaMI", "#bar{#Lambda}^{0} MI associated;Radius (cm);#Delta r (cm)",8,radius, 50, -0.25, 0.25);
 481   fListHist->Add(fHistAsMcResrVsRadiusAntiLambdaMI);
 482   fHistAsMcReszVsRadiusAntiLambdaMI = new TH2F("h2AsMcReszVsRadiusAntiLambdaMI", "#bar{#Lambda}^{0} MI associated;Radius (cm);#Delta z (cm)",8,radius, 50, -0.25, 0.25);
 483   fListHist->Add(fHistAsMcReszVsRadiusAntiLambdaMI);
 484
 485
 486 }
 487
 488 //________________________________________________________________________
 489 void AliAnalysisTaskESDStrangeMC::Exec(Option_t *)
 490 {
 491   AliLog::SetGlobalLogLevel(AliLog::kError);
 492
 493   //**********************************************
 494   // Check Monte Carlo information access first:
 495   //**********************************************
 496   // Process MC truth, therefore we receive the AliAnalysisManager and ask it for the AliMCEventHandler
 497   // This handler can return the current MC event
 498
 499   AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 500   if (!eventHandler) {
 501     Printf("ERROR: Could not retrieve MC event handler");
 502     return;
 503   }
 504
 505   AliMCEvent* mcEvent = eventHandler->MCEvent();
 506   if (!mcEvent) {
 507      Printf("ERROR: Could not retrieve MC event");
 508      return;
 509   }
 510
 511   AliStack* stack = mcEvent->Stack();
 512
 513   //**********************************************
 514   // MC loop
 515   // Called for each event
 516   //**********************************************
 517
 518   //Printf("MC particles: %d", mcEvent->GetNumberOfTracks());
 519   //Printf("MC primary particles: %d", stack->GetNprimary());
 520
 521   for (Int_t iTracksMC = 0; iTracksMC < mcEvent->GetNumberOfTracks(); iTracksMC++) {
 522     AliMCParticle* trackMC = mcEvent->GetTrack(iTracksMC);
 523     if (!trackMC) {
 524       Printf("ERROR: Could not receive track %d (mc loop)", iTracksMC);
 525       continue;
 526     }
 527     fHistMCMultiplicity->Fill(mcEvent->GetNumberOfTracks());
 528     fHistPtMC->Fill(trackMC->Pt());
 529   } //end MC tracks loop
 530
 531   Double_t y =999.9;
 532   for (Int_t iMc = 0; iMc < stack->GetNprimary(); ++iMc) {
 533       TParticle *p0 = stack->Particle(iMc);
 534       if (!p0) {
 535         Printf("ERROR: particle with label %d not found in stack (mc loop)", iMc);
 536         continue;
 537       }
 538       y=myRap(p0->Energy(),p0->Pz());
 539       if (p0->GetPdgCode()==310) fHistMCPtVsYK0s->Fill(p0->Pt(),y);
 540       else if (p0->GetPdgCode()==3122)  fHistMCPtVsYLambda->Fill(p0->Pt(),y);
 541       else if (p0->GetPdgCode()==-3122) fHistMCPtVsYAntiLambda->Fill(p0->Pt(),y);
 542
 543     }// end loop MC primary particles
 544
 545
 546   //*********************************************
 547   // ESD loop
 548   // Called for each event
 549   //*********************************************
 550
 551   if (!fESD) {
 552     Printf("ERROR: fESD not available");
 553     return;
 554   }
 555
 556   //Printf("There are %d tracks in this event", fESD->GetNumberOfTracks());
 557   if (!(fESD->GetNumberOfTracks())) {
 558     //Printf("No ESD track in the event");
 559   return;
 560   }
 561   fHistTrackPerEvent->Fill(fESD->GetNumberOfTracks());
 562
 563   // Track loop to fill a pT spectrum
 564   for (Int_t iesdTracks = 0; iesdTracks < fESD->GetNumberOfTracks(); iesdTracks++) {
 565     AliESDtrack* esdtrack = fESD->GetTrack(iesdTracks);
 566     if (!esdtrack) {
 567       Printf("ERROR: Could not receive track %d", iesdTracks);
 568       continue;
 569     }
 570   }
 571
 572   // Primary Vertex
 573   Double_t  PrimaryVtxPosition[3];
 574   Double_t  PrimaryVtxCov[6];
 575
 576   const AliESDVertex *primaryVtx = fESD->GetPrimaryVertex();
 577
 578   primaryVtx->GetXYZ(PrimaryVtxPosition);
 579   primaryVtx->GetCovMatrix(PrimaryVtxCov);
 580
 581   AliAODVertex *primary = new AliAODVertex(PrimaryVtxPosition, PrimaryVtxCov, primaryVtx->GetChi2toNDF(), NULL, -1, AliAODVertex::kPrimary);
 582
 583   fHistPrimaryVertexX->Fill(primary->GetX());
 584   fHistPrimaryVertexY->Fill(primary->GetY());
 585   fHistPrimaryVertexZ->Fill(primary->GetZ());
 586
 587
 588   // V0 variables:
 589   // to get info from ESD files and fill AliAODVertex:
 590   Float_t   tdcaPosToPrimVertexXYZ[2], tdcaNegToPrimVertexXYZ[2]; // ..[0] = Impact parameter in XY plane and ..[1] = Impact parameter in Z
 591   Double_t  tdcaDaughterToPrimVertex[2];                          // ..[0] = Pos and ..[1] = Neg
 592   Double_t  tdcaV0Daughters     = 0, tdcaV0ToPrimVertex   = 0;
 593   Double_t  tMomPos[3];
 594   Double_t  tMomNeg[3];
 595   Double_t  V0Position[3];
 596   Double_t  V0Cov[6];
 597
 598   // to fill AliAODtrack:
 599   Double_t  TrackP[3];
 600   Double_t  TrackPosition[3];
 601   Double_t  TrackcovTr[21];
 602   Double_t  Trackpid[10];
 603
 604   Int_t    lIndexTrackPos       = 0, lIndexTrackNeg       = 0;
 605   UInt_t   lLabelTrackPos       = 0, lLabelTrackNeg       = 0;
 606   Int_t    lCheckPIdK0Short     = 0, lCheckMcK0Short      = 0;
 607   Int_t    lCheckPIdLambda      = 0, lCheckMcLambda       = 0;
 608   Int_t    lCheckPIdAntiLambda  = 0, lCheckMcAntiLambda   = 0;
 609
 610   Double_t mcPosX        = 0,  mcPosY  = 0, mcPosZ  = 0;
 611   Double_t rcPosX        = 0,  rcPosY  = 0, rcPosZ  = 0;
 612   Double_t rcPosR        = 0,  mcPosR  = 0;
 613   Double_t cosPointAngle = 0;
 614   Double_t deltaPos2     = 0;
 615   Double_t deltaPos[3];
 616
 617   Int_t    myStatus  = 0;
 618   Double_t pt        = 0;
 619
 620   AliAODTrack  *myPosAodTrack  = new AliAODTrack();
 621   AliAODTrack  *myNegAodTrack  = new AliAODTrack();
 622   AliAODVertex *myAODVertex    = new AliAODVertex();
 623   AliAODv0     *myAODv0        = new AliAODv0();
 624
 625   // V0 loop
 626   for (Int_t iV0 = 0; iV0 < fESD->GetNumberOfV0s(); iV0++) {
 627
 628     AliESDv0 *v0 = fESD->GetV0(iV0);
 629     if (!v0) continue;
 630     myAODv0->ResetV0();
 631
 632     // AliAODVertex
 633     v0->GetXYZ(V0Position[0], V0Position[1], V0Position[2]);
 634     v0->GetPosCov(V0Cov);
 635     myAODVertex->SetPosition(V0Position[0],V0Position[1],V0Position[2]);
 636     myAODVertex->SetCovMatrix(V0Cov);
 637     myAODVertex->SetChi2perNDF(v0->GetChi2V0());
 638     myAODVertex->SetID((Short_t)iV0);
 639     myAODVertex->SetParent(primary);
 640     myAODVertex->SetType(AliAODVertex::kV0);
 641
 642     // AliAODtrack (V0 Daughters)
 643     lIndexTrackPos = TMath::Abs(v0->GetPindex());
 644     lIndexTrackNeg = TMath::Abs(v0->GetNindex());
 645     AliESDtrack *TrackPos = fESD->GetTrack(lIndexTrackPos);
 646     AliESDtrack *TrackNeg = fESD->GetTrack(lIndexTrackNeg);
 647     if (!TrackPos || !TrackNeg) {
 648       Printf("ERROR: Could not retreive one of the daughter track");
 649       continue;
 650     }
 651     lLabelTrackPos = (UInt_t)TMath::Abs(TrackPos->GetLabel());
 652     lLabelTrackNeg = (UInt_t)TMath::Abs(TrackNeg->GetLabel());
 653
 654     myPosAodTrack->SetID((Short_t)(TrackPos->GetID()));
 655     myPosAodTrack->SetLabel(lLabelTrackPos);
 656     TrackPos->GetPxPyPz(TrackP);
 657     myPosAodTrack->SetP(TrackP);
 658     TrackPos->GetXYZ(TrackPosition);
 659     myPosAodTrack->SetPosition(TrackPosition,kFALSE);
 660     TrackPos->GetCovarianceXYZPxPyPz(TrackcovTr);
 661     myPosAodTrack->SetCovMatrix(TrackcovTr);
 662     TrackPos->GetESDpid(Trackpid);
 663     myPosAodTrack->SetPID(Trackpid);
 664     myPosAodTrack->SetCharge((Short_t)(TrackPos->Charge()));
 665     myPosAodTrack->SetITSClusterMap(TrackPos->GetITSClusterMap());
 666     myPosAodTrack->SetProdVertex(myAODVertex);
 667     myPosAodTrack->SetUsedForVtxFit(kTRUE);
 668     myPosAodTrack->SetUsedForPrimVtxFit(kFALSE);
 669     myPosAodTrack->SetType(AliAODTrack::kSecondary);
 670     myPosAodTrack->ConvertAliPIDtoAODPID();
 671
 672     myNegAodTrack->SetID((Short_t)(TrackNeg->GetID()));
 673     myNegAodTrack->SetLabel(lLabelTrackNeg);
 674     TrackNeg->GetPxPyPz(TrackP);
 675     myNegAodTrack->SetP(TrackP);
 676     TrackNeg->GetXYZ(TrackPosition);
 677     myNegAodTrack->SetPosition(TrackPosition,kFALSE);
 678     TrackNeg->GetCovarianceXYZPxPyPz(TrackcovTr);
 679     myNegAodTrack->SetCovMatrix(TrackcovTr);
 680     TrackNeg->GetESDpid(Trackpid);
 681     myNegAodTrack->SetPID(Trackpid);
 682     myNegAodTrack->SetCharge((Short_t)(TrackNeg->Charge()));
 683     myNegAodTrack->SetITSClusterMap(TrackPos->GetITSClusterMap());
 684     myNegAodTrack->SetProdVertex(myAODVertex);
 685     myNegAodTrack->SetUsedForVtxFit(kTRUE);
 686     myNegAodTrack->SetUsedForPrimVtxFit(kFALSE);
 687     myNegAodTrack->SetType(AliAODTrack::kSecondary);
 688     myNegAodTrack->ConvertAliPIDtoAODPID();
 689
 690     myAODVertex->AddDaughter(myPosAodTrack);
 691     myAODVertex->AddDaughter(myNegAodTrack);
 692
 693     // filling myAODv0
 694     tdcaV0Daughters    = v0->GetDcaV0Daughters();
 695     tdcaV0ToPrimVertex = v0->GetD(primary->GetX(),primary->GetY(),primary->GetZ());
 696
 697     if (TrackPos) TrackPos->GetImpactParameters(tdcaPosToPrimVertexXYZ[0],tdcaPosToPrimVertexXYZ[1]);
 698     if (TrackNeg) TrackNeg->GetImpactParameters(tdcaNegToPrimVertexXYZ[0],tdcaNegToPrimVertexXYZ[1]);
 699     tdcaDaughterToPrimVertex[0] = TMath::Sqrt(tdcaPosToPrimVertexXYZ[0]*tdcaPosToPrimVertexXYZ[0]+tdcaPosToPrimVertexXYZ[1]*tdcaPosToPrimVertexXYZ[1]);
 700     tdcaDaughterToPrimVertex[1] = TMath::Sqrt(tdcaNegToPrimVertexXYZ[0]*tdcaNegToPrimVertexXYZ[0]+tdcaNegToPrimVertexXYZ[1]*tdcaNegToPrimVertexXYZ[1]);
 701
 702     v0->GetPPxPyPz(tMomPos[0],tMomPos[1],tMomPos[2]);
 703     v0->GetNPxPyPz(tMomNeg[0],tMomNeg[1],tMomNeg[2]);
 704
 705     myAODv0->Fill(myAODVertex, tdcaV0Daughters, tdcaV0ToPrimVertex, tMomPos, tMomNeg, tdcaDaughterToPrimVertex);
 706
 707     // Checking if the v0s is associated with a Monte Carlo particle
 708     //AliMCParticle  *lPartPos= mcEvent->GetTrack(lLabelTrackPos);
 709     //AliMCParticle  *lPartNeg= mcEvent->GetTrack(lLabelTrackNeg);
 710     TParticle  *lPartPos = stack->Particle(lLabelTrackPos);
 711     TParticle  *lPartNeg = stack->Particle(lLabelTrackNeg);
 712     Int_t lPartPosMother = lPartPos->GetFirstMother();
 713     Int_t lPartNegMother = lPartNeg->GetFirstMother();
 714
 715     lCheckPIdK0Short    = 0; lCheckMcK0Short     = 0;
 716     lCheckPIdLambda     = 0; lCheckMcLambda      = 0;
 717     lCheckPIdAntiLambda = 0; lCheckMcAntiLambda  = 0;
 718
 719     if( (lPartPosMother==-1) || (lPartNegMother==-1) ) {
 720       fHistMCDaughterTrack->Fill(1);
 721     }
 722     else if( (lPartPos->GetPdgCode()==+211) &&
 723              (lPartNeg->GetPdgCode()==-211) ) {
 724       lCheckPIdK0Short    = 1;
 725       fHistMCDaughterTrack->Fill(3);
 726       if ( (lPartPosMother==lPartNegMother) &&
 727            (stack->Particle(lPartPosMother)->GetPdgCode()==310) ){
 728         lCheckMcK0Short  = 1;
 729         mcPosX = lPartPos->Vx();
 730         mcPosY = lPartPos->Vy();
 731         mcPosZ = lPartPos->Vz();
 732         mcPosR = TMath::Sqrt(mcPosX*mcPosX+mcPosY*mcPosY);
 733       }
 734     }
 735     else if( (lPartPos->GetPdgCode()==+2212) &&
 736              (lPartNeg->GetPdgCode()==-211) ) {
 737       lCheckPIdLambda     = 1;
 738       fHistMCDaughterTrack->Fill(5);
 739       if ( (lPartPosMother==lPartNegMother) &&
 740            (stack->Particle(lPartPosMother)->GetPdgCode()==3122) ){
 741         lCheckMcLambda  = 1;
 742         mcPosX = lPartPos->Vx();
 743         mcPosY = lPartPos->Vy();
 744         mcPosZ = lPartPos->Vz();
 745         mcPosR = TMath::Sqrt(mcPosX*mcPosX+mcPosY*mcPosY);
 746       }
 747     }
 748     else if( (lPartPos->GetPdgCode()==211) &&
 749              (lPartNeg->GetPdgCode()==-2212) ) {
 750       lCheckPIdAntiLambda = 1;
 751       fHistMCDaughterTrack->Fill(7);
 752       if ( (lPartPosMother==lPartNegMother) &&
 753            (stack->Particle(lPartPosMother)->GetPdgCode()==-3122) ){
 754         lCheckMcAntiLambda  = 1;
 755         mcPosX = lPartPos->Vx();
 756         mcPosY = lPartPos->Vy();
 757         mcPosZ = lPartPos->Vz();
 758         mcPosR = TMath::Sqrt(mcPosX*mcPosX+mcPosY*mcPosY);
 759       }
 760     }
 761
 762     // Reconstructed V0 position and Cos pointing angle:
 763     rcPosX   = myAODv0->DecayVertexV0X();
 764     rcPosY   = myAODv0->DecayVertexV0Y();
 765     rcPosZ   = myAODv0->DecayVertexV0Z();
 766     rcPosR   = TMath::Sqrt(rcPosX*rcPosX+rcPosY*rcPosY);
 767
 768     deltaPos[0]   = rcPosX - PrimaryVtxPosition[0];
 769     deltaPos[1]   = rcPosY - PrimaryVtxPosition[1];
 770     deltaPos[2]   = rcPosZ - PrimaryVtxPosition[2];
 771     deltaPos2     = deltaPos[0]*deltaPos[0] + deltaPos[1]*deltaPos[1] + deltaPos[2]*deltaPos[2];
 772     cosPointAngle = (deltaPos[0]*(myAODv0->MomV0X()) + deltaPos[1]*(myAODv0->MomV0Y()) + deltaPos[2]*(myAODv0->MomV0Z()))/TMath::Sqrt((myAODv0->Ptot2V0())*deltaPos2);
 773
 774
 775
 776     myStatus = v0->GetOnFlyStatus();
 777     pt     = TMath::Sqrt(myAODv0->Ptot2V0());
 778
 779     // filling histograms
 780     fHistDcaPosToPrimVertex->Fill(myAODv0->DcaPosToPrimVertex(),myStatus);
 781     fHistDcaNegToPrimVertex->Fill(myAODv0->DcaNegToPrimVertex(),myStatus);
 782     fHistDcaPosToPrimVertexZoom->Fill(myAODv0->DcaPosToPrimVertex(),myStatus);
 783     fHistDcaNegToPrimVertexZoom->Fill(myAODv0->DcaNegToPrimVertex(),myStatus);
 784     fHistRadiusV0->Fill(myAODv0->RadiusV0(),myStatus);
 785     fHistDecayLengthV0->Fill(myAODv0->DecayLengthV0(PrimaryVtxPosition),myStatus);
 786     fHistDcaV0Daughters->Fill(myAODv0->DcaV0Daughters(),myStatus);
 787     fHistChi2->Fill(myAODv0->Chi2V0(),myStatus);
 788     fHistCosPointAngle->Fill(cosPointAngle,myStatus);
 789     if (cosPointAngle >= 0.9) fHistCosPointAngleZoom->Fill(cosPointAngle,myStatus);
 790     if (!myStatus) {
 791       fHistPtVsYK0s->Fill(pt,myAODv0->RapK0Short());
 792       fHistPtVsYLambda->Fill(pt,myAODv0->RapLambda());
 793       fHistPtVsYAntiLambda->Fill(pt,myAODv0->RapLambda());
 794       fHistArmenterosPodolanski->Fill(myAODv0->AlphaV0(),myAODv0->PtArmV0());
 795     }
 796     else {
 797       fHistPtVsYK0sMI->Fill(pt,myAODv0->RapK0Short());
 798       fHistPtVsYLambdaMI->Fill(pt,myAODv0->RapLambda());
 799       fHistPtVsYAntiLambdaMI->Fill(pt,myAODv0->RapLambda());
 800       fHistArmenterosPodolanskiMI->Fill(myAODv0->AlphaV0(),myAODv0->PtArmV0());
 801     }
 802     // K0s associated histograms:
 803     if (TMath::Abs(myAODv0->RapK0Short()) < 1) {
 804       switch (myStatus){
 805       case 0 :
 806         fHistMassK0->Fill(myAODv0->MassK0Short());
 807         fHistMassVsRadiusK0->Fill(myAODv0->RadiusV0(),myAODv0->MassK0Short());
 808         if(lCheckPIdK0Short) fHistPidMcMassK0->Fill(myAODv0->MassK0Short());
 809         if(lCheckMcK0Short) {
 810           fHistAsMcMassK0->Fill(myAODv0->MassK0Short());
 811           fHistAsMcPtK0->Fill(pt);
 812           if (pt <= 1) fHistAsMcPtZoomK0->Fill(pt);
 813           fHistAsMcResxK0->Fill(rcPosX-mcPosX);
 814           fHistAsMcResyK0->Fill(rcPosY-mcPosY);
 815           fHistAsMcReszK0->Fill(rcPosZ-mcPosZ);
 816           fHistAsMcResrVsRadiusK0->Fill(rcPosR,rcPosR-mcPosR);
 817           fHistAsMcReszVsRadiusK0->Fill(rcPosR,rcPosZ-mcPosZ);
 818         }
 819         break;
 820
 821       case 1 :
 822         fHistMassK0MI->Fill(myAODv0->MassK0Short());
 823         fHistMassVsRadiusK0MI->Fill(myAODv0->RadiusV0(),myAODv0->MassK0Short());
 824         if(lCheckPIdK0Short) fHistPidMcMassK0MI->Fill(myAODv0->MassK0Short());
 825         if(lCheckMcK0Short) {
 826           fHistAsMcMassK0MI->Fill(myAODv0->MassK0Short());
 827           fHistAsMcPtK0MI->Fill(pt);
 828           if (pt <= 1) fHistAsMcPtZoomK0MI->Fill(pt);
 829           fHistAsMcResxK0MI->Fill(rcPosX-mcPosX);
 830           fHistAsMcResyK0MI->Fill(rcPosY-mcPosY);
 831           fHistAsMcReszK0MI->Fill(rcPosZ-mcPosZ);
 832           fHistAsMcResrVsRadiusK0MI->Fill(rcPosR,rcPosR-mcPosR);
 833           fHistAsMcReszVsRadiusK0MI->Fill(rcPosR,rcPosZ-mcPosZ);
 834         }
 835         break;
 836
 837       }
 838     }
 839     // Lambda and AntiLambda associated histograms:
 840     if (TMath::Abs(myAODv0->RapLambda()) < 1) {
 841       switch (myStatus){
 842       case 0 :
 843         fHistMassLambda->Fill(myAODv0->MassLambda());
 844         fHistMassAntiLambda->Fill(myAODv0->MassAntiLambda());
 845         fHistMassVsRadiusLambda->Fill(myAODv0->RadiusV0(),myAODv0->MassLambda());
 846         fHistMassVsRadiusAntiLambda->Fill(myAODv0->RadiusV0(),myAODv0->MassAntiLambda());
 847         if(lCheckPIdLambda) fHistPidMcMassLambda->Fill(myAODv0->MassLambda());
 848         else if (lCheckPIdAntiLambda) fHistPidMcMassAntiLambda->Fill(myAODv0->MassAntiLambda());
 849         if(lCheckMcLambda) {
 850           fHistAsMcMassLambda->Fill(myAODv0->MassLambda());
 851           fHistAsMcPtLambda->Fill(pt);
 852           if (pt <= 1) fHistAsMcPtZoomLambda->Fill(pt);
 853           fHistAsMcMassVsRadiusLambda->Fill(myAODv0->RadiusV0(),myAODv0->MassLambda());
 854           fHistAsMcResxLambda->Fill(rcPosX-mcPosX);
 855           fHistAsMcResyLambda->Fill(rcPosY-mcPosY);
 856           fHistAsMcReszLambda->Fill(rcPosZ-mcPosZ);
 857           fHistAsMcResrVsRadiusLambda->Fill(rcPosR,rcPosR-mcPosR);
 858           fHistAsMcReszVsRadiusLambda->Fill(rcPosR,rcPosZ-mcPosZ);
 859
 860         }
 861         else if(lCheckMcAntiLambda) {
 862           fHistAsMcMassAntiLambda->Fill(myAODv0->MassAntiLambda());
 863           fHistAsMcPtAntiLambda->Fill(pt);
 864           fHistAsMcMassVsRadiusAntiLambda->Fill(myAODv0->RadiusV0(),myAODv0->MassAntiLambda());
 865           fHistAsMcResxAntiLambda->Fill(rcPosX-mcPosX);
 866           fHistAsMcResyAntiLambda->Fill(rcPosY-mcPosY);
 867           fHistAsMcReszAntiLambda->Fill(rcPosZ-mcPosZ);
 868           fHistAsMcResrVsRadiusAntiLambda->Fill(rcPosR,rcPosR-mcPosR);
 869           fHistAsMcReszVsRadiusAntiLambda->Fill(rcPosR,rcPosZ-mcPosZ);
 870         }
 871         break;
 872
 873       case 1 :
 874         fHistMassLambdaMI->Fill(myAODv0->MassLambda());
 875         fHistMassAntiLambdaMI->Fill(myAODv0->MassAntiLambda());
 876         fHistMassVsRadiusLambdaMI->Fill(myAODv0->RadiusV0(),myAODv0->MassLambda());
 877         fHistMassVsRadiusAntiLambdaMI->Fill(myAODv0->RadiusV0(),myAODv0->MassAntiLambda());
 878         if(lCheckPIdLambda) fHistPidMcMassLambdaMI->Fill(myAODv0->MassLambda());
 879         else if (lCheckPIdAntiLambda) fHistPidMcMassAntiLambdaMI->Fill(myAODv0->MassAntiLambda());
 880         if(lCheckMcLambda) {
 881           fHistAsMcMassLambdaMI->Fill(myAODv0->MassLambda());
 882           fHistAsMcPtLambdaMI->Fill(pt);
 883           fHistAsMcMassVsRadiusLambdaMI->Fill(myAODv0->RadiusV0(),myAODv0->MassLambda());
 884           fHistAsMcResxLambdaMI->Fill(rcPosX-mcPosX);
 885           fHistAsMcResyLambdaMI->Fill(rcPosY-mcPosY);
 886           fHistAsMcReszLambdaMI->Fill(rcPosZ-mcPosZ);
 887           fHistAsMcResrVsRadiusLambdaMI->Fill(rcPosR,rcPosR-mcPosR);
 888           fHistAsMcReszVsRadiusLambdaMI->Fill(rcPosR,rcPosZ-mcPosZ);
 889         }
 890         else if(lCheckMcAntiLambda) {
 891           fHistAsMcMassAntiLambdaMI->Fill(myAODv0->MassAntiLambda());
 892           fHistAsMcPtAntiLambdaMI->Fill(pt);
 893           fHistAsMcMassVsRadiusAntiLambdaMI->Fill(myAODv0->RadiusV0(),myAODv0->MassAntiLambda());
 894           fHistAsMcResxAntiLambdaMI->Fill(rcPosX-mcPosX);
 895           fHistAsMcResyAntiLambdaMI->Fill(rcPosY-mcPosY);
 896           fHistAsMcReszAntiLambdaMI->Fill(rcPosZ-mcPosZ);
 897           fHistAsMcResrVsRadiusAntiLambdaMI->Fill(rcPosR,rcPosR-mcPosR);
 898           fHistAsMcReszVsRadiusAntiLambdaMI->Fill(rcPosR,rcPosZ-mcPosZ);
 899         }
 900         break;
 901
 902       }
 903     }
 904
 905   } // end V0 loop
 906
 907   if (primary) delete primary;
 908   if (myPosAodTrack) delete myPosAodTrack;
 909   if (myNegAodTrack) delete myNegAodTrack;
 910   //if (myAODVertex) delete myAODVertex; // don't !!!
 911   if (myAODv0) delete myAODv0;
 912
 913
 914   // Post output data
 915   PostData(0, fListHist);
 916 }
 917
 918 //________________________________________________________________________
 919 void AliAnalysisTaskESDStrangeMC::Terminate(Option_t *)
 920 {
 921   // Draw result to the screen
 922   // Called once at the end of the query
 923
 924   /*
 925   fHistPtMC = dynamic_cast<TH1F*> (((TList*)GetOutputData(0))->FindObject("h1PtMC"));
 926   if (!fHistPtMC) {
 927     Printf("ERROR: fHistPtMC not available");
 928     return;
 929   }
 930
 931   TCanvas *c1 = new TCanvas("AliAnalysisTaskESDStrangeMC","Pt MC",10,10,510,510);
 932   c1->cd(1)->SetLogy();
 933   fHistPtMC->DrawCopy("E");
 934 */
 935 }
 936
 937 //----------------------------------------------------------------------------
 938
 939 Double_t myRap(Double_t rE, Double_t rPz)
 940 {
 941   Double_t lRapidity = 1.e30;
 942   if(rPz && rE && (rPz != rE) && (1.+(2./((rE/rPz)-1.))>0))
 943     lRapidity = 0.5*log(1.+(2./((rE/rPz)-1.)));
 944   return lRapidity;
 945   }
 946
 947 //----------------------------------------------------------------------------
 948