PWGHF/vertexingHF/AliAnalysisTaskSELc2V0bachelorTMVA.cxx

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