PWGHF/vertexingHF/AliAnalysisTaskSELc2V0bachelorTMVA.cxx

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