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