PWG2/SPECTRA/LambdaK0PbPb/AliAnalysisTaskPerformanceStrange.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2009, 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 //-----------------------------------------------------------------
  17 //              AliAnalysisTaskPerformanceStrange class
  18 //    This task is for a performance study of V0 identification.
  19 //                It works with MC info and ESD tree.
  20 //                 Author: H.Ricaud, H.Ricaud@gsi.de
  21 //-----------------------------------------------------------------
  22
  23 #include <Riostream.h>
  24
  25 #include <stdio.h>
  26 #include <iostream>
  27 #include "TChain.h"
  28 #include "TTree.h"
  29 #include "TH1F.h"
  30 #include "TH2F.h"
  31 #include "TF1.h"
  32 #include "TList.h"
  33 #include "TMath.h"
  34 #include "TCanvas.h"
  35
  36 #include "AliAnalysisManager.h"
  37
  38 #include "AliPhysicsSelection.h"
  39 #include "AliBackgroundSelection.h"
  40
  41 #include "AliESDVertex.h"
  42 #include "AliESDEvent.h"
  43 #include "AliESDInputHandler.h"
  44 #include "AliESDtrack.h"
  45 #include "AliESDv0.h"
  46 #include "AliESDtrackCuts.h"
  47 #include "AliESDpid.h"
  48 #include "AliMultiplicity.h"
  49
  50 #include "AliAODEvent.h"
  51 #include "AliAODVertex.h"
  52 #include "AliAODTrack.h"
  53 #include "AliAODv0.h"
  54 #include "AliAODMCHeader.h"
  55 #include "AliAODInputHandler.h"
  56
  57 #include "AliAODMCParticle.h"
  58
  59 #include "AliMCEventHandler.h"
  60 #include "AliMCEvent.h"
  61 #include "AliStack.h"
  62 #include "AliGenEventHeader.h"
  63
  64 #include "AliLog.h"
  65
  66 #include "AliKFVertex.h"
  67 #include "AliVertexerTracks.h"
  68
  69 #include "AliAnalysisTaskPerformanceStrange.h"
  70 #include "AliAnalysisCentralitySelector.h"
  71
  72
  73 ClassImp(AliAnalysisTaskPerformanceStrange)
  74
  75
  76 //________________________________________________________________________
  77 AliAnalysisTaskPerformanceStrange::AliAnalysisTaskPerformanceStrange()
  78   : AliAnalysisTaskSE("TaskStrange"), fAnalysisMC(0), fAnalysisType("infoType"),  fCollidingSystems(0), fUsePID("infoPID"), fUseCut("infoCut"), fUseOnTheFly(0), fESD(0), fListHist(0), fCentrSelector(0),
  79
  80     // MC histograms  ---------------------------------------
  81     fHistMCPrimaryVertexX(0),
  82     fHistMCPrimaryVertexY(0),
  83     fHistMCPrimaryVertexZ(0),
  84
  85     fHistMCMultiplicityPrimary(0),
  86     fHistMCMultiplicityTracks(0),
  87
  88     fHistMCtracksProdRadiusK0s(0),
  89     fHistMCtracksProdRadiusLambda(0),
  90     fHistMCtracksProdRadiusAntiLambda(0),
  91
  92     fHistMCtracksDecayRadiusK0s(0),
  93     fHistMCtracksDecayRadiusLambda(0),
  94     fHistMCtracksDecayRadiusAntiLambda(0),
  95
  96     fHistMCPtAllK0s(0),
  97     fHistMCPtAllLambda(0),
  98     fHistMCPtAllAntiLambda(0),
  99
 100     fHistMCProdRadiusK0s(0),
 101     fHistMCProdRadiusLambda(0),
 102     fHistMCProdRadiusAntiLambda(0),
 103
 104     fHistMCRapK0s(0),
 105     fHistMCRapInPtRangeK0s(0),
 106     fHistMCRapLambda(0),
 107     fHistMCRapInPtRangeLambda(0),
 108     fHistMCRapAntiLambda(0),
 109     fHistMCRapInPtRangeAntiLambda(0),
 110     fHistMCRapXi(0),
 111     fHistMCRapInPtRangeXi(0),
 112     fHistMCRapPhi(0),
 113     fHistMCRapInPtRangePhi(0),
 114
 115     fHistMCPtK0s(0),
 116     fHistMCPtLambda(0),
 117
 118     fHistMCPtLambdaFromSigma(0),
 119     fHistMCPtAntiLambdaFromSigma(0),
 120
 121     fHistNTimesRecK0s(0),
 122     fHistNTimesRecLambda(0),
 123     fHistNTimesRecAntiLambda(0),
 124     fHistNTimesRecK0sVsPt(0),
 125     fHistNTimesRecLambdaVsPt(0),
 126     fHistNTimesRecAntiLambdaVsPt(0),
 127     // ------------------------------------------------------
 128
 129     // Reconstructed particle histograms  -------------------
 130     fHistNumberEvents(0),
 131     fHistTrackPerEvent(0),
 132     fHistTrackletPerEvent(0),
 133     fHistMCDaughterTrack(0),
 134
 135     fHistSPDPrimaryVertexZ(0),
 136
 137     fHistPrimaryVertexX(0),
 138     fHistPrimaryVertexY(0),
 139     fHistPrimaryVertexZ(0),
 140
 141     fHistPrimaryVertexResX(0),
 142     fHistPrimaryVertexResY(0),
 143     fHistPrimaryVertexResZ(0),
 144
 145     fHistPrimaryVertexPosXV0events(0),
 146     fHistPrimaryVertexPosYV0events(0),
 147     fHistPrimaryVertexPosZV0events(0),
 148
 149     fHistDaughterPt(0),
 150
 151     fHistDcaPosToPrimVertex(0),
 152     fHistDcaNegToPrimVertex(0),
 153     fHistDcaPosToPrimVertexZoom(0),
 154     fHistDcaNegToPrimVertexZoom(0),
 155     fHistRadiusV0(0),
 156     fHistDecayLengthV0(0),
 157     fHistDcaV0Daughters(0),
 158     fHistChi2(0),
 159     fHistCosPointAngle(0),
 160     fHistCosPointAngleZoom(0),
 161     fHistProdRadius(0),
 162
 163     fHistV0Multiplicity(0),
 164
 165     fHistChi2KFBeforeCutK0s(0),
 166     fHistChi2KFBeforeCutLambda(0),
 167     fHistChi2KFBeforeCutAntiLambda(0),
 168     fHistChi2KFAfterCutK0s(0),
 169     fHistChi2KFAfterCutLambda(0),
 170     fHistChi2KFAfterCutAntiLambda(0),
 171
 172     fHistMassK0(0),
 173     fHistMassLambda(0),
 174     fHistMassAntiLambda(0),
 175     fHistMassVsRadiusK0(0),
 176     fHistMassVsRadiusLambda(0),
 177     fHistMassVsRadiusAntiLambda(0),
 178
 179     fHistPtVsMassK0(0),
 180     fHistPtVsMassLambda(0),
 181     fHistArmenterosPodolanski(0),
 182     // ------------------------------------------------------
 183
 184
 185     // PID histograms  --------------------------------------
 186     fHistNsigmaPosPionAntiLambda(0),
 187     fHistNsigmaNegProtonAntiLambda(0),
 188     fHistNsigmaPosProtonLambda(0),
 189     fHistNsigmaNegPionLambda(0),
 190     fHistNsigmaPosPionK0(0),
 191     fHistNsigmaNegPionK0(0),
 192     // ------------------------------------------------------
 193
 194     // Associated particles ---------------------------------
 195     fHistAsMcRapK0(0),
 196     fHistAsMcRapLambda(0),
 197     fHistAsMcRapAntiLambda(0),
 198
 199     fHistAsMcPtK0(0),
 200     fHistAsMcPtLambda(0),
 201
 202     fHistAsMcPtZoomK0(0),
 203     fHistAsMcPtZoomLambda(0),
 204
 205     fHistAsMcProdRadiusK0(0),
 206     fHistAsMcProdRadiusLambda(0),
 207     fHistAsMcProdRadiusAntiLambda(0),
 208
 209     fHistAsMcProdRadiusXvsYK0s(0),
 210     fHistAsMcProdRadiusXvsYLambda(0),
 211     fHistAsMcProdRadiusXvsYAntiLambda(0),
 212
 213     fHistPidMcMassK0(0),
 214     fHistPidMcMassLambda(0),
 215     fHistPidMcMassAntiLambda(0),
 216     fHistAsMcMassK0(0),
 217     fHistAsMcMassLambda(0),
 218     fHistAsMcMassAntiLambda(0),
 219
 220     fHistAsMcPtVsMassK0(0),
 221     fHistAsMcPtVsMassLambda(0),
 222     fHistAsMcPtVsMassAntiLambda(0),
 223
 224     fHistAsMcMassVsRadiusK0(0),
 225     fHistAsMcMassVsRadiusLambda(0),
 226     fHistAsMcMassVsRadiusAntiLambda(0),
 227
 228     fHistAsMcResxK0(0),
 229     fHistAsMcResyK0(0),
 230     fHistAsMcReszK0(0),
 231
 232     fHistAsMcResrVsRadiusK0(0),
 233     fHistAsMcReszVsRadiusK0(0),
 234
 235     fHistAsMcResxLambda(0),
 236     fHistAsMcResyLambda(0),
 237     fHistAsMcReszLambda(0),
 238
 239     fHistAsMcResrVsRadiusLambda(0),
 240     fHistAsMcReszVsRadiusLambda(0),
 241
 242     fHistAsMcResxAntiLambda(0),
 243     fHistAsMcResyAntiLambda(0),
 244     fHistAsMcReszAntiLambda(0),
 245
 246     fHistAsMcResrVsRadiusAntiLambda(0),
 247     fHistAsMcReszVsRadiusAntiLambda(0),
 248
 249     fHistAsMcResPtK0(0),
 250     fHistAsMcResPtLambda(0),
 251     fHistAsMcResPtAntiLambda(0),
 252
 253     fHistAsMcResPtVsRapK0(0),
 254     fHistAsMcResPtVsRapLambda(0),
 255     fHistAsMcResPtVsRapAntiLambda(0),
 256     fHistAsMcResPtVsPtK0(0),
 257     fHistAsMcResPtVsPtLambda(0),
 258     fHistAsMcResPtVsPtAntiLambda(0),
 259
 260     fHistAsMcMotherPdgCodeK0s(0),
 261     fHistAsMcMotherPdgCodeLambda(0),
 262     fHistAsMcMotherPdgCodeAntiLambda(0),
 263
 264     fHistAsMcPtLambdaFromSigma(0),
 265     fHistAsMcPtAntiLambdaFromSigma(0),
 266     // ------------------------------------------------------
 267
 268     // Associated secondary particle histograms -------------
 269     fHistAsMcSecondaryPtVsRapK0s(0),
 270     fHistAsMcSecondaryPtVsRapLambda(0),
 271     fHistAsMcSecondaryPtVsRapAntiLambda(0),
 272
 273     fHistAsMcSecondaryProdRadiusK0s(0),
 274     fHistAsMcSecondaryProdRadiusLambda(0),
 275     fHistAsMcSecondaryProdRadiusAntiLambda(0),
 276
 277     fHistAsMcSecondaryProdRadiusXvsYK0s(0),
 278     fHistAsMcSecondaryProdRadiusXvsYLambda(0),
 279     fHistAsMcSecondaryProdRadiusXvsYAntiLambda(0),
 280
 281     fHistAsMcSecondaryMotherPdgCodeK0s(0),
 282     fHistAsMcSecondaryMotherPdgCodeLambda(0),
 283     fHistAsMcSecondaryMotherPdgCodeAntiLambda(0),
 284
 285     fHistAsMcSecondaryPtLambdaFromSigma(0),
 286     fHistAsMcSecondaryPtAntiLambdaFromSigma(0)
 287 {
 288   // Constructor
 289
 290   // New V0 cuts
 291 /*  fCuts[0]=33;    // max allowed chi2
 292   fCuts[1]=0.05;  // min allowed impact parameter for the 1st daughter
 293   fCuts[2]=0.05;  // min allowed impact parameter for the 2nd daughter
 294   fCuts[3]=0.5;   // max allowed DCA between the daughter tracks
 295   fCuts[4]=0.00;  // max allowed cosine of V0's pointing angle
 296   fCuts[5]=0.2;   // min radius of the fiducial volume
 297   fCuts[6]=100;   // max radius of the fiducial volume
 298 */
 299 }
 300
 301 //________________________________________________________________________
 302 AliAnalysisTaskPerformanceStrange::AliAnalysisTaskPerformanceStrange(const char *name)
 303   : AliAnalysisTaskSE(name), fAnalysisMC(0), fAnalysisType("infoType"), fCollidingSystems(0), fUsePID("infoPID"), fUseCut("infoCut"), fUseOnTheFly(0), fESD(0), fListHist(),fCentrSelector(0),
 304
 305     // MC histograms  ---------------------------------------
 306     fHistMCPrimaryVertexX(0),
 307     fHistMCPrimaryVertexY(0),
 308     fHistMCPrimaryVertexZ(0),
 309
 310     fHistMCMultiplicityPrimary(0),
 311     fHistMCMultiplicityTracks(0),
 312
 313     fHistMCtracksProdRadiusK0s(0),
 314     fHistMCtracksProdRadiusLambda(0),
 315     fHistMCtracksProdRadiusAntiLambda(0),
 316
 317     fHistMCtracksDecayRadiusK0s(0),
 318     fHistMCtracksDecayRadiusLambda(0),
 319     fHistMCtracksDecayRadiusAntiLambda(0),
 320
 321     fHistMCPtAllK0s(0),
 322     fHistMCPtAllLambda(0),
 323     fHistMCPtAllAntiLambda(0),
 324
 325     fHistMCProdRadiusK0s(0),
 326     fHistMCProdRadiusLambda(0),
 327     fHistMCProdRadiusAntiLambda(0),
 328
 329     fHistMCRapK0s(0),
 330     fHistMCRapInPtRangeK0s(0),
 331     fHistMCRapLambda(0),
 332     fHistMCRapInPtRangeLambda(0),
 333     fHistMCRapAntiLambda(0),
 334     fHistMCRapInPtRangeAntiLambda(0),
 335     fHistMCRapXi(0),
 336     fHistMCRapInPtRangeXi(0),
 337     fHistMCRapPhi(0),
 338     fHistMCRapInPtRangePhi(0),
 339
 340     fHistMCPtK0s(0),
 341     fHistMCPtLambda(0),
 342
 343     fHistMCPtLambdaFromSigma(0),
 344     fHistMCPtAntiLambdaFromSigma(0),
 345
 346     fHistNTimesRecK0s(0),
 347     fHistNTimesRecLambda(0),
 348     fHistNTimesRecAntiLambda(0),
 349     fHistNTimesRecK0sVsPt(0),
 350     fHistNTimesRecLambdaVsPt(0),
 351     fHistNTimesRecAntiLambdaVsPt(0),
 352     // ------------------------------------------------------
 353
 354     // Reconstructed particle histograms  -------------------
 355     fHistNumberEvents(0),
 356     fHistTrackPerEvent(0),
 357     fHistTrackletPerEvent(0),
 358     fHistMCDaughterTrack(0),
 359
 360     fHistSPDPrimaryVertexZ(0),
 361
 362     fHistPrimaryVertexX(0),
 363     fHistPrimaryVertexY(0),
 364     fHistPrimaryVertexZ(0),
 365
 366     fHistPrimaryVertexResX(0),
 367     fHistPrimaryVertexResY(0),
 368     fHistPrimaryVertexResZ(0),
 369
 370     fHistPrimaryVertexPosXV0events(0),
 371     fHistPrimaryVertexPosYV0events(0),
 372     fHistPrimaryVertexPosZV0events(0),
 373
 374     fHistDaughterPt(0),
 375
 376     fHistDcaPosToPrimVertex(0),
 377     fHistDcaNegToPrimVertex(0),
 378     fHistDcaPosToPrimVertexZoom(0),
 379     fHistDcaNegToPrimVertexZoom(0),
 380     fHistRadiusV0(0),
 381     fHistDecayLengthV0(0),
 382     fHistDcaV0Daughters(0),
 383     fHistChi2(0),
 384     fHistCosPointAngle(0),
 385     fHistCosPointAngleZoom(0),
 386     fHistProdRadius(0),
 387
 388     fHistV0Multiplicity(0),
 389
 390     fHistChi2KFBeforeCutK0s(0),
 391     fHistChi2KFBeforeCutLambda(0),
 392     fHistChi2KFBeforeCutAntiLambda(0),
 393     fHistChi2KFAfterCutK0s(0),
 394     fHistChi2KFAfterCutLambda(0),
 395     fHistChi2KFAfterCutAntiLambda(0),
 396
 397     fHistMassK0(0),
 398     fHistMassLambda(0),
 399     fHistMassAntiLambda(0),
 400     fHistMassVsRadiusK0(0),
 401     fHistMassVsRadiusLambda(0),
 402     fHistMassVsRadiusAntiLambda(0),
 403
 404     fHistPtVsMassK0(0),
 405     fHistPtVsMassLambda(0),
 406     fHistArmenterosPodolanski(0),
 407     // ------------------------------------------------------
 408
 409
 410     // PID histograms  --------------------------------------
 411     fHistNsigmaPosPionAntiLambda(0),
 412     fHistNsigmaNegProtonAntiLambda(0),
 413     fHistNsigmaPosProtonLambda(0),
 414     fHistNsigmaNegPionLambda(0),
 415     fHistNsigmaPosPionK0(0),
 416     fHistNsigmaNegPionK0(0),
 417     // ------------------------------------------------------
 418
 419     // Associated particles ---------------------------------
 420     fHistAsMcRapK0(0),
 421     fHistAsMcRapLambda(0),
 422     fHistAsMcRapAntiLambda(0),
 423
 424     fHistAsMcPtK0(0),
 425     fHistAsMcPtLambda(0),
 426
 427     fHistAsMcPtZoomK0(0),
 428     fHistAsMcPtZoomLambda(0),
 429
 430     fHistAsMcProdRadiusK0(0),
 431     fHistAsMcProdRadiusLambda(0),
 432     fHistAsMcProdRadiusAntiLambda(0),
 433
 434     fHistAsMcProdRadiusXvsYK0s(0),
 435     fHistAsMcProdRadiusXvsYLambda(0),
 436     fHistAsMcProdRadiusXvsYAntiLambda(0),
 437
 438     fHistPidMcMassK0(0),
 439     fHistPidMcMassLambda(0),
 440     fHistPidMcMassAntiLambda(0),
 441     fHistAsMcMassK0(0),
 442     fHistAsMcMassLambda(0),
 443     fHistAsMcMassAntiLambda(0),
 444
 445     fHistAsMcPtVsMassK0(0),
 446     fHistAsMcPtVsMassLambda(0),
 447     fHistAsMcPtVsMassAntiLambda(0),
 448
 449     fHistAsMcMassVsRadiusK0(0),
 450     fHistAsMcMassVsRadiusLambda(0),
 451     fHistAsMcMassVsRadiusAntiLambda(0),
 452
 453     fHistAsMcResxK0(0),
 454     fHistAsMcResyK0(0),
 455     fHistAsMcReszK0(0),
 456
 457     fHistAsMcResrVsRadiusK0(0),
 458     fHistAsMcReszVsRadiusK0(0),
 459
 460     fHistAsMcResxLambda(0),
 461     fHistAsMcResyLambda(0),
 462     fHistAsMcReszLambda(0),
 463
 464     fHistAsMcResrVsRadiusLambda(0),
 465     fHistAsMcReszVsRadiusLambda(0),
 466
 467     fHistAsMcResxAntiLambda(0),
 468     fHistAsMcResyAntiLambda(0),
 469     fHistAsMcReszAntiLambda(0),
 470
 471     fHistAsMcResrVsRadiusAntiLambda(0),
 472     fHistAsMcReszVsRadiusAntiLambda(0),
 473
 474     fHistAsMcResPtK0(0),
 475     fHistAsMcResPtLambda(0),
 476     fHistAsMcResPtAntiLambda(0),
 477
 478     fHistAsMcResPtVsRapK0(0),
 479     fHistAsMcResPtVsRapLambda(0),
 480     fHistAsMcResPtVsRapAntiLambda(0),
 481     fHistAsMcResPtVsPtK0(0),
 482     fHistAsMcResPtVsPtLambda(0),
 483     fHistAsMcResPtVsPtAntiLambda(0),
 484
 485     fHistAsMcMotherPdgCodeK0s(0),
 486     fHistAsMcMotherPdgCodeLambda(0),
 487     fHistAsMcMotherPdgCodeAntiLambda(0),
 488
 489     fHistAsMcPtLambdaFromSigma(0),
 490     fHistAsMcPtAntiLambdaFromSigma(0),
 491     // ------------------------------------------------------
 492
 493     // Associated secondary particle histograms -------------
 494     fHistAsMcSecondaryPtVsRapK0s(0),
 495     fHistAsMcSecondaryPtVsRapLambda(0),
 496     fHistAsMcSecondaryPtVsRapAntiLambda(0),
 497
 498     fHistAsMcSecondaryProdRadiusK0s(0),
 499     fHistAsMcSecondaryProdRadiusLambda(0),
 500     fHistAsMcSecondaryProdRadiusAntiLambda(0),
 501
 502     fHistAsMcSecondaryProdRadiusXvsYK0s(0),
 503     fHistAsMcSecondaryProdRadiusXvsYLambda(0),
 504     fHistAsMcSecondaryProdRadiusXvsYAntiLambda(0),
 505
 506     fHistAsMcSecondaryMotherPdgCodeK0s(0),
 507     fHistAsMcSecondaryMotherPdgCodeLambda(0),
 508     fHistAsMcSecondaryMotherPdgCodeAntiLambda(0),
 509
 510     fHistAsMcSecondaryPtLambdaFromSigma(0),
 511     fHistAsMcSecondaryPtAntiLambdaFromSigma(0)
 512 {
 513   // Constructor
 514
 515   //New V0 cuts
 516 /*  fCuts[0]=33;    // max allowed chi2
 517   fCuts[1]=0.05;  // min allowed impact parameter for the 1st daughter
 518   fCuts[2]=0.05;  // min allowed impact parameter for the 2nd daughter
 519   fCuts[3]=0.5;   // max allowed DCA between the daughter tracks
 520   fCuts[4]=0.00;  // max allowed cosine of V0's pointing angle
 521   fCuts[5]=0.2;   // min radius of the fiducial volume
 522   fCuts[6]=100;   // max radius of the fiducial volume
 523 */
 524   // Define output slots only here
 525   // Output slot #1 writes into a TList container
 526   DefineOutput(1, TList::Class());
 527   DefineOutput(2, AliAnalysisCentralitySelector::Class());
 528
 529 }
 530
 531 //________________________________________________________________________
 532 void AliAnalysisTaskPerformanceStrange::UserCreateOutputObjects()
 533 {
 534   //******************
 535   // Create histograms
 536   //*******************
 537   fListHist = new TList();
 538   fListHist->SetOwner();
 539
 540   // Bo: tbd: condition before allocation (i.e. if (!fHistMCMultiplicityPrimary){...} for each histo...
 541
 542   //***************
 543   // MC histograms
 544   //***************
 545
 546   // Primary Vertex X,Y,Z:
 547   fHistMCPrimaryVertexX          = new TH1F("h1MCPrimaryVertexX", "MC Primary Vertex Position X;Primary Vertex Position X (cm);Events",100,-0.5,0.5);
 548   fListHist->Add(fHistMCPrimaryVertexX);
 549
 550   fHistMCPrimaryVertexY          = new TH1F("h1MCPrimaryVertexY", "MC Primary Vertex Position Y;Primary Vertex Position Y (cm);Events",100,-0.5,0.5);
 551   fListHist->Add(fHistMCPrimaryVertexY);
 552
 553   fHistMCPrimaryVertexZ          = new TH1F("h1MCPrimaryVertexZ", "MC Primary Vertex Position Z;Primary Vertex Position Z (cm);Events",200,-20,20);
 554   fListHist->Add(fHistMCPrimaryVertexZ);
 555
 556   // Multiplicity:
 557   fHistMCMultiplicityPrimary           = new TH1F("h1MCMultiplicityPrimary", "MC Primary Particles;NPrimary;Count", 201, -0.5, 200.5);
 558   fListHist->Add(fHistMCMultiplicityPrimary);
 559
 560   fHistMCMultiplicityTracks            = new TH1F("h1MCMultiplicityTracks", "MC Tracks;Ntracks;Count", 201, -0.5, 200.5);
 561   fListHist->Add(fHistMCMultiplicityTracks);
 562
 563   // Production Radius of non-primary particles:
 564   fHistMCtracksProdRadiusK0s           = new TH2F("h2MCtracksProdRadiusK0s","Non-primary MC K^{0} Production Radius;x (cm); y (cm)",200,-50,50,200,-50,50);
 565   fListHist->Add(fHistMCtracksProdRadiusK0s);
 566
 567   fHistMCtracksProdRadiusLambda        = new TH2F("h2MCtracksProdRadiusLambda","Non-primary MC #Lambda^{0} Production Radius;x (cm); y (cm)",200,-50,50,200,-50,50);
 568   fListHist->Add(fHistMCtracksProdRadiusLambda);
 569
 570   fHistMCtracksProdRadiusAntiLambda    = new TH2F("h2MCtracksProdRadiusAntiLambda","Non-primary MC #bar{#Lambda}^{0} Production Radius;x (cm); y (cm)",200,-50,50,200,-50,50);
 571   fListHist->Add(fHistMCtracksProdRadiusAntiLambda);
 572
 573   // Decay Radius of non-primary particles:
 574   fHistMCtracksDecayRadiusK0s          = new TH1F("h1MCtracksDecayRadiusK0s","Non-primary MC K^{0} Decay Radius;r (cm)",101,-1,100);
 575   fListHist->Add(fHistMCtracksDecayRadiusK0s);
 576
 577   fHistMCtracksDecayRadiusLambda       = new TH1F("h1MCtracksDecayRadiusLambda","Non-primary MC #Lambda^{0} Decay Radius;r (cm)",101,-1,100);
 578   fListHist->Add(fHistMCtracksDecayRadiusLambda);
 579
 580   fHistMCtracksDecayRadiusAntiLambda   = new TH1F("h1MCtracksDecayRadiusAntiLambda","Non-primary #bar{#Lambda}^{0} Decay Radius;r (cm)",100,1,101);
 581   fListHist->Add(fHistMCtracksDecayRadiusAntiLambda);
 582
 583   // Pt Distribution of non-primary particles:
 584   fHistMCPtAllK0s                      = new TH1F("h1MCPtAllK0s", "Non-primary MC K^{0};p_{t} (GeV/c);Counts",240,0,12);
 585   fListHist->Add(fHistMCPtAllK0s);
 586
 587   fHistMCPtAllLambda                   = new TH1F("h1MCPtAllLambda", "Non-primary MC #Lambda^{0};p_{t} (GeV/c);Counts",240,0,12);
 588   fListHist->Add(fHistMCPtAllLambda);
 589
 590   fHistMCPtAllAntiLambda               = new TH1F("h1MCPtAllAntiLambda", "Non-primary MC #bar{#Lambda}^{0};p_{t} (GeV/c);Counts",240,0,12);
 591   fListHist->Add(fHistMCPtAllAntiLambda);
 592
 593   // Production Radius
 594   fHistMCProdRadiusK0s                 = new TH1F("h1MCProdRadiusK0s", "MC K^{0} Production Radius;r (cm);Count", 400, -2, 2);
 595   fListHist->Add(fHistMCProdRadiusK0s);
 596
 597   fHistMCProdRadiusLambda              = new TH1F("h1MCProdRadiusLambda", "MC #Lambda^{0} Production Radius;r (cm);Count", 400, -2, 2);
 598   fListHist->Add(fHistMCProdRadiusLambda);
 599
 600    fHistMCProdRadiusAntiLambda         = new TH1F("h1MCProdRadiusAntiLambda", "MC #bar{#Lambda}^{0} Production Radius;r (cm);Count", 400, -2, 2);
 601   fListHist->Add(fHistMCProdRadiusAntiLambda);
 602
 603
 604   // Rapidity distribution:
 605   fHistMCRapK0s                 = new TH1F("h1MCRapK0s", "K^{0};y",160,-4,4);
 606   fListHist->Add(fHistMCRapK0s);
 607
 608   fHistMCRapInPtRangeK0s        = new TH1F("h1MCRapInPtRangeK0s", "K^{0};y",160,-4,4);
 609   fListHist->Add(fHistMCRapInPtRangeK0s);
 610
 611   fHistMCRapLambda              = new TH1F("h1MCRapLambda", "#Lambda;y",160,-4,4);
 612   fListHist->Add(fHistMCRapLambda);
 613
 614   fHistMCRapInPtRangeLambda     = new TH1F("h1MCRapInPtRangeLambda", "#Lambda;y",160,-4,4);
 615   fListHist->Add(fHistMCRapInPtRangeLambda);
 616
 617   fHistMCRapAntiLambda          = new TH1F("h1MCRapAntiLambda", "#bar{#Lambda};y",160,-4,4);
 618   fListHist->Add(fHistMCRapAntiLambda);
 619
 620   fHistMCRapInPtRangeAntiLambda = new TH1F("h1MCRapInPtRangeAntiLambda", "#bar{#Lambda};y",160,-4,4);
 621   fListHist->Add(fHistMCRapInPtRangeAntiLambda);
 622
 623   fHistMCRapXi                  = new TH1F("h1MCRapXi", "#Xi;y",160,-4,4);
 624   fListHist->Add(fHistMCRapXi);
 625
 626   fHistMCRapInPtRangeXi         = new TH1F("h1MCRapInPtRangeXi", "#Xi;y",160,-4,4);
 627   fListHist->Add(fHistMCRapInPtRangeXi);
 628
 629   fHistMCRapPhi                  = new TH1F("h1MCRapPhi", "#phi;y",160,-4,4);
 630   fListHist->Add(fHistMCRapPhi);
 631
 632   fHistMCRapInPtRangePhi         = new TH1F("h1MCRapInPtRangePhi", "#phi;y",160,-4,4);
 633   fListHist->Add(fHistMCRapInPtRangePhi);
 634
 635
 636   // Pt distribution:
 637   fHistMCPtK0s               = new TH1F("h1MCPtK0s", "K^{0};p_{t} (GeV/c)",240,0,12);
 638   fListHist->Add(fHistMCPtK0s);
 639
 640   fHistMCPtLambda            = new TH1F("h1MCPtLambda", "#Lambda^{0};p_{t} (GeV/c)",240,0,12);
 641   fListHist->Add(fHistMCPtLambda);
 642
 643   // Pt distribution of Lambda coming from Sigma decay:
 644   fHistMCPtLambdaFromSigma      = new TH1F("h1MCPtLambdaFromSigma", "#Lambda^{0};p_{t} (GeV/c)",240,0,12);
 645   fListHist->Add(fHistMCPtLambdaFromSigma);
 646
 647   fHistMCPtAntiLambdaFromSigma  = new TH1F("h1MCPtAntiLambdaFromSigma", "#Lambda^{0};p_{t} (GeV/c)",240,0,12);
 648   fListHist->Add(fHistMCPtAntiLambdaFromSigma);
 649
 650   // Multiple reconstruction studies:
 651   fHistNTimesRecK0s             = new TH1F("h1NTimesRecK0s","number of times a K0s is reconstructed in -1<y<1;number of times;counts",500,-0.5,4.5);
 652   fListHist->Add(fHistNTimesRecK0s);
 653
 654   fHistNTimesRecLambda          = new TH1F("h1NTimesRecLambda","number of times a Lambda is reconstructed in -1<y<1;number of times;counts",500,-0.5,4.5);
 655   fListHist->Add(fHistNTimesRecLambda);
 656
 657   fHistNTimesRecAntiLambda      = new TH1F("h1NTimesRecAntiLambda","number of times an AntiLambda is reconstructed in -1<y<1;number of times;counts",500,-0.5,4.5);
 658   fListHist->Add(fHistNTimesRecAntiLambda);
 659
 660   fHistNTimesRecK0sVsPt         = new TH2F("h2NTimesRecK0sVsPt","NTimes versus Pt, K^{0} in -1<y<1;p_{t} (GeV/c);number of times",75,0,15,5,-0.5,4.5);
 661   fListHist->Add(fHistNTimesRecK0sVsPt);
 662
 663   fHistNTimesRecLambdaVsPt      = new TH2F("h2NTimesRecLambdaVsPt","NTimes versus Pt, #Lambda^{0} in -1<y<1;p_{t} (GeV/c);number of times",75,0,15,5,-0.5,4.5);
 664   fListHist->Add(fHistNTimesRecLambdaVsPt);
 665
 666   fHistNTimesRecAntiLambdaVsPt  = new TH2F("h2NTimesRecAntiLambdaVsPt","NTimes versus Pt, #bar{#Lambda}^{0} in -1<y<1;p_{t} (GeV/c);number of times",75,0,15,5,-0.5,4.5);
 667   fListHist->Add(fHistNTimesRecAntiLambdaVsPt);
 668
 669   //***********************************
 670   // Reconstructed particles histograms
 671   //***********************************
 672
 673   // Number of events:
 674   fHistNumberEvents           = new TH1F("h1NumberEvents", "Number of events; index;Number of Events",10,0,10);
 675   fListHist->Add(fHistNumberEvents);
 676
 677   // Multiplicity:
 678   fHistTrackPerEvent           = new TH1F("h1TrackPerEvent", "Tracks per event;Number of Tracks;Number of Events",20000,0,20000);
 679   fListHist->Add(fHistTrackPerEvent);
 680
 681   fHistTrackletPerEvent       = new TH1F("h1TrackletPerEvent", "Number of tracklets;Number of tracklets per events;Number of events",1000,0,1000);
 682   fListHist->Add(fHistTrackletPerEvent);
 683
 684   fHistMCDaughterTrack         = new TH1F("h1MCDaughterTrack","Distribution of mc id for daughters;id tags;Counts",15,0,15);
 685   fListHist->Add(fHistMCDaughterTrack);
 686
 687   // Primary Vertex:
 688   fHistSPDPrimaryVertexZ          = new TH1F("h1SPDPrimaryVertexZ", "SPD Primary Vertex Position Z;Primary Vertex Position Z (cm);Events",200,-20,20);
 689   fListHist->Add(fHistSPDPrimaryVertexZ);
 690
 691   fHistPrimaryVertexX          = new TH1F("h1PrimaryVertexX", "Primary Vertex Position X;Primary Vertex Position X (cm);Events",100,-0.5,0.5);
 692   fListHist->Add(fHistPrimaryVertexX);
 693
 694   fHistPrimaryVertexY          = new TH1F("h1PrimaryVertexY", "Primary Vertex Position Y;Primary Vertex Position Y (cm);Events",100,-0.5,0.5);
 695   fListHist->Add(fHistPrimaryVertexY);
 696
 697   fHistPrimaryVertexZ          = new TH1F("h1PrimaryVertexZ", "Primary Vertex Position Z;Primary Vertex Position Z (cm);Events",200,-20,20);
 698   fListHist->Add(fHistPrimaryVertexZ);
 699
 700
 701   // Primary vertex resolution:
 702   fHistPrimaryVertexResX          = new TH1F("h1PrimaryVertexResX", "Primary Vertex Resolution X;Primary Vertex Resolution X (cm);Events",100,-0.25,0.25);
 703   fListHist->Add(fHistPrimaryVertexResX);
 704
 705   fHistPrimaryVertexResY          = new TH1F("h1PrimaryVertexResY", "Primary Vertex Resolution Y;Primary Vertex Resolution Y (cm);Events",100,-0.25,0.25);
 706   fListHist->Add(fHistPrimaryVertexResY);
 707
 708   fHistPrimaryVertexResZ          = new TH1F("h1PrimaryVertexResZ", "Primary Vertex Resolution Z;Primary Vertex Resolution Z (cm);Events",200,-0.25,0.25);
 709   fListHist->Add(fHistPrimaryVertexResZ);
 710
 711
 712   // Primary Vertex in events with V0 candidates:
 713   fHistPrimaryVertexPosXV0events       = new TH1F("h1PrimaryVertexPosXV0events", "Primary Vertex Position X;Primary Vertex Position X (cm);Events",100,-0.5,0.5);
 714   fListHist->Add(fHistPrimaryVertexPosXV0events);
 715   fHistPrimaryVertexPosYV0events       = new TH1F("h1PrimaryVertexPosYV0events", "Primary Vertex Position Y;Primary Vertex Position Y (cm);Events",100,-0.5,0.5);
 716   fListHist->Add(fHistPrimaryVertexPosYV0events);
 717   fHistPrimaryVertexPosZV0events       = new TH1F("h1PrimaryVertexPosZV0events", "Primary Vertex Position Z;Primary Vertex Position Z (cm);Events",200,-20.0,20.0);
 718   fListHist->Add(fHistPrimaryVertexPosZV0events);
 719
 720   // Daughters Pt:
 721   fHistDaughterPt              = new TH2F("h2DaughterPt", "Daughter Pt;Positive Daughter Pt; Negative Daughter Pt",200,0,2,200,0,2);
 722   fListHist->Add(fHistDaughterPt);
 723
 724   // Cut checks:
 725   fHistDcaPosToPrimVertex      = new TH2F("h2DcaPosToPrimVertex", "Positive V0 daughter;dca(cm);Status",500,0,5,2,-0.5,1.5);
 726   fListHist->Add(fHistDcaPosToPrimVertex);
 727
 728   fHistDcaNegToPrimVertex      = new TH2F("h2DcaNegToPrimVertex", "Negative V0 daughter;dca(cm);Status",500,0,5,2,-0.5,1.5);
 729   fListHist->Add(fHistDcaNegToPrimVertex);
 730
 731   fHistDcaPosToPrimVertexZoom  = new TH2F("h2DcaPosToPrimVertexZoom", "Positive V0 daughter;dca(cm);Status",100,0,0.1,2,-0.5,1.5);
 732   fListHist->Add(fHistDcaPosToPrimVertexZoom);
 733
 734   fHistDcaNegToPrimVertexZoom  = new TH2F("h2DcaNegToPrimVertexZoom", "Negative V0 daughter;dca(cm);Status",100,0,0.1,2,-0.5,1.5);
 735   fListHist->Add(fHistDcaNegToPrimVertexZoom);
 736
 737   fHistRadiusV0                = new TH2F("h2RadiusV0", "Radius;Radius(cm);Status",5000,0,500,2,-0.5,1.5);
 738   fListHist->Add(fHistRadiusV0);
 739
 740   fHistDecayLengthV0           = new TH2F("h2DecayLengthV0", "V0s decay Length;decay length(cm);Status", 240, 0, 120,2,-0.5,1.5);
 741   fListHist->Add(fHistDecayLengthV0);
 742
 743   fHistDcaV0Daughters          = new TH2F("h2DcaV0Daughters", "DCA between daughters;dca(cm);Status", 400, 0, 4,2,-0.5,1.5);
 744   fListHist->Add(fHistDcaV0Daughters);
 745
 746   fHistChi2                    = new TH2F("h2Chi2", "V0s chi2;chi2;Status", 33, 0, 33,2,-0.5,1.5);
 747   fListHist->Add(fHistChi2);
 748
 749   fHistCosPointAngle           = new TH2F("h2CosPointAngle", "Cosine of V0's pointing angle", 100,0,1,2,-0.5,1.5);
 750   fListHist->Add(fHistCosPointAngle);
 751
 752   fHistCosPointAngleZoom       = new TH2F("h2CosPointAngleZoom", "Cosine of V0's pointing angle", 1000,0.9,1,2,-0.5,1.5);
 753   fListHist->Add(fHistCosPointAngleZoom);
 754
 755   fHistProdRadius              = new TH2F("h2ProdRadius", "Production position;x (cm);y (cm)", 100,-50,50,100,-50,50);
 756   fListHist->Add(fHistProdRadius);
 757
 758   // V0 Multiplicity:
 759   if (!fHistV0Multiplicity) {
 760     if (fCollidingSystems)
 761       fHistV0Multiplicity = new TH1F("fHistV0Multiplicity", "Multiplicity distribution;Number of Offline V0s;Events", 200, 0, 40000);
 762     else
 763       fHistV0Multiplicity = new TH1F("fHistV0Multiplicity", "Multiplicity distribution;Number of Offline V0s;Events", 10, 0, 10);
 764     fListHist->Add(fHistV0Multiplicity);
 765   }
 766
 767   // Kalman Filter Chi2:
 768   fHistChi2KFBeforeCutK0s               = new TH2F("h1Chi2KFBeforeCutK0s", "K^{0}  candidates;#Chi^{2});Counts", 250, 0, 50, 2,-0.5,1.5);
 769   fListHist->Add(fHistChi2KFBeforeCutK0s);
 770   fHistChi2KFBeforeCutLambda            = new TH2F("h1Chi2KFBeforeCutLambda", "#Lambda^{0}  candidates;#Chi^{2};Counts", 250, 0, 50, 2,-0.5,1.5);
 771   fListHist->Add(fHistChi2KFBeforeCutLambda);
 772   fHistChi2KFBeforeCutAntiLambda        = new TH2F("h1Chi2KFBeforeCutAntiLambda", "#bar{#Lambda}^{0}  candidates;#Chi^{2};Counts", 250, 0, 50, 2,-0.5,1.5);
 773   fListHist->Add(fHistChi2KFBeforeCutAntiLambda);
 774
 775   fHistChi2KFAfterCutK0s               = new TH2F("h1Chi2KFAfterCutK0s", "K^{0}  candidates;#Chi^{2});Counts", 250, 0, 50, 2,-0.5,1.5);
 776   fListHist->Add(fHistChi2KFAfterCutK0s);
 777   fHistChi2KFAfterCutLambda            = new TH2F("h1Chi2KFAfterCutLambda", "#Lambda^{0}  candidates;#Chi^{2};Counts", 250, 0, 50, 2,-0.5,1.5);
 778   fListHist->Add(fHistChi2KFAfterCutLambda);
 779   fHistChi2KFAfterCutAntiLambda        = new TH2F("h1Chi2KFAfterCutAntiLambda", "#bar{#Lambda}^{0} candidates;#Chi^{2};Counts", 250, 0, 50, 2,-0.5,1.5);
 780   fListHist->Add(fHistChi2KFAfterCutAntiLambda);
 781
 782   // Invariant mass:
 783   fHistMassK0                   = new TH1F("h1MassK0", "K^{0} candidates;M(#pi^{+}#pi^{-}) (GeV/c^{2});Counts", 100, 0.4, 0.6);
 784   fListHist->Add(fHistMassK0);
 785
 786   fHistMassLambda               = new TH1F("h1MassLambda", "#Lambda^{0} candidates;M(p#pi^{-}) (GeV/c^{2});Counts", 75, 1.05, 1.2);
 787   fListHist->Add(fHistMassLambda);
 788
 789   fHistMassAntiLambda           = new TH1F("h1MassAntiLambda", "#bar{#Lambda}^{0} candidates;M(#bar{p}#pi^{+}) (GeV/c^{2});Counts", 75, 1.05, 1.2);
 790   fListHist->Add(fHistMassAntiLambda);
 791
 792   // Invariant mass vs. radius:
 793   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};
 794   Int_t lNbinRadius        = 9;
 795   Int_t lNbinInvMassLambda = 300;
 796
 797   fHistMassVsRadiusK0           = new TH2F("h2MassVsRadiusK0", "K^{0} candidates;radius (cm);M(#pi^{+}#pi^{-}) (GeV/c^{2})",lNbinRadius,radius, 200, 0.4, 0.6);
 798   fListHist->Add(fHistMassVsRadiusK0);
 799
 800   fHistMassVsRadiusLambda       = new TH2F("h2MassVsRadiusLambda", "#Lambda candidates;radius (cm);M(p#pi^{-}) (GeV/c^{2})",lNbinRadius,radius, 140, 1.06, 1.2);
 801   fListHist->Add(fHistMassVsRadiusLambda);
 802
 803   fHistMassVsRadiusAntiLambda   = new TH2F("h2MassVsRadiusAntiLambda", "#bar{#Lambda} candidates;radius (cm);M(#bar{p}#pi^{+}) (GeV/c^{2})",lNbinRadius,radius, 140, 1.06, 1.2);
 804   fListHist->Add(fHistMassVsRadiusAntiLambda);
 805
 806   // Pt vs. mass:
 807   fHistPtVsMassK0               = new TH2F("h2PtVsMassK0","K^{0} candidates;M(#pi^{+}#pi^{-}) (GeV/c^{2});p_{t} (GeV/c)",200, 0.4, 0.6,240,0,12);
 808   fListHist->Add(fHistPtVsMassK0);
 809
 810   fHistPtVsMassLambda           = new TH2F("h2PtVsMassLambda","#Lambda^{0} candidates;M(p#pi^{-}) (GeV/c^{2});p_{t} (GeV/c)",140, 1.06, 1.2,240,0,12);
 811   fListHist->Add(fHistPtVsMassLambda);
 812
 813   fHistArmenterosPodolanski     = new TH2F("h2ArmenterosPodolanski","Armenteros-Podolanski phase space;#alpha;p_{t} arm",100,-1.0,1.0,50,0,0.5);
 814   fListHist->Add(fHistArmenterosPodolanski);
 815
 816
 817   // PID histograms:
 818   fHistNsigmaPosPionAntiLambda   = new TH1F("h1NsigmaPosPionAntiLambda", "Positive daughter of Antilambda;NsigmaPion;Counts",25,0,5);
 819   fListHist->Add(fHistNsigmaPosPionAntiLambda);
 820
 821   fHistNsigmaNegProtonAntiLambda = new TH1F("h1NsigmaNegProtonAntiLambda", "Negative daughter of Antilambda;NsigmaProton;Counts",25,0,5);
 822   fListHist->Add(fHistNsigmaNegProtonAntiLambda);
 823
 824   fHistNsigmaPosProtonLambda     = new TH1F("h1NsigmaPosProtonLambda", "Positive daughter of Lambda;NsigmaProton;Counts",25,0,5);
 825   fListHist->Add(fHistNsigmaPosProtonLambda);
 826
 827   fHistNsigmaNegPionLambda       = new TH1F("h1NsigmaNegPionLambda", "Negative daughter of Lambda;NsigmaPion;Counts",25,0,5);
 828   fListHist->Add(fHistNsigmaNegPionLambda);
 829
 830   fHistNsigmaPosPionK0           = new TH1F("h1NsigmaPosPionK0", "Positive daughter of K0s;NsigmaPion;Counts",25,0,5);
 831   fListHist->Add(fHistNsigmaPosPionK0);
 832
 833   fHistNsigmaNegPionK0           = new TH1F("h1NsigmaNegPionK0", "Negative daughter of K0s;NsigmaPion;Counts",25,0,5);
 834   fListHist->Add(fHistNsigmaNegPionK0);
 835
 836   //********************************
 837   // Associated particle histograms
 838   //********************************
 839
 840   // Rapidity distribution:
 841   fHistAsMcRapK0                = new TH1F("h1AsMcRapK0", "K^{0} associated;eta;Counts", 60, -1.5, 1.5);
 842   fListHist->Add(fHistAsMcRapK0);
 843
 844   fHistAsMcRapLambda            = new TH1F("h1AsMcRapLambda", "#Lambda^{0} associated;eta;Counts", 60, -1.5, 1.5);
 845   fListHist->Add(fHistAsMcRapLambda);
 846
 847   fHistAsMcRapAntiLambda        = new TH1F("h1AsMcRapAntiLambda", "#bar{#Lambda}^{0} associated;eta;Counts", 60, -1.5, 1.5);
 848   fListHist->Add(fHistAsMcRapAntiLambda);
 849
 850   //Pt distribution:
 851   fHistAsMcPtK0                = new TH1F("h1AsMcPtK0", "K^{0} associated;p_{t} (GeV/c);Counts", 240,0,12);
 852   fListHist->Add(fHistAsMcPtK0);
 853
 854   fHistAsMcPtLambda            = new TH1F("h1AsMcPtLambda", "#Lambda^{0} associated;p_{t} (GeV/c);Counts", 240,0,12);
 855   fListHist->Add(fHistAsMcPtLambda);
 856
 857   fHistAsMcPtZoomK0            = new TH1F("h1AsMcPtZoomK0", "K^{0} candidates in -1 <y<1;p_{t} (GeV/c);Counts",20,0,1);
 858   fListHist->Add(fHistAsMcPtZoomK0);
 859
 860   fHistAsMcPtZoomLambda        = new TH1F("h1AsMcPtZoomLambda", "#Lambda^{0} candidates in -1 <y<1;p_{t} (GeV/c);Counts",20,0,1);
 861   fListHist->Add(fHistAsMcPtZoomLambda);
 862
 863   // Radius distribution:
 864   fHistAsMcProdRadiusK0               = new TH1F("h1AsMcProdRadiusK0", "K^{0} associated;r (cm);Counts", 500, 0, 100);
 865   fListHist->Add(fHistAsMcProdRadiusK0);
 866
 867   fHistAsMcProdRadiusLambda           = new TH1F("h1AsMcProdRadiusLambda", "#Lambda^{0} associated;r (cm);Counts", 500, 0, 100);
 868   fListHist->Add(fHistAsMcProdRadiusLambda);
 869
 870   fHistAsMcProdRadiusAntiLambda       = new TH1F("h1AsMcProdRadiusAntiLambda", "#bar{#Lambda}^{0} associated;r (cm);Counts", 500, 0, 100);
 871   fListHist->Add(fHistAsMcProdRadiusAntiLambda);
 872
 873   // Radius distribution vs. rapidity:
 874   fHistAsMcProdRadiusXvsYK0s          = new TH2F("h2AsMcProdRadiusXvsYK0s","Associated Secondary K^{0} Production Radius;x (cm); y (cm)",200,-50,50,200,-50,50);
 875   fListHist->Add(fHistAsMcProdRadiusXvsYK0s);
 876
 877   fHistAsMcProdRadiusXvsYLambda       = new TH2F("h2AsMcProdRadiusXvsYLambda","Associated Secondary #Lambda^{0} Production Radius;x (cm); y (cm)",200,-50,50,200,-50,50);
 878   fListHist->Add(fHistAsMcProdRadiusXvsYLambda);
 879
 880   fHistAsMcProdRadiusXvsYAntiLambda   = new TH2F("h2AsMcProdRadiusXvsYAntiLambda","Associated Secondary #bar{#Lambda}^{0} Production Radius;x (cm); y (cm)",200,-50,50,200,-50,50);
 881   fListHist->Add(fHistAsMcProdRadiusXvsYAntiLambda);
 882
 883   // Invariant mass distribution with PID checked:
 884   fHistPidMcMassK0             = new TH1F("h1PidMcMassK0", "K^{0} MC PId checked;M(#pi^{+}#pi^{-}) (GeV/c^{2});Counts", 100, 0.4, 0.6);
 885   fListHist->Add(fHistPidMcMassK0);
 886
 887   fHistPidMcMassLambda         = new TH1F("h1PidMcMassLambda", "#Lambda^{0} MC PId checked;M(p#pi^{-}) (GeV/c^{2});Counts", 75, 1.05, 1.2);
 888   fListHist->Add(fHistPidMcMassLambda);
 889
 890   fHistPidMcMassAntiLambda     = new TH1F("h1PidMcMassAntiLambda", "#bar{#Lambda}^{0} MC PId checked;M(#bar{p}#pi^{+}) (GeV/c^{2});Counts", 75, 1.05, 1.2);
 891   fListHist->Add(fHistPidMcMassAntiLambda);
 892
 893   // Invariant mass distribution:
 894   fHistAsMcMassK0              = new TH1F("h1AsMcMassK0", "K^{0} associated;M(#pi^{+}#pi^{-}) (GeV/c^{2});Counts", 100, 0.4, 0.6);
 895   fListHist->Add(fHistAsMcMassK0);
 896
 897   fHistAsMcMassLambda          = new TH1F("h1AsMcMassLambda", "#Lambda^{0} associated;M(p#pi^{-}) (GeV/c^{2});Counts", 75, 1.05, 1.2);
 898   fListHist->Add(fHistAsMcMassLambda);
 899
 900   fHistAsMcMassAntiLambda      = new TH1F("h1AsMcMassAntiLambda", "#bar{#Lambda}^{0} associated;M(#bar{p}#pi^{+}) (GeV/c^{2});Counts", 75, 1.05, 1.2);
 901   fListHist->Add(fHistAsMcMassAntiLambda);
 902
 903   // Pt vs. invariant mass:
 904   fHistAsMcPtVsMassK0               = new TH2F("h2AsMcPtVsMassK0","K^{0} associated;M(#pi^{+}#pi^{-}) (GeV/c^{2});p_{t} (GeV/c)",200, 0.4, 0.6,240,0,12);
 905   fListHist->Add(fHistAsMcPtVsMassK0);
 906
 907   fHistAsMcPtVsMassLambda           = new TH2F("h2AsMcPtVsMassLambda","#Lambda^{0} associated;M(p#pi^{-}) (GeV/c^{2});p_{t} (GeV/c)",140, 1.06, 1.2,240,0,12);
 908   fListHist->Add(fHistAsMcPtVsMassLambda);
 909
 910   fHistAsMcPtVsMassAntiLambda       = new TH2F("h2AsMcPtVsMassAntiLambda","#bar{#Lambda}^{0} associated;M(#bar{p}#pi^{+}) (GeV/c^{2});p_{t} (GeV/c)",140, 1.06, 1.2,240,0,12);
 911   fListHist->Add(fHistAsMcPtVsMassAntiLambda);
 912
 913
 914   // Invariant mass vs. radius:
 915   fHistAsMcMassVsRadiusK0             = new TH2F("h2AsMcMassVsRadiusK0", "K^{0} associated;radius (cm);M(#pi^{+}#pi^{-}) (GeV/c^{2})",lNbinRadius,radius, 500, 0.47, 0.52);
 916   fListHist->Add(fHistAsMcMassVsRadiusK0);
 917
 918   fHistAsMcMassVsRadiusLambda         = new TH2F("h2AsMcMassVsRadiusLambda", "#Lambda associated;radius (cm);M(p#pi^{-}) (GeV/c^{2})",lNbinRadius,radius, lNbinInvMassLambda, 1.10, 1.13);
 919   fListHist->Add(fHistAsMcMassVsRadiusLambda);
 920
 921   fHistAsMcMassVsRadiusAntiLambda     = new TH2F("h2AsMcMassVsRadiusAntiLambda", "#bar{#Lambda} associated;radius (cm);M(#bar{p}#pi^{+}) (GeV/c^{2})",lNbinRadius,radius,lNbinInvMassLambda , 1.10, 1.13);
 922   fListHist->Add(fHistAsMcMassVsRadiusAntiLambda);
 923
 924   // Position resolution for K0s:
 925   fHistAsMcResxK0                     = new TH1F("h1AsMcResxK0", "K^{0} associated;#Delta x (cm);Counts", 50, -0.25, 0.25);
 926   fListHist->Add(fHistAsMcResxK0);
 927   fHistAsMcResyK0                     = new TH1F("h1AsMcResyK0", "K^{0} associated;#Delta y (cm);Counts", 50, -0.25, 0.25);
 928   fListHist->Add(fHistAsMcResyK0);
 929   fHistAsMcReszK0                     = new TH1F("h1AsMcReszK0", "K^{0} associated;#Delta z (cm);Counts", 50, -0.25, 0.25);
 930   fListHist->Add(fHistAsMcReszK0);
 931
 932   // Position resolution vs. radius for K0s:
 933   fHistAsMcResrVsRadiusK0             = new TH2F("h2AsMcResrVsRadiusK0", "K^{0} associated;Radius (cm);#Delta r (cm)",200,0.0,50., 50, -0.25, 0.25);
 934   fListHist->Add(fHistAsMcResrVsRadiusK0);
 935   fHistAsMcReszVsRadiusK0             = new TH2F("h2AsMcReszVsRadiusK0", "K^{0} associated;Radius (cm);#Delta z (cm)",200,0.0,50.0, 50, -0.25, 0.25);
 936   fListHist->Add(fHistAsMcReszVsRadiusK0);
 937
 938   // Position resolution for Lambda:
 939   fHistAsMcResxLambda                 = new TH1F("h1AsMcResxLambda", "#Lambda^{0} associated;#Delta x (cm);Counts", 50, -0.25, 0.25);
 940   fListHist->Add(fHistAsMcResxLambda);
 941   fHistAsMcResyLambda                 = new TH1F("h1AsMcResyLambda", "#Lambda^{0} associated;#Delta y (cm);Counts", 50, -0.25, 0.25);
 942   fListHist->Add(fHistAsMcResyLambda);
 943   fHistAsMcReszLambda                 = new TH1F("h1AsMcReszLambda", "#Lambda^{0} associated;#Delta z (cm);Counts", 50, -0.25, 0.25);
 944   fListHist->Add(fHistAsMcReszLambda);
 945
 946   // Position resolution vs. radius for Lambda:
 947   fHistAsMcResrVsRadiusLambda         = new TH2F("h2AsMcResrVsRadiusLambda", "#Lambda^{0} associated;Radius (cm);#Delta r (cm)",200,0.0,50.0, 50, -0.25, 0.25);
 948   fListHist->Add(fHistAsMcResrVsRadiusLambda);
 949   fHistAsMcReszVsRadiusLambda         = new TH2F("h2AsMcReszVsRadiusLambda", "#Lambda^{0} associated;Radius (cm);#Delta z (cm)",200,0.0,50.0, 50, -0.25, 0.25);
 950   fListHist->Add(fHistAsMcReszVsRadiusLambda);
 951
 952   // Position resolution for anti-Lambda:
 953   fHistAsMcResxAntiLambda             = new TH1F("h1AsMcResxAntiLambda", "#bar{#Lambda}^{0} associated;#Delta x (cm);Counts", 50, -0.25, 0.25);
 954   fListHist->Add(fHistAsMcResxAntiLambda);
 955   fHistAsMcResyAntiLambda             = new TH1F("h1AsMcResyAntiLambda", "#bar{#Lambda}^{0} associated;#Delta y (cm);Counts", 50, -0.25, 0.25);
 956   fListHist->Add(fHistAsMcResyAntiLambda);
 957   fHistAsMcReszAntiLambda             = new TH1F("h1AsMcReszAntiLambda", "#bar{#Lambda}^{0} associated;#Delta z (cm);Counts", 50, -0.25, 0.25);
 958   fListHist->Add(fHistAsMcReszAntiLambda);
 959
 960   // Position resolution vs. radius for anti-Lambda:
 961   fHistAsMcResrVsRadiusAntiLambda     = new TH2F("h2AsMcResrVsRadiusAntiLambda", "#bar{#Lambda}^{0} associated;Radius (cm);#Delta r (cm)",200,0.0,50.0, 50, -0.25, 0.25);
 962   fListHist->Add(fHistAsMcResrVsRadiusAntiLambda);
 963   fHistAsMcReszVsRadiusAntiLambda     = new TH2F("h2AsMcReszVsRadiusAntiLambda", "#bar{#Lambda}^{0} associated;Radius (cm);#Delta z (cm)",200,0.0,50.0, 50, -0.25, 0.25);
 964   fListHist->Add(fHistAsMcReszVsRadiusAntiLambda);
 965
 966   // Pt Resolution:
 967   fHistAsMcResPtK0                   = new TH1F("h1AsMcResPtK0","Pt Resolution K^{0};#Delta Pt;Counts",200,-1,1);
 968   fListHist->Add(fHistAsMcResPtK0);
 969   fHistAsMcResPtLambda               = new TH1F("h1AsMcResPtLambda","Pt Resolution #Lambda^{0};#Delta Pt;Counts",200,-1,1);
 970   fListHist->Add(fHistAsMcResPtLambda);
 971   fHistAsMcResPtAntiLambda           = new TH1F("h1AsMcResPtAntiLambda","Pt Resolution #bar{#Lambda}^{0};#Delta Pt;Counts",200,-1,1);
 972   fListHist->Add(fHistAsMcResPtAntiLambda);
 973
 974   // Pt Resolution vs. rapidity:
 975   fHistAsMcResPtVsRapK0              = new TH2F("h2AsMcResPtVsRapK0","Pt Resolution K^{0};#Delta Pt;Rap",200,-1,1,20,-1,1);
 976   fListHist->Add(fHistAsMcResPtVsRapK0);
 977   fHistAsMcResPtVsRapLambda          = new TH2F("h2AsMcResPtVsRapLambda","Pt Resolution #Lambda^{0};#Delta Pt;Rap",200,-1,1,20,-1,1);
 978   fListHist->Add(fHistAsMcResPtVsRapLambda);
 979   fHistAsMcResPtVsRapAntiLambda      = new TH2F("h2AsMcResPtVsRapAntiLambda","Pt Resolution #bar{#Lambda}^{0};#Delta Pt;Rap",200,-1,1,20,-1,1);
 980   fListHist->Add(fHistAsMcResPtVsRapAntiLambda);
 981
 982   // Pt Resolution vs. Pt:
 983   fHistAsMcResPtVsPtK0               = new TH2F("h2AsMcResPtVsPtK0","Pt Resolution K^{0};#Delta Pt;Pt",600,-0.15,0.15,240,0,12);
 984   fListHist->Add(fHistAsMcResPtVsPtK0);
 985   fHistAsMcResPtVsPtLambda           = new TH2F("h2AsMcResPtVsPtLambda","Pt Resolution #Lambda^{0};#Delta Pt;Pt",600,-0.15,0.15,240,0,12);
 986   fListHist->Add(fHistAsMcResPtVsPtLambda);
 987   fHistAsMcResPtVsPtAntiLambda       = new TH2F("h2AsMcResPtVsPtAntiLambda","Pt Resolution #bar{#Lambda}^{0};#Delta Pt;Pt",300,-0.15,0.15,240,0,12);
 988   fListHist->Add(fHistAsMcResPtVsPtAntiLambda);
 989
 990   // Pdg code of mother particle:
 991   fHistAsMcMotherPdgCodeK0s           = new TH1F("h1AsMcMotherPdgCodeK0s","Mother of Associated K^{0};mother;counts",11,0,11);
 992   fListHist->Add(fHistAsMcMotherPdgCodeK0s);
 993   fHistAsMcMotherPdgCodeLambda        = new TH1F("h1AsMcMotherPdgCodeLambda","Mother of Associated #Lambda^{0};mother;counts",11,0,11);
 994   fListHist->Add(fHistAsMcMotherPdgCodeLambda);
 995   fHistAsMcMotherPdgCodeAntiLambda    = new TH1F("h1AsMcMotherPdgCodeAntiLambda","Mother of Associated #bar{#Lambda}^{0};mother;counts",11,0,11);
 996   fListHist->Add(fHistAsMcMotherPdgCodeAntiLambda);
 997
 998   // Pt distribution of Lambda <- Sigma decay
 999   fHistAsMcPtLambdaFromSigma          = new TH1F("h1AsMcPtLambdaFromSigma","#Lambda}^{0} associated from Sigma;p_{t} (GeV/c);Count",240,0,12);
