PWGHF/vertexingHF/AliAnalysisTaskSELc2V0bachelorTMVA.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  * appeuear in the supporting documentation. The authors make no claims   *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id: AliAnalysisTaskSELc2V0bachelorTMVA.cxx 62231 2013-04-29 09:47:26Z fprino $ */
  17
  18 //
  19 //
  20 //               Base class for Lc2V0 Analysis to be used with TMVA
  21 //
  22 //
  23
  24 //------------------------------------------------------------------------------------------
  25 //
  26 //  Author: C. Zampolli (from AliAnalysisTaskSELc2V0bachelor class by A.De Caro, P. Pagano)
  27 //
  28 //------------------------------------------------------------------------------------------
  29
  30 #include <TSystem.h>
  31 #include <TParticle.h>
  32 #include <TParticlePDG.h>
  33 #include <TH1F.h>
  34 #include <TH1I.h>
  35 #include <TH2F.h>
  36 #include <TTree.h>
  37 #include "TROOT.h"
  38 #include <TDatabasePDG.h>
  39 #include "TVector3.h"
  40
  41 #include <AliAnalysisDataSlot.h>
  42 #include <AliAnalysisDataContainer.h>
  43 #include "AliStack.h"
  44 #include "AliMCEvent.h"
  45 #include "AliAnalysisManager.h"
  46 #include "AliAODMCHeader.h"
  47 #include "AliAODHandler.h"
  48 #include "AliLog.h"
  49 #include "AliAODVertex.h"
  50 #include "AliAODRecoDecay.h"
  51 #include "AliAODRecoDecayHF.h"
  52 #include "AliAODRecoCascadeHF.h"
  53 #include "AliAnalysisVertexingHF.h"
  54 #include "AliESDtrack.h"
  55 #include "AliAODTrack.h"
  56 #include "AliAODv0.h"
  57 #include "AliAODMCParticle.h"
  58 #include "AliAnalysisTaskSE.h"
  59 #include "AliAnalysisTaskSELc2V0bachelorTMVA.h"
  60 #include "AliNormalizationCounter.h"
  61 #include "AliAODPidHF.h"
  62 #include "AliPIDResponse.h"
  63 #include "AliPIDCombined.h"
  64 #include "AliTOFPIDResponse.h"
  65 #include "AliInputEventHandler.h"
  66 #include "AliAODRecoDecayHF3Prong.h"
  67 #include "AliKFParticle.h"
  68 #include "AliKFVertex.h"
  69 #include "AliExternalTrackParam.h"
  70 #include "AliESDtrack.h"
  71
  72 using std::cout;
  73 using std::endl;
  74
  75 ClassImp(AliAnalysisTaskSELc2V0bachelorTMVA)
  76
  77 //__________________________________________________________________________
  78 AliAnalysisTaskSELc2V0bachelorTMVA::AliAnalysisTaskSELc2V0bachelorTMVA():
  79 AliAnalysisTaskSE(),
  80   fUseMCInfo(kFALSE),
  81   fOutput(0),
  82   fCEvents(0),
  83   fPIDResponse(0),
  84   fPIDCombined(0),
  85   fIsK0sAnalysis(kFALSE),
  86   fCounter(0),
  87   fAnalCuts(0),
  88   fListCuts(0),
  89   fUseOnTheFlyV0(kFALSE),
  90   fIsEventSelected(kFALSE),
  91   fVariablesTreeSgn(0),
  92   fVariablesTreeBkg(0),
  93   fCandidateVariables(),
  94   fIspA(kFALSE),
  95   fHistoEvents(0),
  96   fHistoLc(0),
  97   fHistoLcOnTheFly(0),
  98   fFillOnlySgn(kFALSE),
  99   fHistoLcBeforeCuts(0),
 100   fHistoFiducialAcceptance(0),
 101   fHistoCodesSgn(0),
 102   fHistoCodesBkg(0),
 103   fHistoLcpKpiBeforeCuts(0),
 104   fVtx1(0),
 105
 106   fHistoDistanceLcToPrimVtx(0),
 107   fHistoDistanceV0ToPrimVtx(0),
 108   fHistoDistanceV0ToLc(0),
 109
 110   fHistoDistanceLcToPrimVtxSgn(0),
 111   fHistoDistanceV0ToPrimVtxSgn(0),
 112   fHistoDistanceV0ToLcSgn(0),
 113
 114   fHistoVtxLcResidualToPrimVtx(0),
 115   fHistoVtxV0ResidualToPrimVtx(0),
 116   fHistoVtxV0ResidualToLc(0),
 117
 118   fHistoMassV0All(0),
 119   fHistoDecayLengthV0All(0),
 120   fHistoLifeTimeV0All(0),
 121
 122   fHistoMassV0True(0),
 123   fHistoDecayLengthV0True(0),
 124   fHistoLifeTimeV0True(0),
 125
 126   fHistoMassV0TrueFromAOD(0),
 127
 128   fHistoMassV0TrueK0S(0),
 129   fHistoDecayLengthV0TrueK0S(0),
 130   fHistoLifeTimeV0TrueK0S(0),
 131
 132   fHistoMassV0TrueK0SFromAOD(0),
 133
 134   fHistoMassLcAll(0),
 135   fHistoDecayLengthLcAll(0),
 136   fHistoLifeTimeLcAll(0),
 137
 138   fHistoMassLcTrue(0),
 139   fHistoDecayLengthLcTrue(0),
 140   fHistoLifeTimeLcTrue(0),
 141
 142   fHistoMassLcTrueFromAOD(0),
 143
 144   fHistoMassV0fromLcAll(0),
 145   fHistoDecayLengthV0fromLcAll(0),
 146   fHistoLifeTimeV0fromLcAll(0),
 147
 148   fHistoMassV0fromLcTrue(0),
 149   fHistoDecayLengthV0fromLcTrue(0),
 150   fHistoLifeTimeV0fromLcTrue(0),
 151
 152   fHistoMassLcSgn(0),
 153   fHistoMassLcSgnFromAOD(0),
 154   fHistoDecayLengthLcSgn(0),
 155   fHistoLifeTimeLcSgn(0),
 156
 157   fHistoMassV0fromLcSgn(0),
 158   fHistoDecayLengthV0fromLcSgn(0),
 159   fHistoLifeTimeV0fromLcSgn(0),
 160
 161   fHistoKF(0),
 162   fHistoKFV0(0),
 163   fHistoKFLc(0),
 164
 165   fHistoMassKFV0(0),
 166   fHistoDecayLengthKFV0(0),
 167   fHistoLifeTimeKFV0(0),
 168
 169   fHistoMassKFLc(0),
 170   fHistoDecayLengthKFLc(0),
 171   fHistoLifeTimeKFLc(0),
 172
 173   fHistoArmenterosPodolanskiV0KF(0),
 174   fHistoArmenterosPodolanskiV0KFSgn(0),
 175   fHistoArmenterosPodolanskiV0AOD(0),
 176   fHistoArmenterosPodolanskiV0AODSgn(0),
 177
 178   fOutputKF(0),
 179   fmcLabelLc(-1),
 180   ftopoConstraint(kTRUE),
 181   fCallKFVertexing(kFALSE),
 182   fKeepingOnlyHIJINGBkg(kFALSE),
 183   fUtils(0),
 184   fHistoBackground(0),
 185   fCutKFChi2NDF(999999.),
 186   fCutKFDeviationFromVtx(999999.),
 187   fCutKFDeviationFromVtxV0(0.),
 188   fCurrentEvent(-1),
 189   fBField(0),
 190   fKeepingOnlyPYTHIABkg(kFALSE),
 191   fHistoMCLcK0SpGen(0x0),
 192   fHistoMCLcK0SpGenAcc(0x0),
 193   fHistoMCLcK0SpGenLimAcc(0x0)
 194 {
 195   //
 196   // Default ctor
 197   //
 198 }
 199 //___________________________________________________________________________
 200 AliAnalysisTaskSELc2V0bachelorTMVA::AliAnalysisTaskSELc2V0bachelorTMVA(const Char_t* name,
 201                                                                        AliRDHFCutsLctoV0* analCuts, Bool_t useOnTheFly) :
 202   AliAnalysisTaskSE(name),
 203   fUseMCInfo(kFALSE),
 204   fOutput(0),
 205   fCEvents(0),
 206   fPIDResponse(0),
 207   fPIDCombined(0),
 208   fIsK0sAnalysis(kFALSE),
 209   fCounter(0),
 210   fAnalCuts(analCuts),
 211   fListCuts(0),
 212   fUseOnTheFlyV0(useOnTheFly),
 213   fIsEventSelected(kFALSE),
 214   fVariablesTreeSgn(0),
 215   fVariablesTreeBkg(0),
 216   fCandidateVariables(),
 217   fIspA(kFALSE),
 218   fHistoEvents(0),
 219   fHistoLc(0),
 220   fHistoLcOnTheFly(0),
 221   fFillOnlySgn(kFALSE),
 222   fHistoLcBeforeCuts(0),
 223   fHistoFiducialAcceptance(0),
 224   fHistoCodesSgn(0),
 225   fHistoCodesBkg(0),
 226   fHistoLcpKpiBeforeCuts(0),
 227   fVtx1(0),
 228
 229   fHistoDistanceLcToPrimVtx(0),
 230   fHistoDistanceV0ToPrimVtx(0),
 231   fHistoDistanceV0ToLc(0),
 232
 233   fHistoDistanceLcToPrimVtxSgn(0),
 234   fHistoDistanceV0ToPrimVtxSgn(0),
 235   fHistoDistanceV0ToLcSgn(0),
 236
 237   fHistoVtxLcResidualToPrimVtx(0),
 238   fHistoVtxV0ResidualToPrimVtx(0),
 239   fHistoVtxV0ResidualToLc(0),
 240
 241   fHistoMassV0All(0),
 242   fHistoDecayLengthV0All(0),
 243   fHistoLifeTimeV0All(0),
 244
 245   fHistoMassV0True(0),
 246   fHistoDecayLengthV0True(0),
 247   fHistoLifeTimeV0True(0),
 248
 249   fHistoMassV0TrueFromAOD(0),
 250
 251   fHistoMassV0TrueK0S(0),
 252   fHistoDecayLengthV0TrueK0S(0),
 253   fHistoLifeTimeV0TrueK0S(0),
 254
 255   fHistoMassV0TrueK0SFromAOD(0),
 256
 257   fHistoMassLcAll(0),
 258   fHistoDecayLengthLcAll(0),
 259   fHistoLifeTimeLcAll(0),
 260
 261   fHistoMassLcTrue(0),
 262   fHistoDecayLengthLcTrue(0),
 263   fHistoLifeTimeLcTrue(0),
 264
 265   fHistoMassLcTrueFromAOD(0),
 266
 267   fHistoMassV0fromLcAll(0),
 268   fHistoDecayLengthV0fromLcAll(0),
 269   fHistoLifeTimeV0fromLcAll(0),
 270
 271   fHistoMassV0fromLcTrue(0),
 272   fHistoDecayLengthV0fromLcTrue(0),
 273   fHistoLifeTimeV0fromLcTrue(0),
 274
 275   fHistoMassLcSgn(0),
 276   fHistoMassLcSgnFromAOD(0),
 277   fHistoDecayLengthLcSgn(0),
 278   fHistoLifeTimeLcSgn(0),
 279
 280   fHistoMassV0fromLcSgn(0),
 281   fHistoDecayLengthV0fromLcSgn(0),
 282   fHistoLifeTimeV0fromLcSgn(0),
 283
 284   fHistoKF(0),
 285   fHistoKFV0(0),
 286   fHistoKFLc(0),
 287
 288   fHistoMassKFV0(0),
 289   fHistoDecayLengthKFV0(0),
 290   fHistoLifeTimeKFV0(0),
 291
 292   fHistoMassKFLc(0),
 293   fHistoDecayLengthKFLc(0),
 294   fHistoLifeTimeKFLc(0),
 295
 296   fHistoArmenterosPodolanskiV0KF(0),
 297   fHistoArmenterosPodolanskiV0KFSgn(0),
 298   fHistoArmenterosPodolanskiV0AOD(0),
 299   fHistoArmenterosPodolanskiV0AODSgn(0),
 300
 301   fOutputKF(0),
 302   fmcLabelLc(-1),
 303   ftopoConstraint(kTRUE),
 304   fCallKFVertexing(kFALSE),
 305   fKeepingOnlyHIJINGBkg(kFALSE),
 306   fUtils(0),
 307   fHistoBackground(0),
 308   fCutKFChi2NDF(999999.),
 309   fCutKFDeviationFromVtx(999999.),
 310   fCutKFDeviationFromVtxV0(0.),
 311   fCurrentEvent(-1),
 312   fBField(0),
 313   fKeepingOnlyPYTHIABkg(kFALSE),
 314   fHistoMCLcK0SpGen(0x0),
 315   fHistoMCLcK0SpGenAcc(0x0),
 316   fHistoMCLcK0SpGenLimAcc(0x0)
 317
 318 {
 319   //
 320   // Constructor. Initialization of Inputs and Outputs
 321   //
 322   Info("AliAnalysisTaskSELc2V0bachelorTMVA","Calling Constructor");
 323
 324   DefineOutput(1, TList::Class());  // Tree signal + Tree Bkg + histoEvents
 325   DefineOutput(2, AliNormalizationCounter::Class()); // normalization counter object
 326   DefineOutput(3, TList::Class());  // Cuts
 327   DefineOutput(4, TTree::Class());  // Tree signal + Tree Bkg + histoEvents
 328   DefineOutput(5, TTree::Class());  // Tree signal + Tree Bkg + histoEvents
 329   DefineOutput(6, TList::Class());  // Tree signal + Tree Bkg + histoEvents
 330
 331 }
 332
 333 //___________________________________________________________________________
 334 AliAnalysisTaskSELc2V0bachelorTMVA::~AliAnalysisTaskSELc2V0bachelorTMVA() {
 335   //
 336   // destructor
 337   //
 338   Info("~AliAnalysisTaskSELc2V0bachelorTMVA","Calling Destructor");
 339
 340   if (fOutput) {
 341     delete fOutput;
 342     fOutput = 0;
 343   }
 344
 345   if (fPIDResponse) {
 346     delete  fPIDResponse;
 347   }
 348
 349   if (fPIDCombined) {
 350     delete  fPIDCombined;
 351   }
 352
 353   if (fCounter) {
 354     delete fCounter;
 355     fCounter = 0;
 356   }
 357
 358   if (fAnalCuts) {
 359     delete fAnalCuts;
 360     fAnalCuts = 0;
 361   }
 362
 363   if (fListCuts) {
 364     delete fListCuts;
 365     fListCuts = 0;
 366   }
 367
 368   if(fVariablesTreeSgn){
 369     delete fVariablesTreeSgn;
 370     fVariablesTreeSgn = 0;
 371   }
 372
 373   if(fVariablesTreeBkg){
 374     delete fVariablesTreeBkg;
 375     fVariablesTreeBkg = 0;
 376   }
 377
 378   if (fOutputKF) {
 379     delete fOutputKF;
 380     fOutputKF = 0;
 381   }
 382
 383   if (fUtils) {
 384     delete fUtils;
 385     fUtils = 0;
 386   }
 387
 388 }
 389 //_________________________________________________
 390 void AliAnalysisTaskSELc2V0bachelorTMVA::Init() {
 391   //
 392   // Initialization
 393   //
 394
 395   fIsEventSelected=kFALSE;
 396
 397   if (fDebug > 1) AliInfo("Init");
 398
 399   fListCuts = new TList();
 400   fListCuts->SetOwner();
 401   fListCuts->Add(new AliRDHFCutsLctoV0(*fAnalCuts));
 402   PostData(3,fListCuts);
 403
 404   if (fUseMCInfo && (fKeepingOnlyHIJINGBkg || fKeepingOnlyPYTHIABkg)) fUtils = new AliVertexingHFUtils();
 405
 406   return;
 407 }
 408
 409 //________________________________________ terminate ___________________________
 410 void AliAnalysisTaskSELc2V0bachelorTMVA::Terminate(Option_t*)
 411 {
 412   // The Terminate() function is the last function to be called during
 413   // a query. It always runs on the client, it can be used to present
 414   // the results graphically or save the results to file.
 415
 416   //AliInfo("Terminate","");
 417   AliAnalysisTaskSE::Terminate();
 418
 419   fOutput = dynamic_cast<TList*> (GetOutputData(1));
 420   if (!fOutput) {
 421     AliError("fOutput not available");
 422     return;
 423   }
 424
 425
 426   AliDebug(2, Form("At MC level, %f Lc --> K0S + p were found", fHistoMCLcK0SpGen->GetEntries()));
 427   AliDebug(2, Form("At MC level, %f Lc --> K0S + p were found in the acceptance", fHistoMCLcK0SpGenAcc->GetEntries()));
 428
 429   fOutputKF = dynamic_cast<TList*> (GetOutputData(6));
 430   if (!fOutputKF) {
 431     AliError("fOutputKF not available");
 432     return;
 433   }
 434
 435   return;
 436 }
 437
 438 //___________________________________________________________________________
 439 void AliAnalysisTaskSELc2V0bachelorTMVA::UserCreateOutputObjects() {
 440   // output
 441   AliInfo(Form("CreateOutputObjects of task %s\n", GetName()));
 442
 443   //slot #1
 444   //OpenFile(1);
 445   fOutput = new TList();
 446   fOutput->SetOwner();
 447   fOutput->SetName("listTrees");
 448
 449   // Output slot 1: list of 2 trees (Sgn + Bkg) of the candidate variables
 450   const char* nameoutput = GetOutputSlot(1)->GetContainer()->GetName();
 451   fVariablesTreeSgn = new TTree(Form("%s_Sgn", nameoutput), "Candidates variables tree, Signal");
 452   fVariablesTreeBkg = new TTree(Form("%s_Bkg", nameoutput), "Candidates variables tree, Background");
 453   Int_t nVar = 86;
 454   fCandidateVariables = new Float_t [nVar];
 455   TString * fCandidateVariableNames = new TString[nVar];
 456   fCandidateVariableNames[0]="massLc2K0Sp";
 457   fCandidateVariableNames[1]="massLc2Lambdapi";
 458   fCandidateVariableNames[2]="massK0S";
 459   fCandidateVariableNames[3]="massLambda";
 460   fCandidateVariableNames[4]="massLambdaBar";
 461   fCandidateVariableNames[5]="cosPAK0S";
 462   fCandidateVariableNames[6]="dcaV0";
 463   fCandidateVariableNames[7]="tImpParBach";
 464   fCandidateVariableNames[8]="tImpParV0";
 465   fCandidateVariableNames[9]="nSigmaTPCpr";
 466   fCandidateVariableNames[10]="nSigmaTPCpi";
 467   fCandidateVariableNames[11]="nSigmaTPCka";
 468   fCandidateVariableNames[12]="nSigmaTOFpr";
 469   fCandidateVariableNames[13]="nSigmaTOFpi";
 470   fCandidateVariableNames[14]="nSigmaTOFka";
 471   fCandidateVariableNames[15]="bachelorPt";
 472   fCandidateVariableNames[16]="V0positivePt";
 473   fCandidateVariableNames[17]="V0negativePt";
 474   fCandidateVariableNames[18]="dcaV0pos";
 475   fCandidateVariableNames[19]="dcaV0neg";
 476   fCandidateVariableNames[20]="v0Pt";
 477   fCandidateVariableNames[21]="massGamma";
 478   fCandidateVariableNames[22]="LcPt";
 479   fCandidateVariableNames[23]="combinedProtonProb";
 480   fCandidateVariableNames[24]="LcEta";
 481   fCandidateVariableNames[25]="V0positiveEta";
 482   fCandidateVariableNames[26]="V0negativeEta";
 483   fCandidateVariableNames[27]="TPCProtonProb";
 484   fCandidateVariableNames[28]="TOFProtonProb";
 485   fCandidateVariableNames[29]="bachelorEta";
 486   fCandidateVariableNames[30]="LcP";
 487   fCandidateVariableNames[31]="bachelorP";
 488   fCandidateVariableNames[32]="v0P";
 489   fCandidateVariableNames[33]="V0positiveP";
 490   fCandidateVariableNames[34]="V0negativeP";
 491   fCandidateVariableNames[35]="LcY";
 492   fCandidateVariableNames[36]="v0Y";
 493   fCandidateVariableNames[37]="bachelorY";
 494   fCandidateVariableNames[38]="V0positiveY";
 495   fCandidateVariableNames[39]="V0negativeY";
 496   fCandidateVariableNames[40]="v0Eta";
 497   fCandidateVariableNames[41]="DecayLengthLc";
 498   fCandidateVariableNames[42]="DecayLengthK0S";
 499   fCandidateVariableNames[43]="CtLc";
 500   fCandidateVariableNames[44]="CtK0S";
 501   fCandidateVariableNames[45]="bachCode";
 502   fCandidateVariableNames[46]="k0SCode";
 503
 504   fCandidateVariableNames[47]="V0KFmass";
 505   fCandidateVariableNames[48]="V0KFdecayLength";
 506   fCandidateVariableNames[49]="V0KFlifeTime";
 507
 508   fCandidateVariableNames[50]="V0KFmassErr";
 509   fCandidateVariableNames[51]="V0KFdecayTimeErr";
 510   fCandidateVariableNames[52]="V0KFlifeTimeErr";
 511
 512   fCandidateVariableNames[53]="LcKFmass";
 513   fCandidateVariableNames[54]="LcKFdecayLength";
 514   fCandidateVariableNames[55]="LcKFlifeTime";
 515
 516   fCandidateVariableNames[56]="LcKFmassErr";
 517   fCandidateVariableNames[57]="LcKFdecayTimeErr";
 518   fCandidateVariableNames[58]="LcKFlifeTimeErr";
 519
 520   fCandidateVariableNames[59]="LcKFDistToPrimVtx";
 521   fCandidateVariableNames[60]="V0KFDistToPrimVtx";
 522   fCandidateVariableNames[61]="V0KFDistToLc";
 523   fCandidateVariableNames[62]="alphaArmKF";
 524   fCandidateVariableNames[63]="ptArmKF";
 525   fCandidateVariableNames[64]="alphaArm";
 526   fCandidateVariableNames[65]="ptArm";
 527
 528   fCandidateVariableNames[66]="ITSrefitV0pos";
 529   fCandidateVariableNames[67]="ITSrefitV0neg";
 530
 531   fCandidateVariableNames[68]="TPCClV0pos";
 532   fCandidateVariableNames[69]="TPCClV0neg";
 533
 534   fCandidateVariableNames[70]="v0Xcoord";
 535   fCandidateVariableNames[71]="v0Ycoord";
 536   fCandidateVariableNames[72]="v0Zcoord";
 537   fCandidateVariableNames[73]="primVtxX";
 538   fCandidateVariableNames[74]="primVtxY";
 539   fCandidateVariableNames[75]="primVtxZ";
 540
 541   fCandidateVariableNames[76]="ITSclBach";
 542   fCandidateVariableNames[77]="SPDclBach";
 543
 544   fCandidateVariableNames[78]="ITSclV0pos";
 545   fCandidateVariableNames[79]="SPDclV0pos";
 546   fCandidateVariableNames[80]="ITSclV0neg";
 547   fCandidateVariableNames[81]="SPDclV0neg";
 548
 549   fCandidateVariableNames[82]="alphaArmLc";
 550   fCandidateVariableNames[83]="alphaArmLcCharge";
 551   fCandidateVariableNames[84]="ptArmLc";
 552
 553   fCandidateVariableNames[85]="CosThetaStar";
 554
 555
 556   for(Int_t ivar=0; ivar<nVar; ivar++){
 557     fVariablesTreeSgn->Branch(fCandidateVariableNames[ivar].Data(),&fCandidateVariables[ivar],Form("%s/f",fCandidateVariableNames[ivar].Data()));
 558     fVariablesTreeBkg->Branch(fCandidateVariableNames[ivar].Data(),&fCandidateVariables[ivar],Form("%s/f",fCandidateVariableNames[ivar].Data()));
 559   }
 560
 561   fHistoEvents = new TH1F("fHistoEvents", "fHistoEvents", 2, -0.5, 1.5);
 562   TString labelEv[2] = {"NotSelected", "Selected"};
 563   for (Int_t ibin = 1; ibin <= fHistoEvents->GetNbinsX(); ibin++){
 564     fHistoEvents->GetXaxis()->SetBinLabel(ibin, labelEv[ibin-1].Data());
 565   }
 566
 567   fHistoLc = new TH1F("fHistoLc", "fHistoLc", 2, -0.5, 1.5);
 568
 569   fHistoLcOnTheFly = new TH1F("fHistoLcOnTheFly", "fHistoLcOnTheFly", 4, -0.5, 3.5);
 570   TString labelOnTheFly[4] = {"OnTheFly-Bkg", "OfflineBkg", "OnTheFly-Sgn", "OfflineSgn"};
 571   for (Int_t ibin = 1; ibin <= fHistoLcOnTheFly->GetNbinsX(); ibin++){
 572     fHistoLcOnTheFly->GetXaxis()->SetBinLabel(ibin, labelOnTheFly[ibin-1].Data());
 573   }
 574
 575   fHistoLcBeforeCuts = new TH1F("fHistoLcBeforeCuts", "fHistoLcBeforeCuts", 2, -0.5, 1.5);
 576   TString labelBeforeCuts[2] = {"Bkg", "Sgn"};
 577   for (Int_t ibin = 1; ibin <= fHistoLcBeforeCuts->GetNbinsX(); ibin++){
 578     fHistoLcBeforeCuts->GetXaxis()->SetBinLabel(ibin, labelBeforeCuts[ibin-1].Data());
 579     fHistoLc->GetXaxis()->SetBinLabel(ibin, labelBeforeCuts[ibin-1].Data());
 580   }
 581
 582   fHistoFiducialAcceptance = new TH1F("fHistoFiducialAcceptance", "fHistoFiducialAcceptance", 4, -0.5, 3.5);
 583   TString labelAcc[4] = {"NotOk-Bkg", "Ok-Bkg", "NotOk-Sgn", "Ok-Sgn"};
 584   for (Int_t ibin = 1; ibin <= fHistoFiducialAcceptance->GetNbinsX(); ibin++){
 585     fHistoFiducialAcceptance->GetXaxis()->SetBinLabel(ibin, labelAcc[ibin-1].Data());
 586   }
 587
 588   fHistoCodesSgn = new TH2F("fHistoCodesSgn", "fHistoCodes for Signal; bachelor; K0S", 7, -1.5, 5.5, 9, -1.5, 7.5);
 589   fHistoCodesBkg = new TH2F("fHistoCodesBkg", "fHistoCodes for Background; bachelor; K0S", 7, -1.5, 5.5, 9, -1.5, 7.5);
 590
 591   TString labelx[7] = { "kBachInvalid", "kBachFake", "kBachNoProton", "kBachPrimary", "kBachNoLambdaMother",
 592                         "kBachDifferentLambdaMother",   "kBachCorrectLambdaMother"};
 593   TString labely[9] = { "kK0SInvalid", "kK0SFake", "kK0SNoK0S", "kK0SWithoutMother", "kK0SNotFromK0",
 594                         "kK0Primary", "kK0NoLambdaMother", "kK0DifferentLambdaMother", "kK0CorrectLambdaMother"};
 595
 596   for (Int_t ibin = 1; ibin <= fHistoCodesSgn->GetNbinsX(); ibin++){
 597     fHistoCodesSgn->GetXaxis()->SetBinLabel(ibin, labelx[ibin-1].Data());
 598     fHistoCodesBkg->GetXaxis()->SetBinLabel(ibin, labelx[ibin-1].Data());
 599   }
 600   for (Int_t ibin = 1; ibin <= fHistoCodesSgn->GetNbinsY(); ibin++){
 601     fHistoCodesSgn->GetYaxis()->SetBinLabel(ibin, labely[ibin-1].Data());
 602     fHistoCodesBkg->GetYaxis()->SetBinLabel(ibin, labely[ibin-1].Data());
 603   }
 604
 605   fHistoLcpKpiBeforeCuts = new TH1F("fHistoLcpKpiBeforeCuts", "fHistoLcpKpiBeforeCuts", 2, -0.5, 1.5);
 606   for (Int_t ibin = 1; ibin <= fHistoLcpKpiBeforeCuts->GetNbinsX(); ibin++){
 607     fHistoLcpKpiBeforeCuts->GetXaxis()->SetBinLabel(ibin, labelBeforeCuts[ibin-1].Data());
 608   }
 609
 610   fHistoBackground = new TH1F("fHistoBackground", "fHistoBackground", 4, -0.5, 3.5);
 611   TString labelBkg[4] = {"Injected", "Non-injected", "Non-PYTHIA", "PYTHIA"};
 612   for (Int_t ibin = 1; ibin <= fHistoBackground->GetNbinsX(); ibin++){
 613     fHistoBackground->GetXaxis()->SetBinLabel(ibin, labelBkg[ibin-1].Data());
 614   }
 615
 616   //fOutput->Add(fVariablesTreeSgn);
 617   //fOutput->Add(fVariablesTreeBkg);
 618
 619   const Float_t ptbins[15] = {0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 12., 17., 25., 35.};
 620
 621   fHistoMCLcK0SpGen = new TH1F("fHistoMCLcK0SpGen", "fHistoMCLcK0SpGen", 14, ptbins);
 622   fHistoMCLcK0SpGenAcc = new TH1F("fHistoMCLcK0SpGenAcc", "fHistoMCLcK0SpGenAcc", 14, ptbins);
 623   fHistoMCLcK0SpGenLimAcc = new TH1F("fHistoMCLcK0SpGenLimAcc", "fHistoMCLcK0SpGenLimAcc", 14, ptbins);
 624
 625   fOutput->Add(fHistoEvents);
 626   fOutput->Add(fHistoLc);
 627   fOutput->Add(fHistoLcOnTheFly);
 628   fOutput->Add(fHistoLcBeforeCuts);
 629   fOutput->Add(fHistoFiducialAcceptance);
 630   fOutput->Add(fHistoCodesSgn);
 631   fOutput->Add(fHistoCodesBkg);
 632   fOutput->Add(fHistoLcpKpiBeforeCuts);
 633   fOutput->Add(fHistoBackground);
 634   fOutput->Add(fHistoMCLcK0SpGen);
 635   fOutput->Add(fHistoMCLcK0SpGenAcc);
 636   fOutput->Add(fHistoMCLcK0SpGenLimAcc);
 637
 638   PostData(1, fOutput);
 639   PostData(4, fVariablesTreeSgn);
 640   PostData(5, fVariablesTreeBkg);
 641   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
 642   AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
 643   fPIDResponse = inputHandler->GetPIDResponse();
 644
 645   if (fAnalCuts->GetIsUsePID()){
 646     /*
 647       fAnalCuts->GetPidHF()->SetPidResponse(fPIDResponse);
 648       fAnalCuts->GetPidV0pos()->SetPidResponse(fPIDResponse);
 649       fAnalCuts->GetPidV0neg()->SetPidResponse(fPIDResponse);
 650       fAnalCuts->GetPidHF()->SetOldPid(kFALSE);
 651       fAnalCuts->GetPidV0pos()->SetOldPid(kFALSE);
 652       fAnalCuts->GetPidV0neg()->SetOldPid(kFALSE);
 653     */
 654     fAnalCuts->SetUsePID(kFALSE); // I don't want to use the PID through the cut object, but I will use the PID response directly!!!
 655   }
 656
 657   // Setting properties of PID
 658   fPIDCombined=new AliPIDCombined;
 659   fPIDCombined->SetDefaultTPCPriors();
 660   fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC+AliPIDResponse::kDetTOF);
 661   //fPIDCombined->SetPriorDistribution((AliPID::EParticleType)ispec,fPriors[ispec]);
 662
 663   fCounter = new AliNormalizationCounter("NormalizationCounter");
 664   fCounter->Init();
 665   PostData(2, fCounter);
 666
 667   // Histograms from KF
 668
 669   if (fCallKFVertexing){
 670     fHistoDistanceLcToPrimVtx    = new TH1D("fHistoDistanceLcToPrimVtx", "Lc distance to Prim Vertex from KF; distance [cm]", 1000, 0., 1);
 671     fHistoDistanceV0ToPrimVtx    = new TH1D("fHistoDistanceV0ToPrimVtx", "V0 distance to Prim Vertex from KF; distance [cm]", 1000, 0., 100.);
 672     fHistoDistanceV0ToLc         = new TH1D("fHistoDistanceV0ToLc", "V0 distance to Lc Vertex from KF; distance [cm]", 1000, 0., 100.);
 673
 674     fHistoDistanceLcToPrimVtxSgn = new TH1D("fHistoDistanceLcToPrimVtxSgn", "Lc Sgn distance to Prim Vertex from KF; distance [cm]", 1000, 0., 1);
 675     fHistoDistanceV0ToPrimVtxSgn = new TH1D("fHistoDistanceV0ToPrimVtxSgn", "V0 Sgn distance to Prim Vertex from KF; distance [cm]", 1000, 0., 100.);
 676     fHistoDistanceV0ToLcSgn      = new TH1D("fHistoDistanceV0ToLcSgn", "V0 Sgn distance to Lc Vertex from KF; distance [cm]", 1000, 0., 100.);
 677
 678     fHistoVtxLcResidualToPrimVtx = new TH1D("fHistoVtxLcResidualToPrimVtx", "Residual between MC and KF (MC - KF): Lc to Prim Vtx; distance [cm]", 1000, -5., 5.);
 679     fHistoVtxV0ResidualToPrimVtx = new TH1D("fHistoVtxV0ResidualToPrimVtx", "Residual between MC and KF (MC - KF): V0 to Prim Vtx; distance [cm]", 1000, -5., 5.);
 680     fHistoVtxV0ResidualToLc      = new TH1D("fHistoVtxV0ResidualToLc", "Residual between MC and KF: V0 to Lc (MC - KF); distance [cm]", 1000, -5., 5.);
 681
 682     fHistoMassV0All              = new TH1D("fHistoMassV0All", "V0 Mass; mass", 500, 0.4, 0.6);
 683     fHistoDecayLengthV0All       = new TH1D("fHistoDecayLengthV0All", "V0 Decay Length; decayLength", 500, -10, 10.0);
 684     fHistoLifeTimeV0All          = new TH1D("fHistoLifeTimeV0All", "V0 Life Time; lifeTime", 500, -10.0, 10.0);
 685
 686     fHistoMassV0True             = new TH1D("fHistoMassV0True", "True V0 Mass; mass", 500, 0.4, 0.6);
 687     fHistoDecayLengthV0True      = new TH1D("fHistoDecayLengthV0True", "True V0 Decay Length; decayLength", 500, -10, 10.0);
 688     fHistoLifeTimeV0True         = new TH1D("fHistoLifeTimeV0True", "True V0 Life Time; lifeTime", 500, -10.0, 10.0);
 689
 690     fHistoMassV0TrueFromAOD      = new TH1D("fHistoMassV0TrueFormAOD", "True V0 Mass (AOD); mass", 500, 0.4, 0.6);
 691
 692     fHistoMassV0TrueK0S          = new TH1D("fHistoMassV0TrueK0S", "True V0-K0S Mass; mass", 500, 0.4, 0.6);
 693     fHistoDecayLengthV0TrueK0S   = new TH1D("fHistoDecayLengthV0TrueK0S", "True V0-K0S Decay Length; decayLength", 500, -10, 10.0);
 694     fHistoLifeTimeV0TrueK0S      = new TH1D("fHistoLifeTimeV0TrueK0S", "True V0-K0S Life Time; lifeTime", 500, -10.0, 10.0);
 695
 696     fHistoMassV0TrueK0SFromAOD   = new TH1D("fHistoMassV0TrueK0SFormAOD", "True V0-K0S Mass (AOD); mass", 500, 0.4, 0.6);
 697
 698     fHistoMassLcAll              = new TH1D("fHistoMassLcAll", "Lc Mass; mass", 500, 2.0, 3.0);
 699     fHistoDecayLengthLcAll       = new TH1D("fHistoDecayLengthLcAll", "Lc Decay Length; decayLenght", 100000, -0.1, 0.1);
 700     fHistoLifeTimeLcAll          = new TH1D("fHistoLifeTimeLcAll", "Lc Life Time; lifeTime", 100000, -0.1, 0.1);
 701
 702     fHistoMassLcTrue             = new TH1D("fHistoMassLcTrue", "True Lc Mass; mass", 500, 2.0, 3.0);
 703     fHistoDecayLengthLcTrue      = new TH1D("fHistoDecayLengthLcTrue", "True Lc Decay Length; decayLength", 100000, -0.1, 0.1);
 704     fHistoLifeTimeLcTrue         = new TH1D("fHistoLifeTimeLcTrue", "True Lc Life Time; lifeTime", 100000, -0.1, 0.1);
 705
 706     fHistoMassLcTrueFromAOD      = new TH1D("fHistoMassLcTrueFromAOD", "True Lc Mass (AOD); mass", 500, 2.0, 3.0);
 707
 708     fHistoMassV0fromLcAll        = new TH1D("fHistoMassV0fromLcAll", "V0 mass from Lc built in KF; mass", 500, 0.4, 0.6);
 709     fHistoDecayLengthV0fromLcAll = new TH1D("fHistoDecayLengthV0fromLcAll", "V0 Decay Length from Lc built in KF; decayLength", 500, 0, 10.0);
 710     fHistoLifeTimeV0fromLcAll    = new TH1D("fHistoLifeTimeV0fromLcAll", "V0 Life Time from Lc built in KF; lifeTime", 500, 0.0, 3.0);
 711
 712     fHistoMassV0fromLcTrue       = new TH1D("fHistoMassV0fromLcTrue", "V0 mass from true Lc built in KF; mass", 500, 0.4, 0.6);
 713     fHistoDecayLengthV0fromLcTrue= new TH1D("fHistoDecayLengthV0fromLcTrue", "V0 Decay Length from true Lc built in KF; decayLength", 500, 0, 10.0);
 714     fHistoLifeTimeV0fromLcTrue   = new TH1D("fHistoLifeTimeV0fromLcTrue", "V0 Life Time from true Lc built in KF; lifeTime", 500, 0.0, 3.0);
 715
 716     fHistoMassLcSgn              = new TH1D("fHistoMassLcSgn", "True Lc Signal Mass; mass", 500, 2.0, 3.0);
 717     fHistoMassLcSgnFromAOD       = new TH1D("fHistoMassLcSgnFromAOD", "True Lc Signal Mass (AOD); mass", 500, 2.0, 3.0);
 718     fHistoDecayLengthLcSgn       = new TH1D("fHistoDecayLengthLcSgn", "True Lc Signal Decay Length; decayLength", 100000, -0.1, 0.1);
 719     fHistoLifeTimeLcSgn          = new TH1D("fHistoLifeTimeLcSgn", "True Lc Signal Life Time; lifeTime", 100000, -0.1, 0.1);
 720
 721     fHistoMassV0fromLcSgn        = new TH1D("fHistoMassV0fromLcSgn", "V0 from True Lc Signal Mass; mass", 500, 0.4, 0.6);
 722     fHistoDecayLengthV0fromLcSgn = new TH1D("fHistoDecayLengthV0fromLcSgn", "V0 True Lc Signal Decay Length; decayLength", 500, 0, 10.0);
 723     fHistoLifeTimeV0fromLcSgn    = new TH1D("fHistoLifeTimeV0fromLcSgn", "V0 True Lc Signal Life Time; lifeTime", 500, 0.0, 3.0);
 724
 725     fHistoKF = new TH2D("fHistoKF", "Summary from KF; V0 KF; Lc KF", 16, -0.5, 15.5, 16, -0.5, 15.5);
 726     fHistoKFV0 = new TH1D("fHistoKFV0", "Summary from KF; V0 KF", 16, -0.5, 15.5);
 727     fHistoKFLc = new TH1D("fHistoKFLc", "Summary from KF; V0 KF", 16, -0.5, 15.5);
 728     TString axisLabel[16] = {"AllOk", "M_NotOk", "Sm_NotOk", "Dl_NotOk",
 729                              "Lt_NotOk", "M_Sm_NotOk", "M_Dl_NotOk", "M_Lt_NotOk",
 730                              "Dl_Sm_NotOk", "Dl_Lt_NotOk", "Sm_Lt_NotOk", "M_Sm_Dl_NotOk",
 731                              "M_Sm_Lt_NotOk", "Sm_Dl_Lt_NotOk", "M_Dl_Lt_NotOk", "All_NotOk"};
 732
 733     for (Int_t ibin = 1; ibin <=16; ibin++){
 734       fHistoKF->GetXaxis()->SetBinLabel(ibin, axisLabel[ibin-1].Data());
 735       fHistoKF->GetYaxis()->SetBinLabel(ibin, axisLabel[ibin-1].Data());
 736       fHistoKFV0->GetXaxis()->SetBinLabel(ibin, axisLabel[ibin-1].Data());
 737       fHistoKFLc->GetXaxis()->SetBinLabel(ibin, axisLabel[ibin-1].Data());
 738     }
 739
 740     fHistoMassKFV0 = new TH2D("fHistoMassKFV0", "mass vs sigmaMass for V0; mass; sigmaMass", 500, 0.4, 0.6, 500, 0., 10);
 741     fHistoDecayLengthKFV0 = new TH2D("fHistoDecayLengthKFV0", "decayLength vs sigmaDecayLength for V0; decayLength; sigmaDecayLength", 500, -10, 10, 500, 0., 10);
 742     fHistoLifeTimeKFV0 = new TH2D("fHistoLifeTimeKFV0", "lifeTime vs sigmalifeTime for V0; lifeTime; sigmaLifeTime", 500, -10, 10, 500, 0., 10);
 743
 744     fHistoMassKFLc = new TH2D("fHistoMassKFLc", "mass vs sigmaMass for Lc; mass; sigmaMass", 500, 0.4, 0.6, 500, 0., 10);
 745     fHistoDecayLengthKFLc = new TH2D("fHistoDecayLengthKFLc", "decayLength vs sigmaDecayLength for Lc; decayLength; sigmaDecayLength", 500, -10, 10, 500, 0., 10);
 746     fHistoLifeTimeKFLc = new TH2D("fHistoLifeTimeKFLc", "lifeTime vs sigmalifeTime for Lc; lifeTime; sigmaLifeTime", 500, -10, 10, 500, 0., 10);
 747
 748     fHistoArmenterosPodolanskiV0KF = new TH2D("fHistoArmenterosPodolanskiV0KF", "V0 ArmenterosPodolanski from KF; #alpha; Qt", 1000, -1, 1, 1000, 0, 1);
 749     fHistoArmenterosPodolanskiV0KFSgn = new TH2D("fHistoArmenterosPodolanskiV0KFSgn", "V0 (signal) ArmenterosPodolanski from KF; #alpha; Qt", 1000, -1, 1, 1000, 0, 1);
 750     fHistoArmenterosPodolanskiV0AOD = new TH2D("fHistoArmenterosPodolanskiV0AOD", "V0 ArmenterosPodolanski from AOD; #alpha; Qt", 1000, -1, 1, 1000, 0, 1);
 751     fHistoArmenterosPodolanskiV0AODSgn = new TH2D("fHistoArmenterosPodolanskiV0AODSgn", "V0 (signal) ArmenterosPodolanski from AOD; #alpha; Qt", 1000, -1, 1, 1000, 0, 1);
 752   }
 753
 754   fOutputKF = new TList();
 755   fOutputKF->SetOwner();
 756   fOutputKF->SetName("listHistoKF");
 757
 758   if (fCallKFVertexing){
 759     fOutputKF->Add(fHistoDistanceLcToPrimVtx);
 760     fOutputKF->Add(fHistoDistanceV0ToPrimVtx);
 761     fOutputKF->Add(fHistoDistanceV0ToLc);
 762
 763     fOutputKF->Add(fHistoDistanceLcToPrimVtxSgn);
 764     fOutputKF->Add(fHistoDistanceV0ToPrimVtxSgn);
 765     fOutputKF->Add(fHistoDistanceV0ToLcSgn);
 766
 767     fOutputKF->Add(fHistoVtxLcResidualToPrimVtx);
 768     fOutputKF->Add(fHistoVtxV0ResidualToPrimVtx);
 769     fOutputKF->Add(fHistoVtxV0ResidualToLc);
 770
 771     fOutputKF->Add(fHistoMassV0All);
 772     fOutputKF->Add(fHistoDecayLengthV0All);
 773     fOutputKF->Add(fHistoLifeTimeV0All);
 774
 775     fOutputKF->Add(fHistoMassV0True);
 776     fOutputKF->Add(fHistoDecayLengthV0True);
 777     fOutputKF->Add(fHistoLifeTimeV0True);
 778
 779     fOutputKF->Add(fHistoMassV0TrueFromAOD);
 780
 781     fOutputKF->Add(fHistoMassV0TrueK0S);
 782     fOutputKF->Add(fHistoDecayLengthV0TrueK0S);
 783     fOutputKF->Add(fHistoLifeTimeV0TrueK0S);
 784
 785     fOutputKF->Add(fHistoMassV0TrueK0SFromAOD);
 786
 787     fOutputKF->Add(fHistoMassLcAll);
 788     fOutputKF->Add(fHistoDecayLengthLcAll);
 789     fOutputKF->Add(fHistoLifeTimeLcAll);
 790
 791     fOutputKF->Add(fHistoMassLcTrue);
 792     fOutputKF->Add(fHistoDecayLengthLcTrue);
 793     fOutputKF->Add(fHistoLifeTimeLcTrue);
 794
 795     fOutputKF->Add(fHistoMassLcTrueFromAOD);
 796
 797     fOutputKF->Add(fHistoMassV0fromLcAll);
 798     fOutputKF->Add(fHistoDecayLengthV0fromLcAll);
 799     fOutputKF->Add(fHistoLifeTimeV0fromLcAll);
 800
 801     fOutputKF->Add(fHistoMassV0fromLcTrue);
 802     fOutputKF->Add(fHistoDecayLengthV0fromLcTrue);
 803     fOutputKF->Add(fHistoLifeTimeV0fromLcTrue);
 804
 805     fOutputKF->Add(fHistoMassLcSgn);
 806     fOutputKF->Add(fHistoMassLcSgnFromAOD);
 807     fOutputKF->Add(fHistoDecayLengthLcSgn);
 808     fOutputKF->Add(fHistoLifeTimeLcSgn);
 809
 810     fOutputKF->Add(fHistoMassV0fromLcSgn);
 811     fOutputKF->Add(fHistoDecayLengthV0fromLcSgn);
 812     fOutputKF->Add(fHistoLifeTimeV0fromLcSgn);
 813
 814     fOutputKF->Add(fHistoKF);
 815     fOutputKF->Add(fHistoKFV0);
 816     fOutputKF->Add(fHistoKFLc);
 817
 818     fOutputKF->Add(fHistoMassKFV0);
 819     fOutputKF->Add(fHistoDecayLengthKFV0);
 820     fOutputKF->Add(fHistoLifeTimeKFV0);
 821
 822     fOutputKF->Add(fHistoMassKFLc);
 823     fOutputKF->Add(fHistoDecayLengthKFLc);
 824     fOutputKF->Add(fHistoLifeTimeKFLc);
 825
 826     fOutputKF->Add(fHistoArmenterosPodolanskiV0KF);
 827     fOutputKF->Add(fHistoArmenterosPodolanskiV0KFSgn);
 828     fOutputKF->Add(fHistoArmenterosPodolanskiV0AOD);
 829     fOutputKF->Add(fHistoArmenterosPodolanskiV0AODSgn);
 830   }
 831
 832   PostData(6, fOutputKF);
 833
 834   return;
 835 }
 836
 837 //_________________________________________________
 838 void AliAnalysisTaskSELc2V0bachelorTMVA::UserExec(Option_t *)
 839 {
 840   // user exec
 841   if (!fInputEvent) {
 842     AliError("NO EVENT FOUND!");
 843     return;
 844   }
 845
 846   fCurrentEvent++;
 847   AliDebug(2, Form("Processing event = %d", fCurrentEvent));
 848   AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
 849   TClonesArray *arrayLctopKos=0;
 850
 851   TClonesArray *array3Prong = 0;
 852
 853   if (!aodEvent && AODEvent() && IsStandardAOD()) {
 854     // In case there is an AOD handler writing a standard AOD, use the AOD
 855     // event in memory rather than the input (ESD) event.
 856     aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
 857     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 858     // have to taken from the AOD event hold by the AliAODExtension
 859     AliAODHandler* aodHandler = (AliAODHandler*)
 860       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 861
 862     if (aodHandler->GetExtensions()) {
 863       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 864       AliAODEvent *aodFromExt = ext->GetAOD();
 865       arrayLctopKos=(TClonesArray*)aodFromExt->GetList()->FindObject("CascadesHF");
 866
 867       array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
 868     }
 869   } else {
 870     arrayLctopKos=(TClonesArray*)aodEvent->GetList()->FindObject("CascadesHF");
 871
 872     array3Prong=(TClonesArray*)aodEvent->GetList()->FindObject("Charm3Prong");
 873   }
 874
 875   if ( !fUseMCInfo && fIspA) {
 876     fAnalCuts->SetTriggerClass("");
 877     fAnalCuts->SetTriggerMask(AliVEvent::kINT7);
 878   }
 879
 880   fCounter->StoreEvent(aodEvent,fAnalCuts,fUseMCInfo);
 881
 882   // mc analysis
 883   TClonesArray *mcArray = 0;
 884   AliAODMCHeader *mcHeader=0;
 885
 886   if (fUseMCInfo) {
 887     // MC array need for matching
 888     mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
 889     if (!mcArray) {
 890       AliError("Could not find Monte-Carlo in AOD");
 891       return;
 892     }
 893     // load MC header
 894     mcHeader = (AliAODMCHeader*)aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 895     if (!mcHeader) {
 896       AliError("AliAnalysisTaskSELc2V0bachelorTMVA::UserExec: MC header branch not found!\n");
 897       return;
 898     }
 899     //Printf("Filling MC histo");
 900     FillMCHisto(mcArray);
 901   }
 902
 903   // AOD primary vertex
 904   fVtx1 = (AliAODVertex*)aodEvent->GetPrimaryVertex();
 905   if (!fVtx1) return;
 906   if (fVtx1->GetNContributors()<1) return;
 907
 908   fIsEventSelected = fAnalCuts->IsEventSelected(aodEvent);
 909
 910   if ( !fIsEventSelected ) {
 911     fHistoEvents->Fill(0);
 912     return; // don't take into account not selected events
 913   }
 914   fHistoEvents->Fill(1);
 915
 916   // Setting magnetic field for KF vertexing
 917   fBField = aodEvent->GetMagneticField();
 918   AliKFParticle::SetField(fBField);
 919
 920   Int_t nSelectedAnal = 0;
 921   if (fIsK0sAnalysis) {
 922     MakeAnalysisForLc2prK0S(arrayLctopKos, mcArray,
 923                             nSelectedAnal, fAnalCuts,
 924                             array3Prong, mcHeader);
 925   }
 926   fCounter->StoreCandidates(aodEvent,nSelectedAnal,kFALSE);
 927
 928   PostData(1, fOutput);
 929   PostData(2, fCounter);
 930   PostData(4, fVariablesTreeSgn);
 931   PostData(5, fVariablesTreeBkg);
 932   PostData(6, fOutputKF);
 933
 934 }
 935 //-------------------------------------------------------------------------------
 936 void AliAnalysisTaskSELc2V0bachelorTMVA::FillMCHisto(TClonesArray *mcArray){
 937
 938   // method to fill MC histo: how many Lc --> K0S + p are there at MC level
 939   for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) {
 940     AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart));
 941     if (!mcPart){
 942       AliError("Failed casting particle from MC array!, Skipping particle");
 943       AliInfo("Failed casting particle from MC array!, Skipping particle");
 944       continue;
 945     }
 946     Int_t pdg = mcPart->GetPdgCode();
 947     if (TMath::Abs(pdg) != 4122){
 948       AliDebug(2, Form("MC particle %d is not a Lc: its pdg code is %d", iPart, pdg));
 949       continue;
 950     }
 951     AliDebug(2, Form("Step 0 ok: MC particle %d is a Lc: its pdg code is %d", iPart, pdg));
 952     Int_t labeldaugh0 = mcPart->GetDaughter(0);
 953     Int_t labeldaugh1 = mcPart->GetDaughter(1);
 954     if (labeldaugh0 <= 0 || labeldaugh1 <= 0){
 955       AliDebug(2, Form("The MC particle doesn't have correct daughters, skipping!!"));
 956       continue;
 957     }
 958     else if (labeldaugh1 - labeldaugh0 == 1){
 959       AliDebug(2, Form("Step 1 ok: The MC particle has correct daughters!!"));
 960       AliAODMCParticle* daugh0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labeldaugh0));
 961       AliAODMCParticle* daugh1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labeldaugh1));
 962       if(!daugh0 || !daugh1){
 963         AliDebug(2,"Particle daughters not properly retrieved!");
 964         return;
 965       }
 966       Int_t pdgCodeDaugh0 = TMath::Abs(daugh0->GetPdgCode());
 967       Int_t pdgCodeDaugh1 = TMath::Abs(daugh1->GetPdgCode());
 968       AliAODMCParticle* bachelorMC = daugh0;
 969       AliAODMCParticle* v0MC = daugh1;
 970       AliDebug(2, Form("pdgCodeDaugh0 = %d, pdgCodeDaugh1 = %d", pdgCodeDaugh0, pdgCodeDaugh1));
 971       if ((pdgCodeDaugh0 == 311 && pdgCodeDaugh1 == 2212) || (pdgCodeDaugh0 == 2212 && pdgCodeDaugh1 == 311)){
 972         // we are in the case of Lc --> K0 + p; now we have to check if the K0 decays in K0S, and if this goes in pi+pi-
 973         /// first, we set the bachelor and the v0: above we assumed first proton and second V0, but we could have to change it:
 974         if (pdgCodeDaugh0 == 311 && pdgCodeDaugh1 == 2212) {
 975           bachelorMC = daugh1;
 976           v0MC = daugh0;
 977         }
 978         AliDebug(2, Form("Number of Daughters of v0 = %d", v0MC->GetNDaughters()));
 979         if (v0MC->GetNDaughters() != 1) {
 980           AliDebug(2, "The K0 does not decay in 1 body only! Impossible... Continuing...");
 981           continue;
 982         }
 983         else { // So far: Lc --> K0 + p, K0 with 1 daughter
 984           AliDebug(2, "Step 2 ok: The K0 does decay in 1 body only! ");
 985           Int_t labelK0daugh = v0MC->GetDaughter(0);
 986           AliAODMCParticle* partK0S = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelK0daugh));
 987           if(!partK0S){
 988             AliError("Error while casting particle! returning a NULL array");
 989             AliInfo("Error while casting particle! returning a NULL array");
 990             continue;
 991           }
 992           else { // So far: Lc --> K0 + p, K0 with 1 daughter that we can access
 993             if (partK0S->GetNDaughters() != 2 || TMath::Abs(partK0S->GetPdgCode() != 310)){
 994               AliDebug(2, "The K0 daughter is not a K0S or does not decay in 2 bodies");
 995               continue;
 996             }
 997             else { // So far: Lc --> K0 + p, K0 --> K0S, K0S in 2 bodies
 998               AliDebug(2, "Step 3 ok: The K0 daughter is a K0S and does decay in 2 bodies");
 999               Int_t labelK0Sdaugh0 = partK0S->GetDaughter(0);
1000               Int_t labelK0Sdaugh1 = partK0S->GetDaughter(1);
1001               AliAODMCParticle* daughK0S0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelK0Sdaugh0));
1002               AliAODMCParticle* daughK0S1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelK0Sdaugh1));
1003               if (!daughK0S0 || ! daughK0S1){
1004                 AliDebug(2, "Could not access K0S daughters, continuing...");
1005                 continue;
1006               }
1007               else { // So far: Lc --> K0 + p, K0 --> K0S, K0S in 2 bodies that we can access
1008                 AliDebug(2, "Step 4 ok: Could access K0S daughters, continuing...");
1009                 Int_t pdgK0Sdaugh0 = daughK0S0->GetPdgCode();
1010                 Int_t pdgK0Sdaugh1 = daughK0S1->GetPdgCode();
1011                 if (TMath::Abs(pdgK0Sdaugh0) != 211 || TMath::Abs(pdgK0Sdaugh1) != 211){
1012                   AliDebug(2, "The K0S does not decay in pi+pi-, continuing");
1013                   //AliInfo("The K0S does not decay in pi+pi-, continuing");
1014                 }
1015                 else { // Full chain: Lc --> K0 + p, K0 --> K0S, K0S --> pi+pi-
1016                   if (fAnalCuts->IsInFiducialAcceptance(mcPart->Pt(), mcPart->Y())) {
1017                     AliDebug(2, Form("----> Filling histo with pt = %f", mcPart->Pt()));
1018                     if(TMath::Abs(mcPart->Y()) < 0.5) fHistoMCLcK0SpGenLimAcc->Fill(mcPart->Pt());
1019                     fHistoMCLcK0SpGen->Fill(mcPart->Pt());
1020                     if(!(TMath::Abs(bachelorMC->Eta()) > 0.9 || bachelorMC->Pt() < 0.1 ||
1021                          TMath::Abs(daughK0S0->Eta()) > 0.9 || daughK0S0->Pt() < 0.1 ||
1022                          TMath::Abs(daughK0S1->Eta()) > 0.9 || daughK0S1->Pt() < 0.1)) {
1023                       fHistoMCLcK0SpGenAcc->Fill(mcPart->Pt());
1024                     }
1025                   }
1026                   else {
1027                     AliDebug(2, "not in fiducial acceptance! Skipping");
1028                     continue;
1029                   }
1030                 }
1031               }
1032             }
1033           }
1034         }
1035       }
1036     }
1037   } // closing loop over mcArray
1038
1039   return;
1040
1041 }
1042
1043 //-------------------------------------------------------------------------------
1044 void AliAnalysisTaskSELc2V0bachelorTMVA::MakeAnalysisForLc2prK0S(TClonesArray *arrayLctopKos,
1045                                                                  TClonesArray *mcArray,
1046                                                                  Int_t &nSelectedAnal,
1047                                                                  AliRDHFCutsLctoV0 *cutsAnal, TClonesArray *array3Prong,
1048                                                                  AliAODMCHeader* aodheader){
1049   //Lc prong needed to MatchToMC method
1050
1051   Int_t pdgCand = 4122;
1052   Int_t pdgDgLctoV0bachelor[2]={2212, 310};
1053   Int_t pdgDgV0toDaughters[2]={211, 211};
1054
1055   Int_t pdgDgLctopKpi[3]={2212, 321, 211};
1056
1057   // loop to search for candidates Lc->p+K+pi
1058   Int_t n3Prong = array3Prong->GetEntriesFast();
1059   Int_t nCascades= arrayLctopKos->GetEntriesFast();
1060
1061   //AliInfo(Form("\n\n\n\n3 prong candidates = %d, ncascades = %d \n\n\n\n\n", n3Prong, nCascades));
1062   for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
1063     AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
1064     //Filling a control histogram with no cuts
1065     if (fUseMCInfo) {
1066
1067       // find associated MC particle for Lc -> p+K+pi
1068       Int_t mcLabel = d->MatchToMC(4122, mcArray, 3, pdgDgLctopKpi);
1069       //Int_t mcLabelTemp = d->MatchToMC(4122, mcArray);
1070       //Printf("mcLabel = %d, mcLabelTemp = %d", mcLabel, mcLabelTemp);
1071       if (mcLabel >= 0) {
1072
1073         AliAODMCParticle *partLcpKpi = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcLabel));
1074         if(partLcpKpi){
1075           Int_t pdgCode = partLcpKpi->GetPdgCode();
1076           AliDebug(2,Form(" ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ MClabel=%d ~~~~~~~~~~ pdgCode=%d", mcLabel, pdgCode));
1077           //AliInfo(Form(" ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ MClabel=%d ~~~~~~~~~~ pdgCode=%d", mcLabel, pdgCode));
1078           fHistoLcpKpiBeforeCuts->Fill(1);
1079
1080         }
1081       }
1082       else {
1083         //AliInfo(Form(" ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ MClabel=%d ~~~~~~~~~~", mcLabel));
1084         fHistoLcpKpiBeforeCuts->Fill(0);
1085       }
1086     }
1087   }
1088
1089   // loop over cascades to search for candidates Lc->p+K0S
1090
1091   for (Int_t iLctopK0s = 0; iLctopK0s < nCascades; iLctopK0s++) {
1092
1093     // Lc candidates and K0s from Lc
1094     AliAODRecoCascadeHF* lcK0spr = dynamic_cast<AliAODRecoCascadeHF*>(arrayLctopKos->At(iLctopK0s));
1095     if (!lcK0spr) {
1096       AliDebug(2,Form("Cascade %d doens't exist, skipping",iLctopK0s));
1097       continue;
1098     }
1099
1100     //Filling a control histogram with no cuts
1101     if (fUseMCInfo) {
1102
1103       Int_t pdgCode=-2;
1104
1105       // find associated MC particle for Lc -> p+K0 and K0S->pi+pi
1106       fmcLabelLc = lcK0spr->MatchToMC(pdgCand, pdgDgLctoV0bachelor[1], pdgDgLctoV0bachelor, pdgDgV0toDaughters, mcArray, kTRUE);
1107       if (fmcLabelLc>=0) {
1108         AliDebug(2, Form("----> cascade number %d (total cascade number = %d) is a Lc!", iLctopK0s, nCascades));
1109
1110         AliAODMCParticle *partLc = dynamic_cast<AliAODMCParticle*>(mcArray->At(fmcLabelLc));
1111         if(partLc){
1112           pdgCode = partLc->GetPdgCode();
1113           if (pdgCode<0) AliDebug(2,Form(" ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ MClabel=%d ~~~~~~~~~~ pdgCode=%d", fmcLabelLc, pdgCode));
1114           pdgCode = TMath::Abs(pdgCode);
1115           fHistoLcBeforeCuts->Fill(1);
1116
1117         }
1118       }
1119       else {
1120         fHistoLcBeforeCuts->Fill(0);
1121         pdgCode=-1;
1122       }
1123     }
1124
1125     //if (!lcK0spr->GetSecondaryVtx()) {
1126     //  AliInfo("No secondary vertex");
1127     //continue;
1128     //}
1129
1130     if (lcK0spr->GetNDaughters()!=2) {
1131       AliDebug(2,Form("Cascade %d has not 2 daughters (nDaughters=%d)",iLctopK0s,lcK0spr->GetNDaughters()));
1132       continue;
1133     }
1134
1135     AliAODv0 * v0part = dynamic_cast<AliAODv0*>(lcK0spr->Getv0());
1136     AliAODTrack * bachPart = dynamic_cast<AliAODTrack*>(lcK0spr->GetBachelor());
1137     if (!v0part || !bachPart) {
1138       AliDebug(2,Form("Cascade %d has no V0 or no bachelor object",iLctopK0s));
1139       continue;
1140     }
1141
1142
1143     if (!v0part->GetSecondaryVtx()) {
1144       AliDebug(2,Form("No secondary vertex for V0 by cascade %d",iLctopK0s));
1145       continue;
1146     }
1147
1148     if (v0part->GetNDaughters()!=2) {
1149       AliDebug(2,Form("current V0 has not 2 daughters (onTheFly=%d, nDaughters=%d)",v0part->GetOnFlyStatus(),v0part->GetNDaughters()));
1150       continue;
1151     }
1152
1153     AliAODTrack * v0Pos = dynamic_cast<AliAODTrack*>(lcK0spr->Getv0PositiveTrack());
1154     AliAODTrack * v0Neg = dynamic_cast<AliAODTrack*>(lcK0spr->Getv0NegativeTrack());
1155     if (!v0Neg || !v0Pos) {
1156       AliDebug(2,Form("V0 by cascade %d has no V0positive of V0negative object",iLctopK0s));
1157       continue;
1158     }
1159
1160
1161     if (v0Pos->Charge() == v0Neg->Charge()) {
1162       AliDebug(2,Form("V0 by cascade %d has daughters with the same sign: IMPOSSIBLE!",iLctopK0s));
1163       continue;
1164     }
1165
1166     Int_t isLc = 0;
1167
1168     if (fUseMCInfo) {
1169
1170       Int_t pdgCode=-2;
1171
1172       // find associated MC particle for Lc -> p+K0 and K0S->pi+pi
1173       Int_t mcLabel = lcK0spr->MatchToMC(pdgCand, pdgDgLctoV0bachelor[1], pdgDgLctoV0bachelor, pdgDgV0toDaughters, mcArray, kTRUE);
1174       if (mcLabel>=0) {
1175         AliDebug(2,Form(" ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~cascade number %d (total cascade number = %d)", iLctopK0s, nCascades));
1176
1177         AliAODMCParticle *partLc = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcLabel));
1178         if(partLc){
1179           pdgCode = partLc->GetPdgCode();
1180           if (pdgCode<0) AliDebug(2,Form(" ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ MClabel=%d ~~~~~~~~~~ pdgCode=%d", mcLabel, pdgCode));
1181           pdgCode = TMath::Abs(pdgCode);
1182           isLc = 1;
1183           fHistoLc->Fill(1);
1184
1185         }
1186       }
1187       else {
1188         fHistoLc->Fill(0);
1189         pdgCode=-1;
1190       }
1191     }
1192     AliDebug(2, Form("\n\n\n Analysing candidate %d\n", iLctopK0s));
1193     AliDebug(2, Form(">>>>>>>>>> Candidate is background, fFillOnlySgn = %d --> SKIPPING", fFillOnlySgn));
1194     if (!isLc) {
1195       if (fFillOnlySgn) { // if it is background, and we want only signal, we do not fill the tree
1196         continue;
1197       }
1198       else { // checking if we want to fill the background
1199         if (fKeepingOnlyHIJINGBkg){
1200           // we have decided to fill the background only when the candidate has the daugthers that all come from HIJING underlying event!
1201           Bool_t isCandidateInjected = fUtils->HasCascadeCandidateAnyDaughInjected(lcK0spr, aodheader, mcArray);
1202           if (!isCandidateInjected){
1203             AliDebug(2, "The candidate is from HIJING (i.e. not injected), keeping it to fill background");
1204             fHistoBackground->Fill(1);
1205           }
1206           else {
1207             AliDebug(2, "The candidate is NOT from HIJING, we skip it when filling background");
1208             fHistoBackground->Fill(0);
1209             continue;
1210           }
1211         }
1212         else if (fKeepingOnlyPYTHIABkg){
1213           // we have decided to fill the background only when the candidate has the daugthers that all come from HIJING underlying event!
1214           AliAODTrack *bachelor = (AliAODTrack*)lcK0spr->GetBachelor();
1215           AliAODTrack *v0pos = (AliAODTrack*)lcK0spr->Getv0PositiveTrack();
1216           AliAODTrack *v0neg = (AliAODTrack*)lcK0spr->Getv0NegativeTrack();
1217           if (!bachelor || !v0pos || !v0neg) {
1218             AliDebug(2, "Cannot retrieve one of the tracks while checking origin, continuing");
1219             continue;
1220           }
1221           else {
1222             Int_t labelbachelor = TMath::Abs(bachelor->GetLabel());
1223             Int_t labelv0pos = TMath::Abs(v0pos->GetLabel());
1224             Int_t labelv0neg = TMath::Abs(v0neg->GetLabel());
1225             AliAODMCParticle* MCbachelor =  (AliAODMCParticle*)mcArray->At(labelbachelor);
1226             AliAODMCParticle* MCv0pos =  (AliAODMCParticle*)mcArray->At(labelv0pos);
1227             AliAODMCParticle* MCv0neg =  (AliAODMCParticle*)mcArray->At(labelv0neg);
1228             if (!MCbachelor || !MCv0pos || !MCv0neg) {
1229               AliDebug(2, "Cannot retrieve MC particle for one of the tracks while checking origin, continuing");
1230               continue;
1231             }
1232             else {
1233               Int_t isBachelorFromPythia = fUtils->CheckOrigin(mcArray, MCbachelor, kTRUE);
1234               Int_t isv0posFromPythia = fUtils->CheckOrigin(mcArray, MCv0pos, kTRUE);
1235               Int_t isv0negFromPythia = fUtils->CheckOrigin(mcArray, MCv0neg, kTRUE);
1236               if (isBachelorFromPythia != 0 && isv0posFromPythia != 0 && isv0negFromPythia != 0){
1237                 AliDebug(2, "The candidate is from PYTHIA (i.e. all daughters originate from a quark), keeping it to fill background");
1238                 fHistoBackground->Fill(2);
1239               }
1240               else {
1241                 AliDebug(2, "The candidate is NOT from PYTHIA, we skip it when filling background");
1242                 fHistoBackground->Fill(3);
1243                 continue;
1244               }
1245             }
1246           }
1247         }
1248       }
1249     }
1250
1251     FillLc2pK0Sspectrum(lcK0spr, isLc, nSelectedAnal, cutsAnal, mcArray);
1252
1253   }
1254
1255   return;
1256
1257 }
1258 //________________________________________________________________________
1259 void AliAnalysisTaskSELc2V0bachelorTMVA::FillLc2pK0Sspectrum(AliAODRecoCascadeHF *part,
1260                                                              Int_t isLc,
1261                                                              Int_t &nSelectedAnal,
1262                                                              AliRDHFCutsLctoV0 *cutsAnal,
1263                                                              TClonesArray *mcArray){
1264   //
1265   // Fill histos for Lc -> K0S+proton
1266   //
1267
1268   /*
1269     if (!part->GetOwnPrimaryVtx()) {
1270     //Printf("No primary vertex for Lc found!!");
1271     part->SetOwnPrimaryVtx(fVtx1);
1272     }
1273     else {
1274     //Printf("Yu-huuuu!!! primary vertex for Lc found!!");
1275     }
1276   */
1277   Double_t invmassLc = part->InvMassLctoK0sP();
1278
1279   AliAODv0 * v0part = part->Getv0();
1280   Bool_t onFlyV0 = v0part->GetOnFlyStatus(); // on-the-flight V0s
1281   if (onFlyV0){ // on-the-fly V0
1282     if (isLc) { // Lc
1283       fHistoLcOnTheFly->Fill(2.);
1284     }
1285     else { // not Lc
1286       fHistoLcOnTheFly->Fill(0.);
1287     }
1288   }
1289   else { // offline V0
1290     if (isLc) { // Lc
1291       fHistoLcOnTheFly->Fill(3.);
1292     }
1293     else { // not Lc
1294       fHistoLcOnTheFly->Fill(1.);
1295     }
1296   }
1297
1298   Double_t dcaV0 = v0part->GetDCA();
1299   Double_t invmassK0s = v0part->MassK0Short();
1300
1301   if ( (cutsAnal->IsInFiducialAcceptance(part->Pt(),part->Y(4122))) ) {
1302     if (isLc) {
1303       fHistoFiducialAcceptance->Fill(3.);
1304     }
1305     else {
1306       fHistoFiducialAcceptance->Fill(1.);
1307     }
1308   }
1309   else {
1310     if (isLc) {
1311       fHistoFiducialAcceptance->Fill(2.);
1312     }
1313     else {
1314       fHistoFiducialAcceptance->Fill(0.);
1315     }
1316   }
1317
1318   Int_t isInV0window = (((cutsAnal->IsSelectedSingleCut(part, AliRDHFCuts::kCandidate, 2)) & (AliRDHFCutsLctoV0::kLcToK0Spr)) == (AliRDHFCutsLctoV0::kLcToK0Spr)); // cut on V0 invMass
1319
1320   if (isInV0window == 0) {
1321     AliDebug(2, "No: The candidate has NOT passed the mass cuts!");
1322     if (isLc) Printf("SIGNAL candidate rejected");
1323     return;
1324   }
1325   else AliDebug(2, "Yes: The candidate has passed the mass cuts!");
1326
1327   Bool_t isInCascadeWindow = (((cutsAnal->IsSelectedSingleCut(part, AliRDHFCuts::kCandidate, 0)) & (AliRDHFCutsLctoV0::kLcToK0Spr)) == (AliRDHFCutsLctoV0::kLcToK0Spr)); // cut on Lc->p+K0S invMass
1328
1329   if (!isInCascadeWindow) {
1330     AliDebug(2, "No: The candidate has NOT passed the cascade window cuts!");
1331     if (isLc) Printf("SIGNAL candidate rejected");
1332     return;
1333   }
1334   else AliDebug(2, "Yes: The candidate has passed the cascade window cuts!");
1335
1336   Bool_t isCandidateSelectedCuts = (((cutsAnal->IsSelected(part, AliRDHFCuts::kCandidate)) & (AliRDHFCutsLctoV0::kLcToK0Spr)) == (AliRDHFCutsLctoV0::kLcToK0Spr)); // kinematic/topological cuts
1337   if (!isCandidateSelectedCuts){
1338     AliDebug(2, "No: Analysis cuts kCandidate level NOT passed");
1339     if (isLc) Printf("SIGNAL candidate rejected");
1340     return;
1341   }
1342   else {
1343     AliDebug(2, "Yes: Analysis cuts kCandidate level passed");
1344   }
1345
1346   AliAODTrack *bachelor = (AliAODTrack*)part->GetBachelor();
1347   if (!bachelor) {
1348     AliDebug(2, Form("Very weird, the bachelor is not there... returning for this candidate"));
1349     return;
1350   }
1351
1352   //Bool_t isBachelorID = (((cutsAnal->IsSelected(part,AliRDHFCuts::kPID))&(AliRDHFCutsLctoV0::kLcToK0Spr))==(AliRDHFCutsLctoV0::kLcToK0Spr)); // ID x bachelor
1353   Double_t probTPCTOF[AliPID::kSPECIES]={-1.};
1354
1355   UInt_t detUsed = fPIDCombined->ComputeProbabilities(bachelor, fPIDResponse, probTPCTOF);
1356   AliDebug(2, Form("detUsed (TPCTOF case) = %d", detUsed));
1357   Double_t probProton = -1.;
1358   //  Double_t probPion = -1.;
1359   //  Double_t probKaon = -1.;
1360   if (detUsed == (UInt_t)fPIDCombined->GetDetectorMask() ) {
1361     AliDebug(2, Form("We have found the detector mask for TOF + TPC: probProton will be set to %f", probTPCTOF[AliPID::kProton]));
1362     probProton = probTPCTOF[AliPID::kProton];
1363     // probPion = probTPCTOF[AliPID::kPion];
1364     // probKaon = probTPCTOF[AliPID::kKaon];
1365   }
1366   else { // if you don't have both TOF and TPC, try only TPC
1367     fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
1368     AliDebug(2, "We did not find the detector mask for TOF + TPC, let's see only TPC");
1369     detUsed = fPIDCombined->ComputeProbabilities(bachelor, fPIDResponse, probTPCTOF);
1370     AliDebug(2,Form(" detUsed (TPC case) = %d", detUsed));
1371     if (detUsed == (UInt_t)fPIDCombined->GetDetectorMask()) {
1372       probProton = probTPCTOF[AliPID::kProton];
1373       // probPion = probTPCTOF[AliPID::kPion];
1374       // probKaon = probTPCTOF[AliPID::kKaon];
1375       AliDebug(2, Form("TPC only worked: probProton will be set to %f", probTPCTOF[AliPID::kProton]));
1376     }
1377     else {
1378       AliDebug(2, "Only TPC did not work...");
1379     }
1380     // resetting mask to ask for both TPC+TOF
1381     fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC+AliPIDResponse::kDetTOF);
1382   }
1383   AliDebug(2, Form("probProton = %f", probProton));
1384
1385   // now we get the TPC and TOF single PID probabilities (only for Proton, or the tree will explode :) )
1386   Double_t probProtonTPC = -1.;
1387   Double_t probProtonTOF = -1.;
1388   Double_t pidTPC[AliPID::kSPECIES]={-1.};
1389   Double_t pidTOF[AliPID::kSPECIES]={-1.};
1390   Int_t respTPC = fPIDResponse->ComputePIDProbability(AliPIDResponse::kDetTPC, bachelor, AliPID::kSPECIES, pidTPC);
1391   Int_t respTOF = fPIDResponse->ComputePIDProbability(AliPIDResponse::kDetTOF, bachelor, AliPID::kSPECIES, pidTOF);
1392   if (respTPC == AliPIDResponse::kDetPidOk) probProtonTPC = pidTPC[AliPID::kProton];
1393   if (respTOF == AliPIDResponse::kDetPidOk) probProtonTOF = pidTOF[AliPID::kProton];
1394
1395   // checking V0 status (on-the-fly vs offline)
1396   if ( !( !onFlyV0 || (onFlyV0 && fUseOnTheFlyV0) ) ) {
1397     AliDebug(2, "On-the-fly discarded");
1398     return;
1399   }
1400
1401   if ( (((cutsAnal->IsSelected(part,AliRDHFCuts::kAll))&(AliRDHFCutsLctoV0::kLcToK0Spr))==(AliRDHFCutsLctoV0::kLcToK0Spr)) ) {
1402     nSelectedAnal++;
1403   }
1404
1405   if ( !(cutsAnal->IsInFiducialAcceptance(part->Pt(),part->Y(4122))) ) return;
1406
1407   if ( !( ( (cutsAnal->IsSelected(part, AliRDHFCuts::kTracks)) & (AliRDHFCutsLctoV0::kLcToK0Spr)) == (AliRDHFCutsLctoV0::kLcToK0Spr) ) ) { // esd track cuts
1408     if (isLc) Printf("SIGNAL candidate rejected");
1409     AliDebug(2, "No: Analysis cuts kTracks level NOT passed");
1410     return;
1411   }
1412   else {
1413     AliDebug(2, "Yes: Analysis cuts kTracks level passed");
1414   }
1415
1416   Int_t pdgCand = 4122;
1417   Int_t pdgDgLctoV0bachelor[2]={211, 3122}; // case of wrong decay! Lc --> L + pi
1418   Int_t pdgDgV0toDaughters[2]={2212, 211}; // case of wrong decay! Lc --> L + pi
1419   Int_t isLc2LBarpi=0, isLc2Lpi=0;
1420   Int_t currentLabel = part->GetLabel();
1421   Int_t mcLabel = 0;
1422   if (fUseMCInfo) {
1423     mcLabel = part->MatchToMC(pdgCand, pdgDgLctoV0bachelor[1], pdgDgLctoV0bachelor, pdgDgV0toDaughters, mcArray, kTRUE);
1424     if (mcLabel>=0) {
1425       if (bachelor->Charge()==-1) isLc2LBarpi=1;
1426       if (bachelor->Charge()==+1) isLc2Lpi=1;
1427     }
1428   }
1429
1430   Int_t pdgDg2prong[2] = {211, 211};
1431   Int_t labelK0S = 0;
1432   Int_t isK0S = 0;
1433   if (fUseMCInfo) {
1434     labelK0S = v0part->MatchToMC(310, mcArray, 2, pdgDg2prong);
1435     if (labelK0S>=0) isK0S = 1;
1436   }
1437
1438   pdgDg2prong[0] = 211;
1439   pdgDg2prong[1] = 2212;
1440   Int_t isLambda = 0;
1441   Int_t isLambdaBar = 0;
1442   Int_t lambdaLabel = 0;
1443   if (fUseMCInfo) {
1444     lambdaLabel = v0part->MatchToMC(3122, mcArray, 2, pdgDg2prong);
1445     if (lambdaLabel>=0) {
1446       AliAODMCParticle *lambdaTrack = (AliAODMCParticle*)mcArray->At(lambdaLabel);
1447       if (lambdaTrack->GetPdgCode()==3122) isLambda = 1;
1448       else if (lambdaTrack->GetPdgCode()==-3122) isLambdaBar = 1;
1449     }
1450   }
1451
1452   pdgDg2prong[0] = 11;
1453   pdgDg2prong[1] = 11;
1454   Int_t isGamma = 0;
1455   Int_t gammaLabel = 0;
1456   if (fUseMCInfo) {
1457     gammaLabel = v0part->MatchToMC(22, mcArray, 2, pdgDg2prong);
1458     if (gammaLabel>=0) {
1459       AliAODMCParticle *gammaTrack = (AliAODMCParticle*)mcArray->At(gammaLabel);
1460       if (gammaTrack->GetPdgCode()==22) isGamma = 1;
1461     }
1462   }
1463
1464   Int_t pdgTemp = -1;
1465   if (currentLabel != -1){
1466     AliAODMCParticle *tempPart = (AliAODMCParticle*)mcArray->At(currentLabel);
1467     pdgTemp = tempPart->GetPdgCode();
1468   }
1469   if (isLc) AliDebug(2, Form("Signal: Candidate is a Lc in K0s+p"));
1470   else if (isLc2LBarpi) AliDebug(2, Form("Background: Candidate is a Lc in Lbar + pi"));
1471   else if (isLc2Lpi) AliDebug(2, Form("Background: Candidate is a Lc in L + pi"));
1472   else AliDebug(2, Form("Pure bkg: Candidate has pdg = %d", pdgTemp));
1473   if (isK0S) AliDebug(2, Form("V0 is a K0S"));
1474   else if (isLambda)  AliDebug(2, Form("V0 is a Lambda"));
1475   else if (isLambdaBar)  AliDebug(2, Form("V0 is a LambdaBar"));
1476   else if (isGamma)  AliDebug(2, Form("V0 is a Gamma"));
1477   //else AliDebug(2, Form("V0 is something else!!"));
1478
1479   Double_t invmassLc2Lpi = part->InvMassLctoLambdaPi();
1480   Double_t invmassLambda = v0part->MassLambda();
1481   Double_t invmassLambdaBar = v0part->MassAntiLambda();
1482
1483   //Double_t nSigmaITSpr=-999.;
1484   Double_t nSigmaTPCpr=-999.;
1485   Double_t nSigmaTOFpr=-999.;
1486
1487   //Double_t nSigmaITSpi=-999.;
1488   Double_t nSigmaTPCpi=-999.;
1489   Double_t nSigmaTOFpi=-999.;
1490
1491   //Double_t nSigmaITSka=-999.;
1492   Double_t nSigmaTPCka=-999.;
1493   Double_t nSigmaTOFka=-999.;
1494
1495   /*
1496     cutsAnal->GetPidHF()->GetnSigmaITS(bachelor,4,nSigmaITSpr);
1497     cutsAnal->GetPidHF()->GetnSigmaTPC(bachelor,4,nSigmaTPCpr);
1498     cutsAnal->GetPidHF()->GetnSigmaTOF(bachelor,4,nSigmaTOFpr);
1499     cutsAnal->GetPidHF()->GetnSigmaITS(bachelor,2,nSigmaITSpi);
1500     cutsAnal->GetPidHF()->GetnSigmaTPC(bachelor,2,nSigmaTPCpi);
1501     cutsAnal->GetPidHF()->GetnSigmaTOF(bachelor,2,nSigmaTOFpi);
1502     cutsAnal->GetPidHF()->GetnSigmaITS(bachelor,3,nSigmaITSka);
1503     cutsAnal->GetPidHF()->GetnSigmaTPC(bachelor,3,nSigmaTPCka);
1504     cutsAnal->GetPidHF()->GetnSigmaTOF(bachelor,3,nSigmaTOFka);
1505   */
1506
1507   nSigmaTPCpi = fPIDResponse->NumberOfSigmasTPC(bachelor,(AliPID::kPion));
1508   nSigmaTPCka = fPIDResponse->NumberOfSigmasTPC(bachelor,(AliPID::kKaon));
1509   nSigmaTPCpr = fPIDResponse->NumberOfSigmasTPC(bachelor,(AliPID::kProton));
1510   nSigmaTOFpi = fPIDResponse->NumberOfSigmasTOF(bachelor,(AliPID::kPion));
1511   nSigmaTOFka = fPIDResponse->NumberOfSigmasTOF(bachelor,(AliPID::kKaon));
1512   nSigmaTOFpr = fPIDResponse->NumberOfSigmasTOF(bachelor,(AliPID::kProton));
1513
1514
1515   // Fill candidate variable Tree (track selection, V0 invMass selection)
1516   if (!onFlyV0 && isInV0window && isInCascadeWindow && part->CosV0PointingAngle()>0.99 && TMath::Abs(nSigmaTPCpr) <= 3 && v0part->Getd0Prong(0) < 20 && v0part->Getd0Prong(1) < 20) {
1517
1518     fCandidateVariables[0] = invmassLc;
1519     fCandidateVariables[1] = invmassLc2Lpi;
1520     fCandidateVariables[2] = invmassK0s;
1521     fCandidateVariables[3] = invmassLambda;
1522     fCandidateVariables[4] = invmassLambdaBar;
1523     fCandidateVariables[5] = part->CosV0PointingAngle();
1524     fCandidateVariables[6] = dcaV0;
1525     fCandidateVariables[7] = part->Getd0Prong(0);
1526     fCandidateVariables[8] = part->Getd0Prong(1);
1527     fCandidateVariables[9] = nSigmaTPCpr;
1528     fCandidateVariables[10] = nSigmaTPCpi;
1529     fCandidateVariables[11] = nSigmaTPCka;
1530     fCandidateVariables[12] = nSigmaTOFpr;
1531     fCandidateVariables[13] = nSigmaTOFpi;
1532     fCandidateVariables[14] = nSigmaTOFka;
1533     fCandidateVariables[15] = bachelor->Pt();
1534     AliAODTrack *v0pos = (AliAODTrack*)part->Getv0PositiveTrack();
1535     fCandidateVariables[16] = v0pos->Pt();
1536     AliAODTrack *v0neg = (AliAODTrack*)part->Getv0NegativeTrack();
1537     fCandidateVariables[17] = v0neg->Pt();
1538     fCandidateVariables[18] = v0part->Getd0Prong(0);
1539     fCandidateVariables[19] = v0part->Getd0Prong(1);
1540     fCandidateVariables[20] = v0part->Pt();
1541     fCandidateVariables[21] = v0part->InvMass2Prongs(0,1,11,11);
1542     fCandidateVariables[22] = part->Pt();
1543     fCandidateVariables[23] = probProton;
1544     fCandidateVariables[24] = part->Eta();
1545     fCandidateVariables[25] = v0pos->Eta();
1546     fCandidateVariables[26] = v0neg->Eta();
1547     fCandidateVariables[27] = probProtonTPC;
1548     fCandidateVariables[28] = probProtonTOF;
1549     fCandidateVariables[29] = bachelor->Eta();
1550
1551     fCandidateVariables[30] = part->P();
1552     fCandidateVariables[31] = bachelor->P();
1553     fCandidateVariables[32] = v0part->P();
1554     fCandidateVariables[33] = v0pos->P();
1555     fCandidateVariables[34] = v0neg->P();
1556
1557     fCandidateVariables[35] = part->Y(4122);
1558     fCandidateVariables[36] = bachelor->Y(2212);
1559     fCandidateVariables[37] = v0part->Y(310);
1560     fCandidateVariables[38] = v0pos->Y(211);
1561     fCandidateVariables[39] = v0neg->Y(211);
1562
1563     fCandidateVariables[40] = v0part->Eta();
1564
1565     fCandidateVariables[41] = part->DecayLength();
1566     fCandidateVariables[42] = part->DecayLengthV0();
1567     fCandidateVariables[43] = part->Ct(4122);
1568     fCandidateVariables[44] = v0part->Ct(310, v0part->GetSecondaryVtx());
1569
1570     EBachelor bachCode = kBachInvalid;
1571     EK0S k0SCode = kK0SInvalid;
1572     if (fUseMCInfo) {
1573       bachCode = CheckBachelor(part, bachelor, mcArray);
1574       k0SCode = CheckK0S(part, v0part, mcArray);
1575     }
1576
1577     fCandidateVariables[45] = bachCode;
1578     fCandidateVariables[46] = k0SCode;
1579
1580     Double_t V0KF[3] = {-999999, -999999, -999999};
1581     Double_t errV0KF[3] = {-999999, -999999, -999999};
1582     Double_t LcKF[3] = {-999999, -999999, -999999};
1583     Double_t errLcKF[3] = {-999999, -999999, -999999};
1584     Double_t distances[3] = {-999999, -999999, -999999};
1585     Double_t armPolKF[2] = {-999999, -999999};
1586
1587     if (fCallKFVertexing){
1588       Int_t kfResult = CallKFVertexing(part, v0part, bachelor, mcArray, &V0KF[0], &errV0KF[0], &LcKF[0], &errLcKF[0], &distances[0], &armPolKF[0]);
1589       AliDebug(2, Form("Result from KF = %d", kfResult));
1590     }
1591
1592     /*
1593       for (Int_t i = 0; i< 3; i++){
1594       Printf("i = %d, V0KF = %f, errV0KF = %f, LcKF = %f, errLcKF = %f", V0KF[i], errV0KF[i], LcKF[i], errLcKF[i]);
1595       }
1596     */
1597
1598     fCandidateVariables[47] = V0KF[0];
1599     fCandidateVariables[48] = V0KF[1];
1600     fCandidateVariables[49] = V0KF[2];
1601
1602     fCandidateVariables[50] = errV0KF[0];
1603     fCandidateVariables[51] = errV0KF[1];
1604     fCandidateVariables[52] = errV0KF[2];
1605
1606     fCandidateVariables[53] = LcKF[0];
1607     fCandidateVariables[54] = LcKF[1];
1608     fCandidateVariables[55] = LcKF[2];
1609
1610     fCandidateVariables[56] = errLcKF[0];
1611     fCandidateVariables[57] = errLcKF[1];
1612     fCandidateVariables[58] = errLcKF[2];
1613
1614     fCandidateVariables[59] = distances[0];
1615     fCandidateVariables[60] = distances[1];
1616     fCandidateVariables[61] = distances[2];
1617     fCandidateVariables[62] = armPolKF[0];
1618     fCandidateVariables[63] = armPolKF[1];
1619     fCandidateVariables[64] = v0part->AlphaV0();
1620     fCandidateVariables[65] = v0part->PtArmV0();
1621
1622     AliDebug(2, Form("v0pos->GetStatus() & AliESDtrack::kITSrefit= %d, v0neg->GetStatus() & AliESDtrack::kITSrefit = %d, v0pos->GetTPCClusterInfo(2, 1)= %f, v0neg->GetTPCClusterInfo(2, 1) = %f", (Int_t)(v0pos->GetStatus() & AliESDtrack::kITSrefit), (Int_t)(v0pos->GetStatus() & AliESDtrack::kITSrefit), v0pos->GetTPCClusterInfo(2, 1), v0neg->GetTPCClusterInfo(2, 1)));
1623     fCandidateVariables[66] = v0pos->GetStatus() & AliESDtrack::kITSrefit;
1624     fCandidateVariables[67] = v0neg->GetStatus() & AliESDtrack::kITSrefit;
1625     fCandidateVariables[68] = v0pos->GetTPCClusterInfo(2, 1);
1626     fCandidateVariables[69] = v0neg->GetTPCClusterInfo(2, 1);
1627
1628     fCandidateVariables[70] = v0part->Xv();
1629     fCandidateVariables[71] = v0part->Yv();
1630     fCandidateVariables[72] = v0part->Zv();
1631
1632     fCandidateVariables[73] = fVtx1->GetX();
1633     fCandidateVariables[74] = fVtx1->GetY();
1634     fCandidateVariables[75] = fVtx1->GetZ();
1635
1636     fCandidateVariables[76] = bachelor->GetITSNcls();
1637     fCandidateVariables[77] = bachelor->HasPointOnITSLayer(0) + bachelor->HasPointOnITSLayer(1);
1638
1639     fCandidateVariables[78] = v0pos->GetITSNcls();
1640     fCandidateVariables[79] = v0pos->HasPointOnITSLayer(0) + v0pos->HasPointOnITSLayer(1);
1641
1642     fCandidateVariables[80] = v0neg->GetITSNcls();
1643     fCandidateVariables[81] = v0neg->HasPointOnITSLayer(0) + v0neg->HasPointOnITSLayer(1);
1644
1645     TVector3 mom1(bachelor->Px(), bachelor->Py(), bachelor->Pz());
1646     TVector3 mom2(v0part->Px(), v0part->Py(), v0part->Pz());
1647     TVector3 momTot(part->Px(), part->Py(), part->Pz());
1648
1649     Double_t Ql1 = mom1.Dot(momTot)/momTot.Mag();
1650     Double_t Ql2 = mom2.Dot(momTot)/momTot.Mag();
1651
1652     Double_t alphaArmLc = (Ql1 - Ql2)/(Ql1 + Ql2);
1653     Double_t alphaArmLcCharge = ( bachelor->Charge() > 0 ? (Ql1 - Ql2)/(Ql1 + Ql2) : (Ql2 - Ql1)/(Ql1 + Ql2) );
1654     Double_t ptArmLc = mom1.Perp(momTot);
1655
1656     fCandidateVariables[82] = alphaArmLc;
1657     fCandidateVariables[83] = alphaArmLcCharge;
1658     fCandidateVariables[84] = ptArmLc;
1659
1660     Double_t massK0SPDG = TDatabasePDG::Instance()->GetParticle(310)->Mass();    // mass K0S PDG
1661     Double_t massPrPDG = TDatabasePDG::Instance()->GetParticle(2212)->Mass();    // mass Proton PDG
1662     Double_t massLcPDG = TDatabasePDG::Instance()->GetParticle(4122)->Mass();    // mass Lc PDG
1663
1664     Double_t pStar = TMath::Sqrt((massLcPDG*massLcPDG-massPrPDG*massPrPDG-massK0SPDG*massK0SPDG)*(massLcPDG*massLcPDG-massPrPDG*massPrPDG-massK0SPDG*massK0SPDG)-4.*massPrPDG*massPrPDG*massK0SPDG*massK0SPDG)/(2.*massLcPDG);
1665     Double_t e = part->E(4122);
1666     Double_t beta = part->P()/e;
1667     Double_t gamma = e/massLcPDG;
1668
1669     Double_t cts = (Ql1/gamma-beta*TMath::Sqrt(pStar*pStar+massPrPDG*massPrPDG))/pStar;
1670
1671     fCandidateVariables[85] = cts;
1672
1673     if (fUseMCInfo) {
1674       if (isLc){
1675         AliDebug(2, "Filling Sgn");
1676         fVariablesTreeSgn->Fill();
1677         fHistoCodesSgn->Fill(bachCode, k0SCode);
1678       }
1679       else {
1680         if (fFillOnlySgn == kFALSE){
1681           AliDebug(2, "Filling Bkg");
1682           fVariablesTreeBkg->Fill();
1683           fHistoCodesBkg->Fill(bachCode, k0SCode);
1684         }
1685       }
1686     }
1687     else {
1688       fVariablesTreeSgn->Fill();
1689     }
1690   }
1691
1692   return;
1693
1694 }
1695
1696 //________________________________________________________________________
1697 Int_t AliAnalysisTaskSELc2V0bachelorTMVA::CallKFVertexing(AliAODRecoCascadeHF *cascade, AliAODv0* v0part, AliAODTrack* bach, TClonesArray *mcArray,
1698                                                           Double_t* V0KF, Double_t* errV0KF, Double_t* LcKF, Double_t* errLcKF,
1699                                                           Double_t* distances, Double_t* armPolKF) {
1700
1701   //
1702   // method to perform KF vertexing
1703   // elements: [0] = mass, [1] = DecayLength, [2] = lifeTime
1704   //
1705
1706   Int_t codeKFV0 = -1, codeKFLc = -1;
1707
1708   AliKFVertex primVtxCopy;
1709   Int_t nt = 0, ntcheck = 0;
1710   Double_t pos[3] = {0., 0., 0.};
1711
1712   fVtx1->GetXYZ(pos);
1713   Int_t contr = fVtx1->GetNContributors();
1714   Double_t covmatrix[6] = {0.};
1715   fVtx1->GetCovarianceMatrix(covmatrix);
1716   Double_t chi2 = fVtx1->GetChi2();
1717   AliESDVertex primaryESDVtxCopy(pos, covmatrix, chi2, contr, "Vertex");
1718
1719   // topological constraint
1720   primVtxCopy = AliKFVertex(primaryESDVtxCopy);
1721   nt = primaryESDVtxCopy.GetNContributors();
1722   ntcheck = nt;
1723
1724   Int_t pdg[2] = {211, -211};
1725   Int_t pdgLc[2] = {2212, 310};
1726
1727   Int_t pdgDgV0toDaughters[2] = {211, 211};
1728
1729   Int_t mcLabelV0 = v0part->MatchToMC(310, mcArray, 2, pdgDgV0toDaughters);
1730
1731   // the KF vertex for the V0 has to be built with the prongs of the V0!
1732   Bool_t isMCokV0 = kTRUE, isBkgV0 = kFALSE;
1733   AliKFParticle  V0, positiveV0KF, negativeV0KF;
1734   Int_t labelsv0daugh[2] = {-1, -1};
1735   Int_t idv0daugh[2] = {-1, -1};
1736   AliExternalTrackParam* esdv0Daugh1 = 0x0;
1737   AliExternalTrackParam* esdv0Daugh2 = 0x0;
1738   for(Int_t ipr= 0; ipr < 2; ipr++){ // 0 is positive, 1 is negative
1739     AliAODTrack *aodTrack = (AliAODTrack*)v0part->GetDaughter(ipr);
1740     if(!aodTrack) {
1741       AliDebug(2, "No V0 daughters available");
1742       return -1;
1743     }
1744     Double_t xyz[3], pxpypz[3], cv[21];
1745     Short_t sign;
1746     aodTrack->GetXYZ(xyz);
1747     aodTrack->PxPyPz(pxpypz);
1748     aodTrack->GetCovarianceXYZPxPyPz(cv);
1749     sign = aodTrack->Charge();
1750     AliExternalTrackParam tmp1( xyz, pxpypz, cv, sign);
1751
1752     if (ipr == 0) esdv0Daugh1 = new AliExternalTrackParam( xyz, pxpypz, cv, sign);
1753     else esdv0Daugh2 = new AliExternalTrackParam( xyz, pxpypz, cv, sign);
1754     labelsv0daugh[ipr] = TMath::Abs(aodTrack->GetLabel());
1755     idv0daugh[ipr] = aodTrack->GetID();
1756     if (labelsv0daugh[ipr] == -1) isBkgV0 = kTRUE;
1757
1758     //Printf("v0 daughter %d has label %d", ipr, labelsv0daugh[ipr]);
1759
1760     AliKFParticle daughterKF(*aodTrack, pdg[ipr]); // we assume that the PDG is correct
1761     if (aodTrack->Charge() > 0) { // assigning positive and negative track to KF V0 for Armenteros-Podolanski plot
1762       positiveV0KF = daughterKF;
1763     }
1764     else {
1765       negativeV0KF = daughterKF;
1766     }
1767   }
1768
1769   Double_t xn=0., xp=0.;//, dca;
1770   AliDebug(2, Form("bField = %f, esdv0Daugh1 = %p, esdv0Daugh2 = %p", fBField, esdv0Daugh1, esdv0Daugh2));
1771   //  dca = esdv0Daugh1->GetDCA(esdv0Daugh2, fBField, xn, xp);
1772
1773   AliExternalTrackParam tr1(*esdv0Daugh1);
1774   AliExternalTrackParam tr2(*esdv0Daugh2);
1775   tr1.PropagateTo(xn, fBField);
1776   tr2.PropagateTo(xp, fBField);
1777
1778   AliKFParticle daughterKF1(tr1, 211);
1779   AliKFParticle daughterKF2(tr2, 211);
1780   V0.AddDaughter(positiveV0KF);
1781   V0.AddDaughter(negativeV0KF);
1782   //V0.AddDaughter(daughterKF1);
1783   //V0.AddDaughter(daughterKF2);
1784
1785   delete esdv0Daugh1;
1786   delete esdv0Daugh2;
1787   esdv0Daugh1=0;
1788   esdv0Daugh2=0;
1789   // Checking the quality of the KF V0 vertex
1790   if( V0.GetNDF() < 1 ) {
1791     //Printf("Number of degrees of freedom < 1, continuing");
1792     return -1;
1793   }
1794   if( TMath::Sqrt(TMath::Abs(V0.GetChi2()/V0.GetNDF())) > fCutKFChi2NDF ) {
1795     //Printf("Chi2 per DOF too big, continuing");
1796     return -1;
1797   }
1798
1799   if(ftopoConstraint && nt > 0){
1800     for(Int_t ipr = 0; ipr < 2; ipr++){ // 0 is positive, 1 is negative
1801       AliAODTrack *aodTrack = (AliAODTrack*)v0part->GetDaughter(ipr);
1802       //* subtruct daughters from primary vertex
1803       if(!aodTrack->GetUsedForPrimVtxFit()) {
1804         //Printf("Track %d was not used for primary vertex, continuing", i);
1805         continue;
1806       }
1807       AliKFParticle daughterKF(*aodTrack, pdg[ipr]); // we assume that the PDG is correct
1808       primVtxCopy -= daughterKF;
1809       ntcheck--;
1810     }
1811   }
1812
1813   // Check V0 Chi^2 deviation from primary vertex // not needed for V0 for Lc decay!!
1814   /*
1815   if( V0.GetDeviationFromVertex( primVtxCopy ) < fCutKFDeviationFromVtxV0) {
1816   //Printf("Deviation from vertex too big, continuing");
1817     return -1;
1818   }
1819   */
1820
1821   //* Get V0 invariant mass
1822   Double_t massV0 = 999999, sigmaMassV0 = 999999;
1823   Int_t retMV0 = V0.GetMass( massV0, sigmaMassV0 );
1824   if( retMV0 ) {
1825     if (massV0 < 0) {
1826       codeKFV0 = 1; // Mass not ok
1827       if (sigmaMassV0 > 1e19) codeKFV0 = 5;  // Mass and SigmaMass not ok
1828     }
1829     else if (sigmaMassV0 > 1e19) codeKFV0 = 2; // SigmaMass not ok
1830   }
1831   fHistoMassKFV0->Fill(massV0, sigmaMassV0);
1832
1833   if (massV0 < 0.4) Printf("\n\n>>>>>>>>>> Found the Funny V0 (mass = %f, sigma = %f, AOD mass = %f): labels of the tracks = %d, %d, id = %d and %d", massV0, sigmaMassV0, v0part->MassK0Short(), labelsv0daugh[0], labelsv0daugh[1], idv0daugh[0], idv0daugh[1]);
1834   if (massV0 > 0.55) Printf("\n\n>>>>>>>>>> Found the Funny V0 (mass = %f, , sigma = %f, AOD mass = %f): labels of the tracks = %d, %d, id = %d and %d", massV0, sigmaMassV0, v0part->MassK0Short(), labelsv0daugh[0], labelsv0daugh[1], idv0daugh[0], idv0daugh[1]);
1835
1836   Printf("Vertices: KF:  x = %f, y = %f, z = %f", V0.GetX(), V0.GetY(), V0.GetZ());
1837   Printf("Vertices: AOD: x = %f, y = %f, z = %f", v0part->Xv(), v0part->Yv(), v0part->Zv());
1838
1839   //Printf("Got MC vtx for V0");
1840   if (fUseMCInfo && TMath::Abs(labelsv0daugh[0] - labelsv0daugh[1]) == 1) {
1841     AliAODMCParticle* tmpdaughv01 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelsv0daugh[0]));
1842     AliAODMCParticle* tmpdaughv02 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelsv0daugh[1]));
1843     if (!tmpdaughv01 && labelsv0daugh[0] > 0){
1844       AliDebug(2, "Could not access MC info for first daughter of V0, continuing");
1845     }
1846     if (!tmpdaughv02 && labelsv0daugh[1] > 0){
1847       AliDebug(2, "Could not access MC info for second daughter of V0, continuing");
1848     }
1849     if(tmpdaughv01){
1850       Double_t xPionMC = tmpdaughv01->Xv(); //Production vertex of Pion --> Where K0S decays
1851       Double_t yPionMC = tmpdaughv01->Yv();
1852       Double_t zPionMC = tmpdaughv01->Zv();
1853       //Printf("Got MC vtx for Pion");
1854       Printf("Vertices: MC:  x = %f, y = %f, z = %f", xPionMC, yPionMC, zPionMC);
1855     }
1856   }
1857   else {
1858     Printf("Not a true V0");
1859   }
1860   //massV0=-1;//return -1;// !!!!
1861
1862   // now use what we just try with the bachelor, to build the Lc
1863
1864   // topological constraint
1865   nt = primVtxCopy.GetNContributors();
1866   ntcheck = nt;
1867
1868   Bool_t isMCokLc = kTRUE, isBkgLc = kFALSE;
1869   AliKFParticle  Lc;
1870   Int_t labelsLcdaugh[2] = {-1, -1};
1871   labelsLcdaugh[0] = TMath::Abs(bach->GetLabel());
1872   labelsLcdaugh[1] = mcLabelV0;
1873
1874   if (bach->Charge() < 0) pdgLc[0] = -pdgLc[0];
1875   AliKFParticle daughterKFLc(*bach, pdgLc[0]);
1876   Lc.AddDaughter(daughterKFLc);
1877   TParticlePDG* particlePDG = TDatabasePDG::Instance()->GetParticle(310);
1878   Double_t massPDGK0S = particlePDG->Mass();
1879   V0.SetMassConstraint(massPDGK0S);
1880   Lc.AddDaughter(V0);
1881   if( Lc.GetNDF() < 1 ) {
1882     AliDebug(2, Form("Lc: Number of degrees of freedom < 1 (%d), continuing", Lc.GetNDF()));
1883     return -1;
1884   }
1885   if( TMath::Sqrt(TMath::Abs(Lc.GetChi2()/Lc.GetNDF())) > fCutKFChi2NDF) {
1886     AliDebug(2, Form("Lc: Chi2 per DOF too big, continuing (%f)", TMath::Sqrt(TMath::Abs(Lc.GetChi2()/Lc.GetNDF()))));
1887     return -1;
1888   }
1889
1890   if(ftopoConstraint && nt > 0){
1891     //* subtruct daughters from primary vertex
1892     if(!bach->GetUsedForPrimVtxFit()) {
1893       AliDebug(3, "Lc: Bachelor was not used for primary vertex, not subtracting it from primary vertex");
1894     }
1895     else{
1896       primVtxCopy -= daughterKFLc;
1897       ntcheck--;
1898     }
1899     /* the V0 was added above, so it is ok to remove it without checking
1900        if(!V0->GetUsedForPrimVtxFit()) {
1901        Printf("Lc: V0 was not used for primary vertex, continuing");
1902        continue;
1903        }
1904     */
1905     //primVtxCopy -= V0;
1906     //ntcheck--;
1907   }
1908
1909   // Check Lc Chi^2 deviation from primary vertex
1910   /*
1911   if( Lc.GetDeviationFromVertex( primVtxCopy ) > fCutKFDeviationFromVtx) {
1912     AliDebug(2, Form("Lc: Deviation from vertex too big, continuing (%f)", Lc.GetDeviationFromVertex( primVtxCopy )));
1913     return -1;
1914   }
1915
1916   if(ftopoConstraint){
1917     if(ntcheck>0) {
1918       // Add Lc to primary vertex to improve the primary vertex resolution
1919       primVtxCopy += Lc;
1920       Lc.SetProductionVertex(primVtxCopy);
1921     }
1922   }
1923   */
1924   //* Check chi^2
1925   if( TMath::Sqrt( TMath::Abs(Lc.GetChi2()/Lc.GetNDF())) > fCutKFChi2NDF) {
1926     AliDebug(2, Form("Lc: Final Chi2 per DOF too big, continuing (%f)", TMath::Sqrt( TMath::Abs(Lc.GetChi2()/Lc.GetNDF()))));
1927     return -1;
1928   }
1929
1930   if(ftopoConstraint){
1931     V0.SetProductionVertex(Lc);
1932   }
1933
1934   // After setting the vertex of the V0, getting/filling some info
1935
1936   //* Get V0 decayLength
1937   Double_t decayLengthV0 = 999999, sigmaDecayLengthV0 = 999999;
1938   Int_t retDLV0 = V0.GetDecayLength( decayLengthV0, sigmaDecayLengthV0 );
1939   if( retDLV0 ) {
1940     if (sigmaDecayLengthV0 > 1e19) {
1941       codeKFV0 = 3; // DecayLength not ok
1942       if (massV0 < 0) {
1943         codeKFV0 = 6; // DecayLength and Mass not ok
1944         if (sigmaMassV0 > 1e19) codeKFV0 = 11;  // DecayLength and Mass and SigmaMass not ok
1945       }
1946       else if (sigmaMassV0 > 1e19) codeKFV0 = 8;  // DecayLength and SigmaMass not ok
1947     }
1948   }
1949   fHistoDecayLengthKFV0->Fill(decayLengthV0, sigmaDecayLengthV0);
1950
1951   //* Get V0 life time
1952   Double_t lifeTimeV0 = 999999, sigmaLifeTimeV0 = 999999;
1953   Int_t retTLV0 = V0.GetLifeTime( lifeTimeV0, sigmaLifeTimeV0 );
1954   if( retTLV0 ) {
1955     if (sigmaLifeTimeV0 > 1e19) {
1956       codeKFV0 = 4;  // LifeTime not ok
1957       if (sigmaDecayLengthV0 > 1e19) {
1958         codeKFV0 = 9;  // LifeTime and DecayLength not ok
1959         if (massV0 < 0) {
1960           codeKFV0 = 14; // LifeTime and DecayLength and Mass not ok
1961           if (sigmaMassV0 > 1e19) codeKFV0 = 15; // LifeTime and DecayLength and Mass and SigmaMass not ok
1962         }
1963         else if (sigmaMassV0 > 1e19) codeKFV0 = 13; // LifeTime and DecayLength and SigmaMass not ok
1964       }
1965       else if (massV0 < 0) { // LifeTime and Mass and SigmaMass not ok
1966         codeKFV0 = 7; // LifeTime and Mass not ok
1967         if (sigmaMassV0 > 1e19) codeKFV0 = 12; // LifeTime and Mass and SigmaMass not ok
1968       }
1969       else if (sigmaMassV0 > 1e19) codeKFV0 = 10;  // LifeTime and SigmaMass not ok
1970     }
1971   }
1972   fHistoLifeTimeKFV0->Fill(lifeTimeV0, sigmaLifeTimeV0);
1973
1974   if (codeKFV0 == -1) codeKFV0 = 0;
1975   fHistoKFV0->Fill(codeKFV0);
1976
1977   AliDebug(2, Form("V0: mass = %f, decay length = %f, life time = %f", massV0, decayLengthV0, lifeTimeV0 ));
1978
1979   fHistoMassV0All->Fill(massV0);
1980   fHistoDecayLengthV0All->Fill(decayLengthV0);
1981   fHistoLifeTimeV0All->Fill(lifeTimeV0);
1982
1983   Double_t qtAlphaV0[2] = {0., 0.};
1984   Double_t vtxV0KF[3] = {V0.GetX(), V0.GetY(), V0.GetZ()};
1985   positiveV0KF.TransportToPoint(vtxV0KF);
1986   negativeV0KF.TransportToPoint(vtxV0KF);
1987   V0.GetArmenterosPodolanski(positiveV0KF, negativeV0KF, qtAlphaV0);
1988   AliDebug(2, Form("Armenteros-Podolanski variables: alpha = %f, qt = %f", qtAlphaV0[1], qtAlphaV0[0]));
1989   fHistoArmenterosPodolanskiV0KF->Fill(qtAlphaV0[1], qtAlphaV0[0]);
1990   fHistoArmenterosPodolanskiV0AOD->Fill(v0part->AlphaV0(), v0part->PtArmV0());
1991   armPolKF[0] = qtAlphaV0[1];
1992   armPolKF[1] = qtAlphaV0[0];
1993
1994   // Checking MC info for V0
1995
1996   AliAODMCParticle *motherV0 = 0x0;
1997   AliAODMCParticle *daughv01 = 0x0;
1998   AliAODMCParticle *daughv02 = 0x0;
1999
2000   if (fUseMCInfo) {
2001     daughv01 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelsv0daugh[0]));
2002     daughv02 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelsv0daugh[1]));
2003     if (!daughv01 && labelsv0daugh[0] > 0){
2004       AliDebug(2, "Could not access MC info for first daughter of V0, continuing");
2005       isMCokV0 = kFALSE;
2006     }
2007     if (!daughv02 && labelsv0daugh[1] > 0){
2008       AliDebug(2, "Could not access MC info for second daughter of V0, continuing");
2009       isMCokV0 = kFALSE;
2010     }
2011     if (isMCokV0){
2012       if( daughv01->GetMother() ==  daughv02->GetMother() && daughv01->GetMother()>=0 ){
2013         AliDebug(3, Form("The mother has label %d", daughv01->GetMother()));
2014         motherV0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughv01->GetMother()));
2015         if( motherV0 && TMath::Abs(motherV0->GetPdgCode()) != 21 ){ // These are all V0 that are truly V0, not only K0S, but no gluons
2016           if( motherV0->GetNDaughters() == 2 ){
2017             fHistoMassV0True->Fill(massV0);
2018             fHistoDecayLengthV0True->Fill(decayLengthV0);
2019             fHistoLifeTimeV0True->Fill(lifeTimeV0);
2020             fHistoMassV0TrueFromAOD->Fill(v0part->MassK0Short());
2021             if (TMath::Abs(motherV0->GetPdgCode()) == 310){ // These are true V0 that are also K0S
2022               fHistoMassV0TrueK0S->Fill(massV0);
2023               fHistoDecayLengthV0TrueK0S->Fill(decayLengthV0);
2024               fHistoLifeTimeV0TrueK0S->Fill(lifeTimeV0);
2025               fHistoMassV0TrueK0SFromAOD->Fill(v0part->MassK0Short());
2026             }
2027           }
2028           AliDebug(2, Form("PDG V0 = %d, pi = %d, pj = %d, ndaughters = %d, mc mass = %f, reco mass = %f, v0 mass = %f", motherV0->GetPdgCode(), daughv01->GetPdgCode(), daughv02->GetPdgCode(), motherV0->GetNDaughters(), motherV0->GetCalcMass(), massV0, v0part->MassK0Short()));
2029         }
2030         else if (!motherV0){
2031           AliDebug(3, "could not access MC info for mother, continuing");
2032         }
2033         else {
2034           AliDebug(3, "MC mother is a gluon, continuing");
2035         }
2036       }
2037       else {
2038         AliDebug(3, "Background V0!");
2039         isBkgV0 = kTRUE;
2040       }
2041     }
2042   }
2043
2044   AliDebug(2, Form("isMCokV0 = %d, isBkgV0 = %d", (Int_t)isMCokV0, (Int_t)isBkgV0));
2045
2046   // Going back to Lc
2047
2048   //* Get Lc invariant mass
2049   Double_t massLc = 999999, sigmaMassLc= 999999;
2050   Int_t retMLc = Lc.GetMass( massLc, sigmaMassLc );
2051   if( retMLc ) {
2052     AliDebug(3, Form("----> Could not get mass (%e), and sigma(%e) for Lc, continuing", massLc, sigmaMassLc));
2053     if (massLc < 0) {
2054       codeKFLc = 1; // Mass not ok
2055       if (sigmaMassLc > 1e19) codeKFLc = 5;  // Mass and SigmaMass not ok
2056     }
2057     else if (sigmaMassLc > 1e19) codeKFLc = 2; // SigmaMass not ok
2058   }
2059   fHistoMassKFLc->Fill(massLc, sigmaMassLc);
2060
2061   //* Get Lc Decay length
2062   Double_t decayLengthLc = 999999, sigmaDecayLengthLc = 999999;
2063   Int_t retDLLc = Lc.GetDecayLength( decayLengthLc, sigmaDecayLengthLc );
2064   if( retDLLc ) {
2065     AliDebug(3, "----> Lc: Could not get decay length, and sigma");
2066     if (sigmaDecayLengthLc > 1e19) {
2067       codeKFLc = 3; // DecayLength not ok
2068       if (massLc < 0) {
2069         codeKFLc = 6; // DecayLength and Mass not ok
2070         if (sigmaMassLc > 1e19) codeKFLc = 11;  // DecayLength and Mass and SigmaMass not ok
2071       }
2072       else if (sigmaMassLc > 1e19) codeKFLc = 8;  // DecayLength and SigmaMass not ok
2073     }
2074   }
2075   AliDebug(3, Form("retDLLc = %d, with decayLength = %f and error = %e", retDLLc, decayLengthLc, sigmaDecayLengthLc));
2076
2077   fHistoDecayLengthKFLc->Fill(decayLengthLc, sigmaDecayLengthLc);
2078
2079   //* Get Lc life time
2080   Double_t lifeTimeLc = 999999, sigmaLifeTimeLc = 999999;
2081   Int_t retTLLc = Lc.GetLifeTime( lifeTimeLc, sigmaLifeTimeLc );
2082   if( retTLLc ) {
2083     AliDebug(3, "----> Lc: Could not get lifeTime, and sigma");
2084     if (sigmaLifeTimeLc > 1e19) {
2085       codeKFLc = 4;  // LifeTime not ok
2086       if (sigmaDecayLengthLc > 1e19) {
2087         codeKFLc = 9;  // LifeTime and DecayLength not ok
2088         if (massLc < 0) {
2089           codeKFLc = 14; // LifeTime and DecayLength and Mass not ok
2090           if (sigmaMassLc > 1e19) codeKFLc = 15; // LifeTime and DecayLength and Mass and SigmaMass not ok
2091         }
2092         else if (sigmaMassLc > 1e19) codeKFLc = 13; // LifeTime and DecayLength and SigmaMass not ok
2093       }
2094       else if (massLc < 0) { // LifeTime and Mass and SigmaMass not ok
2095         codeKFLc = 7; // LifeTime and Mass not ok
2096         if (sigmaMassLc > 1e19) codeKFLc = 12; // LifeTime and Mass and SigmaMass not ok
2097       }
2098       else if (sigmaMassLc > 1e19) codeKFLc = 10;  // LifeTime and SigmaMass not ok
2099     }
2100   }
2101
2102   fHistoLifeTimeKFLc->Fill(lifeTimeLc, sigmaLifeTimeLc);
2103
2104   AliDebug(2, Form("Lc: mass = %f (error = %e), decay length = %f (error = %e), life time = %f (error = %e) --> codeKFLc = %d", massLc, sigmaMassLc, decayLengthLc, sigmaDecayLengthLc, lifeTimeLc, sigmaLifeTimeLc, codeKFLc));
2105
2106   if (codeKFLc == -1) codeKFLc = 0;
2107   fHistoKFLc->Fill(codeKFLc);
2108
2109   fHistoKF->Fill(codeKFV0, codeKFLc);
2110
2111   // here we fill the histgrams for all the reconstructed KF vertices for the cascade
2112   fHistoMassLcAll->Fill(massLc);
2113   fHistoDecayLengthLcAll->Fill(decayLengthLc);
2114   fHistoLifeTimeLcAll->Fill(lifeTimeLc);
2115
2116   fHistoMassV0fromLcAll->Fill(massV0);
2117   fHistoDecayLengthV0fromLcAll->Fill(decayLengthV0);
2118   fHistoLifeTimeV0fromLcAll->Fill(lifeTimeV0);
2119
2120   Double_t xV0 = V0.GetX();
2121   Double_t yV0 = V0.GetY();
2122   Double_t zV0 = V0.GetZ();
2123
2124   Double_t xLc = Lc.GetX();
2125   Double_t yLc = Lc.GetY();
2126   Double_t zLc = Lc.GetZ();
2127
2128   Double_t xPrimVtx = primVtxCopy.GetX();
2129   Double_t yPrimVtx = primVtxCopy.GetY();
2130   Double_t zPrimVtx = primVtxCopy.GetZ();
2131
2132   Double_t distanceLcToPrimVtx = TMath::Sqrt((xPrimVtx - xLc) * (xPrimVtx - xLc) +
2133                                              (yPrimVtx - yLc) * (yPrimVtx - yLc) +
2134                                              (zPrimVtx - zLc) * (zPrimVtx - zLc));
2135
2136   Double_t distanceV0ToPrimVtx = TMath::Sqrt((xPrimVtx - xV0) * (xPrimVtx - xV0) +
2137                                              (yPrimVtx - yV0) * (yPrimVtx - yV0) +
2138                                              (zPrimVtx - zV0) * (zPrimVtx - zV0));
2139
2140   Double_t distanceV0ToLc = TMath::Sqrt((xLc - xV0)*(xLc - xV0) +
2141                                         (yLc - yV0)*(yLc - yV0) +
2142                                         (zLc - zV0)*(zLc - zV0));
2143
2144   //Printf("distanceLcToPrimVtx = %e, distanceV0ToPrimVtx= %f, distanceV0ToLc = %f", distanceLcToPrimVtx, distanceV0ToPrimVtx, distanceV0ToLc);
2145
2146   fHistoDistanceLcToPrimVtx->Fill(distanceLcToPrimVtx);
2147   fHistoDistanceV0ToPrimVtx->Fill(distanceV0ToPrimVtx);
2148   fHistoDistanceV0ToLc->Fill(distanceV0ToLc);
2149
2150   distances[0] = distanceLcToPrimVtx;
2151   distances[1] = distanceV0ToPrimVtx;
2152   distances[2] = distanceV0ToLc;
2153
2154   if (fUseMCInfo) {
2155
2156     AliAODMCParticle *daughv01Lc = 0x0;
2157     AliAODMCParticle *K0S = 0x0;
2158     AliAODMCParticle *daughv02Lc = 0x0;
2159
2160     if (labelsLcdaugh[0] >= 0) {
2161       //      Printf("Getting Bachelor from label %d", labelsLcdaugh[1]);
2162       daughv01Lc = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelsLcdaugh[0]));
2163       if (!daughv01Lc){
2164         AliDebug(3, "Could not access MC info for first daughter of Lc");
2165         isMCokLc = kFALSE;
2166       }
2167       else {
2168         AliDebug(2, Form("The bachelor has label = %d", daughv01Lc->GetLabel()));
2169         if (TMath::Abs(daughv01Lc->GetPdgCode()) != 2212) isBkgLc = kTRUE;
2170       }
2171     }
2172     else { // The bachelor is a fake
2173       isBkgLc = kTRUE;
2174     }
2175
2176     if (labelsLcdaugh[1] >= 0) {
2177       //Printf("Getting K0S");
2178       K0S = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelsLcdaugh[1]));
2179       if (!K0S) {
2180         AliDebug(3, "Could not access MC info for second daughter of Lc");
2181         isMCokLc = kFALSE;
2182       }
2183       else{
2184         if (TMath::Abs(K0S->GetPdgCode()) != 310) isBkgLc = kTRUE;
2185       }
2186     }
2187     else{
2188       AliDebug(2, "The K0S is not true --> it does not have a label, continuing...");
2189       isBkgLc = kTRUE;
2190     }
2191
2192     if (!isBkgLc){ // so far, we only checked that the V0 and the bachelor are not fake, and in particular, we know that the V0 is a K0S since we used the MatchToMC method
2193       if (isMCokLc) { // We can then access its MC info, and it might then be that also the Lc is a true Lc
2194         Int_t iK0 = K0S->GetMother();
2195         if (iK0 < 0) {
2196           Printf("The K0S has no mother... IMPOSSIBLE"); // the K0S MUST have a mother!
2197         }
2198         else { // The K0S has a mother
2199           daughv02Lc = dynamic_cast<AliAODMCParticle*>(mcArray->At(iK0));
2200           if (!daughv02Lc){
2201             AliDebug(3, "Could not access MC info for second daughter of Lc");
2202           }
2203           else { // we can access safely the K0S mother in the MC
2204             if( daughv01Lc && (daughv01Lc->GetMother() ==  daughv02Lc->GetMother()) && (daughv01Lc->GetMother()>=0) ){  // This is a true cascade! bachelor and V0 come from the same mother
2205               //Printf("Lc: The mother has label %d", daughv01Lc->GetMother());
2206               AliAODMCParticle *motherLc = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughv01Lc->GetMother()));
2207               Int_t pdgMum = 0, pdgBach = 0, pdgV0 = 0;
2208               if (motherLc) pdgMum = motherLc->GetPdgCode();
2209               if (daughv01Lc) pdgBach = daughv01Lc->GetPdgCode();
2210               if (daughv02Lc) pdgV0 = daughv02Lc->GetPdgCode();
2211               AliDebug(2, Form("pdgMum = %d, pdgBach = %d, pdgV0 = %d", pdgMum, pdgBach, pdgV0));
2212
2213               if( motherLc && TMath::Abs(motherLc->GetPdgCode()) != 21 ){ // These are true cascades, we don't know yet if they are Lc
2214                 fHistoMassLcTrue->Fill(massLc);
2215                 fHistoDecayLengthLcTrue->Fill(decayLengthLc);
2216                 fHistoLifeTimeLcTrue->Fill(lifeTimeLc);
2217                 fHistoMassLcTrueFromAOD->Fill(cascade->InvMassLctoK0sP());
2218
2219                 fHistoMassV0fromLcTrue->Fill(massV0);
2220                 fHistoDecayLengthV0fromLcTrue->Fill(decayLengthV0);
2221                 fHistoLifeTimeV0fromLcTrue->Fill(lifeTimeV0);
2222
2223                 if (TMath::Abs(motherLc->GetPdgCode()) == 4122 && TMath::Abs(motherV0->GetPdgCode()) == 310 && TMath::Abs(daughv01Lc->GetPdgCode()) == 2212){  // This is Lc --> K0S + p (the check on the PDG code of the V0 is useless, since we used MathcToMC with it, but fine...
2224                   AliDebug(2, Form("IT IS SIGNAL!!! with label = %d", motherLc->GetLabel()));
2225
2226                   fHistoArmenterosPodolanskiV0KFSgn->Fill(qtAlphaV0[1], qtAlphaV0[0]);
2227                   fHistoArmenterosPodolanskiV0AODSgn->Fill(v0part->AlphaV0(), v0part->PtArmV0());
2228
2229                   fHistoDistanceLcToPrimVtxSgn->Fill(distanceLcToPrimVtx);
2230                   fHistoDistanceV0ToPrimVtxSgn->Fill(distanceV0ToPrimVtx);
2231                   fHistoDistanceV0ToLcSgn->Fill(distanceV0ToLc);
2232
2233                   fHistoMassLcSgn->Fill(massLc);
2234                   fHistoMassLcSgnFromAOD->Fill(cascade->InvMassLctoK0sP());
2235
2236                   fHistoDecayLengthLcSgn->Fill(decayLengthLc);
2237                   fHistoLifeTimeLcSgn->Fill(lifeTimeLc);
2238
2239                   fHistoMassV0fromLcSgn->Fill(massV0);
2240                   fHistoDecayLengthV0fromLcSgn->Fill(decayLengthV0);
2241                   fHistoLifeTimeV0fromLcSgn->Fill(lifeTimeV0);
2242                 }
2243                 else {
2244                   isBkgLc = kTRUE;  // it is not a Lc, since the pdg != 4122
2245                 }
2246
2247                 // if we got here, we can compare with MC information; this part is done also for cases where the cascade is not a Lc!
2248                 // First, we compare the vtx of the Lc
2249                 Double_t xLcMC = motherLc->Xv();
2250                 Double_t yLcMC = motherLc->Yv();
2251                 Double_t zLcMC = motherLc->Zv();
2252                 //Printf("Got MC vtx for Lc");
2253                 Double_t xProtonMC = daughv01Lc->Xv();
2254                 Double_t yProtonMC = daughv01Lc->Yv();
2255                 Double_t zProtonMC = daughv01Lc->Zv();
2256                 //Printf("Got MC vtx for Proton");
2257
2258                 Double_t vtxLcResidualToPrimVtx = TMath::Sqrt((xLcMC - xProtonMC) * (xLcMC - xProtonMC) +
2259                                                               (yLcMC - yProtonMC) * (yLcMC - yProtonMC) +
2260                                                               (zLcMC - zProtonMC) * (zLcMC - zProtonMC)) - distanceLcToPrimVtx;
2261
2262                 // Then, we compare the vtx of the K0s
2263                 Double_t xV0MC = motherV0->Xv();     //Production vertex of K0S --> Where Lc decays
2264                 Double_t yV0MC = motherV0->Yv();
2265                 Double_t zV0MC = motherV0->Zv();
2266
2267                 //Printf("Got MC vtx for V0");
2268                 Double_t xPionMC = daughv01->Xv(); //Production vertex of Pion --> Where K0S decays
2269                 Double_t yPionMC = daughv01->Yv();
2270                 Double_t zPionMC = daughv01->Zv();
2271                 //Printf("Got MC vtx for Pion");
2272                 Printf("Vertices: MC:  x = %f, y = %f, z = %f", xPionMC, yPionMC, zPionMC);
2273
2274                 Double_t vtxV0ResidualToLc = TMath::Sqrt((xV0MC - xPionMC) * (xV0MC - xPionMC) +
2275                                                          (yV0MC - yPionMC) * (yV0MC - yPionMC) +
2276                                                          (zV0MC - zPionMC) * (zV0MC - zPionMC)) - distanceV0ToLc;
2277
2278                 // Then, the K0S vertex wrt primary vertex
2279
2280                 Double_t vtxV0ResidualToPrimVtx = TMath::Sqrt((xPionMC - xLcMC) * (xPionMC - xLcMC) +
2281                                                               (yPionMC - yLcMC) * (yPionMC - yLcMC) +
2282                                                               (zPionMC - zLcMC) * (zPionMC - zLcMC)) - distanceV0ToPrimVtx;
2283
2284                 fHistoVtxLcResidualToPrimVtx->Fill(vtxLcResidualToPrimVtx);
2285                 fHistoVtxV0ResidualToLc->Fill(vtxV0ResidualToLc);
2286                 fHistoVtxV0ResidualToPrimVtx->Fill(vtxV0ResidualToPrimVtx);
2287
2288               } //closing: if( motherLc && TMath::Abs(motherLc->GetPdgCode()) != 21 )
2289               else if (!motherLc){
2290                 AliDebug(2, "We could not access MC info for Lc mother, so we did nothing");
2291               }
2292               else {
2293                 AliDebug(2, "MC Lc mother is a gluon, so we did nothing");
2294               }
2295             } //closing: if( daughv01Lc->GetMother() ==  daughv02Lc->GetMother() && daughv01Lc->GetMother()>=0 )
2296             else {
2297               isBkgLc = kTRUE; // it cannot be a Lc, since the daughters do not have the same mother
2298             }
2299           } // closing: else { // we can access safely the K0S mother in the MC
2300         } // closing: else { // The K0S has a mother
2301       } // closing isMCLcok
2302     }  // closing !isBkgLc
2303   } // closing fUseMCInfo
2304
2305   //Printf("retMV0 = %d, retMLc = %d", retMV0, retMLc);
2306   if ( retMV0 == 0 && retMLc == 0){
2307     V0KF[0] = massV0;
2308     errV0KF[0] = sigmaMassV0;
2309     V0KF[1] = decayLengthV0;
2310     errV0KF[1] = sigmaDecayLengthV0;
2311     V0KF[2] = lifeTimeV0;
2312     errV0KF[2] = sigmaLifeTimeV0;
2313     LcKF[0] = massLc;
2314     errLcKF[0] = sigmaMassLc;
2315     LcKF[1] = decayLengthLc;
2316     errLcKF[1] = sigmaDecayLengthLc;
2317     LcKF[2] = lifeTimeLc;
2318     errLcKF[2] = sigmaLifeTimeLc;
2319   }
2320
2321   return 1;
2322
2323 }
2324 //________________________________________________________________________
2325 AliAnalysisTaskSELc2V0bachelorTMVA::EBachelor AliAnalysisTaskSELc2V0bachelorTMVA::CheckBachelor( AliAODRecoCascadeHF *part,
2326                                                                                                  AliAODTrack* bachelor,
2327                                                                                                  TClonesArray *mcArray ){
2328
2329   //Printf("In CheckBachelor");
2330
2331   // function to check if the bachelor is a p, if it is a p but not from Lc
2332   // to be filled for background candidates
2333
2334   Int_t label = bachelor->GetLabel();
2335   if (label == -1) {
2336     return kBachFake;
2337   }
2338
2339   AliAODMCParticle *mcpart = dynamic_cast<AliAODMCParticle*>(mcArray->At(TMath::Abs(label)));
2340   if (!mcpart) return kBachInvalid;
2341   Int_t pdg = mcpart->PdgCode();
2342   if (TMath::Abs(pdg) != 2212) {
2343     AliDebug(2, Form("Bachelor is not a p, but a particle with pdg code =  %d", pdg));
2344     return kBachNoProton;
2345   }
2346   else { // it is a proton
2347     //Int_t labelLc = part->GetLabel();
2348     Int_t labelLc = FindLcLabel(part, mcArray);
2349     //Printf(">>>>>>>>>>>>> label for Lc = %d", labelLc);
2350     Int_t bachelorMotherLabel = mcpart->GetMother();
2351     //Printf(">>>>>>>>>>>>> label for bachelorMother = %d", bachelorMotherLabel);
2352     if (bachelorMotherLabel == -1) {
2353       return kBachPrimary;
2354     }
2355     AliAODMCParticle *bachelorMother = dynamic_cast<AliAODMCParticle*>(mcArray->At(bachelorMotherLabel));
2356     if (!bachelorMother) return kBachInvalid;
2357     Int_t pdgMother = bachelorMother->PdgCode();
2358     if (TMath::Abs(pdgMother) != 4122) {
2359       AliDebug(2, Form("The proton does not come from a Lc, but from a particle with pdgcode = %d", pdgMother));
2360       return kBachNoLambdaMother;
2361     }
2362     else { // it comes from Lc
2363       if (labelLc != bachelorMotherLabel){
2364         //AliInfo(Form("The proton comes from a Lc, but it is not the candidate we are analyzing (label Lc = %d, label p mother = %d", labelLc, bachelorMotherLabel));
2365         AliDebug(2, Form("The proton comes from a Lc, but it is not the candidate we are analyzing"));
2366         return kBachDifferentLambdaMother;
2367       }
2368       else { // it comes from the correct Lc
2369         AliDebug(2, Form("The proton comes from a Lc, and it is EXACTLY the candidate we are analyzing"));
2370         return kBachCorrectLambdaMother;
2371       }
2372     }
2373   }
2374
2375   return kBachInvalid;
2376
2377 }
2378
2379 //________________________________________________________________________
2380 AliAnalysisTaskSELc2V0bachelorTMVA::EK0S AliAnalysisTaskSELc2V0bachelorTMVA::CheckK0S( AliAODRecoCascadeHF *part,
2381                                                                                        AliAODv0* v0part,
2382                                                                                        //AliAODTrack* v0part,
2383                                                                                        TClonesArray *mcArray ){
2384
2385   // function to check if the K0Spart is a p, if it is a p but not from Lc
2386   // to be filled for background candidates
2387
2388   //Printf(" CheckK0S");
2389
2390   //Int_t labelMatchToMC = v0part->MatchToMC(310, mcArray);
2391   //Int_t label = v0part->GetLabel();
2392   Int_t labelFind = FindV0Label(v0part, mcArray);
2393   //Printf("\n >>>>>>>>>>>>> label for V0 = %d, from MatchToMC = %d, from Find = %d", label, labelMatchToMC, labelFind);
2394   if (labelFind == -1) {
2395     return kK0SFake;
2396   }
2397
2398   AliAODMCParticle *mcpart = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelFind));
2399   if (!mcpart) return kK0SInvalid;
2400   Int_t pdg = mcpart->PdgCode();
2401   if (TMath::Abs(pdg) != 310) {
2402     AliDebug(2, Form("K0Spart is not a K0S, but a particle with pdg code =  %d", pdg));
2403     //AliInfo(Form("K0Spart is not a K0S, but a particle with pdg code =  %d", pdg));
2404     return kK0SNoK0S;
2405   }
2406   else { // it is a K0S
2407     //Int_t labelLc = part->GetLabel();
2408     Int_t labelLc = FindLcLabel(part, mcArray);
2409     Int_t K0SpartMotherLabel = mcpart->GetMother();
2410     if (K0SpartMotherLabel == -1) {
2411       return kK0SWithoutMother;
2412     }
2413     AliAODMCParticle *K0SpartMother = dynamic_cast<AliAODMCParticle*>(mcArray->At(K0SpartMotherLabel));  // should be a K0 in case of signal
2414     if (!K0SpartMother) return kK0SInvalid;
2415     Int_t pdgMotherK0S = K0SpartMother->PdgCode();
2416     if (TMath::Abs(pdgMotherK0S) != 311) {
2417       AliDebug(2, Form("The K0S does not come from a K0, but from a particle with pdgcode = %d", pdgMotherK0S));
2418       //AliInfo(Form("The K0S does not come from a K0, but from a particle with pdgcode = %d", pdgMotherK0S));
2419       return kK0SNotFromK0;
2420     }
2421     else { // the K0S comes from a K0
2422       Int_t K0MotherLabel = K0SpartMother->GetMother(); // mother of K0 --> Lc in case of signal
2423       if (K0MotherLabel == -1) {
2424         return kK0Primary;
2425       }
2426       AliAODMCParticle *K0Mother = dynamic_cast<AliAODMCParticle*>(mcArray->At(K0MotherLabel));
2427       if (!K0Mother) return kK0SInvalid;
2428       Int_t pdgK0Mother = K0Mother->PdgCode();
2429       if (TMath::Abs(pdgK0Mother) != 4122) { // the K0 does not come from a Lc
2430         AliDebug(2, Form("The K0 does not come from a Lc, but from a particle with pdgcode = %d", pdgK0Mother));
2431         //AliInfo(Form("The K0 does not come from a Lc, but from a particle with pdgcode = %d", pdgK0Mother));
2432         return kK0NoLambdaMother;
2433       }
2434       else { // the K0 comes from Lc
2435         if (labelLc != K0MotherLabel){ // The K0 comes from a different Lc
2436           AliDebug(2, Form("The K0S comes from a Lc, but it is not the candidate we are analyzing"));
2437           //AliInfo(Form("The K0S comes from a Lc, but it is not the candidate we are analyzing"));
2438           return kK0DifferentLambdaMother;
2439         }
2440         else { // it comes from the correct Lc
2441           AliDebug(2, Form("The K0S comes from a Lc, and it is EXACTLY the candidate we are analyzing"));
2442           //AliInfo(Form("The K0S comes from a Lc, and it is EXACTLY the candidate we are analyzing"));
2443           return kK0CorrectLambdaMother;
2444         }
2445       }
2446     }
2447   }
2448
2449   return kK0SInvalid;
2450
2451 }
2452
2453 //----------------------------------------------------------------------------
2454 Int_t AliAnalysisTaskSELc2V0bachelorTMVA::FindV0Label(AliAODRecoDecay* v0part, TClonesArray *mcArray) const
2455 {
2456
2457   //Printf(" FindV0Label");
2458
2459   // finding the label of teh V0; inspired from AliAODRecoDecay::MatchToMC
2460
2461   Int_t labMother[2]={-1, -1};
2462   AliAODMCParticle *part=0;
2463   AliAODMCParticle *mother=0;
2464   Int_t dgLabels = -1;
2465
2466   Int_t ndg = v0part->GetNDaughters();
2467   if (ndg != 2) {
2468     AliFatal(Form("IMPOSSIBLE!! There are %d daughters, but they should be 2!", ndg));
2469   }
2470
2471   for(Int_t i = 0; i < 2; i++) {
2472     AliAODTrack *trk = (AliAODTrack*)v0part->GetDaughter(i);
2473     dgLabels = trk->GetLabel();
2474     if (dgLabels == -1) {
2475       //printf("daughter with negative label %d\n", dgLabels);
2476       return -1;
2477     }
2478     part = (AliAODMCParticle*)mcArray->At(TMath::Abs(dgLabels));
2479     if (!part) {
2480       //printf("no MC particle\n");
2481       return -1;
2482     }
2483     labMother[i] = part->GetMother();
2484     if (labMother[i] != -1){
2485       mother = (AliAODMCParticle*)mcArray->At(labMother[i]);
2486       if(!mother) {
2487         //printf("no MC mother particle\n");
2488         return -1;
2489       }
2490     }
2491     else {
2492       return -1;
2493     }
2494   }
2495
2496   if (labMother[0] == labMother[1]) return labMother[0];
2497
2498   else return -1;
2499
2500 }
2501
2502 //----------------------------------------------------------------------------
2503 Int_t AliAnalysisTaskSELc2V0bachelorTMVA::FindLcLabel(AliAODRecoCascadeHF* cascade, TClonesArray *mcArray) const
2504 {
2505
2506   // finding the label of teh V0; inspired from AliAODRecoDecay::MatchToMC
2507
2508   //Printf(" FindLcLabel");
2509
2510   AliAODMCParticle *part=0;
2511   AliAODMCParticle *mother=0;
2512   AliAODMCParticle *grandmother=0;
2513   Int_t dgLabels = -1;
2514
2515   Int_t ndg = cascade->GetNDaughters();
2516   if (ndg != 2) {
2517     AliFatal(Form("IMPOSSIBLE!! There are %d daughters, but they should be 2!", ndg));
2518   }
2519
2520   // daughter 0 --> bachelor, daughter 1 --> V0
2521   AliAODTrack *trk = (AliAODTrack*)cascade->GetDaughter(0); // bachelor
2522   dgLabels = trk->GetLabel();
2523   if (dgLabels == -1 ) {
2524     //printf("daughter with negative label %d\n", dgLabels);
2525     return -1;
2526   }
2527   part = (AliAODMCParticle*)mcArray->At(TMath::Abs(dgLabels));
2528   if (!part) {
2529     //printf("no MC particle\n");
2530     return -1;
2531   }
2532   Int_t labMotherBach = part->GetMother();
2533   if (labMotherBach == -1){
2534     return -1;
2535   }
2536   mother = (AliAODMCParticle*)mcArray->At(labMotherBach);
2537   if(!mother) {
2538     //printf("no MC mother particle\n");
2539     return -1;
2540   }
2541
2542   AliAODv0 *v0 = (AliAODv0*)cascade->GetDaughter(1); // V0
2543   dgLabels = FindV0Label(v0, mcArray);
2544   if (dgLabels == -1 ) {
2545     //printf("daughter with negative label (v0 was a fake) %d\n", dgLabels);
2546     return -1;
2547   }
2548   part = (AliAODMCParticle*)mcArray->At(TMath::Abs(dgLabels));
2549   if (!part) {
2550     //printf("no MC particle\n");
2551     return -1;
2552   }
2553   Int_t labMotherv0 = part->GetMother();
2554   if (labMotherv0 == -1){
2555     return -1;
2556   }
2557   mother = (AliAODMCParticle*)mcArray->At(labMotherv0);
2558   if(!mother) {
2559     //printf("no MC mother particle\n");
2560     return -1;
2561   }
2562   Int_t labGrandMotherv0 = mother->GetMother();
2563   if (labGrandMotherv0 == -1){
2564     return -1;
2565   }
2566   grandmother = (AliAODMCParticle*)mcArray->At(labGrandMotherv0);
2567   if(!grandmother) {
2568     //printf("no MC mother particle\n");
2569     return -1;
2570   }
2571
2572   //Printf("labMotherBach = %d, labMotherv0 = %d, labGrandMotherv0 = %d", labMotherBach, labMotherv0, labGrandMotherv0);
2573   if (labGrandMotherv0 == labMotherBach) return labMotherBach;
2574
2575   else return -1;
2576
2577 }
2578
2579
2580
2581
2582
2583