Task for Lc->K0s+p with TMVA (Chiara Z.)
authorfprino <prino@to.infn.it>
Wed, 12 Mar 2014 17:07:44 +0000 (18:07 +0100)
committerfprino <prino@to.infn.it>
Wed, 12 Mar 2014 17:07:59 +0000 (18:07 +0100)
PWGHF/CMakelibPWGHFvertexingHF.pkg
PWGHF/PWGHFvertexingHFLinkDef.h
PWGHF/vertexingHF/AliAnalysisTaskSELc2V0bachelorTMVA.cxx [new file with mode: 0644]
PWGHF/vertexingHF/AliAnalysisTaskSELc2V0bachelorTMVA.h [new file with mode: 0644]
PWGHF/vertexingHF/macros/AddTaskLc2V0bachpA_TMVA_KF_ON_OnlyHIJINGBkg.C [new file with mode: 0644]

index 07dc2ba..b7da08f 100644 (file)
@@ -84,6 +84,7 @@ set ( SRCS
     vertexingHF/charmFlow/AliAnalysisTaskFlowD2H.cxx
     vertexingHF/upgrade/AliAnalysisTaskSELambdacUp.cxx
     vertexingHF/upgrade/AliAnalysisTaskCountLcEta.cxx
+    vertexingHF/AliAnalysisTaskSELc2V0bachelorTMVA.cxx 
     )
 
 string ( REPLACE ".cxx" ".h" HDRS "${SRCS}" )
index 82429dd..0a968b8 100644 (file)
@@ -62,6 +62,7 @@
 #pragma link C++ class AliHFAfterBurner+;
 #pragma link C++ class AliAnalysisTaskSELambdacUp+;
 #pragma link C++ class AliAnalysisTaskCountLcEta+;
+#pragma link C++ class AliAnalysisTaskSELc2V0bachelorTMVA+;
 
 
 #endif
