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