PWGHF/vertexingHF/AliAnalysisTaskSELc2V0bachelorTMVA.cxx

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