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