PWGHF/vertexingHF/AliAnalysisTaskSELc2V0bachelorTMVA.cxx

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