1000   fListHist->Add(fHistAsMcPtLambdaFromSigma);
1001   fHistAsMcPtAntiLambdaFromSigma      = new TH1F("h1AsMcPtAntiLambdaFromSigma","#bar{#Lambda}^{0} associated from Sigma;p_{t} (GeV/c);Count",240,0,12);
1002   fListHist->Add(fHistAsMcPtAntiLambdaFromSigma);
1003
1004   //*******************************************
1005   // Associated secondary particles histograms
1006   //*******************************************
1007
1008   // Pt vs. rapidity distribution:
1009   fHistAsMcSecondaryPtVsRapK0s          = new TH2F("h2AsMcSecondaryPtVsRapK0s", "K^{0} associated secondary;p_{t} (GeV/c);rapidity",240,0,12,30,-1.5,1.5);
1010   fListHist->Add(fHistAsMcSecondaryPtVsRapK0s);
1011   fHistAsMcSecondaryPtVsRapLambda       = new TH2F("h2AsMcSecondaryPtVsRapLambda", "#Lambda^{0} associated secondary;p_{t} (GeV/c);rapidity",240,0,12,30,-1.5,1.5);
1012   fListHist->Add(fHistAsMcSecondaryPtVsRapLambda);
1013
1014   fHistAsMcSecondaryPtVsRapAntiLambda   = new TH2F("h2AsMcSecondaryPtVsRapAntiLambda", "#bar{#Lambda}^{0} associated secondary;p_{t} (GeV/c);rapidity",240,0,12,30,-1.5,1.5);
1015   fListHist->Add(fHistAsMcSecondaryPtVsRapAntiLambda);
1016
1017   // Production radius
1018   fHistAsMcSecondaryProdRadiusK0s              = new TH1F("h1AsMcSecondaryProdRadiusK0s", "K^{0} Production Radius;r (cm);Count", 170, -2, 15);
1019   fListHist->Add(fHistAsMcSecondaryProdRadiusK0s);
1020
1021   fHistAsMcSecondaryProdRadiusLambda           = new TH1F("h1AsMcSecondaryProdRadiusLambda", "#Lambda^{0} Production Radius;r (cm);Count", 170, -2, 15);
1022   fListHist->Add(fHistAsMcSecondaryProdRadiusLambda);
1023
1024   fHistAsMcSecondaryProdRadiusAntiLambda       = new TH1F("h1AsMcSecondaryProdRadiusAntiLambda", "#bar{#Lambda}^{0} Production Radius;r (cm);Count", 170, -2, 15);
1025   fListHist->Add(fHistAsMcSecondaryProdRadiusAntiLambda);
1026
1027   // Production radius vs. rapidity:
1028   fHistAsMcSecondaryProdRadiusXvsYK0s          = new TH2F("h2AsMcSecondaryProdRadiusXvsYK0s","Associated Secondary K^{0} Production Radius;x (cm); y (cm)",200,-20,20,200,-20,20);
1029   fListHist->Add(fHistAsMcSecondaryProdRadiusXvsYK0s);
1030   fHistAsMcSecondaryProdRadiusXvsYLambda       = new TH2F("h2AsMcSecondaryProdRadiusXvsYLambda","Associated Secondary #Lambda^{0} Production Radius;x (cm); y (cm)",200,-20,20,200,-20,20);
1031   fListHist->Add(fHistAsMcSecondaryProdRadiusXvsYLambda);
1032   fHistAsMcSecondaryProdRadiusXvsYAntiLambda   = new TH2F("h2AsMcSecondaryProdRadiusXvsYAntiLambda","Associated Secondary #bar{#Lambda}^{0} Production Radius;x (cm); y (cm)",200,-20,20,200,-20,20);
1033   fListHist->Add(fHistAsMcSecondaryProdRadiusXvsYAntiLambda);
1034
1035   // Pdg code of mother particle for secondary V0s:
1036   fHistAsMcSecondaryMotherPdgCodeK0s           = new TH1F("h1AsMcSecondaryMotherPdgCodeK0s","Mother of Associated Secondary K^{0};mother;counts",11,0,11);
1037   fListHist->Add(fHistAsMcSecondaryMotherPdgCodeK0s);
1038   fHistAsMcSecondaryMotherPdgCodeLambda        = new TH1F("h1AsMcSecondaryMotherPdgCodeLambda","Mother of Associated Secondary #Lambda^{0};mother;counts",11,0,11);
1039   fListHist->Add(fHistAsMcSecondaryMotherPdgCodeLambda);
1040   fHistAsMcSecondaryMotherPdgCodeAntiLambda    = new TH1F("h1AsMcSecondaryMotherPdgCodeAntiLambda","Mother of Associated Secondary #bar{#Lambda}^{0};mother;counts",11,0,11);
1041   fListHist->Add(fHistAsMcSecondaryMotherPdgCodeAntiLambda);
1042
1043   // Pt distribution of secondary Lambda <- Sigma decay:
1044   fHistAsMcSecondaryPtLambdaFromSigma          = new TH1F("h1AsMcSecondaryPtLambdaFromSigma","#Lambda}^{0} associated from Sigma;p_{t} (GeV/c);Count",240,0,12);
1045   fListHist->Add(fHistAsMcSecondaryPtLambdaFromSigma);
1046   fHistAsMcSecondaryPtAntiLambdaFromSigma      = new TH1F("h1AsMcSecondaryPtAntiLambdaFromSigma","#bar{#Lambda}^{0} associated from Sigma;p_{t} (GeV/c);Count",240,0,12);
1047   fListHist->Add(fHistAsMcSecondaryPtAntiLambdaFromSigma);
1048
1049   PostData(1, fListHist);
1050   PostData(2, fCentrSelector);
1051 }
1052
1053 //________________________________________________________________________
1054 void AliAnalysisTaskPerformanceStrange::UserExec(Option_t *)
1055 {
1056   // Main loop
1057   // Called for each event
1058
1059   AliStack* stack = NULL;
1060   TClonesArray *mcArray = NULL;
1061   TArrayF mcPrimaryVtx;
1062
1063   fESD=(AliESDEvent *)InputEvent();
1064
1065   if (!fESD) {
1066     Printf("ERROR: fESD not available");
1067     return;
1068   }
1069
1070   // FIXME: levent not used. Can I remove it?
1071   AliVEvent* lEvent = InputEvent();
1072
1073   if (!lEvent) {
1074     Printf("ERROR: Event not available");
1075     return;
1076   }
1077
1078   fHistNumberEvents->Fill(0.5);
1079
1080   //******************
1081   // Trigger Selection ! Warning Works only for ESD, add protection in case of AOD loop
1082   //******************
1083
1084   Bool_t isSelected =
1085     (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected()
1086      & AliVEvent::kMB);
1087   if (!isSelected) return;
1088
1089   // Centrality selection
1090   static AliESDtrackCuts * trackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(); // FIXME: make it a data member
1091   Bool_t isCentralitySelected = fCentrSelector->IsCentralityBinSelected(fESD,trackCuts);
1092   if(!isCentralitySelected) return;
1093   // FIXME: add to hist number events another entry for centrality.
1094
1095  // Done by the AliPhysicsSelection Task ! Only the selected events are passed to this task
1096
1097   fHistNumberEvents->Fill(1.5); // FIXME: use enum here
1098
1099
1100   //*************************
1101   //End track multiplicity
1102   //*************************
1103
1104
1105   // Remove Events with no tracks
1106   //if (!(fESD->GetNumberOfTracks()))  return;
1107
1108   fHistNumberEvents->Fill(2.5);
1109   fHistTrackPerEvent->Fill(fESD->GetNumberOfTracks());
1110
1111   //*************************************
1112   // Cut used:
1113   //*************************************
1114   // FIXME: Create a cut object, to be configured in the steering macro and to be streamed in the output to reference those cuts
1115   // Cut Rapidity:
1116   Double_t lCutRap  = 0.75;
1117
1118   // Cut AliKF Chi2 for Reconstructed particles
1119   Double_t cutChi2KF  = 1E3;
1120
1121   // If PID is used:
1122   Double_t lLimitPPID    = 0.7;
1123   Float_t cutNSigmaLowP  = 1E3;
1124   Float_t cutNSigmaHighP = 1E3;
1125   if (fUsePID.Contains("withPID")) {
1126     cutNSigmaLowP  = 5.0;
1127     cutNSigmaHighP = 3.0;
1128   }
1129
1130
1131   // Cut Daughters pt (GeV/c):
1132   Double_t cutMinPtDaughter = 0.160;
1133
1134   // Cut primary vertex:
1135   Double_t cutPrimVertex = 10.0;
1136
1137   // Min number of TPC clusters:
1138   Int_t nbMinTPCclusters = 80;
1139
1140   //*******************
1141   // PID parameters:
1142   //*******************
1143   // FIXME: OADB or momber TFormula?
1144   Double_t fAlephParameters[5] = {0,0,0,0,0,};
1145
1146   fAlephParameters[0] = 0.0283086;
1147   fAlephParameters[1] = 2.63394e+01;
1148   fAlephParameters[2] = 5.04114e-11;
1149   fAlephParameters[3] = 2.12543e+00;
1150   fAlephParameters[4] = 4.88663e+00;
1151
1152
1153   //*******************
1154   // Access MC:
1155   //*******************
1156   // FIXME:: move this two branches directly in the loops below
1157   if (fAnalysisMC) {
1158     if(fAnalysisType == "ESD") {
1159       AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
1160       if (!eventHandler) {
1161         Printf("ERROR: Could not retrieve MC event handler");
1162         return;
1163       }
1164       AliMCEvent* mcEvent = eventHandler->MCEvent();
1165       if (!mcEvent) {
1166         Printf("ERROR: Could not retrieve MC event");
1167         return;
1168       }
1169       stack = mcEvent->Stack();
1170       if (!stack) {
1171         Printf("ERROR: Could not retrieve stack");
1172         return;
1173       }
1174
1175       AliGenEventHeader* mcHeader=mcEvent->GenEventHeader();
1176       if(!mcHeader) return;
1177       mcHeader->PrimaryVertex(mcPrimaryVtx);
1178
1179     }
1180
1181     else if(fAnalysisType == "AOD") {
1182
1183       // load MC particles
1184       mcArray = (TClonesArray*)fESD->GetList()->FindObject(AliAODMCParticle::StdBranchName());
1185       if(!mcArray) {
1186         Printf("strange analysis::UserExec: MC particles branch not found!\n");
1187         return;
1188       }
1189
1190       // load MC header
1191       AliAODMCHeader *mcHeader =
1192         (AliAODMCHeader*)fESD->GetList()->FindObject(AliAODMCHeader::StdBranchName());
1193       if(!mcHeader) {
1194         Printf("strange analysis::UserExec: MC header branch not found!\n");
1195         return;
1196       }
1197     }
1198
1199     // PID parameters for MC simulations:
1200     // FIXME: set above, with the others
1201     fAlephParameters[0] = 2.15898e+00/50.;
1202     fAlephParameters[1] = 1.75295e+01;
1203     fAlephParameters[2] = 3.40030e-09;
1204     fAlephParameters[3] = 1.96178e+00;
1205     fAlephParameters[4] = 3.91720e+00;
1206   }
1207
1208
1209   //**********************************************
1210   // MC loop
1211   //**********************************************
1212
1213   Double_t lmcPrimVtxR      = 0;
1214
1215   Int_t lNbMCPrimary        = 0;
1216   Int_t lNbMCPart           = 0;
1217
1218   Int_t lPdgcodeCurrentPart = 0;
1219   Double_t lRapCurrentPart  = 0;
1220   Double_t lPtCurrentPart   = 0;
1221
1222   Int_t lComeFromSigma      = 0;
1223
1224
1225   // Production Radius
1226   Double_t mcPosX     = 0.0,  mcPosY      = 0.0,  mcPosZ      = 0.0;
1227   Double_t mcPosR     = 0.0;
1228
1229   // Decay Radius
1230   Double_t mcDecayPosX = 0, mcDecayPosY = 0, mcDecayPosR = 0;
1231
1232   // current mc particle 's mother
1233   Int_t iCurrentMother  = 0, lPdgCurrentMother    = 0;
1234   Bool_t lCurrentMotherIsPrimary;
1235
1236   // current mc particles 's daughter:
1237   Int_t lPdgCurrentDaughter0 = 0, lPdgCurrentDaughter1 = 0;
1238
1239   // variables for multiple reconstruction studies:
1240   Int_t id0           = 0, id1          = 0;
1241   //Int_t lLabelTrackN  = 0, lLabelTrackP = 0;
1242   //Int_t lPartNMother  = 0, lPartPMother = 0;
1243   //Int_t lPartPMotherPDGcode      = 0;
1244   Int_t lNtimesReconstructedK0s   = 0, lNtimesReconstructedLambda   = 0, lNtimesReconstructedAntiLambda   = 0;
1245  // Int_t lNtimesReconstructedK0sMI = 0, lNtimesReconstructedLambdaMI = 0, lNtimesReconstructedAntiLambdaMI = 0;
1246
1247   //****************************
1248   // Start loop over MC particles
1249   if (fAnalysisMC) {
1250
1251     // Primary vertex
1252     fHistMCPrimaryVertexX->Fill(mcPrimaryVtx.At(0));
1253     fHistMCPrimaryVertexY->Fill(mcPrimaryVtx.At(1));
1254     fHistMCPrimaryVertexZ->Fill(mcPrimaryVtx.At(2));
1255
1256     lmcPrimVtxR = TMath::Sqrt(mcPrimaryVtx.At(0)*mcPrimaryVtx.At(0)+mcPrimaryVtx.At(1)*mcPrimaryVtx.At(1));
1257
1258     // FIXME: move these loops to other functions, for better readibility?
1259     if(fAnalysisType == "ESD") {
1260
1261       lNbMCPrimary = stack->GetNprimary(); // FIXME: This does not correspond to our definition of primaries, but maybe it is ok for strange particles
1262       lNbMCPart    = stack->GetNtrack();
1263
1264       fHistMCMultiplicityPrimary->Fill(lNbMCPrimary);
1265       fHistMCMultiplicityTracks->Fill(lNbMCPart);
1266
1267
1268       for (Int_t iMc = 0; iMc < (stack->GetNtrack()); iMc++) {
1269         TParticle *p0 = stack->Particle(iMc);
1270         if (!p0) {
1271           //Printf("ERROR: particle with label %d not found in stack (mc loop)", iMc);
1272           continue;
1273         }
1274         lPdgcodeCurrentPart = p0->GetPdgCode();
1275
1276         // Keep only K0s, Lambda and AntiLambda, Xi and Phi:
1277         if ( (lPdgcodeCurrentPart != 310 ) && (lPdgcodeCurrentPart != 3122 ) && (lPdgcodeCurrentPart != -3122 ) && (lPdgcodeCurrentPart != 3312 ) && (lPdgcodeCurrentPart != -3312) && (lPdgcodeCurrentPart != -333) ) continue;
1278
1279         lRapCurrentPart   = MyRapidity(p0->Energy(),p0->Pz());
1280         //lEtaCurrentPart   = p0->Eta();
1281         lPtCurrentPart    = p0->Pt();
1282
1283         iCurrentMother    = p0->GetFirstMother();
1284
1285 //      lPdgCurrentMother = stack->Particle(iCurrentMother)->GetPdgCode();
1286         if (iCurrentMother == -1){lPdgCurrentMother=0; } else {lPdgCurrentMother = stack->Particle(iCurrentMother)->GetPdgCode();}
1287
1288         mcPosX = p0->Vx();
1289         mcPosY = p0->Vy();
1290         mcPosZ = p0->Vz();
1291         mcPosR = TMath::Sqrt(mcPosX*mcPosX+mcPosY*mcPosY);
1292
1293         id0  = p0->GetDaughter(0);
1294         id1  = p0->GetDaughter(1);
1295
1296         // Decay Radius
1297         if ( id0 <= lNbMCPart && id0 > 0 && id1 <= lNbMCPart && id1 > 0) {
1298           TParticle *pDaughter0 = stack->Particle(id0);
1299           TParticle *pDaughter1 = stack->Particle(id1);
1300           lPdgCurrentDaughter0 = pDaughter0->GetPdgCode();
1301           lPdgCurrentDaughter1 = pDaughter1->GetPdgCode();
1302
1303           mcDecayPosX = pDaughter0->Vx();
1304           mcDecayPosY = pDaughter0->Vy();
1305           mcDecayPosR = TMath::Sqrt(mcDecayPosX*mcDecayPosX+mcDecayPosY*mcDecayPosY);
1306         }
1307         else  {
1308           //FIXME: shouldn't this be a fatal?
1309           //Printf("ERROR: particle with label %d and/or %d not found in stack (mc loop)", id0,id1);
1310           mcDecayPosR = -1.0;
1311         }
1312         // FIXME using array of histos and conversion PDGCode -> enum would make this much easier to read
1313         // We could also have a function FillMcHistos (pos, radius, rap...) which we call from both the AOD and ESD loops
1314         if (lPdgcodeCurrentPart==310)   {
1315           fHistMCtracksProdRadiusK0s->Fill(mcPosX,mcPosY);
1316           fHistMCtracksDecayRadiusK0s->Fill(mcDecayPosR);
1317           if (TMath::Abs(lRapCurrentPart) < lCutRap) fHistMCPtAllK0s->Fill(lPtCurrentPart);
1318         }
1319         else if (lPdgcodeCurrentPart==3122)  {
1320           fHistMCtracksProdRadiusLambda->Fill(mcPosX,mcPosY);
1321           fHistMCtracksDecayRadiusLambda->Fill(mcDecayPosR);
1322           if (TMath::Abs(lRapCurrentPart) < lCutRap) fHistMCPtAllLambda->Fill(lPtCurrentPart);
1323         }
1324         else if (lPdgcodeCurrentPart==-3122) {
1325           fHistMCtracksProdRadiusAntiLambda->Fill(mcPosX,mcPosY);
1326           fHistMCtracksDecayRadiusAntiLambda->Fill(mcDecayPosR);
1327           if (TMath::Abs(lRapCurrentPart) < lCutRap) fHistMCPtAllAntiLambda->Fill(lPtCurrentPart);
1328         }
1329
1330           // FIXME: not sure if I understand this: is it correct? (definition of primaries)
1331         if ( ( ( TMath::Abs(lPdgCurrentMother) == 3212)  ||
1332                ( TMath::Abs(lPdgCurrentMother) == 3224)  ||
1333                ( TMath::Abs(lPdgCurrentMother) == 3214)  ||
1334                ( TMath::Abs(lPdgCurrentMother) == 3114) )
1335              && ( iCurrentMother <= lNbMCPrimary )
1336              ) lComeFromSigma = 1;
1337         else lComeFromSigma = 0;
1338
1339         //*********************************************
1340         // Now keep only primary particles
1341         if ( ( iMc > lNbMCPrimary ) && (!lComeFromSigma) ) continue;
1342
1343         //********************************************
1344         //check if V0 is reconstructed several times
1345
1346         lNtimesReconstructedK0s   = 0; lNtimesReconstructedLambda   = 0; lNtimesReconstructedAntiLambda   = 0;
1347
1348         //for (Int_t jV0 = 0; jV0 < fESD->GetNumberOfV0s(); jV0++) {
1349
1350         //lLabelTrackN  = 0; lLabelTrackP = 0;
1351         //lPartNMother  = 0; lPartPMother = 0;
1352
1353         //AliESDv0    *vertexESD = ((AliESDEvent*)fESD)->GetV0(jV0);
1354         //if (!vertexESD) continue;
1355
1356         //AliESDtrack *trackNESD = ((AliESDEvent*)fESD)->GetTrack(TMath::Abs(vertexESD->GetNindex()));
1357         //lLabelTrackN = (UInt_t)TMath::Abs(trackNESD->GetLabel());
1358         //if (lLabelTrackN!=id0 && lLabelTrackN!=id1) continue;
1359
1360         //AliESDtrack *trackPESD = ((AliESDEvent*)fESD)->GetTrack(TMath::Abs(vertexESD->GetPindex()));
1361         //lLabelTrackP = (UInt_t)TMath::Abs(trackPESD->GetLabel());
1362         //if (lLabelTrackP!=id0 && lLabelTrackP!=id1) continue;
1363
1364         //TParticle   *lPartNESD = stack->Particle(lLabelTrackN);
1365         //TParticle   *lPartPESD = stack->Particle(lLabelTrackP);
1366         //lPartNMother = lPartNESD->GetFirstMother();
1367         //lPartPMother = lPartPESD->GetFirstMother();
1368
1369         //lPartPMotherPDGcode = stack->Particle(lPartPMother)->GetPdgCode();
1370
1371         //switch (vertexESD->GetOnFlyStatus()){
1372
1373         //case 0 :
1374         //if ( (lPartPMother==lPartNMother) && (lPartPMotherPDGcode==310) ) lNtimesReconstructedK0s++;
1375         //else if ( (lPartPMother==lPartNMother) && (lPartPMotherPDGcode==3122) ) lNtimesReconstructedLambda++;
1376         //else if ( (lPartPMother==lPartNMother) && (lPartPMotherPDGcode==-3122) ) lNtimesReconstructedAntiLambda++;
1377         //break;
1378
1379         //case 1 :
1380         //if ( (lPartPMother==lPartNMother) && (lPartPMotherPDGcode==310) ) lNtimesReconstructedK0sMI++;
1381         //else if ( (lPartPMother==lPartNMother) && (lPartPMotherPDGcode==3122) ) lNtimesReconstructedLambdaMI++;
1382         //else if ( (lPartPMother==lPartNMother) && (lPartPMotherPDGcode==-3122) ) lNtimesReconstructedAntiLambdaMI++;
1383         //break;
1384
1385         //}
1386         //} // end loop over reconstructed V0s inside MC loop
1387
1388         // FIXME: same comemtn for array of histos
1389         // Rap distribution
1390         if (lPdgcodeCurrentPart==310) {
1391           fHistMCRapK0s->Fill(lRapCurrentPart);
1392           if (lPtCurrentPart < 0.2 && lPtCurrentPart < 3.0)
1393             fHistMCRapInPtRangeK0s->Fill(lRapCurrentPart);
1394         }
1395
1396         if (lPdgcodeCurrentPart==3122) {
1397           fHistMCRapLambda->Fill(lRapCurrentPart);
1398           if (lPtCurrentPart < 0.6 && lPtCurrentPart < 3.5)
1399             fHistMCRapInPtRangeLambda->Fill(lRapCurrentPart);
1400         }
1401
1402         if (lPdgcodeCurrentPart==-3122) {
1403           fHistMCRapAntiLambda->Fill(lRapCurrentPart);
1404           if (lPtCurrentPart < 0.6 && lPtCurrentPart < 3.5)
1405             fHistMCRapInPtRangeAntiLambda->Fill(lRapCurrentPart);
1406         }
1407
1408         if (lPdgcodeCurrentPart==3312 || lPdgcodeCurrentPart==-3312) {
1409           fHistMCRapXi->Fill(lRapCurrentPart);
1410           if (lPtCurrentPart < 0.6 && lPtCurrentPart < 3.0)
1411             fHistMCRapInPtRangeXi->Fill(lRapCurrentPart);
1412         }
1413
1414         if (lPdgcodeCurrentPart==333) {
1415           fHistMCRapPhi->Fill(lRapCurrentPart);
1416           if (lPtCurrentPart < 0.7 && lPtCurrentPart < 3.0)
1417             fHistMCRapInPtRangePhi->Fill(lRapCurrentPart);
1418         }
1419
1420         // Rapidity Cut
1421         if (TMath::Abs(lRapCurrentPart) > lCutRap) continue;
1422
1423         if (lPdgcodeCurrentPart==310) {
1424           fHistMCProdRadiusK0s->Fill(mcPosR);
1425
1426           fHistMCPtK0s->Fill(lPtCurrentPart);
1427
1428           fHistNTimesRecK0s->Fill(lNtimesReconstructedK0s);
1429           fHistNTimesRecK0sVsPt->Fill(lPtCurrentPart,lNtimesReconstructedK0s);
1430         }
1431         else
1432         if (lPdgcodeCurrentPart==3122) {
1433           fHistMCProdRadiusLambda->Fill(mcPosR);
1434
1435           fHistMCPtLambda->Fill(lPtCurrentPart);
1436
1437
1438           fHistNTimesRecLambda->Fill(lNtimesReconstructedLambda);
1439           fHistNTimesRecLambdaVsPt->Fill(lPtCurrentPart,lNtimesReconstructedLambda);
1440           if (lComeFromSigma) fHistMCPtLambdaFromSigma->Fill(lPtCurrentPart);
1441           //printf("found Lambda MC pT=%e\n",lPtCurrentPart);
1442           //printf("found Lambda MC Plabel=%d PPDGcode=%d Nlabel=%d NPDGcode=%d\n\n",id0,lPdgCurrentDaughter0,id1,lPdgCurrentDaughter1);
1443
1444         }
1445
1446       } // end loop ESD MC
1447
1448     } // end ESD condition
1449
1450     // FIXME: I skipped the AOD loop
1451     else if(fAnalysisType == "AOD") {
1452       lNbMCPart = mcArray->GetEntriesFast();
1453       lNbMCPrimary = 0;
1454
1455       fHistMCMultiplicityTracks->Fill(lNbMCPart);
1456
1457       for (Int_t iMc = 0; iMc < lNbMCPart; iMc++) {
1458
1459         // Primary vertex TO DO !!
1460         //
1461
1462         AliAODMCParticle *mcAODPart = (AliAODMCParticle*)mcArray->At(iMc);
1463         if (!mcAODPart) {
1464           //Printf("Strange analysis task (mc loop): particle with label %d not found", iMc);
1465           continue;
1466         }
1467         lPdgcodeCurrentPart = mcAODPart->GetPdgCode();
1468         if (mcAODPart->IsPhysicalPrimary()) {lNbMCPrimary = lNbMCPrimary +1;}
1469
1470         // Keep only K0s, Lambda and AntiLambda:
1471         if ( (lPdgcodeCurrentPart != 310 ) && (lPdgcodeCurrentPart != 3122 ) && (lPdgcodeCurrentPart != -3122 ) ) continue;
1472
1473         //lEtaCurrentPart   = mcAODPart->Eta();
1474         lRapCurrentPart   = mcAODPart->Y();
1475         lPtCurrentPart    = mcAODPart->Pt();
1476         iCurrentMother    = mcAODPart->GetMother();
1477         lPdgCurrentMother = ((AliAODMCParticle*)mcArray->At(iCurrentMother))->GetPdgCode();
1478         lCurrentMotherIsPrimary = ((AliAODMCParticle*)mcArray->At(iCurrentMother))->IsPhysicalPrimary();
1479
1480         mcPosX = mcAODPart->Xv();
1481         mcPosY = mcAODPart->Yv();
1482         mcPosZ = mcAODPart->Zv();
1483         mcPosR = TMath::Sqrt(mcPosX*mcPosX+mcPosY*mcPosY);
1484
1485         id0  = mcAODPart->GetDaughter(0);
1486         id1  = mcAODPart->GetDaughter(1);
1487
1488         // Decay Radius and Production Radius
1489         if ( id0 <= lNbMCPart && id0 > 0 && id1 <= lNbMCPart && id1 > 0) {
1490           AliAODMCParticle *mcAODDaughter1 = (AliAODMCParticle*)mcArray->At(id1);
1491           if (!mcAODPart) {
1492             //Printf("Strange analysis task (mc loop): daughter not found");
1493             continue;
1494           }
1495           mcDecayPosX = mcAODDaughter1->Xv();
1496           mcDecayPosY = mcAODDaughter1->Yv();
1497           mcDecayPosR = TMath::Sqrt(mcDecayPosX*mcDecayPosX+mcDecayPosY*mcDecayPosY);
1498         }
1499         else  {
1500           //Printf("ERROR: particle with label %d and/or %d not found in stack (mc loop)", id0,id1);
1501           mcDecayPosR = -1.0;
1502         }
1503
1504         if (lPdgcodeCurrentPart==310)   {
1505           fHistMCtracksProdRadiusK0s->Fill(mcPosX,mcPosY);
1506           fHistMCtracksDecayRadiusK0s->Fill(mcDecayPosR);
1507           if (TMath::Abs(lRapCurrentPart) < lCutRap) fHistMCPtAllK0s->Fill(lPtCurrentPart);
1508         }
1509         else if (lPdgcodeCurrentPart==3122)  {
1510           fHistMCtracksProdRadiusLambda->Fill(mcPosX,mcPosY);
1511           fHistMCtracksDecayRadiusLambda->Fill(mcDecayPosR);
1512           if (TMath::Abs(lRapCurrentPart) < lCutRap) fHistMCPtAllLambda->Fill(lPtCurrentPart);
1513         }
1514         else if (lPdgcodeCurrentPart==-3122) {
1515           fHistMCtracksProdRadiusAntiLambda->Fill(mcPosX,mcPosY);
1516           fHistMCtracksDecayRadiusAntiLambda->Fill(mcDecayPosR);
1517           if (TMath::Abs(lRapCurrentPart) < lCutRap) fHistMCPtAllAntiLambda->Fill(lPtCurrentPart);
1518         }
1519
1520         if ( ( ( TMath::Abs(lPdgCurrentMother) == 3212)  ||
1521                ( TMath::Abs(lPdgCurrentMother) == 3224)  ||
1522                ( TMath::Abs(lPdgCurrentMother) == 3214)  ||
1523                ( TMath::Abs(lPdgCurrentMother) == 3114) )
1524              && (lCurrentMotherIsPrimary)
1525              ) lComeFromSigma = 1;
1526         else lComeFromSigma = 0;
1527
1528         //*********************************************
1529         // Now keep only primary particles
1530
1531         // FIX IT !!!!    iMC is not defined !!!! FIX IT also in ESD/AOD loop !!
1532         if ( ( iMc > lNbMCPrimary ) && (!lComeFromSigma) ) continue;
1533
1534         //********************************************
1535         // check if V0 is reconstructed several times
1536
1537         //lNtimesReconstructedK0s   = 0; lNtimesReconstructedLambda   = 0; lNtimesReconstructedAntiLambda   = 0;
1538         //lNtimesReconstructedK0sMI = 0; lNtimesReconstructedLambdaMI = 0; lNtimesReconstructedAntiLambdaMI = 0;
1539
1540         //for (Int_t jV0 = 0; jV0 < fESD->GetNumberOfV0s(); jV0++) {
1541
1542         //lLabelTrackN  = 0; lLabelTrackP = 0;
1543         //lPartNMother  = 0; lPartPMother = 0;
1544
1545         //AliAODv0    *vertexAOD= ((AliAODEvent*)fESD)->GetV0(jV0);
1546         //if (!vertexAOD) continue;
1547         //printf("enter!!");
1548         //AliVParticle  *trackP  = ((AliVEvent*)fESD)->GetTrack(vertexAOD->GetPosID());
1549         //if (!trackP) continue;
1550         //lLabelTrackP = TMath::Abs(trackP->GetLabel());
1551         //if (lLabelTrackP!=id0 && lLabelTrackP!=id1) continue;
1552
1553         //AliVParticle  *trackN  = ((AliVEvent*)fESD)->GetTrack(vertexAOD->GetNegID());
1554         //if (!trackN) continue;
1555         //lLabelTrackN = TMath::Abs(trackN->GetLabel());
1556         //if (lLabelTrackN!=id0 && lLabelTrackN!=id1) continue;
1557
1558         //AliAODMCParticle *lPartNAOD = (AliAODMCParticle*)mcArray->At(lLabelTrackN);
1559         //if (!lPartNAOD) continue;
1560         //AliAODMCParticle *lPartPAOD = (AliAODMCParticle*)mcArray->At(lLabelTrackP);
1561         //if (!lPartPAOD) continue;
1562
1563         //lPartNMother = lPartNAOD->GetMother();
1564         //lPartPMother = lPartPAOD->GetMother();
1565
1566         //lPartPMotherPDGcode = ((AliAODMCParticle*)mcArray->At(lPartPMother))->GetPdgCode();
1567
1568         //switch (vertexAOD->GetOnFlyStatus()){
1569
1570         //case 0 :
1571           //if ( (lPartPMother==lPartNMother) && (lPartPMotherPDGcode==310) ) lNtimesReconstructedK0s++;
1572           //else if ( (lPartPMother==lPartNMother) && (lPartPMotherPDGcode==3122) ) lNtimesReconstructedLambda++;
1573           //else if ( (lPartPMother==lPartNMother) && (lPartPMotherPDGcode==-3122) ) lNtimesReconstructedAntiLambda++;
1574           //break;
1575
1576         //case 1 :
1577           //if ( (lPartPMother==lPartNMother) && (lPartPMotherPDGcode==310) ) lNtimesReconstructedK0sMI++;
1578           //else if ( (lPartPMother==lPartNMother) && (lPartPMotherPDGcode==3122) ) lNtimesReconstructedLambdaMI++;
1579           //else if ( (lPartPMother==lPartNMother) && (lPartPMotherPDGcode==-3122) ) lNtimesReconstructedAntiLambdaMI++;
1580           //break;
1581
1582        ///}
1583       //} // end loop over reconstructed V0s inside MC loop
1584
1585         if (TMath::Abs(lRapCurrentPart) > lCutRap) continue;
1586
1587         if (lPdgcodeCurrentPart==310) {
1588           fHistMCProdRadiusK0s->Fill(mcPosR);
1589           fHistMCPtK0s->Fill(lPtCurrentPart);
1590           fHistNTimesRecK0s->Fill(lNtimesReconstructedK0s);
1591           fHistNTimesRecK0sVsPt->Fill(lPtCurrentPart,lNtimesReconstructedK0s);
1592         }
1593         else if (lPdgcodeCurrentPart==3122) {
1594           fHistMCProdRadiusLambda->Fill(mcPosR);
1595           fHistMCPtLambda->Fill(lPtCurrentPart);
1596           fHistNTimesRecLambda->Fill(lNtimesReconstructedLambda);
1597           fHistNTimesRecLambdaVsPt->Fill(lPtCurrentPart,lNtimesReconstructedLambda);
1598           if (lComeFromSigma) fHistMCPtLambdaFromSigma->Fill(lPtCurrentPart);
1599         }
1600         else if (lPdgcodeCurrentPart==-3122) {
1601           fHistMCProdRadiusAntiLambda->Fill(mcPosR);
1602           fHistNTimesRecAntiLambda->Fill(lNtimesReconstructedAntiLambda);
1603           fHistNTimesRecAntiLambdaVsPt->Fill(lPtCurrentPart,lNtimesReconstructedAntiLambda);
1604           if (lComeFromSigma) fHistMCPtAntiLambdaFromSigma->Fill(lPtCurrentPart);
1605         }
1606
1607       } // end loop over AODMC particles
1608       fHistMCMultiplicityPrimary->Fill(lNbMCPrimary);
1609
1610     } // end AOD condition
1611
1612   } // End Loop over MC condition
1613
1614
1615
1616
1617   //************************************
1618   // ESD or AOD loop
1619   //************************************
1620
1621   Double_t lMagneticField = 999;
1622
1623   //Multiplcity:
1624   Int_t    nv0sTot= 0, nv0s = 0;
1625 //  Int_t nv0sMI =0;
1626   // Variables:
1627   Double_t  lV0Position[3];
1628
1629   Double_t lDcaPosToPrimVertex = 0;
1630   Double_t lDcaNegToPrimVertex = 0;
1631   Double_t lDcaV0Daughters     = 0;
1632   Double_t lV0cosPointAngle    = 0;
1633   Double_t lChi2V0             = 0;
1634   Double_t lV0DecayLength      = 0;
1635   Double_t lV0Radius           = 0;
1636   Double_t lDcaV0ToPrimVertex  = 0;
1637
1638   Int_t    lOnFlyStatus        = 0;
1639   //Float_t   tdcaPosToPrimVertexXYZ[2], tdcaNegToPrimVertexXYZ[2]; // ..[0] = Impact parameter in XY plane and ..[1] = Impact parameter in Z
1640   //Double_t  tdcaDaughterToPrimVertex[2];                          // ..[0] = Pos and ..[1] = Neg
1641
1642
1643
1644   Double_t lInvMassK0s = 0, lInvMassLambda = 0, lInvMassAntiLambda = 0;
1645   Double_t lPtK0s      = 0, lPtLambda      = 0, lPtAntiLambda      = 0;
1646   Double_t lRapK0s     = 0, lRapLambda     = 0, lRapAntiLambda     = 0;
1647   Double_t lEtaK0s     = 0, lEtaLambda     = 0, lEtaAntiLambda     = 0;
1648   Double_t lAlphaV0      = 0, lPtArmV0       = 0;
1649
1650   Double_t lPzK0s      = 0, lPzLambda      = 0, lPzAntiLambda      = 0;
1651
1652
1653   Double_t lV0Eta = 999;
1654
1655   // to study Associated V0s:
1656   Int_t    lIndexTrackPos       = 0, lIndexTrackNeg         = 0;
1657   UInt_t   lLabelTrackPos       = 0, lLabelTrackNeg         = 0;
1658   Int_t    lCheckPIdK0Short     = 0, lCheckMcK0Short        = 0;
1659   Int_t    lCheckPIdLambda      = 0, lCheckMcLambda         = 0;
1660   Int_t    lCheckPIdAntiLambda  = 0, lCheckMcAntiLambda     = 0;
1661   Int_t    lCheckSecondaryK0s   = 0, lCheckSecondaryLambda  = 0, lCheckSecondaryAntiLambda  = 0;
1662   Int_t    lCheckGamma          = 0;
1663   Double_t mcPosMotherX         = 0, mcPosMotherY           = 0, mcPosMotherZ  = 0;
1664   Double_t mcPosMotherR         = 0;
1665   Double_t mcMotherPt           = 0;
1666
1667   Int_t lIndexPosMother        = 0;
1668   Int_t lIndexNegMother        = 0;
1669   Int_t lIndexMotherOfMother   = 0;
1670   Int_t lPDGCodePosDaughter    = 0;
1671   Int_t lPDGCodeNegDaughter    = 0;
1672   Int_t lPdgcodeMother         = 0;
1673   Int_t lPdgcodeMotherOfMother = 0;
1674
1675   // Reconstructed position
1676   Double_t rcPosXK0s        = 0,  rcPosYK0s        = 0, rcPosZK0s        = 0;
1677   Double_t rcPosRK0s        = 0;
1678   Double_t rcPosXLambda     = 0,  rcPosYLambda     = 0, rcPosZLambda     = 0;
1679   Double_t rcPosRLambda     = 0;
1680   Double_t rcPosXAntiLambda = 0,  rcPosYAntiLambda = 0, rcPosZAntiLambda = 0;
1681   Double_t rcPosRAntiLambda = 0;
1682
1683   // Pt resolution
1684   Double_t deltaPtK0s  = 0, deltaPtLambda  = 0, deltaPtAntiLambda  = 0;
1685
1686   // Daughters
1687   AliESDtrack  *myTrackPos  = NULL;
1688   AliESDtrack  *myTrackNeg  = NULL;
1689   AliVParticle *lVPartPos   = NULL;
1690   AliVParticle *lVPartNeg   = NULL;
1691
1692   // Daughters' momentum:
1693   Double_t  lMomPos[3] = {999,999,999};
1694   Double_t  lMomNeg[3] = {999,999,999};
1695   Double_t  lPtPos = 999, lPtNeg = 999;
1696   Double_t  lPPos = 999, lPNeg = 999;
1697
1698   // Inner Wall parameters:
1699   Double_t  lMomInnerWallPos =999, lMomInnerWallNeg = 999;
1700
1701   // AliKF Chi2 and Armenteros variables
1702   Double_t lChi2KFK0s  = 0, lChi2KFLambda = 0,  lChi2KFAntiLambda = 0;
1703   Double_t lAlphaV0K0s = 0, lAlphaV0Lambda = 0,  lAlphaV0AntiLambda = 0;
1704   Double_t lPtArmV0K0s = 0, lPtArmV0Lambda = 0,  lPtArmV0AntiLambda = 0;
1705   Double_t lQlPos   = 0, lQlNeg   = 0;
1706
1707
1708   // PID
1709   Float_t nSigmaPosPion   = 0;
1710   Float_t nSigmaNegPion   = 0;
1711
1712   Float_t nSigmaPosProton = 0;
1713   Float_t nSigmaNegProton = 0;
1714
1715   Int_t lCheckPIDK0sPosDaughter        = 0, lCheckPIDK0sNegDaughter        = 0;
1716   Int_t lCheckPIDLambdaPosDaughter     = 0, lCheckPIDLambdaNegDaughter     = 0;
1717   Int_t lCheckPIDAntiLambdaPosDaughter = 0, lCheckPIDAntiLambdaNegDaughter = 0;
1718
1719
1720
1721   //***********************
1722   // Primary Vertex cuts &
1723   // Magnetic field and Quality tracks cuts
1724
1725   Double_t  lPrimaryVtxPosition[3];
1726   Double_t  lPrimaryVtxCov[6];
1727   Double_t  lPrimaryVtxChi2 = 999;
1728   Double_t  lResPrimaryVtxX = 999;
1729   Double_t  lResPrimaryVtxY = 999;
1730   Double_t  lResPrimaryVtxZ = 999;
1731
1732   AliAODVertex *myPrimaryVertex = NULL;
1733   //const AliVVertex *mySPDPrimaryVertex = NULL;
1734
1735   AliESDtrackCuts *myTracksCuts = NULL;
1736
1737   const AliMultiplicity *myMultiplicty = ((AliESDEvent*)fESD)->GetMultiplicity();
1738
1739   if(fAnalysisType == "ESD") {
1740
1741     // Best Primary Vertex:
1742     const AliESDVertex *myBestPrimaryVertex = ((AliESDEvent*)fESD)->GetPrimaryVertex();
1743     myBestPrimaryVertex = ((AliESDEvent*)fESD)->GetPrimaryVertex();
1744     if (!myBestPrimaryVertex) return;
1745     if (!myBestPrimaryVertex->GetStatus()) return;
1746     fHistNumberEvents->Fill(3.5);
1747     myBestPrimaryVertex->GetXYZ(lPrimaryVtxPosition);
1748     myBestPrimaryVertex->GetCovMatrix(lPrimaryVtxCov);
1749     if ( ( TMath::Abs(lPrimaryVtxPosition[2]) ) > cutPrimVertex) return ;
1750     fHistNumberEvents->Fill(4.5);
1751     lPrimaryVtxChi2 = myBestPrimaryVertex->GetChi2toNDF();
1752     lResPrimaryVtxX = myBestPrimaryVertex->GetXRes();
1753     lResPrimaryVtxY = myBestPrimaryVertex->GetYRes();
1754     lResPrimaryVtxZ = myBestPrimaryVertex->GetZRes();
1755
1756     // remove TPC-only primary vertex : retain only events with tracking + SPD vertex
1757     const AliESDVertex *mySPDPrimaryVertex = ((AliESDEvent*)fESD)->GetPrimaryVertexSPD();
1758     if (!mySPDPrimaryVertex) return;
1759     fHistSPDPrimaryVertexZ->Fill(mySPDPrimaryVertex->GetZ());
1760     const AliESDVertex *myPrimaryVertexTracking = ((AliESDEvent*)fESD)->GetPrimaryVertexTracks();
1761     if (!myPrimaryVertexTracking) return;
1762     if (!mySPDPrimaryVertex->GetStatus() && !myPrimaryVertexTracking->GetStatus() ) return;
1763     fHistNumberEvents->Fill(5.5);
1764
1765
1766     myPrimaryVertex = new AliAODVertex(lPrimaryVtxPosition, lPrimaryVtxCov, lPrimaryVtxChi2, NULL, -1, AliAODVertex::kPrimary);
1767     if (!myPrimaryVertex) return;
1768
1769
1770      // Number of Tracklets:
1771     //const AliMultiplicity *myMultiplicty = ((AliESDEvent*)fESD)->GetMultiplicity();
1772     //if (myMultiplicty->GetNumberOfTracklets() < 10) return;
1773     fHistTrackletPerEvent->Fill(myMultiplicty->GetNumberOfTracklets());
1774
1775     lMagneticField = ((AliESDEvent*)fESD)->GetMagneticField();
1776
1777     myTracksCuts = new AliESDtrackCuts();
1778     // require TPC refit
1779     myTracksCuts->SetRequireTPCRefit(kTRUE);
1780     // minimum number of clusters in TPC
1781     myTracksCuts->SetMinNClustersTPC(nbMinTPCclusters);
1782
1783   }
1784
1785   else if(fAnalysisType == "AOD") {
1786     printf("enter AOD!!");
1787     myPrimaryVertex = ((AliAODEvent*)fESD)->GetPrimaryVertex();
1788     if (!myPrimaryVertex) return;
1789
1790     lPrimaryVtxPosition[0] = myPrimaryVertex->GetX();
1791     lPrimaryVtxPosition[1] = myPrimaryVertex->GetY();
1792     lPrimaryVtxPosition[2] = myPrimaryVertex->GetZ();
1793
1794     // Cut on SPD vertex and fill histo Nevents: FIX it !
1795
1796      // Tracks cuts FIX IT !
1797
1798     // FIX it !!!
1799     lMagneticField = 999;
1800
1801   }
1802
1803
1804   fHistPrimaryVertexX->Fill(lPrimaryVtxPosition[0]);
1805   fHistPrimaryVertexY->Fill(lPrimaryVtxPosition[1]);
1806   fHistPrimaryVertexZ->Fill(lPrimaryVtxPosition[2]);
1807   //Double_t lrcPrimVtxR = TMath::Sqrt(lPrimaryVtxPosition[0]*lPrimaryVtxPosition[0]+lPrimaryVtxPosition[0]*lPrimaryVtxPosition[0]);
1808
1809   fHistPrimaryVertexResX->Fill(lResPrimaryVtxX);
1810   fHistPrimaryVertexResY->Fill(lResPrimaryVtxY);
1811   fHistPrimaryVertexResZ->Fill(lResPrimaryVtxZ);
1812
1813
1814   //***********************
1815   // AliKF Primary Vertex
1816
1817   AliKFVertex primaryVtxKF( *myPrimaryVertex );
1818   AliKFParticle::SetField(lMagneticField);
1819
1820
1821   //************************************
1822   // PID
1823
1824   AliESDpid *fESDpid = new AliESDpid(); // FIXME delete
1825   fESDpid->GetTPCResponse().SetBetheBlochParameters(fAlephParameters[0],fAlephParameters[1],fAlephParameters[2],fAlephParameters[3],fAlephParameters[4]);
1826
1827
1828
1829   //***Rerun the V0 finder
1830
1831 //  fESD->ResetV0s();
1832 //  AliV0vertexer v0Vertexer;
1833 //  v0Vertexer.SetCuts(fCuts);
1834 //  v0Vertexer.Tracks2V0vertices(fESD);
1835
1836   //*************************
1837   // V0 loop
1838
1839   nv0sTot = fESD->GetNumberOfV0s();
1840   if (!nv0sTot) fHistNumberEvents->Fill(6.5);
1841
1842   for (Int_t iV0 = 0; iV0 < nv0sTot; iV0++) {
1843
1844     // ALiKF
1845     AliKFParticle* negPiKF = NULL;
1846     AliKFParticle* posPiKF = NULL;
1847     AliKFParticle* posPKF  = NULL;
1848     AliKFParticle* negAPKF = NULL;
1849
1850
1851     lIndexPosMother     = 0; lIndexNegMother     = 0; lIndexMotherOfMother       = 0;
1852     lCheckPIdK0Short    = 0; lCheckMcK0Short     = 0; lCheckSecondaryK0s         = 0;
1853     lCheckPIdLambda     = 0; lCheckMcLambda      = 0; lCheckSecondaryLambda      = 0;
1854     lCheckPIdAntiLambda = 0; lCheckMcAntiLambda  = 0; lCheckSecondaryAntiLambda  = 0;
1855     lComeFromSigma      = -1;
1856
1857
1858     if(fAnalysisType == "ESD") {
1859
1860
1861       AliESDv0 *v0 = ((AliESDEvent*)fESD)->GetV0(iV0);
1862       if (!v0) continue;
1863
1864       // Primary vertex:
1865       fHistPrimaryVertexPosXV0events->Fill(lPrimaryVtxPosition[0]);
1866       fHistPrimaryVertexPosYV0events->Fill(lPrimaryVtxPosition[1]);
1867       fHistPrimaryVertexPosZV0events->Fill(lPrimaryVtxPosition[2]);
1868
1869       // V0's Daughters
1870       lIndexTrackPos = TMath::Abs(v0->GetPindex());
1871       lIndexTrackNeg = TMath::Abs(v0->GetNindex());
1872       AliESDtrack *myTrackPosTest = ((AliESDEvent*)fESD)->GetTrack(lIndexTrackPos);
1873       AliESDtrack *myTrackNegTest = ((AliESDEvent*)fESD)->GetTrack(lIndexTrackNeg);
1874       if (!myTrackPosTest || !myTrackNegTest) {
1875         // FIXME: shouldn't this be fatal?
1876         Printf("strange analysis::UserExec:: Error:Could not retreive one of the daughter track\n");
1877         continue;
1878       }
1879       // Remove like-sign
1880       if ( myTrackPosTest->GetSign() == myTrackNegTest->GetSign()){
1881         // FIXME: how can this happen?
1882         continue;
1883       }
1884
1885       //    FIXME: are the GetParamN/GetParamP reliable? If so, why do you need to check the sign below?
1886       // VO's main characteristics to check the reconstruction cuts
1887       lOnFlyStatus       = v0->GetOnFlyStatus();
1888       lChi2V0            = v0->GetChi2V0();
1889       lDcaV0Daughters    = v0->GetDcaV0Daughters();
1890       lDcaV0ToPrimVertex = v0->GetD(lPrimaryVtxPosition[0],lPrimaryVtxPosition[1],lPrimaryVtxPosition[2]);
1891       lV0cosPointAngle   = v0->GetV0CosineOfPointingAngle(lPrimaryVtxPosition[0],lPrimaryVtxPosition[1], lPrimaryVtxPosition[2]);
1892
1893       v0->GetXYZ(lV0Position[0], lV0Position[1], lV0Position[2]);
1894
1895       lV0Radius      = TMath::Sqrt(lV0Position[0]*lV0Position[0]+lV0Position[1]*lV0Position[1]);
1896       lV0DecayLength = TMath::Sqrt(TMath::Power(lV0Position[0] - lPrimaryVtxPosition[0],2) +
1897                                    TMath::Power(lV0Position[1] - lPrimaryVtxPosition[1],2) +
1898                                    TMath::Power(lV0Position[2] - lPrimaryVtxPosition[2],2 ));
1899
1900
1901       if( myTrackPosTest->GetSign() ==1){
1902
1903         myTrackPos = ((AliESDEvent*)fESD)->GetTrack(lIndexTrackPos);
1904         myTrackNeg = ((AliESDEvent*)fESD)->GetTrack(lIndexTrackNeg);
1905
1906         // Daughters' momentum;
1907         v0->GetPPxPyPz(lMomPos[0],lMomPos[1],lMomPos[2]);
1908         v0->GetNPxPyPz(lMomNeg[0],lMomNeg[1],lMomNeg[2]);
1909
1910         if (negPiKF) delete negPiKF; negPiKF=NULL;
1911         if (posPiKF) delete posPiKF; posPiKF=NULL;
1912         if (posPKF)  delete posPKF;  posPKF=NULL;
1913         if (negAPKF) delete negAPKF; negAPKF=NULL;
1914
1915         negPiKF = new AliKFParticle( *(v0->GetParamN()) ,-211);
1916         posPiKF = new AliKFParticle( *(v0->GetParamP()) ,211);
1917         posPKF  = new AliKFParticle( *(v0->GetParamP()) ,2212);
1918         negAPKF = new AliKFParticle( *(v0->GetParamN()) ,-2212);
1919
1920       }
1921
1922       if( myTrackPosTest->GetSign() ==-1){
1923
1924         myTrackPos = ((AliESDEvent*)fESD)->GetTrack(lIndexTrackNeg);
1925         myTrackNeg = ((AliESDEvent*)fESD)->GetTrack(lIndexTrackPos);
1926
1927         // Daughters' momentum;
1928         v0->GetPPxPyPz(lMomNeg[0],lMomNeg[1],lMomNeg[2]);
1929         v0->GetNPxPyPz(lMomPos[0],lMomPos[1],lMomPos[2]);
1930
1931         if (negPiKF) delete negPiKF; negPiKF=NULL;
1932         if (posPiKF) delete posPiKF; posPiKF=NULL;
1933         if (posPKF)  delete posPKF;  posPKF=NULL;
1934         if (negAPKF) delete negAPKF; negAPKF=NULL;
1935
1936         negPiKF = new AliKFParticle( *(v0->GetParamP()) ,-211);
1937         posPiKF = new AliKFParticle( *(v0->GetParamN()) ,211);
1938         posPKF  = new AliKFParticle( *(v0->GetParamN()) ,2212);
1939         negAPKF = new AliKFParticle( *(v0->GetParamP()) ,-2212);
1940
1941       }
1942
1943       lLabelTrackPos = (UInt_t)TMath::Abs(myTrackPos->GetLabel());
1944       lLabelTrackNeg = (UInt_t)TMath::Abs(myTrackNeg->GetLabel());
1945
1946       // Daughters Pt and P:
1947       lPtPos = TMath::Sqrt(lMomPos[0]*lMomPos[0] + lMomPos[1]*lMomPos[1]);
1948       lPtNeg = TMath::Sqrt(lMomNeg[0]*lMomNeg[0] + lMomNeg[1]*lMomNeg[1]);
1949
1950       lPPos = TMath::Sqrt(lMomPos[0]*lMomPos[0] + lMomPos[1]*lMomPos[1] + lMomPos[2]*lMomPos[2]);
1951       lPNeg = TMath::Sqrt(lMomNeg[0]*lMomNeg[0] + lMomNeg[1]*lMomNeg[1] + lMomNeg[2]*lMomNeg[2]);
1952
1953       // Inner Wall parameter (used for pid):
1954       const AliExternalTrackParam *myInnerWallTrackPos = myTrackPos->GetInnerParam();
1955       if(myInnerWallTrackPos) lMomInnerWallPos = myInnerWallTrackPos->GetP();
1956       const AliExternalTrackParam *myInnerWallTrackNeg = myTrackNeg->GetInnerParam();
1957       if(myInnerWallTrackNeg) lMomInnerWallNeg = myInnerWallTrackNeg->GetP();
1958
1959       // DCA between daughter and Primary Vertex:
1960       if (myTrackPos) lDcaPosToPrimVertex = TMath::Abs(myTrackPos->GetD(lPrimaryVtxPosition[0],lPrimaryVtxPosition[1],lMagneticField) );
1961
1962       if (myTrackNeg) lDcaNegToPrimVertex = TMath::Abs(myTrackNeg->GetD(lPrimaryVtxPosition[0],lPrimaryVtxPosition[1],lMagneticField) );
1963
1964       // Quality tracks cuts:
1965       if ( !(myTracksCuts->IsSelected(myTrackPos)) || !(myTracksCuts->IsSelected(myTrackNeg)) ) {
1966         if (negPiKF) delete negPiKF; negPiKF=NULL;
1967         if (posPiKF) delete posPiKF; posPiKF=NULL;
1968         if (posPKF)  delete posPKF;  posPKF=NULL;
1969         if (negAPKF) delete negAPKF; negAPKF=NULL;
1970         continue;
1971       }
1972       // Armenteros variables:
1973       lAlphaV0      =  v0->AlphaV0();
1974       lPtArmV0      =  v0->PtArmV0();
1975
1976       // Pseudorapidity:
1977       lV0Eta = v0->Eta();
1978
1979       // PID
1980       if (fUsePID.Contains("withPID")) {
1981         nSigmaPosPion   = TMath::Abs(fESDpid->NumberOfSigmasTPC(myTrackPos,AliPID::kPion));
1982
1983         nSigmaNegPion   = TMath::Abs(fESDpid->NumberOfSigmasTPC(myTrackNeg,AliPID::kPion));
1984
1985         nSigmaPosProton = TMath::Abs(fESDpid->NumberOfSigmasTPC(myTrackPos,AliPID::kProton));
1986
1987         nSigmaNegProton = TMath::Abs(fESDpid->NumberOfSigmasTPC(myTrackNeg,AliPID::kProton));
1988       }
1989       else {
1990         nSigmaPosPion = 0; nSigmaNegPion =0; nSigmaPosProton = 0; nSigmaNegProton= 0;
1991       }
1992
1993       // Monte-Carlo particle associated to reconstructed particles:
1994       if (fAnalysisMC) {
1995         //if (lLabelTrackPos < 0 || lLabelTrackNeg < 0) continue;
1996         TParticle  *lMCESDPartPos  = stack->Particle(lLabelTrackPos);
1997         if(!lMCESDPartPos) {
1998           Printf("no MC particle for positive and/or negative daughter\n");
1999
2000         if (negPiKF) delete negPiKF; negPiKF=NULL;
2001         if (posPiKF) delete posPiKF; posPiKF=NULL;
2002         if (posPKF)  delete posPKF;  posPKF=NULL;
2003         if (negAPKF) delete negAPKF; negAPKF=NULL;
2004         continue;
2005
2006         }
2007         TParticle  *lMCESDPartNeg  = stack->Particle(lLabelTrackNeg);
2008         if (!lMCESDPartNeg) {
2009         if (negPiKF) delete negPiKF; negPiKF=NULL;
2010         if (posPiKF) delete posPiKF; posPiKF=NULL;
2011         if (posPKF)  delete posPKF;  posPKF=NULL;
2012         if (negAPKF) delete negAPKF; negAPKF=NULL;
2013         continue;
2014     }
2015         lPDGCodePosDaughter = lMCESDPartPos->GetPdgCode();
2016         lPDGCodeNegDaughter = lMCESDPartNeg->GetPdgCode();
2017         lIndexPosMother = lMCESDPartPos->GetFirstMother();
2018         lIndexNegMother = lMCESDPartNeg->GetFirstMother();
2019
2020         if (lIndexPosMother == -1) {
2021
2022
2023                                         lPdgcodeMother = 0;
2024                                         lIndexMotherOfMother = 0;
2025                                         mcPosX = 0;
2026                                         mcPosY = 0;
2027                                         mcPosZ = 0;
2028                                         mcPosR = 0;
2029                                         mcPosMotherX = 0;
2030                                         mcPosMotherY = 0;
2031                                         mcPosMotherZ = 0;
2032                                         mcPosMotherR = 0;
2033                                         mcMotherPt = 1;
2034                                         }
2035
2036                 else {
2037
2038         TParticle  *lMCESDMother    = stack->Particle(lIndexPosMother);
2039         if (!lMCESDMother) {
2040         if (negPiKF) delete negPiKF; negPiKF=NULL;
2041         if (posPiKF) delete posPiKF; posPiKF=NULL;
2042         if (posPKF)  delete posPKF;  posPKF=NULL;
2043         if (negAPKF) delete negAPKF; negAPKF=NULL;
2044         continue;
2045     }
2046         lPdgcodeMother         = lMCESDMother->GetPdgCode();
2047         lIndexMotherOfMother   = lMCESDMother->GetFirstMother();
2048         if (lIndexMotherOfMother ==-1) lPdgcodeMotherOfMother = 0;
2049         else {
2050           TParticle  *lMCESDMotherOfMother    = stack->Particle(lIndexMotherOfMother);
2051           if (!lMCESDMotherOfMother) {
2052         if (negPiKF) delete negPiKF; negPiKF=NULL;
2053         if (posPiKF) delete posPiKF; posPiKF=NULL;
2054         if (posPKF)  delete posPKF;  posPKF=NULL;
2055         if (negAPKF) delete negAPKF; negAPKF=NULL;
2056         continue;
2057     }
2058           lPdgcodeMotherOfMother = lMCESDMotherOfMother->GetPdgCode();
2059         }
2060
2061         mcPosX = lMCESDPartPos->Vx();
2062         mcPosY = lMCESDPartPos->Vy();
2063         mcPosZ = lMCESDPartPos->Vz();
2064         mcPosR = TMath::Sqrt(mcPosX*mcPosX+mcPosY*mcPosY);
2065         mcPosMotherX = lMCESDMother->Vx();
2066         mcPosMotherY = lMCESDMother->Vy();
2067         mcPosMotherZ = lMCESDMother->Vz();
2068         mcPosMotherR = TMath::Sqrt(mcPosMotherX*mcPosMotherX+mcPosMotherY*mcPosMotherY);
2069
2070         mcMotherPt   = lMCESDMother->Pt();
2071       }
2072      }
2073     } // end ESD condition
2074
2075
2076
2077     else if(fAnalysisType == "AOD") {
2078
2079       AliAODv0     *myAODv0 = ((AliAODEvent*)fESD)->GetV0(iV0);
2080       if (!myAODv0) continue;
2081
2082       // Primary vertex:
2083       fHistPrimaryVertexPosXV0events->Fill(lPrimaryVtxPosition[0]);
2084       fHistPrimaryVertexPosYV0events->Fill(lPrimaryVtxPosition[1]);
2085       fHistPrimaryVertexPosZV0events->Fill(lPrimaryVtxPosition[2]);
2086
2087
2088       //Multiplicity:
2089       if(lOnFlyStatus == fUseOnTheFly) nv0s++;
2090
2091       // V0's Daughters
2092       lIndexTrackPos = TMath::Abs(myAODv0->GetPosID());
2093       lIndexTrackNeg = TMath::Abs(myAODv0->GetNegID());
2094
2095       AliVParticle  *lVPartPosTest  = ((AliVEvent*)fESD)->GetTrack(lIndexTrackPos);
2096       AliVParticle  *lVPartNegTest  = ((AliVEvent*)fESD)->GetTrack(lIndexTrackNeg);
2097       //AliAODTrack  *lVPartPos  = ((AliAODEvent*)fESD)->GetTrack(lIndexTrackPos);
2098       //AliAODTrack  *lVPartNeg  = ((AliAODEvent*)fESD)->GetTrack(lIndexTrackNeg);
2099
2100       if (!lVPartPosTest ||(!lVPartNegTest )) {
2101         Printf("strange analysis::UserExec:: Could not retreive one of the daughter track\n");
2102         continue;
2103       }
2104
2105       // Quality cuts:
2106       // TO DO !!!!!!!
2107
2108       // TPC refit condition (done during reconstruction for Offline but not for On-the-fly)
2109       //if( !(lVPartPosTest->GetStatus() & AliAODTrack::kTPCrefit)) continue;
2110       //if( !(lVPartNegTest->GetStatus() & AliAODTrack::kTPCrefit)) continue;
2111
2112
2113       lDcaPosToPrimVertex = myAODv0->DcaPosToPrimVertex();
2114       lDcaNegToPrimVertex = myAODv0->DcaNegToPrimVertex();
2115       lOnFlyStatus        = myAODv0->GetOnFlyStatus();
2116       lChi2V0             = myAODv0->Chi2V0();
2117       lDcaV0Daughters     = myAODv0->DcaV0Daughters();
2118       lDcaV0ToPrimVertex  = myAODv0->DcaV0ToPrimVertex();
2119       lV0DecayLength      = myAODv0->DecayLengthV0(lPrimaryVtxPosition);
2120       lV0cosPointAngle    = myAODv0->CosPointingAngle(lPrimaryVtxPosition);
2121       lV0Radius           = myAODv0->RadiusV0();
2122
2123       if( lVPartPosTest->Charge() ==1){
2124
2125         lVPartPos = ((AliVEvent*)fESD)->GetTrack(lIndexTrackPos);
2126         lVPartNeg = ((AliVEvent*)fESD)->GetTrack(lIndexTrackNeg);
2127
2128
2129         if (negPiKF) delete negPiKF; negPiKF=NULL;
2130         if (posPiKF) delete posPiKF; posPiKF=NULL;
2131         if (posPKF)  delete posPKF; posPKF=NULL;
2132         if (negAPKF) delete negAPKF; negAPKF=NULL;
2133
2134         //negPiKF = new AliKFParticle( *(myAODv0->GetParamN()) ,-211);
2135         //posPiKF = new AliKFParticle( *(myAODv0->GetParamP()) ,211);
2136         //posPKF  = new AliKFParticle( *(myAODv0->GetParamP()) ,2212);
2137         //negAPKF = new AliKFParticle( *(myAODv0->GetParamN()) ,-2212);
2138         // TO DO !!!!!!
2139         negPiKF = NULL;
2140         posPiKF = NULL;
2141         posPKF  = NULL;
2142         negAPKF = NULL;
2143
2144       }
2145
2146       if( lVPartPosTest->Charge() ==-1){
2147
2148         lVPartPos = ((AliVEvent*)fESD)->GetTrack(lIndexTrackNeg);
2149         lVPartNeg = ((AliVEvent*)fESD)->GetTrack(lIndexTrackPos);
2150
2151         if (negPiKF) delete negPiKF; negPiKF=NULL;
2152         if (posPiKF) delete posPiKF; posPiKF=NULL;
2153         if (posPKF)  delete posPKF; posPKF=NULL;
2154         if (negAPKF) delete negAPKF; negAPKF=NULL;
2155
2156         //negPiKF = new AliKFParticle( *(myAODv0->GetParamP()) ,-211);
2157         //posPiKF = new AliKFParticle( *(myAODv0->GetParamN()) ,211);
2158         //posPKF  = new AliKFParticle( *(myAODv0->GetParamN()) ,2212);
2159         //negAPKF = new AliKFParticle( *(myAODv0->GetParamP()) ,-2212);
2160         negPiKF = NULL;
2161         posPiKF = NULL;
2162         posPKF  = NULL;
2163         negAPKF = NULL;
2164       }
2165
2166       lLabelTrackPos  = TMath::Abs(lVPartPos->GetLabel());
2167       lLabelTrackNeg  = TMath::Abs(lVPartNeg->GetLabel());
2168
2169       // Armenteros variables:
2170       lAlphaV0   = myAODv0->AlphaV0();
2171       lPtArmV0   = myAODv0->PtArmV0();
2172
2173       // Pseudorapidity:
2174       lV0Eta = myAODv0->PseudoRapV0();
2175
2176       // PID not accessible with AOD !
2177       nSigmaPosPion = 0; nSigmaNegPion =0; nSigmaPosProton = 0; nSigmaNegProton= 0;
2178
2179
2180       // Monte-Carlo particle associated to reconstructed particles:
2181       if (fAnalysisMC) {
2182         AliAODMCParticle *lMCAODPartPos = (AliAODMCParticle*)mcArray->At(lLabelTrackPos);
2183         if (!lMCAODPartPos) {
2184         if (negPiKF) delete negPiKF; negPiKF=NULL;
2185         if (posPiKF) delete posPiKF; posPiKF=NULL;
2186         if (posPKF)  delete posPKF;  posPKF=NULL;
2187         if (negAPKF) delete negAPKF; negAPKF=NULL;
2188         continue;
2189     }
2190         AliAODMCParticle *lMCAODPartNeg = (AliAODMCParticle*)mcArray->At(lLabelTrackNeg);
2191         if(!lMCAODPartNeg)
2192          // Printf("strange analysis::UserExec:no MC particle for negative daughter\n");
2193           {
2194         if (negPiKF) delete negPiKF; negPiKF=NULL;
2195         if (posPiKF) delete posPiKF; posPiKF=NULL;
2196         if (posPKF)  delete posPKF;  posPKF=NULL;
2197         if (negAPKF) delete negAPKF; negAPKF=NULL;
2198         continue;
2199
2200         }
2201         lPDGCodePosDaughter = lMCAODPartPos->GetPdgCode();
2202         lPDGCodeNegDaughter = lMCAODPartNeg->GetPdgCode();
2203         lIndexPosMother = lMCAODPartPos->GetMother();
2204         lIndexNegMother = lMCAODPartNeg->GetMother();
2205
2206         AliAODMCParticle *lMCAODMother = (AliAODMCParticle*)mcArray->At(lIndexPosMother);
2207         lPdgcodeMother = lMCAODMother->GetPdgCode();
2208         lIndexMotherOfMother  = lMCAODMother->GetMother();
2209         if (lIndexMotherOfMother ==-1) lPdgcodeMotherOfMother = 0;
2210         else {
2211           lPdgcodeMotherOfMother   = ((AliAODMCParticle*)mcArray->At(lIndexMotherOfMother))->GetPdgCode();
2212         }
2213
2214         mcPosX = lMCAODPartPos->Xv();
2215         mcPosY = lMCAODPartPos->Yv();
2216         mcPosZ = lMCAODPartPos->Zv();
2217         mcPosR = TMath::Sqrt(mcPosX*mcPosX+mcPosY*mcPosY);
2218         mcPosMotherX = lMCAODMother->Xv();
2219         mcPosMotherY = lMCAODMother->Yv();
2220         mcPosMotherZ = lMCAODMother->Zv();
2221         mcPosMotherR = TMath::Sqrt(mcPosMotherX*mcPosMotherX+mcPosMotherY*mcPosMotherY);
2222         mcMotherPt   = lMCAODMother->Pt();
2223       }
2224
2225     } // end AOD condition
2226
2227
2228     // Multiplicity:
2229     if(lOnFlyStatus == fUseOnTheFly) nv0s++;
2230
2231     // Daughter momentum cut: ! FIX it in case of AOD !
2232     if ( (lPtPos  < cutMinPtDaughter ) ||
2233          (lPtNeg  < cutMinPtDaughter )
2234          ) {
2235       if (negPiKF) delete negPiKF; negPiKF=NULL;
2236       if (posPiKF) delete posPiKF; posPiKF=NULL;
2237       if (posPKF)  delete posPKF;  posPKF=NULL;
2238       if (negAPKF) delete negAPKF; negAPKF=NULL;
2239       continue;
2240     }
2241
2242     AliKFParticle v0K0sKF;
2243     v0K0sKF+=(*negPiKF);
2244     v0K0sKF+=(*posPiKF);
2245     v0K0sKF.SetProductionVertex(primaryVtxKF);
2246
2247     AliKFParticle v0LambdaKF;
2248     v0LambdaKF+=(*negPiKF);
2249     v0LambdaKF+=(*posPKF);
2250     v0LambdaKF.SetProductionVertex(primaryVtxKF);
2251
2252     AliKFParticle v0AntiLambdaKF;
2253     v0AntiLambdaKF+=(*posPiKF);
2254     v0AntiLambdaKF+=(*negAPKF);
2255     v0AntiLambdaKF.SetProductionVertex(primaryVtxKF);
2256
2257     // Invariant mass
2258     lInvMassK0s        = v0K0sKF.GetMass();
2259     lInvMassLambda     = v0LambdaKF.GetMass();
2260     lInvMassAntiLambda = v0AntiLambdaKF.GetMass();
2261
2262     // Rapidity:
2263     lRapK0s        = 0.5*TMath::Log((v0K0sKF.GetE()+v0K0sKF.GetPz())/(v0K0sKF.GetE()-v0K0sKF.GetPz()+1.e-13));
2264     lRapLambda     = 0.5*TMath::Log((v0LambdaKF.GetE()+v0LambdaKF.GetPz())/(v0LambdaKF.GetE()-v0LambdaKF.GetPz()+1.e-13));
2265     lRapAntiLambda = 0.5*TMath::Log((v0AntiLambdaKF.GetE()+v0AntiLambdaKF.GetPz())/(v0AntiLambdaKF.GetE()-v0AntiLambdaKF.GetPz()+1.e-13));
2266
2267     // Pseudo-rapidity
2268     lEtaK0s     = v0K0sKF.GetEta();
2269     lEtaLambda  = v0LambdaKF.GetEta();
2270     lEtaAntiLambda  = v0AntiLambdaKF.GetEta();
2271
2272     // Pz:
2273     lPzK0s        = v0K0sKF.GetPz();
2274     lPzLambda     = v0LambdaKF.GetPz();
2275     lPzAntiLambda = v0AntiLambdaKF.GetPz();
2276
2277     // Pt:
2278     lPtK0s        = v0K0sKF.GetPt();
2279     lPtLambda     = v0LambdaKF.GetPt();
2280     lPtAntiLambda = v0AntiLambdaKF.GetPt();
2281
2282     if (lPtK0s==0) {
2283         if (negPiKF) delete negPiKF; negPiKF=NULL;
2284         if (posPiKF) delete posPiKF; posPiKF=NULL;
2285         if (posPKF)  delete posPKF;  posPKF=NULL;
2286         if (negAPKF) delete negAPKF; negAPKF=NULL;
2287         // FIXME: should you really continue here and below?
2288         continue;
2289     }
2290     if (lPtLambda==0) {
2291       if (negPiKF) delete negPiKF; negPiKF=NULL;
2292       if (posPiKF) delete posPiKF; posPiKF=NULL;
2293       if (posPKF)  delete posPKF;  posPKF=NULL;
2294       if (negAPKF) delete negAPKF; negAPKF=NULL;
2295       continue;
2296     }
2297     // Pt Resolution
2298     deltaPtK0s        = (lPtK0s - mcMotherPt)/mcMotherPt;
2299     deltaPtLambda     = (lPtLambda - mcMotherPt)/mcMotherPt;
2300     deltaPtAntiLambda = (lPtAntiLambda - mcMotherPt)/mcMotherPt;
2301
2302     // KF Chi2
2303     lChi2KFK0s        = v0K0sKF.GetChi2();
2304     lChi2KFLambda     = v0LambdaKF.GetChi2();
2305     lChi2KFAntiLambda = v0AntiLambdaKF.GetChi2();
2306
2307     // Reconstructed Position
2308     rcPosXK0s = v0K0sKF.GetX();
2309     rcPosYK0s = v0K0sKF.GetY();
2310     rcPosZK0s = v0K0sKF.GetZ();
2311     rcPosRK0s = TMath::Sqrt(rcPosXK0s*rcPosXK0s+rcPosYK0s*rcPosYK0s);
2312
2313     rcPosXLambda = v0LambdaKF.GetX();
2314     rcPosYLambda = v0LambdaKF.GetY();
2315     rcPosZLambda = v0LambdaKF.GetZ();
2316     rcPosRLambda = TMath::Sqrt(rcPosXLambda*rcPosXLambda+rcPosYLambda*rcPosYLambda);
2317
2318     rcPosXAntiLambda = v0AntiLambdaKF.GetX();
2319     rcPosYAntiLambda = v0AntiLambdaKF.GetY();
2320     rcPosZAntiLambda = v0AntiLambdaKF.GetZ();
2321     rcPosRAntiLambda = TMath::Sqrt(rcPosXAntiLambda*rcPosXAntiLambda+rcPosYAntiLambda*rcPosYAntiLambda);
2322
2323     TVector3 momPos(lMomPos[0],lMomPos[1],lMomPos[2]);
2324     TVector3 momNeg(lMomNeg[0],lMomNeg[1],lMomNeg[2]);
2325     TVector3 momTotK0s(v0K0sKF.GetPx(),v0K0sKF.GetPy(),v0K0sKF.GetPz());
2326     TVector3 momTotLambda(v0LambdaKF.GetPx(),v0LambdaKF.GetPy(),v0LambdaKF.GetPz());
2327     TVector3 momTotAntiLambda(v0AntiLambdaKF.GetPx(),v0AntiLambdaKF.GetPy(),v0AntiLambdaKF.GetPz());
2328
2329     lQlPos = momPos.Dot(momTotK0s)/momTotK0s.Mag();
2330     lQlNeg = momNeg.Dot(momTotK0s)/momTotK0s.Mag();
2331     lAlphaV0K0s = 1.-2./(1.+lQlPos/lQlNeg);
2332     lQlPos = momPos.Dot(momTotLambda)/momTotLambda.Mag();
2333     lQlNeg = momNeg.Dot(momTotLambda)/momTotLambda.Mag();
2334     lAlphaV0Lambda = 1.-2./(1.+lQlPos/lQlNeg);
2335     lQlPos = momPos.Dot(momTotAntiLambda)/momTotAntiLambda.Mag();
2336     lQlNeg = momNeg.Dot(momTotAntiLambda)/momTotAntiLambda.Mag();
2337     lAlphaV0AntiLambda = 1.-2./(1.+lQlPos/lQlNeg);
2338
2339     lPtArmV0K0s = momPos.Perp(momTotK0s);
2340     lPtArmV0Lambda = momPos.Perp(momTotLambda);
2341     lPtArmV0AntiLambda = momPos.Perp(momTotAntiLambda);
2342
2343     // Look for associated particles:
2344     if (fAnalysisMC) {
2345       if( (lIndexPosMother==-1) || (lIndexNegMother==-1) ) {
2346         fHistMCDaughterTrack->Fill(1);
2347       }
2348
2349       else if( ( (lPDGCodePosDaughter==+211) && (lPDGCodeNegDaughter==-211) )
2350                ) {
2351         lCheckPIdK0Short    = 1;
2352         fHistMCDaughterTrack->Fill(3);
2353         if ( (lIndexPosMother==lIndexNegMother) &&
2354              (lPdgcodeMother==310) ) {
2355           if (lIndexPosMother <= lNbMCPrimary) lCheckMcK0Short  = 1;
2356           else lCheckSecondaryK0s = 1;
2357         }
2358       }
2359       else if( ( (lPDGCodePosDaughter==+2212) && (lPDGCodeNegDaughter==-211)  )
2360                ) {
2361         lCheckPIdLambda     = 1;
2362         fHistMCDaughterTrack->Fill(5);
2363         if ( (lIndexPosMother==lIndexNegMother) &&
2364              (lPdgcodeMother==3122)  ){
2365           if ( ( TMath::Abs(lPdgcodeMotherOfMother) == 3212) ||
2366                ( TMath::Abs(lPdgcodeMotherOfMother) == 3224) ||
2367                ( TMath::Abs(lPdgcodeMotherOfMother) == 3214) ||
2368                ( TMath::Abs(lPdgcodeMotherOfMother) == 3114)
2369                ) lComeFromSigma = 1;
2370           else lComeFromSigma = 0;
2371           if ( (lIndexPosMother <= lNbMCPrimary) ||
2372              ( ( lIndexPosMother > lNbMCPrimary) && (lComeFromSigma) )
2373              ) lCheckMcLambda  = 1;
2374           else lCheckSecondaryLambda    = 1;
2375         }
2376       }
2377       else if( ( (lPDGCodePosDaughter==211)   && (lPDGCodeNegDaughter==-2212) )
2378                ) {
2379         lCheckPIdAntiLambda = 1;
2380         fHistMCDaughterTrack->Fill(7);
2381         if ( (lIndexPosMother==lIndexNegMother) &&
2382              (lPdgcodeMother==-3122) ) {
2383           if ( ( TMath::Abs(lPdgcodeMotherOfMother) == 3212) ||
2384                ( TMath::Abs(lPdgcodeMotherOfMother) == 3224) ||
2385                ( TMath::Abs(lPdgcodeMotherOfMother) == 3214) ||
2386                ( TMath::Abs(lPdgcodeMotherOfMother) == 3114)
2387                ) lComeFromSigma = 1;
2388           else lComeFromSigma = 0;
2389           if ( (lIndexPosMother <= lNbMCPrimary) ||
2390              ( ( lIndexPosMother > lNbMCPrimary) && (lComeFromSigma) )
2391              ) lCheckMcAntiLambda  = 1;
2392           else lCheckSecondaryAntiLambda = 1;
2393         }
2394       }
2395
2396       // Gamma conversion
2397       else if ( (lPDGCodePosDaughter==11) &&
2398                 (lPDGCodeNegDaughter==-11) &&
2399                 (lPdgcodeMother==22 ) )
2400         lCheckGamma = 1;
2401     } // end "look for associated particles
2402
2403
2404     // Cuts:
2405 /*    if (fUseCut.Contains("yes")) {
2406       if ( (lDcaPosToPrimVertex < 0.036 ) ||
2407            (lDcaNegToPrimVertex < 0.036 ) ||
2408            (lDcaV0Daughters     > 0.5   ) ||
2409            (lV0cosPointAngle    < 0.999 )
2410            )
2411         continue;
2412     }
2413 */
2414
2415 /*
2416       if ( (lDcaV0Daughters     > 0.3   ) ||
2417            (lV0cosPointAngle    < 0.998 )
2418
2419            )    continue;
2420 */
2421     // PID condition:
2422     lCheckPIDK0sPosDaughter        = 0, lCheckPIDK0sNegDaughter        = 0;
2423     lCheckPIDLambdaPosDaughter     = 0, lCheckPIDLambdaNegDaughter     = 0;
2424     lCheckPIDAntiLambdaPosDaughter = 0, lCheckPIDAntiLambdaNegDaughter = 0;
2425
2426     if (lMomInnerWallPos < lLimitPPID) {
2427       if (nSigmaPosPion < cutNSigmaLowP)   {
2428         lCheckPIDK0sPosDaughter        = 1;
2429         lCheckPIDAntiLambdaPosDaughter = 1;
2430       }
2431       if (nSigmaPosProton < cutNSigmaLowP) lCheckPIDLambdaPosDaughter    = 1;
2432     }
2433
2434     else if (lMomInnerWallPos > lLimitPPID) {
2435       if (nSigmaPosPion < cutNSigmaHighP)   {
2436         lCheckPIDK0sPosDaughter        = 1;
2437         lCheckPIDAntiLambdaPosDaughter = 1;
2438       }
2439       if (nSigmaPosProton < cutNSigmaHighP) lCheckPIDLambdaPosDaughter    = 1;
2440     }
2441
2442     if (lMomInnerWallNeg < lLimitPPID) {
2443       if (nSigmaNegPion < cutNSigmaLowP)    {
2444         lCheckPIDK0sNegDaughter       = 1;
2445         lCheckPIDLambdaNegDaughter    = 1;
2446       }
2447       if (nSigmaNegProton < cutNSigmaLowP)  lCheckPIDAntiLambdaNegDaughter = 1;
2448
2449     }
2450     else if (lMomInnerWallNeg > lLimitPPID) {
2451       if (nSigmaNegPion < cutNSigmaHighP)   {
2452         lCheckPIDK0sNegDaughter       = 1;
2453         lCheckPIDLambdaNegDaughter    = 1;
2454       }
2455       if (nSigmaNegProton < cutNSigmaHighP) lCheckPIDAntiLambdaNegDaughter = 1;
2456     }
2457
2458
2459     //*****************************
2460     // filling histograms
2461
2462     fHistDcaPosToPrimVertex->Fill(lDcaPosToPrimVertex,lOnFlyStatus);
2463     fHistDcaNegToPrimVertex->Fill(lDcaNegToPrimVertex,lOnFlyStatus);
2464     fHistDcaPosToPrimVertexZoom->Fill(lDcaPosToPrimVertex,lOnFlyStatus);
2465     fHistDcaNegToPrimVertexZoom->Fill(lDcaNegToPrimVertex,lOnFlyStatus);
2466     fHistRadiusV0->Fill(lV0Radius,lOnFlyStatus);
2467     fHistDecayLengthV0->Fill(lV0DecayLength,lOnFlyStatus);
2468     fHistDcaV0Daughters->Fill(lDcaV0Daughters,lOnFlyStatus);
2469     fHistChi2->Fill(lChi2V0,lOnFlyStatus);
2470     fHistCosPointAngle->Fill(lV0cosPointAngle,lOnFlyStatus);
2471     if (lV0cosPointAngle >= 0.9) fHistCosPointAngleZoom->Fill(lV0cosPointAngle,lOnFlyStatus);
2472     fHistChi2KFBeforeCutK0s->Fill(lChi2KFK0s,lOnFlyStatus);
2473     fHistChi2KFBeforeCutLambda->Fill(lChi2KFLambda,lOnFlyStatus);
2474     fHistChi2KFBeforeCutAntiLambda->Fill(lChi2KFAntiLambda,lOnFlyStatus);
2475
2476
2477     // Histo versus Rap and armenteros plot
2478     if (lOnFlyStatus == fUseOnTheFly){
2479       if (lCheckMcK0Short) fHistAsMcRapK0->Fill(lRapK0s);
2480       if (lCheckMcLambda) fHistAsMcRapLambda->Fill(lRapLambda);
2481       if (lCheckMcAntiLambda) fHistAsMcRapLambda->Fill(lRapAntiLambda);
2482 //      fHistArmenterosPodolanski->Fill(lAlphaV0,lPtArmV0);
2483 //      fHistDaughterPt->Fill(lPtPos,lPtNeg);
2484     }
2485
2486     // FIXME: associated histos, what are they used for?
2487
2488     // K0s associated histograms in |rap| < lCutRap:
2489
2490
2491     if ( lCheckPIDK0sPosDaughter && lCheckPIDK0sNegDaughter
2492          && (lChi2KFK0s < cutChi2KF) && (TMath::Abs(lPzK0s/lPtK0s)<0.7) )     {
2493
2494
2495       fHistChi2KFAfterCutK0s->Fill(lChi2KFK0s,lOnFlyStatus);
2496
2497       if (TMath::Abs(lRapK0s) < lCutRap) {
2498
2499         fHistNsigmaPosPionK0->Fill(nSigmaPosPion);
2500         fHistNsigmaNegPionK0->Fill(nSigmaNegPion);
2501
2502         if (lOnFlyStatus == fUseOnTheFly){
2503           fHistMassK0->Fill(lInvMassK0s);
2504           fHistMassVsRadiusK0->Fill(rcPosRK0s,lInvMassK0s);
2505           fHistPtVsMassK0->Fill(lInvMassK0s,lPtK0s);
2506
2507
2508 //        fHistMultVsPtVsMassK0->Fill(multiplicity ,lInvMassK0s,lPtK0s);
2509           if(lCheckPIdK0Short) fHistPidMcMassK0->Fill(lInvMassK0s);
2510           if(lCheckMcK0Short) {
2511             fHistAsMcMassK0->Fill(lInvMassK0s);
2512             fHistAsMcPtK0->Fill(lPtK0s);
2513
2514
2515             fHistAsMcPtVsMassK0->Fill(lInvMassK0s,lPtK0s);
2516             if (lPtK0s <= 1) fHistAsMcPtZoomK0->Fill(lPtK0s);
2517             fHistAsMcMassVsRadiusK0->Fill(rcPosRK0s,lInvMassK0s);
2518             fHistAsMcResxK0->Fill(rcPosXK0s-mcPosX);
2519             fHistAsMcResyK0->Fill(rcPosYK0s-mcPosY);
2520             fHistAsMcReszK0->Fill(rcPosZK0s-mcPosZ);
2521             fHistAsMcResrVsRadiusK0->Fill(rcPosRK0s,rcPosRK0s-mcPosR);
2522             fHistAsMcReszVsRadiusK0->Fill(rcPosZK0s,rcPosZK0s-mcPosZ);
2523             fHistAsMcProdRadiusK0->Fill(mcPosMotherR);
2524             fHistAsMcProdRadiusXvsYK0s->Fill(mcPosMotherX,mcPosMotherY);
2525             fHistAsMcResPtK0->Fill(deltaPtK0s);
2526             fHistAsMcResPtVsRapK0->Fill(deltaPtK0s,lRapK0s);
2527             fHistAsMcResPtVsPtK0->Fill(deltaPtK0s,lPtK0s);
2528           }
2529           else if (lCheckSecondaryK0s) {
2530             fHistAsMcSecondaryPtVsRapK0s->Fill(lPtK0s,lRapK0s);
2531             fHistAsMcSecondaryProdRadiusK0s->Fill(mcPosMotherR);
2532             fHistAsMcSecondaryProdRadiusXvsYK0s->Fill(mcPosMotherX,mcPosMotherY);
2533             switch (lPdgcodeMotherOfMother) {
2534             case 130   : fHistAsMcSecondaryMotherPdgCodeK0s->Fill(0.5);break; // K0L
2535             case 321   : fHistAsMcSecondaryMotherPdgCodeK0s->Fill(1.5);break; // K+
2536             case -321  : fHistAsMcSecondaryMotherPdgCodeK0s->Fill(2.5);break; // K-
2537             case -3122 : fHistAsMcSecondaryMotherPdgCodeK0s->Fill(3.5);break; //AntiLambda
2538             default    : fHistAsMcSecondaryMotherPdgCodeK0s->Fill(6.5);break;
2539             }
2540           }
2541         }
2542       } // end rapidity condition
2543     } // end nsigma condition
2544
2545
2546     // Associated Lambda histograms in |rap| < lCutRap
2547
2548
2549     if ( lCheckPIDLambdaPosDaughter && lCheckPIDLambdaNegDaughter
2550          && (lChi2KFLambda < cutChi2KF) && (TMath::Abs(lPzLambda/lPtLambda)<0.7) )  {
2551
2552
2553
2554       fHistChi2KFAfterCutLambda->Fill(lChi2KFLambda,lOnFlyStatus);
2555
2556       if (TMath::Abs(lRapLambda) < lCutRap) {
2557
2558         fHistNsigmaPosProtonLambda->Fill(nSigmaPosProton);
2559         fHistNsigmaNegPionLambda->Fill(nSigmaNegPion);
2560         if (lOnFlyStatus == fUseOnTheFly){
2561           fHistMassLambda->Fill(lInvMassLambda);
2562           fHistMassVsRadiusLambda->Fill(rcPosRLambda,lInvMassLambda);
2563           fHistPtVsMassLambda->Fill(lInvMassLambda,lPtLambda);
2564
2565 //          fHistMultVsPtVsMassLambda->Fill(multiplicity ,lInvMassLambda,lPtLambda);
2566           if(lCheckPIdLambda) fHistPidMcMassLambda->Fill(lInvMassLambda);
2567
2568           if(lCheckMcLambda) {
2569             fHistAsMcMassLambda->Fill(lInvMassLambda);
2570             fHistAsMcPtLambda->Fill(lPtLambda);
2571
2572
2573             fHistAsMcPtVsMassLambda->Fill(lInvMassLambda,lPtLambda);
2574             if (lPtLambda <= 1) fHistAsMcPtZoomLambda->Fill(lPtLambda);
2575             fHistAsMcMassVsRadiusLambda->Fill(rcPosRLambda,lInvMassLambda);
2576             fHistAsMcResxLambda->Fill(rcPosXLambda-mcPosX);
2577             fHistAsMcResyLambda->Fill(rcPosYLambda-mcPosY);
2578             fHistAsMcReszLambda->Fill(rcPosZLambda-mcPosZ);
2579             fHistAsMcResrVsRadiusLambda->Fill(rcPosRLambda,rcPosRLambda-mcPosR);
2580             fHistAsMcReszVsRadiusLambda->Fill(rcPosZLambda,rcPosZLambda-mcPosZ);
2581             fHistAsMcProdRadiusLambda->Fill(mcPosMotherR);
2582             fHistAsMcProdRadiusXvsYLambda->Fill(mcPosMotherX,mcPosMotherY);
2583             fHistAsMcResPtLambda->Fill(deltaPtLambda);
2584             fHistAsMcResPtVsRapLambda->Fill(deltaPtLambda,lRapLambda);
2585             fHistAsMcResPtVsPtLambda->Fill(deltaPtLambda,lPtLambda);
2586             if (lComeFromSigma) fHistAsMcPtLambdaFromSigma->Fill(lPtLambda);
2587             switch (lPdgcodeMotherOfMother) {
2588             case 3222 : fHistAsMcMotherPdgCodeLambda->Fill(0.5); break; // Sigma +
2589             case 3212 : fHistAsMcMotherPdgCodeLambda->Fill(1.5); break; // Sigma 0
2590             case 3112 : fHistAsMcMotherPdgCodeLambda->Fill(2.5); break;// Sigma -
2591             case 3224 : fHistAsMcMotherPdgCodeLambda->Fill(3.5); break;// Sigma(1385) +
2592             case 3214 : fHistAsMcMotherPdgCodeLambda->Fill(4.5); break;// Sigma(1385) 0
2593             case 3114 : fHistAsMcMotherPdgCodeLambda->Fill(5.5); break;// Sigma(1385) -
2594             case 3322 : fHistAsMcMotherPdgCodeLambda->Fill(6.5); break; // Xi 0
2595             case 3312 : fHistAsMcMotherPdgCodeLambda->Fill(7.5); break; // Xi -
2596             case 3334 : fHistAsMcMotherPdgCodeLambda->Fill(8.5); break; // Omega
2597             case -1   : fHistAsMcMotherPdgCodeLambda->Fill(9.5); break;
2598             default   : fHistAsMcMotherPdgCodeLambda->Fill(10.5);break;
2599             }
2600
2601
2602             //printf("found Lambda RC dcaPos=%e dcaNeg=%e dcaDau=%e cosP=%e pT=%e mass=%e\n",lDcaPosToPrimVertex ,lDcaNegToPrimVertex ,lDcaV0Daughters,lV0cosPointAngle,lPtLambda,lInvMassLambda);
2603             //printf("found Lambda RC Pindex=%d  Nindex=%d  Plabel=%d  Nlabel=%d\n\n",lIndexTrackPos,lIndexTrackNeg,lLabelTrackPos,lLabelTrackNeg);
2604
2605           }
2606
2607           else if (lCheckSecondaryLambda) {
2608             fHistAsMcSecondaryPtVsRapLambda->Fill(lPtLambda,lRapLambda);
2609             fHistAsMcSecondaryProdRadiusLambda->Fill(mcPosMotherR);
2610             fHistAsMcSecondaryProdRadiusXvsYLambda->Fill(mcPosMotherX,mcPosMotherY);
2611             if (lComeFromSigma) fHistAsMcSecondaryPtLambdaFromSigma->Fill(lPtLambda);
2612             printf(" lPdgcodeMotherOfMother= %d",lPdgcodeMotherOfMother);
2613             switch (lPdgcodeMotherOfMother) {
2614             case 3222 : fHistAsMcSecondaryMotherPdgCodeLambda->Fill(0.5); break;// Sigma +
2615             case 3212 : fHistAsMcSecondaryMotherPdgCodeLambda->Fill(1.5); break;// Sigma 0
2616             case 3112 : fHistAsMcSecondaryMotherPdgCodeLambda->Fill(2.5); break;// Sigma -
2617             case 3224 : fHistAsMcSecondaryMotherPdgCodeLambda->Fill(3.5); break;// Sigma(1385) +
2618             case 3214 : fHistAsMcSecondaryMotherPdgCodeLambda->Fill(4.5); break;// Sigma(1385) 0
2619             case 3114 : fHistAsMcSecondaryMotherPdgCodeLambda->Fill(5.5); break;// Sigma(1385) -
2620             case 3322 : fHistAsMcSecondaryMotherPdgCodeLambda->Fill(6.5); break; // Xi 0
2621             case 3312 : fHistAsMcSecondaryMotherPdgCodeLambda->Fill(7.5); break; // Xi -
2622             case 3334 : fHistAsMcSecondaryMotherPdgCodeLambda->Fill(8.5); break; // Omega
2623             case -1   : fHistAsMcSecondaryMotherPdgCodeLambda->Fill(9.5); break;
2624             default   : fHistAsMcSecondaryMotherPdgCodeLambda->Fill(10.5);break;
2625             }
2626           }
2627         }
2628       } // end rapidity condition
2629     } //end nsigma condition - lambda
2630
2631
2632     if (negPiKF) delete negPiKF; negPiKF= NULL;
2633     if (posPiKF) delete posPiKF; posPiKF= NULL;
2634     if (posPKF)  delete posPKF;  posPKF = NULL;
2635     if (negAPKF) delete negAPKF; negAPKF= NULL;
2636
2637   } // end V0 loop
2638
2639   fHistV0Multiplicity->Fill(nv0s);
2640   if (fAnalysisType == "ESD") { if(myPrimaryVertex) delete myPrimaryVertex; }
2641
2642   if(myTracksCuts) delete myTracksCuts;
2643
2644   // Post output data
2645   PostData(1, fListHist);
2646   PostData(2, fCentrSelector);
2647 }
2648
2649 //________________________________________________________________________
2650 void AliAnalysisTaskPerformanceStrange::Terminate(Option_t *)
2651 {/*
2652   // Draw result to the screen
2653   // Called once at the end of the query
2654
2655   TList *cRetrievedList = 0x0;
2656   cRetrievedList = (TList*)GetOutputData(1);
2657
2658   if(!cRetrievedList){
2659     AliWarning("ERROR - AliAnalysisTaskPerformanceStrange: output data container list not available\n"); return;
2660   }
2661
2662
2663   fHistV0Multiplicity = dynamic_cast<TH1F*> ( cRetrievedList->FindObject("fHistV0Multiplicity"));
2664   if (!fHistV0Multiplicity) {
2665     Printf("ERROR: fHistV0Multiplicity not available");
2666     return;
2667   }
2668
2669   fHistV0MultiplicityMI = dynamic_cast<TH1F*> ( cRetrievedList->FindObject("fHistV0MultiplicityMI"));
2670   if (!fHistV0MultiplicityMI) {
2671     Printf("ERROR: fHistV0MultiplicityMI not available");
2672     return;
2673   }
2674
2675   TCanvas *canPerformanceStrange = new TCanvas("AliAnalysisTaskCheckV0","Multiplicity",10,10,510,510);
2676   canPerformanceStrange->Divide(2,1);
2677   if (fHistV0Multiplicity->GetMaximum() > 0.) canPerformanceStrange->cd(1)->SetLogy();
2678   fHistV0Multiplicity->SetMarkerStyle(25);
2679   fHistV0Multiplicity->DrawCopy("E");
2680   if (fHistV0MultiplicityMI->GetMaximum() > 0.) canPerformanceStrange->cd(2)->SetLogy();
2681   fHistV0MultiplicityMI->SetMarkerStyle(24);
2682   fHistV0MultiplicityMI->DrawCopy("E");
2683
2684
2685 */
2686 }
2687
2688 //----------------------------------------------------------------------------
2689
2690 Double_t AliAnalysisTaskPerformanceStrange::MyRapidity(Double_t rE, Double_t rPz) const
2691 {
2692   // Local calculation for rapidity
2693   return 0.5*TMath::Log((rE+rPz)/(rE-rPz+1.e-13));
2694 }
2695 //----------------------------------------------------------------------------
2696