PWGHF/vertexingHF/AliAnalysisTaskSELc2V0bachelorTMVA.cxx

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