diff --git a/PWGHF/vertexingHF/AliAnalysisTaskSELc2V0bachelorTMVA.cxx b/PWGHF/vertexingHF/AliAnalysisTaskSELc2V0bachelorTMVA.cxx
new file mode 100644 (file)
index 0000000..b13c1a8
--- /dev/null
@@ -0,0 +1,2249 @@
+/**************************************************************************
+ * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appeuear in the supporting documentation. The authors make no claims   *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/* $Id: AliAnalysisTaskSELc2V0bachelorTMVA.cxx 62231 2013-04-29 09:47:26Z fprino $ */
+
+//
+//
+//               Base class for Lc2V0 Analysis to be used with TMVA
+//
+//
+
+//------------------------------------------------------------------------------------------
+//
+//  Author: C. Zampolli (from AliAnalysisTaskSELc2V0bachelor class by A.De Caro, P. Pagano)
+//
+//------------------------------------------------------------------------------------------
+
+#include <TSystem.h>
+#include <TParticle.h>
+#include <TParticlePDG.h>
+#include <TH1F.h>
+#include <TH1I.h>
+#include <TH2F.h>
+#include <TTree.h>
+#include "TROOT.h"
+#include <TDatabasePDG.h>
+#include <AliAnalysisDataSlot.h>
+#include <AliAnalysisDataContainer.h>
+#include "AliStack.h"
+#include "AliMCEvent.h"
+#include "AliAnalysisManager.h"
+#include "AliAODMCHeader.h"
+#include "AliAODHandler.h"
+#include "AliLog.h"
+#include "AliAODVertex.h"
+#include "AliAODRecoDecay.h"
+#include "AliAODRecoDecayHF.h"
+#include "AliAODRecoCascadeHF.h"
+#include "AliAnalysisVertexingHF.h"
+#include "AliESDtrack.h"
+#include "AliAODTrack.h"
+#include "AliAODv0.h"
+#include "AliAODMCParticle.h"
+#include "AliAnalysisTaskSE.h"
+#include "AliAnalysisTaskSELc2V0bachelorTMVA.h"
+#include "AliNormalizationCounter.h"
+#include "AliAODPidHF.h"
+#include "AliPIDResponse.h"
+#include "AliPIDCombined.h"
+#include "AliTOFPIDResponse.h"
+#include "AliInputEventHandler.h"
+#include "AliAODRecoDecayHF3Prong.h"
+#include "AliKFParticle.h"
+#include "AliKFVertex.h"
+#include "AliVertexingHFUtils.h"
+
+using std::cout;
+using std::endl;
+
+ClassImp(AliAnalysisTaskSELc2V0bachelorTMVA)
+
+//__________________________________________________________________________
+AliAnalysisTaskSELc2V0bachelorTMVA::AliAnalysisTaskSELc2V0bachelorTMVA():  
+AliAnalysisTaskSE(),
+  fUseMCInfo(kFALSE),
+  fOutput(0),
+  fCEvents(0),
+  fPIDResponse(0),
+  fPIDCombined(0),
+  fIsK0sAnalysis(kFALSE),
+  fCounter(0),
+  fAnalCuts(0),
+  fListCuts(0),
+  fUseOnTheFlyV0(kFALSE),
+  fIsEventSelected(kFALSE),
+  fVariablesTreeSgn(0),
+  fVariablesTreeBkg(0),
+  fCandidateVariables(),
+  fIspA(kFALSE),
+  fHistoEvents(0),
+  fHistoLc(0),
+  fHistoLcOnTheFly(0),
+  fFillOnlySgn(kFALSE),
+fHistoLcBeforeCuts(0),
+fHistoFiducialAcceptance(0), 
+fHistoCodesSgn(0),
+fHistoCodesBkg(0),
+  fHistoLcpKpiBeforeCuts(0),
+  fVtx1(0),
+
+  fHistoDistanceLcToPrimVtx(0),       
+  fHistoDistanceV0ToPrimVtx(0),      
+  fHistoDistanceV0ToLc(0),            
+
+  fHistoDistanceLcToPrimVtxSgn(0),    
+  fHistoDistanceV0ToPrimVtxSgn(0),    
+  fHistoDistanceV0ToLcSgn(0),
+
+  fHistoVtxLcResidualToPrimVtx(0),    
+  fHistoVtxV0ResidualToPrimVtx(0),    
+  fHistoVtxV0ResidualToLc(0),         
+
+  fHistoMassV0All(0),             
+  fHistoDecayLengthV0All(0),      
+  fHistoLifeTimeV0All(0),         
+
+  fHistoMassV0True(0),             
+  fHistoDecayLengthV0True(0),      
+  fHistoLifeTimeV0True(0),         
+
+  fHistoMassV0TrueFromAOD(0),             
+
+  fHistoMassV0TrueK0S(0),             
+  fHistoDecayLengthV0TrueK0S(0),      
+  fHistoLifeTimeV0TrueK0S(0),         
+
+  fHistoMassV0TrueK0SFromAOD(0),             
+
+  fHistoMassLcAll(0),             
+  fHistoDecayLengthLcAll(0),      
+  fHistoLifeTimeLcAll(0),         
+
+  fHistoMassLcTrue(0),             
+  fHistoDecayLengthLcTrue(0),      
+  fHistoLifeTimeLcTrue(0),         
+
+  fHistoMassLcTrueFromAOD(0),             
+
+  fHistoMassV0fromLcAll(0),       
+  fHistoDecayLengthV0fromLcAll(0),
+  fHistoLifeTimeV0fromLcAll(0),   
+
+  fHistoMassV0fromLcTrue(0),       
+  fHistoDecayLengthV0fromLcTrue(0),
+  fHistoLifeTimeV0fromLcTrue(0),   
+
+  fHistoMassLcSgn(0),             
+  fHistoMassLcSgnFromAOD(0),             
+  fHistoDecayLengthLcSgn(0),      
+  fHistoLifeTimeLcSgn(0),         
+
+  fHistoMassV0fromLcSgn(0),       
+  fHistoDecayLengthV0fromLcSgn(0),
+  fHistoLifeTimeV0fromLcSgn(0),
+
+  fHistoKF(0),
+  fHistoKFV0(0),
+  fHistoKFLc(0),
+
+  fHistoMassKFV0(0), 
+  fHistoDecayLengthKFV0(0), 
+  fHistoLifeTimeKFV0(0), 
+
+  fHistoMassKFLc(0), 
+  fHistoDecayLengthKFLc(0), 
+  fHistoLifeTimeKFLc(0), 
+
+  fHistoArmenterosPodolanskiV0KF(0),
+  fHistoArmenterosPodolanskiV0KFSgn(0),
+  fHistoArmenterosPodolanskiV0AOD(0),
+  fHistoArmenterosPodolanskiV0AODSgn(0),
+
+  fOutputKF(0),
+  fmcLabelLc(-1),
+  ftopoConstraint(kTRUE),
+  fCallKFVertexing(kFALSE),
+  fKeepingOnlyHIJINGBkg(kFALSE),
+  fUtils(0),
+  fHistoBackground(0) 
+
+{
+  //
+  // Default ctor
+  //
+}
+//___________________________________________________________________________
+AliAnalysisTaskSELc2V0bachelorTMVA::AliAnalysisTaskSELc2V0bachelorTMVA(const Char_t* name,
+                                                                      AliRDHFCutsLctoV0* analCuts, Bool_t useOnTheFly) :
+  AliAnalysisTaskSE(name),
+  fUseMCInfo(kFALSE),
+  fOutput(0),
+  fCEvents(0),
+  fPIDResponse(0),
+  fPIDCombined(0),
+  fIsK0sAnalysis(kFALSE),
+  fCounter(0),
+  fAnalCuts(analCuts),
+  fListCuts(0),
+  fUseOnTheFlyV0(useOnTheFly),
+  fIsEventSelected(kFALSE),
+  fVariablesTreeSgn(0),
+  fVariablesTreeBkg(0),
+  fCandidateVariables(),
+  fIspA(kFALSE),
+  fHistoEvents(0),
+  fHistoLc(0),
+  fHistoLcOnTheFly(0),
+  fFillOnlySgn(kFALSE),
+  fHistoLcBeforeCuts(0),
+  fHistoFiducialAcceptance(0), 
+  fHistoCodesSgn(0),
+  fHistoCodesBkg(0),
+  fHistoLcpKpiBeforeCuts(0),
+  fVtx1(0),
+
+  fHistoDistanceLcToPrimVtx(0),       
+  fHistoDistanceV0ToPrimVtx(0),      
+  fHistoDistanceV0ToLc(0),            
+
+  fHistoDistanceLcToPrimVtxSgn(0),    
+  fHistoDistanceV0ToPrimVtxSgn(0),    
+  fHistoDistanceV0ToLcSgn(0),
+
+  fHistoVtxLcResidualToPrimVtx(0),    
+  fHistoVtxV0ResidualToPrimVtx(0),    
+  fHistoVtxV0ResidualToLc(0),         
+
+  fHistoMassV0All(0),             
+  fHistoDecayLengthV0All(0),      
+  fHistoLifeTimeV0All(0),         
+
+  fHistoMassV0True(0),             
+  fHistoDecayLengthV0True(0),      
+  fHistoLifeTimeV0True(0),   
+
+  fHistoMassV0TrueFromAOD(0),                   
+
+  fHistoMassV0TrueK0S(0),             
+  fHistoDecayLengthV0TrueK0S(0),      
+  fHistoLifeTimeV0TrueK0S(0),   
+
+  fHistoMassV0TrueK0SFromAOD(0),                   
+
+  fHistoMassLcAll(0),             
+  fHistoDecayLengthLcAll(0),      
+  fHistoLifeTimeLcAll(0),         
+
+  fHistoMassLcTrue(0),             
+  fHistoDecayLengthLcTrue(0),      
+  fHistoLifeTimeLcTrue(0),         
+
+  fHistoMassLcTrueFromAOD(0),             
+
+  fHistoMassV0fromLcAll(0),       
+  fHistoDecayLengthV0fromLcAll(0),
+  fHistoLifeTimeV0fromLcAll(0),   
+
+  fHistoMassV0fromLcTrue(0),       
+  fHistoDecayLengthV0fromLcTrue(0),
+  fHistoLifeTimeV0fromLcTrue(0),   
+
+  fHistoMassLcSgn(0),             
+  fHistoMassLcSgnFromAOD(0),             
+  fHistoDecayLengthLcSgn(0),      
+  fHistoLifeTimeLcSgn(0),         
+
+  fHistoMassV0fromLcSgn(0),       
+  fHistoDecayLengthV0fromLcSgn(0),
+  fHistoLifeTimeV0fromLcSgn(0),
+  fHistoKF(0),
+  fHistoKFV0(0),
+  fHistoKFLc(0),
+
+  fHistoMassKFV0(0), 
+  fHistoDecayLengthKFV0(0), 
+  fHistoLifeTimeKFV0(0), 
+
+  fHistoMassKFLc(0), 
+  fHistoDecayLengthKFLc(0), 
+  fHistoLifeTimeKFLc(0), 
+
+  fHistoArmenterosPodolanskiV0KF(0),
+  fHistoArmenterosPodolanskiV0KFSgn(0),
+  fHistoArmenterosPodolanskiV0AOD(0),
+  fHistoArmenterosPodolanskiV0AODSgn(0),
+
+  fOutputKF(0),
+  fmcLabelLc(-1),
+  ftopoConstraint(kTRUE),
+  fCallKFVertexing(kFALSE),
+  fKeepingOnlyHIJINGBkg(kFALSE),
+  fUtils(0),
+  fHistoBackground(0) 
+
+
+{
+  //
+  // Constructor. Initialization of Inputs and Outputs
+  //
+  Info("AliAnalysisTaskSELc2V0bachelorTMVA","Calling Constructor");
+
+  DefineOutput(1, TList::Class());  // Tree signal + Tree Bkg + histoEvents
+  DefineOutput(2, AliNormalizationCounter::Class()); // normalization counter object
+  DefineOutput(3, TList::Class());  // Cuts
+  DefineOutput(4, TTree::Class());  // Tree signal + Tree Bkg + histoEvents
+  DefineOutput(5, TTree::Class());  // Tree signal + Tree Bkg + histoEvents
+  DefineOutput(6, TList::Class());  // Tree signal + Tree Bkg + histoEvents
+
+}
+
+//___________________________________________________________________________
+AliAnalysisTaskSELc2V0bachelorTMVA::~AliAnalysisTaskSELc2V0bachelorTMVA() {
+  //
+  // destructor
+  //
+  Info("~AliAnalysisTaskSELc2V0bachelorTMVA","Calling Destructor");
+  
+  if (fOutput) {
+    delete fOutput;
+    fOutput = 0;
+  }
+
+  if (fPIDResponse) {
+    delete  fPIDResponse;
+  }
+
+  if (fPIDCombined) {
+    delete  fPIDCombined;
+  }
+
+  if (fCounter) {
+    delete fCounter;
+    fCounter = 0;
+  }
+
+  if (fAnalCuts) {
+    delete fAnalCuts;
+    fAnalCuts = 0;
+  }
+
+  if (fListCuts) {
+    delete fListCuts;
+    fListCuts = 0;
+  }
+
+  if(fVariablesTreeSgn){
+    delete fVariablesTreeSgn;
+    fVariablesTreeSgn = 0;
+  }
+
+  if(fVariablesTreeBkg){
+    delete fVariablesTreeBkg;
+    fVariablesTreeBkg = 0;
+  }
+
+  if (fOutputKF) {
+    delete fOutputKF;
+    fOutputKF = 0;
+  }
+
+  if (fUtils) {
+    delete fUtils;
+    fUtils = 0;
+  }
+
+}
+//_________________________________________________
+void AliAnalysisTaskSELc2V0bachelorTMVA::Init() {
+  //
+  // Initialization
+  //
+
+  fIsEventSelected=kFALSE;
+
+  if (fDebug > 1) AliInfo("Init");
+
+  fListCuts = new TList();
+  fListCuts->SetOwner();
+  fListCuts->Add(new AliRDHFCutsLctoV0(*fAnalCuts));
+  PostData(3,fListCuts);
+
+  if (fUseMCInfo && fKeepingOnlyHIJINGBkg) fUtils = new AliVertexingHFUtils();
+
+  return;
+}
+
+//________________________________________ terminate ___________________________
+void AliAnalysisTaskSELc2V0bachelorTMVA::Terminate(Option_t*)
+{    
+  // The Terminate() function is the last function to be called during
+  // a query. It always runs on the client, it can be used to present
+  // the results graphically or save the results to file.
+  
+  //AliInfo("Terminate","");
+  AliAnalysisTaskSE::Terminate();
+  
+  fOutput = dynamic_cast<TList*> (GetOutputData(1));
+  if (!fOutput) {     
+    AliError("fOutput not available");
+    return;
+  }
+  fOutputKF = dynamic_cast<TList*> (GetOutputData(6));
+  if (!fOutputKF) {     
+    AliError("fOutputKF not available");
+    return;
+  }
+  
+  return;
+}
+
+//___________________________________________________________________________
+void AliAnalysisTaskSELc2V0bachelorTMVA::UserCreateOutputObjects() { 
+  // output
+  AliInfo(Form("CreateOutputObjects of task %s\n", GetName()));
+  
+  //slot #1  
+  //OpenFile(1);
+  fOutput = new TList();
+  fOutput->SetOwner();
+  fOutput->SetName("listTrees");
+
+  // Output slot 1: list of 2 trees (Sgn + Bkg) of the candidate variables
+  const char* nameoutput = GetOutputSlot(1)->GetContainer()->GetName();
+  fVariablesTreeSgn = new TTree(Form("%s_Sgn", nameoutput), "Candidates variables tree, Signal");
+  fVariablesTreeBkg = new TTree(Form("%s_Bkg", nameoutput), "Candidates variables tree, Background");
+  Int_t nVar = 66;
+  fCandidateVariables = new Float_t [nVar];
+  TString * fCandidateVariableNames = new TString[nVar];
+  fCandidateVariableNames[0]="massLc2K0Sp";
+  fCandidateVariableNames[1]="massLc2Lambdapi";
+  fCandidateVariableNames[2]="massK0S";
+  fCandidateVariableNames[3]="massLambda";
+  fCandidateVariableNames[4]="massLambdaBar";
+  fCandidateVariableNames[5]="cosPAK0S";
+  fCandidateVariableNames[6]="dcaV0";
+  fCandidateVariableNames[7]="tImpParBach";
+  fCandidateVariableNames[8]="tImpParV0";
+  fCandidateVariableNames[9]="nSigmaTPCpr";
+  fCandidateVariableNames[10]="nSigmaTPCpi";
+  fCandidateVariableNames[11]="nSigmaTPCka";
+  fCandidateVariableNames[12]="nSigmaTOFpr";
+  fCandidateVariableNames[13]="nSigmaTOFpi";
+  fCandidateVariableNames[14]="nSigmaTOFka";
+  fCandidateVariableNames[15]="bachelorPt";
+  fCandidateVariableNames[16]="V0positivePt";
+  fCandidateVariableNames[17]="V0negativePt";
+  fCandidateVariableNames[18]="dcaV0pos";
+  fCandidateVariableNames[19]="dcaV0neg";
+  fCandidateVariableNames[20]="v0Pt";
+  fCandidateVariableNames[21]="massGamma";
+  fCandidateVariableNames[22]="LcPt";
+  fCandidateVariableNames[23]="combinedProtonProb";
+  fCandidateVariableNames[24]="LcEta";
+  fCandidateVariableNames[25]="V0positiveEta";
+  fCandidateVariableNames[26]="V0negativeEta";
+  fCandidateVariableNames[27]="combinedPionProb";
+  fCandidateVariableNames[28]="combinedKaonProb";
+  fCandidateVariableNames[29]="bachelorEta";
+  fCandidateVariableNames[30]="LcP";
+  fCandidateVariableNames[31]="bachelorP";
+  fCandidateVariableNames[32]="v0P";
+  fCandidateVariableNames[33]="V0positiveP";
+  fCandidateVariableNames[34]="V0negativeP";
+  fCandidateVariableNames[35]="LcY";
+  fCandidateVariableNames[36]="v0Y";
+  fCandidateVariableNames[37]="bachelorY";
+  fCandidateVariableNames[38]="V0positiveY";
+  fCandidateVariableNames[39]="V0negativeY";
+  fCandidateVariableNames[40]="v0Eta";
+  fCandidateVariableNames[41]="DecayLengthLc";
+  fCandidateVariableNames[42]="DecayLengthK0S";
+  fCandidateVariableNames[43]="CtLc";
+  fCandidateVariableNames[44]="CtK0S";
+  fCandidateVariableNames[45]="bachCode";
+  fCandidateVariableNames[46]="k0SCode";
+
+  fCandidateVariableNames[47]="V0KFmass";
+  fCandidateVariableNames[48]="V0KFdecayLength";
+  fCandidateVariableNames[49]="V0KFlifeTime";
+
+  fCandidateVariableNames[50]="V0KFmassErr";
+  fCandidateVariableNames[51]="V0KFdecayTimeErr";
+  fCandidateVariableNames[52]="V0KFlifeTimeErr";
+
+  fCandidateVariableNames[53]="LcKFmass";
+  fCandidateVariableNames[54]="LcKFdecayLength";
+  fCandidateVariableNames[55]="LcKFlifeTime";
+
+  fCandidateVariableNames[56]="LcKFmassErr";
+  fCandidateVariableNames[57]="LcKFdecayTimeErr";
+  fCandidateVariableNames[58]="LcKFlifeTimeErr";
+
+  fCandidateVariableNames[59]="LcKFDistToPrimVtx";
+  fCandidateVariableNames[60]="V0KFDistToPrimVtx";
+  fCandidateVariableNames[61]="V0KFDistToLc";
+  fCandidateVariableNames[62]="alphaArmKF";
+  fCandidateVariableNames[63]="ptArmKF";
+  fCandidateVariableNames[64]="alphaArm";
+  fCandidateVariableNames[65]="ptArm";
+
+  for(Int_t ivar=0; ivar<nVar; ivar++){
+    fVariablesTreeSgn->Branch(fCandidateVariableNames[ivar].Data(),&fCandidateVariables[ivar],Form("%s/f",fCandidateVariableNames[ivar].Data()));
+    fVariablesTreeBkg->Branch(fCandidateVariableNames[ivar].Data(),&fCandidateVariables[ivar],Form("%s/f",fCandidateVariableNames[ivar].Data()));
+  }
+       
+  fHistoEvents = new TH1F("fHistoEvents", "fHistoEvents", 2, -0.5, 1.5);
+  TString labelEv[2] = {"NotSelected", "Selected"};
+  for (Int_t ibin = 1; ibin <= fHistoEvents->GetNbinsX(); ibin++){
+    fHistoEvents->GetXaxis()->SetBinLabel(ibin, labelEv[ibin-1].Data());
+  }  
+
+  fHistoLc = new TH1F("fHistoLc", "fHistoLc", 2, -0.5, 1.5);
+
+  fHistoLcOnTheFly = new TH1F("fHistoLcOnTheFly", "fHistoLcOnTheFly", 4, -0.5, 3.5);
+  TString labelOnTheFly[4] = {"OnTheFly-Bkg", "OfflineBkg", "OnTheFly-Sgn", "OfflineSgn"};
+  for (Int_t ibin = 1; ibin <= fHistoLcOnTheFly->GetNbinsX(); ibin++){
+    fHistoLcOnTheFly->GetXaxis()->SetBinLabel(ibin, labelOnTheFly[ibin-1].Data());
+  }
+
+  fHistoLcBeforeCuts = new TH1F("fHistoLcBeforeCuts", "fHistoLcBeforeCuts", 2, -0.5, 1.5);
+  TString labelBeforeCuts[2] = {"Bkg", "Sgn"};
+  for (Int_t ibin = 1; ibin <= fHistoLcBeforeCuts->GetNbinsX(); ibin++){
+    fHistoLcBeforeCuts->GetXaxis()->SetBinLabel(ibin, labelBeforeCuts[ibin-1].Data());
+    fHistoLc->GetXaxis()->SetBinLabel(ibin, labelBeforeCuts[ibin-1].Data());
+  }
+
+  fHistoFiducialAcceptance = new TH1F("fHistoFiducialAcceptance", "fHistoFiducialAcceptance", 4, -0.5, 3.5);
+  TString labelAcc[4] = {"NotOk-Bkg", "Ok-Bkg", "NotOk-Sgn", "Ok-Sgn"};
+  for (Int_t ibin = 1; ibin <= fHistoFiducialAcceptance->GetNbinsX(); ibin++){
+    fHistoFiducialAcceptance->GetXaxis()->SetBinLabel(ibin, labelAcc[ibin-1].Data());
+  }
+
+  fHistoCodesSgn = new TH2F("fHistoCodesSgn", "fHistoCodes for Signal; bachelor; K0S", 7, -1.5, 5.5, 9, -1.5, 7.5);
+  fHistoCodesBkg = new TH2F("fHistoCodesBkg", "fHistoCodes for Background; bachelor; K0S", 7, -1.5, 5.5, 9, -1.5, 7.5);
+
+  TString labelx[7] = { "kBachInvalid", "kBachFake", "kBachNoProton", "kBachPrimary", "kBachNoLambdaMother",
+                       "kBachDifferentLambdaMother",   "kBachCorrectLambdaMother"};
+  TString labely[9] = { "kK0SInvalid", "kK0SFake", "kK0SNoK0S", "kK0SWithoutMother", "kK0SNotFromK0", 
+                       "kK0Primary", "kK0NoLambdaMother", "kK0DifferentLambdaMother", "kK0CorrectLambdaMother"};  
+  
+  for (Int_t ibin = 1; ibin <= fHistoCodesSgn->GetNbinsX(); ibin++){
+    fHistoCodesSgn->GetXaxis()->SetBinLabel(ibin, labelx[ibin-1].Data());
+    fHistoCodesBkg->GetXaxis()->SetBinLabel(ibin, labelx[ibin-1].Data());
+  }
+  for (Int_t ibin = 1; ibin <= fHistoCodesSgn->GetNbinsY(); ibin++){
+    fHistoCodesSgn->GetYaxis()->SetBinLabel(ibin, labely[ibin-1].Data());
+    fHistoCodesBkg->GetYaxis()->SetBinLabel(ibin, labely[ibin-1].Data());
+  }
+
+  fHistoLcpKpiBeforeCuts = new TH1F("fHistoLcpKpiBeforeCuts", "fHistoLcpKpiBeforeCuts", 2, -0.5, 1.5);
+  for (Int_t ibin = 1; ibin <= fHistoLcpKpiBeforeCuts->GetNbinsX(); ibin++){
+    fHistoLcpKpiBeforeCuts->GetXaxis()->SetBinLabel(ibin, labelBeforeCuts[ibin-1].Data());
+  }
+
+  fHistoBackground = new TH1F("fHistoBackground", "fHistoBackground", 2, -0.5, 1.5);
+  TString labelBkg[2] = {"Injected", "Non-injected"};
+  for (Int_t ibin = 1; ibin <= fHistoBackground->GetNbinsX(); ibin++){
+    fHistoBackground->GetXaxis()->SetBinLabel(ibin, labelBkg[ibin-1].Data());
+  }
+
+  //fOutput->Add(fVariablesTreeSgn);
+  //fOutput->Add(fVariablesTreeBkg);
+
+  fOutput->Add(fHistoEvents);
+  fOutput->Add(fHistoLc);
+  fOutput->Add(fHistoLcOnTheFly);
+  fOutput->Add(fHistoLcBeforeCuts);
+  fOutput->Add(fHistoFiducialAcceptance);
+  fOutput->Add(fHistoCodesSgn);
+  fOutput->Add(fHistoCodesBkg);
+  fOutput->Add(fHistoLcpKpiBeforeCuts);
+  fOutput->Add(fHistoBackground);
+
+  PostData(1, fOutput);
+  PostData(4, fVariablesTreeSgn);
+  PostData(5, fVariablesTreeBkg);
+  AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
+  AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
+  fPIDResponse = inputHandler->GetPIDResponse();
+
+  if (fAnalCuts->GetIsUsePID()){
+    /*
+      fAnalCuts->GetPidHF()->SetPidResponse(fPIDResponse);
+      fAnalCuts->GetPidV0pos()->SetPidResponse(fPIDResponse);
+      fAnalCuts->GetPidV0neg()->SetPidResponse(fPIDResponse);
+      fAnalCuts->GetPidHF()->SetOldPid(kFALSE);
+      fAnalCuts->GetPidV0pos()->SetOldPid(kFALSE);
+      fAnalCuts->GetPidV0neg()->SetOldPid(kFALSE);
+    */
+    fAnalCuts->SetUsePID(kFALSE); // I don't want to use the PID through the cut object, but I will use the PID response directly!!!
+  }
+
+  // Setting properties of PID
+  fPIDCombined=new AliPIDCombined;
+  fPIDCombined->SetDefaultTPCPriors();
+  fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC+AliPIDResponse::kDetTOF);
+  //fPIDCombined->SetPriorDistribution((AliPID::EParticleType)ispec,fPriors[ispec]);
+
+  //Printf("Here ok!");
+
+  fCounter = new AliNormalizationCounter("NormalizationCounter");
+  fCounter->Init();
+  PostData(2, fCounter);
+
+  // Histograms from KF
+
+  fHistoDistanceLcToPrimVtx    = new TH1D("fHistoDistanceLcToPrimVtx", "Lc distance to Prim Vertex from KF; distance [cm]", 1000, 0., 1);
+  fHistoDistanceV0ToPrimVtx    = new TH1D("fHistoDistanceV0ToPrimVtx", "V0 distance to Prim Vertex from KF; distance [cm]", 1000, 0., 100.);
+  fHistoDistanceV0ToLc         = new TH1D("fHistoDistanceV0ToLc", "V0 distance to Lc Vertex from KF; distance [cm]", 1000, 0., 100.);
+
+  fHistoDistanceLcToPrimVtxSgn = new TH1D("fHistoDistanceLcToPrimVtxSgn", "Lc Sgn distance to Prim Vertex from KF; distance [cm]", 1000, 0., 1);
+  fHistoDistanceV0ToPrimVtxSgn = new TH1D("fHistoDistanceV0ToPrimVtxSgn", "V0 Sgn distance to Prim Vertex from KF; distance [cm]", 1000, 0., 100.);
+  fHistoDistanceV0ToLcSgn      = new TH1D("fHistoDistanceV0ToLcSgn", "V0 Sgn distance to Lc Vertex from KF; distance [cm]", 1000, 0., 100.);
+
+  fHistoVtxLcResidualToPrimVtx = new TH1D("fHistoVtxLcResidualToPrimVtx", "Residual between MC and KF (MC - KF): Lc to Prim Vtx; distance [cm]", 1000, -5., 5.);
+  fHistoVtxV0ResidualToPrimVtx = new TH1D("fHistoVtxV0ResidualToPrimVtx", "Residual between MC and KF (MC - KF): V0 to Prim Vtx; distance [cm]", 1000, -5., 5.);
+  fHistoVtxV0ResidualToLc      = new TH1D("fHistoVtxV0ResidualToLc", "Residual between MC and KF: V0 to Lc (MC - KF); distance [cm]", 1000, -5., 5.);
+
+  fHistoMassV0All              = new TH1D("fHistoMassV0All", "V0 Mass; mass", 500, 0.4, 0.6);
+  fHistoDecayLengthV0All       = new TH1D("fHistoDecayLengthV0All", "V0 Decay Length; decayLength", 500, -10, 10.0);
+  fHistoLifeTimeV0All          = new TH1D("fHistoLifeTimeV0All", "V0 Life Time; lifeTime", 500, -10.0, 10.0);
+
+  fHistoMassV0True             = new TH1D("fHistoMassV0True", "True V0 Mass; mass", 500, 0.4, 0.6);
+  fHistoDecayLengthV0True      = new TH1D("fHistoDecayLengthV0True", "True V0 Decay Length; decayLength", 500, -10, 10.0);
+  fHistoLifeTimeV0True         = new TH1D("fHistoLifeTimeV0True", "True V0 Life Time; lifeTime", 500, -10.0, 10.0);
+
+  fHistoMassV0TrueFromAOD      = new TH1D("fHistoMassV0TrueFormAOD", "True V0 Mass (AOD); mass", 500, 0.4, 0.6);
+
+  fHistoMassV0TrueK0S          = new TH1D("fHistoMassV0TrueK0S", "True V0-K0S Mass; mass", 500, 0.4, 0.6);
+  fHistoDecayLengthV0TrueK0S   = new TH1D("fHistoDecayLengthV0TrueK0S", "True V0-K0S Decay Length; decayLength", 500, -10, 10.0);
+  fHistoLifeTimeV0TrueK0S      = new TH1D("fHistoLifeTimeV0TrueK0S", "True V0-K0S Life Time; lifeTime", 500, -10.0, 10.0);
+
+  fHistoMassV0TrueK0SFromAOD   = new TH1D("fHistoMassV0TrueK0SFormAOD", "True V0-K0S Mass (AOD); mass", 500, 0.4, 0.6);
+
+  fHistoMassLcAll              = new TH1D("fHistoMassLcAll", "Lc Mass; mass", 500, 2.0, 3.0);
+  fHistoDecayLengthLcAll       = new TH1D("fHistoDecayLengthLcAll", "Lc Decay Length; decayLenght", 100000, -0.1, 0.1);
+  fHistoLifeTimeLcAll          = new TH1D("fHistoLifeTimeLcAll", "Lc Life Time; lifeTime", 100000, -0.1, 0.1);
+
+  fHistoMassLcTrue             = new TH1D("fHistoMassLcTrue", "True Lc Mass; mass", 500, 2.0, 3.0);
+  fHistoDecayLengthLcTrue      = new TH1D("fHistoDecayLengthLcTrue", "True Lc Decay Length; decayLength", 100000, -0.1, 0.1);
+  fHistoLifeTimeLcTrue         = new TH1D("fHistoLifeTimeLcTrue", "True Lc Life Time; lifeTime", 100000, -0.1, 0.1);
+
+  fHistoMassLcTrueFromAOD      = new TH1D("fHistoMassLcTrueFromAOD", "True Lc Mass (AOD); mass", 500, 2.0, 3.0);
+
+  fHistoMassV0fromLcAll        = new TH1D("fHistoMassV0fromLcAll", "V0 mass from Lc built in KF; mass", 500, 0.4, 0.6);
+  fHistoDecayLengthV0fromLcAll = new TH1D("fHistoDecayLengthV0fromLcAll", "V0 Decay Length from Lc built in KF; decayLength", 500, 0, 10.0);
+  fHistoLifeTimeV0fromLcAll    = new TH1D("fHistoLifeTimeV0fromLcAll", "V0 Life Time from Lc built in KF; lifeTime", 500, 0.0, 3.0);
+
+  fHistoMassV0fromLcTrue       = new TH1D("fHistoMassV0fromLcTrue", "V0 mass from true Lc built in KF; mass", 500, 0.4, 0.6);
+  fHistoDecayLengthV0fromLcTrue= new TH1D("fHistoDecayLengthV0fromLcTrue", "V0 Decay Length from true Lc built in KF; decayLength", 500, 0, 10.0);
+  fHistoLifeTimeV0fromLcTrue   = new TH1D("fHistoLifeTimeV0fromLcTrue", "V0 Life Time from true Lc built in KF; lifeTime", 500, 0.0, 3.0);
+
+  fHistoMassLcSgn              = new TH1D("fHistoMassLcSgn", "True Lc Signal Mass; mass", 500, 2.0, 3.0);
+  fHistoMassLcSgnFromAOD       = new TH1D("fHistoMassLcSgnFromAOD", "True Lc Signal Mass (AOD); mass", 500, 2.0, 3.0);
+  fHistoDecayLengthLcSgn       = new TH1D("fHistoDecayLengthLcSgn", "True Lc Signal Decay Length; decayLength", 100000, -0.1, 0.1);
+  fHistoLifeTimeLcSgn          = new TH1D("fHistoLifeTimeLcSgn", "True Lc Signal Life Time; lifeTime", 100000, -0.1, 0.1);
+
+  fHistoMassV0fromLcSgn        = new TH1D("fHistoMassV0fromLcSgn", "V0 from True Lc Signal Mass; mass", 500, 0.4, 0.6);
+  fHistoDecayLengthV0fromLcSgn = new TH1D("fHistoDecayLengthV0fromLcSgn", "V0 True Lc Signal Decay Length; decayLength", 500, 0, 10.0);
+  fHistoLifeTimeV0fromLcSgn    = new TH1D("fHistoLifeTimeV0fromLcSgn", "V0 True Lc Signal Life Time; lifeTime", 500, 0.0, 3.0);
+
+  fHistoKF = new TH2D("fHistoKF", "Summary from KF; V0 KF; Lc KF", 16, -0.5, 15.5, 16, -0.5, 15.5);
+  fHistoKFV0 = new TH1D("fHistoKFV0", "Summary from KF; V0 KF", 16, -0.5, 15.5);
+  fHistoKFLc = new TH1D("fHistoKFLc", "Summary from KF; V0 KF", 16, -0.5, 15.5);
+  TString axisLabel[16] = {"AllOk", "M_NotOk", "Sm_NotOk", "Dl_NotOk", 
+                          "Lt_NotOk", "M_Sm_NotOk", "M_Dl_NotOk", "M_Lt_NotOk", 
+                          "Dl_Sm_NotOk", "Dl_Lt_NotOk", "Sm_Lt_NotOk", "M_Sm_Dl_NotOk", 
+                          "M_Sm_Lt_NotOk", "Sm_Dl_Lt_NotOk", "M_Dl_Lt_NotOk", "All_NotOk"};
+
+  for (Int_t ibin = 1; ibin <=16; ibin++){
+    fHistoKF->GetXaxis()->SetBinLabel(ibin, axisLabel[ibin-1].Data());
+    fHistoKF->GetYaxis()->SetBinLabel(ibin, axisLabel[ibin-1].Data());
+    fHistoKFV0->GetXaxis()->SetBinLabel(ibin, axisLabel[ibin-1].Data());
+    fHistoKFLc->GetXaxis()->SetBinLabel(ibin, axisLabel[ibin-1].Data());
+  }
+
+  fHistoMassKFV0 = new TH2D("fHistoMassKFV0", "mass vs sigmaMass for V0; mass; sigmaMass", 500, 0.4, 0.6, 500, 0., 10);
+  fHistoDecayLengthKFV0 = new TH2D("fHistoDecayLengthKFV0", "decayLength vs sigmaDecayLength for V0; decayLength; sigmaDecayLength", 500, -10, 10, 500, 0., 10);
+  fHistoLifeTimeKFV0 = new TH2D("fHistoLifeTimeKFV0", "lifeTime vs sigmalifeTime for V0; lifeTime; sigmaLifeTime", 500, -10, 10, 500, 0., 10);
+
+  fHistoMassKFLc = new TH2D("fHistoMassKFLc", "mass vs sigmaMass for Lc; mass; sigmaMass", 500, 0.4, 0.6, 500, 0., 10);
+  fHistoDecayLengthKFLc = new TH2D("fHistoDecayLengthKFLc", "decayLength vs sigmaDecayLength for Lc; decayLength; sigmaDecayLength", 500, -10, 10, 500, 0., 10);
+  fHistoLifeTimeKFLc = new TH2D("fHistoLifeTimeKFLc", "lifeTime vs sigmalifeTime for Lc; lifeTime; sigmaLifeTime", 500, -10, 10, 500, 0., 10);
+
+  fHistoArmenterosPodolanskiV0KF = new TH2D("fHistoArmenterosPodolanskiV0KF", "V0 ArmenterosPodolanski from KF; #alpha; Qt", 1000, -1, 1, 1000, 0, 1); 
+  fHistoArmenterosPodolanskiV0KFSgn = new TH2D("fHistoArmenterosPodolanskiV0KFSgn", "V0 (signal) ArmenterosPodolanski from KF; #alpha; Qt", 1000, -1, 1, 1000, 0, 1); 
+  fHistoArmenterosPodolanskiV0AOD = new TH2D("fHistoArmenterosPodolanskiV0AOD", "V0 ArmenterosPodolanski from AOD; #alpha; Qt", 1000, -1, 1, 1000, 0, 1); 
+  fHistoArmenterosPodolanskiV0AODSgn = new TH2D("fHistoArmenterosPodolanskiV0AODSgn", "V0 (signal) ArmenterosPodolanski from AOD; #alpha; Qt", 1000, -1, 1, 1000, 0, 1); 
+
+  fOutputKF = new TList();
+  fOutputKF->SetOwner();
+  fOutputKF->SetName("listHistoKF");
+
+  fOutputKF->Add(fHistoDistanceLcToPrimVtx);
+  fOutputKF->Add(fHistoDistanceV0ToPrimVtx);
+  fOutputKF->Add(fHistoDistanceV0ToLc);
+
+  fOutputKF->Add(fHistoDistanceLcToPrimVtxSgn);
+  fOutputKF->Add(fHistoDistanceV0ToPrimVtxSgn);
+  fOutputKF->Add(fHistoDistanceV0ToLcSgn);
+
+  fOutputKF->Add(fHistoVtxLcResidualToPrimVtx);
+  fOutputKF->Add(fHistoVtxV0ResidualToPrimVtx);
+  fOutputKF->Add(fHistoVtxV0ResidualToLc);
+
+  fOutputKF->Add(fHistoMassV0All);
+  fOutputKF->Add(fHistoDecayLengthV0All);
+  fOutputKF->Add(fHistoLifeTimeV0All);
+
+  fOutputKF->Add(fHistoMassV0True);
+  fOutputKF->Add(fHistoDecayLengthV0True);
+  fOutputKF->Add(fHistoLifeTimeV0True);
+
+  fOutputKF->Add(fHistoMassV0TrueFromAOD);
+
+  fOutputKF->Add(fHistoMassV0TrueK0S);
+  fOutputKF->Add(fHistoDecayLengthV0TrueK0S);
+  fOutputKF->Add(fHistoLifeTimeV0TrueK0S);
+
+  fOutputKF->Add(fHistoMassV0TrueK0SFromAOD);
+
+  fOutputKF->Add(fHistoMassLcAll);
+  fOutputKF->Add(fHistoDecayLengthLcAll);
+  fOutputKF->Add(fHistoLifeTimeLcAll);
+
+  fOutputKF->Add(fHistoMassLcTrue);
+  fOutputKF->Add(fHistoDecayLengthLcTrue);
+  fOutputKF->Add(fHistoLifeTimeLcTrue);
+
+  fOutputKF->Add(fHistoMassLcTrueFromAOD);
+
+  fOutputKF->Add(fHistoMassV0fromLcAll);
+  fOutputKF->Add(fHistoDecayLengthV0fromLcAll);
+  fOutputKF->Add(fHistoLifeTimeV0fromLcAll);
+
+  fOutputKF->Add(fHistoMassV0fromLcTrue);
+  fOutputKF->Add(fHistoDecayLengthV0fromLcTrue);
+  fOutputKF->Add(fHistoLifeTimeV0fromLcTrue);
+
+  fOutputKF->Add(fHistoMassLcSgn);
+  fOutputKF->Add(fHistoDecayLengthLcSgn);
+  fOutputKF->Add(fHistoLifeTimeLcSgn);
+
+  fOutputKF->Add(fHistoMassV0fromLcSgn);
+  fOutputKF->Add(fHistoDecayLengthV0fromLcSgn);
+  fOutputKF->Add(fHistoLifeTimeV0fromLcSgn);
+  
+  fOutputKF->Add(fHistoKF);
+  fOutputKF->Add(fHistoKFV0);
+  fOutputKF->Add(fHistoKFLc);
+
+  fOutputKF->Add(fHistoMassKFV0);
+  fOutputKF->Add(fHistoDecayLengthKFV0);
+  fOutputKF->Add(fHistoLifeTimeKFV0);
+
+  fOutputKF->Add(fHistoMassKFLc);
+  fOutputKF->Add(fHistoDecayLengthKFLc);
+  fOutputKF->Add(fHistoLifeTimeKFLc);
+
+  fOutputKF->Add(fHistoArmenterosPodolanskiV0KF);
+  fOutputKF->Add(fHistoArmenterosPodolanskiV0KFSgn);
+  fOutputKF->Add(fHistoArmenterosPodolanskiV0AOD);
+  fOutputKF->Add(fHistoArmenterosPodolanskiV0AODSgn);
+
+  PostData(6, fOutputKF);
+
+  //Printf("Here ok 1!");
+
+  return;
+}
+
+//_________________________________________________
+void AliAnalysisTaskSELc2V0bachelorTMVA::UserExec(Option_t *)
+{
+  // user exec
+  if (!fInputEvent) {
+    AliError("NO EVENT FOUND!");
+    return;
+  }
+
+  AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
+  TClonesArray *arrayLctopKos=0;
+
+  TClonesArray *array3Prong = 0;
+
+  if (!aodEvent && AODEvent() && IsStandardAOD()) {
+    // In case there is an AOD handler writing a standard AOD, use the AOD 
+    // event in memory rather than the input (ESD) event.    
+    aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
+    // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
+    // have to taken from the AOD event hold by the AliAODExtension
+    AliAODHandler* aodHandler = (AliAODHandler*) 
+      ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
+   
+    if (aodHandler->GetExtensions()) {
+      AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
+      AliAODEvent *aodFromExt = ext->GetAOD();
+      arrayLctopKos=(TClonesArray*)aodFromExt->GetList()->FindObject("CascadesHF");
+
+      array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
+    }
+  } else {
+    arrayLctopKos=(TClonesArray*)aodEvent->GetList()->FindObject("CascadesHF");
+
+    array3Prong=(TClonesArray*)aodEvent->GetList()->FindObject("Charm3Prong");
+  }
+  
+  if ( !fUseMCInfo && fIspA) {
+    fAnalCuts->SetTriggerClass("");
+    fAnalCuts->SetTriggerMask(AliVEvent::kINT7);
+  }
+
+  fCounter->StoreEvent(aodEvent,fAnalCuts,fUseMCInfo);
+
+  // AOD primary vertex
+  fVtx1 = (AliAODVertex*)aodEvent->GetPrimaryVertex();
+  if (!fVtx1) return;
+  if (fVtx1->GetNContributors()<1) return;
+
+  fIsEventSelected = fAnalCuts->IsEventSelected(aodEvent);
+
+  if ( !fIsEventSelected ) {
+    fHistoEvents->Fill(0);
+    return; // don't take into account not selected events 
+  }
+  fHistoEvents->Fill(1);
+
+  // Setting magnetic field for KF vertexing
+  AliKFParticle::SetField(aodEvent->GetMagneticField());
+
+  // mc analysis 
+  TClonesArray *mcArray = 0;
+  AliAODMCHeader *mcHeader=0;
+
+  if (fUseMCInfo) {
+    // MC array need for maching
+    mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
+    if (!mcArray) {
+      AliError("Could not find Monte-Carlo in AOD");
+      return;
+    }
+    // load MC header
+    mcHeader = (AliAODMCHeader*)aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName());
+    if (!mcHeader) {
+      AliError("AliAnalysisTaskSELc2V0bachelorTMVA::UserExec: MC header branch not found!\n");
+      return;
+    }
+  }
+
+  Int_t nSelectedAnal = 0;
+  if (fIsK0sAnalysis) {
+    MakeAnalysisForLc2prK0S(arrayLctopKos, mcArray,
+                           nSelectedAnal, fAnalCuts, 
+                           array3Prong, mcHeader);
+  }
+  fCounter->StoreCandidates(aodEvent,nSelectedAnal,kFALSE);
+
+  PostData(1, fOutput);
+  PostData(2, fCounter);
+  PostData(4, fVariablesTreeSgn);
+  PostData(5, fVariablesTreeBkg);
+  PostData(6, fOutputKF);
+
+}
+
+//-------------------------------------------------------------------------------
+void AliAnalysisTaskSELc2V0bachelorTMVA::MakeAnalysisForLc2prK0S(TClonesArray *arrayLctopKos,
+                                                                TClonesArray *mcArray,
+                                                                Int_t &nSelectedAnal,
+                                                                AliRDHFCutsLctoV0 *cutsAnal, TClonesArray *array3Prong,
+                                                                AliAODMCHeader* aodheader){
+  //Lc prong needed to MatchToMC method
+
+  Int_t pdgCand = 4122;
+  Int_t pdgDgLctoV0bachelor[2]={2212, 310};
+  Int_t pdgDgV0toDaughters[2]={211, 211};
+  
+  Int_t pdgDgLctopKpi[3]={2212, 321, 211};
+
+  // loop to search for candidates Lc->p+K+pi
+  Int_t n3Prong = array3Prong->GetEntriesFast();
+  Int_t nCascades= arrayLctopKos->GetEntriesFast();
+  
+  //AliInfo(Form("\n\n\n\n3 prong candidates = %d, ncascades = %d \n\n\n\n\n", n3Prong, nCascades));
+  for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
+    AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
+    //Filling a control histogram with no cuts
+    if (fUseMCInfo) {
+      
+      // find associated MC particle for Lc -> p+K+pi
+      Int_t mcLabel = d->MatchToMC(4122, mcArray, 3, pdgDgLctopKpi);
+      //Int_t mcLabelTemp = d->MatchToMC(4122, mcArray);
+      //Printf("mcLabel = %d, mcLabelTemp = %d", mcLabel, mcLabelTemp);
+      if (mcLabel >= 0) {
+       
+       AliAODMCParticle *partLcpKpi = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcLabel));
+       if(partLcpKpi){
+         Int_t pdgCode = partLcpKpi->GetPdgCode();
+         AliDebug(2,Form(" ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ MClabel=%d ~~~~~~~~~~ pdgCode=%d", mcLabel, pdgCode));
+         //AliInfo(Form(" ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ MClabel=%d ~~~~~~~~~~ pdgCode=%d", mcLabel, pdgCode));
+         fHistoLcpKpiBeforeCuts->Fill(1);
+                                       
+       }
+      } 
+      else {
+       //AliInfo(Form(" ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ MClabel=%d ~~~~~~~~~~", mcLabel));
+       fHistoLcpKpiBeforeCuts->Fill(0);
+      }
+    }  
+  }
+  
+  // loop over cascades to search for candidates Lc->p+K0S
+
+  for (Int_t iLctopK0s = 0; iLctopK0s < nCascades; iLctopK0s++) {
+
+    // Lc candidates and K0s from Lc
+    AliAODRecoCascadeHF* lcK0spr = dynamic_cast<AliAODRecoCascadeHF*>(arrayLctopKos->At(iLctopK0s));
+    if (!lcK0spr) {
+      AliDebug(2,Form("Cascade %d doens't exist, skipping",iLctopK0s));
+      continue;
+    }
+
+    //Filling a control histogram with no cuts
+    if (fUseMCInfo) {
+
+      Int_t pdgCode=-2;
+
+      // find associated MC particle for Lc -> p+K0 and K0S->pi+pi
+      fmcLabelLc = lcK0spr->MatchToMC(pdgCand, pdgDgLctoV0bachelor[1], pdgDgLctoV0bachelor, pdgDgV0toDaughters, mcArray, kTRUE);
+      if (fmcLabelLc>=0) {
+       AliDebug(2,Form(" ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~cascade number %d (total cascade number = %d)", iLctopK0s, nCascades));
+
+       AliAODMCParticle *partLc = dynamic_cast<AliAODMCParticle*>(mcArray->At(fmcLabelLc));
+       if(partLc){
+         pdgCode = partLc->GetPdgCode();
+         if (pdgCode<0) AliDebug(2,Form(" ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ MClabel=%d ~~~~~~~~~~ pdgCode=%d", fmcLabelLc, pdgCode));
+         pdgCode = TMath::Abs(pdgCode);
+         fHistoLcBeforeCuts->Fill(1);
+                                       
+       }
+      } 
+      else {
+       fHistoLcBeforeCuts->Fill(0);
+       pdgCode=-1;
+      }
+    }
+
+    //if (!lcK0spr->GetSecondaryVtx()) {
+    // AliInfo("No secondary vertex");
+    //continue;
+    //}
+
+    if (lcK0spr->GetNDaughters()!=2) {
+      AliDebug(2,Form("Cascade %d has not 2 daughters (nDaughters=%d)",iLctopK0s,lcK0spr->GetNDaughters()));
+      continue;
+    }
+
+    AliAODv0 * v0part = dynamic_cast<AliAODv0*>(lcK0spr->Getv0());
+    AliAODTrack * bachPart = dynamic_cast<AliAODTrack*>(lcK0spr->GetBachelor());
+    if (!v0part || !bachPart) {
+      AliDebug(2,Form("Cascade %d has no V0 or no bachelor object",iLctopK0s));
+      continue;
+    }
+
+
+    if (!v0part->GetSecondaryVtx()) {
+      AliDebug(2,Form("No secondary vertex for V0 by cascade %d",iLctopK0s));
+      continue;
+    }
+
+    if (v0part->GetNDaughters()!=2) {
+      AliDebug(2,Form("current V0 has not 2 daughters (onTheFly=%d, nDaughters=%d)",v0part->GetOnFlyStatus(),v0part->GetNDaughters()));
+      continue;
+    }
+
+    AliAODTrack * v0Pos = dynamic_cast<AliAODTrack*>(lcK0spr->Getv0PositiveTrack());
+    AliAODTrack * v0Neg = dynamic_cast<AliAODTrack*>(lcK0spr->Getv0NegativeTrack());
+    if (!v0Neg || !v0Neg) {
+      AliDebug(2,Form("V0 by cascade %d has no V0positive of V0negative object",iLctopK0s));
+      continue;
+    }
+
+
+    if (v0Pos->Charge() == v0Neg->Charge()) {
+      AliDebug(2,Form("V0 by cascade %d has daughters with the same sign: IMPOSSIBLE!",iLctopK0s));
+      continue;
+    }
+
+    Int_t isLc = 0;
+
+    if (fUseMCInfo) {
+
+      Int_t pdgCode=-2;
+
+      // find associated MC particle for Lc -> p+K0 and K0S->pi+pi
+      Int_t mcLabel = lcK0spr->MatchToMC(pdgCand, pdgDgLctoV0bachelor[1], pdgDgLctoV0bachelor, pdgDgV0toDaughters, mcArray, kTRUE);
+      if (mcLabel>=0) {
+       AliDebug(2,Form(" ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~cascade number %d (total cascade number = %d)", iLctopK0s, nCascades));
+
+       AliAODMCParticle *partLc = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcLabel));
+       if(partLc){
+         pdgCode = partLc->GetPdgCode();
+         if (pdgCode<0) AliDebug(2,Form(" ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ MClabel=%d ~~~~~~~~~~ pdgCode=%d", mcLabel, pdgCode));
+         pdgCode = TMath::Abs(pdgCode);
+         isLc = 1;
+         fHistoLc->Fill(1);
+                                       
+       }
+      } 
+      else {
+       fHistoLc->Fill(0);
+       pdgCode=-1;
+      }
+    }
+    AliDebug(2, Form("\n\n\n Analysing candidate %d\n", iLctopK0s));
+    AliDebug(2, Form(">>>>>>>>>> Candidate is background, fFillOnlySgn = %d --> SKIPPING", fFillOnlySgn));
+    if (!isLc) {
+      if (fFillOnlySgn) { // if it is background, and we want only signal, we do not fill the tree
+       continue;
+      }
+      else { // checking if we want to fill the background
+       if (fKeepingOnlyHIJINGBkg){
+         // we have decided to fill the background only when the candidate has the daugthers that all come from HIJING underlying event!
+         //Bool_t isCandidateInjected = fUtils->HasCascadeCandidateAnyDaughInjected(part, aodheader, mcArray);
+         Bool_t isCandidateInjected = CheckInjection(lcK0spr, aodheader, mcArray);
+         if (!isCandidateInjected){
+           AliDebug(2, "The candidate is from HIJING (i.e. not injected), keeping it to fill background");
+           fHistoBackground->Fill(1);
+         }
+         else {
+           AliDebug(2, "The candidate is NOT from HIJING, we skip it when filling background");
+           fHistoBackground->Fill(0);
+           continue; 
+         }
+       }
+      }
+    }
+
+    FillLc2pK0Sspectrum(lcK0spr, isLc, nSelectedAnal, cutsAnal, mcArray);
+    
+  }
+  
+  return;
+
+}
+//________________________________________________________________________
+void AliAnalysisTaskSELc2V0bachelorTMVA::FillLc2pK0Sspectrum(AliAODRecoCascadeHF *part,
+                                                            Int_t isLc,
+                                                            Int_t &nSelectedAnal,
+                                                            AliRDHFCutsLctoV0 *cutsAnal,
+                                                            TClonesArray *mcArray){
+  //
+  // Fill histos for Lc -> K0S+proton
+  //
+
+  /*
+    if (!part->GetOwnPrimaryVtx()) {
+    //Printf("No primary vertex for Lc found!!");
+    part->SetOwnPrimaryVtx(fVtx1);
+    }
+    else {
+    //Printf("Yu-huuuu!!! primary vertex for Lc found!!");
+    }
+  */
+  Double_t invmassLc = part->InvMassLctoK0sP();
+
+  AliAODv0 * v0part = part->Getv0();
+  Bool_t onFlyV0 = v0part->GetOnFlyStatus(); // on-the-flight V0s
+  if (onFlyV0){ // on-the-fly V0
+    if (isLc) { // Lc
+      fHistoLcOnTheFly->Fill(2.);
+    }
+    else { // not Lc
+      fHistoLcOnTheFly->Fill(0.);
+    }
+  }
+  else { // offline V0
+    if (isLc) { // Lc
+      fHistoLcOnTheFly->Fill(3.);
+    }
+    else { // not Lc
+      fHistoLcOnTheFly->Fill(1.);
+    }
+  }
+  
+  Double_t dcaV0 = v0part->GetDCA();
+  Double_t invmassK0s = v0part->MassK0Short();
+
+  if ( (cutsAnal->IsInFiducialAcceptance(part->Pt(),part->Y(4122))) ) {
+    if (isLc) {
+      fHistoFiducialAcceptance->Fill(3.);
+    }
+    else {
+      fHistoFiducialAcceptance->Fill(1.);
+    }
+  }
+  else {
+    if (isLc) {
+      fHistoFiducialAcceptance->Fill(2.);
+    }
+    else {
+      fHistoFiducialAcceptance->Fill(0.);
+    }
+  }
+
+  Int_t isInV0window = (((cutsAnal->IsSelectedSingleCut(part, AliRDHFCuts::kCandidate, 2)) & (AliRDHFCutsLctoV0::kLcToK0Spr)) == (AliRDHFCutsLctoV0::kLcToK0Spr)); // cut on V0 invMass
+
+  if (isInV0window == 0) {
+    AliDebug(2, "No: The candidate has NOT passed the mass cuts!");
+    if (isLc) Printf("SIGNAL candidate rejected");
+    return;           
+  }
+  else AliDebug(2, "Yes: The candidate has passed the mass cuts!");
+
+  //   Printf("ciccio!!!!!!!");
+
+  Bool_t isInCascadeWindow = (((cutsAnal->IsSelectedSingleCut(part, AliRDHFCuts::kCandidate, 0)) & (AliRDHFCutsLctoV0::kLcToK0Spr)) == (AliRDHFCutsLctoV0::kLcToK0Spr)); // cut on Lc->p+K0S invMass
+
+  if (!isInCascadeWindow) {
+    AliDebug(2, "No: The candidate has NOT passed the cascade window cuts!");
+    if (isLc) Printf("SIGNAL candidate rejected");
+    return;
+  }
+  else AliDebug(2, "Yes: The candidate has passed the cascade window cuts!");
+
+  Bool_t isCandidateSelectedCuts = (((cutsAnal->IsSelected(part, AliRDHFCuts::kCandidate)) & (AliRDHFCutsLctoV0::kLcToK0Spr)) == (AliRDHFCutsLctoV0::kLcToK0Spr)); // kinematic/topological cuts
+  if (!isCandidateSelectedCuts){
+    AliDebug(2, "No: Analysis cuts kCandidate level NOT passed");
+    if (isLc) Printf("SIGNAL candidate rejected");
+    return;
+  }
+  else {
+    AliDebug(2, "Yes: Analysis cuts kCandidate level passed");
+  }
+
+  AliAODTrack *bachelor = (AliAODTrack*)part->GetBachelor();
+  if (!bachelor) {
+    AliDebug(2, Form("Very weird, the bachelor is not there... returning for this candidate"));
+    return;
+  }
+
+  //Bool_t isBachelorID = (((cutsAnal->IsSelected(part,AliRDHFCuts::kPID))&(AliRDHFCutsLctoV0::kLcToK0Spr))==(AliRDHFCutsLctoV0::kLcToK0Spr)); // ID x bachelor
+  Double_t probTPCTOF[AliPID::kSPECIES]={0.};
+       
+  UInt_t detUsed = fPIDCombined->ComputeProbabilities(bachelor, fPIDResponse, probTPCTOF);
+  Double_t probProton = 0.;
+  Double_t probPion = 0.;
+  Double_t probKaon = 0.;
+  if (detUsed == (UInt_t)fPIDCombined->GetDetectorMask() ) {
+    probProton = probTPCTOF[AliPID::kProton];
+    probPion = probTPCTOF[AliPID::kPion];
+    probKaon = probTPCTOF[AliPID::kKaon];
+  }
+
+  if ( !( !onFlyV0 || (onFlyV0 && fUseOnTheFlyV0) ) ) {
+    AliDebug(2, "On-the-fly discarded");
+    return;            
+  }
+
+  if ( (((cutsAnal->IsSelected(part,AliRDHFCuts::kAll))&(AliRDHFCutsLctoV0::kLcToK0Spr))==(AliRDHFCutsLctoV0::kLcToK0Spr)) ) {
+    nSelectedAnal++;
+  }
+
+  if ( !(cutsAnal->IsInFiducialAcceptance(part->Pt(),part->Y(4122))) ) return;
+
+  if ( !( ( (cutsAnal->IsSelected(part, AliRDHFCuts::kTracks)) & (AliRDHFCutsLctoV0::kLcToK0Spr)) == (AliRDHFCutsLctoV0::kLcToK0Spr) ) ) { // esd track cuts
+    if (isLc) Printf("SIGNAL candidate rejected");
+    AliDebug(2, "No: Analysis cuts kTracks level NOT passed");
+    return;
+  }
+  else {
+    AliDebug(2, "Yes: Analysis cuts kTracks level passed");
+  }
+
+  Int_t pdgCand = 4122;
+  Int_t pdgDgLctoV0bachelor[2]={211, 3122}; // case of wrong decay! Lc --> L + pi
+  Int_t pdgDgV0toDaughters[2]={2212, 211}; // case of wrong decay! Lc --> L + pi
+  Int_t isLc2LBarpi=0, isLc2Lpi=0;
+  Int_t currentLabel = part->GetLabel();
+  Int_t mcLabel = 0;
+  if (fUseMCInfo) {
+    mcLabel = part->MatchToMC(pdgCand, pdgDgLctoV0bachelor[1], pdgDgLctoV0bachelor, pdgDgV0toDaughters, mcArray, kTRUE);
+    if (mcLabel>=0) {
+      if (bachelor->Charge()==-1) isLc2LBarpi=1;
+      if (bachelor->Charge()==+1) isLc2Lpi=1;
+    }
+  }
+
+  Int_t pdgDg2prong[2] = {211, 211};
+  Int_t labelK0S = 0;
+  Int_t isK0S = 0;
+  if (fUseMCInfo) {
+    labelK0S = v0part->MatchToMC(310, mcArray, 2, pdgDg2prong);
+    if (labelK0S>=0) isK0S = 1;
+  }
+
+  pdgDg2prong[0] = 211;
+  pdgDg2prong[1] = 2212;
+  Int_t isLambda = 0;
+  Int_t isLambdaBar = 0;
+  Int_t lambdaLabel = 0;
+  if (fUseMCInfo) {
+    lambdaLabel = v0part->MatchToMC(3122, mcArray, 2, pdgDg2prong);
+    if (lambdaLabel>=0) {
+      AliAODMCParticle *lambdaTrack = (AliAODMCParticle*)mcArray->At(lambdaLabel);
+      if (lambdaTrack->GetPdgCode()==3122) isLambda = 1;
+      else if (lambdaTrack->GetPdgCode()==-3122) isLambdaBar = 1;
+    }
+  }
+
+  pdgDg2prong[0] = 11;
+  pdgDg2prong[1] = 11;
+  Int_t isGamma = 0;
+  Int_t gammaLabel = 0;
+  if (fUseMCInfo) {
+    gammaLabel = v0part->MatchToMC(22, mcArray, 2, pdgDg2prong);
+    if (gammaLabel>=0) {
+      AliAODMCParticle *gammaTrack = (AliAODMCParticle*)mcArray->At(gammaLabel);
+      if (gammaTrack->GetPdgCode()==22) isGamma = 1;
+    }
+  }
+               
+  Int_t pdgTemp = -1;
+  if (currentLabel != -1){
+    AliAODMCParticle *tempPart = (AliAODMCParticle*)mcArray->At(currentLabel);
+    pdgTemp = tempPart->GetPdgCode();
+  }
+  if (isLc) AliDebug(2, Form("Signal: Candidate is a Lc in K0s+p"));
+  else if (isLc2LBarpi) AliDebug(2, Form("Background: Candidate is a Lc in Lbar + pi"));
+  else if (isLc2Lpi) AliDebug(2, Form("Background: Candidate is a Lc in L + pi"));
+  else AliDebug(2, Form("Pure bkg: Candidate has pdg = %d", pdgTemp));
+  if (isK0S) AliDebug(2, Form("V0 is a K0S"));
+  else if (isLambda)  AliDebug(2, Form("V0 is a Lambda"));
+  else if (isLambdaBar)  AliDebug(2, Form("V0 is a LambdaBar"));
+  else if (isGamma)  AliDebug(2, Form("V0 is a Gamma"));
+  //else AliDebug(2, Form("V0 is something else!!"));
+
+  Double_t invmassLc2Lpi = part->InvMassLctoLambdaPi();
+  Double_t invmassLambda = v0part->MassLambda();
+  Double_t invmassLambdaBar = v0part->MassAntiLambda();
+
+  //Double_t nSigmaITSpr=-999.;
+  Double_t nSigmaTPCpr=-999.;
+  Double_t nSigmaTOFpr=-999.;
+
+  //Double_t nSigmaITSpi=-999.;
+  Double_t nSigmaTPCpi=-999.;
+  Double_t nSigmaTOFpi=-999.;
+
+  //Double_t nSigmaITSka=-999.;
+  Double_t nSigmaTPCka=-999.;
+  Double_t nSigmaTOFka=-999.;
+
+  /*
+    cutsAnal->GetPidHF()->GetnSigmaITS(bachelor,4,nSigmaITSpr);
+    cutsAnal->GetPidHF()->GetnSigmaTPC(bachelor,4,nSigmaTPCpr);
+    cutsAnal->GetPidHF()->GetnSigmaTOF(bachelor,4,nSigmaTOFpr);
+    cutsAnal->GetPidHF()->GetnSigmaITS(bachelor,2,nSigmaITSpi);
+    cutsAnal->GetPidHF()->GetnSigmaTPC(bachelor,2,nSigmaTPCpi);
+    cutsAnal->GetPidHF()->GetnSigmaTOF(bachelor,2,nSigmaTOFpi);
+    cutsAnal->GetPidHF()->GetnSigmaITS(bachelor,3,nSigmaITSka);
+    cutsAnal->GetPidHF()->GetnSigmaTPC(bachelor,3,nSigmaTPCka);
+    cutsAnal->GetPidHF()->GetnSigmaTOF(bachelor,3,nSigmaTOFka);
+  */
+
+  nSigmaTPCpi = fPIDResponse->NumberOfSigmasTPC(bachelor,(AliPID::kPion));
+  nSigmaTPCka = fPIDResponse->NumberOfSigmasTPC(bachelor,(AliPID::kKaon));
+  nSigmaTPCpr = fPIDResponse->NumberOfSigmasTPC(bachelor,(AliPID::kProton));
+  nSigmaTOFpi = fPIDResponse->NumberOfSigmasTOF(bachelor,(AliPID::kPion));
+  nSigmaTOFka = fPIDResponse->NumberOfSigmasTOF(bachelor,(AliPID::kKaon));
+  nSigmaTOFpr = fPIDResponse->NumberOfSigmasTOF(bachelor,(AliPID::kProton));
+
+
+  // Fill candidate variable Tree (track selection, V0 invMass selection)
+  if (!onFlyV0 && isInV0window && isInCascadeWindow && part->CosV0PointingAngle()>0.99 && TMath::Abs(nSigmaTPCpr) <= 3 && v0part->Getd0Prong(0) < 20 && v0part->Getd0Prong(1) < 20) {
+
+    fCandidateVariables[0] = invmassLc;
+    fCandidateVariables[1] = invmassLc2Lpi;
+    fCandidateVariables[2] = invmassK0s;
+    fCandidateVariables[3] = invmassLambda;
+    fCandidateVariables[4] = invmassLambdaBar;
+    fCandidateVariables[5] = part->CosV0PointingAngle();
+    fCandidateVariables[6] = dcaV0;
+    fCandidateVariables[7] = part->Getd0Prong(0);
+    fCandidateVariables[8] = part->Getd0Prong(1);
+    fCandidateVariables[9] = nSigmaTPCpr;
+    fCandidateVariables[10] = nSigmaTPCpi;
+    fCandidateVariables[11] = nSigmaTPCka;
+    fCandidateVariables[12] = nSigmaTOFpr;
+    fCandidateVariables[13] = nSigmaTOFpi;
+    fCandidateVariables[14] = nSigmaTOFka;
+    fCandidateVariables[15] = bachelor->Pt();
+    AliAODTrack *v0pos = (AliAODTrack*)part->Getv0PositiveTrack();
+    fCandidateVariables[16] = v0pos->Pt();
+    AliAODTrack *v0neg = (AliAODTrack*)part->Getv0NegativeTrack(); 
+    fCandidateVariables[17] = v0neg->Pt();
+    fCandidateVariables[18] = v0part->Getd0Prong(0);
+    fCandidateVariables[19] = v0part->Getd0Prong(1);
+    fCandidateVariables[20] = v0part->Pt();
+    fCandidateVariables[21] = v0part->InvMass2Prongs(0,1,11,11);
+    fCandidateVariables[22] = part->Pt();
+    fCandidateVariables[23] = probProton;
+    fCandidateVariables[24] = part->Eta();
+    fCandidateVariables[25] = v0pos->Eta();
+    fCandidateVariables[26] = v0neg->Eta();
+    fCandidateVariables[27] = probPion;
+    fCandidateVariables[28] = probKaon;
+    fCandidateVariables[29] = bachelor->Eta();
+
+    fCandidateVariables[30] = part->P();
+    fCandidateVariables[31] = bachelor->P();
+    fCandidateVariables[32] = v0part->P();
+    fCandidateVariables[33] = v0pos->P();
+    fCandidateVariables[34] = v0neg->P();
+
+    fCandidateVariables[35] = part->Y(4122);
+    fCandidateVariables[36] = bachelor->Y(2212);
+    fCandidateVariables[37] = v0part->Y(310);
+    fCandidateVariables[38] = v0pos->Y(211);
+    fCandidateVariables[39] = v0neg->Y(211);
+
+    fCandidateVariables[40] = v0part->Eta();
+
+    fCandidateVariables[41] = part->DecayLength();
+    fCandidateVariables[42] = part->DecayLengthV0();
+    fCandidateVariables[43] = part->Ct(4122);
+    fCandidateVariables[44] = v0part->Ct(310, v0part->GetSecondaryVtx());
+
+    EBachelor bachCode = CheckBachelor(part, bachelor, mcArray);         
+    EK0S k0SCode = CheckK0S(part, v0part, mcArray);
+
+    fCandidateVariables[45] = bachCode;
+    fCandidateVariables[46] = k0SCode;
+
+    Double_t V0KF[3] = {-999999, -999999, -999999};
+    Double_t errV0KF[3] = {-999999, -999999, -999999};
+    Double_t LcKF[3] = {-999999, -999999, -999999};
+    Double_t errLcKF[3] = {-999999, -999999, -999999};
+    Double_t distances[3] = {-999999, -999999, -999999};
+    Double_t armPolKF[2] = {-999999, -999999};
+
+    if (fCallKFVertexing){
+      Int_t kfResult = CallKFVertexing(part, v0part, bachelor, mcArray, &V0KF[0], &errV0KF[0], &LcKF[0], &errLcKF[0], &distances[0], &armPolKF[0]);
+      AliDebug(2, Form("Result from KF = %d", kfResult));
+    }
+
+    /*
+      for (Int_t i = 0; i< 3; i++){
+      Printf("i = %d, V0KF = %f, errV0KF = %f, LcKF = %f, errLcKF = %f", V0KF[i], errV0KF[i], LcKF[i], errLcKF[i]);
+      }
+    */
+
+    fCandidateVariables[47] = V0KF[0];
+    fCandidateVariables[48] = V0KF[1];
+    fCandidateVariables[49] = V0KF[2];
+
+    fCandidateVariables[50] = errV0KF[0];
+    fCandidateVariables[51] = errV0KF[1];
+    fCandidateVariables[52] = errV0KF[2];
+
+    fCandidateVariables[53] = LcKF[0];
+    fCandidateVariables[54] = LcKF[1];
+    fCandidateVariables[55] = LcKF[2];
+
+    fCandidateVariables[56] = errLcKF[0];
+    fCandidateVariables[57] = errLcKF[1];
+    fCandidateVariables[58] = errLcKF[2];
+
+    fCandidateVariables[59] = distances[0];
+    fCandidateVariables[60] = distances[1];
+    fCandidateVariables[61] = distances[2];
+    fCandidateVariables[62] = armPolKF[0];
+    fCandidateVariables[63] = armPolKF[1];
+    fCandidateVariables[64] = v0part->AlphaV0();
+    fCandidateVariables[65] = v0part->PtArmV0();
+
+    if (fUseMCInfo) {
+      if (isLc){
+       AliDebug(2, "Filling Sgn");
+       fVariablesTreeSgn->Fill();
+       fHistoCodesSgn->Fill(bachCode, k0SCode);          
+      }
+      else {
+       if (fFillOnlySgn == kFALSE){
+         AliDebug(2, "Filling Bkg");
+         fVariablesTreeBkg->Fill();
+         fHistoCodesBkg->Fill(bachCode, k0SCode);        
+       }
+      }
+    }
+    else {
+      fVariablesTreeSgn->Fill();
+    }
+  }
+       
+  return;
+
+}
+
+//________________________________________________________________________
+Int_t AliAnalysisTaskSELc2V0bachelorTMVA::CallKFVertexing(AliAODRecoCascadeHF *cascade, AliAODv0* v0part, AliAODTrack* bach, TClonesArray *mcArray,
+                                                         Double_t* V0KF, Double_t* errV0KF, Double_t* LcKF, Double_t* errLcKF, 
+                                                         Double_t* distances, Double_t* armPolKF) {
+
+  //
+  // method to perform KF vertexing 
+  // elements: [0] = mass, [1] = DecayLength, [2] = lifeTime
+  //
+  Int_t codeKFV0 = -1, codeKFLc = -1;
+
+  AliKFVertex primVtxCopy;
+  Int_t nt = 0, ntcheck = 0;
+  Double_t pos[3] = {0., 0., 0.};
+
+  fVtx1->GetXYZ(pos);
+  Int_t contr = fVtx1->GetNContributors();
+  Double_t covmatrix[6] = {0.};
+  fVtx1->GetCovarianceMatrix(covmatrix);
+  Double_t chi2 = fVtx1->GetChi2();
+  AliESDVertex primaryESDVtxCopy(pos, covmatrix, chi2, contr, "Vertex");
+
+  // topological constraint
+  primVtxCopy = AliKFVertex(primaryESDVtxCopy);
+  nt = primaryESDVtxCopy.GetNContributors();
+  ntcheck = nt;
+  
+  Int_t pdg[2] = {211, -211};
+  Int_t pdgLc[2] = {2212, 310};
+  
+  Int_t pdgDgV0toDaughters[2] = {211, 211};
+
+  Int_t mcLabelV0 = v0part->MatchToMC(310, mcArray, 2, pdgDgV0toDaughters);
+
+  // the KF vertex for the V0 has to be built with the prongs of the V0!
+  Bool_t isMCokV0 = kTRUE, isBkgV0 = kFALSE;
+  AliKFParticle  V0, positiveV0KF, negativeV0KF;;
+  Int_t labelsv0daugh[2] = {-1, -1};
+  for(Int_t ipr= 0; ipr < 2; ipr++){ // 0 is positive, 1 is negative
+    AliAODTrack *aodTrack = (AliAODTrack*)v0part->GetDaughter(ipr);
+    if(!aodTrack) {
+      AliDebug(2, "No V0 daughters available");
+      return -1;
+    }
+    labelsv0daugh[ipr] = TMath::Abs(aodTrack->GetLabel());
+    if (labelsv0daugh[ipr] == -1) isBkgV0 = kTRUE;
+
+    //Printf("v0 daughter %d has label %d", ipr, labelsv0daugh[ipr]);
+    AliKFParticle daughterKF(*aodTrack, pdg[ipr]); // we assume that the PDG is correct
+    if (aodTrack->Charge() > 0) { // assigning positive and negative track to KF V0 for Armenteros-Podolanski plot 
+      positiveV0KF = daughterKF;
+    }
+    else { 
+      negativeV0KF = daughterKF;
+    }  
+    V0.AddDaughter(daughterKF);
+  }
+
+  // Checking the quality of the KF V0 vertex
+  if( V0.GetNDF() < 1 ) {
+    //Printf("Number of degrees of freedom < 1, continuing");
+    return -1;
+  }
+  if( TMath::Sqrt(TMath::Abs(V0.GetChi2()/V0.GetNDF())) > 3. ) {
+    //Printf("Chi2 per DOF too big, continuing");
+    return -1;
+  }
+
+  if(ftopoConstraint && nt > 0){
+    for(Int_t ipr = 0; ipr < 2; ipr++){ // 0 is positive, 1 is negative
+      AliAODTrack *aodTrack = (AliAODTrack*)v0part->GetDaughter(ipr);
+      //* subtruct daughters from primary vertex 
+      if(!aodTrack->GetUsedForPrimVtxFit()) {
+       //Printf("Track %d was not used for primary vertex, continuing", i);
+       continue;
+      }
+      AliKFParticle daughterKF(*aodTrack, pdg[ipr]); // we assume that the PDG is correct
+      primVtxCopy -= daughterKF;
+      ntcheck--;
+    }
+  }
+
+  /*
+  // Check V0 Chi^2 deviation from primary vertex // not needed for V0 for Lc decay!!
+  if( V0.GetDeviationFromVertex( primVtxCopy ) >3. ) {
+  //Printf("Deviation from vertex too big, continuing");
+  return -1;
+  }
+  */
+  
+  if(ftopoConstraint){
+    if(ntcheck>0) { // number of constraints left with which the vertex was built, after removing V0 daughters
+      //* Add V0 to primary vertex to improve the primary vertex resolution
+      primVtxCopy += V0;
+      V0.SetProductionVertex(primVtxCopy); // FIXME: is this correct here? 
+    }
+  }
+
+  //* Check chi^2 
+  if( TMath::Sqrt( TMath::Abs(V0.GetChi2()/V0.GetNDF())) >3. ) {
+    AliDebug(2, "Final Chi2 per DOF too big, continuing");
+    return -1;
+  }
+
+  //* Get V0 invariant mass
+  Double_t massV0 = 999999, sigmaMassV0 = 999999;
+  Int_t retMV0 = V0.GetMass( massV0, sigmaMassV0 ); 
+  if( retMV0 ) {
+    if (massV0 < 0) {
+      codeKFV0 = 1; // Mass not ok
+      if (sigmaMassV0 > 1e19) codeKFV0 = 5;  // Mass and SigmaMass not ok
+    }
+    else if (sigmaMassV0 > 1e19) codeKFV0 = 2; // SigmaMass not ok
+  }
+  fHistoMassKFV0->Fill(massV0, sigmaMassV0);   
+
+  //* Get V0 decayLength
+  Double_t decayLengthV0 = 999999, sigmaDecayLengthV0 = 999999;
+  Int_t retDLV0 = V0.GetDecayLength( decayLengthV0, sigmaDecayLengthV0 );
+  if( retDLV0 ) {
+    if (sigmaDecayLengthV0 > 1e19) {  
+      codeKFV0 = 3; // DecayLength not ok
+      if (massV0 < 0) {
+       codeKFV0 = 6; // DecayLength and Mass not ok      
+       if (sigmaMassV0 > 1e19) codeKFV0 = 11;  // DecayLength and Mass and SigmaMass not ok
+      }
+      else if (sigmaMassV0 > 1e19) codeKFV0 = 8;  // DecayLength and SigmaMass not ok
+    }
+  }
+  fHistoDecayLengthKFV0->Fill(decayLengthV0, sigmaDecayLengthV0);      
+
+  //* Get V0 life time
+  Double_t lifeTimeV0 = 999999, sigmaLifeTimeV0 = 999999;
+  Int_t retTLV0 = V0.GetLifeTime( lifeTimeV0, sigmaLifeTimeV0 );
+  if( retTLV0 ) {
+    if (sigmaLifeTimeV0 > 1e19) {
+      codeKFV0 = 4;  // LifeTime not ok
+      if (sigmaDecayLengthV0 > 1e19) {
+       codeKFV0 = 9;  // LifeTime and DecayLength not ok
+       if (massV0 < 0) {
+         codeKFV0 = 14; // LifeTime and DecayLength and Mass not ok
+         if (sigmaMassV0 > 1e19) codeKFV0 = 15; // LifeTime and DecayLength and Mass and SigmaMass not ok
+       }
+       else if (sigmaMassV0 > 1e19) codeKFV0 = 13; // LifeTime and DecayLength and SigmaMass not ok
+      }
+      else if (massV0 < 0) { // LifeTime and Mass and SigmaMass not ok
+       codeKFV0 = 7; // LifeTime and Mass not ok
+       if (sigmaMassV0 > 1e19) codeKFV0 = 12; // LifeTime and Mass and SigmaMass not ok
+      }
+      else if (sigmaMassV0 > 1e19) codeKFV0 = 10;  // LifeTime and SigmaMass not ok      
+    }
+  }  
+  fHistoLifeTimeKFV0->Fill(lifeTimeV0, sigmaLifeTimeV0);       
+
+  if (codeKFV0 == -1) codeKFV0 = 0;
+  fHistoKFV0->Fill(codeKFV0);
+
+  AliDebug(2, Form("V0: mass = %f, decay length = %f, life time = %f", massV0, decayLengthV0, lifeTimeV0 ));
+  
+  fHistoMassV0All->Fill(massV0); 
+  fHistoDecayLengthV0All->Fill(decayLengthV0);
+  fHistoLifeTimeV0All->Fill(lifeTimeV0);
+
+  Double_t qtAlphaV0[2] = {0., 0.};
+  Double_t vtxV0KF[3] = {V0.GetX(), V0.GetY(), V0.GetZ()}; 
+  positiveV0KF.TransportToPoint(vtxV0KF);
+  negativeV0KF.TransportToPoint(vtxV0KF);
+  V0.GetArmenterosPodolanski(positiveV0KF, negativeV0KF, qtAlphaV0);
+  AliDebug(2, Form("Armenteros-Podolanski variables: alpha = %f, qt = %f", qtAlphaV0[1], qtAlphaV0[0])); 
+  fHistoArmenterosPodolanskiV0KF->Fill(qtAlphaV0[1], qtAlphaV0[0]);
+  fHistoArmenterosPodolanskiV0AOD->Fill(v0part->AlphaV0(), v0part->PtArmV0());
+  armPolKF[0] = qtAlphaV0[1];
+  armPolKF[1] = qtAlphaV0[0];
+
+  // Checking MC info for V0
+
+  AliAODMCParticle *motherV0 = 0x0;
+  AliAODMCParticle *daughv01 = 0x0;
+  AliAODMCParticle *daughv02 = 0x0;
+
+  if (fUseMCInfo) {
+    daughv01 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelsv0daugh[0]));
+    daughv02 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelsv0daugh[1]));
+    if (!daughv01 && labelsv0daugh[0] > 0){
+      AliDebug(2, "Could not access MC info for first daughter of V0, continuing");
+      isMCokV0 = kFALSE;
+    }
+    if (!daughv02 && labelsv0daugh[1] > 0){
+      AliDebug(2, "Could not access MC info for second daughter of V0, continuing");
+      isMCokV0 = kFALSE;
+    }
+    if (isMCokV0){
+      if( daughv01->GetMother() ==  daughv02->GetMother() && daughv01->GetMother()>=0 ){
+       AliDebug(3, Form("The mother has label %d", daughv01->GetMother()));
+       motherV0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughv01->GetMother()));
+       if( motherV0 && TMath::Abs(motherV0->GetPdgCode()) != 21 ){ // These are all V0 that are truly V0, not only K0S, but no gluons 
+         if( motherV0->GetNDaughters() == 2 ){
+           fHistoMassV0True->Fill(massV0);     
+           fHistoDecayLengthV0True->Fill(decayLengthV0);
+           fHistoLifeTimeV0True->Fill(lifeTimeV0);
+           fHistoMassV0TrueFromAOD->Fill(v0part->MassK0Short());
+           if (TMath::Abs(motherV0->GetPdgCode()) == 310){ // These are true V0 that are also K0S
+             fHistoMassV0TrueK0S->Fill(massV0);        
+             fHistoDecayLengthV0TrueK0S->Fill(decayLengthV0);
+             fHistoLifeTimeV0TrueK0S->Fill(lifeTimeV0);
+             fHistoMassV0TrueK0SFromAOD->Fill(v0part->MassK0Short());
+           }
+         }
+         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()));    
+       }
+       else if (!motherV0){
+         AliDebug(3, "could not access MC info for mother, continuing");
+       }
+       else {
+         AliDebug(3, "MC mother is a gluon, continuing");
+       }
+      }
+      else {
+       AliDebug(3, "Background V0!");
+       isBkgV0 = kTRUE; 
+      }
+    }
+  }
+      
+  AliDebug(2, Form("isMCokV0 = %d, isBkgV0 = %d", (Int_t)isMCokV0, (Int_t)isBkgV0));
+
+  // now use what we just try with the bachelor, to build the Lc
+  
+  // topological constraint
+  nt = primVtxCopy.GetNContributors();
+  ntcheck = nt;
+  
+  Bool_t isMCokLc = kTRUE, isBkgLc = kFALSE;
+  AliKFParticle  Lc;
+  Int_t labelsLcdaugh[2] = {-1, -1};
+  labelsLcdaugh[0] = TMath::Abs(bach->GetLabel());
+  labelsLcdaugh[1] = mcLabelV0;
+
+  if (bach->Charge() < 0) pdgLc[0] = -pdgLc[0];
+  AliKFParticle daughterKFLc(*bach, pdgLc[0]);
+  Lc.AddDaughter(daughterKFLc);
+  Lc.AddDaughter(V0);
+  if( Lc.GetNDF() < 1 ) {
+    AliDebug(2, Form("Lc: Number of degrees of freedom < 1 (%d), continuing", Lc.GetNDF()));
+    return -1;
+  }
+  if( TMath::Sqrt(TMath::Abs(Lc.GetChi2()/Lc.GetNDF())) >3. ) {
+    AliDebug(2, Form("Lc: Chi2 per DOF too big, continuing (%f)", TMath::Sqrt(TMath::Abs(Lc.GetChi2()/Lc.GetNDF()))));
+    return -1;
+  }
+  
+  if(ftopoConstraint && nt > 0){
+    //* subtruct daughters from primary vertex 
+    if(!bach->GetUsedForPrimVtxFit()) {
+      AliDebug(3, "Lc: Bachelor was not used for primary vertex, not subtracting it from primary vertex");
+    }
+    else{
+      primVtxCopy -= daughterKFLc;
+      ntcheck--;
+    }
+    /* the V0 was added above, so it is ok to remove it without checking
+       if(!V0->GetUsedForPrimVtxFit()) {
+       Printf("Lc: V0 was not used for primary vertex, continuing");
+       continue;
+       }
+    */
+    primVtxCopy -= V0;
+    ntcheck--;
+  }
+  
+  //* Check V0 Chi^2 deviation from primary vertex 
+  if( Lc.GetDeviationFromVertex( primVtxCopy ) >3. ) {
+    AliDebug(2, Form("Lc: Deviation from vertex too big, continuing (%f)", Lc.GetDeviationFromVertex( primVtxCopy )));
+    return -1;
+  }
+  
+  if(ftopoConstraint){
+    if(ntcheck>0) {
+      //* Add Lc to primary vertex to improve the primary vertex resolution
+      primVtxCopy += Lc;
+      Lc.SetProductionVertex(primVtxCopy);
+    }
+  }
+  
+  //* Check chi^2 
+  if( TMath::Sqrt( TMath::Abs(Lc.GetChi2()/Lc.GetNDF())) >3. ) {
+    AliDebug(2, Form("Lc: Final Chi2 per DOF too big, continuing (%f)", TMath::Sqrt( TMath::Abs(Lc.GetChi2()/Lc.GetNDF()))));
+    return -1;
+  }
+
+  //* Get Lc invariant mass
+  Double_t massLc = 999999, sigmaMassLc= 999999;
+  Int_t retMLc = Lc.GetMass( massLc, sigmaMassLc ); 
+  if( retMLc ) {
+    AliDebug(3, Form("----> Could not get mass (%e), and sigma(%e) for Lc, continuing", massLc, sigmaMassLc));
+    if (massLc < 0) {
+      codeKFLc = 1; // Mass not ok
+      if (sigmaMassLc > 1e19) codeKFLc = 5;  // Mass and SigmaMass not ok
+    }
+    else if (sigmaMassLc > 1e19) codeKFLc = 2; // SigmaMass not ok
+  }
+  fHistoMassKFLc->Fill(massLc, sigmaMassLc);   
+
+  //* Get Lc Decay length
+  Double_t decayLengthLc = 999999, sigmaDecayLengthLc = 999999;
+  Int_t retDLLc = Lc.GetDecayLength( decayLengthLc, sigmaDecayLengthLc );
+  if( retDLLc ) {
+    AliDebug(3, "----> Lc: Could not get decay length, and sigma");
+    if (sigmaDecayLengthLc > 1e19) {  
+      codeKFLc = 3; // DecayLength not ok
+      if (massLc < 0) {
+       codeKFLc = 6; // DecayLength and Mass not ok      
+       if (sigmaMassLc > 1e19) codeKFLc = 11;  // DecayLength and Mass and SigmaMass not ok
+      }
+      else if (sigmaMassLc > 1e19) codeKFLc = 8;  // DecayLength and SigmaMass not ok
+    }
+  }
+  AliDebug(3, Form("retDLLc = %d, with decayLength = %f and error = %e", retDLLc, decayLengthLc, sigmaDecayLengthLc));
+
+  fHistoDecayLengthKFLc->Fill(decayLengthLc, sigmaDecayLengthLc);      
+
+  //* Get Lc life time
+  Double_t lifeTimeLc = 999999, sigmaLifeTimeLc = 999999;
+  Int_t retTLLc = Lc.GetLifeTime( lifeTimeLc, sigmaLifeTimeLc );
+  if( retTLLc ) {
+    AliDebug(3, "----> Lc: Could not get lifeTime, and sigma");
+    if (sigmaLifeTimeLc > 1e19) {
+      codeKFLc = 4;  // LifeTime not ok
+      if (sigmaDecayLengthLc > 1e19) {
+       codeKFLc = 9;  // LifeTime and DecayLength not ok
+       if (massLc < 0) {
+         codeKFLc = 14; // LifeTime and DecayLength and Mass not ok
+         if (sigmaMassLc > 1e19) codeKFLc = 15; // LifeTime and DecayLength and Mass and SigmaMass not ok
+       }
+       else if (sigmaMassLc > 1e19) codeKFLc = 13; // LifeTime and DecayLength and SigmaMass not ok
+      }
+      else if (massLc < 0) { // LifeTime and Mass and SigmaMass not ok
+       codeKFLc = 7; // LifeTime and Mass not ok
+       if (sigmaMassLc > 1e19) codeKFLc = 12; // LifeTime and Mass and SigmaMass not ok
+      }
+      else if (sigmaMassLc > 1e19) codeKFLc = 10;  // LifeTime and SigmaMass not ok      
+    }
+  }
+
+  fHistoLifeTimeKFLc->Fill(lifeTimeLc, sigmaLifeTimeLc);       
+
+  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));
+  
+  if (codeKFLc == -1) codeKFLc = 0;
+  fHistoKFLc->Fill(codeKFLc);
+  
+  fHistoKF->Fill(codeKFV0, codeKFLc);
+
+  // here we fill the histgrams for all the reconstructed KF vertices for the cascade
+  fHistoMassLcAll->Fill(massLc);            
+  fHistoDecayLengthLcAll->Fill(decayLengthLc);            
+  fHistoLifeTimeLcAll->Fill(lifeTimeLc);            
+
+  fHistoMassV0fromLcAll->Fill(massV0); 
+  fHistoDecayLengthV0fromLcAll->Fill(decayLengthV0);
+  fHistoLifeTimeV0fromLcAll->Fill(lifeTimeV0);
+  
+  Double_t xV0 = V0.GetX();
+  Double_t yV0 = V0.GetY();
+  Double_t zV0 = V0.GetZ();
+  
+  Double_t xLc = Lc.GetX();
+  Double_t yLc = Lc.GetY();
+  Double_t zLc = Lc.GetZ();
+  
+  Double_t xPrimVtx = primVtxCopy.GetX();
+  Double_t yPrimVtx = primVtxCopy.GetY();
+  Double_t zPrimVtx = primVtxCopy.GetZ();
+  
+  Double_t distanceLcToPrimVtx = TMath::Sqrt((xPrimVtx - xLc) * (xPrimVtx - xLc) +
+                                            (yPrimVtx - yLc) * (yPrimVtx - yLc) +
+                                            (zPrimVtx - zLc) * (zPrimVtx - zLc));
+  
+  Double_t distanceV0ToPrimVtx = TMath::Sqrt((xPrimVtx - xV0) * (xPrimVtx - xV0) +
+                                            (yPrimVtx - yV0) * (yPrimVtx - yV0) +
+                                            (zPrimVtx - zV0) * (zPrimVtx - zV0));
+  
+  Double_t distanceV0ToLc = TMath::Sqrt((xLc - xV0)*(xLc - xV0) +
+                                       (yLc - yV0)*(yLc - yV0) +
+                                       (zLc - zV0)*(zLc - zV0));
+  
+  //Printf("distanceLcToPrimVtx = %e, distanceV0ToPrimVtx= %f, distanceV0ToLc = %f", distanceLcToPrimVtx, distanceV0ToPrimVtx, distanceV0ToLc);
+  
+  fHistoDistanceLcToPrimVtx->Fill(distanceLcToPrimVtx);
+  fHistoDistanceV0ToPrimVtx->Fill(distanceV0ToPrimVtx);
+  fHistoDistanceV0ToLc->Fill(distanceV0ToLc);
+  
+  distances[0] = distanceLcToPrimVtx;
+  distances[1] = distanceV0ToPrimVtx;
+  distances[2] = distanceV0ToLc;
+
+  if (fUseMCInfo) {
+
+    AliAODMCParticle *daughv01Lc = 0x0;
+    AliAODMCParticle *K0S = 0x0;
+    AliAODMCParticle *daughv02Lc = 0x0;
+    
+    if (labelsLcdaugh[0] >= 0) {
+      //      Printf("Getting Bachelor from label %d", labelsLcdaugh[1]);
+      daughv01Lc = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelsLcdaugh[0]));
+      if (!daughv01Lc){
+       AliDebug(3, "Could not access MC info for first daughter of Lc");
+       isMCokLc = kFALSE;
+      }
+      else {
+       AliDebug(2, Form("The bachelor has label = %d", daughv01Lc->GetLabel()));
+       if (TMath::Abs(daughv01Lc->GetPdgCode()) != 2212) isBkgLc = kTRUE;      
+      }
+    }  
+    else { // The bachelor is a fake
+      isBkgLc = kTRUE;
+    }
+
+    if (labelsLcdaugh[1] >= 0) {
+      //Printf("Getting K0S");
+      K0S = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelsLcdaugh[1]));    
+      if (!K0S) {
+       AliDebug(3, "Could not access MC info for second daughter of Lc");
+       isMCokLc = kFALSE;
+      }
+      else{
+       if (TMath::Abs(K0S->GetPdgCode()) != 310) isBkgLc = kTRUE; 
+      }      
+    }
+    else{
+      AliDebug(2, "The K0S is not true --> it does not have a label, continuing...");
+      isBkgLc = kTRUE;
+    }
+  
+    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
+      if (isMCokLc) { // We can then access its MC info, and it might then be that also the Lc is a true Lc
+       Int_t iK0 = K0S->GetMother();
+       if (iK0 < 0) {
+         Printf("The K0S has no mother... IMPOSSIBLE"); // the K0S MUST have a mother!
+       }
+       else { // The K0S has a mother
+         daughv02Lc = dynamic_cast<AliAODMCParticle*>(mcArray->At(iK0));
+         if (!daughv02Lc){
+           AliDebug(3, "Could not access MC info for second daughter of Lc");
+         }
+         else { // we can access safely the K0S mother in the MC
+           if( daughv01Lc->GetMother() ==  daughv02Lc->GetMother() && daughv01Lc->GetMother()>=0 ){  // This is a true cascade! bachelor and V0 come from the same mother
+             //Printf("Lc: The mother has label %d", daughv01Lc->GetMother());
+             AliAODMCParticle *motherLc = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughv01Lc->GetMother()));
+             Int_t pdgMum = 0, pdgBach = 0, pdgV0 = 0;
+             if (motherLc) pdgMum = motherLc->GetPdgCode();
+             if (daughv01Lc) pdgBach = daughv01Lc->GetPdgCode();
+             if (daughv02Lc) pdgV0 = daughv02Lc->GetPdgCode();
+             AliDebug(2, Form("pdgMum = %d, pdgBach = %d, pdgV0 = %d", pdgMum, pdgBach, pdgV0));
+             
+             if( motherLc && TMath::Abs(motherLc->GetPdgCode()) != 21 ){ // These are true cascades, we don't know yet if they are Lc    
+               fHistoMassLcTrue->Fill(massLc); 
+               fHistoDecayLengthLcTrue->Fill(decayLengthLc);            
+               fHistoLifeTimeLcTrue->Fill(lifeTimeLc);            
+               fHistoMassLcTrueFromAOD->Fill(cascade->InvMassLctoK0sP());
+               
+               fHistoMassV0fromLcTrue->Fill(massV0);   
+               fHistoDecayLengthV0fromLcTrue->Fill(decayLengthV0);            
+               fHistoLifeTimeV0fromLcTrue->Fill(lifeTimeV0);            
+               
+               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...
+                 AliDebug(2, Form("IT IS SIGNAL!!! with label = %d", motherLc->GetLabel()));
+
+                 fHistoArmenterosPodolanskiV0KFSgn->Fill(qtAlphaV0[1], qtAlphaV0[0]);
+                 fHistoArmenterosPodolanskiV0AODSgn->Fill(v0part->AlphaV0(), v0part->PtArmV0());
+
+                 fHistoDistanceLcToPrimVtxSgn->Fill(distanceLcToPrimVtx);
+                 fHistoDistanceV0ToPrimVtxSgn->Fill(distanceV0ToPrimVtx);
+                 fHistoDistanceV0ToLcSgn->Fill(distanceV0ToLc);
+
+                 fHistoMassLcSgn->Fill(massLc);        
+                 fHistoMassLcSgnFromAOD->Fill(cascade->InvMassLctoK0sP());
+
+                 fHistoDecayLengthLcSgn->Fill(decayLengthLc);            
+                 fHistoLifeTimeLcSgn->Fill(lifeTimeLc);            
+
+                 fHistoMassV0fromLcSgn->Fill(massV0);  
+                 fHistoDecayLengthV0fromLcSgn->Fill(decayLengthV0);            
+                 fHistoLifeTimeV0fromLcSgn->Fill(lifeTimeV0);    
+               }        
+               else {
+                 isBkgLc = kTRUE;  // it is not a Lc, since the pdg != 4122
+               }
+               
+               // if we got here, we can compare with MC information; this part is done also for cases where the cascade is not a Lc!
+               // First, we compare the vtx of the Lc
+               Double_t xLcMC = motherLc->Xv();
+               Double_t yLcMC = motherLc->Yv();
+               Double_t zLcMC = motherLc->Zv();
+               //Printf("Got MC vtx for Lc");
+               Double_t xProtonMC = daughv01Lc->Xv();
+               Double_t yProtonMC = daughv01Lc->Yv();
+               Double_t zProtonMC = daughv01Lc->Zv();
+               //Printf("Got MC vtx for Proton");
+               
+               Double_t vtxLcResidualToPrimVtx = TMath::Sqrt((xLcMC - xProtonMC) * (xLcMC - xProtonMC) + 
+                                                             (yLcMC - yProtonMC) * (yLcMC - yProtonMC) + 
+                                                             (zLcMC - zProtonMC) * (zLcMC - zProtonMC)) - distanceLcToPrimVtx;
+               
+               // Then, we compare the vtx of the K0s
+               Double_t xV0MC = motherV0->Xv();     //Production vertex of K0S --> Where Lc decays
+               Double_t yV0MC = motherV0->Yv();
+               Double_t zV0MC = motherV0->Zv();
+               //Printf("Got MC vtx for V0");
+               Double_t xPionMC = daughv01->Xv(); //Production vertex of Pion --> Where K0S decays
+               Double_t yPionMC = daughv01->Yv();
+               Double_t zPionMC = daughv01->Zv();
+               //Printf("Got MC vtx for Pion");
+               
+               Double_t vtxV0ResidualToLc = TMath::Sqrt((xV0MC - xPionMC) * (xV0MC - xPionMC) + 
+                                                        (yV0MC - yPionMC) * (yV0MC - yPionMC) + 
+                                                        (zV0MC - zPionMC) * (zV0MC - zPionMC)) - distanceV0ToLc;
+               
+               // Then, the K0S vertex wrt primary vertex
+               
+               Double_t vtxV0ResidualToPrimVtx = TMath::Sqrt((xPionMC - xLcMC) * (xPionMC - xLcMC) + 
+                                                             (yPionMC - yLcMC) * (yPionMC - yLcMC) + 
+                                                             (zPionMC - zLcMC) * (zPionMC - zLcMC)) - distanceV0ToPrimVtx;
+               
+               fHistoVtxLcResidualToPrimVtx->Fill(vtxLcResidualToPrimVtx);
+               fHistoVtxV0ResidualToLc->Fill(vtxV0ResidualToLc);
+               fHistoVtxV0ResidualToPrimVtx->Fill(vtxV0ResidualToPrimVtx);
+               
+             } //closing: if( motherLc && TMath::Abs(motherLc->GetPdgCode()) != 21 )
+             else if (!motherLc){
+               AliDebug(2, "We could not access MC info for Lc mother, so we did nothing");
+             }
+             else {
+               AliDebug(2, "MC Lc mother is a gluon, so we did nothing");
+             }
+           } //closing: if( daughv01Lc->GetMother() ==  daughv02Lc->GetMother() && daughv01Lc->GetMother()>=0 )
+           else {
+             isBkgLc = kTRUE; // it cannot be a Lc, since the daughters do not have the same mother
+           }
+         } // closing: else { // we can access safely the K0S mother in the MC
+       } // closing: else { // The K0S has a mother
+      } // closing isMCLcok
+    }  // closing !isBkgLc
+  } // closing fUseMCInfo
+
+  //Printf("retMV0 = %d, retMLc = %d", retMV0, retMLc); 
+  if ( retMV0 == 0 && retMLc == 0){
+    V0KF[0] = massV0;
+    errV0KF[0] = sigmaMassV0;
+    V0KF[1] = decayLengthV0;
+    errV0KF[1] = sigmaDecayLengthV0;
+    V0KF[2] = lifeTimeV0;
+    errV0KF[2] = sigmaLifeTimeV0;
+    LcKF[0] = massLc;
+    errLcKF[0] = sigmaMassLc;
+    LcKF[1] = decayLengthLc;
+    errLcKF[1] = sigmaDecayLengthLc;
+    LcKF[2] = lifeTimeLc;
+    errLcKF[2] = sigmaLifeTimeLc;
+  }
+
+  return 1;
+
+}
+//________________________________________________________________________
+AliAnalysisTaskSELc2V0bachelorTMVA::EBachelor AliAnalysisTaskSELc2V0bachelorTMVA::CheckBachelor( AliAODRecoCascadeHF *part,
+                                                                                                AliAODTrack* bachelor,
+                                                                                                TClonesArray *mcArray ){
+  
+  //Printf("In CheckBachelor");
+
+  // function to check if the bachelor is a p, if it is a p but not from Lc
+  // to be filled for background candidates
+
+  Int_t label = bachelor->GetLabel();
+  if (label == -1) {
+    return kBachFake;
+  }
+
+  AliAODMCParticle *mcpart = dynamic_cast<AliAODMCParticle*>(mcArray->At(TMath::Abs(label)));
+  if (!mcpart) return kBachInvalid;
+  Int_t pdg = mcpart->PdgCode();
+  if (TMath::Abs(pdg) != 2212) {
+    AliDebug(2, Form("Bachelor is not a p, but a particle with pdg code =  %d", pdg));
+    return kBachNoProton;
+  }
+  else { // it is a proton
+    //Int_t labelLc = part->GetLabel();
+    Int_t labelLc = FindLcLabel(part, mcArray);
+    //Printf(">>>>>>>>>>>>> label for Lc = %d", labelLc);
+    Int_t bachelorMotherLabel = mcpart->GetMother();
+    //Printf(">>>>>>>>>>>>> label for bachelorMother = %d", bachelorMotherLabel);
+    if (bachelorMotherLabel == -1) {
+      return kBachPrimary;
+    }
+    AliAODMCParticle *bachelorMother = dynamic_cast<AliAODMCParticle*>(mcArray->At(bachelorMotherLabel));
+    if (!bachelorMother) return kBachInvalid;
+    Int_t pdgMother = bachelorMother->PdgCode();
+    if (TMath::Abs(pdgMother) != 4122) {
+      AliDebug(2, Form("The proton does not come from a Lc, but from a particle with pdgcode = %d", pdgMother));
+      return kBachNoLambdaMother;
+    }
+    else { // it comes from Lc
+      if (labelLc != bachelorMotherLabel){
+       //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));
+       AliDebug(2, Form("The proton comes from a Lc, but it is not the candidate we are analyzing"));
+       return kBachDifferentLambdaMother;
+      }
+      else { // it comes from the correct Lc
+       AliDebug(2, Form("The proton comes from a Lc, and it is EXACTLY the candidate we are analyzing"));
+       return kBachCorrectLambdaMother;
+      }      
+    }
+  }
+
+  return kBachInvalid;
+
+}
+
+//________________________________________________________________________
+AliAnalysisTaskSELc2V0bachelorTMVA::EK0S AliAnalysisTaskSELc2V0bachelorTMVA::CheckK0S( AliAODRecoCascadeHF *part,
+                                                                                      AliAODv0* v0part,
+                                                                                      //AliAODTrack* v0part,
+                                                                                      TClonesArray *mcArray ){
+  
+  // function to check if the K0Spart is a p, if it is a p but not from Lc
+  // to be filled for background candidates
+
+  //Printf(" CheckK0S");
+
+  //Int_t labelMatchToMC = v0part->MatchToMC(310, mcArray);
+  //Int_t label = v0part->GetLabel();
+  Int_t labelFind = FindV0Label(v0part, mcArray);
+  //Printf("\n >>>>>>>>>>>>> label for V0 = %d, from MatchToMC = %d, from Find = %d", label, labelMatchToMC, labelFind);
+  if (labelFind == -1) {
+    return kK0SFake;
+  }
+
+  AliAODMCParticle *mcpart = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelFind));
+  if (!mcpart) return kK0SInvalid;
+  Int_t pdg = mcpart->PdgCode();
+  if (TMath::Abs(pdg) != 310) {
+    AliDebug(2, Form("K0Spart is not a K0S, but a particle with pdg code =  %d", pdg));
+    //AliInfo(Form("K0Spart is not a K0S, but a particle with pdg code =  %d", pdg));
+    return kK0SNoK0S;
+  }
+  else { // it is a K0S
+    //Int_t labelLc = part->GetLabel();
+    Int_t labelLc = FindLcLabel(part, mcArray);
+    Int_t K0SpartMotherLabel = mcpart->GetMother();
+    if (K0SpartMotherLabel == -1) {
+      return kK0SWithoutMother;
+    }
+    AliAODMCParticle *K0SpartMother = dynamic_cast<AliAODMCParticle*>(mcArray->At(K0SpartMotherLabel));  // should be a K0 in case of signal
+    if (!K0SpartMother) return kK0SInvalid;
+    Int_t pdgMotherK0S = K0SpartMother->PdgCode();
+    if (TMath::Abs(pdgMotherK0S) != 311) {
+      AliDebug(2, Form("The K0S does not come from a K0, but from a particle with pdgcode = %d", pdgMotherK0S));
+      //AliInfo(Form("The K0S does not come from a K0, but from a particle with pdgcode = %d", pdgMotherK0S));
+      return kK0SNotFromK0;
+    }
+    else { // the K0S comes from a K0
+      Int_t K0MotherLabel = K0SpartMother->GetMother(); // mother of K0 --> Lc in case of signal
+      if (K0MotherLabel == -1) {
+       return kK0Primary;
+      }
+      AliAODMCParticle *K0Mother = dynamic_cast<AliAODMCParticle*>(mcArray->At(K0MotherLabel)); 
+      if (!K0Mother) return kK0SInvalid;
+      Int_t pdgK0Mother = K0Mother->PdgCode();
+      if (TMath::Abs(pdgK0Mother) != 4122) { // the K0 does not come from a Lc
+       AliDebug(2, Form("The K0 does not come from a Lc, but from a particle with pdgcode = %d", pdgK0Mother));
+       //AliInfo(Form("The K0 does not come from a Lc, but from a particle with pdgcode = %d", pdgK0Mother));
+       return kK0NoLambdaMother;
+      }
+      else { // the K0 comes from Lc
+       if (labelLc != K0MotherLabel){ // The K0 comes from a different Lc
+         AliDebug(2, Form("The K0S comes from a Lc, but it is not the candidate we are analyzing"));
+         //AliInfo(Form("The K0S comes from a Lc, but it is not the candidate we are analyzing"));
+         return kK0DifferentLambdaMother;
+       }
+       else { // it comes from the correct Lc
+         AliDebug(2, Form("The K0S comes from a Lc, and it is EXACTLY the candidate we are analyzing"));
+         //AliInfo(Form("The K0S comes from a Lc, and it is EXACTLY the candidate we are analyzing"));
+         return kK0CorrectLambdaMother;
+       }      
+      }
+    }
+  }
+
+  return kK0SInvalid;
+
+}
+  
+//----------------------------------------------------------------------------
+Int_t AliAnalysisTaskSELc2V0bachelorTMVA::FindV0Label(AliAODRecoDecay* v0part, TClonesArray *mcArray) const
+{
+
+  //Printf(" FindV0Label");
+
+  // finding the label of teh V0; inspired from AliAODRecoDecay::MatchToMC
+
+  Int_t labMother[2]={-1, -1};
+  AliAODMCParticle *part=0;
+  AliAODMCParticle *mother=0;
+  Int_t dgLabels = -1;
+
+  Int_t ndg = v0part->GetNDaughters();
+  if (ndg != 2) {
+    AliFatal(Form("IMPOSSIBLE!! There are %d daughters, but they should be 2!", ndg));
+  }
+
+  for(Int_t i = 0; i < 2; i++) {
+    AliAODTrack *trk = (AliAODTrack*)v0part->GetDaughter(i);
+    dgLabels = trk->GetLabel();
+    if (dgLabels == -1) {
+      //printf("daughter with negative label %d\n", dgLabels);
+      return -1;
+    }
+    part = (AliAODMCParticle*)mcArray->At(TMath::Abs(dgLabels));
+    if (!part) { 
+      //printf("no MC particle\n");
+      return -1;
+    }
+    labMother[i] = part->GetMother();
+    if (labMother[i] != -1){
+      mother = (AliAODMCParticle*)mcArray->At(labMother[i]);
+      if(!mother) {
+       //printf("no MC mother particle\n");
+       return -1;
+      }
+    }
+    else {
+      return -1;
+    }
+  }
+
+  if (labMother[0] == labMother[1]) return labMother[0];
+
+  else return -1;
+
+}
+
+//----------------------------------------------------------------------------
+Int_t AliAnalysisTaskSELc2V0bachelorTMVA::FindLcLabel(AliAODRecoCascadeHF* cascade, TClonesArray *mcArray) const
+{
+
+  // finding the label of teh V0; inspired from AliAODRecoDecay::MatchToMC
+
+  //Printf(" FindLcLabel");
+
+  AliAODMCParticle *part=0;
+  AliAODMCParticle *mother=0;
+  AliAODMCParticle *grandmother=0;
+  Int_t dgLabels = -1;
+
+  Int_t ndg = cascade->GetNDaughters();
+  if (ndg != 2) {
+    AliFatal(Form("IMPOSSIBLE!! There are %d daughters, but they should be 2!", ndg));
+  }
+
+  // daughter 0 --> bachelor, daughter 1 --> V0
+  AliAODTrack *trk = (AliAODTrack*)cascade->GetDaughter(0); // bachelor
+  dgLabels = trk->GetLabel();
+  if (dgLabels == -1 ) {
+    //printf("daughter with negative label %d\n", dgLabels);
+    return -1;
+  }
+  part = (AliAODMCParticle*)mcArray->At(TMath::Abs(dgLabels));
+  if (!part) { 
+    //printf("no MC particle\n");
+    return -1;
+  }
+  Int_t labMotherBach = part->GetMother();
+  if (labMotherBach == -1){
+    return -1;
+  }
+  mother = (AliAODMCParticle*)mcArray->At(labMotherBach);
+  if(!mother) {
+    //printf("no MC mother particle\n");
+    return -1;
+  }
+  
+  AliAODv0 *v0 = (AliAODv0*)cascade->GetDaughter(1); // V0
+  dgLabels = FindV0Label(v0, mcArray);
+  if (dgLabels == -1 ) {
+    //printf("daughter with negative label (v0 was a fake) %d\n", dgLabels);
+    return -1;
+  }
+  part = (AliAODMCParticle*)mcArray->At(TMath::Abs(dgLabels));
+  if (!part) { 
+    //printf("no MC particle\n");
+    return -1;
+  }
+  Int_t labMotherv0 = part->GetMother();
+  if (labMotherv0 == -1){
+    return -1;
+  }
+  mother = (AliAODMCParticle*)mcArray->At(labMotherv0);
+  if(!mother) {
+    //printf("no MC mother particle\n");
+    return -1;
+  }
+  Int_t labGrandMotherv0 = mother->GetMother();
+  if (labGrandMotherv0 == -1){
+    return -1;
+  }
+  grandmother = (AliAODMCParticle*)mcArray->At(labGrandMotherv0);
+  if(!grandmother) {
+    //printf("no MC mother particle\n");
+    return -1;
+  }
+
+  //Printf("labMotherBach = %d, labMotherv0 = %d, labGrandMotherv0 = %d", labMotherBach, labMotherv0, labGrandMotherv0);
+  if (labGrandMotherv0 == labMotherBach) return labMotherBach;
+
+  else return -1;
+
+}
+  
+//----------------------------------------------------------------------------
+Bool_t AliAnalysisTaskSELc2V0bachelorTMVA::CheckInjection(AliAODRecoCascadeHF *cand, AliAODMCHeader *header,TClonesArray *arrayMC){
+
+  // temporary method until the same is in AliVertexingHFUtils
+
+  AliAODTrack* bach = cand->GetBachelor();
+  if(fUtils->IsTrackInjected(bach, header, arrayMC)) {
+    AliDebug(2, "Bachelor is injected, the whole candidate is then injected");
+    return kTRUE;
+  }
+  AliAODv0* v0 = cand->Getv0();
+  Int_t nprongs = v0->GetNProngs();
+  for(Int_t i = 0; i < nprongs; i++){
+    AliAODTrack *daugh = (AliAODTrack*)v0->GetDaughter(i);
+    if(fUtils->IsTrackInjected(daugh,header,arrayMC)) {
+      AliDebug(2, Form("V0 daughter number %d is injected, the whole candidate is then injected", i));
+      return kTRUE;
+    }
+  }
+  return kFALSE;
+}
+
+
+
+
+
+
diff --git a/PWGHF/vertexingHF/AliAnalysisTaskSELc2V0bachelorTMVA.h b/PWGHF/vertexingHF/AliAnalysisTaskSELc2V0bachelorTMVA.h
new file mode 100644 (file)
index 0000000..9415a79
--- /dev/null
@@ -0,0 +1,244 @@
+#ifndef ALIANALYSISTASKSELC2V0BACHELORTMVA_H
+#define ALIANALYSISTASKSELC2V0BACHELORTMVA_H
+/**************************************************************************
+ * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/* $Id: AliAnalysisTaskSELc2V0bachelorTMVA.h 61835 2013-04-05 23:07:23Z fprino $ */ 
+
+#include "TROOT.h"
+#include "TSystem.h"
+
+#include "AliAnalysisTaskSE.h"
+#include "AliAODEvent.h"
+#include "AliAODMCHeader.h"
+#include "AliPID.h"
+#include "AliAODTrack.h"
+#include "AliAODVertex.h"
+#include "AliAODRecoDecay.h"
+#include "AliPIDResponse.h"
+#include "AliPIDCombined.h"
+#include "AliTPCPIDResponse.h"
+#include "AliRDHFCutsLctoV0.h"
+#include "AliNormalizationCounter.h"
+#include "AliVertexingHFUtils.h"
+#include "AliAODRecoCascadeHF.h"
+
+class TH1F;
+class TH1D;
+
+class AliAnalysisTaskSELc2V0bachelorTMVA : public AliAnalysisTaskSE 
+{
+  
+ public:
+
+  enum EBachelor {
+    kBachInvalid = -1,
+    kBachFake = 0,
+    kBachNoProton = 1,
+    kBachPrimary = 2,
+    kBachNoLambdaMother = 3,
+    kBachDifferentLambdaMother = 4,
+    kBachCorrectLambdaMother = 5 };
+
+  enum EK0S {
+    kK0SInvalid = -1,
+    kK0SFake = 0,
+    kK0SNoK0S = 1,
+    kK0SWithoutMother = 2,
+    kK0SNotFromK0 = 3,
+    kK0Primary = 4, 
+    kK0NoLambdaMother = 5,
+    kK0DifferentLambdaMother = 6,
+    kK0CorrectLambdaMother = 7 };    
+  
+  AliAnalysisTaskSELc2V0bachelorTMVA();
+  AliAnalysisTaskSELc2V0bachelorTMVA(const Char_t* name, AliRDHFCutsLctoV0* cutsA,
+                                Bool_t useOnTheFly=kFALSE);
+  virtual ~AliAnalysisTaskSELc2V0bachelorTMVA();
+
+  // Implementation of interface methods  
+  virtual void UserCreateOutputObjects();
+  virtual void Init();
+  virtual void LocalInit() {Init();}
+  virtual void UserExec(Option_t *option);
+  virtual void Terminate(Option_t *option);
+  // histos
+  void FillLc2pK0Sspectrum(AliAODRecoCascadeHF *part, Int_t isLc,
+                          Int_t &nSelectedAnal, AliRDHFCutsLctoV0 *cutsAnal,
+                          TClonesArray *mcArray);
+
+  void MakeAnalysisForLc2prK0S(TClonesArray *arrayLctopK0s,
+                              TClonesArray *mcArray,
+                              Int_t &nSelectedAnal, AliRDHFCutsLctoV0 *cutsAnal, 
+                              TClonesArray *array3Prong, AliAODMCHeader *aodheader);
+  // set MC usage
+  void SetMC(Bool_t theMCon) {fUseMCInfo = theMCon;}
+  Bool_t GetMC() const {return fUseMCInfo;}
+
+  void SetK0sAnalysis(Bool_t a) {fIsK0sAnalysis=a;}
+  Bool_t GetK0sAnalysis() const {return fIsK0sAnalysis;}
+
+  void SetUseOnTheFlyV0(Bool_t a) { fUseOnTheFlyV0=a; }
+  Bool_t GetUseOnTheFlyV0() { return fUseOnTheFlyV0; }
+
+  void SetIspA(Bool_t a) { fIspA=a; }
+  Bool_t GetIspA() { return fIspA; }
+
+  void SetFillOnlySgn(Bool_t a) { fFillOnlySgn=a; }
+  Bool_t GetFillOnlySgn() { return fFillOnlySgn; }
+
+  void SetTopoConstraint(Bool_t a) { ftopoConstraint=a; }
+  Bool_t GetTopoConstraint() { return ftopoConstraint; }
+
+  void SetCallKFVertexing(Bool_t a) { fCallKFVertexing=a; }
+  Bool_t GetCallKFVertexing() { return fCallKFVertexing; }
+
+  void SetKeepingKeepingOnlyHIJINGBkg(Bool_t a) { fKeepingOnlyHIJINGBkg = a;}
+  Bool_t GetKeepingOnlyHIJINGBkg() {return fKeepingOnlyHIJINGBkg;}
+
+ private:
+  
+  EBachelor CheckBachelor(AliAODRecoCascadeHF *part, AliAODTrack* bachelor, TClonesArray *mcArray);
+  EK0S CheckK0S(AliAODRecoCascadeHF *part, AliAODv0* v0part, TClonesArray *mcArray);
+  //EK0S CheckK0S(AliAODRecoCascadeHF *part, AliAODTrack* v0part, TClonesArray *mcArray );
+  Int_t FindV0Label(AliAODRecoDecay* v0part, TClonesArray *mcArray) const;
+  Int_t FindLcLabel(AliAODRecoCascadeHF* cascade, TClonesArray *mcArray) const;
+  Int_t CallKFVertexing(AliAODRecoCascadeHF *cascade, AliAODv0* v0part, AliAODTrack* bach, TClonesArray *mcArray, 
+                       Double_t* V0KF, Double_t* errV0KF, Double_t* LcKF, Double_t* errLcKF,
+                       Double_t* distances, Double_t* armPolKF);
+  Bool_t CheckInjection(AliAODRecoCascadeHF *cand, AliAODMCHeader *header,TClonesArray *arrayMC);
+
+  AliAnalysisTaskSELc2V0bachelorTMVA(const AliAnalysisTaskSELc2V0bachelorTMVA &source);
+  AliAnalysisTaskSELc2V0bachelorTMVA& operator=(const AliAnalysisTaskSELc2V0bachelorTMVA& source); 
+  
+  Bool_t fUseMCInfo;          // Use MC info
+  TList *fOutput;             // User output1: list of trees
+
+  // define the histograms
+  TH1F *fCEvents;                    // Histogram to check selected events
+  AliPIDResponse *fPIDResponse;      //! PID response object
+  AliPIDCombined *fPIDCombined;      //! combined PID response object
+  Bool_t fIsK0sAnalysis;             // switch between Lpi and K0sp
+  AliNormalizationCounter *fCounter; // AliNormalizationCounter on output slot 4
+  AliRDHFCutsLctoV0 *fAnalCuts;      // Cuts - sent to output slot 5
+  TList *fListCuts;                  // list of cuts
+  Bool_t fUseOnTheFlyV0;             // flag to analyze also on-the-fly V0 candidates
+  Bool_t fIsEventSelected;           // flag for event selected
+
+  TTree   *fVariablesTreeSgn;        //! tree of the candidate variables after track selection (Signal)
+  TTree   *fVariablesTreeBkg;        //! tree of the candidate variables after track selection (Background)
+  Float_t *fCandidateVariables;      //! variables to be written to the tree
+
+  Bool_t fIspA;                       // flag for running on pA
+
+  TH1F* fHistoEvents;                 // histogram with number of events analyzed
+  TH1F* fHistoLc;                     // histogram with number of Lc
+  TH1F* fHistoLcOnTheFly;             // histogram with number of Lc with on-the-fly V0
+  Bool_t fFillOnlySgn;                // flag to fill only signal (speeding up processing)
+  TH1F* fHistoLcBeforeCuts;           // histogram with number of Lc before any cut 
+  TH1F* fHistoFiducialAcceptance;     // histogram to check FiducialAcceptance cut
+  TH2F* fHistoCodesSgn;               // histogram with codes for bachelor and V0 for signal
+  TH2F* fHistoCodesBkg;               // histogram with codes for bachelor and V0 for background
+  TH1F* fHistoLcpKpiBeforeCuts;       // histogram number of true Lc-->pKpi (3 prong) before any cut
+  AliAODVertex *fVtx1;                // primary vertex
+
+  TH1D* fHistoDistanceLcToPrimVtx;    // KF: distance Lc vertex from primary vertex   
+  TH1D* fHistoDistanceV0ToPrimVtx;    // KF: distance V0 vertex from primary vertex   
+  TH1D* fHistoDistanceV0ToLc;         // KF: distance V0 vertex from Lc vertex    
+
+  TH1D* fHistoDistanceLcToPrimVtxSgn; // KF: distance of signal Lc vertex from primary vertex    
+  TH1D* fHistoDistanceV0ToPrimVtxSgn; // KF: distance for signal Lc of V0 vertex from primary vertex   
+  TH1D* fHistoDistanceV0ToLcSgn;      // KF: distance for signal Lc of V0 vertex from Lc vertex 
+         
+  TH1D* fHistoVtxLcResidualToPrimVtx; // KF: residual wrt MC of distance Lc vertex from primary vertex (MC - KF)
+  TH1D* fHistoVtxV0ResidualToPrimVtx; // KF: residual wrt MC of distance V0 vertex from primary vertex (MC - KF)
+  TH1D* fHistoVtxV0ResidualToLc;      // KF: residual wrt MC of distance V0 vertex from Lc vertex (MC - KF)
+
+  TH1D* fHistoMassV0All;              // KF: mass for all V0 reconstructed with KF
+  TH1D* fHistoDecayLengthV0All;       // KF: decay length for all V0 reconstructed with KF
+  TH1D* fHistoLifeTimeV0All;          // KF: life time for all V0 reconstructed with KF
+
+  TH1D* fHistoMassV0True;             // KF: mass for true V0 reconstructed with KF
+  TH1D* fHistoDecayLengthV0True;      // KF: decay length for true V0 reconstructed with KF
+  TH1D* fHistoLifeTimeV0True;         // KF: life time for true V0 reconstructed with KF
+
+  TH1D* fHistoMassV0TrueFromAOD;      // KF: AOD mass for true V0 reconstructed with KF
+
+  TH1D* fHistoMassV0TrueK0S;          // KF: mass for true V0 which are really K0S reconstructed with KF
+  TH1D* fHistoDecayLengthV0TrueK0S;   // KF: decay length for true V0 which are really K0S reconstructed with KF
+  TH1D* fHistoLifeTimeV0TrueK0S;      // KF: life time for true V0 which are really K0S reconstructed with KF
+
+  TH1D* fHistoMassV0TrueK0SFromAOD;   // KF: AOD mass for true V0 which are really K0S reconstructed with KF
+
+  TH1D* fHistoMassLcAll;              // KF: mass for all Lc reconstructed with KF
+  TH1D* fHistoDecayLengthLcAll;       // KF: decay length for all Lc reconstructed with KF
+  TH1D* fHistoLifeTimeLcAll;          // KF: life time for all Lc reconstructed with KF
+
+  TH1D* fHistoMassLcTrue;             // KF: mass for true cascades reconstructed with KF
+  TH1D* fHistoDecayLengthLcTrue;      // KF: decay length for true cascades reconstructed with KF
+  TH1D* fHistoLifeTimeLcTrue;         // KF: life time for true cascades reconstructed with KF
+
+  TH1D* fHistoMassLcTrueFromAOD;      // KF: AOD mass for true cascades reconstructed with KF
+
+  TH1D* fHistoMassV0fromLcAll;        // KF: mass of V0 for all cascades reconstructed with KF
+  TH1D* fHistoDecayLengthV0fromLcAll; // KF: decay length of V0 for all cascades reconstructed with KF
+  TH1D* fHistoLifeTimeV0fromLcAll;    // KF: life time of V0 for all cascades reconstructed with KF
+
+  TH1D* fHistoMassV0fromLcTrue;       // KF: mass of V0 for true cascades reconstructed with KF
+  TH1D* fHistoDecayLengthV0fromLcTrue;// KF: decay length of V0 for true cascades reconstructed with KF
+  TH1D* fHistoLifeTimeV0fromLcTrue;   // KF: life time of V0 for true cascades reconstructed with KF
+
+  TH1D* fHistoMassLcSgn;              // KF: mass of signal Lc reconstructed with KF
+  TH1D* fHistoMassLcSgnFromAOD;       // KF: AOD mass of signal Lc reconstructed with KF
+  TH1D* fHistoDecayLengthLcSgn;       // KF: decay length of signal Lc reconstructed with KF
+  TH1D* fHistoLifeTimeLcSgn;          // KF: life time of signal Lc reconstructed with KF
+
+  TH1D* fHistoMassV0fromLcSgn;        // KF: mass of V0 for signal Lc reconstructed with KF
+  TH1D* fHistoDecayLengthV0fromLcSgn; // KF: decay length of V0 for signal Lc reconstructed with KF
+  TH1D* fHistoLifeTimeV0fromLcSgn;    // KF: life time of V0 for signal Lc reconstructed with KF
+
+  TH2D* fHistoKF;                     // KF: V0 code vs Lc code from KF (mass, decaylength, lifetime considered) 
+  TH1D* fHistoKFV0;                   // KF: V0 code from KF (mass, decaylength, lifetime considered) 
+  TH1D* fHistoKFLc;                   // KF: Lc code from KF (mass, decaylength, lifetime considered) 
+
+  TH2D* fHistoMassKFV0;               // KF: mass vs mass error for V0 from KF  
+  TH2D* fHistoDecayLengthKFV0;        // KF: decay length vs decay length error for V0 from KF
+  TH2D* fHistoLifeTimeKFV0;           // KF: life time vs life time error for V0 from KF
+
+  TH2D* fHistoMassKFLc;               // KF: mass vs mass error for Lc from KF
+  TH2D* fHistoDecayLengthKFLc;        // KF: decay length vs decay length error for Lc from KF
+  TH2D* fHistoLifeTimeKFLc;           // KF: life time vs life time error for Lc from KF
+
+  TH2D* fHistoArmenterosPodolanskiV0KF;      // KF: Armeteros-Podolanski plot for all V0 from KF
+  TH2D* fHistoArmenterosPodolanskiV0KFSgn;   // KF: Armeteros-Podolanski plot for V0 from signal Lc from KF
+  TH2D* fHistoArmenterosPodolanskiV0AOD;     // KF: AOD Armeteros-Podolanski plot for all V0 from KF
+  TH2D* fHistoArmenterosPodolanskiV0AODSgn;  // KF: AOD Armeteros-Podolanski plot for V0 from signal Lc from KF
+
+  TList *fOutputKF;             // User output1: list of histograms from KF
+
+  Int_t fmcLabelLc;             // label of candidate    
+  Bool_t ftopoConstraint;       // flag to use topological constraints in KF
+  Bool_t fCallKFVertexing;      // flag to decide whether to call or not KF
+  Bool_t fKeepingOnlyHIJINGBkg; // flag to fill bkg with only candidates that have daughters generated by HIJING (to be used for enriched MC)
+  AliVertexingHFUtils* fUtils;  // AliVertexingHFUtils used to check the generator of a specific candidate
+  TH1F* fHistoBackground;       // histo to check the number of candidates with at least one daughter for the injected signal
+
+  ClassDef(AliAnalysisTaskSELc2V0bachelorTMVA,2); // class for Lc->p K0
+};
+
+#endif
+
diff --git a/PWGHF/vertexingHF/macros/AddTaskLc2V0bachpA_TMVA_KF_ON_OnlyHIJINGBkg.C b/PWGHF/vertexingHF/macros/AddTaskLc2V0bachpA_TMVA_KF_ON_OnlyHIJINGBkg.C
new file mode 100644 (file)
index 0000000..6790a41
--- /dev/null
@@ -0,0 +1,82 @@
+AliAnalysisTaskSELc2V0bachelorTMVA* AddTaskLc2V0bachpA_TMVA_KF_ON_OnlyHIJINGBkg(TString finname="Lc2V0bachelorCuts.root",
+                                                                                 Bool_t theMCon=kTRUE,
+                                                                                 Bool_t onTheFly=kFALSE){
+
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+  if (!mgr) {
+    ::Error("AddTaskLc2V0bachelor", "No analysis manager to connect to.");
+    return NULL;
+  }  
+  
+  // cuts are stored in a TFile generated by makeTFile4CutsLc2V0bachelor.C in ./macros/
+  // set there the cuts!!!!!
+
+  Bool_t stdcuts=kFALSE;
+  TFile* filecuts;
+  if( finname.EqualTo("") ) {
+    stdcuts=kTRUE; 
+  } else {
+    filecuts=TFile::Open(finname.Data());
+    if(!filecuts ||(filecuts&& !filecuts->IsOpen())){
+      AliFatal("Input file not found : check your cut object");
+    }
+  }
+
+  AliRDHFCutsLctoV0* RDHFCutsLctoV0anal = new AliRDHFCutsLctoV0();
+  if (stdcuts) RDHFCutsLctoV0anal->SetStandardCutsPP2010();
+  else RDHFCutsLctoV0anal = (AliRDHFCutsLctoV0*)filecuts->Get("LctoV0AnalysisCuts");
+  RDHFCutsLctoV0anal->SetName("LctoV0AnalysisCuts");
+  RDHFCutsLctoV0anal->SetMinPtCandidate(-1.);
+  RDHFCutsLctoV0anal->SetMaxPtCandidate(10000.);
+
+
+  // mm let's see if everything is ok
+  if (!RDHFCutsLctoV0anal) {
+    cout << "Specific AliRDHFCutsLctoV0 not found\n";
+    return;
+  }
+
+
+  //CREATE THE TASK
+
+  printf("CREATE TASK\n");
+  AliAnalysisTaskSELc2V0bachelorTMVA *task = new AliAnalysisTaskSELc2V0bachelorTMVA("AliAnalysisTaskSELc2V0bachelorTMVA", RDHFCutsLctoV0anal, onTheFly);
+  task->SetIspA(kTRUE);
+  task->SetMC(theMCon);
+  task->SetK0sAnalysis(kTRUE);
+  task->SetDebugLevel(0);
+  task->SetCallKFVertexing(kTRUE);
+  task->SetKeepingKeepingOnlyHIJINGBkg(kTRUE);
+  mgr->AddTask(task);
+
+  // Create and connect containers for input/output  
+  TString outputfile = AliAnalysisManager::GetCommonFileName();
+  outputfile += ":PWG3_D2H_Lc2pK0S";
+
+  mgr->ConnectInput(task,0,mgr->GetCommonInputContainer());
+  AliAnalysisDataContainer *coutput1   = mgr->CreateContainer(Form("treeList"),TList::Class(),AliAnalysisManager::kOutputContainer,"Lc2K0Sp_tree_pA.root"); // trees
+
+  mgr->ConnectOutput(task, 1, coutput1);
+
+  AliAnalysisDataContainer *coutputLc2 = mgr->CreateContainer(Form("Lc2pK0Scounter"),AliNormalizationCounter::Class(),AliAnalysisManager::kOutputContainer,"Lc2K0Sp_tree_pA.root"); //counter
+  mgr->ConnectOutput(task, 2, coutputLc2);
+
+  AliAnalysisDataContainer *coutputLc3 = mgr->CreateContainer(Form("Lc2pK0SCuts"),TList::Class(),AliAnalysisManager::kOutputContainer,"Lc2K0Sp_tree_pA.root"); // cuts
+  mgr->ConnectOutput(task, 3, coutputLc3);
+
+  AliAnalysisDataContainer *coutput4   = mgr->CreateContainer(Form("treeSgn"),TTree::Class(),AliAnalysisManager::kOutputContainer,"Lc2K0Sp_tree_pA.root"); // trees
+
+  mgr->ConnectOutput(task, 4, coutput4);
+
+  AliAnalysisDataContainer *coutput5   = mgr->CreateContainer(Form("treeBkg"),TTree::Class(),AliAnalysisManager::kOutputContainer,"Lc2K0Sp_tree_pA.root"); // trees
+
+  mgr->ConnectOutput(task, 5, coutput5);
+
+  AliAnalysisDataContainer *coutput6   = mgr->CreateContainer(Form("listHistoKF"), TList::Class(), AliAnalysisManager::kOutputContainer,"Lc2K0Sp_tree_pA.root"); // trees
+
+  mgr->ConnectOutput(task, 6, coutput6);
+
+
+  return task;
+
+}