AliPhysics  master (3d17d9d)
AliAnalysisTaskSELc2V0bachelorTMVAApp.cxx
Go to the documentation of this file.
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: AliAnalysisTaskSELc2V0bachelorTMVAApp.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, A. Alici
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 <TProfile.h>
36 #include <TH2F.h>
37 #include <TTree.h>
38 #include "TROOT.h"
39 #include <TDatabasePDG.h>
40 #include <TMath.h>
41 #include <TVector3.h>
42 #include <TLorentzVector.h>
43 
44 #include <AliAnalysisDataSlot.h>
45 #include <AliAnalysisDataContainer.h>
46 #include "AliMCEvent.h"
47 #include "AliAnalysisManager.h"
48 #include "AliAODMCHeader.h"
49 #include "AliAODHandler.h"
50 #include "AliLog.h"
51 #include "AliAODVertex.h"
52 #include "AliAODRecoDecay.h"
53 #include "AliAODRecoDecayHF.h"
54 #include "AliAODRecoCascadeHF.h"
55 #include "AliAnalysisVertexingHF.h"
56 #include "AliESDtrack.h"
57 #include "AliAODTrack.h"
58 #include "AliAODv0.h"
59 #include "AliAODMCParticle.h"
60 #include "AliAnalysisTaskSE.h"
63 #include "AliAODPidHF.h"
64 #include "AliPIDResponse.h"
65 #include "AliPIDCombined.h"
66 #include "AliTOFPIDResponse.h"
67 #include "AliInputEventHandler.h"
69 #include "AliKFParticle.h"
70 #include "AliKFVertex.h"
71 #include "AliExternalTrackParam.h"
72 #include "AliESDUtils.h"
73 #include "AliDataFile.h"
74 #include <TMVA/Tools.h>
75 #include <TMVA/Reader.h>
76 #include <TMVA/MethodCuts.h>
77 
78 #include "IClassifierReader.h"
79 
80 using std::cout;
81 using std::endl;
82 
83 #include <dlfcn.h>
84 
88 
89 //__________________________________________________________________________
92  fUseMCInfo(kFALSE),
93  fOutput(0),
94  fPIDResponse(0),
95  fPIDCombined(0),
96  fIsK0sAnalysis(kFALSE),
97  fCounter(0),
98  fCounterC(0),
99  fAnalCuts(0),
100  fListCuts(0),
101  fListWeight(0),
102  fListCounters(0),
103  fListProfiles(0),
104  fUseOnTheFlyV0(kFALSE),
105  fIsEventSelected(kFALSE),
106  fVariablesTreeSgn(0),
107  fVariablesTreeBkg(0),
108  fCandidateVariables(),
109  fHistoCentrality(0),
110  fHistoEvents(0),
111  fHistoTracklets_1(0),
112  fHistoTracklets_1_cent(0),
113  fHistoTracklets_All(0),
114  fHistoTracklets_All_cent(0),
115  fHistoLc(0),
116  fHistoLcOnTheFly(0),
117  fFillOnlySgn(kFALSE),
118  fHistoLcBeforeCuts(0),
119  fHistoFiducialAcceptance(0),
120  fHistoCodesSgn(0),
121  fHistoCodesBkg(0),
122  fHistoLcpKpiBeforeCuts(0),
123  fVtx1(0),
124  fHistoDistanceLcToPrimVtx(0),
125  fHistoDistanceV0ToPrimVtx(0),
126  fHistoDistanceV0ToLc(0),
127  fHistoDistanceLcToPrimVtxSgn(0),
128  fHistoDistanceV0ToPrimVtxSgn(0),
129  fHistoDistanceV0ToLcSgn(0),
130  fHistoVtxLcResidualToPrimVtx(0),
131  fHistoVtxV0ResidualToPrimVtx(0),
132  fHistoVtxV0ResidualToLc(0),
133  fHistoMassV0All(0),
134  fHistoDecayLengthV0All(0),
135  fHistoLifeTimeV0All(0),
136  fHistoMassV0True(0),
137  fHistoDecayLengthV0True(0),
138  fHistoLifeTimeV0True(0),
139  fHistoMassV0TrueFromAOD(0),
140  fHistoMassV0TrueK0S(0),
141  fHistoDecayLengthV0TrueK0S(0),
142  fHistoLifeTimeV0TrueK0S(0),
143  fHistoMassV0TrueK0SFromAOD(0),
144  fHistoMassLcAll(0),
145  fHistoDecayLengthLcAll(0),
146  fHistoLifeTimeLcAll(0),
147  fHistoMassLcTrue(0),
148  fHistoDecayLengthLcTrue(0),
149  fHistoLifeTimeLcTrue(0),
150  fHistoMassLcTrueFromAOD(0),
151  fHistoMassV0fromLcAll(0),
152  fHistoDecayLengthV0fromLcAll(0),
153  fHistoLifeTimeV0fromLcAll(0),
154  fHistoMassV0fromLcTrue(0),
155  fHistoDecayLengthV0fromLcTrue(0),
156  fHistoLifeTimeV0fromLcTrue(0),
157  fHistoMassLcSgn(0),
158  fHistoMassLcSgnFromAOD(0),
159  fHistoDecayLengthLcSgn(0),
160  fHistoLifeTimeLcSgn(0),
161  fHistoMassV0fromLcSgn(0),
162  fHistoDecayLengthV0fromLcSgn(0),
163  fHistoLifeTimeV0fromLcSgn(0),
164  fHistoKF(0),
165  fHistoKFV0(0),
166  fHistoKFLc(0),
167  fHistoMassKFV0(0),
168  fHistoDecayLengthKFV0(0),
169  fHistoLifeTimeKFV0(0),
170  fHistoMassKFLc(0),
171  fHistoDecayLengthKFLc(0),
172  fHistoLifeTimeKFLc(0),
173  fHistoArmenterosPodolanskiV0KF(0),
174  fHistoArmenterosPodolanskiV0KFSgn(0),
175  fHistoArmenterosPodolanskiV0AOD(0),
176  fHistoArmenterosPodolanskiV0AODSgn(0),
177  fOutputKF(0),
178  fmcLabelLc(-1),
179  ftopoConstraint(kTRUE),
180  fCallKFVertexing(kFALSE),
181  fKeepingOnlyHIJINGBkg(kFALSE),
182  fUtils(0),
183  fHistoBackground(0),
184  fCutKFChi2NDF(999999.),
185  fCutKFDeviationFromVtx(999999.),
186  fCutKFDeviationFromVtxV0(0.),
187  fCurrentEvent(-1),
188  fBField(0),
189  fKeepingOnlyPYTHIABkg(kFALSE),
190  fHistoMCLcK0SpGen(0),
191  fHistoMCLcK0SpGenAcc(0),
192  fHistoMCLcK0SpGenLimAcc(0),
193  fTriggerMask(0),
194  fFuncWeightPythia(0),
195  fFuncWeightFONLL5overLHC13d3(0),
196  fFuncWeightFONLL5overLHC13d3Lc(0),
197  fHistoMCNch(0),
198  fNTracklets_1(0),
199  fNTracklets_All(0),
200  fCentrality(0),
201  fFillTree(0),
202  fUseWeightsLibrary(kFALSE),
203  fBDTReader(0),
204  fTMVAlibName(""),
205  fTMVAlibPtBin(""),
206  fNamesTMVAVar(""),
207  fBDTHisto(0),
208  fBDTHistoVsMassK0S(0),
209  fBDTHistoVstImpParBach(0),
210  fBDTHistoVstImpParV0(0),
211  fBDTHistoVsBachelorPt(0),
212  fBDTHistoVsCombinedProtonProb(0),
213  fBDTHistoVsCtau(0),
214  fBDTHistoVsCosPAK0S(0),
215  fBDTHistoVsSignd0(0),
216  fBDTHistoVsCosThetaStar(0),
217  fBDTHistoVsnSigmaTPCpr(0),
218  fBDTHistoVsnSigmaTOFpr(0),
219  fBDTHistoVsnSigmaTPCpi(0),
220  fBDTHistoVsnSigmaTPCka(0),
221  fBDTHistoVsBachelorP(0),
222  fBDTHistoVsBachelorTPCP(0),
223  fHistoNsigmaTPC(0),
224  fHistoNsigmaTOF(0),
225  fDebugHistograms(kFALSE),
226  fAODProtection(1),
227  fUsePIDresponseForNsigma(kFALSE),
228  fNVars(14),
229  fTimestampCut(0),
230  fUseXmlWeightsFile(kTRUE),
231  fReader(0),
232  fVarsTMVA(0),
233  fNVarsSpectators(0),
234  fVarsTMVASpectators(0),
235  fXmlWeightsFile(""),
236  fBDTHistoTMVA(0),
237  fRefMult(9.26),
238  fYearNumber(16),
239  fHistoNtrUnCorr(0),
240  fHistoNtrCorr(0),
241  fHistoVzVsNtrUnCorr(0),
242  fHistoVzVsNtrCorr(0),
243  fUseMultCorrection(kFALSE),
244  fMultiplicityEstimator(kNtrk10),
245  fDoVZER0ParamVertexCorr(1),
246  fUseMultiplicityCut(kFALSE),
247  fMultiplicityCutMin(0.),
248  fMultiplicityCutMax(99999.),
249  fUseXmlFileFromCVMFS(kFALSE),
250  fXmlFileFromCVMFS("")
251 {
253  //
254  for(Int_t i=0; i<14; i++) fMultEstimatorAvg[i]=0;
255 }
256 //___________________________________________________________________________
258  AliRDHFCutsLctoV0* analCuts, Bool_t useOnTheFly) :
259  AliAnalysisTaskSE(name),
260  fUseMCInfo(kFALSE),
261  fOutput(0),
262  fPIDResponse(0),
263  fPIDCombined(0),
264  fIsK0sAnalysis(kFALSE),
265  fCounter(0),
266  fCounterC(0),
267  fAnalCuts(analCuts),
268  fListCuts(0),
269  fListWeight(0),
270  fListCounters(0),
271  fListProfiles(0),
272  fUseOnTheFlyV0(useOnTheFly),
273  fIsEventSelected(kFALSE),
277  fHistoCentrality(0),
278  fHistoEvents(0),
283  fHistoLc(0),
284  fHistoLcOnTheFly(0),
285  fFillOnlySgn(kFALSE),
288  fHistoCodesSgn(0),
289  fHistoCodesBkg(0),
291  fVtx1(0),
301  fHistoMassV0All(0),
304  fHistoMassV0True(0),
312  fHistoMassLcAll(0),
315  fHistoMassLcTrue(0),
325  fHistoMassLcSgn(0),
332  fHistoKF(0),
333  fHistoKFV0(0),
334  fHistoKFLc(0),
335  fHistoMassKFV0(0),
338  fHistoMassKFLc(0),
345  fOutputKF(0),
346  fmcLabelLc(-1),
347  ftopoConstraint(kTRUE),
348  fCallKFVertexing(kFALSE),
349  fKeepingOnlyHIJINGBkg(kFALSE),
350  fUtils(0),
351  fHistoBackground(0),
352  fCutKFChi2NDF(999999.),
353  fCutKFDeviationFromVtx(999999.),
355  fCurrentEvent(-1),
356  fBField(0),
357  fKeepingOnlyPYTHIABkg(kFALSE),
361  fTriggerMask(0),
365  fHistoMCNch(0),
366  fNTracklets_1(0),
367  fNTracklets_All(0),
368  fCentrality(0),
369  fFillTree(0),
370  fUseWeightsLibrary(kFALSE),
371  fBDTReader(0),
372  fTMVAlibName(""),
373  fTMVAlibPtBin(""),
374  fNamesTMVAVar(""),
375  fBDTHisto(0),
381  fBDTHistoVsCtau(0),
391  fHistoNsigmaTPC(0),
392  fHistoNsigmaTOF(0),
393  fDebugHistograms(kFALSE),
394  fAODProtection(1),
395  fUsePIDresponseForNsigma(kFALSE),
396  fNVars(14),
397  fTimestampCut(0),
398  fUseXmlWeightsFile(kTRUE),
399  fReader(0),
400  fVarsTMVA(0),
401  fNVarsSpectators(0),
404  fXmlWeightsFile(""),
405  fBDTHistoTMVA(0),
406  fRefMult(9.26),
407  fYearNumber(16),
408  fHistoNtrUnCorr(0),
409  fHistoNtrCorr(0),
412  fUseMultCorrection(kFALSE),
415  fUseMultiplicityCut(kFALSE),
417  fMultiplicityCutMax(99999.),
418  fUseXmlFileFromCVMFS(kFALSE),
420 {
421  //
423  //
424  Info("AliAnalysisTaskSELc2V0bachelorTMVAApp","Calling Constructor");
425 
426  for(Int_t i=0; i<14; i++) fMultEstimatorAvg[i]=0;
427 
428  DefineOutput(1, TList::Class()); // Tree signal + Tree Bkg + histoEvents
429  DefineOutput(2, TList::Class()); // normalization counter object
430  DefineOutput(3, TList::Class()); // Cuts
431  DefineOutput(4, TTree::Class()); // Tree signal + Tree Bkg + histoEvents
432  DefineOutput(5, TTree::Class()); // Tree signal + Tree Bkg + histoEvents
433  DefineOutput(6, TList::Class()); // Tree signal + Tree Bkg + histoEvents
434  DefineOutput(7, TList::Class()); // weights
435 }
436 
437 //___________________________________________________________________________
439  //
441  //
442  Info("~AliAnalysisTaskSELc2V0bachelorTMVAApp","Calling Destructor");
443 
444  if (fOutput) {
445  delete fOutput;
446  fOutput = 0;
447  }
448 
449  if (fPIDResponse) {
450  delete fPIDResponse;
451  }
452 
453  if (fPIDCombined) {
454  delete fPIDCombined;
455  }
456 
457  if (fCounter) {
458  delete fCounter;
459  fCounter = 0;
460  }
461 
462  if (fCounterC) {
463  delete fCounterC;
464  fCounterC = 0;
465  }
466 
467  if (fAnalCuts) {
468  delete fAnalCuts;
469  fAnalCuts = 0;
470  }
471 
472  if (fListCuts) {
473  delete fListCuts;
474  fListCuts = 0;
475  }
476 
477  if (fListCounters) {
478  delete fListCounters;
479  fListCounters = 0;
480  }
481 
482  if (fListWeight) {
483  delete fListWeight;
484  fListWeight = 0;
485  }
486 
487  if(fVariablesTreeSgn){
488  delete fVariablesTreeSgn;
489  fVariablesTreeSgn = 0;
490  }
491 
492  if(fVariablesTreeBkg){
493  delete fVariablesTreeBkg;
494  fVariablesTreeBkg = 0;
495  }
496 
497  if (fOutputKF) {
498  delete fOutputKF;
499  fOutputKF = 0;
500  }
501 
502  if (fUtils) {
503  delete fUtils;
504  fUtils = 0;
505  }
506 
507  if (fBDTReader) {
508  //delete fBDTReader;
509  fBDTReader = 0;
510  }
511 
512  if (fReader) {
513  delete fReader;
514  fReader = 0;
515  }
516 
517  if (fVarsTMVA) {
518  delete fVarsTMVA;
519  fVarsTMVA = 0;
520  }
521 
522  if (fVarsTMVASpectators) {
523  delete fVarsTMVASpectators;
525  }
526 
527  for(Int_t i=0; i<14; i++){
528  if (fMultEstimatorAvg[i]) delete fMultEstimatorAvg[i];
529  }
530 }
531 //_________________________________________________
533  //
535  //
536 
537  fIsEventSelected=kFALSE;
538 
539  if (fDebug > 1) AliInfo("Init");
540 
541  fListCuts = new TList();
542  fListCuts->SetOwner();
544  PostData(3, fListCuts);
545 
546  // Save the weight functions or histograms
547  fListWeight = new TList();
548  fListWeight->SetOwner();
549  fListWeight->Add(fHistoMCNch);
550  PostData(7, fListWeight);
551 
553 
554  return;
555 }
556 
557 //________________________________________ terminate ___________________________
559 {
563 
564  AliInfo("Terminate");
565  AliAnalysisTaskSE::Terminate();
566 
567  fOutput = dynamic_cast<TList*> (GetOutputData(1));
568  if (!fOutput) {
569  AliError("fOutput not available");
570  return;
571  }
572 
573  if(fHistoMCLcK0SpGen) {
574  AliInfo(Form("At MC level, %f Lc --> K0S + p were found", fHistoMCLcK0SpGen->GetEntries()));
575  } else {
576  AliInfo("fHistoMCLcK0SpGen not available");
577  }
579  AliInfo(Form("At MC level, %f Lc --> K0S + p were found in the acceptance", fHistoMCLcK0SpGenAcc->GetEntries()));
580  } else {
581  AliInfo("fHistoMCLcK0SpGenAcc not available");
582  }
583  if(fVariablesTreeSgn) {
584  AliInfo(Form("At Reco level, %lld Lc --> K0S + p were found", fVariablesTreeSgn->GetEntries()));
585  } else {
586  AliInfo("fVariablesTreeSgn not available");
587  }
588 
589  fOutputKF = dynamic_cast<TList*> (GetOutputData(6));
590  if (!fOutputKF) {
591  AliError("fOutputKF not available");
592  return;
593  }
594 
595  return;
596 }
597 
598 //___________________________________________________________________________
601  AliInfo(Form("CreateOutputObjects of task %s\n", GetName()));
602 
603  fOutput = new TList();
604  fOutput->SetOwner();
605  fOutput->SetName("listTrees");
606 
607  // Output slot 1: list of 2 trees (Sgn + Bkg) of the candidate variables
608  const char* nameoutput = GetOutputSlot(1)->GetContainer()->GetName();
609  fVariablesTreeSgn = new TTree(Form("%s_Sgn", nameoutput), "Candidates variables tree, Signal");
610  fVariablesTreeBkg = new TTree(Form("%s_Bkg", nameoutput), "Candidates variables tree, Background");
611 
612  Int_t nVar;
613  if (fUseMCInfo) nVar = 52; //"full" tree if MC
614  else nVar = 35; //"reduced" tree if data
615 
616  fCandidateVariables = new Float_t [nVar];
617  TString * fCandidateVariableNames = new TString[nVar];
618 
619  if (fUseMCInfo) { // "full tree" for MC
620  fCandidateVariableNames[0] = "massLc2K0Sp";
621  fCandidateVariableNames[1] = "massLc2Lambdapi";
622  fCandidateVariableNames[2] = "massK0S";
623  fCandidateVariableNames[3] = "massLambda";
624  fCandidateVariableNames[4] = "massLambdaBar";
625  fCandidateVariableNames[5] = "cosPAK0S";
626  fCandidateVariableNames[6] = "dcaV0";
627  fCandidateVariableNames[7] = "tImpParBach";
628  fCandidateVariableNames[8] = "tImpParV0";
629  fCandidateVariableNames[9] = "nSigmaTPCpr";
630  fCandidateVariableNames[10] = "nSigmaTOFpr";
631  fCandidateVariableNames[11] = "bachelorPt";
632  fCandidateVariableNames[12] = "V0positivePt";
633  fCandidateVariableNames[13] = "V0negativePt";
634  fCandidateVariableNames[14] = "dcaV0pos";
635  fCandidateVariableNames[15] = "dcaV0neg";
636  fCandidateVariableNames[16] = "v0Pt";
637  fCandidateVariableNames[17] = "massGamma";
638  fCandidateVariableNames[18] = "LcPt";
639  fCandidateVariableNames[19] = "combinedProtonProb";
640  fCandidateVariableNames[20] = "LcEta";
641  fCandidateVariableNames[21] = "V0positiveEta";
642  fCandidateVariableNames[22] = "V0negativeEta";
643  fCandidateVariableNames[23] = "TPCProtonProb";
644  fCandidateVariableNames[24] = "TOFProtonProb";
645  fCandidateVariableNames[25] = "bachelorEta";
646  fCandidateVariableNames[26] = "LcP";
647  fCandidateVariableNames[27] = "bachelorP";
648  fCandidateVariableNames[28] = "v0P";
649  fCandidateVariableNames[29] = "V0positiveP";
650  fCandidateVariableNames[30] = "V0negativeP";
651  fCandidateVariableNames[31] = "v0Eta";
652  fCandidateVariableNames[32] = "LcPtMC";
653  fCandidateVariableNames[33] = "DecayLengthK0S";
654  fCandidateVariableNames[34] = "bachCode";
655  fCandidateVariableNames[35] = "k0SCode";
656  fCandidateVariableNames[36] = "alphaArm";
657  fCandidateVariableNames[37] = "ptArm";
658  fCandidateVariableNames[38] = "CosThetaStar";
659  fCandidateVariableNames[39] = "weightPtFlat";
660  fCandidateVariableNames[40] = "weightFONLL5overLHC13d3";
661  fCandidateVariableNames[41] = "weightFONLL5overLHC13d3Lc";
662  fCandidateVariableNames[42] = "weightNch";
663  fCandidateVariableNames[43] = "NtrkRaw";
664  fCandidateVariableNames[44] = "NtrkCorr";
665  fCandidateVariableNames[45] = "signd0";
666  fCandidateVariableNames[46] = "centrality";
667  fCandidateVariableNames[47] = "NtrkAll";
668  fCandidateVariableNames[48] = "origin";
669  fCandidateVariableNames[49] = "nSigmaTPCpi";
670  fCandidateVariableNames[50] = "nSigmaTPCka";
671  fCandidateVariableNames[51] = "bachTPCmom";
672  }
673  else { // "light mode"
674  fCandidateVariableNames[0] = "massLc2K0Sp";
675  fCandidateVariableNames[1] = "alphaArm";
676  fCandidateVariableNames[2] = "massK0S";
677  fCandidateVariableNames[3] = "massLambda";
678  fCandidateVariableNames[4] = "massLambdaBar";
679  fCandidateVariableNames[5] = "cosPAK0S";
680  fCandidateVariableNames[6] = "dcaV0";
681  fCandidateVariableNames[7] = "tImpParBach";
682  fCandidateVariableNames[8] = "tImpParV0";
683  fCandidateVariableNames[9] = "nSigmaTPCpr";
684  fCandidateVariableNames[10] = "nSigmaTOFpr";
685  fCandidateVariableNames[11] = "bachelorPt";
686  fCandidateVariableNames[12] = "V0positivePt";
687  fCandidateVariableNames[13] = "V0negativePt";
688  fCandidateVariableNames[14] = "dcaV0pos";
689  fCandidateVariableNames[15] = "dcaV0neg";
690  fCandidateVariableNames[16] = "v0Pt";
691  fCandidateVariableNames[17] = "bachTPCmom";
692  fCandidateVariableNames[18] = "LcPt";
693  fCandidateVariableNames[19] = "combinedProtonProb";
694  fCandidateVariableNames[20] = "V0positiveEta";
695  fCandidateVariableNames[21] = "bachelorP"; // we replaced the V0negativeEta with the bachelor P as this is more useful (for PID) while the V0 daughters' eta we don't use... And are practically the same (positive and negative)
696  fCandidateVariableNames[22] = "bachelorEta";
697  fCandidateVariableNames[23] = "v0P";
698  fCandidateVariableNames[24] = "DecayLengthK0S";
699  fCandidateVariableNames[25] = "nSigmaTPCpi";
700  fCandidateVariableNames[26] = "nSigmaTPCka";
701  fCandidateVariableNames[27] = "NtrkRaw";
702  fCandidateVariableNames[28] = "NtrkCorr";
703  fCandidateVariableNames[29] = "CosThetaStar";
704  fCandidateVariableNames[30] = "signd0";
705  fCandidateVariableNames[31] = "centrality";
706  fCandidateVariableNames[32] = "NtrkAll";
707  fCandidateVariableNames[33] = "origin";
708  fCandidateVariableNames[34] = "ptArm";
709  }
710 
711  for(Int_t ivar=0; ivar < nVar; ivar++){
712  fVariablesTreeSgn->Branch(fCandidateVariableNames[ivar].Data(), &fCandidateVariables[ivar], Form("%s/f",fCandidateVariableNames[ivar].Data()));
713  fVariablesTreeBkg->Branch(fCandidateVariableNames[ivar].Data(), &fCandidateVariables[ivar], Form("%s/f",fCandidateVariableNames[ivar].Data()));
714  }
715 
716  fHistoCentrality = new TH1F("fHistoCentrality", "fHistoCentrality", 100, 0., 100.);
717 
718  fHistoEvents = new TH1F("fHistoEvents", "fHistoEvents", 6, -0.5, 5.5);
719  TString labelEv[6] = {"RejectedDeltaMismatch", "AcceptedDeltaMismatch", "NotSelected", "TimeStampCut", "Selected","AcceptedMultCut"};
720  for (Int_t ibin = 1; ibin <= fHistoEvents->GetNbinsX(); ibin++){
721  fHistoEvents->GetXaxis()->SetBinLabel(ibin, labelEv[ibin-1].Data());
722  }
723 
724  fHistoTracklets_1 = new TH1F("fHistoTracklets_1", "fHistoTracklets_1", 1000, 0, 5000);
725  fHistoTracklets_1_cent = new TH2F("fHistoTracklets_1_cent", "fHistoTracklets_1_cent; centrality; SPD tracklets [-1, 1]", 100, 0., 100., 1000, 0, 5000);
726  fHistoTracklets_All = new TH1F("fHistoTracklets_All", "fHistoTracklets_All", 1000, 0, 5000);
727  fHistoTracklets_All_cent = new TH2F("fHistoTracklets_All_cent", "fHistoTracklets_All_cent; centrality; SPD tracklets [-999, 999]", 100, 0., 100., 1000, 0, 5000);
728 
729  fHistoLc = new TH1F("fHistoLc", "fHistoLc", 2, -0.5, 1.5);
730 
731  fHistoLcOnTheFly = new TH1F("fHistoLcOnTheFly", "fHistoLcOnTheFly", 4, -0.5, 3.5);
732  TString labelOnTheFly[4] = {"OnTheFly-Bkg", "OfflineBkg", "OnTheFly-Sgn", "OfflineSgn"};
733  for (Int_t ibin = 1; ibin <= fHistoLcOnTheFly->GetNbinsX(); ibin++){
734  fHistoLcOnTheFly->GetXaxis()->SetBinLabel(ibin, labelOnTheFly[ibin-1].Data());
735  }
736 
737  fHistoLcBeforeCuts = new TH1F("fHistoLcBeforeCuts", "fHistoLcBeforeCuts", 2, -0.5, 1.5);
738  TString labelBeforeCuts[2] = {"Bkg", "Sgn"};
739  for (Int_t ibin = 1; ibin <= fHistoLcBeforeCuts->GetNbinsX(); ibin++){
740  fHistoLcBeforeCuts->GetXaxis()->SetBinLabel(ibin, labelBeforeCuts[ibin-1].Data());
741  fHistoLc->GetXaxis()->SetBinLabel(ibin, labelBeforeCuts[ibin-1].Data());
742  }
743 
744  fHistoFiducialAcceptance = new TH1F("fHistoFiducialAcceptance", "fHistoFiducialAcceptance", 4, -0.5, 3.5);
745  TString labelAcc[4] = {"NotOk-Bkg", "Ok-Bkg", "NotOk-Sgn", "Ok-Sgn"};
746  for (Int_t ibin = 1; ibin <= fHistoFiducialAcceptance->GetNbinsX(); ibin++){
747  fHistoFiducialAcceptance->GetXaxis()->SetBinLabel(ibin, labelAcc[ibin-1].Data());
748  }
749 
750  fHistoCodesSgn = new TH2F("fHistoCodesSgn", "fHistoCodes for Signal; bachelor; K0S", 7, -1.5, 5.5, 9, -1.5, 7.5);
751  fHistoCodesBkg = new TH2F("fHistoCodesBkg", "fHistoCodes for Background; bachelor; K0S", 7, -1.5, 5.5, 9, -1.5, 7.5);
752 
753  TString labelx[7] = { "kBachInvalid", "kBachFake", "kBachNoProton", "kBachPrimary", "kBachNoLambdaMother",
754  "kBachDifferentLambdaMother", "kBachCorrectLambdaMother"};
755  TString labely[9] = { "kK0SInvalid", "kK0SFake", "kK0SNoK0S", "kK0SWithoutMother", "kK0SNotFromK0",
756  "kK0Primary", "kK0NoLambdaMother", "kK0DifferentLambdaMother", "kK0CorrectLambdaMother"};
757 
758  for (Int_t ibin = 1; ibin <= fHistoCodesSgn->GetNbinsX(); ibin++){
759  fHistoCodesSgn->GetXaxis()->SetBinLabel(ibin, labelx[ibin-1].Data());
760  fHistoCodesBkg->GetXaxis()->SetBinLabel(ibin, labelx[ibin-1].Data());
761  }
762  for (Int_t ibin = 1; ibin <= fHistoCodesSgn->GetNbinsY(); ibin++){
763  fHistoCodesSgn->GetYaxis()->SetBinLabel(ibin, labely[ibin-1].Data());
764  fHistoCodesBkg->GetYaxis()->SetBinLabel(ibin, labely[ibin-1].Data());
765  }
766 
767  fHistoLcpKpiBeforeCuts = new TH1F("fHistoLcpKpiBeforeCuts", "fHistoLcpKpiBeforeCuts", 2, -0.5, 1.5);
768  for (Int_t ibin = 1; ibin <= fHistoLcpKpiBeforeCuts->GetNbinsX(); ibin++){
769  fHistoLcpKpiBeforeCuts->GetXaxis()->SetBinLabel(ibin, labelBeforeCuts[ibin-1].Data());
770  }
771 
772  fHistoBackground = new TH1F("fHistoBackground", "fHistoBackground", 4, -0.5, 3.5);
773  TString labelBkg[4] = {"Injected", "Non-injected", "Non-PYTHIA", "PYTHIA"};
774  for (Int_t ibin = 1; ibin <= fHistoBackground->GetNbinsX(); ibin++){
775  fHistoBackground->GetXaxis()->SetBinLabel(ibin, labelBkg[ibin-1].Data());
776  }
777 
778  const Float_t ptbins[15] = {0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 12., 17., 25., 35.};
779 
780  fHistoMCLcK0SpGen = new TH1F("fHistoMCLcK0SpGen", "fHistoMCLcK0SpGen", 14, ptbins);
781  fHistoMCLcK0SpGenAcc = new TH1F("fHistoMCLcK0SpGenAcc", "fHistoMCLcK0SpGenAcc", 14, ptbins);
782  fHistoMCLcK0SpGenLimAcc = new TH1F("fHistoMCLcK0SpGenLimAcc", "fHistoMCLcK0SpGenLimAcc", 14, ptbins);
783 
784 
785  // //code to run the BDT Application on-the-fly
786  // TString inputVariablesBDT = "massK0S,tImpParBach,tImpParV0,bachelorPt,combinedProtonProb,DecayLengthK0S*0.497/v0P,cosPAK0S,CosThetaStar,signd0";
787  // TObjArray *tokens = inputVariablesBDT.Tokenize(",");
788  // tokens->Print();
789  // std::vector<std::string> inputNamesVec;
790  // for(Int_t i=0; i<tokens->GetEntries(); i++){
791  // TString variable = ((TObjString*)(tokens->At(i)))->String();
792  // string tmpvar = variable.Data();
793  // inputNamesVec.push_back(tmpvar);
794  // }
795  // Printf("************************************************ fBDTReader = %p", fBDTReader);
796  // //fBDTReader = new ReadBDT_Default(inputNamesVec);
797 
798 
799  fBDTHisto = new TH2D("fBDTHisto", "Lc inv mass vs bdt output; bdt; m_{inv}(pK^{0}_{S})[GeV/#it{c}^{2}]", 10000, -1, 1, 1000, 2.05, 2.55);
800  fBDTHistoTMVA = new TH2D("fBDTHistoTMVA", "Lc inv mass vs bdt output; bdt; m_{inv}(pK^{0}_{S})[GeV/#it{c}^{2}]", 10000, -1, 1, 1000, 2.05, 2.55);
801  if (fDebugHistograms) {
802  fBDTHistoVsMassK0S = new TH2D("fBDTHistoVsMassK0S", "K0S inv mass vs bdt output; bdt; m_{inv}(#pi^{+}#pi^{#minus})[GeV/#it{c}^{2}]", 1000, -1, 1, 1000, 0.485, 0.51);
803  fBDTHistoVstImpParBach = new TH2D("fBDTHistoVstImpParBach", "d0 bachelor vs bdt output; bdt; d_{0, bachelor}[cm]", 1000, -1, 1, 100, -1, 1);
804  fBDTHistoVstImpParV0 = new TH2D("fBDTHistoVstImpParV0", "d0 K0S vs bdt output; bdt; d_{0, V0}[cm]", 1000, -1, 1, 100, -1, 1);
805  fBDTHistoVsBachelorPt = new TH2D("fBDTHistoVsBachelorPt", "bachelor pT vs bdt output; bdt; p_{T, bachelor}[GeV/#it{c}]", 1000, -1, 1, 100, 0, 20);
806  fBDTHistoVsCombinedProtonProb = new TH2D("fBDTHistoVsCombinedProtonProb", "combined proton probability vs bdt output; bdt; Bayesian PID_{bachelor}", 1000, -1, 1, 100, 0, 1);
807  fBDTHistoVsCtau = new TH2D("fBDTHistoVsCtau", "K0S ctau vs bdt output; bdt; c#tau_{V0}[cm]", 1000, -1, 1, 1000, 0, 100);
808  fBDTHistoVsCosPAK0S = new TH2D("fBDTHistoVsCosPAK0S", "V0 cosine pointing angle vs bdt output; bdt; CosPAK^{0}_{S}", 1000, -1, 1, 100, 0.9, 1);
809  fBDTHistoVsCosThetaStar = new TH2D("fBDTHistoVsCosThetaStar", "proton emission angle in pK0s pair rest frame vs bdt output; bdt; Cos#Theta*", 1000, -1, 1, 100, -1, 1);
810  fBDTHistoVsSignd0 = new TH2D("fBDTHistoVsSignd0", "signed d0 bachelor vs bdt output; bdt; signd_{0, bachelor}[cm]", 1000, -1, 1, 100, -1, 1);
811  fBDTHistoVsnSigmaTPCpr = new TH2D("fBDTHistoVsnSigmaTPCpr", "nSigmaTPCpr vs bdt output; bdt; n_{#sigma}^{TPC}_{pr}", 1000, -1, 1, 1000, -10, 10);
812  fBDTHistoVsnSigmaTOFpr = new TH2D("fBDTHistoVsnSigmaTOFpr", "nSigmaTOFpr vs bdt output; bdt; n_{#sigma}^{TOF}_{pr}", 1000, -1, 1, 1000, -10, 10);
813  fBDTHistoVsnSigmaTPCpi = new TH2D("fBDTHistoVsnSigmaTPCpi", "nSigmaTPCpi vs bdt output; bdt; n_{#sigma}^{TPC}_{pi}", 1000, -1, 1, 1000, -10, 10);
814  fBDTHistoVsnSigmaTPCka = new TH2D("fBDTHistoVsnSigmaTPCka", "nSigmaTPCka vs bdt output; bdt; n_{#sigma}^{TPC}_{ka}", 1000, -1, 1, 1000, -10, 10);
815  fBDTHistoVsBachelorP = new TH2D("fBDTHistoVsBachelorP", "bachelor p vs bdt output; bdt; p_{bachelor}[GeV/#it{c}]", 1000, -1, 1, 100, 0, 20);
816  fBDTHistoVsBachelorTPCP = new TH2D("fBDTHistoVsBachelorTPCP", "bachelor TPC momentum vs bdt output; bdt; p_{TPC, bachelor}[GeV/#it{c}]", 1000, -1, 1, 100, 0, 20);
817  fHistoNsigmaTPC = new TH2D("fHistoNsigmaTPC", "; #it{p} (GeV/#it{c}); n_{#sigma}^{TPC} (proton hypothesis)", 500, 0, 5, 1000, -5, 5);
818  fHistoNsigmaTOF = new TH2D("fHistoNsigmaTOF", "; #it{p} (GeV/#it{c}); n_{#sigma}^{TOF} (proton hypothesis)", 500, 0, 5, 1000, -5, 5);
819  }
820 
821  fOutput->Add(fHistoEvents);
826  fOutput->Add(fHistoLc);
830  fOutput->Add(fHistoCodesSgn);
831  fOutput->Add(fHistoCodesBkg);
838  fOutput->Add(fBDTHisto);
839  fOutput->Add(fBDTHistoTMVA);
840  if (fDebugHistograms) {
846  fOutput->Add(fBDTHistoVsCtau);
856  fOutput->Add(fHistoNsigmaTPC);
857  fOutput->Add(fHistoNsigmaTOF);
858  }
859 
860  if(fUseMultCorrection){
861 
862  Int_t nMultBins = 200;
863  Float_t firstMultBin = -0.5;
864  Float_t lastMultBin = 199.5;
865  const char *estimatorName="tracklets";
867  nMultBins = 400;
868  lastMultBin = 799.5;
869  estimatorName = "vzero";
870  }
871 
872  fHistoNtrUnCorr = new TH1F("fHistoNtrUnCorr", Form("Uncorrected %s mulitplicity; %s; Entries",estimatorName,estimatorName),nMultBins, firstMultBin, lastMultBin);
873  fHistoNtrCorr = new TH1F("fHistoNtrCorr", Form("Corrected %s mulitplicity; %s; Entries",estimatorName,estimatorName),nMultBins, firstMultBin, lastMultBin);
874 
875  fHistoVzVsNtrUnCorr = new TH2F("fHistoVzVsNtrUnCorr", Form("VtxZ vs Uncorrected %s mulitplicity; VtxZ; %s; Entries",estimatorName,estimatorName), 300,-15.,15., nMultBins, firstMultBin, lastMultBin);
876  fHistoVzVsNtrCorr = new TH2F("fHistoVzVsNtrCorr", Form("VtxZ vs Corrected %s mulitplicity; VtxZ; %s; Entries",estimatorName,estimatorName), 300,-15.,15., nMultBins, firstMultBin, lastMultBin);
877 
878  fOutput->Add(fHistoNtrUnCorr);
879  fOutput->Add(fHistoNtrCorr);
882  }
883 
884  PostData(1, fOutput);
885  PostData(4, fVariablesTreeSgn);
886  PostData(5, fVariablesTreeBkg);
887  AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
888  AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
889  fPIDResponse = inputHandler->GetPIDResponse();
890 
891  // fAnalCuts->SetUsePID(kTRUE); // this forces the PID to be used, despite the settings in the cut object
892  // if (fAnalCuts->GetIsUsePID()){
893  // /*
894  // fAnalCuts->GetPidHF()->SetPidResponse(fPIDResponse);
895  // fAnalCuts->GetPidV0pos()->SetPidResponse(fPIDResponse);
896  // fAnalCuts->GetPidV0neg()->SetPidResponse(fPIDResponse);
897  // fAnalCuts->GetPidHF()->SetOldPid(kFALSE);
898  // fAnalCuts->GetPidV0pos()->SetOldPid(kFALSE);
899  // fAnalCuts->GetPidV0neg()->SetOldPid(kFALSE);
900  // */
901  // fAnalCuts->SetUsePID(kFALSE); // I don't want to use the PID through the cut object, but I will use the PID response directly!!!
902  // }
903 
904  // Setting properties of PID
905  fPIDCombined = new AliPIDCombined;
906  fPIDCombined->SetDefaultTPCPriors();
907  fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC+AliPIDResponse::kDetTOF);
908  //fPIDCombined->SetPriorDistribution((AliPID::EParticleType)ispec,fPriors[ispec]);
909 
910  fCounter = new AliNormalizationCounter("NormalizationCounter");
911  fCounter->Init();
912 
913  fCounterC = new AliNormalizationCounter("NormalizationCounterCorrMult");
914  fCounterC->SetStudyMultiplicity(kTRUE,1.);
915  fCounterC->Init();
916 
917  fListCounters = new TList();
918  fListCounters->SetOwner();
919  fListCounters->SetName("ListCounters");
920  fListCounters->Add(fCounter);
921  fListCounters->Add(fCounterC);
922 
923  PostData(2, fListCounters);
924 
925 
926  // Histograms from KF
927 
928  if (fCallKFVertexing){
929  fHistoDistanceLcToPrimVtx = new TH1D("fHistoDistanceLcToPrimVtx", "Lc distance to Prim Vertex from KF; distance [cm]", 1000, 0., 1);
930  fHistoDistanceV0ToPrimVtx = new TH1D("fHistoDistanceV0ToPrimVtx", "V0 distance to Prim Vertex from KF; distance [cm]", 1000, 0., 100.);
931  fHistoDistanceV0ToLc = new TH1D("fHistoDistanceV0ToLc", "V0 distance to Lc Vertex from KF; distance [cm]", 1000, 0., 100.);
932 
933  fHistoDistanceLcToPrimVtxSgn = new TH1D("fHistoDistanceLcToPrimVtxSgn", "Lc Sgn distance to Prim Vertex from KF; distance [cm]", 1000, 0., 1);
934  fHistoDistanceV0ToPrimVtxSgn = new TH1D("fHistoDistanceV0ToPrimVtxSgn", "V0 Sgn distance to Prim Vertex from KF; distance [cm]", 1000, 0., 100.);
935  fHistoDistanceV0ToLcSgn = new TH1D("fHistoDistanceV0ToLcSgn", "V0 Sgn distance to Lc Vertex from KF; distance [cm]", 1000, 0., 100.);
936 
937  fHistoVtxLcResidualToPrimVtx = new TH1D("fHistoVtxLcResidualToPrimVtx", "Residual between MC and KF (MC - KF): Lc to Prim Vtx; distance [cm]", 1000, -5., 5.);
938  fHistoVtxV0ResidualToPrimVtx = new TH1D("fHistoVtxV0ResidualToPrimVtx", "Residual between MC and KF (MC - KF): V0 to Prim Vtx; distance [cm]", 1000, -5., 5.);
939  fHistoVtxV0ResidualToLc = new TH1D("fHistoVtxV0ResidualToLc", "Residual between MC and KF: V0 to Lc (MC - KF); distance [cm]", 1000, -5., 5.);
940 
941  fHistoMassV0All = new TH1D("fHistoMassV0All", "V0 Mass; mass", 500, 0.4, 0.6);
942  fHistoDecayLengthV0All = new TH1D("fHistoDecayLengthV0All", "V0 Decay Length; decayLength", 500, -10, 10.0);
943  fHistoLifeTimeV0All = new TH1D("fHistoLifeTimeV0All", "V0 Life Time; lifeTime", 500, -10.0, 10.0);
944 
945  fHistoMassV0True = new TH1D("fHistoMassV0True", "True V0 Mass; mass", 500, 0.4, 0.6);
946  fHistoDecayLengthV0True = new TH1D("fHistoDecayLengthV0True", "True V0 Decay Length; decayLength", 500, -10, 10.0);
947  fHistoLifeTimeV0True = new TH1D("fHistoLifeTimeV0True", "True V0 Life Time; lifeTime", 500, -10.0, 10.0);
948 
949  fHistoMassV0TrueFromAOD = new TH1D("fHistoMassV0TrueFormAOD", "True V0 Mass (AOD); mass", 500, 0.4, 0.6);
950 
951  fHistoMassV0TrueK0S = new TH1D("fHistoMassV0TrueK0S", "True V0-K0S Mass; mass", 500, 0.4, 0.6);
952  fHistoDecayLengthV0TrueK0S = new TH1D("fHistoDecayLengthV0TrueK0S", "True V0-K0S Decay Length; decayLength", 500, -10, 10.0);
953  fHistoLifeTimeV0TrueK0S = new TH1D("fHistoLifeTimeV0TrueK0S", "True V0-K0S Life Time; lifeTime", 500, -10.0, 10.0);
954 
955  fHistoMassV0TrueK0SFromAOD = new TH1D("fHistoMassV0TrueK0SFormAOD", "True V0-K0S Mass (AOD); mass", 500, 0.4, 0.6);
956 
957  fHistoMassLcAll = new TH1D("fHistoMassLcAll", "Lc Mass; mass", 500, 2.0, 3.0);
958  fHistoDecayLengthLcAll = new TH1D("fHistoDecayLengthLcAll", "Lc Decay Length; decayLenght", 100000, -0.1, 0.1);
959  fHistoLifeTimeLcAll = new TH1D("fHistoLifeTimeLcAll", "Lc Life Time; lifeTime", 100000, -0.1, 0.1);
960 
961  fHistoMassLcTrue = new TH1D("fHistoMassLcTrue", "True Lc Mass; mass", 500, 2.0, 3.0);
962  fHistoDecayLengthLcTrue = new TH1D("fHistoDecayLengthLcTrue", "True Lc Decay Length; decayLength", 100000, -0.1, 0.1);
963  fHistoLifeTimeLcTrue = new TH1D("fHistoLifeTimeLcTrue", "True Lc Life Time; lifeTime", 100000, -0.1, 0.1);
964 
965  fHistoMassLcTrueFromAOD = new TH1D("fHistoMassLcTrueFromAOD", "True Lc Mass (AOD); mass", 500, 2.0, 3.0);
966 
967  fHistoMassV0fromLcAll = new TH1D("fHistoMassV0fromLcAll", "V0 mass from Lc built in KF; mass", 500, 0.4, 0.6);
968  fHistoDecayLengthV0fromLcAll = new TH1D("fHistoDecayLengthV0fromLcAll", "V0 Decay Length from Lc built in KF; decayLength", 500, 0, 10.0);
969  fHistoLifeTimeV0fromLcAll = new TH1D("fHistoLifeTimeV0fromLcAll", "V0 Life Time from Lc built in KF; lifeTime", 500, 0.0, 3.0);
970 
971  fHistoMassV0fromLcTrue = new TH1D("fHistoMassV0fromLcTrue", "V0 mass from true Lc built in KF; mass", 500, 0.4, 0.6);
972  fHistoDecayLengthV0fromLcTrue= new TH1D("fHistoDecayLengthV0fromLcTrue", "V0 Decay Length from true Lc built in KF; decayLength", 500, 0, 10.0);
973  fHistoLifeTimeV0fromLcTrue = new TH1D("fHistoLifeTimeV0fromLcTrue", "V0 Life Time from true Lc built in KF; lifeTime", 500, 0.0, 3.0);
974 
975  fHistoMassLcSgn = new TH1D("fHistoMassLcSgn", "True Lc Signal Mass; mass", 500, 2.0, 3.0);
976  fHistoMassLcSgnFromAOD = new TH1D("fHistoMassLcSgnFromAOD", "True Lc Signal Mass (AOD); mass", 500, 2.0, 3.0);
977  fHistoDecayLengthLcSgn = new TH1D("fHistoDecayLengthLcSgn", "True Lc Signal Decay Length; decayLength", 100000, -0.1, 0.1);
978  fHistoLifeTimeLcSgn = new TH1D("fHistoLifeTimeLcSgn", "True Lc Signal Life Time; lifeTime", 100000, -0.1, 0.1);
979 
980  fHistoMassV0fromLcSgn = new TH1D("fHistoMassV0fromLcSgn", "V0 from True Lc Signal Mass; mass", 500, 0.4, 0.6);
981  fHistoDecayLengthV0fromLcSgn = new TH1D("fHistoDecayLengthV0fromLcSgn", "V0 True Lc Signal Decay Length; decayLength", 500, 0, 10.0);
982  fHistoLifeTimeV0fromLcSgn = new TH1D("fHistoLifeTimeV0fromLcSgn", "V0 True Lc Signal Life Time; lifeTime", 500, 0.0, 3.0);
983 
984  fHistoKF = new TH2D("fHistoKF", "Summary from KF; V0 KF; Lc KF", 16, -0.5, 15.5, 16, -0.5, 15.5);
985  fHistoKFV0 = new TH1D("fHistoKFV0", "Summary from KF; V0 KF", 16, -0.5, 15.5);
986  fHistoKFLc = new TH1D("fHistoKFLc", "Summary from KF; V0 KF", 16, -0.5, 15.5);
987  TString axisLabel[16] = {"AllOk", "M_NotOk", "Sm_NotOk", "Dl_NotOk",
988  "Lt_NotOk", "M_Sm_NotOk", "M_Dl_NotOk", "M_Lt_NotOk",
989  "Dl_Sm_NotOk", "Dl_Lt_NotOk", "Sm_Lt_NotOk", "M_Sm_Dl_NotOk",
990  "M_Sm_Lt_NotOk", "Sm_Dl_Lt_NotOk", "M_Dl_Lt_NotOk", "All_NotOk"};
991 
992  for (Int_t ibin = 1; ibin <= 16; ibin++){
993  fHistoKF->GetXaxis()->SetBinLabel(ibin, axisLabel[ibin-1].Data());
994  fHistoKF->GetYaxis()->SetBinLabel(ibin, axisLabel[ibin-1].Data());
995  fHistoKFV0->GetXaxis()->SetBinLabel(ibin, axisLabel[ibin-1].Data());
996  fHistoKFLc->GetXaxis()->SetBinLabel(ibin, axisLabel[ibin-1].Data());
997  }
998 
999  fHistoMassKFV0 = new TH2D("fHistoMassKFV0", "mass vs sigmaMass for V0; mass; sigmaMass", 500, 0.4, 0.6, 500, 0., 10);
1000  fHistoDecayLengthKFV0 = new TH2D("fHistoDecayLengthKFV0", "decayLength vs sigmaDecayLength for V0; decayLength; sigmaDecayLength", 500, -10, 10, 500, 0., 10);
1001  fHistoLifeTimeKFV0 = new TH2D("fHistoLifeTimeKFV0", "lifeTime vs sigmalifeTime for V0; lifeTime; sigmaLifeTime", 500, -10, 10, 500, 0., 10);
1002 
1003  fHistoMassKFLc = new TH2D("fHistoMassKFLc", "mass vs sigmaMass for Lc; mass; sigmaMass", 500, 0.4, 0.6, 500, 0., 10);
1004  fHistoDecayLengthKFLc = new TH2D("fHistoDecayLengthKFLc", "decayLength vs sigmaDecayLength for Lc; decayLength; sigmaDecayLength", 500, -10, 10, 500, 0., 10);
1005  fHistoLifeTimeKFLc = new TH2D("fHistoLifeTimeKFLc", "lifeTime vs sigmalifeTime for Lc; lifeTime; sigmaLifeTime", 500, -10, 10, 500, 0., 10);
1006 
1007  fHistoArmenterosPodolanskiV0KF = new TH2D("fHistoArmenterosPodolanskiV0KF", "V0 ArmenterosPodolanski from KF; #alpha; Qt", 1000, -1, 1, 1000, 0, 1);
1008  fHistoArmenterosPodolanskiV0KFSgn = new TH2D("fHistoArmenterosPodolanskiV0KFSgn", "V0 (signal) ArmenterosPodolanski from KF; #alpha; Qt", 1000, -1, 1, 1000, 0, 1);
1009  fHistoArmenterosPodolanskiV0AOD = new TH2D("fHistoArmenterosPodolanskiV0AOD", "V0 ArmenterosPodolanski from AOD; #alpha; Qt", 1000, -1, 1, 1000, 0, 1);
1010  fHistoArmenterosPodolanskiV0AODSgn = new TH2D("fHistoArmenterosPodolanskiV0AODSgn", "V0 (signal) ArmenterosPodolanski from AOD; #alpha; Qt", 1000, -1, 1, 1000, 0, 1);
1011  }
1012 
1013  fOutputKF = new TList();
1014  fOutputKF->SetOwner();
1015  fOutputKF->SetName("listHistoKF");
1016 
1017  if (fCallKFVertexing){
1027  fOutputKF->Add(fHistoMassV0All);
1038  fOutputKF->Add(fHistoMassLcAll);
1051  fOutputKF->Add(fHistoMassLcSgn);
1058  fOutputKF->Add(fHistoKF);
1059  fOutputKF->Add(fHistoKFV0);
1060  fOutputKF->Add(fHistoKFLc);
1061  fOutputKF->Add(fHistoMassKFV0);
1064  fOutputKF->Add(fHistoMassKFLc);
1071  }
1072 
1073  // weight function from ratio of flat value (1/30) to pythia
1074  // use to normalise to flat distribution (should lead to flat pT distribution)
1075  fFuncWeightPythia = new TF1("funcWeightPythia","1./(30. *[0]*x/TMath::Power(1.+(TMath::Power((x/[1]),[3])),[2]))", 0.15, 30);
1076  fFuncWeightPythia->SetParameters(0.36609, 1.94635, 1.40463,2.5);
1077 
1078  // weight function from the ratio of the LHC13d3 MC
1079  // and FONLL calculations for pp data
1080  fFuncWeightFONLL5overLHC13d3 = new TF1("funcWeightFONLL5overLHC13d3","([0]*x)/TMath::Power([2],(1+TMath::Power([3],x/[1])))+[4]*TMath::Exp([5]+[6]*x)+[7]*TMath::Exp([8]*x)", 0.15, 30.);
1081  fFuncWeightFONLL5overLHC13d3->SetParameters(2.94999e+00, 3.47032e+00, 2.81278e+00, 2.5, 1.93370e-02, 3.86865e+00, -1.54113e-01, 8.86944e-02, 2.56267e-02);
1082 
1083  fFuncWeightFONLL5overLHC13d3Lc = new TF1("funcWeightFONLL5overLHC13d3Lc","([0]*x)/TMath::Power([2],(1+TMath::Power([3],x/[1])))+[4]*TMath::Exp([5]+[6]*x)+[7]*TMath::Exp([8]*x)", 0.15, 20.);
1084  fFuncWeightFONLL5overLHC13d3Lc->SetParameters(5.94428e+01, 1.63585e+01, 9.65555e+00, 6.71944e+00, 8.88338e-02, 2.40477e+00, -4.88649e-02, -6.78599e-01, -2.10951e-01);
1085 
1086  PostData(6, fOutputKF);
1087 
1088  PostData(7, fListWeight);
1089 
1090  if (!fFillTree) {
1091  Printf("Booking methods");
1092  // creating the BDT and TMVA reader
1093  fVarsTMVA = new Float_t[fNVars];
1095  fReader = new TMVA::Reader( "!Color:!Silent" );
1096  std::vector<std::string> inputNamesVec;
1097  TObjArray *tokens = fNamesTMVAVar.Tokenize(",");
1098  for(Int_t i = 0; i < tokens->GetEntries(); i++){
1099  TString variable = ((TObjString*)(tokens->At(i)))->String();
1100  std::string tmpvar = variable.Data();
1101  inputNamesVec.push_back(tmpvar);
1102  if (fUseXmlWeightsFile || fUseXmlFileFromCVMFS) fReader->AddVariable(variable.Data(), &fVarsTMVA[i]);
1103  }
1104  delete tokens;
1105  TObjArray *tokensSpectators = fNamesTMVAVarSpectators.Tokenize(",");
1106  for(Int_t i = 0; i < tokensSpectators->GetEntries(); i++){
1107  TString variable = ((TObjString*)(tokensSpectators->At(i)))->String();
1108  if (fUseXmlWeightsFile || fUseXmlFileFromCVMFS) fReader->AddSpectator(variable.Data(), &fVarsTMVASpectators[i]);
1109  }
1110  delete tokensSpectators;
1111  if (fUseWeightsLibrary) {
1112  void* lib = dlopen(fTMVAlibName.Data(), RTLD_NOW);
1113  void* p = dlsym(lib, Form("%s", fTMVAlibPtBin.Data()));
1114  IClassifierReader* (*maker1)(std::vector<std::string>&) = (IClassifierReader* (*)(std::vector<std::string>&)) p;
1115  fBDTReader = maker1(inputNamesVec);
1116  }
1117 
1118  if (fUseXmlWeightsFile) fReader->BookMVA("BDT method", fXmlWeightsFile);
1119 
1120  if (fUseXmlFileFromCVMFS){
1121  TString pathToFileCVMFS = AliDataFile::GetFileName(fXmlFileFromCVMFS.Data());
1122  if (pathToFileCVMFS.IsNull()){
1123  AliFatal("Cannot access data files from CVMFS");
1124  }
1125  fReader->BookMVA("BDT method", pathToFileCVMFS);
1126  }
1127  }
1128 
1129  return;
1130 }
1131 
1132 //_________________________________________________
1134 {
1136  if (!fInputEvent) {
1137  AliError("NO EVENT FOUND!");
1138  return;
1139  }
1140 
1141  fCurrentEvent++;
1142  AliDebug(2, Form("Processing event = %d", fCurrentEvent));
1143  AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
1144 
1145  if(fAODProtection >= 0){
1146  // Protection against different number of events in the AOD and deltaAOD
1147  // In case of discrepancy the event is rejected.
1148  Int_t matchingAODdeltaAODlevel = AliRDHFCuts::CheckMatchingAODdeltaAODevents();
1149  if (matchingAODdeltaAODlevel < 0 || (matchingAODdeltaAODlevel == 0 && fAODProtection == 1)) {
1150  // AOD/deltaAOD trees have different number of entries || TProcessID do not match while it was required
1151  fHistoEvents->Fill(0);
1152  return;
1153  }
1154  fHistoEvents->Fill(1);
1155  }
1156 
1157  TClonesArray *arrayLctopKos=0;
1158 
1159  TClonesArray *array3Prong = 0;
1160 
1161  if (!aodEvent && AODEvent() && IsStandardAOD()) {
1162  // In case there is an AOD handler writing a standard AOD, use the AOD
1163  // event in memory rather than the input (ESD) event.
1164  aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
1165  // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
1166  // have to taken from the AOD event hold by the AliAODExtension
1167  AliAODHandler* aodHandler = (AliAODHandler*)
1168  ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
1169 
1170  if (aodHandler->GetExtensions()) {
1171  AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
1172  AliAODEvent *aodFromExt = ext->GetAOD();
1173  arrayLctopKos = (TClonesArray*)aodFromExt->GetList()->FindObject("CascadesHF");
1174 
1175  array3Prong = (TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
1176  }
1177  } else {
1178  arrayLctopKos = (TClonesArray*)aodEvent->GetList()->FindObject("CascadesHF");
1179 
1180  array3Prong = (TClonesArray*)aodEvent->GetList()->FindObject("Charm3Prong");
1181  }
1182 
1183  if (!fUseMCInfo) {
1186  }
1187 
1188  Int_t runnumber = aodEvent->GetRunNumber();
1189  if (aodEvent->GetTriggerMask() == 0 && (runnumber >= 195344 && runnumber <= 195677)){
1190  AliDebug(3, "Event rejected because of null trigger mask");
1191  return;
1192  }
1193 
1194  fCounter->StoreEvent(aodEvent, fAnalCuts, fUseMCInfo);
1195 
1196  // mc analysis
1197  TClonesArray *mcArray = 0;
1198  AliAODMCHeader *mcHeader = 0;
1199 
1200  // multiplicity definition with tracklets
1201  fNTracklets_1 = static_cast<Int_t>(AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aodEvent, -1., 1.));
1202  fNTracklets_All = static_cast<Int_t>(AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aodEvent, -999., 999.));
1203  if (fUseMCInfo) {
1204  // MC array need for matching
1205  mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
1206  if (!mcArray) {
1207  AliError("Could not find Monte-Carlo in AOD");
1208  return;
1209  }
1210  // load MC header
1211  mcHeader = (AliAODMCHeader*)aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName());
1212  if (!mcHeader) {
1213  AliError("AliAnalysisTaskSELc2V0bachelorTMVAApp::UserExec: MC header branch not found!\n");
1214  return;
1215  }
1216 
1217  Double_t zMCVertex = mcHeader->GetVtxZ();
1218  if (TMath::Abs(zMCVertex) > fAnalCuts->GetMaxVtxZ()){
1219  AliDebug(3, Form("z coordinate of MC vertex = %f, it was required to be within [-%f, +%f], skipping event", zMCVertex, fAnalCuts->GetMaxVtxZ(), fAnalCuts->GetMaxVtxZ()));
1220  return;
1221  }
1222 
1223  FillMCHisto(mcArray);
1224 
1225  }
1226 
1227  //centrality
1229 
1230  // AOD primary vertex
1231  fVtx1 = (AliAODVertex*)aodEvent->GetPrimaryVertex();
1232  if (!fVtx1) {
1233  AliDebug(2, "No primary vertex found, returning");
1234  return;
1235  }
1236  if (fVtx1->GetNContributors()<1) {
1237  AliDebug(2, "Number of contributors to vertex < 1, returning");
1238  return;
1239  }
1240 
1241  Int_t fNTracklets_03=AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aodEvent,-0.3,0.3);
1242  Int_t fNTracklets_05=AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aodEvent,-0.5,0.5);
1243  Int_t fNTracklets_16=AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aodEvent,-1.6,1.6);
1244 
1245  Int_t vzeroMult=0, vzeroMultA=0, vzeroMultC=0;
1246  Int_t vzeroMultEq=0, vzeroMultAEq=0, vzeroMultCEq=0;
1247  AliAODVZERO *vzeroAOD = (AliAODVZERO*)aodEvent->GetVZEROData();
1248  if(vzeroAOD) {
1249  vzeroMultA = static_cast<Int_t>(vzeroAOD->GetMTotV0A());
1250  vzeroMultC = static_cast<Int_t>(vzeroAOD->GetMTotV0C());
1251  vzeroMult = vzeroMultA + vzeroMultC;
1252  vzeroMultAEq = static_cast<Int_t>(AliVertexingHFUtils::GetVZEROAEqualizedMultiplicity(aodEvent));
1253  vzeroMultCEq = static_cast<Int_t>(AliVertexingHFUtils::GetVZEROCEqualizedMultiplicity(aodEvent));
1254  vzeroMultEq = vzeroMultAEq + vzeroMultCEq;
1255  }
1256 
1257  Int_t countMult = fNTracklets_1;
1258  if(fMultiplicityEstimator==kNtrk03) { countMult = fNTracklets_03; }
1259  else if(fMultiplicityEstimator==kNtrk05) { countMult = fNTracklets_05; }
1260  else if(fMultiplicityEstimator==kNtrk10to16) { countMult = fNTracklets_16 - fNTracklets_1; }
1261  else if(fMultiplicityEstimator==kVZERO) { countMult = vzeroMult; }
1262  else if(fMultiplicityEstimator==kVZEROA) { countMult = vzeroMultA; }
1263  else if(fMultiplicityEstimator==kVZEROEq) { countMult = vzeroMultEq; }
1264  else if(fMultiplicityEstimator==kVZEROAEq) { countMult = vzeroMultAEq; }
1265 
1266  // Double_t countTreta1corr = fNTracklets_1;
1267  Double_t countCorr=countMult;
1268 
1269  if(fUseMultCorrection){
1270  // In case of VZERO multiplicity, consider the zvtx correction flag
1271  // fDoVZER0ParamVertexCorr: 0= none, 1= usual d2h, 2=AliESDUtils
1272  Bool_t isDataDrivenZvtxCorr=kTRUE;
1273  Int_t vzeroMultACorr=vzeroMultA, vzeroMultCCorr=vzeroMultC, vzeroMultCorr=vzeroMult;
1274  Int_t vzeroMultAEqCorr=vzeroMultAEq, vzeroMultCEqCorr=vzeroMultCEq, vzeroMultEqCorr=vzeroMultEq;
1275 
1278  {
1279  if(fDoVZER0ParamVertexCorr==0){
1280  // do not correct
1281  isDataDrivenZvtxCorr=kFALSE;
1282  }
1283  else if (fDoVZER0ParamVertexCorr==2){
1284  // use AliESDUtils correction
1285  Float_t zvtx = fVtx1->GetZ();
1286  isDataDrivenZvtxCorr=kFALSE;
1287  vzeroMultACorr = static_cast<Int_t>(AliESDUtils::GetCorrV0A(vzeroMultA,zvtx));
1288  vzeroMultCCorr = static_cast<Int_t>(AliESDUtils::GetCorrV0C(vzeroMultC,zvtx));
1289  vzeroMultCorr = vzeroMultACorr + vzeroMultCCorr;
1290  vzeroMultAEqCorr = static_cast<Int_t>(AliESDUtils::GetCorrV0A(vzeroMultAEq,zvtx));
1291  vzeroMultCEqCorr =static_cast<Int_t>( AliESDUtils::GetCorrV0C(vzeroMultCEq,zvtx));
1292  vzeroMultEqCorr = vzeroMultAEqCorr + vzeroMultCEqCorr;
1293  if(fMultiplicityEstimator==kVZERO) { countCorr = vzeroMultCorr; }
1294  else if(fMultiplicityEstimator==kVZEROA) { countCorr = vzeroMultACorr; }
1295  else if(fMultiplicityEstimator==kVZEROEq) { countCorr = vzeroMultEqCorr; }
1296  else if(fMultiplicityEstimator==kVZEROAEq) { countCorr = vzeroMultAEqCorr; }
1297  }
1298  }
1299  if(isDataDrivenZvtxCorr){
1300  TProfile* estimatorAvg = GetEstimatorHistogram(aodEvent);
1301  if(estimatorAvg){
1302  // countTreta1corr = static_cast<Int_t>(AliVertexingHFUtils::GetCorrectedNtracklets(estimatorAvg,fNTracklets_1,vtx1->GetZ(),fRefMult));
1303  countCorr = static_cast<Int_t>(AliVertexingHFUtils::GetCorrectedNtracklets(estimatorAvg,countMult,fVtx1->GetZ(),fRefMult));
1304  }
1305  }
1306  }
1307 
1308  fCounterC->StoreEvent(aodEvent, fAnalCuts, fUseMCInfo,countCorr);
1309 
1310  Bool_t isSelectedMultCut = kTRUE;
1311  if(fUseMultiplicityCut){
1312  // Multiplicity cut used
1313  if(countCorr >= fMultiplicityCutMin && countCorr < fMultiplicityCutMax){
1314  // Within multiplicity range
1315  fHistoEvents->Fill(5);
1316  }
1317  else{
1318  // Outside multiplicity range
1319  isSelectedMultCut = kFALSE;
1320  }
1321  }
1322 
1324 
1325  if ( !fIsEventSelected || !isSelectedMultCut) {
1326  fHistoEvents->Fill(2);
1327  return; // don't take into account not selected events
1328  }
1329 
1330  // check on the timestamp
1331  AliVHeader* h = aodEvent->GetHeader();
1332  UInt_t timestamp = h->GetTimeStamp();
1333  //Printf("timestamp = %d, cut = %u", timestamp, fTimestampCut);
1334  if (fTimestampCut != 0) {
1335  //Printf("timestamp = %d, cut = %u", timestamp, fTimestampCut);
1336  if (timestamp > fTimestampCut) {
1337  fHistoEvents->Fill(3);
1338  return;
1339  }
1340  }
1341 
1342  fHistoEvents->Fill(4);
1343 
1348 
1350 
1351  if(fUseMultCorrection){
1352  fHistoNtrUnCorr->Fill(countMult);
1353  fHistoNtrCorr->Fill(countCorr);
1354  fHistoVzVsNtrUnCorr->Fill(fVtx1->GetZ(),countMult);
1355  fHistoVzVsNtrCorr->Fill(fVtx1->GetZ(),countCorr);
1356  }
1357 
1358  //Setting magnetic field for KF vertexing
1359  fBField = aodEvent->GetMagneticField();
1360  AliKFParticle::SetField(fBField);
1361 
1362  Int_t nSelectedAnal = 0;
1363  if (fIsK0sAnalysis) {
1364  MakeAnalysisForLc2prK0S(aodEvent, arrayLctopKos, mcArray,
1365  nSelectedAnal, fAnalCuts,
1366  array3Prong, mcHeader);
1367  }
1368  fCounter->StoreCandidates(aodEvent, nSelectedAnal, kFALSE);
1369 
1370  Double_t nchWeight = 1.;
1371  if (fNTracklets_1 > 0) {
1372 
1373  if(!fHistoMCNch) AliInfo("Input histos to evaluate Nch weights missing");
1374  if(fHistoMCNch) nchWeight *= fHistoMCNch->GetBinContent(fHistoMCNch->FindBin(fNTracklets_1));
1375  }
1376 
1377 
1378  PostData(1, fOutput);
1379  PostData(2, fListCounters);
1380  PostData(4, fVariablesTreeSgn);
1381  PostData(5, fVariablesTreeBkg);
1382  PostData(6, fOutputKF);
1383  PostData(7, fListWeight);
1384 }
1385 //-------------------------------------------------------------------------------
1387 
1389 
1390  for (Int_t iPart = 0; iPart < mcArray->GetEntriesFast(); iPart++) {
1391  AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart));
1392  if (!mcPart){
1393  AliError("Failed casting particle from MC array!, Skipping particle");
1394  continue;
1395  }
1396  Int_t pdg = mcPart->GetPdgCode();
1397  if (TMath::Abs(pdg) != 4122){
1398  AliDebug(2, Form("MC particle %d is not a Lc: its pdg code is %d", iPart, pdg));
1399  continue;
1400  }
1401  AliDebug(2, Form("Step 0 ok: MC particle %d is a Lc: its pdg code is %d", iPart, pdg));
1402  Int_t labeldaugh0 = mcPart->GetDaughterLabel(0);
1403  Int_t labeldaugh1 = mcPart->GetDaughterLabel(1);
1404  if (labeldaugh0 <= 0 || labeldaugh1 <= 0){
1405  AliDebug(2, Form("The MC particle doesn't have correct daughters, skipping!!"));
1406  continue;
1407  }
1408  else if (labeldaugh1 - labeldaugh0 == 1){
1409  AliDebug(2, Form("Step 1 ok: The MC particle has correct daughters!!"));
1410  AliAODMCParticle* daugh0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labeldaugh0));
1411  AliAODMCParticle* daugh1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labeldaugh1));
1412  if(!daugh0 || !daugh1){
1413  AliDebug(2,"Particle daughters not properly retrieved!");
1414  return;
1415  }
1416  Int_t pdgCodeDaugh0 = TMath::Abs(daugh0->GetPdgCode());
1417  Int_t pdgCodeDaugh1 = TMath::Abs(daugh1->GetPdgCode());
1418  AliAODMCParticle* bachelorMC = daugh0;
1419  AliAODMCParticle* v0MC = daugh1;
1420  AliDebug(2, Form("pdgCodeDaugh0 = %d, pdgCodeDaugh1 = %d", pdgCodeDaugh0, pdgCodeDaugh1));
1421  if ((pdgCodeDaugh0 == 311 && pdgCodeDaugh1 == 2212) || (pdgCodeDaugh0 == 2212 && pdgCodeDaugh1 == 311)){
1422  // 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-
1424  if (pdgCodeDaugh0 == 311 && pdgCodeDaugh1 == 2212) {
1425  bachelorMC = daugh1;
1426  v0MC = daugh0;
1427  }
1428  AliDebug(2, Form("Number of Daughters of v0 = %d", v0MC->GetNDaughters()));
1429  if (v0MC->GetNDaughters() != 1) {
1430  AliDebug(2, "The K0 does not decay in 1 body only! Impossible... Continuing...");
1431  continue;
1432  }
1433  else { // So far: Lc --> K0 + p, K0 with 1 daughter
1434  AliDebug(2, "Step 2 ok: The K0 does decay in 1 body only! ");
1435  Int_t labelK0daugh = v0MC->GetDaughterLabel(0);
1436  AliAODMCParticle* partK0S = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelK0daugh));
1437  if(!partK0S){
1438  AliError("Error while casting particle! returning a NULL array");
1439  continue;
1440  }
1441  else { // So far: Lc --> K0 + p, K0 with 1 daughter that we can access
1442  if (partK0S->GetNDaughters() != 2 || TMath::Abs(partK0S->GetPdgCode() != 310)){
1443  AliDebug(2, "The K0 daughter is not a K0S or does not decay in 2 bodies");
1444  continue;
1445  }
1446  else { // So far: Lc --> K0 + p, K0 --> K0S, K0S in 2 bodies
1447  AliDebug(2, "Step 3 ok: The K0 daughter is a K0S and does decay in 2 bodies");
1448  Int_t labelK0Sdaugh0 = partK0S->GetDaughterLabel(0);
1449  Int_t labelK0Sdaugh1 = partK0S->GetDaughterLabel(1);
1450  AliAODMCParticle* daughK0S0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelK0Sdaugh0));
1451  AliAODMCParticle* daughK0S1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelK0Sdaugh1));
1452  if (!daughK0S0 || ! daughK0S1){
1453  AliDebug(2, "Could not access K0S daughters, continuing...");
1454  continue;
1455  }
1456  else { // So far: Lc --> K0 + p, K0 --> K0S, K0S in 2 bodies that we can access
1457  AliDebug(2, "Step 4 ok: Could access K0S daughters, continuing...");
1458  Int_t pdgK0Sdaugh0 = daughK0S0->GetPdgCode();
1459  Int_t pdgK0Sdaugh1 = daughK0S1->GetPdgCode();
1460  if (TMath::Abs(pdgK0Sdaugh0) != 211 || TMath::Abs(pdgK0Sdaugh1) != 211){
1461  AliDebug(2, "The K0S does not decay in pi+pi-, continuing");
1462  //AliInfo("The K0S does not decay in pi+pi-, continuing");
1463  }
1464  else { // Full chain: Lc --> K0 + p, K0 --> K0S, K0S --> pi+pi-
1465  if (fAnalCuts->IsInFiducialAcceptance(mcPart->Pt(), mcPart->Y())) {
1466  AliDebug(2, Form("----> Filling histo with pt = %f", mcPart->Pt()));
1467  if(TMath::Abs(mcPart->Y()) < 0.5) fHistoMCLcK0SpGenLimAcc->Fill(mcPart->Pt());
1468  //AliInfo(Form("\nparticle = %d, Filling MC Gen histo\n", iPart));
1469  fHistoMCLcK0SpGen->Fill(mcPart->Pt());
1470  if(!(TMath::Abs(bachelorMC->Eta()) > 0.9 || bachelorMC->Pt() < 0.1 ||
1471  TMath::Abs(daughK0S0->Eta()) > 0.9 || daughK0S0->Pt() < 0.1 ||
1472  TMath::Abs(daughK0S1->Eta()) > 0.9 || daughK0S1->Pt() < 0.1)) {
1473  fHistoMCLcK0SpGenAcc->Fill(mcPart->Pt());
1474  }
1475  }
1476  else {
1477  AliDebug(2, "not in fiducial acceptance! Skipping");
1478  continue;
1479  }
1480  }
1481  }
1482  }
1483  }
1484  }
1485  }
1486  }
1487  } // closing loop over mcArray
1488 
1489  return;
1490 
1491 }
1492 
1493 //-------------------------------------------------------------------------------
1495  TClonesArray *arrayLctopKos,
1496  TClonesArray *mcArray,
1497  Int_t &nSelectedAnal,
1498  AliRDHFCutsLctoV0 *cutsAnal, TClonesArray *array3Prong,
1499  AliAODMCHeader* aodheader){
1500 
1502 
1503  Int_t pdgCand = 4122;
1504  Int_t pdgDgLctoV0bachelor[2] = {2212, 310};
1505  Int_t pdgDgV0toDaughters[2] = {211, 211};
1506 
1507  // loop to search for candidates Lc->K0sp
1508  Int_t nCascades= arrayLctopKos->GetEntriesFast();
1509 
1510  // loop over cascades to search for candidates Lc->p+K0S
1511 
1512  Int_t mcLabel = -1;
1513 
1515  for (Int_t iLctopK0s = 0; iLctopK0s < nCascades; iLctopK0s++) {
1516 
1517  // Lc candidates and K0s from Lc
1518  AliAODRecoCascadeHF* lcK0spr = dynamic_cast<AliAODRecoCascadeHF*>(arrayLctopKos->At(iLctopK0s));
1519  if (!lcK0spr) {
1520  AliDebug(2, Form("Cascade %d doens't exist, skipping",iLctopK0s));
1521  continue;
1522  }
1523 
1524  if (!(lcK0spr->CheckCascadeFlags())) {
1525  AliDebug(2, Form("Cascade %d is not flagged as Lc candidate",iLctopK0s));
1526  continue;
1527  }
1528 
1529  // use Preselect to filter out the tracks according to the pT
1530  if (cutsAnal->GetUsePreselect()){
1531  TObjArray arrTracks(2);
1532  for(Int_t ipr = 0; ipr < 2; ipr++){
1533  AliAODTrack *tr;
1534  if (ipr == 0) tr = vHF->GetProng(aodEvent, lcK0spr, ipr);
1535  else tr = (AliAODTrack*)(aodEvent->GetV0(lcK0spr->GetProngID(1)));
1536  arrTracks.AddAt(tr, ipr);
1537  }
1538  Int_t preSelectLc = cutsAnal->PreSelect(arrTracks);
1539  if (preSelectLc == 0) continue;
1540  }
1541 
1542  // fill the cascade candidate
1543  if(!vHF->FillRecoCasc(aodEvent, lcK0spr, kFALSE)){ //Fill the data members of the candidate only if they are empty.
1544  continue;
1545  }
1546  //if (!(vHF->RecoSecondaryVertexForCascades(aodEvent, lcK0spr))) continue;
1547 
1548 
1549  if (!lcK0spr->GetSecondaryVtx()) {
1550  AliInfo("No secondary vertex");
1551  continue;
1552  }
1553 
1554  if (lcK0spr->GetNDaughters()!=2) {
1555  AliDebug(2, Form("Cascade %d has not 2 daughters (nDaughters=%d)", iLctopK0s, lcK0spr->GetNDaughters()));
1556  continue;
1557  }
1558 
1559  AliAODv0 * v0part = dynamic_cast<AliAODv0*>(lcK0spr->Getv0());
1560  AliAODTrack * bachPart = dynamic_cast<AliAODTrack*>(lcK0spr->GetBachelor());
1561  if (!v0part || !bachPart) {
1562  AliDebug(2, Form("Cascade %d has no V0 or no bachelor object", iLctopK0s));
1563  continue;
1564  }
1565 
1566  if (!v0part->GetSecondaryVtx()) {
1567  AliDebug(2, Form("No secondary vertex for V0 by cascade %d", iLctopK0s));
1568  continue;
1569  }
1570 
1571  if (v0part->GetNDaughters()!=2) {
1572  AliDebug(2, Form("current V0 has not 2 daughters (onTheFly=%d, nDaughters=%d)", v0part->GetOnFlyStatus(), v0part->GetNDaughters()));
1573  continue;
1574  }
1575 
1576  AliAODTrack * v0Pos = dynamic_cast<AliAODTrack*>(lcK0spr->Getv0PositiveTrack());
1577  AliAODTrack * v0Neg = dynamic_cast<AliAODTrack*>(lcK0spr->Getv0NegativeTrack());
1578  if (!v0Neg || !v0Pos) {
1579  AliDebug(2, Form("V0 by cascade %d has no V0positive of V0negative object", iLctopK0s));
1580  continue;
1581  }
1582 
1583  if (v0Pos->Charge() == v0Neg->Charge()) {
1584  AliDebug(2, Form("V0 by cascade %d has daughters with the same sign: IMPOSSIBLE!", iLctopK0s));
1585  continue;
1586  }
1587 
1588  //Filling a control histogram with no cuts
1589  if (fUseMCInfo) {
1590 
1591  Int_t pdgCode=-2;
1592 
1593  // find associated MC particle for Lc -> p+K0 and K0S->pi+pi
1594  fmcLabelLc = lcK0spr->MatchToMC(pdgCand, pdgDgLctoV0bachelor[1], pdgDgLctoV0bachelor, pdgDgV0toDaughters, mcArray, kTRUE);
1595  //if (fmcLabelLc>=0) {
1596  if (fmcLabelLc != -1) {
1597  AliDebug(2, Form("----> cascade number %d (total cascade number = %d) is a Lc!", iLctopK0s, nCascades));
1598 
1599  AliAODMCParticle *partLc = dynamic_cast<AliAODMCParticle*>(mcArray->At(fmcLabelLc));
1600  if(partLc){
1601  pdgCode = partLc->GetPdgCode();
1602  if (pdgCode<0) AliDebug(2,Form(" ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ MClabel=%d ~~~~~~~~~~ pdgCode=%d", fmcLabelLc, pdgCode));
1603  pdgCode = TMath::Abs(pdgCode);
1604  fHistoLcBeforeCuts->Fill(1);
1605  }
1606  }
1607  else {
1608  fHistoLcBeforeCuts->Fill(0);
1609  pdgCode = -1;
1610  }
1611 
1612  }
1613 
1614  Int_t isLc = 0;
1615 
1616  if (fUseMCInfo) {
1617 
1618  Int_t pdgCode = -2;
1619 
1620  // find associated MC particle for Lc -> p+K0 and K0S->pi+pi
1621  mcLabel = lcK0spr->MatchToMC(pdgCand, pdgDgLctoV0bachelor[1], pdgDgLctoV0bachelor, pdgDgV0toDaughters, mcArray, kTRUE);
1622  if (mcLabel >= 0) {
1623  AliDebug(2, Form(" ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~cascade number %d (total cascade number = %d)", iLctopK0s, nCascades));
1624 
1625  AliAODMCParticle *partLc = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcLabel));
1626  if(partLc){
1627  pdgCode = partLc->GetPdgCode();
1628  if (pdgCode < 0) AliDebug(2, Form(" ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ MClabel=%d ~~~~~~~~~~ pdgCode=%d", mcLabel, pdgCode));
1629  pdgCode = TMath::Abs(pdgCode);
1630  isLc = 1;
1631  fHistoLc->Fill(1);
1632  }
1633  }
1634  else {
1635  fHistoLc->Fill(0);
1636  pdgCode = -1;
1637  }
1638  }
1639  AliDebug(2, Form("Analysing candidate %d\n", iLctopK0s));
1640  AliDebug(2, Form(">>>>>>>>>> Candidate is background, fFillOnlySgn = %d --> SKIPPING", fFillOnlySgn));
1641  if (!isLc) {
1642  if (fFillOnlySgn) { // if it is background, and we want only signal, we do not fill the tree
1643  continue;
1644  }
1645  else { // checking if we want to fill the background
1646  if (fKeepingOnlyHIJINGBkg){
1647  // we have decided to fill the background only when the candidate has the daugthers that all come from HIJING underlying event!
1648  Bool_t isCandidateInjected = fUtils->HasCascadeCandidateAnyDaughInjected(lcK0spr, aodheader, mcArray);
1649  if (!isCandidateInjected){
1650  AliDebug(2, "The candidate is from HIJING (i.e. not injected), keeping it to fill background");
1651  fHistoBackground->Fill(1);
1652  }
1653  else {
1654  AliDebug(2, "The candidate is NOT from HIJING, we skip it when filling background");
1655  fHistoBackground->Fill(0);
1656  continue;
1657  }
1658  }
1659  else if (fKeepingOnlyPYTHIABkg){
1660  // we have decided to fill the background only when the candidate has the daugthers that all come from PYTHIA underlying event!
1661  AliAODTrack *bachelor = (AliAODTrack*)lcK0spr->GetBachelor();
1662  AliAODTrack *v0pos = (AliAODTrack*)lcK0spr->Getv0PositiveTrack();
1663  AliAODTrack *v0neg = (AliAODTrack*)lcK0spr->Getv0NegativeTrack();
1664  if (!bachelor || !v0pos || !v0neg) {
1665  AliDebug(2, "Cannot retrieve one of the tracks while checking origin, continuing");
1666  continue;
1667  }
1668  else {
1669  Int_t labelbachelor = TMath::Abs(bachelor->GetLabel());
1670  Int_t labelv0pos = TMath::Abs(v0pos->GetLabel());
1671  Int_t labelv0neg = TMath::Abs(v0neg->GetLabel());
1672  AliAODMCParticle* MCbachelor = (AliAODMCParticle*)mcArray->At(labelbachelor);
1673  AliAODMCParticle* MCv0pos = (AliAODMCParticle*)mcArray->At(labelv0pos);
1674  AliAODMCParticle* MCv0neg = (AliAODMCParticle*)mcArray->At(labelv0neg);
1675  if (!MCbachelor || !MCv0pos || !MCv0neg) {
1676  AliDebug(2, "Cannot retrieve MC particle for one of the tracks while checking origin, continuing");
1677  continue;
1678  }
1679  else {
1680  Int_t isBachelorFromPythia = fUtils->CheckOrigin(mcArray, MCbachelor, kTRUE);
1681  Int_t isv0posFromPythia = fUtils->CheckOrigin(mcArray, MCv0pos, kTRUE);
1682  Int_t isv0negFromPythia = fUtils->CheckOrigin(mcArray, MCv0neg, kTRUE);
1683  if (isBachelorFromPythia != 0 && isv0posFromPythia != 0 && isv0negFromPythia != 0){
1684  AliDebug(2, "The candidate is from PYTHIA (i.e. all daughters originate from a quark), keeping it to fill background");
1685  fHistoBackground->Fill(2);
1686  }
1687  else {
1688  AliDebug(2, "The candidate is NOT from PYTHIA, we skip it when filling background");
1689  fHistoBackground->Fill(3);
1690  continue;
1691  }
1692  }
1693  }
1694  }
1695  }
1696  }
1697 
1698  //FillLc2pK0Sspectrum(lcK0spr, isLc, nSelectedAnal, cutsAnal, mcArray, iLctopK0s);
1699  FillLc2pK0Sspectrum(lcK0spr, isLc, nSelectedAnal, cutsAnal, mcArray, mcLabel);
1700  }
1701 
1702  delete vHF;
1703 
1704  return;
1705 
1706 }
1707 //________________________________________________________________________
1709  Int_t isLc,
1710  Int_t &nSelectedAnal,
1711  AliRDHFCutsLctoV0 *cutsAnal,
1712  TClonesArray *mcArray, Int_t iLctopK0s){
1713  //
1715  //
1716  /*
1717  if (!part->GetOwnPrimaryVtx()) {
1718  //Printf("No primary vertex for Lc found!!");
1719  part->SetOwnPrimaryVtx(fVtx1);
1720  }
1721  else {
1722  //Printf("Yu-huuuu!!! primary vertex for Lc found!!");
1723  }
1724  */
1725  Double_t invmassLc = part->InvMassLctoK0sP();
1726 
1727  AliAODv0 * v0part = part->Getv0();
1728  Bool_t onFlyV0 = v0part->GetOnFlyStatus(); // on-the-flight V0s
1729  if (onFlyV0){ // on-the-fly V0
1730  if (isLc) { // Lc
1731  fHistoLcOnTheFly->Fill(2.);
1732  }
1733  else { // not Lc
1734  fHistoLcOnTheFly->Fill(0.);
1735  }
1736  }
1737  else { // offline V0
1738  if (isLc) { // Lc
1739  fHistoLcOnTheFly->Fill(3.);
1740  }
1741  else { // not Lc
1742  fHistoLcOnTheFly->Fill(1.);
1743  }
1744  }
1745 
1746  Double_t dcaV0 = v0part->GetDCA();
1747  Double_t invmassK0s = v0part->MassK0Short();
1748 
1749  if ( (cutsAnal->IsInFiducialAcceptance(part->Pt(), part->Y(4122))) ) {
1750  if (isLc) {
1751  fHistoFiducialAcceptance->Fill(3.);
1752  }
1753  else {
1754  fHistoFiducialAcceptance->Fill(1.);
1755  }
1756  }
1757  else {
1758  if (isLc) {
1759  fHistoFiducialAcceptance->Fill(2.);
1760  }
1761  else {
1762  fHistoFiducialAcceptance->Fill(0.);
1763  }
1764  }
1765 
1766  Int_t isInV0window = (((cutsAnal->IsSelectedSingleCut(part, AliRDHFCuts::kCandidate, 2)) & (AliRDHFCutsLctoV0::kLcToK0Spr)) == (AliRDHFCutsLctoV0::kLcToK0Spr)); // cut on V0 invMass
1767 
1768  if (isInV0window == 0) {
1769  AliDebug(2, "No: The candidate has NOT passed the V0 window cuts!");
1770  if (isLc) AliDebug(2, "SIGNAL candidate rejected: V0 window cuts");
1771  return;
1772  }
1773  else AliDebug(2, "Yes: The candidate has passed the mass cuts!");
1774 
1775  Bool_t isInCascadeWindow = (((cutsAnal->IsSelectedSingleCut(part, AliRDHFCuts::kCandidate, 0)) & (AliRDHFCutsLctoV0::kLcToK0Spr)) == (AliRDHFCutsLctoV0::kLcToK0Spr)); // cut on Lc->p+K0S invMass
1776 
1777  if (!isInCascadeWindow) {
1778  AliDebug(2, "No: The candidate has NOT passed the cascade window cuts!");
1779  if (isLc) AliDebug(2, "SIGNAL candidate rejected: cascade window cuts");
1780  return;
1781  }
1782  else AliDebug(2, "Yes: The candidate has passed the cascade window cuts!");
1783 
1784  Bool_t isCandidateSelectedCuts = (((cutsAnal->IsSelected(part, AliRDHFCuts::kCandidate)) & (AliRDHFCutsLctoV0::kLcToK0Spr)) == (AliRDHFCutsLctoV0::kLcToK0Spr)); // kinematic/topological cuts
1785  AliDebug(2, Form("recoAnalysisCuts = %d", cutsAnal->IsSelected(part, AliRDHFCuts::kCandidate) & (AliRDHFCutsLctoV0::kLcToK0Spr)));
1786  if (!isCandidateSelectedCuts){
1787  AliDebug(2, "No: Analysis cuts kCandidate level NOT passed");
1788  if (isLc) AliDebug(2, "SIGNAL candidate rejected");
1789  return;
1790  }
1791  else {
1792  AliDebug(2, "Yes: Analysis cuts kCandidate level passed");
1793  }
1794 
1795  AliAODTrack *bachelor = (AliAODTrack*)part->GetBachelor();
1796  if (!bachelor) {
1797  AliDebug(2, Form("Very weird, the bachelor is not there... returning for this candidate"));
1798  return;
1799  }
1800 
1801  //Bool_t isBachelorID = (((cutsAnal->IsSelected(part,AliRDHFCuts::kPID))&(AliRDHFCutsLctoV0::kLcToK0Spr))==(AliRDHFCutsLctoV0::kLcToK0Spr)); // ID x bachelor
1802  Double_t probTPCTOF[AliPID::kSPECIES] = {-1.};
1803 
1804  UInt_t detUsed = fPIDCombined->ComputeProbabilities(bachelor, fPIDResponse, probTPCTOF);
1805  AliDebug(2, Form("detUsed (TPCTOF case) = %d", detUsed));
1806  Double_t probProton = -1.;
1807  // Double_t probPion = -1.;
1808  // Double_t probKaon = -1.;
1809  if (detUsed == (UInt_t)fPIDCombined->GetDetectorMask() ) {
1810  AliDebug(2, Form("We have found the detector mask for TOF + TPC: probProton will be set to %f", probTPCTOF[AliPID::kProton]));
1811  probProton = probTPCTOF[AliPID::kProton];
1812  // probPion = probTPCTOF[AliPID::kPion];
1813  // probKaon = probTPCTOF[AliPID::kKaon];
1814  }
1815  else { // if you don't have both TOF and TPC, try only TPC
1816  fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
1817  AliDebug(2, "We did not find the detector mask for TOF + TPC, let's see only TPC");
1818  detUsed = fPIDCombined->ComputeProbabilities(bachelor, fPIDResponse, probTPCTOF);
1819  AliDebug(2, Form(" detUsed (TPC case) = %d", detUsed));
1820  if (detUsed == (UInt_t)fPIDCombined->GetDetectorMask()) {
1821  probProton = probTPCTOF[AliPID::kProton];
1822  // probPion = probTPCTOF[AliPID::kPion];
1823  // probKaon = probTPCTOF[AliPID::kKaon];
1824  AliDebug(2, Form("TPC only worked: probProton will be set to %f", probTPCTOF[AliPID::kProton]));
1825  }
1826  else {
1827  AliDebug(2, "Only TPC did not work...");
1828  }
1829  // resetting mask to ask for both TPC+TOF
1830  fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC+AliPIDResponse::kDetTOF);
1831  }
1832  AliDebug(2, Form("probProton = %f", probProton));
1833 
1834  // now we get the TPC and TOF single PID probabilities (only for Proton, or the tree will explode :) )
1835  Double_t probProtonTPC = -1.;
1836  Double_t probProtonTOF = -1.;
1837  Double_t pidTPC[AliPID::kSPECIES]={-1.};
1838  Double_t pidTOF[AliPID::kSPECIES]={-1.};
1839  Int_t respTPC = fPIDResponse->ComputePIDProbability(AliPIDResponse::kDetTPC, bachelor, AliPID::kSPECIES, pidTPC);
1840  Int_t respTOF = fPIDResponse->ComputePIDProbability(AliPIDResponse::kDetTOF, bachelor, AliPID::kSPECIES, pidTOF);
1841  if (respTPC == AliPIDResponse::kDetPidOk) probProtonTPC = pidTPC[AliPID::kProton];
1842  if (respTOF == AliPIDResponse::kDetPidOk) probProtonTOF = pidTOF[AliPID::kProton];
1843 
1844  // checking V0 status (on-the-fly vs offline)
1845  if ( !( !onFlyV0 || (onFlyV0 && fUseOnTheFlyV0) ) ) {
1846  AliDebug(2, "On-the-fly discarded");
1847  return;
1848  }
1849 
1851  nSelectedAnal++;
1852  }
1853 
1854  if ( !(cutsAnal->IsInFiducialAcceptance(part->Pt(), part->Y(4122))) ) return;
1855 
1856  if ( !( ( (cutsAnal->IsSelected(part, AliRDHFCuts::kTracks)) & (AliRDHFCutsLctoV0::kLcToK0Spr)) == (AliRDHFCutsLctoV0::kLcToK0Spr) ) ) { // esd track cuts
1857  if (isLc) AliDebug(2, "SIGNAL candidate rejected");
1858  AliDebug(2, "No: Analysis cuts kTracks level NOT passed");
1859  return;
1860  }
1861  else {
1862  AliDebug(2, "Yes: Analysis cuts kTracks level passed");
1863  }
1864 
1865  Int_t pdgCand = 4122;
1866  Int_t pdgDgLctoV0bachelor[2] = {211, 3122}; // case of wrong decay! Lc --> L + pi
1867  Int_t pdgDgV0toDaughters[2] = {2212, 211}; // case of wrong decay! Lc --> L + pi
1868  Int_t isLc2LBarpi = 0, isLc2Lpi = 0;
1869  Int_t currentLabel = part->GetLabel();
1870  Int_t mcLabel = 0;
1871  if (fUseMCInfo) {
1872  mcLabel = part->MatchToMC(pdgCand, pdgDgLctoV0bachelor[1], pdgDgLctoV0bachelor, pdgDgV0toDaughters, mcArray, kTRUE);
1873  if (mcLabel >= 0) {
1874  if (bachelor->Charge() == -1) isLc2LBarpi=1;
1875  if (bachelor->Charge() == +1) isLc2Lpi=1;
1876  }
1877  }
1878 
1879  Int_t pdgDg2prong[2] = {211, 211};
1880  Int_t labelK0S = 0;
1881  Int_t isK0S = 0;
1882  if (fUseMCInfo) {
1883  labelK0S = v0part->MatchToMC(310, mcArray, 2, pdgDg2prong);
1884  if (labelK0S>=0) isK0S = 1;
1885  }
1886 
1887  pdgDg2prong[0] = 211;
1888  pdgDg2prong[1] = 2212;
1889  Int_t isLambda = 0;
1890  Int_t isLambdaBar = 0;
1891  Int_t lambdaLabel = 0;
1892  if (fUseMCInfo) {
1893  lambdaLabel = v0part->MatchToMC(3122, mcArray, 2, pdgDg2prong);
1894  if (lambdaLabel >= 0) {
1895  AliAODMCParticle *lambdaTrack = (AliAODMCParticle*)mcArray->At(lambdaLabel);
1896  if (lambdaTrack->GetPdgCode() == 3122) isLambda = 1;
1897  else if (lambdaTrack->GetPdgCode() == -3122) isLambdaBar = 1;
1898  }
1899  }
1900 
1901  pdgDg2prong[0] = 11;
1902  pdgDg2prong[1] = 11;
1903  Int_t isGamma = 0;
1904  Int_t gammaLabel = 0;
1905  if (fUseMCInfo) {
1906  gammaLabel = v0part->MatchToMC(22, mcArray, 2, pdgDg2prong);
1907  if (gammaLabel>=0) {
1908  AliAODMCParticle *gammaTrack = (AliAODMCParticle*)mcArray->At(gammaLabel);
1909  if (gammaTrack->GetPdgCode() == 22) isGamma = 1;
1910  }
1911  }
1912 
1913  Int_t pdgTemp = -1;
1914  if (currentLabel != -1){
1915  AliAODMCParticle *tempPart = (AliAODMCParticle*)mcArray->At(currentLabel);
1916  pdgTemp = tempPart->GetPdgCode();
1917  }
1918  if (isLc) AliDebug(2, Form("Signal: Candidate is a Lc in K0s+p"));
1919  else if (isLc2LBarpi) AliDebug(2, Form("Background: Candidate is a Lc in Lbar + pi"));
1920  else if (isLc2Lpi) AliDebug(2, Form("Background: Candidate is a Lc in L + pi"));
1921  else AliDebug(2, Form("Pure bkg: Candidate has pdg = %d", pdgTemp));
1922  if (isK0S) AliDebug(2, Form("V0 is a K0S"));
1923  else if (isLambda) AliDebug(2, Form("V0 is a Lambda"));
1924  else if (isLambdaBar) AliDebug(2, Form("V0 is a LambdaBar"));
1925  else if (isGamma) AliDebug(2, Form("V0 is a Gamma"));
1926  //else AliDebug(2, Form("V0 is something else!!"));
1927 
1928  Double_t invmassLc2Lpi = part->InvMassLctoLambdaPi();
1929  Double_t invmassLambda = v0part->MassLambda();
1930  Double_t invmassLambdaBar = v0part->MassAntiLambda();
1931 
1932  //Double_t nSigmaITSpr = -999.;
1933  Double_t nSigmaTPCpr = -999.;
1934  Double_t nSigmaTOFpr = -999.;
1935 
1936  //Double_t nSigmaITSpi = -999.;
1937  Double_t nSigmaTPCpi = -999.;
1938  Double_t nSigmaTOFpi = -999.;
1939 
1940  //Double_t nSigmaITSka = -999.;
1941  Double_t nSigmaTPCka = -999.;
1942  Double_t nSigmaTOFka = -999.;
1943 
1944  /*
1945  cutsAnal->GetPidHF()->GetnSigmaITS(bachelor, 4, nSigmaITSpr);
1946  cutsAnal->GetPidHF()->GetnSigmaTPC(bachelor, 4, nSigmaTPCpr);
1947  cutsAnal->GetPidHF()->GetnSigmaTOF(bachelor, 4, nSigmaTOFpr);
1948  cutsAnal->GetPidHF()->GetnSigmaITS(bachelor, 2, nSigmaITSpi);
1949  cutsAnal->GetPidHF()->GetnSigmaTPC(bachelor, 2, nSigmaTPCpi);
1950  cutsAnal->GetPidHF()->GetnSigmaTOF(bachelor, 2, nSigmaTOFpi);
1951  cutsAnal->GetPidHF()->GetnSigmaITS(bachelor, 3, nSigmaITSka);
1952  cutsAnal->GetPidHF()->GetnSigmaTPC(bachelor, 3, nSigmaTPCka);
1953  cutsAnal->GetPidHF()->GetnSigmaTOF(bachelor, 3, nSigmaTOFka);
1954  */
1955 
1957  nSigmaTPCpi = fPIDResponse->NumberOfSigmasTPC(bachelor, (AliPID::kPion));
1958  nSigmaTPCka = fPIDResponse->NumberOfSigmasTPC(bachelor, (AliPID::kKaon));
1959  nSigmaTPCpr = fPIDResponse->NumberOfSigmasTPC(bachelor, (AliPID::kProton));
1960  nSigmaTOFpi = fPIDResponse->NumberOfSigmasTOF(bachelor, (AliPID::kPion));
1961  nSigmaTOFka = fPIDResponse->NumberOfSigmasTOF(bachelor, (AliPID::kKaon));
1962  nSigmaTOFpr = fPIDResponse->NumberOfSigmasTOF(bachelor, (AliPID::kProton));
1963  }
1964  else {
1965  cutsAnal->GetPidHF()->GetnSigmaTPC(bachelor, (AliPID::kPion), nSigmaTPCpi);
1966  cutsAnal->GetPidHF()->GetnSigmaTPC(bachelor, (AliPID::kKaon), nSigmaTPCka);
1967  cutsAnal->GetPidHF()->GetnSigmaTPC(bachelor, (AliPID::kProton), nSigmaTPCpr);
1968  cutsAnal->GetPidHF()->GetnSigmaTOF(bachelor, (AliPID::kPion), nSigmaTOFpi);
1969  cutsAnal->GetPidHF()->GetnSigmaTOF(bachelor, (AliPID::kKaon), nSigmaTOFka);
1970  cutsAnal->GetPidHF()->GetnSigmaTOF(bachelor, (AliPID::kProton), nSigmaTOFpr);
1971  }
1972 
1973  Double_t ptLcMC = -1;
1974  Double_t weightPythia = -1, weight5LHC13d3 = -1, weight5LHC13d3Lc = -1;
1975 
1976  AliAODMCParticle *partLcMC = 0x0;
1977 
1978  if (fUseMCInfo) {
1979  if (iLctopK0s >= 0) {
1980  partLcMC = (AliAODMCParticle*)mcArray->At(iLctopK0s);
1981  ptLcMC = partLcMC->Pt();
1982  //Printf("--------------------- Reco pt = %f, MC particle pt = %f", part->Pt(), ptLcMC);
1983  weightPythia = fFuncWeightPythia->Eval(ptLcMC);
1984  weight5LHC13d3 = fFuncWeightFONLL5overLHC13d3->Eval(ptLcMC);
1985  weight5LHC13d3Lc = fFuncWeightFONLL5overLHC13d3Lc->Eval(ptLcMC);
1986  }
1987  }
1988 
1989  Double_t weightNch = 1;
1990  if (fUseMCInfo) {
1991  //Int_t nChargedMCPhysicalPrimary=AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(mcArray,-1.0,1.0);
1992  // if(nChargedMCPhysicalPrimary > 0)
1993 
1994  if(fNTracklets_1 > 0){
1995  if(!fHistoMCNch) AliDebug(2, "Input histos to evaluate Nch weights missing");
1996  if(fHistoMCNch) weightNch *= fHistoMCNch->GetBinContent(fHistoMCNch->FindBin(fNTracklets_1));
1997  }
1998  }
1999 
2000 
2001  // Fill candidate variable Tree (track selection, V0 invMass selection)
2002  // if (!onFlyV0 && isInV0window && isInCascadeWindow && part->CosV0PointingAngle()>0.99 && TMath::Abs(nSigmaTPCpr) <= 3 &&
2003  // TMath::Abs(nSigmaTOFpr) <= 3 && v0part->Getd0Prong(0) < 20 && v0part->Getd0Prong(1) < 20) {
2004  if (!onFlyV0 && isInV0window && isInCascadeWindow){
2005 
2006  EBachelor bachCode = kBachInvalid;
2007  EK0S k0SCode = kK0SInvalid;
2008  if (fUseMCInfo) {
2009  bachCode = CheckBachelor(part, bachelor, mcArray);
2010  k0SCode = CheckK0S(part, v0part, mcArray);
2011  }
2012 
2013  Double_t V0KF[3] = {-999999, -999999, -999999};
2014  Double_t errV0KF[3] = {-999999, -999999, -999999};
2015  Double_t LcKF[3] = {-999999, -999999, -999999};
2016  Double_t errLcKF[3] = {-999999, -999999, -999999};
2017  Double_t distances[3] = {-999999, -999999, -999999};
2018  Double_t armPolKF[2] = {-999999, -999999};
2019 
2020  AliAODTrack *v0pos = (AliAODTrack*)part->Getv0PositiveTrack();
2021  AliAODTrack *v0neg = (AliAODTrack*)part->Getv0NegativeTrack();
2022  Int_t kfResult;
2023  TVector3 mom1(bachelor->Px(), bachelor->Py(), bachelor->Pz());
2024  TVector3 mom2(v0part->Px(), v0part->Py(), v0part->Pz());
2025  TVector3 momTot(part->Px(), part->Py(), part->Pz());
2026 
2027  Double_t Ql1 = mom1.Dot(momTot)/momTot.Mag();
2028  Double_t Ql2 = mom2.Dot(momTot)/momTot.Mag();
2029 
2030  Double_t alphaArmLc = (Ql1 - Ql2)/(Ql1 + Ql2);
2031  Double_t alphaArmLcCharge = ( bachelor->Charge() > 0 ? (Ql1 - Ql2)/(Ql1 + Ql2) : (Ql2 - Ql1)/(Ql1 + Ql2) );
2032  Double_t ptArmLc = mom1.Perp(momTot);
2033 
2034  Double_t massK0SPDG = TDatabasePDG::Instance()->GetParticle(310)->Mass(); // mass K0S PDG
2035  Double_t massPrPDG = TDatabasePDG::Instance()->GetParticle(2212)->Mass(); // mass Proton PDG
2036  Double_t massLcPDG = TDatabasePDG::Instance()->GetParticle(4122)->Mass(); // mass Lc PDG
2037 
2038  // Cosine of proton emission angle (theta*) in the rest frame of the mother particle
2039  // for prong ip (0 or 1) with mass hypotheses massLcPDG for mother particle (from AliAODRecoDecay)
2040  Double_t pStar = TMath::Sqrt((massLcPDG*massLcPDG-massPrPDG*massPrPDG-massK0SPDG*massK0SPDG)*(massLcPDG*massLcPDG-massPrPDG*massPrPDG-massK0SPDG*massK0SPDG)-4.*massPrPDG*massPrPDG*massK0SPDG*massK0SPDG)/(2.*massLcPDG);
2041  Double_t e = part->E(4122);
2042  Double_t beta = part->P()/e;
2043  Double_t gamma = e/massLcPDG;
2044  //Double_t cts = (Ql1/gamma-beta*TMath::Sqrt(pStar*pStar+massPrPDG*massPrPDG))/pStar;
2045 
2046  // Cosine of proton emission angle (theta*) in the rest frame of the mother particle
2047  // (from AliRDHFCutsLctoV0)
2048  TLorentzVector vpr, vk0s,vlc;
2049  vpr.SetXYZM(part->PxProng(0), part->PyProng(0), part->PzProng(0), massPrPDG);
2050  vk0s.SetXYZM(part->PxProng(1), part->PyProng(1), part->PzProng(1), massK0SPDG);
2051  vlc = vpr + vk0s;
2052  TVector3 vboost = vlc.BoostVector();
2053  vpr.Boost(-vboost);
2054  Double_t cts = TMath::Cos(vpr.Angle(vlc.Vect()));
2055 
2056  if (fCallKFVertexing){
2057  kfResult = CallKFVertexing(part, v0part, bachelor, mcArray, &V0KF[0], &errV0KF[0], &LcKF[0], &errLcKF[0], &distances[0], &armPolKF[0]);
2058  AliDebug(2, Form("Result from KF = %d", kfResult));
2059  }
2060 
2061  /*
2062  for (Int_t i = 0; i< 3; i++){
2063  Printf("i = %d, V0KF = %f, errV0KF = %f, LcKF = %f, errLcKF = %f", V0KF[i], errV0KF[i], LcKF[i], errLcKF[i]);
2064  }
2065  */
2066 
2067  Double_t countTreta1corr = fNTracklets_1;
2068 
2069  AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
2070 
2071  if(fUseMultCorrection){
2072  TProfile* estimatorAvg = GetEstimatorHistogram(aodEvent);
2073  if(estimatorAvg) {
2074  countTreta1corr = static_cast<Int_t>(AliVertexingHFUtils::GetCorrectedNtracklets(estimatorAvg, fNTracklets_1, fVtx1->GetZ(), fRefMult));
2075  }
2076  }
2077 
2078  AliAODVertex *primvert = dynamic_cast<AliAODVertex*>(part->GetPrimaryVtx());
2079 
2080  Double_t d0z0bach[2], covd0z0bach[3];
2081  bachelor->PropagateToDCA(primvert, fBField, kVeryBig, d0z0bach, covd0z0bach);
2082  Double_t tx[3];
2083  bachelor->GetXYZ(tx);
2084  tx[0] -= primvert->GetX();
2085  tx[1] -= primvert->GetY();
2086  tx[2] -= primvert->GetZ();
2087  Double_t innerpro = tx[0]*part->Px()+tx[1]*part->Py();
2088  Double_t signd0 = 1.;
2089  if(innerpro<0.) signd0 = -1.;
2090 
2091  signd0 = signd0*TMath::Abs(d0z0bach[0]);
2092 
2093  if (fUseMCInfo) { // save full tree if on MC
2094  fCandidateVariables[0] = invmassLc;
2095  fCandidateVariables[1] = invmassLc2Lpi;
2096  fCandidateVariables[2] = invmassK0s;
2097  fCandidateVariables[3] = invmassLambda;
2098  fCandidateVariables[4] = invmassLambdaBar;
2100  fCandidateVariables[6] = dcaV0;
2101  fCandidateVariables[7] = part->Getd0Prong(0);
2102  fCandidateVariables[8] = part->Getd0Prong(1);
2103  fCandidateVariables[9] = nSigmaTPCpr;
2104  fCandidateVariables[10] = nSigmaTOFpr;
2105  fCandidateVariables[11] = bachelor->Pt();
2106  fCandidateVariables[12] = v0pos->Pt();
2107  fCandidateVariables[13] = v0neg->Pt();
2108  fCandidateVariables[14] = v0part->Getd0Prong(0);
2109  fCandidateVariables[15] = v0part->Getd0Prong(1);
2110  fCandidateVariables[16] = v0part->Pt();
2111  fCandidateVariables[17] = v0part->InvMass2Prongs(0,1,11,11);
2112  fCandidateVariables[18] = part->Pt();
2113  fCandidateVariables[19] = probProton;
2114  fCandidateVariables[20] = part->Eta();
2115  fCandidateVariables[21] = v0pos->Eta();
2116  fCandidateVariables[22] = v0neg->Eta();
2117  fCandidateVariables[23] = probProtonTPC;
2118  fCandidateVariables[24] = probProtonTOF;
2119  fCandidateVariables[25] = bachelor->Eta();
2120  fCandidateVariables[26] = part->P();
2121  fCandidateVariables[27] = bachelor->P();
2122  fCandidateVariables[28] = v0part->P();
2123  fCandidateVariables[29] = v0pos->P();
2124  fCandidateVariables[30] = v0neg->P();
2125  fCandidateVariables[31] = v0part->Eta();
2126  fCandidateVariables[32] = ptLcMC;
2127  fCandidateVariables[33] = part->DecayLengthV0();
2128  fCandidateVariables[34] = bachCode;
2129  fCandidateVariables[35] = k0SCode;
2130  fCandidateVariables[36] = v0part->AlphaV0();
2131  fCandidateVariables[37] = v0part->PtArmV0();
2132 
2133  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)));
2134  fCandidateVariables[38] = cts;
2135  fCandidateVariables[39] = weightPythia;
2136  fCandidateVariables[40] = weight5LHC13d3;
2137  fCandidateVariables[41] = weight5LHC13d3Lc;
2138  fCandidateVariables[42] = weightNch;
2140  fCandidateVariables[44] = countTreta1corr;
2141  fCandidateVariables[45] = signd0;
2144  if (partLcMC)
2145  fCandidateVariables[48] = fUtils->CheckOrigin(mcArray, partLcMC, kTRUE);
2146  else
2147  fCandidateVariables[48] = -1;
2148  fCandidateVariables[49] = nSigmaTPCpi;
2149  fCandidateVariables[50] = nSigmaTPCka;
2150  fCandidateVariables[51] = bachelor->GetTPCmomentum();
2151  }
2152  else { //remove MC-only variables from tree if data
2153  fCandidateVariables[0] = invmassLc;
2154  fCandidateVariables[1] = v0part->AlphaV0();
2155  fCandidateVariables[2] = invmassK0s;
2156  fCandidateVariables[3] = invmassLambda;
2157  fCandidateVariables[4] = invmassLambdaBar;
2159  fCandidateVariables[6] = dcaV0;
2160  fCandidateVariables[7] = part->Getd0Prong(0);
2161  fCandidateVariables[8] = part->Getd0Prong(1);
2162  fCandidateVariables[9] = nSigmaTPCpr;
2163  fCandidateVariables[10] = nSigmaTOFpr;
2164  fCandidateVariables[11] = bachelor->Pt();
2165  fCandidateVariables[12] = v0pos->Pt();
2166  fCandidateVariables[13] = v0neg->Pt();
2167  fCandidateVariables[14] = v0part->Getd0Prong(0);
2168  fCandidateVariables[15] = v0part->Getd0Prong(1);
2169  fCandidateVariables[16] = v0part->Pt();
2170  fCandidateVariables[17] = bachelor->GetTPCmomentum();
2171  fCandidateVariables[18] = part->Pt();
2172  fCandidateVariables[19] = probProton;
2173  fCandidateVariables[20] = v0pos->Eta();
2174  fCandidateVariables[21] = bachelor->P();
2175  fCandidateVariables[22] = bachelor->Eta();
2176  fCandidateVariables[23] = v0part->P();
2177  fCandidateVariables[24] = part->DecayLengthV0();
2178  fCandidateVariables[25] = nSigmaTPCpi;
2179  fCandidateVariables[26] = nSigmaTPCka;
2181  fCandidateVariables[28] = countTreta1corr;
2182  fCandidateVariables[29] = cts;
2183  fCandidateVariables[30] = signd0;
2186  fCandidateVariables[33] = -1;
2187  fCandidateVariables[34] = v0part->PtArmV0();
2188  }
2189 
2190  // fill multiplicity histograms for events with a candidate
2191  //fHistNtrUnCorrEvWithCand->Fill(fNTracklets_1, weightNch);
2192  //fHistNtrCorrEvWithCand->Fill(countTreta1corr, weightNch);
2193  if (fUseMCInfo) {
2194  if (isLc){
2195  AliDebug(2, Form("Reco particle %d --> Filling Sgn", iLctopK0s));
2196  if(fFillTree) fVariablesTreeSgn->Fill();
2197  fHistoCodesSgn->Fill(bachCode, k0SCode);
2198  }
2199  else {
2200  if (fFillOnlySgn == kFALSE){
2201  AliDebug(2, "Filling Bkg");
2202  fVariablesTreeBkg->Fill();
2203  if(fFillTree) fHistoCodesBkg->Fill(bachCode, k0SCode);
2204  }
2205  }
2206  }
2207  else {
2208  if(fFillTree) fVariablesTreeSgn->Fill();
2209  }
2210 
2211 
2212  if(!fFillTree){
2213  std::vector<Double_t> inputVars(fNVars);
2214  if (fNVars == 14) {
2215  inputVars[0] = invmassK0s;
2216  inputVars[1] = part->Getd0Prong(0);
2217  inputVars[2] = part->Getd0Prong(1);
2218  inputVars[3] = bachelor->Pt();
2219  inputVars[4] = (part->DecayLengthV0())*0.497/(v0part->P());
2220  inputVars[5] = part->CosV0PointingAngle();
2221  inputVars[6] = cts;
2222  inputVars[7] = signd0;
2223  inputVars[8] = bachelor->P();
2224  inputVars[9] = nSigmaTOFpr;
2225  inputVars[10] = nSigmaTPCpr;
2226  inputVars[11] = nSigmaTPCpi;
2227  inputVars[12] = nSigmaTPCka;
2228  inputVars[13] = bachelor->GetTPCmomentum();
2229  }
2230  else if (fNVars == 11) {
2231  inputVars[0] = invmassK0s;
2232  inputVars[1] = part->Getd0Prong(0);
2233  inputVars[2] = part->Getd0Prong(1);
2234  inputVars[3] = (part->DecayLengthV0())*0.497/(v0part->P());
2235  inputVars[4] = part->CosV0PointingAngle();
2236  inputVars[5] = cts;
2237  inputVars[6] = signd0;
2238  inputVars[7] = nSigmaTOFpr;
2239  inputVars[8] = nSigmaTPCpr;
2240  inputVars[9] = nSigmaTPCpi;
2241  inputVars[10] = nSigmaTPCka;
2242  }
2243  else if (fNVars == 10) {
2244  inputVars[0] = invmassK0s;
2245  inputVars[1] = part->Getd0Prong(0);
2246  inputVars[2] = part->Getd0Prong(1);
2247  inputVars[3] = (part->DecayLengthV0())*0.497/(v0part->P());
2248  inputVars[4] = part->CosV0PointingAngle();
2249  inputVars[5] = signd0;
2250  inputVars[6] = nSigmaTOFpr;
2251  inputVars[7] = nSigmaTPCpr;
2252  inputVars[8] = nSigmaTPCpi;
2253  inputVars[9] = nSigmaTPCka;
2254  }
2255  else if (fNVars == 7) {
2256  inputVars[0] = invmassK0s;
2257  inputVars[1] = part->Getd0Prong(0);
2258  inputVars[2] = part->Getd0Prong(1);
2259  inputVars[3] = (part->DecayLengthV0())*0.497/(v0part->P());
2260  inputVars[4] = part->CosV0PointingAngle();
2261  inputVars[5] = cts;
2262  inputVars[6] = signd0;
2263  }
2264 
2265  for (Int_t i = 0; i < fNVars; i++) {
2266  fVarsTMVA[i] = inputVars[i];
2267  }
2268 
2269  Double_t BDTResponse = -1;
2270  Double_t tmva = -1;
2271  if (fUseXmlWeightsFile || fUseXmlFileFromCVMFS) tmva = fReader->EvaluateMVA("BDT method");
2272  if (fUseWeightsLibrary) BDTResponse = fBDTReader->GetMvaValue(inputVars);
2273  //Printf("BDTResponse = %f, invmassLc = %f", BDTResponse, invmassLc);
2274  //Printf("tmva = %f", tmva);
2275  fBDTHisto->Fill(BDTResponse, invmassLc);
2276  fBDTHistoTMVA->Fill(tmva, invmassLc);
2277  if (fDebugHistograms) {
2278  if (fUseXmlWeightsFile || fUseXmlFileFromCVMFS) BDTResponse = tmva; // we fill the debug histogram with the output from the xml file
2279  fBDTHistoVsMassK0S->Fill(BDTResponse, invmassK0s);
2280  fBDTHistoVstImpParBach->Fill(BDTResponse, part->Getd0Prong(0));
2281  fBDTHistoVstImpParV0->Fill(BDTResponse, part->Getd0Prong(1));
2282  fBDTHistoVsBachelorPt->Fill(BDTResponse, bachelor->Pt());
2283  fBDTHistoVsCombinedProtonProb->Fill(BDTResponse, probProton);
2284  fBDTHistoVsCtau->Fill(BDTResponse, (part->DecayLengthV0())*0.497/(v0part->P()));
2285  fBDTHistoVsCosPAK0S->Fill(BDTResponse, part->CosV0PointingAngle());
2286  fBDTHistoVsSignd0->Fill(BDTResponse, signd0);
2287  fBDTHistoVsCosThetaStar->Fill(BDTResponse, cts);
2288  fBDTHistoVsnSigmaTPCpr->Fill(BDTResponse, nSigmaTPCpr);
2289  fBDTHistoVsnSigmaTOFpr->Fill(BDTResponse, nSigmaTOFpr);
2290  fBDTHistoVsnSigmaTPCpi->Fill(BDTResponse, nSigmaTPCpi);
2291  fBDTHistoVsnSigmaTPCka->Fill(BDTResponse, nSigmaTPCka);
2292  fBDTHistoVsBachelorP->Fill(BDTResponse, bachelor->P());
2293  fBDTHistoVsBachelorTPCP->Fill(BDTResponse, bachelor->GetTPCmomentum());
2294  fHistoNsigmaTPC->Fill(bachelor->P(), nSigmaTPCpr);
2295  fHistoNsigmaTOF->Fill(bachelor->P(), nSigmaTOFpr);
2296  }
2297  }
2298 
2299  }
2300 
2301  return;
2302 
2303 
2304 }
2305 
2306 //________________________________________________________________________
2307 Int_t AliAnalysisTaskSELc2V0bachelorTMVAApp::CallKFVertexing(AliAODRecoCascadeHF *cascade, AliAODv0* v0part, AliAODTrack* bach, TClonesArray *mcArray,
2308  Double_t* V0KF, Double_t* errV0KF, Double_t* LcKF, Double_t* errLcKF,
2309  Double_t* distances, Double_t* armPolKF) {
2310 
2311  //
2314  //
2315 
2316  Int_t codeKFV0 = -1, codeKFLc = -1;
2317 
2318  AliKFVertex primVtxCopy;
2319  Int_t nt = 0, ntcheck = 0;
2320  Double_t pos[3] = {0., 0., 0.};
2321 
2322  fVtx1->GetXYZ(pos);
2323  Int_t contr = fVtx1->GetNContributors();
2324  Double_t covmatrix[6] = {0.};
2325  fVtx1->GetCovarianceMatrix(covmatrix);
2326  Double_t chi2 = fVtx1->GetChi2();
2327  AliESDVertex primaryESDVtxCopy(pos, covmatrix, chi2, contr, "Vertex");
2328 
2329  // topological constraint
2330  primVtxCopy = AliKFVertex(primaryESDVtxCopy);
2331  nt = primaryESDVtxCopy.GetNContributors();
2332  ntcheck = nt;
2333 
2334  Int_t pdg[2] = {211, -211};
2335  Int_t pdgLc[2] = {2212, 310};
2336 
2337  Int_t pdgDgV0toDaughters[2] = {211, 211};
2338 
2339  Int_t mcLabelV0 = v0part->MatchToMC(310, mcArray, 2, pdgDgV0toDaughters);
2340 
2341  // the KF vertex for the V0 has to be built with the prongs of the V0!
2342  Bool_t isMCokV0 = kTRUE, isBkgV0 = kFALSE;
2343  AliKFParticle V0, positiveV0KF, negativeV0KF;
2344  Int_t labelsv0daugh[2] = {-1, -1};
2345  Int_t idv0daugh[2] = {-1, -1};
2346  AliExternalTrackParam* esdv0Daugh1 = 0x0;
2347  AliExternalTrackParam* esdv0Daugh2 = 0x0;
2348  for(Int_t ipr = 0; ipr < 2; ipr++){ // 0 is positive, 1 is negative
2349  AliAODTrack *aodTrack = (AliAODTrack*)v0part->GetDaughter(ipr);
2350  if(!aodTrack) {
2351  AliDebug(2, "No V0 daughters available");
2352  return -1;
2353  }
2354  Double_t xyz[3], pxpypz[3], cv[21];
2355  Short_t sign;
2356  aodTrack->GetXYZ(xyz);
2357  aodTrack->PxPyPz(pxpypz);
2358  aodTrack->GetCovarianceXYZPxPyPz(cv);
2359  sign = aodTrack->Charge();
2360  AliExternalTrackParam tmp1( xyz, pxpypz, cv, sign);
2361 
2362  if (ipr == 0) esdv0Daugh1 = new AliExternalTrackParam( xyz, pxpypz, cv, sign);
2363  else esdv0Daugh2 = new AliExternalTrackParam( xyz, pxpypz, cv, sign);
2364  labelsv0daugh[ipr] = TMath::Abs(aodTrack->GetLabel());
2365  idv0daugh[ipr] = aodTrack->GetID();
2366  if (labelsv0daugh[ipr] == -1) isBkgV0 = kTRUE;
2367 
2368  //Printf("v0 daughter %d has label %d", ipr, labelsv0daugh[ipr]);
2369 
2370  AliKFParticle daughterKF(*aodTrack, pdg[ipr]); // we assume that the PDG is correct
2371  if (aodTrack->Charge() > 0) { // assigning positive and negative track to KF V0 for Armenteros-Podolanski plot
2372  positiveV0KF = daughterKF;
2373  }
2374  else {
2375  negativeV0KF = daughterKF;
2376  }
2377  }
2378 
2379  Double_t xn = 0., xp = 0.;//, dca;
2380  AliDebug(2, Form("bField = %f, esdv0Daugh1 = %p, esdv0Daugh2 = %p", fBField, esdv0Daugh1, esdv0Daugh2));
2381  // dca = esdv0Daugh1->GetDCA(esdv0Daugh2, fBField, xn, xp);
2382 
2383  AliExternalTrackParam tr1(*esdv0Daugh1);
2384  AliExternalTrackParam tr2(*esdv0Daugh2);
2385  tr1.PropagateTo(xn, fBField);
2386  tr2.PropagateTo(xp, fBField);
2387 
2388  AliKFParticle daughterKF1(tr1, 211);
2389  AliKFParticle daughterKF2(tr2, 211);
2390  V0.AddDaughter(positiveV0KF);
2391  V0.AddDaughter(negativeV0KF);
2392  //V0.AddDaughter(daughterKF1);
2393  //V0.AddDaughter(daughterKF2);
2394 
2395  delete esdv0Daugh1;
2396  delete esdv0Daugh2;
2397  esdv0Daugh1=0;
2398  esdv0Daugh2=0;
2399  // Checking the quality of the KF V0 vertex
2400  if( V0.GetNDF() < 1 ) {
2401  //Printf("Number of degrees of freedom < 1, continuing");
2402  return -1;
2403  }
2404  if( TMath::Sqrt(TMath::Abs(V0.GetChi2()/V0.GetNDF())) > fCutKFChi2NDF ) {
2405  //Printf("Chi2 per DOF too big, continuing");
2406  return -1;
2407  }
2408 
2409  if(ftopoConstraint && nt > 0){
2410  for(Int_t ipr = 0; ipr < 2; ipr++){ // 0 is positive, 1 is negative
2411  AliAODTrack *aodTrack = (AliAODTrack*)v0part->GetDaughter(ipr);
2412  //* subtruct daughters from primary vertex
2413  if(!aodTrack->GetUsedForPrimVtxFit()) {
2414  //Printf("Track %d was not used for primary vertex, continuing", i);
2415  continue;
2416  }
2417  AliKFParticle daughterKF(*aodTrack, pdg[ipr]); // we assume that the PDG is correct
2418  primVtxCopy -= daughterKF;
2419  ntcheck--;
2420  }
2421  }
2422 
2423  // Check V0 Chi^2 deviation from primary vertex // not needed for V0 for Lc decay!!
2424  /*
2425  if( V0.GetDeviationFromVertex( primVtxCopy ) < fCutKFDeviationFromVtxV0) {
2426  //Printf("Deviation from vertex too big, continuing");
2427  return -1;
2428  }
2429  */
2430 
2431  //* Get V0 invariant mass
2432  Double_t massV0 = 999999, sigmaMassV0 = 999999;
2433  Int_t retMV0 = V0.GetMass( massV0, sigmaMassV0 );
2434  if( retMV0 ) {
2435  if (massV0 < 0) {
2436  codeKFV0 = 1; // Mass not ok
2437  if (sigmaMassV0 > 1e19) codeKFV0 = 5; // Mass and SigmaMass not ok
2438  }
2439  else if (sigmaMassV0 > 1e19) codeKFV0 = 2; // SigmaMass not ok
2440  }
2441  fHistoMassKFV0->Fill(massV0, sigmaMassV0);
2442 
2443  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]);
2444  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]);
2445 
2446  Printf("Vertices: KF: x = %f, y = %f, z = %f", V0.GetX(), V0.GetY(), V0.GetZ());
2447  Printf("Vertices: AOD: x = %f, y = %f, z = %f", v0part->Xv(), v0part->Yv(), v0part->Zv());
2448 
2449  //Printf("Got MC vtx for V0");
2450  if (fUseMCInfo && TMath::Abs(labelsv0daugh[0] - labelsv0daugh[1]) == 1) {
2451  AliAODMCParticle* tmpdaughv01 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelsv0daugh[0]));
2452  AliAODMCParticle* tmpdaughv02 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelsv0daugh[1]));
2453  if (!tmpdaughv01 && labelsv0daugh[0] > 0){
2454  AliDebug(2, "Could not access MC info for first daughter of V0, continuing");
2455  }
2456  if (!tmpdaughv02 && labelsv0daugh[1] > 0){
2457  AliDebug(2, "Could not access MC info for second daughter of V0, continuing");
2458  }
2459  if(tmpdaughv01){
2460  Double_t xPionMC = tmpdaughv01->Xv(); //Production vertex of Pion --> Where K0S decays
2461  Double_t yPionMC = tmpdaughv01->Yv();
2462  Double_t zPionMC = tmpdaughv01->Zv();
2463  //Printf("Got MC vtx for Pion");
2464  Printf("Vertices: MC: x = %f, y = %f, z = %f", xPionMC, yPionMC, zPionMC);
2465  }
2466  }
2467  else {
2468  Printf("Not a true V0");
2469  }
2470  //massV0=-1;//return -1;// !!!!
2471 
2472  // now use what we just try with the bachelor, to build the Lc
2473 
2474  // topological constraint
2475  nt = primVtxCopy.GetNContributors();
2476  ntcheck = nt;
2477 
2478  Bool_t isMCokLc = kTRUE, isBkgLc = kFALSE;
2479  AliKFParticle Lc;
2480  Int_t labelsLcdaugh[2] = {-1, -1};
2481  labelsLcdaugh[0] = TMath::Abs(bach->GetLabel());
2482  labelsLcdaugh[1] = mcLabelV0;
2483 
2484  if (bach->Charge() < 0) pdgLc[0] = -pdgLc[0];
2485  AliKFParticle daughterKFLc(*bach, pdgLc[0]);
2486  Lc.AddDaughter(daughterKFLc);
2487  TParticlePDG* particlePDG = TDatabasePDG::Instance()->GetParticle(310);
2488  Double_t massPDGK0S = particlePDG->Mass();
2489  V0.SetMassConstraint(massPDGK0S);
2490  Lc.AddDaughter(V0);
2491  if( Lc.GetNDF() < 1 ) {
2492  AliDebug(2, Form("Lc: Number of degrees of freedom < 1 (%d), continuing", Lc.GetNDF()));
2493  return -1;
2494  }
2495  if( TMath::Sqrt(TMath::Abs(Lc.GetChi2()/Lc.GetNDF())) > fCutKFChi2NDF) {
2496  AliDebug(2, Form("Lc: Chi2 per DOF too big, continuing (%f)", TMath::Sqrt(TMath::Abs(Lc.GetChi2()/Lc.GetNDF()))));
2497  return -1;
2498  }
2499 
2500  if(ftopoConstraint && nt > 0){
2501  //* subtruct daughters from primary vertex
2502  if(!bach->GetUsedForPrimVtxFit()) {
2503  AliDebug(3, "Lc: Bachelor was not used for primary vertex, not subtracting it from primary vertex");
2504  }
2505  else{
2506  primVtxCopy -= daughterKFLc;
2507  ntcheck--;
2508  }
2509  /* the V0 was added above, so it is ok to remove it without checking
2510  if(!V0->GetUsedForPrimVtxFit()) {
2511  Printf("Lc: V0 was not used for primary vertex, continuing");
2512  continue;
2513  }
2514  */
2515  //primVtxCopy -= V0;
2516  //ntcheck--;
2517  }
2518 
2519  // Check Lc Chi^2 deviation from primary vertex
2520  /*
2521  if( Lc.GetDeviationFromVertex( primVtxCopy ) > fCutKFDeviationFromVtx) {
2522  AliDebug(2, Form("Lc: Deviation from vertex too big, continuing (%f)", Lc.GetDeviationFromVertex( primVtxCopy )));
2523  return -1;
2524  }
2525  if(ftopoConstraint){
2526  if(ntcheck>0) {
2527  // Add Lc to primary vertex to improve the primary vertex resolution
2528  primVtxCopy += Lc;
2529  Lc.SetProductionVertex(primVtxCopy);
2530  }
2531  }
2532  */
2533  //* Check chi^2
2534  if( TMath::Sqrt( TMath::Abs(Lc.GetChi2()/Lc.GetNDF())) > fCutKFChi2NDF) {
2535  AliDebug(2, Form("Lc: Final Chi2 per DOF too big, continuing (%f)", TMath::Sqrt( TMath::Abs(Lc.GetChi2()/Lc.GetNDF()))));
2536  return -1;
2537  }
2538 
2539  if(ftopoConstraint){
2540  V0.SetProductionVertex(Lc);
2541  }
2542 
2543  // After setting the vertex of the V0, getting/filling some info
2544 
2545  //* Get V0 decayLength
2546  Double_t decayLengthV0 = 999999, sigmaDecayLengthV0 = 999999;
2547  Int_t retDLV0 = V0.GetDecayLength( decayLengthV0, sigmaDecayLengthV0 );
2548  if( retDLV0 ) {
2549  if (sigmaDecayLengthV0 > 1e19) {
2550  codeKFV0 = 3; // DecayLength not ok
2551  if (massV0 < 0) {
2552  codeKFV0 = 6; // DecayLength and Mass not ok
2553  if (sigmaMassV0 > 1e19) codeKFV0 = 11; // DecayLength and Mass and SigmaMass not ok
2554  }
2555  else if (sigmaMassV0 > 1e19) codeKFV0 = 8; // DecayLength and SigmaMass not ok
2556  }
2557  }
2558  fHistoDecayLengthKFV0->Fill(decayLengthV0, sigmaDecayLengthV0);
2559 
2560  //* Get V0 life time
2561  Double_t lifeTimeV0 = 999999, sigmaLifeTimeV0 = 999999;
2562  Int_t retTLV0 = V0.GetLifeTime( lifeTimeV0, sigmaLifeTimeV0 );
2563  if( retTLV0 ) {
2564  if (sigmaLifeTimeV0 > 1e19) {
2565  codeKFV0 = 4; // LifeTime not ok
2566  if (sigmaDecayLengthV0 > 1e19) {
2567  codeKFV0 = 9; // LifeTime and DecayLength not ok
2568  if (massV0 < 0) {
2569  codeKFV0 = 14; // LifeTime and DecayLength and Mass not ok
2570  if (sigmaMassV0 > 1e19) codeKFV0 = 15; // LifeTime and DecayLength and Mass and SigmaMass not ok
2571  }
2572  else if (sigmaMassV0 > 1e19) codeKFV0 = 13; // LifeTime and DecayLength and SigmaMass not ok
2573  }
2574  else if (massV0 < 0) { // LifeTime and Mass and SigmaMass not ok
2575  codeKFV0 = 7; // LifeTime and Mass not ok
2576  if (sigmaMassV0 > 1e19) codeKFV0 = 12; // LifeTime and Mass and SigmaMass not ok
2577  }
2578  else if (sigmaMassV0 > 1e19) codeKFV0 = 10; // LifeTime and SigmaMass not ok
2579  }
2580  }
2581  fHistoLifeTimeKFV0->Fill(lifeTimeV0, sigmaLifeTimeV0);
2582 
2583  if (codeKFV0 == -1) codeKFV0 = 0;
2584  fHistoKFV0->Fill(codeKFV0);
2585 
2586  AliDebug(2, Form("V0: mass = %f, decay length = %f, life time = %f", massV0, decayLengthV0, lifeTimeV0 ));
2587 
2588  fHistoMassV0All->Fill(massV0);
2589  fHistoDecayLengthV0All->Fill(decayLengthV0);
2590  fHistoLifeTimeV0All->Fill(lifeTimeV0);
2591 
2592  Double_t qtAlphaV0[2] = {0., 0.};
2593  Double_t vtxV0KF[3] = {V0.GetX(), V0.GetY(), V0.GetZ()};
2594  positiveV0KF.TransportToPoint(vtxV0KF);
2595  negativeV0KF.TransportToPoint(vtxV0KF);
2596  V0.GetArmenterosPodolanski(positiveV0KF, negativeV0KF, qtAlphaV0);
2597  AliDebug(2, Form("Armenteros-Podolanski variables: alpha = %f, qt = %f", qtAlphaV0[1], qtAlphaV0[0]));
2598  fHistoArmenterosPodolanskiV0KF->Fill(qtAlphaV0[1], qtAlphaV0[0]);
2599  fHistoArmenterosPodolanskiV0AOD->Fill(v0part->AlphaV0(), v0part->PtArmV0());
2600  armPolKF[0] = qtAlphaV0[1];
2601  armPolKF[1] = qtAlphaV0[0];
2602 
2603  // Checking MC info for V0
2604 
2605  AliAODMCParticle *motherV0 = 0x0;
2606  AliAODMCParticle *daughv01 = 0x0;
2607  AliAODMCParticle *daughv02 = 0x0;
2608 
2609  if (fUseMCInfo) {
2610  daughv01 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelsv0daugh[0]));
2611  daughv02 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelsv0daugh[1]));
2612  if (!daughv01 && labelsv0daugh[0] > 0){
2613  AliDebug(2, "Could not access MC info for first daughter of V0, continuing");
2614  isMCokV0 = kFALSE;
2615  }
2616  if (!daughv02 && labelsv0daugh[1] > 0){
2617  AliDebug(2, "Could not access MC info for second daughter of V0, continuing");
2618  isMCokV0 = kFALSE;
2619  }
2620  if (isMCokV0){
2621  if( daughv01->GetMother() == daughv02->GetMother() && daughv01->GetMother()>=0 ){
2622  AliDebug(3, Form("The mother has label %d", daughv01->GetMother()));
2623  motherV0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughv01->GetMother()));
2624  if( motherV0 && TMath::Abs(motherV0->GetPdgCode()) != 21 ){ // These are all V0 that are truly V0, not only K0S, but no gluons
2625  if( motherV0->GetNDaughters() == 2 ){
2626  fHistoMassV0True->Fill(massV0);
2627  fHistoDecayLengthV0True->Fill(decayLengthV0);
2628  fHistoLifeTimeV0True->Fill(lifeTimeV0);
2629  fHistoMassV0TrueFromAOD->Fill(v0part->MassK0Short());
2630  if (TMath::Abs(motherV0->GetPdgCode()) == 310){ // These are true V0 that are also K0S
2631  fHistoMassV0TrueK0S->Fill(massV0);
2632  fHistoDecayLengthV0TrueK0S->Fill(decayLengthV0);
2633  fHistoLifeTimeV0TrueK0S->Fill(lifeTimeV0);
2634  fHistoMassV0TrueK0SFromAOD->Fill(v0part->MassK0Short());
2635  }
2636  }
2637  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()));
2638  }
2639  else if (!motherV0){
2640  AliDebug(3, "could not access MC info for mother, continuing");
2641  }
2642  else {
2643  AliDebug(3, "MC mother is a gluon, continuing");
2644  }
2645  }
2646  else {
2647  AliDebug(3, "Background V0!");
2648  isBkgV0 = kTRUE;
2649  }
2650  }
2651  }
2652 
2653  AliDebug(2, Form("isMCokV0 = %d, isBkgV0 = %d", (Int_t)isMCokV0, (Int_t)isBkgV0));
2654 
2655  // Going back to Lc
2656 
2657  //* Get Lc invariant mass
2658  Double_t massLc = 999999, sigmaMassLc= 999999;
2659  Int_t retMLc = Lc.GetMass( massLc, sigmaMassLc );
2660  if( retMLc ) {
2661  AliDebug(3, Form("----> Could not get mass (%e), and sigma(%e) for Lc, continuing", massLc, sigmaMassLc));
2662  if (massLc < 0) {
2663  codeKFLc = 1; // Mass not ok
2664  if (sigmaMassLc > 1e19) codeKFLc = 5; // Mass and SigmaMass not ok
2665  }
2666  else if (sigmaMassLc > 1e19) codeKFLc = 2; // SigmaMass not ok
2667  }
2668  fHistoMassKFLc->Fill(massLc, sigmaMassLc);
2669 
2670  //* Get Lc Decay length
2671  Double_t decayLengthLc = 999999, sigmaDecayLengthLc = 999999;
2672  Int_t retDLLc = Lc.GetDecayLength( decayLengthLc, sigmaDecayLengthLc );
2673  if( retDLLc ) {
2674  AliDebug(3, "----> Lc: Could not get decay length, and sigma");
2675  if (sigmaDecayLengthLc > 1e19) {
2676  codeKFLc = 3; // DecayLength not ok
2677  if (massLc < 0) {
2678  codeKFLc = 6; // DecayLength and Mass not ok
2679  if (sigmaMassLc > 1e19) codeKFLc = 11; // DecayLength and Mass and SigmaMass not ok
2680  }
2681  else if (sigmaMassLc > 1e19) codeKFLc = 8; // DecayLength and SigmaMass not ok
2682  }
2683  }
2684  AliDebug(3, Form("retDLLc = %d, with decayLength = %f and error = %e", retDLLc, decayLengthLc, sigmaDecayLengthLc));
2685 
2686  fHistoDecayLengthKFLc->Fill(decayLengthLc, sigmaDecayLengthLc);
2687 
2688  //* Get Lc life time
2689  Double_t lifeTimeLc = 999999, sigmaLifeTimeLc = 999999;
2690  Int_t retTLLc = Lc.GetLifeTime( lifeTimeLc, sigmaLifeTimeLc );
2691  if( retTLLc ) {
2692  AliDebug(3, "----> Lc: Could not get lifeTime, and sigma");
2693  if (sigmaLifeTimeLc > 1e19) {
2694  codeKFLc = 4; // LifeTime not ok
2695  if (sigmaDecayLengthLc > 1e19) {
2696  codeKFLc = 9; // LifeTime and DecayLength not ok
2697  if (massLc < 0) {
2698  codeKFLc = 14; // LifeTime and DecayLength and Mass not ok
2699  if (sigmaMassLc > 1e19) codeKFLc = 15; // LifeTime and DecayLength and Mass and SigmaMass not ok
2700  }
2701  else if (sigmaMassLc > 1e19) codeKFLc = 13; // LifeTime and DecayLength and SigmaMass not ok
2702  }
2703  else if (massLc < 0) { // LifeTime and Mass and SigmaMass not ok
2704  codeKFLc = 7; // LifeTime and Mass not ok
2705  if (sigmaMassLc > 1e19) codeKFLc = 12; // LifeTime and Mass and SigmaMass not ok
2706  }
2707  else if (sigmaMassLc > 1e19) codeKFLc = 10; // LifeTime and SigmaMass not ok
2708  }
2709  }
2710 
2711  fHistoLifeTimeKFLc->Fill(lifeTimeLc, sigmaLifeTimeLc);
2712 
2713  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));
2714 
2715  if (codeKFLc == -1) codeKFLc = 0;
2716  fHistoKFLc->Fill(codeKFLc);
2717 
2718  fHistoKF->Fill(codeKFV0, codeKFLc);
2719 
2720  // here we fill the histgrams for all the reconstructed KF vertices for the cascade
2721  fHistoMassLcAll->Fill(massLc);
2722  fHistoDecayLengthLcAll->Fill(decayLengthLc);
2723  fHistoLifeTimeLcAll->Fill(lifeTimeLc);
2724 
2725  fHistoMassV0fromLcAll->Fill(massV0);
2726  fHistoDecayLengthV0fromLcAll->Fill(decayLengthV0);
2727  fHistoLifeTimeV0fromLcAll->Fill(lifeTimeV0);
2728 
2729  Double_t xV0 = V0.GetX();
2730  Double_t yV0 = V0.GetY();
2731  Double_t zV0 = V0.GetZ();
2732 
2733  Double_t xLc = Lc.GetX();
2734  Double_t yLc = Lc.GetY();
2735  Double_t zLc = Lc.GetZ();
2736 
2737  Double_t xPrimVtx = primVtxCopy.GetX();
2738  Double_t yPrimVtx = primVtxCopy.GetY();
2739  Double_t zPrimVtx = primVtxCopy.GetZ();
2740 
2741  Double_t distanceLcToPrimVtx = TMath::Sqrt((xPrimVtx - xLc) * (xPrimVtx - xLc) +
2742  (yPrimVtx - yLc) * (yPrimVtx - yLc) +
2743  (zPrimVtx - zLc) * (zPrimVtx - zLc));
2744 
2745  Double_t distanceV0ToPrimVtx = TMath::Sqrt((xPrimVtx - xV0) * (xPrimVtx - xV0) +
2746  (yPrimVtx - yV0) * (yPrimVtx - yV0) +
2747  (zPrimVtx - zV0) * (zPrimVtx - zV0));
2748 
2749  Double_t distanceV0ToLc = TMath::Sqrt((xLc - xV0)*(xLc - xV0) +
2750  (yLc - yV0)*(yLc - yV0) +
2751  (zLc - zV0)*(zLc - zV0));
2752 
2753  //Printf("distanceLcToPrimVtx = %e, distanceV0ToPrimVtx= %f, distanceV0ToLc = %f", distanceLcToPrimVtx, distanceV0ToPrimVtx, distanceV0ToLc);
2754 
2755  fHistoDistanceLcToPrimVtx->Fill(distanceLcToPrimVtx);
2756  fHistoDistanceV0ToPrimVtx->Fill(distanceV0ToPrimVtx);
2757  fHistoDistanceV0ToLc->Fill(distanceV0ToLc);
2758 
2759  distances[0] = distanceLcToPrimVtx;
2760  distances[1] = distanceV0ToPrimVtx;
2761  distances[2] = distanceV0ToLc;
2762 
2763  if (fUseMCInfo) {
2764 
2765  AliAODMCParticle *daughv01Lc = 0x0;
2766  AliAODMCParticle *K0S = 0x0;
2767  AliAODMCParticle *daughv02Lc = 0x0;
2768 
2769  if (labelsLcdaugh[0] >= 0) {
2770  // Printf("Getting Bachelor from label %d", labelsLcdaugh[1]);
2771  daughv01Lc = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelsLcdaugh[0]));
2772  if (!daughv01Lc){
2773  AliDebug(3, "Could not access MC info for first daughter of Lc");
2774  isMCokLc = kFALSE;
2775  }
2776  else {
2777  AliDebug(2, Form("The bachelor has label = %d", daughv01Lc->GetLabel()));
2778  if (TMath::Abs(daughv01Lc->GetPdgCode()) != 2212) isBkgLc = kTRUE;
2779  }
2780  }
2781  else { // The bachelor is a fake
2782  isBkgLc = kTRUE;
2783  }
2784 
2785  if (labelsLcdaugh[1] >= 0) {
2786  //Printf("Getting K0S");
2787  K0S = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelsLcdaugh[1]));
2788  if (!K0S) {
2789  AliDebug(3, "Could not access MC info for second daughter of Lc");
2790  isMCokLc = kFALSE;
2791  }
2792  else{
2793  if (TMath::Abs(K0S->GetPdgCode()) != 310) isBkgLc = kTRUE;
2794  }
2795  }
2796  else{
2797  AliDebug(2, "The K0S is not true --> it does not have a label, continuing...");
2798  isBkgLc = kTRUE;
2799  }
2800 
2801  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
2802  if (isMCokLc) { // We can then access its MC info, and it might then be that also the Lc is a true Lc
2803  Int_t iK0 = K0S->GetMother();
2804  if (iK0 < 0) {
2805  Printf("The K0S has no mother... IMPOSSIBLE"); // the K0S MUST have a mother!
2806  }
2807  else { // The K0S has a mother
2808  daughv02Lc = dynamic_cast<AliAODMCParticle*>(mcArray->At(iK0));
2809  if (!daughv02Lc){
2810  AliDebug(3, "Could not access MC info for second daughter of Lc");
2811  }
2812  else { // we can access safely the K0S mother in the MC
2813  if( daughv01Lc && (daughv01Lc->GetMother() == daughv02Lc->GetMother()) && (daughv01Lc->GetMother()>=0) ){ // This is a true cascade! bachelor and V0 come from the same mother
2814  //Printf("Lc: The mother has label %d", daughv01Lc->GetMother());
2815  AliAODMCParticle *motherLc = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughv01Lc->GetMother()));
2816  Int_t pdgMum = 0, pdgBach = 0, pdgV0 = 0;
2817  if (motherLc) pdgMum = motherLc->GetPdgCode();
2818  if (daughv01Lc) pdgBach = daughv01Lc->GetPdgCode();
2819  if (daughv02Lc) pdgV0 = daughv02Lc->GetPdgCode();
2820  AliDebug(2, Form("pdgMum = %d, pdgBach = %d, pdgV0 = %d", pdgMum, pdgBach, pdgV0));
2821 
2822  if( motherLc && TMath::Abs(motherLc->GetPdgCode()) != 21 ){ // These are true cascades, we don't know yet if they are Lc
2823  fHistoMassLcTrue->Fill(massLc);
2824  fHistoDecayLengthLcTrue->Fill(decayLengthLc);
2825  fHistoLifeTimeLcTrue->Fill(lifeTimeLc);
2826  fHistoMassLcTrueFromAOD->Fill(cascade->InvMassLctoK0sP());
2827 
2828  fHistoMassV0fromLcTrue->Fill(massV0);
2829  fHistoDecayLengthV0fromLcTrue->Fill(decayLengthV0);
2830  fHistoLifeTimeV0fromLcTrue->Fill(lifeTimeV0);
2831 
2832  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...
2833  AliDebug(2, Form("IT IS SIGNAL!!! with label = %d", motherLc->GetLabel()));
2834 
2835  fHistoArmenterosPodolanskiV0KFSgn->Fill(qtAlphaV0[1], qtAlphaV0[0]);
2836  fHistoArmenterosPodolanskiV0AODSgn->Fill(v0part->AlphaV0(), v0part->PtArmV0());
2837 
2838  fHistoDistanceLcToPrimVtxSgn->Fill(distanceLcToPrimVtx);
2839  fHistoDistanceV0ToPrimVtxSgn->Fill(distanceV0ToPrimVtx);
2840  fHistoDistanceV0ToLcSgn->Fill(distanceV0ToLc);
2841 
2842  fHistoMassLcSgn->Fill(massLc);
2843  fHistoMassLcSgnFromAOD->Fill(cascade->InvMassLctoK0sP());
2844 
2845  fHistoDecayLengthLcSgn->Fill(decayLengthLc);
2846  fHistoLifeTimeLcSgn->Fill(lifeTimeLc);
2847 
2848  fHistoMassV0fromLcSgn->Fill(massV0);
2849  fHistoDecayLengthV0fromLcSgn->Fill(decayLengthV0);
2850  fHistoLifeTimeV0fromLcSgn->Fill(lifeTimeV0);
2851  }
2852  else {
2853  isBkgLc = kTRUE; // it is not a Lc, since the pdg != 4122
2854  }
2855 
2856  // if we got here, we can compare with MC information; this part is done also for cases where the cascade is not a Lc!
2857  // First, we compare the vtx of the Lc
2858  Double_t xLcMC = motherLc->Xv();
2859  Double_t yLcMC = motherLc->Yv();
2860  Double_t zLcMC = motherLc->Zv();
2861  //Printf("Got MC vtx for Lc");
2862  Double_t xProtonMC = daughv01Lc->Xv();
2863  Double_t yProtonMC = daughv01Lc->Yv();
2864  Double_t zProtonMC = daughv01Lc->Zv();
2865  //Printf("Got MC vtx for Proton");
2866 
2867  Double_t vtxLcResidualToPrimVtx = TMath::Sqrt((xLcMC - xProtonMC) * (xLcMC - xProtonMC) +
2868  (yLcMC - yProtonMC) * (yLcMC - yProtonMC) +
2869  (zLcMC - zProtonMC) * (zLcMC - zProtonMC)) - distanceLcToPrimVtx;
2870 
2871  // Then, we compare the vtx of the K0s
2872  Double_t xV0MC = motherV0->Xv(); //Production vertex of K0S --> Where Lc decays
2873  Double_t yV0MC = motherV0->Yv();
2874  Double_t zV0MC = motherV0->Zv();
2875 
2876  //Printf("Got MC vtx for V0");
2877  Double_t xPionMC = daughv01->Xv(); //Production vertex of Pion --> Where K0S decays
2878  Double_t yPionMC = daughv01->Yv();
2879  Double_t zPionMC = daughv01->Zv();
2880  //Printf("Got MC vtx for Pion");
2881  Printf("Vertices: MC: x = %f, y = %f, z = %f", xPionMC, yPionMC, zPionMC);
2882 
2883  Double_t vtxV0ResidualToLc = TMath::Sqrt((xV0MC - xPionMC) * (xV0MC - xPionMC) +
2884  (yV0MC - yPionMC) * (yV0MC - yPionMC) +
2885  (zV0MC - zPionMC) * (zV0MC - zPionMC)) - distanceV0ToLc;
2886 
2887  // Then, the K0S vertex wrt primary vertex
2888 
2889  Double_t vtxV0ResidualToPrimVtx = TMath::Sqrt((xPionMC - xLcMC) * (xPionMC - xLcMC) +
2890  (yPionMC - yLcMC) * (yPionMC - yLcMC) +
2891  (zPionMC - zLcMC) * (zPionMC - zLcMC)) - distanceV0ToPrimVtx;
2892 
2893  fHistoVtxLcResidualToPrimVtx->Fill(vtxLcResidualToPrimVtx);
2894  fHistoVtxV0ResidualToLc->Fill(vtxV0ResidualToLc);
2895  fHistoVtxV0ResidualToPrimVtx->Fill(vtxV0ResidualToPrimVtx);
2896 
2897  } //closing: if( motherLc && TMath::Abs(motherLc->GetPdgCode()) != 21 )
2898  else if (!motherLc){
2899  AliDebug(2, "We could not access MC info for Lc mother, so we did nothing");
2900  }
2901  else {
2902  AliDebug(2, "MC Lc mother is a gluon, so we did nothing");
2903  }
2904  } //closing: if( daughv01Lc->GetMother() == daughv02Lc->GetMother() && daughv01Lc->GetMother()>=0 )
2905  else {
2906  isBkgLc = kTRUE; // it cannot be a Lc, since the daughters do not have the same mother
2907  }
2908  } // closing: else { // we can access safely the K0S mother in the MC
2909  } // closing: else { // The K0S has a mother
2910  } // closing isMCLcok
2911  } // closing !isBkgLc
2912  } // closing fUseMCInfo
2913 
2914  //Printf("retMV0 = %d, retMLc = %d", retMV0, retMLc);
2915  if ( retMV0 == 0 && retMLc == 0){
2916  V0KF[0] = massV0;
2917  errV0KF[0] = sigmaMassV0;
2918  V0KF[1] = decayLengthV0;
2919  errV0KF[1] = sigmaDecayLengthV0;
2920  V0KF[2] = lifeTimeV0;
2921  errV0KF[2] = sigmaLifeTimeV0;
2922  LcKF[0] = massLc;
2923  errLcKF[0] = sigmaMassLc;
2924  LcKF[1] = decayLengthLc;
2925  errLcKF[1] = sigmaDecayLengthLc;
2926  LcKF[2] = lifeTimeLc;
2927  errLcKF[2] = sigmaLifeTimeLc;
2928  }
2929 
2930  return 1;
2931 
2932 }
2933 //________________________________________________________________________
2935  AliAODTrack* bachelor,
2936  TClonesArray *mcArray ){
2937 
2938  //Printf("In CheckBachelor");
2939 
2942 
2943  Int_t label = bachelor->GetLabel();
2944  if (label == -1) {
2945  return kBachFake;
2946  }
2947 
2948  AliAODMCParticle *mcpart = dynamic_cast<AliAODMCParticle*>(mcArray->At(TMath::Abs(label)));
2949  if (!mcpart) return kBachInvalid;
2950  Int_t pdg = mcpart->PdgCode();
2951  if (TMath::Abs(pdg) != 2212) {
2952  AliDebug(2, Form("Bachelor is not a p, but a particle with pdg code = %d", pdg));
2953  return kBachNoProton;
2954  }
2955  else { // it is a proton
2956  //Int_t labelLc = part->GetLabel();
2957  Int_t labelLc = FindLcLabel(part, mcArray);
2958  //Printf(">>>>>>>>>>>>> label for Lc = %d", labelLc);
2959  Int_t bachelorMotherLabel = mcpart->GetMother();
2960  //Printf(">>>>>>>>>>>>> label for bachelorMother = %d", bachelorMotherLabel);
2961  if (bachelorMotherLabel == -1) {
2962  return kBachPrimary;
2963  }
2964  AliAODMCParticle *bachelorMother = dynamic_cast<AliAODMCParticle*>(mcArray->At(bachelorMotherLabel));
2965  if (!bachelorMother) return kBachInvalid;
2966  Int_t pdgMother = bachelorMother->PdgCode();
2967  if (TMath::Abs(pdgMother) != 4122) {
2968  AliDebug(2, Form("The proton does not come from a Lc, but from a particle with pdgcode = %d", pdgMother));
2969  return kBachNoLambdaMother;
2970  }
2971  else { // it comes from Lc
2972  if (labelLc != bachelorMotherLabel){
2973  //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));
2974  AliDebug(2, Form("The proton comes from a Lc, but it is not the candidate we are analyzing"));
2976  }
2977  else { // it comes from the correct Lc
2978  AliDebug(2, Form("The proton comes from a Lc, and it is EXACTLY the candidate we are analyzing"));
2979  return kBachCorrectLambdaMother;
2980  }
2981  }
2982  }
2983 
2984  return kBachInvalid;
2985 
2986 }
2987 
2988 //________________________________________________________________________
2990  AliAODv0* v0part,
2991  //AliAODTrack* v0part,
2992  TClonesArray *mcArray ){
2993 
2996 
2997  //Printf(" CheckK0S");
2998 
2999  //Int_t labelMatchToMC = v0part->MatchToMC(310, mcArray);
3000  //Int_t label = v0part->GetLabel();
3001  Int_t labelFind = FindV0Label(v0part, mcArray);
3002  //Printf("\n >>>>>>>>>>>>> label for V0 = %d, from MatchToMC = %d, from Find = %d", label, labelMatchToMC, labelFind);
3003  if (labelFind == -1) {
3004  return kK0SFake;
3005  }
3006 
3007  AliAODMCParticle *mcpart = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelFind));
3008  if (!mcpart) return kK0SInvalid;
3009  Int_t pdg = mcpart->PdgCode();
3010  if (TMath::Abs(pdg) != 310) {
3011  AliDebug(2, Form("K0Spart is not a K0S, but a particle with pdg code = %d", pdg));
3012  //AliInfo(Form("K0Spart is not a K0S, but a particle with pdg code = %d", pdg));
3013  return kK0SNoK0S;
3014  }
3015  else { // it is a K0S
3016  //Int_t labelLc = part->GetLabel();
3017  Int_t labelLc = FindLcLabel(part, mcArray);
3018  Int_t K0SpartMotherLabel = mcpart->GetMother();
3019  if (K0SpartMotherLabel == -1) {
3020  return kK0SWithoutMother;
3021  }
3022  AliAODMCParticle *K0SpartMother = dynamic_cast<AliAODMCParticle*>(mcArray->At(K0SpartMotherLabel)); // should be a K0 in case of signal
3023  if (!K0SpartMother) return kK0SInvalid;
3024  Int_t pdgMotherK0S = K0SpartMother->PdgCode();
3025  if (TMath::Abs(pdgMotherK0S) != 311) {
3026  AliDebug(2, Form("The K0S does not come from a K0, but from a particle with pdgcode = %d", pdgMotherK0S));
3027  //AliInfo(Form("The K0S does not come from a K0, but from a particle with pdgcode = %d", pdgMotherK0S));
3028  return kK0SNotFromK0;
3029  }
3030  else { // the K0S comes from a K0
3031  Int_t K0MotherLabel = K0SpartMother->GetMother(); // mother of K0 --> Lc in case of signal
3032  if (K0MotherLabel == -1) {
3033  return kK0Primary;
3034  }
3035  AliAODMCParticle *K0Mother = dynamic_cast<AliAODMCParticle*>(mcArray->At(K0MotherLabel));
3036  if (!K0Mother) return kK0SInvalid;
3037  Int_t pdgK0Mother = K0Mother->PdgCode();
3038  if (TMath::Abs(pdgK0Mother) != 4122) { // the K0 does not come from a Lc
3039  AliDebug(2, Form("The K0 does not come from a Lc, but from a particle with pdgcode = %d", pdgK0Mother));
3040  //AliInfo(Form("The K0 does not come from a Lc, but from a particle with pdgcode = %d", pdgK0Mother));
3041  return kK0NoLambdaMother;
3042  }
3043  else { // the K0 comes from Lc
3044  if (labelLc != K0MotherLabel){ // The K0 comes from a different Lc
3045  AliDebug(2, Form("The K0S comes from a Lc, but it is not the candidate we are analyzing"));
3046  //AliInfo(Form("The K0S comes from a Lc, but it is not the candidate we are analyzing"));
3047  return kK0DifferentLambdaMother;
3048  }
3049  else { // it comes from the correct Lc
3050  AliDebug(2, Form("The K0S comes from a Lc, and it is EXACTLY the candidate we are analyzing"));
3051  //AliInfo(Form("The K0S comes from a Lc, and it is EXACTLY the candidate we are analyzing"));
3052  return kK0CorrectLambdaMother;
3053  }
3054  }
3055  }
3056  }
3057 
3058  return kK0SInvalid;
3059 
3060 }
3061 
3062 //----------------------------------------------------------------------------
3063 Int_t AliAnalysisTaskSELc2V0bachelorTMVAApp::FindV0Label(AliAODRecoDecay* v0part, TClonesArray *mcArray) const
3064 {
3065 
3066  //Printf(" FindV0Label");
3067 
3069 
3070  Int_t labMother[2] = {-1, -1};
3071  AliAODMCParticle *part = 0;
3072  AliAODMCParticle *mother = 0;
3073  Int_t dgLabels = -1;
3074 
3075  Int_t ndg = v0part->GetNDaughters();
3076  if (ndg != 2) {
3077  AliFatal(Form("IMPOSSIBLE!! There are %d daughters, but they should be 2!", ndg));
3078  }
3079 
3080  for(Int_t i = 0; i < 2; i++) {
3081  AliAODTrack *trk = (AliAODTrack*)v0part->GetDaughter(i);
3082  dgLabels = trk->GetLabel();
3083  if (dgLabels == -1) {
3084  //printf("daughter with negative label %d\n", dgLabels);
3085  return -1;
3086  }
3087  part = (AliAODMCParticle*)mcArray->At(TMath::Abs(dgLabels));
3088  if (!part) {
3089  //printf("no MC particle\n");
3090  return -1;
3091  }
3092  labMother[i] = part->GetMother();
3093  if (labMother[i] != -1){
3094  mother = (AliAODMCParticle*)mcArray->At(labMother[i]);
3095  if(!mother) {
3096  //printf("no MC mother particle\n");
3097  return -1;
3098  }
3099  }
3100  else {
3101  return -1;
3102  }
3103  }
3104 
3105  if (labMother[0] == labMother[1]) return labMother[0];
3106 
3107  else return -1;
3108 
3109 }
3110 
3111 //----------------------------------------------------------------------------
3113 {
3114 
3116 
3117  //Printf(" FindLcLabel");
3118 
3119  AliAODMCParticle *part = 0;
3120  AliAODMCParticle *mother = 0;
3121  AliAODMCParticle *grandmother = 0;
3122  Int_t dgLabels = -1;
3123 
3124  Int_t ndg = cascade->GetNDaughters();
3125  if (ndg != 2) {
3126  AliFatal(Form("IMPOSSIBLE!! There are %d daughters, but they should be 2!", ndg));
3127  }
3128 
3129  // daughter 0 --> bachelor, daughter 1 --> V0
3130  AliAODTrack *trk = (AliAODTrack*)cascade->GetDaughter(0); // bachelor
3131  dgLabels = trk->GetLabel();
3132  if (dgLabels == -1 ) {
3133  //printf("daughter with negative label %d\n", dgLabels);
3134  return -1;
3135  }
3136  part = (AliAODMCParticle*)mcArray->At(TMath::Abs(dgLabels));
3137  if (!part) {
3138  //printf("no MC particle\n");
3139  return -1;
3140  }
3141  Int_t labMotherBach = part->GetMother();
3142  if (labMotherBach == -1){
3143  return -1;
3144  }
3145  mother = (AliAODMCParticle*)mcArray->At(labMotherBach);
3146  if(!mother) {
3147  //printf("no MC mother particle\n");
3148  return -1;
3149  }
3150 
3151  AliAODv0 *v0 = (AliAODv0*)cascade->GetDaughter(1); // V0
3152  dgLabels = FindV0Label(v0, mcArray);
3153  if (dgLabels == -1 ) {
3154  //printf("daughter with negative label (v0 was a fake) %d\n", dgLabels);
3155  return -1;
3156  }
3157  part = (AliAODMCParticle*)mcArray->At(TMath::Abs(dgLabels));
3158  if (!part) {
3159  //printf("no MC particle\n");
3160  return -1;
3161  }
3162  Int_t labMotherv0 = part->GetMother();
3163  if (labMotherv0 == -1){
3164  return -1;
3165  }
3166  mother = (AliAODMCParticle*)mcArray->At(labMotherv0);
3167  if(!mother) {
3168  //printf("no MC mother particle\n");
3169  return -1;
3170  }
3171  Int_t labGrandMotherv0 = mother->GetMother();
3172  if (labGrandMotherv0 == -1){
3173  return -1;
3174  }
3175  grandmother = (AliAODMCParticle*)mcArray->At(labGrandMotherv0);
3176  if(!grandmother) {
3177  //printf("no MC mother particle\n");
3178  return -1;
3179  }
3180 
3181  //Printf("labMotherBach = %d, labMotherv0 = %d, labGrandMotherv0 = %d", labMotherBach, labMotherv0, labGrandMotherv0);
3182  if (labGrandMotherv0 == labMotherBach) return labMotherBach;
3183 
3184  else return -1;
3185 
3186 }
3187 
3188 //____________________________________________________________________________
3193 
3194  Int_t runNo = event->GetRunNumber();
3195  Int_t period = -1; // pp: 0-LHC10b, 1-LHC10c, 2-LHC10d, 3-LHC10e
3196 
3197  if(fYearNumber==10){
3198  if(runNo>114930 && runNo<117223) period = 0;//10b
3199  if(runNo>119158 && runNo<120830) period = 1;//10c
3200  if(runNo>122373 && runNo<126438) period = 2;//10d
3201  if(runNo>127711 && runNo<130851) period = 3;//10e
3202  if(period<0 || period>3) return 0;
3203  }else if(fYearNumber==16){
3204  if(runNo>=252235 && runNo<=252375)period = 0;//16d
3205  if(runNo>=252603 && runNo<=253591)period = 1;//16e
3206  if(runNo>=254124 && runNo<=254332)period = 2;//16g
3207  if(runNo>=254378 && runNo<=255469 )period = 3;//16h_1
3208  if(runNo>=254418 && runNo<=254422 )period = 4;//16h_2 negative mag
3209  if(runNo>=256146 && runNo<=256420 )period = 5;//16j
3210  if(runNo>=256504 && runNo<=258537 )period = 6;//16k
3211  if(runNo>=258883 && runNo<=260187)period = 7;//16l
3212  if(runNo>=262395 && runNo<=264035 )period = 8;//16o
3213  if(runNo>=264076 && runNo<=264347 )period = 9;//16p
3214  }else if(fYearNumber==17){
3215  if(runNo>=270822 && runNo<=270830)period = 0;//17e
3216  if(runNo>=270854 && runNo<=270865)period = 1;//17f
3217  if(runNo>=271868 && runNo<=273103)period = 2;//17h
3218  if(runNo>=273591 && runNo<=274442)period = 3;//17i
3219  if(runNo>=274593 && runNo<=274671)period = 4;//17j
3220  if(runNo>=274690 && runNo<=276508)period = 5;//17k
3221  if(runNo>=276551 && runNo<=278216)period = 6;//17l
3222  if(runNo>=278914 && runNo<=280140)period = 7;//17m
3223  if(runNo>=280282 && runNo<=281961)period = 8;//17o
3224  if(runNo>=282504 && runNo<=282704)period = 9;//17r
3225  }else if(fYearNumber==18){
3226  if(runNo>=285008 && runNo<=285447)period = 0;//18b
3227  if(runNo>=285978 && runNo<=286350)period = 1;//18d
3228  if(runNo>=286380 && runNo<=286937)period = 2;//18e
3229  if(runNo>=287000 && runNo<=287977)period = 3;//18f
3230  if(runNo>=288619 && runNo<=288750)period = 4;//18g
3231  if(runNo>=288804 && runNo<=288806)period = 5;//18h
3232  if(runNo>=288861 && runNo<=288909 )period = 6;//18i
3233  if(runNo==288943)period = 7;//18j
3234  if(runNo>=289165 && runNo<=289201)period = 8;//18k
3235  if(runNo>=289240 && runNo<=289971)period = 9;//18l
3236  if(runNo>=290222 && runNo<=292839)period = 10;//18m
3237  if(runNo>=293357 && runNo<=293359)period = 11;//18n
3238  if(runNo>=293368 && runNo<=293898)period = 12;//18o
3239  if(runNo>=294009 && runNo<=294925)period = 13;//18p
3240  }
3241  return fMultEstimatorAvg[period];
3242 
3243 }
TH1F * fHistoMCLcK0SpGenLimAcc
! histo with MC Lc –> K0S + p
TH2D * fHistoArmenterosPodolanskiV0AODSgn
! KF: AOD Armeteros-Podolanski plot for V0 from signal Lc from KF
Int_t fCurrentEvent
cut for KF on distance to primary vtx for V0
Int_t fNVars
flag to decide if to take the nSigma from the PIDresponse or from AliAODPidHF
Int_t pdg
AliAODTrack * GetProng(AliVEvent *event, AliAODRecoDecayHF *rd, Int_t iprong)
TH1D * fHistoVtxV0ResidualToLc
! KF: residual wrt MC of distance V0 vertex from Lc vertex (MC - KF)
TH2D * fHistoDecayLengthKFV0
! KF: decay length vs decay length error for V0 from KF
EBachelor CheckBachelor(AliAODRecoCascadeHF *part, AliAODTrack *bachelor, TClonesArray *mcArray)
double Double_t
Definition: External.C:58
Float_t fCentrality
tracklet multiplicity in event without eta cut
void StoreCandidates(AliVEvent *, Int_t nCand=0, Bool_t flagFilter=kTRUE)
TH2D * fHistoArmenterosPodolanskiV0KFSgn
! KF: Armeteros-Podolanski plot for V0 from signal Lc from KF
Definition: External.C:236
Int_t fNTracklets_All
tracklet multiplicity in event in [-1. 1]
Int_t GetnSigmaTOF(AliAODTrack *track, Int_t species, Double_t &sigma) const
virtual void UserCreateOutputObjects()
Implementation of interface methods.
TH1D * fHistoDecayLengthLcTrue
! KF: decay length for true cascades reconstructed with KF
TH1F * fHistoMCNch
! histogram with Nch distribution from MC production
Int_t FindLcLabel(AliAODRecoCascadeHF *cascade, TClonesArray *mcArray) const
Int_t FindV0Label(AliAODRecoDecay *v0part, TClonesArray *mcArray) const
AliPIDCombined * fPIDCombined
! combined PID response object
TH2F * fHistoCodesBkg
! histogram with codes for bachelor and V0 for background
Int_t MatchToMC(Int_t pdgabs, Int_t pdgabs2prong, Int_t *pdgDg, Int_t *pdgDg2prong, TClonesArray *mcArray, Bool_t isV0=kFALSE) const
AliAODTrack * Getv0NegativeTrack() const
TH1D * fHistoMassV0TrueFromAOD
! KF: AOD mass for true V0 reconstructed with KF
TH1D * fHistoDistanceV0ToLc
! KF: distance V0 vertex from Lc vertex
Int_t GetnSigmaTPC(AliAODTrack *track, Int_t species, Double_t &sigma) const
TH2F * fHistoTracklets_All_cent
! histogram with number of tracklets in the event in eta [-999, 999] vs centrality ...
TH2D * fBDTHistoVsSignd0
! BDT classifier vs V0 proton signed d0
TH1D * fHistoMassV0fromLcTrue
! KF: mass of V0 for true cascades reconstructed with KF
AliAODv0 * Getv0() const
TH1D * fHistoLifeTimeV0TrueK0S
! KF: life time for true V0 which are really K0S reconstructed with KF
TH1D * fHistoLifeTimeV0True
! KF: life time for true V0 reconstructed with KF
TH2D * fBDTHistoVsBachelorTPCP
! BDT classifier vs bachelor p at TPC wall
char Char_t
Definition: External.C:18
Double_t InvMassLctoLambdaPi() const
static Int_t CheckMatchingAODdeltaAODevents()
TH2D * fBDTHistoVsMassK0S
! BDT classifier vs mass (pi+pi-) pairs
TH1D * fHistoDistanceV0ToPrimVtx
! KF: distance V0 vertex from primary vertex
Int_t CallKFVertexing(AliAODRecoCascadeHF *cascade, AliAODv0 *v0part, AliAODTrack *bach, TClonesArray *mcArray, Double_t *V0KF, Double_t *errV0KF, Double_t *LcKF, Double_t *errLcKF, Double_t *distances, Double_t *armPolKF)
TList * fOutputKF
! User output1: list of histograms from KF
AliNormalizationCounter * fCounterC
! AliNormalizationCounter on output slot 4, corrected with multiplicity dependence ...
TH1D * fHistoDecayLengthV0fromLcSgn
! KF: decay length of V0 for signal Lc reconstructed with KF
TH1D * fHistoMassV0fromLcAll
! KF: mass of V0 for all cascades reconstructed with KF
TH2D * fBDTHistoVsCosPAK0S
! BDT classifier vs V0 cosine of pointing angle
virtual Int_t PreSelect(TObjArray aodtracks)
TH1D * fHistoMassV0All
! KF: mass for all V0 reconstructed with KF
TH2D * fHistoMassKFV0
! KF: mass vs mass error for V0 from KF
TH2D * fBDTHistoVsnSigmaTOFpr
! BDT classifier vs nSigmaTOFpr
TH1F * fHistoEvents
! histogram with number of events analyzed
TH1D * fHistoDistanceLcToPrimVtxSgn
! KF: distance of signal Lc vertex from primary vertex
static Int_t CheckOrigin(TClonesArray *arrayMC, AliAODMCParticle *mcPart, Bool_t searchUpToQuark=kTRUE)
Functions to check the decay tree.
void FillLc2pK0Sspectrum(AliAODRecoCascadeHF *part, Int_t isLc, Int_t &nSelectedAnal, AliRDHFCutsLctoV0 *cutsAnal, TClonesArray *mcArray, Int_t iLctopK0s)
histos
TH2D * fHistoMassKFLc
! KF: mass vs mass error for Lc from KF
Double_t GetMaxVtxZ() const
Definition: AliRDHFCuts.h:274
TH2F * fHistoVzVsNtrCorr
! hist. Vz vs corrected tracklets
Bool_t FillRecoCasc(AliVEvent *event, AliAODRecoCascadeHF *rc, Bool_t isDStar, Bool_t recoSecVtx=kFALSE)
Double_t InvMassLctoK0sP() const
TH2F * fHistoTracklets_1_cent
! histogram with number of tracklets in the event in eta [-1, 1] vs centrality
AliAODPidHF * GetPidHF() const
Definition: AliRDHFCuts.h:264
TH2D * fBDTHistoVsCombinedProtonProb
! BDT classifier vs combined proton probability
TH2F * fHistoCodesSgn
! histogram with codes for bachelor and V0 for signal
AliAODTrack * Getv0PositiveTrack() const
TH1F * fHistoCentrality
! histogram with centrality from AliRDHFCuts
void MakeAnalysisForLc2prK0S(AliAODEvent *aodEvent, TClonesArray *arrayLctopK0s, TClonesArray *mcArray, Int_t &nSelectedAnal, AliRDHFCutsLctoV0 *cutsAnal, TClonesArray *array3Prong, AliAODMCHeader *aodheader)
TH1D * fHistoMassLcTrue
! KF: mass for true cascades reconstructed with KF
static Double_t GetVZEROCEqualizedMultiplicity(AliAODEvent *ev)
TH2D * fBDTHistoVsCosThetaStar
! BDT classifier vs proton emission angle in pK0s pair rest frame
void SetStudyMultiplicity(Bool_t flag, Float_t etaRange)
TH1D * fHistoDecayLengthV0All
! KF: decay length for all V0 reconstructed with KF
static Int_t GetNumberOfTrackletsInEtaRange(AliAODEvent *ev, Double_t mineta, Double_t maxeta)
int Int_t
Definition: External.C:63
TH2D * fHistoLifeTimeKFV0
! KF: life time vs life time error for V0 from KF
TH1D * fHistoVtxLcResidualToPrimVtx
! KF: residual wrt MC of distance Lc vertex from primary vertex (MC - KF)
unsigned int UInt_t
Definition: External.C:33
TH1D * fHistoLifeTimeV0All
! KF: life time for all V0 reconstructed with KF
TString fNamesTMVAVar
Pt bin that will be in the library to be loaded for the TMVA.
TH1D * fHistoKFV0
! KF: V0 code from KF (mass, decaylength, lifetime considered)
float Float_t
Definition: External.C:68
TH1D * fHistoMassV0True
! KF: mass for true V0 reconstructed with KF
TH1F * fHistoTracklets_All
! histogram with number of tracklets in the event in eta [-999, 999]
TH1F * fHistoNtrCorr
! hist. with number of corrected tracklets
TF1 * fFuncWeightFONLL5overLHC13d3
! weight function for FONLL vs pPb prod.
TH1D * fHistoDecayLengthV0fromLcAll
! KF: decay length of V0 for all cascades reconstructed with KF
Double_t fBField
current event number - for debug purposes
TH1F * fHistoFiducialAcceptance
! histogram to check FiducialAcceptance cut
TTree * fVariablesTreeBkg
! tree of the candidate variables after track selection (Background)
Definition: External.C:228
Definition: External.C:212
AliAODTrack * GetBachelor() const
TH1D * fHistoLifeTimeV0fromLcSgn
! KF: life time of V0 for signal Lc reconstructed with KF
TH2D * fHistoArmenterosPodolanskiV0KF
! KF: Armeteros-Podolanski plot for all V0 from KF
Int_t GetUsePreselect()
Definition: AliRDHFCuts.h:398
Float_t GetCentrality(AliAODEvent *aodEvent)
Definition: AliRDHFCuts.h:275
TH1D * fHistoMassV0TrueK0SFromAOD
! KF: AOD mass for true V0 which are really K0S reconstructed with KF
UShort_t GetProngID(Int_t ip) const
TH1D * fHistoMassV0fromLcSgn
! KF: mass of V0 for signal Lc reconstructed with KF
static Bool_t HasCascadeCandidateAnyDaughInjected(AliAODRecoCascadeHF *cand, AliAODMCHeader *header, TClonesArray *arrayMC)
Bool_t fUsePIDresponseForNsigma
-1: no protection, 0: check AOD/dAOD nEvents only, 1: check AOD/dAOD nEvents + TProcessID names ...
TH1F * fHistoLcpKpiBeforeCuts
! histogram number of true Lc–>pKpi (3 prong) before any cut
TH1F * fHistoBackground
AliVertexingHFUtils used to check the generator of a specific candidate.
TH1D * fHistoMassLcAll
! KF: mass for all Lc reconstructed with KF
TH1D * fHistoLifeTimeV0fromLcTrue
! KF: life time of V0 for true cascades reconstructed with KF
Bool_t fUseOnTheFlyV0
list of profiles for z-vtx correction of multiplicity
TF1 * fFuncWeightPythia
mask to the trigger word returned by the physics selection
TString fTMVAlibPtBin
Name of the library to load to have the TMVA weights.
TH2D * fBDTHistoVsBachelorP
! BDT classifier vs bachelor p
TH2F * fHistoVzVsNtrUnCorr
! hist. Vz vs UNCORRECTED tracklets
Double_t CosV0PointingAngle() const
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
short Short_t
Definition: External.C:23
TH1F * fHistoNtrUnCorr
! hist. with number of uncorrected tracklets
TH1D * fHistoMassV0TrueK0S
! KF: mass for true V0 which are really K0S reconstructed with KF
TH1D * fHistoMassLcSgnFromAOD
! KF: AOD mass of signal Lc reconstructed with KF
TH1D * fHistoDecayLengthV0fromLcTrue
! KF: decay length of V0 for true cascades reconstructed with KF
static Double_t GetCorrectedNtracklets(TProfile *estimatorAvg, Double_t uncorrectedNacc, Double_t vtxZ, Double_t refMult)
TH1D * fHistoDecayLengthLcSgn
! KF: decay length of signal Lc reconstructed with KF
Bool_t IsEventSelected(AliVEvent *event)
TH1D * fHistoDecayLengthV0TrueK0S
! KF: decay length for true V0 which are really K0S reconstructed with KF
void StoreEvent(AliVEvent *, AliRDHFCuts *, Bool_t mc=kFALSE, Int_t multiplicity=-9999, Double_t spherocity=-99.)
Int_t fAODProtection
flag to decide whether or not to have extra histograms (useful mainly for debug)
AliVertexingHFUtils * fUtils
flag to fill bkg with only candidates that have daughters generated by HIJING (to be used for enriche...
TH2D * fBDTHistoVsnSigmaTPCpi
! BDT classifier vs nSigmaTPCpi
TH2D * fHistoKF
! KF: V0 code vs Lc code from KF (mass, decaylength, lifetime considered)
TH2D * fHistoArmenterosPodolanskiV0AOD
! KF: AOD Armeteros-Podolanski plot for all V0 from KF
TH1D * fHistoKFLc
! KF: Lc code from KF (mass, decaylength, lifetime considered)
AliPIDResponse * fPIDResponse
! PID response object
TH1D * fHistoMassLcSgn
! KF: mass of signal Lc reconstructed with KF
Float_t fCutKFDeviationFromVtxV0
cut for KF on distance to primary vtx
Bool_t CheckCascadeFlags(AliRDHFCuts::ESele selFlag=AliRDHFCuts::kLctoV0Cuts)
Bool_t fIsEventSelected
flag to analyze also on-the-fly V0 candidates
AliNormalizationCounter * fCounter
switch between Lpi and K0sp
TH1F * fHistoMCLcK0SpGen
flag to allow to use only PYTHIA tracks for background
TH1D * fHistoMassLcTrueFromAOD
! KF: AOD mass for true cascades reconstructed with KF
Float_t * fCandidateVariables
! variables to be written to the tree
TH1D * fHistoDistanceV0ToLcSgn
! KF: distance for signal Lc of V0 vertex from Lc vertex
TH1D * fHistoLifeTimeLcSgn
! KF: life time of signal Lc reconstructed with KF
Bool_t fCallKFVertexing
flag to use topological constraints in KF
TH1D * fHistoLifeTimeLcAll
! KF: life time for all Lc reconstructed with KF
TH2D * fBDTHistoTMVA
! BDT histo file for the case in which the xml file is used
TF1 * fFuncWeightFONLL5overLHC13d3Lc
! weight function for FONLL vs pPb prod.
const char Option_t
Definition: External.C:48
TH1D * fHistoLifeTimeLcTrue
! KF: life time for true cascades reconstructed with KF
TH2D * fBDTHistoVsnSigmaTPCpr
! BDT classifier vs nSigmaTPCpr
TH1F * fHistoLcOnTheFly
! histogram with number of Lc with on-the-fly V0
EK0S CheckK0S(AliAODRecoCascadeHF *part, AliAODv0 *v0part, TClonesArray *mcArray)
AliAODVertex * GetPrimaryVtx() const
TH2D * fBDTHistoVstImpParBach
! BDT classifier vs proton d0
static Double_t GetVZEROAEqualizedMultiplicity(AliAODEvent *ev)
Utilities for V0 multiplicity checks.
TH1D * fHistoDecayLengthV0True
! KF: decay length for true V0 reconstructed with KF
bool Bool_t
Definition: External.C:53
TH2D * fBDTHistoVsBachelorPt
! BDT classifier vs proton pT
Double_t DecayLengthV0() const
void SetTriggerClass(TString trclass0, TString trclass1="")
Definition: AliRDHFCuts.h:208
virtual Bool_t IsInFiducialAcceptance(Double_t pt, Double_t y) const
TH1D * fHistoLifeTimeV0fromLcAll
! KF: life time of V0 for all cascades reconstructed with KF
TH1F * fHistoLcBeforeCuts
flag to fill only signal (speeding up processing)
Bool_t fKeepingOnlyPYTHIABkg
magnetic field of current event
virtual double GetMvaValue(const std::vector< double > &inputValues) const =0
TH2D * fHistoDecayLengthKFLc
! KF: decay length vs decay length error for Lc from KF
Bool_t fUseWeightsLibrary
flag to decide whether to fill the sgn and bkg trees
TH2D * fHistoLifeTimeKFLc
! KF: life time vs life time error for Lc from KF
Int_t IsSelectedSingleCut(TObject *obj, Int_t selectionLevel, Int_t cutIndex, AliAODEvent *aod=0x0)
Bool_t fKeepingOnlyHIJINGBkg
flag to decide whether to call or not KF
TH1F * fHistoMCLcK0SpGenAcc
! histo with MC Lc –> K0S + p
TH1F * fHistoTracklets_1
! histogram with number of tracklets in the event in eta [-1, 1]
virtual Int_t IsSelected(TObject *obj, Int_t selectionLevel)
TH2D * fBDTHistoVsnSigmaTPCka
! BDT classifier vs nSigmaTPCka
void SetTriggerMask(ULong64_t mask=0)
Definition: AliRDHFCuts.h:71
TH1D * fHistoVtxV0ResidualToPrimVtx
! KF: residual wrt MC of distance V0 vertex from primary vertex (MC - KF)
TH2D * fBDTHisto
vector of the names of the input variables
TH1D * fHistoDistanceV0ToPrimVtxSgn
! KF: distance for signal Lc of V0 vertex from primary vertex
IClassifierReader * fBDTReader
! BDT reader using BDT class
TH1D * fHistoDecayLengthLcAll
! KF: decay length for all Lc reconstructed with KF
Class with functions useful for different D2H analyses //.