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