AliPhysics  608b256 (608b256)
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 
74 #include "IClassifierReader.h"
75 
76 using std::cout;
77 using std::endl;
78 
79 #include <dlfcn.h>
80 
84 
85 //__________________________________________________________________________
88  fUseMCInfo(kFALSE),
89  fOutput(0),
90  fCEvents(0),
91  fPIDResponse(0),
92  fPIDCombined(0),
93  fIsK0sAnalysis(kFALSE),
94  fCounter(0),
95  fAnalCuts(0),
96  fListCuts(0),
97  fListWeight(0),
98  fListCounters(0),
99  fListProfiles(0),
100  fUseOnTheFlyV0(kFALSE),
101  fIsEventSelected(kFALSE),
102  fVariablesTreeSgn(0),
103  fVariablesTreeBkg(0),
104  fCandidateVariables(),
105  fHistoCentrality(0),
106  fHistoEvents(0),
107  fHistoTracklets_1(0),
108  fHistoTracklets_1_cent(0),
109  fHistoTracklets_All(0),
110  fHistoTracklets_All_cent(0),
111  fHistoLc(0),
112  fHistoLcOnTheFly(0),
113  fFillOnlySgn(kFALSE),
114  fHistoLcBeforeCuts(0),
115  fHistoFiducialAcceptance(0),
116  fHistoCodesSgn(0),
117  fHistoCodesBkg(0),
118  fHistoLcpKpiBeforeCuts(0),
119  fVtx1(0),
120  fHistoDistanceLcToPrimVtx(0),
121  fHistoDistanceV0ToPrimVtx(0),
122  fHistoDistanceV0ToLc(0),
123  fHistoDistanceLcToPrimVtxSgn(0),
124  fHistoDistanceV0ToPrimVtxSgn(0),
125  fHistoDistanceV0ToLcSgn(0),
126  fHistoVtxLcResidualToPrimVtx(0),
127  fHistoVtxV0ResidualToPrimVtx(0),
128  fHistoVtxV0ResidualToLc(0),
129  fHistoMassV0All(0),
130  fHistoDecayLengthV0All(0),
131  fHistoLifeTimeV0All(0),
132  fHistoMassV0True(0),
133  fHistoDecayLengthV0True(0),
134  fHistoLifeTimeV0True(0),
135  fHistoMassV0TrueFromAOD(0),
136  fHistoMassV0TrueK0S(0),
137  fHistoDecayLengthV0TrueK0S(0),
138  fHistoLifeTimeV0TrueK0S(0),
139  fHistoMassV0TrueK0SFromAOD(0),
140  fHistoMassLcAll(0),
141  fHistoDecayLengthLcAll(0),
142  fHistoLifeTimeLcAll(0),
143  fHistoMassLcTrue(0),
144  fHistoDecayLengthLcTrue(0),
145  fHistoLifeTimeLcTrue(0),
146  fHistoMassLcTrueFromAOD(0),
147  fHistoMassV0fromLcAll(0),
148  fHistoDecayLengthV0fromLcAll(0),
149  fHistoLifeTimeV0fromLcAll(0),
150  fHistoMassV0fromLcTrue(0),
151  fHistoDecayLengthV0fromLcTrue(0),
152  fHistoLifeTimeV0fromLcTrue(0),
153  fHistoMassLcSgn(0),
154  fHistoMassLcSgnFromAOD(0),
155  fHistoDecayLengthLcSgn(0),
156  fHistoLifeTimeLcSgn(0),
157  fHistoMassV0fromLcSgn(0),
158  fHistoDecayLengthV0fromLcSgn(0),
159  fHistoLifeTimeV0fromLcSgn(0),
160  fHistoKF(0),
161  fHistoKFV0(0),
162  fHistoKFLc(0),
163  fHistoMassKFV0(0),
164  fHistoDecayLengthKFV0(0),
165  fHistoLifeTimeKFV0(0),
166  fHistoMassKFLc(0),
167  fHistoDecayLengthKFLc(0),
168  fHistoLifeTimeKFLc(0),
169  fHistoArmenterosPodolanskiV0KF(0),
170  fHistoArmenterosPodolanskiV0KFSgn(0),
171  fHistoArmenterosPodolanskiV0AOD(0),
172  fHistoArmenterosPodolanskiV0AODSgn(0),
173  fOutputKF(0),
174  fmcLabelLc(-1),
175  ftopoConstraint(kTRUE),
176  fCallKFVertexing(kFALSE),
177  fKeepingOnlyHIJINGBkg(kFALSE),
178  fUtils(0),
179  fHistoBackground(0),
180  fCutKFChi2NDF(999999.),
181  fCutKFDeviationFromVtx(999999.),
182  fCutKFDeviationFromVtxV0(0.),
183  fCurrentEvent(-1),
184  fBField(0),
185  fKeepingOnlyPYTHIABkg(kFALSE),
186  fHistoMCLcK0SpGen(0),
187  fHistoMCLcK0SpGenAcc(0),
188  fHistoMCLcK0SpGenLimAcc(0),
189  fTriggerMask(0),
190  fFuncWeightPythia(0),
191  fFuncWeightFONLL5overLHC13d3(0),
192  fFuncWeightFONLL5overLHC13d3Lc(0),
193  fHistoMCNch(0),
194  fNTracklets_1(0),
195  fNTracklets_All(0),
196  fCentrality(0),
197  fFillTree(0),
198  fBDTReader(0),
199  fTMVAlibName(""),
200  fTMVAlibPtBin(""),
201  fNamesTMVAVar(""),
202  fBDTHisto(0),
203  fBDTHistoVsMassK0S(0),
204  fBDTHistoVstImpParBach(0),
205  fBDTHistoVstImpParV0(0),
206  fBDTHistoVsBachelorPt(0),
207  fBDTHistoVsCombinedProtonProb(0),
208  fBDTHistoVsCtau(0),
209  fBDTHistoVsCosPAK0S(0),
210  fBDTHistoVsSignd0(0),
211  fBDTHistoVsCosThetaStar(0),
212  fHistoNsigmaTPC(0),
213  fHistoNsigmaTOF(0),
214  fDebugHistograms(kFALSE)
215 {
217  //
218 }
219 //___________________________________________________________________________
221  AliRDHFCutsLctoV0* analCuts, Bool_t useOnTheFly) :
222  AliAnalysisTaskSE(name),
223  fUseMCInfo(kFALSE),
224  fOutput(0),
225  fCEvents(0),
226  fPIDResponse(0),
227  fPIDCombined(0),
228  fIsK0sAnalysis(kFALSE),
229  fCounter(0),
230  fAnalCuts(analCuts),
231  fListCuts(0),
232  fListWeight(0),
233  fListCounters(0),
234  fListProfiles(0),
235  fUseOnTheFlyV0(useOnTheFly),
236  fIsEventSelected(kFALSE),
240  fHistoCentrality(0),
241  fHistoEvents(0),
246  fHistoLc(0),
247  fHistoLcOnTheFly(0),
248  fFillOnlySgn(kFALSE),
251  fHistoCodesSgn(0),
252  fHistoCodesBkg(0),
254  fVtx1(0),
264  fHistoMassV0All(0),
267  fHistoMassV0True(0),
275  fHistoMassLcAll(0),
278  fHistoMassLcTrue(0),
288  fHistoMassLcSgn(0),
295  fHistoKF(0),
296  fHistoKFV0(0),
297  fHistoKFLc(0),
298  fHistoMassKFV0(0),
301  fHistoMassKFLc(0),
308  fOutputKF(0),
309  fmcLabelLc(-1),
310  ftopoConstraint(kTRUE),
311  fCallKFVertexing(kFALSE),
312  fKeepingOnlyHIJINGBkg(kFALSE),
313  fUtils(0),
314  fHistoBackground(0),
315  fCutKFChi2NDF(999999.),
316  fCutKFDeviationFromVtx(999999.),
318  fCurrentEvent(-1),
319  fBField(0),
320  fKeepingOnlyPYTHIABkg(kFALSE),
324  fTriggerMask(0),
328  fHistoMCNch(0),
329  fNTracklets_1(0),
330  fNTracklets_All(0),
331  fCentrality(0),
332  fFillTree(0),
333  fBDTReader(0),
334  fTMVAlibName(""),
335  fTMVAlibPtBin(""),
336  fNamesTMVAVar(""),
337  fBDTHisto(0),
343  fBDTHistoVsCtau(0),
347  fHistoNsigmaTPC(0),
348  fHistoNsigmaTOF(0),
349  fDebugHistograms(kFALSE)
350 {
351  //
353  //
354  Info("AliAnalysisTaskSELc2V0bachelorTMVAApp","Calling Constructor");
355 
356  DefineOutput(1, TList::Class()); // Tree signal + Tree Bkg + histoEvents
357  DefineOutput(2, TList::Class()); // normalization counter object
358  DefineOutput(3, TList::Class()); // Cuts
359  DefineOutput(4, TTree::Class()); // Tree signal + Tree Bkg + histoEvents
360  DefineOutput(5, TTree::Class()); // Tree signal + Tree Bkg + histoEvents
361  DefineOutput(6, TList::Class()); // Tree signal + Tree Bkg + histoEvents
362  DefineOutput(7, TList::Class()); // weights
363 }
364 
365 //___________________________________________________________________________
367  //
369  //
370  Info("~AliAnalysisTaskSELc2V0bachelorTMVAApp","Calling Destructor");
371 
372  if (fOutput) {
373  delete fOutput;
374  fOutput = 0;
375  }
376 
377  if (fPIDResponse) {
378  delete fPIDResponse;
379  }
380 
381  if (fPIDCombined) {
382  delete fPIDCombined;
383  }
384 
385  if (fCounter) {
386  delete fCounter;
387  fCounter = 0;
388  }
389 
390  if (fAnalCuts) {
391  delete fAnalCuts;
392  fAnalCuts = 0;
393  }
394 
395  if (fListCuts) {
396  delete fListCuts;
397  fListCuts = 0;
398  }
399 
400  if (fListCounters) {
401  delete fListCounters;
402  fListCounters = 0;
403  }
404 
405  if (fListWeight) {
406  delete fListWeight;
407  fListWeight = 0;
408  }
409 
410  if(fVariablesTreeSgn){
411  delete fVariablesTreeSgn;
412  fVariablesTreeSgn = 0;
413  }
414 
415  if(fVariablesTreeBkg){
416  delete fVariablesTreeBkg;
417  fVariablesTreeBkg = 0;
418  }
419 
420  if (fOutputKF) {
421  delete fOutputKF;
422  fOutputKF = 0;
423  }
424 
425  if (fUtils) {
426  delete fUtils;
427  fUtils = 0;
428  }
429 
430  if (fBDTReader) {
431  //delete fBDTReader;
432  fBDTReader = 0;
433  }
434 
435 }
436 //_________________________________________________
438  //
440  //
441 
442  fIsEventSelected=kFALSE;
443 
444  if (fDebug > 1) AliInfo("Init");
445 
446  fListCuts = new TList();
447  fListCuts->SetOwner();
449  PostData(3, fListCuts);
450 
451  // Save the weight functions or histograms
452  fListWeight = new TList();
453  fListWeight->SetOwner();
454  fListWeight->Add(fHistoMCNch);
455  PostData(7, fListWeight);
456 
458 
459  return;
460 }
461 
462 //________________________________________ terminate ___________________________
464 {
468 
469  AliInfo("Terminate");
470  AliAnalysisTaskSE::Terminate();
471 
472  fOutput = dynamic_cast<TList*> (GetOutputData(1));
473  if (!fOutput) {
474  AliError("fOutput not available");
475  return;
476  }
477 
478  if(fHistoMCLcK0SpGen) {
479  AliInfo(Form("At MC level, %f Lc --> K0S + p were found", fHistoMCLcK0SpGen->GetEntries()));
480  } else {
481  AliInfo("fHistoMCLcK0SpGen not available");
482  }
484  AliInfo(Form("At MC level, %f Lc --> K0S + p were found in the acceptance", fHistoMCLcK0SpGenAcc->GetEntries()));
485  } else {
486  AliInfo("fHistoMCLcK0SpGenAcc not available");
487  }
488  if(fVariablesTreeSgn) {
489  AliInfo(Form("At Reco level, %lld Lc --> K0S + p were found", fVariablesTreeSgn->GetEntries()));
490  } else {
491  AliInfo("fVariablesTreeSgn not available");
492  }
493 
494  fOutputKF = dynamic_cast<TList*> (GetOutputData(6));
495  if (!fOutputKF) {
496  AliError("fOutputKF not available");
497  return;
498  }
499 
500  return;
501 }
502 
503 //___________________________________________________________________________
506  AliInfo(Form("CreateOutputObjects of task %s\n", GetName()));
507 
508  fOutput = new TList();
509  fOutput->SetOwner();
510  fOutput->SetName("listTrees");
511 
512  // Output slot 1: list of 2 trees (Sgn + Bkg) of the candidate variables
513  const char* nameoutput = GetOutputSlot(1)->GetContainer()->GetName();
514  fVariablesTreeSgn = new TTree(Form("%s_Sgn", nameoutput), "Candidates variables tree, Signal");
515  fVariablesTreeBkg = new TTree(Form("%s_Bkg", nameoutput), "Candidates variables tree, Background");
516 
517  Int_t nVar;
518  if (fUseMCInfo) nVar = 48; //"full" tree if MC
519  else nVar = 33; //"reduced" tree if data
520 
521  fCandidateVariables = new Float_t [nVar];
522  TString * fCandidateVariableNames = new TString[nVar];
523 
524  if (fUseMCInfo) { // "full tree" for MC
525  fCandidateVariableNames[0] = "massLc2K0Sp";
526  fCandidateVariableNames[1] = "massLc2Lambdapi";
527  fCandidateVariableNames[2] = "massK0S";
528  fCandidateVariableNames[3] = "massLambda";
529  fCandidateVariableNames[4] = "massLambdaBar";
530  fCandidateVariableNames[5] = "cosPAK0S";
531  fCandidateVariableNames[6] = "dcaV0";
532  fCandidateVariableNames[7] = "tImpParBach";
533  fCandidateVariableNames[8] = "tImpParV0";
534  fCandidateVariableNames[9] = "nSigmaTPCpr";
535  fCandidateVariableNames[10] = "nSigmaTOFpr";
536  fCandidateVariableNames[11] = "bachelorPt";
537  fCandidateVariableNames[12] = "V0positivePt";
538  fCandidateVariableNames[13] = "V0negativePt";
539  fCandidateVariableNames[14] = "dcaV0pos";
540  fCandidateVariableNames[15] = "dcaV0neg";
541  fCandidateVariableNames[16] = "v0Pt";
542  fCandidateVariableNames[17] = "massGamma";
543  fCandidateVariableNames[18] = "LcPt";
544  fCandidateVariableNames[19] = "combinedProtonProb";
545  fCandidateVariableNames[20] = "LcEta";
546  fCandidateVariableNames[21] = "V0positiveEta";
547  fCandidateVariableNames[22] = "V0negativeEta";
548  fCandidateVariableNames[23] = "TPCProtonProb";
549  fCandidateVariableNames[24] = "TOFProtonProb";
550  fCandidateVariableNames[25] = "bachelorEta";
551  fCandidateVariableNames[26] = "LcP";
552  fCandidateVariableNames[27] = "bachelorP";
553  fCandidateVariableNames[28] = "v0P";
554  fCandidateVariableNames[29] = "V0positiveP";
555  fCandidateVariableNames[30] = "V0negativeP";
556  fCandidateVariableNames[31] = "v0Eta";
557  fCandidateVariableNames[32] = "DecayLengthLc";
558  fCandidateVariableNames[33] = "DecayLengthK0S";
559  fCandidateVariableNames[34] = "bachCode";
560  fCandidateVariableNames[35] = "k0SCode";
561  fCandidateVariableNames[36] = "alphaArm";
562  fCandidateVariableNames[37] = "ptArm";
563  fCandidateVariableNames[38] = "CosThetaStar";
564  fCandidateVariableNames[39] = "weightPtFlat";
565  fCandidateVariableNames[40] = "weightFONLL5overLHC13d3";
566  fCandidateVariableNames[41] = "weightFONLL5overLHC13d3Lc";
567  fCandidateVariableNames[42] = "weightNch";
568  fCandidateVariableNames[43] = "NtrkRaw";
569  fCandidateVariableNames[44] = "NtrkCorr";
570  fCandidateVariableNames[45] = "signd0";
571  fCandidateVariableNames[46] = "centrality";
572  fCandidateVariableNames[47] = "NtrkAll";
573  }
574  else { // "light mode"
575  fCandidateVariableNames[0] = "massLc2K0Sp";
576  fCandidateVariableNames[1] = "massLc2Lambdapi";
577  fCandidateVariableNames[2] = "massK0S";
578  fCandidateVariableNames[3] = "massLambda";
579  fCandidateVariableNames[4] = "massLambdaBar";
580  fCandidateVariableNames[5] = "cosPAK0S";
581  fCandidateVariableNames[6] = "dcaV0";
582  fCandidateVariableNames[7] = "tImpParBach";
583  fCandidateVariableNames[8] = "tImpParV0";
584  fCandidateVariableNames[9] = "nSigmaTPCpr";
585  fCandidateVariableNames[10] = "nSigmaTOFpr";
586  fCandidateVariableNames[11] = "bachelorPt";
587  fCandidateVariableNames[12] = "V0positivePt";
588  fCandidateVariableNames[13] = "V0negativePt";
589  fCandidateVariableNames[14] = "dcaV0pos";
590  fCandidateVariableNames[15] = "dcaV0neg";
591  fCandidateVariableNames[16] = "v0Pt";
592  fCandidateVariableNames[17] = "massGamma";
593  fCandidateVariableNames[18] = "LcPt";
594  fCandidateVariableNames[19] = "combinedProtonProb";
595  fCandidateVariableNames[20] = "V0positiveEta";
596  fCandidateVariableNames[21] = "V0negativeEta";
597  fCandidateVariableNames[22] = "bachelorEta";
598  fCandidateVariableNames[23] = "v0P";
599  fCandidateVariableNames[24] = "DecayLengthK0S";
600  fCandidateVariableNames[25] = "alphaArm";
601  fCandidateVariableNames[26] = "ptArm";
602  fCandidateVariableNames[27] = "NtrkRaw";
603  fCandidateVariableNames[28] = "NtrkCorr";
604  fCandidateVariableNames[29] = "CosThetaStar";
605  fCandidateVariableNames[30] = "signd0";
606  fCandidateVariableNames[31] = "centrality";
607  fCandidateVariableNames[32] = "NtrkAll";
608  }
609 
610  for(Int_t ivar=0; ivar < nVar; ivar++){
611  fVariablesTreeSgn->Branch(fCandidateVariableNames[ivar].Data(), &fCandidateVariables[ivar], Form("%s/f",fCandidateVariableNames[ivar].Data()));
612  fVariablesTreeBkg->Branch(fCandidateVariableNames[ivar].Data(), &fCandidateVariables[ivar], Form("%s/f",fCandidateVariableNames[ivar].Data()));
613  }
614 
615  fHistoCentrality = new TH1F("fHistoCentrality", "fHistoCentrality", 100, 0., 100.);
616 
617  fHistoEvents = new TH1F("fHistoEvents", "fHistoEvents", 2, -0.5, 1.5);
618  TString labelEv[2] = {"NotSelected", "Selected"};
619  for (Int_t ibin = 1; ibin <= fHistoEvents->GetNbinsX(); ibin++){
620  fHistoEvents->GetXaxis()->SetBinLabel(ibin, labelEv[ibin-1].Data());
621  }
622 
623  fHistoTracklets_1 = new TH1F("fHistoTracklets_1", "fHistoTracklets_1", 1000, 0, 5000);
624  fHistoTracklets_1_cent = new TH2F("fHistoTracklets_1_cent", "fHistoTracklets_1_cent; centrality; SPD tracklets [-1, 1]", 100, 0., 100., 1000, 0, 5000);
625  fHistoTracklets_All = new TH1F("fHistoTracklets_All", "fHistoTracklets_All", 1000, 0, 5000);
626  fHistoTracklets_All_cent = new TH2F("fHistoTracklets_All_cent", "fHistoTracklets_All_cent; centrality; SPD tracklets [-999, 999]", 100, 0., 100., 1000, 0, 5000);
627 
628  fHistoLc = new TH1F("fHistoLc", "fHistoLc", 2, -0.5, 1.5);
629 
630  fHistoLcOnTheFly = new TH1F("fHistoLcOnTheFly", "fHistoLcOnTheFly", 4, -0.5, 3.5);
631  TString labelOnTheFly[4] = {"OnTheFly-Bkg", "OfflineBkg", "OnTheFly-Sgn", "OfflineSgn"};
632  for (Int_t ibin = 1; ibin <= fHistoLcOnTheFly->GetNbinsX(); ibin++){
633  fHistoLcOnTheFly->GetXaxis()->SetBinLabel(ibin, labelOnTheFly[ibin-1].Data());
634  }
635 
636  fHistoLcBeforeCuts = new TH1F("fHistoLcBeforeCuts", "fHistoLcBeforeCuts", 2, -0.5, 1.5);
637  TString labelBeforeCuts[2] = {"Bkg", "Sgn"};
638  for (Int_t ibin = 1; ibin <= fHistoLcBeforeCuts->GetNbinsX(); ibin++){
639  fHistoLcBeforeCuts->GetXaxis()->SetBinLabel(ibin, labelBeforeCuts[ibin-1].Data());
640  fHistoLc->GetXaxis()->SetBinLabel(ibin, labelBeforeCuts[ibin-1].Data());
641  }
642 
643  fHistoFiducialAcceptance = new TH1F("fHistoFiducialAcceptance", "fHistoFiducialAcceptance", 4, -0.5, 3.5);
644  TString labelAcc[4] = {"NotOk-Bkg", "Ok-Bkg", "NotOk-Sgn", "Ok-Sgn"};
645  for (Int_t ibin = 1; ibin <= fHistoFiducialAcceptance->GetNbinsX(); ibin++){
646  fHistoFiducialAcceptance->GetXaxis()->SetBinLabel(ibin, labelAcc[ibin-1].Data());
647  }
648 
649  fHistoCodesSgn = new TH2F("fHistoCodesSgn", "fHistoCodes for Signal; bachelor; K0S", 7, -1.5, 5.5, 9, -1.5, 7.5);
650  fHistoCodesBkg = new TH2F("fHistoCodesBkg", "fHistoCodes for Background; bachelor; K0S", 7, -1.5, 5.5, 9, -1.5, 7.5);
651 
652  TString labelx[7] = { "kBachInvalid", "kBachFake", "kBachNoProton", "kBachPrimary", "kBachNoLambdaMother",
653  "kBachDifferentLambdaMother", "kBachCorrectLambdaMother"};
654  TString labely[9] = { "kK0SInvalid", "kK0SFake", "kK0SNoK0S", "kK0SWithoutMother", "kK0SNotFromK0",
655  "kK0Primary", "kK0NoLambdaMother", "kK0DifferentLambdaMother", "kK0CorrectLambdaMother"};
656 
657  for (Int_t ibin = 1; ibin <= fHistoCodesSgn->GetNbinsX(); ibin++){
658  fHistoCodesSgn->GetXaxis()->SetBinLabel(ibin, labelx[ibin-1].Data());
659  fHistoCodesBkg->GetXaxis()->SetBinLabel(ibin, labelx[ibin-1].Data());
660  }
661  for (Int_t ibin = 1; ibin <= fHistoCodesSgn->GetNbinsY(); ibin++){
662  fHistoCodesSgn->GetYaxis()->SetBinLabel(ibin, labely[ibin-1].Data());
663  fHistoCodesBkg->GetYaxis()->SetBinLabel(ibin, labely[ibin-1].Data());
664  }
665 
666  fHistoLcpKpiBeforeCuts = new TH1F("fHistoLcpKpiBeforeCuts", "fHistoLcpKpiBeforeCuts", 2, -0.5, 1.5);
667  for (Int_t ibin = 1; ibin <= fHistoLcpKpiBeforeCuts->GetNbinsX(); ibin++){
668  fHistoLcpKpiBeforeCuts->GetXaxis()->SetBinLabel(ibin, labelBeforeCuts[ibin-1].Data());
669  }
670 
671  fHistoBackground = new TH1F("fHistoBackground", "fHistoBackground", 4, -0.5, 3.5);
672  TString labelBkg[4] = {"Injected", "Non-injected", "Non-PYTHIA", "PYTHIA"};
673  for (Int_t ibin = 1; ibin <= fHistoBackground->GetNbinsX(); ibin++){
674  fHistoBackground->GetXaxis()->SetBinLabel(ibin, labelBkg[ibin-1].Data());
675  }
676 
677  const Float_t ptbins[15] = {0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 12., 17., 25., 35.};
678 
679  fHistoMCLcK0SpGen = new TH1F("fHistoMCLcK0SpGen", "fHistoMCLcK0SpGen", 14, ptbins);
680  fHistoMCLcK0SpGenAcc = new TH1F("fHistoMCLcK0SpGenAcc", "fHistoMCLcK0SpGenAcc", 14, ptbins);
681  fHistoMCLcK0SpGenLimAcc = new TH1F("fHistoMCLcK0SpGenLimAcc", "fHistoMCLcK0SpGenLimAcc", 14, ptbins);
682 
683 
684  // //code to run the BDT Application on-the-fly
685  // TString inputVariablesBDT = "massK0S,tImpParBach,tImpParV0,bachelorPt,combinedProtonProb,DecayLengthK0S*0.497/v0P,cosPAK0S,CosThetaStar,signd0";
686  // TObjArray *tokens = inputVariablesBDT.Tokenize(",");
687  // tokens->Print();
688  // std::vector<std::string> inputNamesVec;
689  // for(Int_t i=0; i<tokens->GetEntries(); i++){
690  // TString variable = ((TObjString*)(tokens->At(i)))->String();
691  // string tmpvar = variable.Data();
692  // inputNamesVec.push_back(tmpvar);
693  // }
694  // Printf("************************************************ fBDTReader = %p", fBDTReader);
695  // //fBDTReader = new ReadBDT_Default(inputNamesVec);
696 
697 
698  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);
699  if (fDebugHistograms) {
700  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);
701  fBDTHistoVstImpParBach = new TH2D("fBDTHistoVstImpParBach", "d0 bachelor vs bdt output; bdt; d_{0, bachelor}[cm]", 1000, -1, 1, 100, -1, 1);
702  fBDTHistoVstImpParV0 = new TH2D("fBDTHistoVstImpParV0", "d0 K0S vs bdt output; bdt; d_{0, V0}[cm]", 1000, -1, 1, 100, -1, 1);
703  fBDTHistoVsBachelorPt = new TH2D("fBDTHistoVsBachelorPt", "bachelor pT vs bdt output; bdt; p_{T, bachelor}[GeV/#it{c}]", 1000, -1, 1, 100, 0, 20);
704  fBDTHistoVsCombinedProtonProb = new TH2D("fBDTHistoVsCombinedProtonProb", "combined proton probability vs bdt output; bdt; Bayesian PID_{bachelor}", 1000, -1, 1, 100, 0, 1);
705  fBDTHistoVsCtau = new TH2D("fBDTHistoVsCtau", "K0S ctau vs bdt output; bdt; c#tau_{V0}[cm]", 1000, -1, 1, 100, -2, 2);
706  fBDTHistoVsCosPAK0S = new TH2D("fBDTHistoVsCosPAK0S", "V0 cosine pointing angle vs bdt output; bdt; CosPAK^{0}_{S}", 1000, -1, 1, 100, 0.9, 1);
707  fBDTHistoVsCosThetaStar = new TH2D("fBDTHistoVsCosThetaStar", "proton emission angle in pK0s pair rest frame vs bdt output; bdt; Cos#Theta*", 1000, -1, 1, 100, -1, 1);
708  fBDTHistoVsSignd0 = new TH2D("fBDTHistoVsSignd0", "signed d0 bachelor vs bdt output; bdt; signd_{0, bachelor}[cm]", 1000, -1, 1, 100, -1, 1);
709  fHistoNsigmaTPC = new TH2D("fHistoNsigmaTPC", "; #it{p} (GeV/#it{c}); n_{#sigma}^{TPC} (proton hypothesis)", 500, 0, 5, 1000, -5, 5);
710  fHistoNsigmaTOF = new TH2D("fHistoNsigmaTOF", "; #it{p} (GeV/#it{c}); n_{#sigma}^{TOF} (proton hypothesis)", 500, 0, 5, 1000, -5, 5);
711  }
712 
713  fOutput->Add(fHistoEvents);
718  fOutput->Add(fHistoLc);
722  fOutput->Add(fHistoCodesSgn);
723  fOutput->Add(fHistoCodesBkg);
730  fOutput->Add(fBDTHisto);
731  if (fDebugHistograms) {
737  fOutput->Add(fBDTHistoVsCtau);
741  fOutput->Add(fHistoNsigmaTPC);
742  fOutput->Add(fHistoNsigmaTOF);
743  }
744 
745  PostData(1, fOutput);
746  PostData(4, fVariablesTreeSgn);
747  PostData(5, fVariablesTreeBkg);
748  AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
749  AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
750  fPIDResponse = inputHandler->GetPIDResponse();
751 
752  // fAnalCuts->SetUsePID(kTRUE); // this forces the PID to be used, despite the settings in the cut object
753  // if (fAnalCuts->GetIsUsePID()){
754  // /*
755  // fAnalCuts->GetPidHF()->SetPidResponse(fPIDResponse);
756  // fAnalCuts->GetPidV0pos()->SetPidResponse(fPIDResponse);
757  // fAnalCuts->GetPidV0neg()->SetPidResponse(fPIDResponse);
758  // fAnalCuts->GetPidHF()->SetOldPid(kFALSE);
759  // fAnalCuts->GetPidV0pos()->SetOldPid(kFALSE);
760  // fAnalCuts->GetPidV0neg()->SetOldPid(kFALSE);
761  // */
762  // fAnalCuts->SetUsePID(kFALSE); // I don't want to use the PID through the cut object, but I will use the PID response directly!!!
763  // }
764 
765  // Setting properties of PID
766  fPIDCombined = new AliPIDCombined;
767  fPIDCombined->SetDefaultTPCPriors();
768  fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC+AliPIDResponse::kDetTOF);
769  //fPIDCombined->SetPriorDistribution((AliPID::EParticleType)ispec,fPriors[ispec]);
770 
771  fCounter = new AliNormalizationCounter("NormalizationCounter");
772  fCounter->Init();
773 
774  fListCounters = new TList();
775  fListCounters->SetOwner();
776  fListCounters->SetName("ListCounters");
777  fListCounters->Add(fCounter);
778 
779  PostData(2, fListCounters);
780 
781 
782  // Histograms from KF
783 
784  if (fCallKFVertexing){
785  fHistoDistanceLcToPrimVtx = new TH1D("fHistoDistanceLcToPrimVtx", "Lc distance to Prim Vertex from KF; distance [cm]", 1000, 0., 1);
786  fHistoDistanceV0ToPrimVtx = new TH1D("fHistoDistanceV0ToPrimVtx", "V0 distance to Prim Vertex from KF; distance [cm]", 1000, 0., 100.);
787  fHistoDistanceV0ToLc = new TH1D("fHistoDistanceV0ToLc", "V0 distance to Lc Vertex from KF; distance [cm]", 1000, 0., 100.);
788 
789  fHistoDistanceLcToPrimVtxSgn = new TH1D("fHistoDistanceLcToPrimVtxSgn", "Lc Sgn distance to Prim Vertex from KF; distance [cm]", 1000, 0., 1);
790  fHistoDistanceV0ToPrimVtxSgn = new TH1D("fHistoDistanceV0ToPrimVtxSgn", "V0 Sgn distance to Prim Vertex from KF; distance [cm]", 1000, 0., 100.);
791  fHistoDistanceV0ToLcSgn = new TH1D("fHistoDistanceV0ToLcSgn", "V0 Sgn distance to Lc Vertex from KF; distance [cm]", 1000, 0., 100.);
792 
793  fHistoVtxLcResidualToPrimVtx = new TH1D("fHistoVtxLcResidualToPrimVtx", "Residual between MC and KF (MC - KF): Lc to Prim Vtx; distance [cm]", 1000, -5., 5.);
794  fHistoVtxV0ResidualToPrimVtx = new TH1D("fHistoVtxV0ResidualToPrimVtx", "Residual between MC and KF (MC - KF): V0 to Prim Vtx; distance [cm]", 1000, -5., 5.);
795  fHistoVtxV0ResidualToLc = new TH1D("fHistoVtxV0ResidualToLc", "Residual between MC and KF: V0 to Lc (MC - KF); distance [cm]", 1000, -5., 5.);
796 
797  fHistoMassV0All = new TH1D("fHistoMassV0All", "V0 Mass; mass", 500, 0.4, 0.6);
798  fHistoDecayLengthV0All = new TH1D("fHistoDecayLengthV0All", "V0 Decay Length; decayLength", 500, -10, 10.0);
799  fHistoLifeTimeV0All = new TH1D("fHistoLifeTimeV0All", "V0 Life Time; lifeTime", 500, -10.0, 10.0);
800 
801  fHistoMassV0True = new TH1D("fHistoMassV0True", "True V0 Mass; mass", 500, 0.4, 0.6);
802  fHistoDecayLengthV0True = new TH1D("fHistoDecayLengthV0True", "True V0 Decay Length; decayLength", 500, -10, 10.0);
803  fHistoLifeTimeV0True = new TH1D("fHistoLifeTimeV0True", "True V0 Life Time; lifeTime", 500, -10.0, 10.0);
804 
805  fHistoMassV0TrueFromAOD = new TH1D("fHistoMassV0TrueFormAOD", "True V0 Mass (AOD); mass", 500, 0.4, 0.6);
806 
807  fHistoMassV0TrueK0S = new TH1D("fHistoMassV0TrueK0S", "True V0-K0S Mass; mass", 500, 0.4, 0.6);
808  fHistoDecayLengthV0TrueK0S = new TH1D("fHistoDecayLengthV0TrueK0S", "True V0-K0S Decay Length; decayLength", 500, -10, 10.0);
809  fHistoLifeTimeV0TrueK0S = new TH1D("fHistoLifeTimeV0TrueK0S", "True V0-K0S Life Time; lifeTime", 500, -10.0, 10.0);
810 
811  fHistoMassV0TrueK0SFromAOD = new TH1D("fHistoMassV0TrueK0SFormAOD", "True V0-K0S Mass (AOD); mass", 500, 0.4, 0.6);
812 
813  fHistoMassLcAll = new TH1D("fHistoMassLcAll", "Lc Mass; mass", 500, 2.0, 3.0);
814  fHistoDecayLengthLcAll = new TH1D("fHistoDecayLengthLcAll", "Lc Decay Length; decayLenght", 100000, -0.1, 0.1);
815  fHistoLifeTimeLcAll = new TH1D("fHistoLifeTimeLcAll", "Lc Life Time; lifeTime", 100000, -0.1, 0.1);
816 
817  fHistoMassLcTrue = new TH1D("fHistoMassLcTrue", "True Lc Mass; mass", 500, 2.0, 3.0);
818  fHistoDecayLengthLcTrue = new TH1D("fHistoDecayLengthLcTrue", "True Lc Decay Length; decayLength", 100000, -0.1, 0.1);
819  fHistoLifeTimeLcTrue = new TH1D("fHistoLifeTimeLcTrue", "True Lc Life Time; lifeTime", 100000, -0.1, 0.1);
820 
821  fHistoMassLcTrueFromAOD = new TH1D("fHistoMassLcTrueFromAOD", "True Lc Mass (AOD); mass", 500, 2.0, 3.0);
822 
823  fHistoMassV0fromLcAll = new TH1D("fHistoMassV0fromLcAll", "V0 mass from Lc built in KF; mass", 500, 0.4, 0.6);
824  fHistoDecayLengthV0fromLcAll = new TH1D("fHistoDecayLengthV0fromLcAll", "V0 Decay Length from Lc built in KF; decayLength", 500, 0, 10.0);
825  fHistoLifeTimeV0fromLcAll = new TH1D("fHistoLifeTimeV0fromLcAll", "V0 Life Time from Lc built in KF; lifeTime", 500, 0.0, 3.0);
826 
827  fHistoMassV0fromLcTrue = new TH1D("fHistoMassV0fromLcTrue", "V0 mass from true Lc built in KF; mass", 500, 0.4, 0.6);
828  fHistoDecayLengthV0fromLcTrue= new TH1D("fHistoDecayLengthV0fromLcTrue", "V0 Decay Length from true Lc built in KF; decayLength", 500, 0, 10.0);
829  fHistoLifeTimeV0fromLcTrue = new TH1D("fHistoLifeTimeV0fromLcTrue", "V0 Life Time from true Lc built in KF; lifeTime", 500, 0.0, 3.0);
830 
831  fHistoMassLcSgn = new TH1D("fHistoMassLcSgn", "True Lc Signal Mass; mass", 500, 2.0, 3.0);
832  fHistoMassLcSgnFromAOD = new TH1D("fHistoMassLcSgnFromAOD", "True Lc Signal Mass (AOD); mass", 500, 2.0, 3.0);
833  fHistoDecayLengthLcSgn = new TH1D("fHistoDecayLengthLcSgn", "True Lc Signal Decay Length; decayLength", 100000, -0.1, 0.1);
834  fHistoLifeTimeLcSgn = new TH1D("fHistoLifeTimeLcSgn", "True Lc Signal Life Time; lifeTime", 100000, -0.1, 0.1);
835 
836  fHistoMassV0fromLcSgn = new TH1D("fHistoMassV0fromLcSgn", "V0 from True Lc Signal Mass; mass", 500, 0.4, 0.6);
837  fHistoDecayLengthV0fromLcSgn = new TH1D("fHistoDecayLengthV0fromLcSgn", "V0 True Lc Signal Decay Length; decayLength", 500, 0, 10.0);
838  fHistoLifeTimeV0fromLcSgn = new TH1D("fHistoLifeTimeV0fromLcSgn", "V0 True Lc Signal Life Time; lifeTime", 500, 0.0, 3.0);
839 
840  fHistoKF = new TH2D("fHistoKF", "Summary from KF; V0 KF; Lc KF", 16, -0.5, 15.5, 16, -0.5, 15.5);
841  fHistoKFV0 = new TH1D("fHistoKFV0", "Summary from KF; V0 KF", 16, -0.5, 15.5);
842  fHistoKFLc = new TH1D("fHistoKFLc", "Summary from KF; V0 KF", 16, -0.5, 15.5);
843  TString axisLabel[16] = {"AllOk", "M_NotOk", "Sm_NotOk", "Dl_NotOk",
844  "Lt_NotOk", "M_Sm_NotOk", "M_Dl_NotOk", "M_Lt_NotOk",
845  "Dl_Sm_NotOk", "Dl_Lt_NotOk", "Sm_Lt_NotOk", "M_Sm_Dl_NotOk",
846  "M_Sm_Lt_NotOk", "Sm_Dl_Lt_NotOk", "M_Dl_Lt_NotOk", "All_NotOk"};
847 
848  for (Int_t ibin = 1; ibin <= 16; ibin++){
849  fHistoKF->GetXaxis()->SetBinLabel(ibin, axisLabel[ibin-1].Data());
850  fHistoKF->GetYaxis()->SetBinLabel(ibin, axisLabel[ibin-1].Data());
851  fHistoKFV0->GetXaxis()->SetBinLabel(ibin, axisLabel[ibin-1].Data());
852  fHistoKFLc->GetXaxis()->SetBinLabel(ibin, axisLabel[ibin-1].Data());
853  }
854 
855  fHistoMassKFV0 = new TH2D("fHistoMassKFV0", "mass vs sigmaMass for V0; mass; sigmaMass", 500, 0.4, 0.6, 500, 0., 10);
856  fHistoDecayLengthKFV0 = new TH2D("fHistoDecayLengthKFV0", "decayLength vs sigmaDecayLength for V0; decayLength; sigmaDecayLength", 500, -10, 10, 500, 0., 10);
857  fHistoLifeTimeKFV0 = new TH2D("fHistoLifeTimeKFV0", "lifeTime vs sigmalifeTime for V0; lifeTime; sigmaLifeTime", 500, -10, 10, 500, 0., 10);
858 
859  fHistoMassKFLc = new TH2D("fHistoMassKFLc", "mass vs sigmaMass for Lc; mass; sigmaMass", 500, 0.4, 0.6, 500, 0., 10);
860  fHistoDecayLengthKFLc = new TH2D("fHistoDecayLengthKFLc", "decayLength vs sigmaDecayLength for Lc; decayLength; sigmaDecayLength", 500, -10, 10, 500, 0., 10);
861  fHistoLifeTimeKFLc = new TH2D("fHistoLifeTimeKFLc", "lifeTime vs sigmalifeTime for Lc; lifeTime; sigmaLifeTime", 500, -10, 10, 500, 0., 10);
862 
863  fHistoArmenterosPodolanskiV0KF = new TH2D("fHistoArmenterosPodolanskiV0KF", "V0 ArmenterosPodolanski from KF; #alpha; Qt", 1000, -1, 1, 1000, 0, 1);
864  fHistoArmenterosPodolanskiV0KFSgn = new TH2D("fHistoArmenterosPodolanskiV0KFSgn", "V0 (signal) ArmenterosPodolanski from KF; #alpha; Qt", 1000, -1, 1, 1000, 0, 1);
865  fHistoArmenterosPodolanskiV0AOD = new TH2D("fHistoArmenterosPodolanskiV0AOD", "V0 ArmenterosPodolanski from AOD; #alpha; Qt", 1000, -1, 1, 1000, 0, 1);
866  fHistoArmenterosPodolanskiV0AODSgn = new TH2D("fHistoArmenterosPodolanskiV0AODSgn", "V0 (signal) ArmenterosPodolanski from AOD; #alpha; Qt", 1000, -1, 1, 1000, 0, 1);
867  }
868 
869  fOutputKF = new TList();
870  fOutputKF->SetOwner();
871  fOutputKF->SetName("listHistoKF");
872 
873  if (fCallKFVertexing){
914  fOutputKF->Add(fHistoKF);
915  fOutputKF->Add(fHistoKFV0);
916  fOutputKF->Add(fHistoKFLc);
927  }
928 
929  // weight function from ratio of flat value (1/30) to pythia
930  // use to normalise to flat distribution (should lead to flat pT distribution)
931  fFuncWeightPythia = new TF1("funcWeightPythia","1./(30. *[0]*x/TMath::Power(1.+(TMath::Power((x/[1]),[3])),[2]))", 0.15, 30);
932  fFuncWeightPythia->SetParameters(0.36609, 1.94635, 1.40463,2.5);
933 
934  // weight function from the ratio of the LHC13d3 MC
935  // and FONLL calculations for pp data
936  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.);
937  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);
938 
939  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.);
940  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);
941 
942  PostData(6, fOutputKF);
943 
944  PostData(7, fListWeight);
945 
946  // creating the BDT reader
947  if(!fFillTree){
948  std::vector<std::string> inputNamesVec;
949  TObjArray *tokens = fNamesTMVAVar.Tokenize(",");
950  for(Int_t i = 0; i < tokens->GetEntries(); i++){
951  TString variable = ((TObjString*)(tokens->At(i)))->String();
952  std::string tmpvar = variable.Data();
953  inputNamesVec.push_back(tmpvar);
954  }
955  void* lib = dlopen(fTMVAlibName.Data(), RTLD_NOW);
956  void* p = dlsym(lib, Form("ReadBDT_Default_maker%s", fTMVAlibPtBin.Data()));
957  IClassifierReader* (*maker1)(std::vector<std::string>&) = (IClassifierReader* (*)(std::vector<std::string>&)) p;
958  fBDTReader = maker1(inputNamesVec);
959  }
960  return;
961 }
962 
963 //_________________________________________________
965 {
967  if (!fInputEvent) {
968  AliError("NO EVENT FOUND!");
969  return;
970  }
971 
972  fCurrentEvent++;
973  AliDebug(2, Form("Processing event = %d", fCurrentEvent));
974  AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
975  TClonesArray *arrayLctopKos=0;
976 
977  TClonesArray *array3Prong = 0;
978 
979  if (!aodEvent && AODEvent() && IsStandardAOD()) {
980  // In case there is an AOD handler writing a standard AOD, use the AOD
981  // event in memory rather than the input (ESD) event.
982  aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
983  // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
984  // have to taken from the AOD event hold by the AliAODExtension
985  AliAODHandler* aodHandler = (AliAODHandler*)
986  ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
987 
988  if (aodHandler->GetExtensions()) {
989  AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
990  AliAODEvent *aodFromExt = ext->GetAOD();
991  arrayLctopKos = (TClonesArray*)aodFromExt->GetList()->FindObject("CascadesHF");
992 
993  array3Prong = (TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
994  }
995  } else {
996  arrayLctopKos = (TClonesArray*)aodEvent->GetList()->FindObject("CascadesHF");
997 
998  array3Prong = (TClonesArray*)aodEvent->GetList()->FindObject("Charm3Prong");
999  }
1000 
1001  if (!fUseMCInfo) {
1004  }
1005 
1006  Int_t runnumber = aodEvent->GetRunNumber();
1007  if (aodEvent->GetTriggerMask() == 0 && (runnumber >= 195344 && runnumber <= 195677)){
1008  AliDebug(3, "Event rejected because of null trigger mask");
1009  return;
1010  }
1011 
1012  fCounter->StoreEvent(aodEvent, fAnalCuts, fUseMCInfo);
1013 
1014  // mc analysis
1015  TClonesArray *mcArray = 0;
1016  AliAODMCHeader *mcHeader = 0;
1017 
1018  // multiplicity definition with tracklets
1019  fNTracklets_1 = static_cast<Int_t>(AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aodEvent, -1., 1.));
1020  fNTracklets_All = static_cast<Int_t>(AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aodEvent, -999., 999.));
1021  if (fUseMCInfo) {
1022  // MC array need for matching
1023  mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
1024  if (!mcArray) {
1025  AliError("Could not find Monte-Carlo in AOD");
1026  return;
1027  }
1028  // load MC header
1029  mcHeader = (AliAODMCHeader*)aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName());
1030  if (!mcHeader) {
1031  AliError("AliAnalysisTaskSELc2V0bachelorTMVAApp::UserExec: MC header branch not found!\n");
1032  return;
1033  }
1034 
1035  Double_t zMCVertex = mcHeader->GetVtxZ();
1036  if (TMath::Abs(zMCVertex) > fAnalCuts->GetMaxVtxZ()){
1037  AliDebug(3, Form("z coordinate of MC vertex = %f, it was required to be within [-%f, +%f], skipping event", zMCVertex, fAnalCuts->GetMaxVtxZ(), fAnalCuts->GetMaxVtxZ()));
1038  return;
1039  }
1040 
1041  FillMCHisto(mcArray);
1042 
1043  }
1044 
1045  //centrality
1047 
1048  // AOD primary vertex
1049  fVtx1 = (AliAODVertex*)aodEvent->GetPrimaryVertex();
1050  if (!fVtx1) {
1051  AliDebug(2, "No primary vertex found, returning");
1052  return;
1053  }
1054  if (fVtx1->GetNContributors()<1) {
1055  AliDebug(2, "Number of contributors to vertex < 1, returning");
1056  return;
1057  }
1058 
1060 
1061  if ( !fIsEventSelected ) {
1062  fHistoEvents->Fill(0);
1063  return; // don't take into account not selected events
1064  }
1065  fHistoEvents->Fill(1);
1066 
1071 
1073 
1074  //Setting magnetic field for KF vertexing
1075  fBField = aodEvent->GetMagneticField();
1076  AliKFParticle::SetField(fBField);
1077 
1078  Int_t nSelectedAnal = 0;
1079  if (fIsK0sAnalysis) {
1080  MakeAnalysisForLc2prK0S(aodEvent, arrayLctopKos, mcArray,
1081  nSelectedAnal, fAnalCuts,
1082  array3Prong, mcHeader);
1083  }
1084  fCounter->StoreCandidates(aodEvent, nSelectedAnal, kFALSE);
1085 
1086 
1087  //Method to get tracklet multiplicity from event
1088 
1089  // Double_t countTreta1corr = fNTracklets_1;
1090 
1091  // //get corrected tracklet multiplicity
1092  // TProfile* estimatorAvg = GetEstimatorHistogram(aodEvent);
1093  // if(estimatorAvg) {
1094  // countTreta1corr = static_cast<Int_t>(AliVertexingHFUtils::GetCorrectedNtracklets(estimatorAvg,fNTracklets_1,fVtx1->GetZ(),fRefMult));
1095  // }
1096 
1097 
1098  Double_t nchWeight = 1.;
1099  if (fNTracklets_1 > 0) {
1100 
1101  if(!fHistoMCNch) AliInfo("Input histos to evaluate Nch weights missing");
1102  if(fHistoMCNch) nchWeight *= fHistoMCNch->GetBinContent(fHistoMCNch->FindBin(fNTracklets_1));
1103  }
1104 
1105 
1106  PostData(1, fOutput);
1107  PostData(2, fListCounters);
1108  PostData(4, fVariablesTreeSgn);
1109  PostData(5, fVariablesTreeBkg);
1110  PostData(6, fOutputKF);
1111  PostData(7, fListWeight);
1112 }
1113 //-------------------------------------------------------------------------------
1115 
1117 
1118  for (Int_t iPart = 0; iPart < mcArray->GetEntriesFast(); iPart++) {
1119  AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart));
1120  if (!mcPart){
1121  AliError("Failed casting particle from MC array!, Skipping particle");
1122  continue;
1123  }
1124  Int_t pdg = mcPart->GetPdgCode();
1125  if (TMath::Abs(pdg) != 4122){
1126  AliDebug(2, Form("MC particle %d is not a Lc: its pdg code is %d", iPart, pdg));
1127  continue;
1128  }
1129  AliDebug(2, Form("Step 0 ok: MC particle %d is a Lc: its pdg code is %d", iPart, pdg));
1130  Int_t labeldaugh0 = mcPart->GetDaughterLabel(0);
1131  Int_t labeldaugh1 = mcPart->GetDaughterLabel(1);
1132  if (labeldaugh0 <= 0 || labeldaugh1 <= 0){
1133  AliDebug(2, Form("The MC particle doesn't have correct daughters, skipping!!"));
1134  continue;
1135  }
1136  else if (labeldaugh1 - labeldaugh0 == 1){
1137  AliDebug(2, Form("Step 1 ok: The MC particle has correct daughters!!"));
1138  AliAODMCParticle* daugh0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labeldaugh0));
1139  AliAODMCParticle* daugh1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labeldaugh1));
1140  if(!daugh0 || !daugh1){
1141  AliDebug(2,"Particle daughters not properly retrieved!");
1142  return;
1143  }
1144  Int_t pdgCodeDaugh0 = TMath::Abs(daugh0->GetPdgCode());
1145  Int_t pdgCodeDaugh1 = TMath::Abs(daugh1->GetPdgCode());
1146  AliAODMCParticle* bachelorMC = daugh0;
1147  AliAODMCParticle* v0MC = daugh1;
1148  AliDebug(2, Form("pdgCodeDaugh0 = %d, pdgCodeDaugh1 = %d", pdgCodeDaugh0, pdgCodeDaugh1));
1149  if ((pdgCodeDaugh0 == 311 && pdgCodeDaugh1 == 2212) || (pdgCodeDaugh0 == 2212 && pdgCodeDaugh1 == 311)){
1150  // 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-
1152  if (pdgCodeDaugh0 == 311 && pdgCodeDaugh1 == 2212) {
1153  bachelorMC = daugh1;
1154  v0MC = daugh0;
1155  }
1156  AliDebug(2, Form("Number of Daughters of v0 = %d", v0MC->GetNDaughters()));
1157  if (v0MC->GetNDaughters() != 1) {
1158  AliDebug(2, "The K0 does not decay in 1 body only! Impossible... Continuing...");
1159  continue;
1160  }
1161  else { // So far: Lc --> K0 + p, K0 with 1 daughter
1162  AliDebug(2, "Step 2 ok: The K0 does decay in 1 body only! ");
1163  Int_t labelK0daugh = v0MC->GetDaughterLabel(0);
1164  AliAODMCParticle* partK0S = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelK0daugh));
1165  if(!partK0S){
1166  AliError("Error while casting particle! returning a NULL array");
1167  continue;
1168  }
1169  else { // So far: Lc --> K0 + p, K0 with 1 daughter that we can access
1170  if (partK0S->GetNDaughters() != 2 || TMath::Abs(partK0S->GetPdgCode() != 310)){
1171  AliDebug(2, "The K0 daughter is not a K0S or does not decay in 2 bodies");
1172  continue;
1173  }
1174  else { // So far: Lc --> K0 + p, K0 --> K0S, K0S in 2 bodies
1175  AliDebug(2, "Step 3 ok: The K0 daughter is a K0S and does decay in 2 bodies");
1176  Int_t labelK0Sdaugh0 = partK0S->GetDaughterLabel(0);
1177  Int_t labelK0Sdaugh1 = partK0S->GetDaughterLabel(1);
1178  AliAODMCParticle* daughK0S0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelK0Sdaugh0));
1179  AliAODMCParticle* daughK0S1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelK0Sdaugh1));
1180  if (!daughK0S0 || ! daughK0S1){
1181  AliDebug(2, "Could not access K0S daughters, continuing...");
1182  continue;
1183  }
1184  else { // So far: Lc --> K0 + p, K0 --> K0S, K0S in 2 bodies that we can access
1185  AliDebug(2, "Step 4 ok: Could access K0S daughters, continuing...");
1186  Int_t pdgK0Sdaugh0 = daughK0S0->GetPdgCode();
1187  Int_t pdgK0Sdaugh1 = daughK0S1->GetPdgCode();
1188  if (TMath::Abs(pdgK0Sdaugh0) != 211 || TMath::Abs(pdgK0Sdaugh1) != 211){
1189  AliDebug(2, "The K0S does not decay in pi+pi-, continuing");
1190  //AliInfo("The K0S does not decay in pi+pi-, continuing");
1191  }
1192  else { // Full chain: Lc --> K0 + p, K0 --> K0S, K0S --> pi+pi-
1193  if (fAnalCuts->IsInFiducialAcceptance(mcPart->Pt(), mcPart->Y())) {
1194  AliDebug(2, Form("----> Filling histo with pt = %f", mcPart->Pt()));
1195  if(TMath::Abs(mcPart->Y()) < 0.5) fHistoMCLcK0SpGenLimAcc->Fill(mcPart->Pt());
1196  //AliInfo(Form("\nparticle = %d, Filling MC Gen histo\n", iPart));
1197  fHistoMCLcK0SpGen->Fill(mcPart->Pt());
1198  if(!(TMath::Abs(bachelorMC->Eta()) > 0.9 || bachelorMC->Pt() < 0.1 ||
1199  TMath::Abs(daughK0S0->Eta()) > 0.9 || daughK0S0->Pt() < 0.1 ||
1200  TMath::Abs(daughK0S1->Eta()) > 0.9 || daughK0S1->Pt() < 0.1)) {
1201  fHistoMCLcK0SpGenAcc->Fill(mcPart->Pt());
1202  }
1203  }
1204  else {
1205  AliDebug(2, "not in fiducial acceptance! Skipping");
1206  continue;
1207  }
1208  }
1209  }
1210  }
1211  }
1212  }
1213  }
1214  }
1215  } // closing loop over mcArray
1216 
1217  return;
1218 
1219 }
1220 
1221 //-------------------------------------------------------------------------------
1223  TClonesArray *arrayLctopKos,
1224  TClonesArray *mcArray,
1225  Int_t &nSelectedAnal,
1226  AliRDHFCutsLctoV0 *cutsAnal, TClonesArray *array3Prong,
1227  AliAODMCHeader* aodheader){
1228 
1230 
1231  Int_t pdgCand = 4122;
1232  Int_t pdgDgLctoV0bachelor[2] = {2212, 310};
1233  Int_t pdgDgV0toDaughters[2] = {211, 211};
1234 
1235  // loop to search for candidates Lc->K0sp
1236  Int_t nCascades= arrayLctopKos->GetEntriesFast();
1237 
1238  // loop over cascades to search for candidates Lc->p+K0S
1239 
1240  Int_t mcLabel = -1;
1241 
1243  for (Int_t iLctopK0s = 0; iLctopK0s < nCascades; iLctopK0s++) {
1244 
1245  // Lc candidates and K0s from Lc
1246  AliAODRecoCascadeHF* lcK0spr = dynamic_cast<AliAODRecoCascadeHF*>(arrayLctopKos->At(iLctopK0s));
1247  if (!lcK0spr) {
1248  AliDebug(2, Form("Cascade %d doens't exist, skipping",iLctopK0s));
1249  continue;
1250  }
1251 
1252  if (!(lcK0spr->CheckCascadeFlags())) {
1253  AliDebug(2, Form("Cascade %d is not flagged as Lc candidate",iLctopK0s));
1254  continue;
1255  }
1256 
1257  if(!vHF->FillRecoCasc(aodEvent, lcK0spr, kFALSE)){ //Fill the data members of the candidate only if they are empty.
1258  //fCEvents->Fill(18);//monitor how often this fails
1259  continue;
1260  }
1261  //if (!(vHF->RecoSecondaryVertexForCascades(aodEvent, lcK0spr))) continue;
1262 
1263 
1264  if (!lcK0spr->GetSecondaryVtx()) {
1265  AliInfo("No secondary vertex");
1266  continue;
1267  }
1268 
1269  if (lcK0spr->GetNDaughters()!=2) {
1270  AliDebug(2, Form("Cascade %d has not 2 daughters (nDaughters=%d)", iLctopK0s, lcK0spr->GetNDaughters()));
1271  continue;
1272  }
1273 
1274  AliAODv0 * v0part = dynamic_cast<AliAODv0*>(lcK0spr->Getv0());
1275  AliAODTrack * bachPart = dynamic_cast<AliAODTrack*>(lcK0spr->GetBachelor());
1276  if (!v0part || !bachPart) {
1277  AliDebug(2, Form("Cascade %d has no V0 or no bachelor object", iLctopK0s));
1278  continue;
1279  }
1280 
1281  if (!v0part->GetSecondaryVtx()) {
1282  AliDebug(2, Form("No secondary vertex for V0 by cascade %d", iLctopK0s));
1283  continue;
1284  }
1285 
1286  if (v0part->GetNDaughters()!=2) {
1287  AliDebug(2, Form("current V0 has not 2 daughters (onTheFly=%d, nDaughters=%d)", v0part->GetOnFlyStatus(), v0part->GetNDaughters()));
1288  continue;
1289  }
1290 
1291  AliAODTrack * v0Pos = dynamic_cast<AliAODTrack*>(lcK0spr->Getv0PositiveTrack());
1292  AliAODTrack * v0Neg = dynamic_cast<AliAODTrack*>(lcK0spr->Getv0NegativeTrack());
1293  if (!v0Neg || !v0Pos) {
1294  AliDebug(2, Form("V0 by cascade %d has no V0positive of V0negative object", iLctopK0s));
1295  continue;
1296  }
1297 
1298  if (v0Pos->Charge() == v0Neg->Charge()) {
1299  AliDebug(2, Form("V0 by cascade %d has daughters with the same sign: IMPOSSIBLE!", iLctopK0s));
1300  continue;
1301  }
1302 
1303  //Filling a control histogram with no cuts
1304  if (fUseMCInfo) {
1305 
1306  Int_t pdgCode=-2;
1307 
1308  // find associated MC particle for Lc -> p+K0 and K0S->pi+pi
1309  fmcLabelLc = lcK0spr->MatchToMC(pdgCand, pdgDgLctoV0bachelor[1], pdgDgLctoV0bachelor, pdgDgV0toDaughters, mcArray, kTRUE);
1310  //if (fmcLabelLc>=0) {
1311  if (fmcLabelLc != -1) {
1312  AliDebug(2, Form("----> cascade number %d (total cascade number = %d) is a Lc!", iLctopK0s, nCascades));
1313 
1314  AliAODMCParticle *partLc = dynamic_cast<AliAODMCParticle*>(mcArray->At(fmcLabelLc));
1315  if(partLc){
1316  pdgCode = partLc->GetPdgCode();
1317  if (pdgCode<0) AliDebug(2,Form(" ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ MClabel=%d ~~~~~~~~~~ pdgCode=%d", fmcLabelLc, pdgCode));
1318  pdgCode = TMath::Abs(pdgCode);
1319  fHistoLcBeforeCuts->Fill(1);
1320  }
1321  }
1322  else {
1323  fHistoLcBeforeCuts->Fill(0);
1324  pdgCode = -1;
1325  }
1326 
1327  }
1328 
1329  Int_t isLc = 0;
1330 
1331  if (fUseMCInfo) {
1332 
1333  Int_t pdgCode = -2;
1334 
1335  // find associated MC particle for Lc -> p+K0 and K0S->pi+pi
1336  mcLabel = lcK0spr->MatchToMC(pdgCand, pdgDgLctoV0bachelor[1], pdgDgLctoV0bachelor, pdgDgV0toDaughters, mcArray, kTRUE);
1337  if (mcLabel >= 0) {
1338  AliDebug(2, Form(" ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~cascade number %d (total cascade number = %d)", iLctopK0s, nCascades));
1339 
1340  AliAODMCParticle *partLc = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcLabel));
1341  if(partLc){
1342  pdgCode = partLc->GetPdgCode();
1343  if (pdgCode < 0) AliDebug(2, Form(" ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ MClabel=%d ~~~~~~~~~~ pdgCode=%d", mcLabel, pdgCode));
1344  pdgCode = TMath::Abs(pdgCode);
1345  isLc = 1;
1346  fHistoLc->Fill(1);
1347  }
1348  }
1349  else {
1350  fHistoLc->Fill(0);
1351  pdgCode = -1;
1352  }
1353  }
1354  AliDebug(2, Form("\n\n\n Analysing candidate %d\n", iLctopK0s));
1355  AliDebug(2, Form(">>>>>>>>>> Candidate is background, fFillOnlySgn = %d --> SKIPPING", fFillOnlySgn));
1356  if (!isLc) {
1357  if (fFillOnlySgn) { // if it is background, and we want only signal, we do not fill the tree
1358  continue;
1359  }
1360  else { // checking if we want to fill the background
1361  if (fKeepingOnlyHIJINGBkg){
1362  // we have decided to fill the background only when the candidate has the daugthers that all come from HIJING underlying event!
1363  Bool_t isCandidateInjected = fUtils->HasCascadeCandidateAnyDaughInjected(lcK0spr, aodheader, mcArray);
1364  if (!isCandidateInjected){
1365  AliDebug(2, "The candidate is from HIJING (i.e. not injected), keeping it to fill background");
1366  fHistoBackground->Fill(1);
1367  }
1368  else {
1369  AliDebug(2, "The candidate is NOT from HIJING, we skip it when filling background");
1370  fHistoBackground->Fill(0);
1371  continue;
1372  }
1373  }
1374  else if (fKeepingOnlyPYTHIABkg){
1375  // we have decided to fill the background only when the candidate has the daugthers that all come from HIJING underlying event!
1376  AliAODTrack *bachelor = (AliAODTrack*)lcK0spr->GetBachelor();
1377  AliAODTrack *v0pos = (AliAODTrack*)lcK0spr->Getv0PositiveTrack();
1378  AliAODTrack *v0neg = (AliAODTrack*)lcK0spr->Getv0NegativeTrack();
1379  if (!bachelor || !v0pos || !v0neg) {
1380  AliDebug(2, "Cannot retrieve one of the tracks while checking origin, continuing");
1381  continue;
1382  }
1383  else {
1384  Int_t labelbachelor = TMath::Abs(bachelor->GetLabel());
1385  Int_t labelv0pos = TMath::Abs(v0pos->GetLabel());
1386  Int_t labelv0neg = TMath::Abs(v0neg->GetLabel());
1387  AliAODMCParticle* MCbachelor = (AliAODMCParticle*)mcArray->At(labelbachelor);
1388  AliAODMCParticle* MCv0pos = (AliAODMCParticle*)mcArray->At(labelv0pos);
1389  AliAODMCParticle* MCv0neg = (AliAODMCParticle*)mcArray->At(labelv0neg);
1390  if (!MCbachelor || !MCv0pos || !MCv0neg) {
1391  AliDebug(2, "Cannot retrieve MC particle for one of the tracks while checking origin, continuing");
1392  continue;
1393  }
1394  else {
1395  Int_t isBachelorFromPythia = fUtils->CheckOrigin(mcArray, MCbachelor, kTRUE);
1396  Int_t isv0posFromPythia = fUtils->CheckOrigin(mcArray, MCv0pos, kTRUE);
1397  Int_t isv0negFromPythia = fUtils->CheckOrigin(mcArray, MCv0neg, kTRUE);
1398  if (isBachelorFromPythia != 0 && isv0posFromPythia != 0 && isv0negFromPythia != 0){
1399  AliDebug(2, "The candidate is from PYTHIA (i.e. all daughters originate from a quark), keeping it to fill background");
1400  fHistoBackground->Fill(2);
1401  }
1402  else {
1403  AliDebug(2, "The candidate is NOT from PYTHIA, we skip it when filling background");
1404  fHistoBackground->Fill(3);
1405  continue;
1406  }
1407  }
1408  }
1409  }
1410  }
1411  }
1412 
1413  //FillLc2pK0Sspectrum(lcK0spr, isLc, nSelectedAnal, cutsAnal, mcArray, iLctopK0s);
1414  FillLc2pK0Sspectrum(lcK0spr, isLc, nSelectedAnal, cutsAnal, mcArray, mcLabel);
1415  }
1416 
1417  delete vHF;
1418 
1419  return;
1420 
1421 }
1422 //________________________________________________________________________
1424  Int_t isLc,
1425  Int_t &nSelectedAnal,
1426  AliRDHFCutsLctoV0 *cutsAnal,
1427  TClonesArray *mcArray, Int_t iLctopK0s){
1428  //
1430  //
1431  /*
1432  if (!part->GetOwnPrimaryVtx()) {
1433  //Printf("No primary vertex for Lc found!!");
1434  part->SetOwnPrimaryVtx(fVtx1);
1435  }
1436  else {
1437  //Printf("Yu-huuuu!!! primary vertex for Lc found!!");
1438  }
1439  */
1440  Double_t invmassLc = part->InvMassLctoK0sP();
1441 
1442  AliAODv0 * v0part = part->Getv0();
1443  Bool_t onFlyV0 = v0part->GetOnFlyStatus(); // on-the-flight V0s
1444  if (onFlyV0){ // on-the-fly V0
1445  if (isLc) { // Lc
1446  fHistoLcOnTheFly->Fill(2.);
1447  }
1448  else { // not Lc
1449  fHistoLcOnTheFly->Fill(0.);
1450  }
1451  }
1452  else { // offline V0
1453  if (isLc) { // Lc
1454  fHistoLcOnTheFly->Fill(3.);
1455  }
1456  else { // not Lc
1457  fHistoLcOnTheFly->Fill(1.);
1458  }
1459  }
1460 
1461  Double_t dcaV0 = v0part->GetDCA();
1462  Double_t invmassK0s = v0part->MassK0Short();
1463 
1464  if ( (cutsAnal->IsInFiducialAcceptance(part->Pt(), part->Y(4122))) ) {
1465  if (isLc) {
1466  fHistoFiducialAcceptance->Fill(3.);
1467  }
1468  else {
1469  fHistoFiducialAcceptance->Fill(1.);
1470  }
1471  }
1472  else {
1473  if (isLc) {
1474  fHistoFiducialAcceptance->Fill(2.);
1475  }
1476  else {
1477  fHistoFiducialAcceptance->Fill(0.);
1478  }
1479  }
1480 
1481  Int_t isInV0window = (((cutsAnal->IsSelectedSingleCut(part, AliRDHFCuts::kCandidate, 2)) & (AliRDHFCutsLctoV0::kLcToK0Spr)) == (AliRDHFCutsLctoV0::kLcToK0Spr)); // cut on V0 invMass
1482 
1483  if (isInV0window == 0) {
1484  AliDebug(2, "No: The candidate has NOT passed the V0 window cuts!");
1485  if (isLc) AliDebug(2, "SIGNAL candidate rejected: V0 window cuts");
1486  return;
1487  }
1488  else AliDebug(2, "Yes: The candidate has passed the mass cuts!");
1489 
1490  Bool_t isInCascadeWindow = (((cutsAnal->IsSelectedSingleCut(part, AliRDHFCuts::kCandidate, 0)) & (AliRDHFCutsLctoV0::kLcToK0Spr)) == (AliRDHFCutsLctoV0::kLcToK0Spr)); // cut on Lc->p+K0S invMass
1491 
1492  if (!isInCascadeWindow) {
1493  AliDebug(2, "No: The candidate has NOT passed the cascade window cuts!");
1494  if (isLc) AliDebug(2, "SIGNAL candidate rejected: cascade window cuts");
1495  return;
1496  }
1497  else AliDebug(2, "Yes: The candidate has passed the cascade window cuts!");
1498 
1499  Bool_t isCandidateSelectedCuts = (((cutsAnal->IsSelected(part, AliRDHFCuts::kCandidate)) & (AliRDHFCutsLctoV0::kLcToK0Spr)) == (AliRDHFCutsLctoV0::kLcToK0Spr)); // kinematic/topological cuts
1500  AliDebug(2, Form("recoAnalysisCuts = %d", cutsAnal->IsSelected(part, AliRDHFCuts::kCandidate) & (AliRDHFCutsLctoV0::kLcToK0Spr)));
1501  if (!isCandidateSelectedCuts){
1502  AliDebug(2, "No: Analysis cuts kCandidate level NOT passed");
1503  if (isLc) AliDebug(2, "SIGNAL candidate rejected");
1504  return;
1505  }
1506  else {
1507  AliDebug(2, "Yes: Analysis cuts kCandidate level passed");
1508  }
1509 
1510  AliAODTrack *bachelor = (AliAODTrack*)part->GetBachelor();
1511  if (!bachelor) {
1512  AliDebug(2, Form("Very weird, the bachelor is not there... returning for this candidate"));
1513  return;
1514  }
1515 
1516  //Bool_t isBachelorID = (((cutsAnal->IsSelected(part,AliRDHFCuts::kPID))&(AliRDHFCutsLctoV0::kLcToK0Spr))==(AliRDHFCutsLctoV0::kLcToK0Spr)); // ID x bachelor
1517  Double_t probTPCTOF[AliPID::kSPECIES] = {-1.};
1518 
1519  UInt_t detUsed = fPIDCombined->ComputeProbabilities(bachelor, fPIDResponse, probTPCTOF);
1520  AliDebug(2, Form("detUsed (TPCTOF case) = %d", detUsed));
1521  Double_t probProton = -1.;
1522  // Double_t probPion = -1.;
1523  // Double_t probKaon = -1.;
1524  if (detUsed == (UInt_t)fPIDCombined->GetDetectorMask() ) {
1525  AliDebug(2, Form("We have found the detector mask for TOF + TPC: probProton will be set to %f", probTPCTOF[AliPID::kProton]));
1526  probProton = probTPCTOF[AliPID::kProton];
1527  // probPion = probTPCTOF[AliPID::kPion];
1528  // probKaon = probTPCTOF[AliPID::kKaon];
1529  }
1530  else { // if you don't have both TOF and TPC, try only TPC
1531  fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
1532  AliDebug(2, "We did not find the detector mask for TOF + TPC, let's see only TPC");
1533  detUsed = fPIDCombined->ComputeProbabilities(bachelor, fPIDResponse, probTPCTOF);
1534  AliDebug(2, Form(" detUsed (TPC case) = %d", detUsed));
1535  if (detUsed == (UInt_t)fPIDCombined->GetDetectorMask()) {
1536  probProton = probTPCTOF[AliPID::kProton];
1537  // probPion = probTPCTOF[AliPID::kPion];
1538  // probKaon = probTPCTOF[AliPID::kKaon];
1539  AliDebug(2, Form("TPC only worked: probProton will be set to %f", probTPCTOF[AliPID::kProton]));
1540  }
1541  else {
1542  AliDebug(2, "Only TPC did not work...");
1543  }
1544  // resetting mask to ask for both TPC+TOF
1545  fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC+AliPIDResponse::kDetTOF);
1546  }
1547  AliDebug(2, Form("probProton = %f", probProton));
1548 
1549  // now we get the TPC and TOF single PID probabilities (only for Proton, or the tree will explode :) )
1550  Double_t probProtonTPC = -1.;
1551  Double_t probProtonTOF = -1.;
1552  Double_t pidTPC[AliPID::kSPECIES]={-1.};
1553  Double_t pidTOF[AliPID::kSPECIES]={-1.};
1554  Int_t respTPC = fPIDResponse->ComputePIDProbability(AliPIDResponse::kDetTPC, bachelor, AliPID::kSPECIES, pidTPC);
1555  Int_t respTOF = fPIDResponse->ComputePIDProbability(AliPIDResponse::kDetTOF, bachelor, AliPID::kSPECIES, pidTOF);
1556  if (respTPC == AliPIDResponse::kDetPidOk) probProtonTPC = pidTPC[AliPID::kProton];
1557  if (respTOF == AliPIDResponse::kDetPidOk) probProtonTOF = pidTOF[AliPID::kProton];
1558 
1559  // checking V0 status (on-the-fly vs offline)
1560  if ( !( !onFlyV0 || (onFlyV0 && fUseOnTheFlyV0) ) ) {
1561  AliDebug(2, "On-the-fly discarded");
1562  return;
1563  }
1564 
1566  nSelectedAnal++;
1567  }
1568 
1569  if ( !(cutsAnal->IsInFiducialAcceptance(part->Pt(), part->Y(4122))) ) return;
1570 
1571  if ( !( ( (cutsAnal->IsSelected(part, AliRDHFCuts::kTracks)) & (AliRDHFCutsLctoV0::kLcToK0Spr)) == (AliRDHFCutsLctoV0::kLcToK0Spr) ) ) { // esd track cuts
1572  if (isLc) AliDebug(2, "SIGNAL candidate rejected");
1573  AliDebug(2, "No: Analysis cuts kTracks level NOT passed");
1574  return;
1575  }
1576  else {
1577  AliDebug(2, "Yes: Analysis cuts kTracks level passed");
1578  }
1579 
1580  Int_t pdgCand = 4122;
1581  Int_t pdgDgLctoV0bachelor[2] = {211, 3122}; // case of wrong decay! Lc --> L + pi
1582  Int_t pdgDgV0toDaughters[2] = {2212, 211}; // case of wrong decay! Lc --> L + pi
1583  Int_t isLc2LBarpi = 0, isLc2Lpi = 0;
1584  Int_t currentLabel = part->GetLabel();
1585  Int_t mcLabel = 0;
1586  if (fUseMCInfo) {
1587  mcLabel = part->MatchToMC(pdgCand, pdgDgLctoV0bachelor[1], pdgDgLctoV0bachelor, pdgDgV0toDaughters, mcArray, kTRUE);
1588  if (mcLabel >= 0) {
1589  if (bachelor->Charge() == -1) isLc2LBarpi=1;
1590  if (bachelor->Charge() == +1) isLc2Lpi=1;
1591  }
1592  }
1593 
1594  Int_t pdgDg2prong[2] = {211, 211};
1595  Int_t labelK0S = 0;
1596  Int_t isK0S = 0;
1597  if (fUseMCInfo) {
1598  labelK0S = v0part->MatchToMC(310, mcArray, 2, pdgDg2prong);
1599  if (labelK0S>=0) isK0S = 1;
1600  }
1601 
1602  pdgDg2prong[0] = 211;
1603  pdgDg2prong[1] = 2212;
1604  Int_t isLambda = 0;
1605  Int_t isLambdaBar = 0;
1606  Int_t lambdaLabel = 0;
1607  if (fUseMCInfo) {
1608  lambdaLabel = v0part->MatchToMC(3122, mcArray, 2, pdgDg2prong);
1609  if (lambdaLabel >= 0) {
1610  AliAODMCParticle *lambdaTrack = (AliAODMCParticle*)mcArray->At(lambdaLabel);
1611  if (lambdaTrack->GetPdgCode() == 3122) isLambda = 1;
1612  else if (lambdaTrack->GetPdgCode() == -3122) isLambdaBar = 1;
1613  }
1614  }
1615 
1616  pdgDg2prong[0] = 11;
1617  pdgDg2prong[1] = 11;
1618  Int_t isGamma = 0;
1619  Int_t gammaLabel = 0;
1620  if (fUseMCInfo) {
1621  gammaLabel = v0part->MatchToMC(22, mcArray, 2, pdgDg2prong);
1622  if (gammaLabel>=0) {
1623  AliAODMCParticle *gammaTrack = (AliAODMCParticle*)mcArray->At(gammaLabel);
1624  if (gammaTrack->GetPdgCode() == 22) isGamma = 1;
1625  }
1626  }
1627 
1628  Int_t pdgTemp = -1;
1629  if (currentLabel != -1){
1630  AliAODMCParticle *tempPart = (AliAODMCParticle*)mcArray->At(currentLabel);
1631  pdgTemp = tempPart->GetPdgCode();
1632  }
1633  if (isLc) AliDebug(2, Form("Signal: Candidate is a Lc in K0s+p"));
1634  else if (isLc2LBarpi) AliDebug(2, Form("Background: Candidate is a Lc in Lbar + pi"));
1635  else if (isLc2Lpi) AliDebug(2, Form("Background: Candidate is a Lc in L + pi"));
1636  else AliDebug(2, Form("Pure bkg: Candidate has pdg = %d", pdgTemp));
1637  if (isK0S) AliDebug(2, Form("V0 is a K0S"));
1638  else if (isLambda) AliDebug(2, Form("V0 is a Lambda"));
1639  else if (isLambdaBar) AliDebug(2, Form("V0 is a LambdaBar"));
1640  else if (isGamma) AliDebug(2, Form("V0 is a Gamma"));
1641  //else AliDebug(2, Form("V0 is something else!!"));
1642 
1643  Double_t invmassLc2Lpi = part->InvMassLctoLambdaPi();
1644  Double_t invmassLambda = v0part->MassLambda();
1645  Double_t invmassLambdaBar = v0part->MassAntiLambda();
1646 
1647  //Double_t nSigmaITSpr = -999.;
1648  Double_t nSigmaTPCpr = -999.;
1649  Double_t nSigmaTOFpr = -999.;
1650 
1651  //Double_t nSigmaITSpi = -999.;
1652  Double_t nSigmaTPCpi = -999.;
1653  Double_t nSigmaTOFpi = -999.;
1654 
1655  //Double_t nSigmaITSka = -999.;
1656  Double_t nSigmaTPCka = -999.;
1657  Double_t nSigmaTOFka = -999.;
1658 
1659  /*
1660  cutsAnal->GetPidHF()->GetnSigmaITS(bachelor, 4, nSigmaITSpr);
1661  cutsAnal->GetPidHF()->GetnSigmaTPC(bachelor, 4, nSigmaTPCpr);
1662  cutsAnal->GetPidHF()->GetnSigmaTOF(bachelor, 4, nSigmaTOFpr);
1663  cutsAnal->GetPidHF()->GetnSigmaITS(bachelor, 2, nSigmaITSpi);
1664  cutsAnal->GetPidHF()->GetnSigmaTPC(bachelor, 2, nSigmaTPCpi);
1665  cutsAnal->GetPidHF()->GetnSigmaTOF(bachelor, 2, nSigmaTOFpi);
1666  cutsAnal->GetPidHF()->GetnSigmaITS(bachelor, 3, nSigmaITSka);
1667  cutsAnal->GetPidHF()->GetnSigmaTPC(bachelor, 3, nSigmaTPCka);
1668  cutsAnal->GetPidHF()->GetnSigmaTOF(bachelor, 3, nSigmaTOFka);
1669  */
1670 
1671  nSigmaTPCpi = fPIDResponse->NumberOfSigmasTPC(bachelor, (AliPID::kPion));
1672  nSigmaTPCka = fPIDResponse->NumberOfSigmasTPC(bachelor, (AliPID::kKaon));
1673  nSigmaTPCpr = fPIDResponse->NumberOfSigmasTPC(bachelor, (AliPID::kProton));
1674  nSigmaTOFpi = fPIDResponse->NumberOfSigmasTOF(bachelor, (AliPID::kPion));
1675  nSigmaTOFka = fPIDResponse->NumberOfSigmasTOF(bachelor, (AliPID::kKaon));
1676  nSigmaTOFpr = fPIDResponse->NumberOfSigmasTOF(bachelor, (AliPID::kProton));
1677 
1678  Double_t ptLcMC = -1;
1679  Double_t weightPythia = -1, weight5LHC13d3 = -1, weight5LHC13d3Lc = -1;
1680 
1681  if (fUseMCInfo) {
1682  if (iLctopK0s >= 0) {
1683  AliAODMCParticle *partLcMC = (AliAODMCParticle*)mcArray->At(iLctopK0s);
1684  ptLcMC = partLcMC->Pt();
1685  //Printf("--------------------- Reco pt = %f, MC particle pt = %f", part->Pt(), ptLcMC);
1686  weightPythia = fFuncWeightPythia->Eval(ptLcMC);
1687  weight5LHC13d3 = fFuncWeightFONLL5overLHC13d3->Eval(ptLcMC);
1688  weight5LHC13d3Lc = fFuncWeightFONLL5overLHC13d3Lc->Eval(ptLcMC);
1689  }
1690  }
1691 
1692  Double_t weightNch = 1;
1693  if (fUseMCInfo) {
1694  //Int_t nChargedMCPhysicalPrimary=AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(mcArray,-1.0,1.0);
1695  // if(nChargedMCPhysicalPrimary > 0)
1696 
1697  if(fNTracklets_1 > 0){
1698  if(!fHistoMCNch) AliDebug(2, "Input histos to evaluate Nch weights missing");
1699  if(fHistoMCNch) weightNch *= fHistoMCNch->GetBinContent(fHistoMCNch->FindBin(fNTracklets_1));
1700  }
1701  }
1702 
1703 
1704  // Fill candidate variable Tree (track selection, V0 invMass selection)
1705  // if (!onFlyV0 && isInV0window && isInCascadeWindow && part->CosV0PointingAngle()>0.99 && TMath::Abs(nSigmaTPCpr) <= 3 &&
1706  // TMath::Abs(nSigmaTOFpr) <= 3 && v0part->Getd0Prong(0) < 20 && v0part->Getd0Prong(1) < 20) {
1707  if (!onFlyV0 && isInV0window && isInCascadeWindow){
1708 
1709  EBachelor bachCode = kBachInvalid;
1710  EK0S k0SCode = kK0SInvalid;
1711  if (fUseMCInfo) {
1712  bachCode = CheckBachelor(part, bachelor, mcArray);
1713  k0SCode = CheckK0S(part, v0part, mcArray);
1714  }
1715 
1716  Double_t V0KF[3] = {-999999, -999999, -999999};
1717  Double_t errV0KF[3] = {-999999, -999999, -999999};
1718  Double_t LcKF[3] = {-999999, -999999, -999999};
1719  Double_t errLcKF[3] = {-999999, -999999, -999999};
1720  Double_t distances[3] = {-999999, -999999, -999999};
1721  Double_t armPolKF[2] = {-999999, -999999};
1722 
1723  AliAODTrack *v0pos = (AliAODTrack*)part->Getv0PositiveTrack();
1724  AliAODTrack *v0neg = (AliAODTrack*)part->Getv0NegativeTrack();
1725  Int_t kfResult;
1726  TVector3 mom1(bachelor->Px(), bachelor->Py(), bachelor->Pz());
1727  TVector3 mom2(v0part->Px(), v0part->Py(), v0part->Pz());
1728  TVector3 momTot(part->Px(), part->Py(), part->Pz());
1729 
1730  Double_t Ql1 = mom1.Dot(momTot)/momTot.Mag();
1731  Double_t Ql2 = mom2.Dot(momTot)/momTot.Mag();
1732 
1733  Double_t alphaArmLc = (Ql1 - Ql2)/(Ql1 + Ql2);
1734  Double_t alphaArmLcCharge = ( bachelor->Charge() > 0 ? (Ql1 - Ql2)/(Ql1 + Ql2) : (Ql2 - Ql1)/(Ql1 + Ql2) );
1735  Double_t ptArmLc = mom1.Perp(momTot);
1736 
1737  Double_t massK0SPDG = TDatabasePDG::Instance()->GetParticle(310)->Mass(); // mass K0S PDG
1738  Double_t massPrPDG = TDatabasePDG::Instance()->GetParticle(2212)->Mass(); // mass Proton PDG
1739  Double_t massLcPDG = TDatabasePDG::Instance()->GetParticle(4122)->Mass(); // mass Lc PDG
1740 
1741  // Cosine of proton emission angle (theta*) in the rest frame of the mother particle
1742  // for prong ip (0 or 1) with mass hypotheses massLcPDG for mother particle (from AliAODRecoDecay)
1743  Double_t pStar = TMath::Sqrt((massLcPDG*massLcPDG-massPrPDG*massPrPDG-massK0SPDG*massK0SPDG)*(massLcPDG*massLcPDG-massPrPDG*massPrPDG-massK0SPDG*massK0SPDG)-4.*massPrPDG*massPrPDG*massK0SPDG*massK0SPDG)/(2.*massLcPDG);
1744  Double_t e = part->E(4122);
1745  Double_t beta = part->P()/e;
1746  Double_t gamma = e/massLcPDG;
1747  //Double_t cts = (Ql1/gamma-beta*TMath::Sqrt(pStar*pStar+massPrPDG*massPrPDG))/pStar;
1748 
1749  // Cosine of proton emission angle (theta*) in the rest frame of the mother particle
1750  // (from AliRDHFCutsLctoV0)
1751  TLorentzVector vpr, vk0s,vlc;
1752  vpr.SetXYZM(part->PxProng(0), part->PyProng(0), part->PzProng(0), massPrPDG);
1753  vk0s.SetXYZM(part->PxProng(1), part->PyProng(1), part->PzProng(1), massK0SPDG);
1754  vlc = vpr + vk0s;
1755  TVector3 vboost = vlc.BoostVector();
1756  vpr.Boost(-vboost);
1757  Double_t cts = TMath::Cos(vpr.Angle(vlc.Vect()));
1758 
1759  if (fCallKFVertexing){
1760  kfResult = CallKFVertexing(part, v0part, bachelor, mcArray, &V0KF[0], &errV0KF[0], &LcKF[0], &errLcKF[0], &distances[0], &armPolKF[0]);
1761  AliDebug(2, Form("Result from KF = %d", kfResult));
1762  }
1763 
1764  /*
1765  for (Int_t i = 0; i< 3; i++){
1766  Printf("i = %d, V0KF = %f, errV0KF = %f, LcKF = %f, errLcKF = %f", V0KF[i], errV0KF[i], LcKF[i], errLcKF[i]);
1767  }
1768  */
1769 
1770  Double_t countTreta1corr = fNTracklets_1;
1771 
1772  AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
1773  // TProfile* estimatorAvg = GetEstimatorHistogram(aodEvent);
1774  // if(estimatorAvg) {
1775  // countTreta1corr = static_cast<Int_t>(AliVertexingHFUtils::GetCorrectedNtracklets(estimatorAvg, fNTracklets_1, fVtx1->GetZ(), fRefMult));
1776  // }
1777 
1778 
1779  AliAODVertex *primvert = dynamic_cast<AliAODVertex*>(part->GetPrimaryVtx());
1780 
1781  Double_t d0z0bach[2], covd0z0bach[3];
1782  bachelor->PropagateToDCA(primvert, fBField, kVeryBig, d0z0bach, covd0z0bach);
1783  Double_t tx[3];
1784  bachelor->GetXYZ(tx);
1785  tx[0] -= primvert->GetX();
1786  tx[1] -= primvert->GetY();
1787  tx[2] -= primvert->GetZ();
1788  Double_t innerpro = tx[0]*part->Px()+tx[1]*part->Py();
1789  Double_t signd0 = 1.;
1790  if(innerpro<0.) signd0 = -1.;
1791 
1792  signd0 = signd0*TMath::Abs(d0z0bach[0]);
1793 
1794 
1795  if (fUseMCInfo) { // save full tree if on MC
1796  fCandidateVariables[0] = invmassLc;
1797  fCandidateVariables[1] = invmassLc2Lpi;
1798  fCandidateVariables[2] = invmassK0s;
1799  fCandidateVariables[3] = invmassLambda;
1800  fCandidateVariables[4] = invmassLambdaBar;
1802  fCandidateVariables[6] = dcaV0;
1803  fCandidateVariables[7] = part->Getd0Prong(0);
1804  fCandidateVariables[8] = part->Getd0Prong(1);
1805  fCandidateVariables[9] = nSigmaTPCpr;
1806  fCandidateVariables[10] = nSigmaTOFpr;
1807  fCandidateVariables[11] = bachelor->Pt();
1808  fCandidateVariables[12] = v0pos->Pt();
1809  fCandidateVariables[13] = v0neg->Pt();
1810  fCandidateVariables[14] = v0part->Getd0Prong(0);
1811  fCandidateVariables[15] = v0part->Getd0Prong(1);
1812  fCandidateVariables[16] = v0part->Pt();
1813  fCandidateVariables[17] = v0part->InvMass2Prongs(0,1,11,11);
1814  fCandidateVariables[18] = part->Pt();
1815  fCandidateVariables[19] = probProton;
1816  fCandidateVariables[20] = part->Eta();
1817  fCandidateVariables[21] = v0pos->Eta();
1818  fCandidateVariables[22] = v0neg->Eta();
1819  fCandidateVariables[23] = probProtonTPC;
1820  fCandidateVariables[24] = probProtonTOF;
1821  fCandidateVariables[25] = bachelor->Eta();
1822  fCandidateVariables[26] = part->P();
1823  fCandidateVariables[27] = bachelor->P();
1824  fCandidateVariables[28] = v0part->P();
1825  fCandidateVariables[29] = v0pos->P();
1826  fCandidateVariables[30] = v0neg->P();
1827  fCandidateVariables[31] = v0part->Eta();
1828  fCandidateVariables[32] = part->DecayLength();
1829  fCandidateVariables[33] = part->DecayLengthV0();
1830  fCandidateVariables[34] = bachCode;
1831  fCandidateVariables[35] = k0SCode;
1832  fCandidateVariables[36] = v0part->AlphaV0();
1833  fCandidateVariables[37] = v0part->PtArmV0();
1834 
1835  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)));
1836  fCandidateVariables[38] = cts;
1837  fCandidateVariables[39] = weightPythia;
1838  fCandidateVariables[40] = weight5LHC13d3;
1839  fCandidateVariables[41] = weight5LHC13d3Lc;
1840  fCandidateVariables[42] = weightNch;
1842  fCandidateVariables[44] = countTreta1corr;
1843  fCandidateVariables[45] = signd0;
1846  }
1847  else { //remove MC-only variables from tree if data
1848  fCandidateVariables[0] = invmassLc;
1849  fCandidateVariables[1] = invmassLc2Lpi;
1850  fCandidateVariables[2] = invmassK0s;
1851  fCandidateVariables[3] = invmassLambda;
1852  fCandidateVariables[4] = invmassLambdaBar;
1854  fCandidateVariables[6] = dcaV0;
1855  fCandidateVariables[7] = part->Getd0Prong(0);
1856  fCandidateVariables[8] = part->Getd0Prong(1);
1857  fCandidateVariables[9] = nSigmaTPCpr;
1858  fCandidateVariables[10] = nSigmaTOFpr;
1859  fCandidateVariables[11] = bachelor->Pt();
1860  fCandidateVariables[12] = v0pos->Pt();
1861  fCandidateVariables[13] = v0neg->Pt();
1862  fCandidateVariables[14] = v0part->Getd0Prong(0);
1863  fCandidateVariables[15] = v0part->Getd0Prong(1);
1864  fCandidateVariables[16] = v0part->Pt();
1865  fCandidateVariables[17] = v0part->InvMass2Prongs(0,1,11,11);
1866  fCandidateVariables[18] = part->Pt();
1867  fCandidateVariables[19] = probProton;
1868  fCandidateVariables[20] = v0pos->Eta();
1869  fCandidateVariables[21] = v0neg->Eta();
1870  fCandidateVariables[22] = bachelor->Eta();
1871  fCandidateVariables[23] = v0part->P();
1872  fCandidateVariables[24] = part->DecayLengthV0();
1873  fCandidateVariables[25] = v0part->AlphaV0();
1874  fCandidateVariables[26] = v0part->PtArmV0();
1876  fCandidateVariables[28] = countTreta1corr;
1877  fCandidateVariables[29] = cts;
1878  fCandidateVariables[30] = signd0;
1881  }
1882 
1883  // fill multiplicity histograms for events with a candidate
1884  //fHistNtrUnCorrEvWithCand->Fill(fNTracklets_1, weightNch);
1885  //fHistNtrCorrEvWithCand->Fill(countTreta1corr, weightNch);
1886  if (fUseMCInfo) {
1887  if (isLc){
1888  AliDebug(2, Form("Reco particle %d --> Filling Sgn", iLctopK0s));
1889  if(fFillTree) fVariablesTreeSgn->Fill();
1890  fHistoCodesSgn->Fill(bachCode, k0SCode);
1891  }
1892  else {
1893  if (fFillOnlySgn == kFALSE){
1894  AliDebug(2, "Filling Bkg");
1895  fVariablesTreeBkg->Fill();
1896  if(fFillTree) fHistoCodesBkg->Fill(bachCode, k0SCode);
1897  }
1898  }
1899  }
1900  else {
1901  if(fFillTree) fVariablesTreeSgn->Fill();
1902  }
1903 
1904 
1905  if(!fFillTree){
1906 
1907  std::vector<Double_t> inputVars(9);
1908  inputVars[0] = invmassK0s;
1909  inputVars[1] = part->Getd0Prong(0);
1910  inputVars[2] = part->Getd0Prong(1);
1911  inputVars[3] = bachelor->Pt();
1912  inputVars[4] = probProton;
1913  inputVars[5] = (part->DecayLengthV0())*0.497/(v0part->P());
1914  inputVars[6] = part->CosV0PointingAngle();
1915  inputVars[7] = cts;
1916  inputVars[8] = signd0;
1917 
1918 
1919  Double_t BDTResponse = -1;
1920  BDTResponse = fBDTReader->GetMvaValue(inputVars);
1921  fBDTHisto->Fill(BDTResponse, invmassLc);
1922  if (fDebugHistograms) {
1923  fBDTHistoVsMassK0S->Fill(BDTResponse, invmassK0s);
1924  fBDTHistoVstImpParBach->Fill(BDTResponse, part->Getd0Prong(0));
1925  fBDTHistoVstImpParV0->Fill(BDTResponse, part->Getd0Prong(1));
1926  fBDTHistoVsBachelorPt->Fill(BDTResponse, bachelor->Pt());
1927  fBDTHistoVsCombinedProtonProb->Fill(BDTResponse, probProton);
1928  fBDTHistoVsCtau->Fill(BDTResponse, (part->DecayLengthV0())*0.497/(v0part->P()));
1929  fBDTHistoVsCosPAK0S->Fill(BDTResponse, part->CosV0PointingAngle());
1930  fBDTHistoVsSignd0->Fill(BDTResponse, signd0);
1931  fBDTHistoVsCosThetaStar->Fill(BDTResponse, cts);
1932 
1933  fHistoNsigmaTPC->Fill(bachelor->P(), nSigmaTPCpr);
1934  fHistoNsigmaTOF->Fill(bachelor->P(), nSigmaTOFpr);
1935  }
1936  }
1937 
1938  }
1939 
1940  return;
1941 
1942 
1943 }
1944 
1945 //________________________________________________________________________
1946 Int_t AliAnalysisTaskSELc2V0bachelorTMVAApp::CallKFVertexing(AliAODRecoCascadeHF *cascade, AliAODv0* v0part, AliAODTrack* bach, TClonesArray *mcArray,
1947  Double_t* V0KF, Double_t* errV0KF, Double_t* LcKF, Double_t* errLcKF,
1948  Double_t* distances, Double_t* armPolKF) {
1949 
1950  //
1953  //
1954 
1955  Int_t codeKFV0 = -1, codeKFLc = -1;
1956 
1957  AliKFVertex primVtxCopy;
1958  Int_t nt = 0, ntcheck = 0;
1959  Double_t pos[3] = {0., 0., 0.};
1960 
1961  fVtx1->GetXYZ(pos);
1962  Int_t contr = fVtx1->GetNContributors();
1963  Double_t covmatrix[6] = {0.};
1964  fVtx1->GetCovarianceMatrix(covmatrix);
1965  Double_t chi2 = fVtx1->GetChi2();
1966  AliESDVertex primaryESDVtxCopy(pos, covmatrix, chi2, contr, "Vertex");
1967 
1968  // topological constraint
1969  primVtxCopy = AliKFVertex(primaryESDVtxCopy);
1970  nt = primaryESDVtxCopy.GetNContributors();
1971  ntcheck = nt;
1972 
1973  Int_t pdg[2] = {211, -211};
1974  Int_t pdgLc[2] = {2212, 310};
1975 
1976  Int_t pdgDgV0toDaughters[2] = {211, 211};
1977 
1978  Int_t mcLabelV0 = v0part->MatchToMC(310, mcArray, 2, pdgDgV0toDaughters);
1979 
1980  // the KF vertex for the V0 has to be built with the prongs of the V0!
1981  Bool_t isMCokV0 = kTRUE, isBkgV0 = kFALSE;
1982  AliKFParticle V0, positiveV0KF, negativeV0KF;
1983  Int_t labelsv0daugh[2] = {-1, -1};
1984  Int_t idv0daugh[2] = {-1, -1};
1985  AliExternalTrackParam* esdv0Daugh1 = 0x0;
1986  AliExternalTrackParam* esdv0Daugh2 = 0x0;
1987  for(Int_t ipr = 0; ipr < 2; ipr++){ // 0 is positive, 1 is negative
1988  AliAODTrack *aodTrack = (AliAODTrack*)v0part->GetDaughter(ipr);
1989  if(!aodTrack) {
1990  AliDebug(2, "No V0 daughters available");
1991  return -1;
1992  }
1993  Double_t xyz[3], pxpypz[3], cv[21];
1994  Short_t sign;
1995  aodTrack->GetXYZ(xyz);
1996  aodTrack->PxPyPz(pxpypz);
1997  aodTrack->GetCovarianceXYZPxPyPz(cv);
1998  sign = aodTrack->Charge();
1999  AliExternalTrackParam tmp1( xyz, pxpypz, cv, sign);
2000 
2001  if (ipr == 0) esdv0Daugh1 = new AliExternalTrackParam( xyz, pxpypz, cv, sign);
2002  else esdv0Daugh2 = new AliExternalTrackParam( xyz, pxpypz, cv, sign);
2003  labelsv0daugh[ipr] = TMath::Abs(aodTrack->GetLabel());
2004  idv0daugh[ipr] = aodTrack->GetID();
2005  if (labelsv0daugh[ipr] == -1) isBkgV0 = kTRUE;
2006 
2007  //Printf("v0 daughter %d has label %d", ipr, labelsv0daugh[ipr]);
2008 
2009  AliKFParticle daughterKF(*aodTrack, pdg[ipr]); // we assume that the PDG is correct
2010  if (aodTrack->Charge() > 0) { // assigning positive and negative track to KF V0 for Armenteros-Podolanski plot
2011  positiveV0KF = daughterKF;
2012  }
2013  else {
2014  negativeV0KF = daughterKF;
2015  }
2016  }
2017 
2018  Double_t xn = 0., xp = 0.;//, dca;
2019  AliDebug(2, Form("bField = %f, esdv0Daugh1 = %p, esdv0Daugh2 = %p", fBField, esdv0Daugh1, esdv0Daugh2));
2020  // dca = esdv0Daugh1->GetDCA(esdv0Daugh2, fBField, xn, xp);
2021 
2022  AliExternalTrackParam tr1(*esdv0Daugh1);
2023  AliExternalTrackParam tr2(*esdv0Daugh2);
2024  tr1.PropagateTo(xn, fBField);
2025  tr2.PropagateTo(xp, fBField);
2026 
2027  AliKFParticle daughterKF1(tr1, 211);
2028  AliKFParticle daughterKF2(tr2, 211);
2029  V0.AddDaughter(positiveV0KF);
2030  V0.AddDaughter(negativeV0KF);
2031  //V0.AddDaughter(daughterKF1);
2032  //V0.AddDaughter(daughterKF2);
2033 
2034  delete esdv0Daugh1;
2035  delete esdv0Daugh2;
2036  esdv0Daugh1=0;
2037  esdv0Daugh2=0;
2038  // Checking the quality of the KF V0 vertex
2039  if( V0.GetNDF() < 1 ) {
2040  //Printf("Number of degrees of freedom < 1, continuing");
2041  return -1;
2042  }
2043  if( TMath::Sqrt(TMath::Abs(V0.GetChi2()/V0.GetNDF())) > fCutKFChi2NDF ) {
2044  //Printf("Chi2 per DOF too big, continuing");
2045  return -1;
2046  }
2047 
2048  if(ftopoConstraint && nt > 0){
2049  for(Int_t ipr = 0; ipr < 2; ipr++){ // 0 is positive, 1 is negative
2050  AliAODTrack *aodTrack = (AliAODTrack*)v0part->GetDaughter(ipr);
2051  //* subtruct daughters from primary vertex
2052  if(!aodTrack->GetUsedForPrimVtxFit()) {
2053  //Printf("Track %d was not used for primary vertex, continuing", i);
2054  continue;
2055  }
2056  AliKFParticle daughterKF(*aodTrack, pdg[ipr]); // we assume that the PDG is correct
2057  primVtxCopy -= daughterKF;
2058  ntcheck--;
2059  }
2060  }
2061 
2062  // Check V0 Chi^2 deviation from primary vertex // not needed for V0 for Lc decay!!
2063  /*
2064  if( V0.GetDeviationFromVertex( primVtxCopy ) < fCutKFDeviationFromVtxV0) {
2065  //Printf("Deviation from vertex too big, continuing");
2066  return -1;
2067  }
2068  */
2069 
2070  //* Get V0 invariant mass
2071  Double_t massV0 = 999999, sigmaMassV0 = 999999;
2072  Int_t retMV0 = V0.GetMass( massV0, sigmaMassV0 );
2073  if( retMV0 ) {
2074  if (massV0 < 0) {
2075  codeKFV0 = 1; // Mass not ok
2076  if (sigmaMassV0 > 1e19) codeKFV0 = 5; // Mass and SigmaMass not ok
2077  }
2078  else if (sigmaMassV0 > 1e19) codeKFV0 = 2; // SigmaMass not ok
2079  }
2080  fHistoMassKFV0->Fill(massV0, sigmaMassV0);
2081 
2082  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]);
2083  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]);
2084 
2085  Printf("Vertices: KF: x = %f, y = %f, z = %f", V0.GetX(), V0.GetY(), V0.GetZ());
2086  Printf("Vertices: AOD: x = %f, y = %f, z = %f", v0part->Xv(), v0part->Yv(), v0part->Zv());
2087 
2088  //Printf("Got MC vtx for V0");
2089  if (fUseMCInfo && TMath::Abs(labelsv0daugh[0] - labelsv0daugh[1]) == 1) {
2090  AliAODMCParticle* tmpdaughv01 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelsv0daugh[0]));
2091  AliAODMCParticle* tmpdaughv02 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelsv0daugh[1]));
2092  if (!tmpdaughv01 && labelsv0daugh[0] > 0){
2093  AliDebug(2, "Could not access MC info for first daughter of V0, continuing");
2094  }
2095  if (!tmpdaughv02 && labelsv0daugh[1] > 0){
2096  AliDebug(2, "Could not access MC info for second daughter of V0, continuing");
2097  }
2098  if(tmpdaughv01){
2099  Double_t xPionMC = tmpdaughv01->Xv(); //Production vertex of Pion --> Where K0S decays
2100  Double_t yPionMC = tmpdaughv01->Yv();
2101  Double_t zPionMC = tmpdaughv01->Zv();
2102  //Printf("Got MC vtx for Pion");
2103  Printf("Vertices: MC: x = %f, y = %f, z = %f", xPionMC, yPionMC, zPionMC);
2104  }
2105  }
2106  else {
2107  Printf("Not a true V0");
2108  }
2109  //massV0=-1;//return -1;// !!!!
2110 
2111  // now use what we just try with the bachelor, to build the Lc
2112 
2113  // topological constraint
2114  nt = primVtxCopy.GetNContributors();
2115  ntcheck = nt;
2116 
2117  Bool_t isMCokLc = kTRUE, isBkgLc = kFALSE;
2118  AliKFParticle Lc;
2119  Int_t labelsLcdaugh[2] = {-1, -1};
2120  labelsLcdaugh[0] = TMath::Abs(bach->GetLabel());
2121  labelsLcdaugh[1] = mcLabelV0;
2122 
2123  if (bach->Charge() < 0) pdgLc[0] = -pdgLc[0];
2124  AliKFParticle daughterKFLc(*bach, pdgLc[0]);
2125  Lc.AddDaughter(daughterKFLc);
2126  TParticlePDG* particlePDG = TDatabasePDG::Instance()->GetParticle(310);
2127  Double_t massPDGK0S = particlePDG->Mass();
2128  V0.SetMassConstraint(massPDGK0S);
2129  Lc.AddDaughter(V0);
2130  if( Lc.GetNDF() < 1 ) {
2131  AliDebug(2, Form("Lc: Number of degrees of freedom < 1 (%d), continuing", Lc.GetNDF()));
2132  return -1;
2133  }
2134  if( TMath::Sqrt(TMath::Abs(Lc.GetChi2()/Lc.GetNDF())) > fCutKFChi2NDF) {
2135  AliDebug(2, Form("Lc: Chi2 per DOF too big, continuing (%f)", TMath::Sqrt(TMath::Abs(Lc.GetChi2()/Lc.GetNDF()))));
2136  return -1;
2137  }
2138 
2139  if(ftopoConstraint && nt > 0){
2140  //* subtruct daughters from primary vertex
2141  if(!bach->GetUsedForPrimVtxFit()) {
2142  AliDebug(3, "Lc: Bachelor was not used for primary vertex, not subtracting it from primary vertex");
2143  }
2144  else{
2145  primVtxCopy -= daughterKFLc;
2146  ntcheck--;
2147  }
2148  /* the V0 was added above, so it is ok to remove it without checking
2149  if(!V0->GetUsedForPrimVtxFit()) {
2150  Printf("Lc: V0 was not used for primary vertex, continuing");
2151  continue;
2152  }
2153  */
2154  //primVtxCopy -= V0;
2155  //ntcheck--;
2156  }
2157 
2158  // Check Lc Chi^2 deviation from primary vertex
2159  /*
2160  if( Lc.GetDeviationFromVertex( primVtxCopy ) > fCutKFDeviationFromVtx) {
2161  AliDebug(2, Form("Lc: Deviation from vertex too big, continuing (%f)", Lc.GetDeviationFromVertex( primVtxCopy )));
2162  return -1;
2163  }
2164  if(ftopoConstraint){
2165  if(ntcheck>0) {
2166  // Add Lc to primary vertex to improve the primary vertex resolution
2167  primVtxCopy += Lc;
2168  Lc.SetProductionVertex(primVtxCopy);
2169  }
2170  }
2171  */
2172  //* Check chi^2
2173  if( TMath::Sqrt( TMath::Abs(Lc.GetChi2()/Lc.GetNDF())) > fCutKFChi2NDF) {
2174  AliDebug(2, Form("Lc: Final Chi2 per DOF too big, continuing (%f)", TMath::Sqrt( TMath::Abs(Lc.GetChi2()/Lc.GetNDF()))));
2175  return -1;
2176  }
2177 
2178  if(ftopoConstraint){
2179  V0.SetProductionVertex(Lc);
2180  }
2181 
2182  // After setting the vertex of the V0, getting/filling some info
2183 
2184  //* Get V0 decayLength
2185  Double_t decayLengthV0 = 999999, sigmaDecayLengthV0 = 999999;
2186  Int_t retDLV0 = V0.GetDecayLength( decayLengthV0, sigmaDecayLengthV0 );
2187  if( retDLV0 ) {
2188  if (sigmaDecayLengthV0 > 1e19) {
2189  codeKFV0 = 3; // DecayLength not ok
2190  if (massV0 < 0) {
2191  codeKFV0 = 6; // DecayLength and Mass not ok
2192  if (sigmaMassV0 > 1e19) codeKFV0 = 11; // DecayLength and Mass and SigmaMass not ok
2193  }
2194  else if (sigmaMassV0 > 1e19) codeKFV0 = 8; // DecayLength and SigmaMass not ok
2195  }
2196  }
2197  fHistoDecayLengthKFV0->Fill(decayLengthV0, sigmaDecayLengthV0);
2198 
2199  //* Get V0 life time
2200  Double_t lifeTimeV0 = 999999, sigmaLifeTimeV0 = 999999;
2201  Int_t retTLV0 = V0.GetLifeTime( lifeTimeV0, sigmaLifeTimeV0 );
2202  if( retTLV0 ) {
2203  if (sigmaLifeTimeV0 > 1e19) {
2204  codeKFV0 = 4; // LifeTime not ok
2205  if (sigmaDecayLengthV0 > 1e19) {
2206  codeKFV0 = 9; // LifeTime and DecayLength not ok
2207  if (massV0 < 0) {
2208  codeKFV0 = 14; // LifeTime and DecayLength and Mass not ok
2209  if (sigmaMassV0 > 1e19) codeKFV0 = 15; // LifeTime and DecayLength and Mass and SigmaMass not ok
2210  }
2211  else if (sigmaMassV0 > 1e19) codeKFV0 = 13; // LifeTime and DecayLength and SigmaMass not ok
2212  }
2213  else if (massV0 < 0) { // LifeTime and Mass and SigmaMass not ok
2214  codeKFV0 = 7; // LifeTime and Mass not ok
2215  if (sigmaMassV0 > 1e19) codeKFV0 = 12; // LifeTime and Mass and SigmaMass not ok
2216  }
2217  else if (sigmaMassV0 > 1e19) codeKFV0 = 10; // LifeTime and SigmaMass not ok
2218  }
2219  }
2220  fHistoLifeTimeKFV0->Fill(lifeTimeV0, sigmaLifeTimeV0);
2221 
2222  if (codeKFV0 == -1) codeKFV0 = 0;
2223  fHistoKFV0->Fill(codeKFV0);
2224 
2225  AliDebug(2, Form("V0: mass = %f, decay length = %f, life time = %f", massV0, decayLengthV0, lifeTimeV0 ));
2226 
2227  fHistoMassV0All->Fill(massV0);
2228  fHistoDecayLengthV0All->Fill(decayLengthV0);
2229  fHistoLifeTimeV0All->Fill(lifeTimeV0);
2230 
2231  Double_t qtAlphaV0[2] = {0., 0.};
2232  Double_t vtxV0KF[3] = {V0.GetX(), V0.GetY(), V0.GetZ()};
2233  positiveV0KF.TransportToPoint(vtxV0KF);
2234  negativeV0KF.TransportToPoint(vtxV0KF);
2235  V0.GetArmenterosPodolanski(positiveV0KF, negativeV0KF, qtAlphaV0);
2236  AliDebug(2, Form("Armenteros-Podolanski variables: alpha = %f, qt = %f", qtAlphaV0[1], qtAlphaV0[0]));
2237  fHistoArmenterosPodolanskiV0KF->Fill(qtAlphaV0[1], qtAlphaV0[0]);
2238  fHistoArmenterosPodolanskiV0AOD->Fill(v0part->AlphaV0(), v0part->PtArmV0());
2239  armPolKF[0] = qtAlphaV0[1];
2240  armPolKF[1] = qtAlphaV0[0];
2241 
2242  // Checking MC info for V0
2243 
2244  AliAODMCParticle *motherV0 = 0x0;
2245  AliAODMCParticle *daughv01 = 0x0;
2246  AliAODMCParticle *daughv02 = 0x0;
2247 
2248  if (fUseMCInfo) {
2249  daughv01 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelsv0daugh[0]));
2250  daughv02 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelsv0daugh[1]));
2251  if (!daughv01 && labelsv0daugh[0] > 0){
2252  AliDebug(2, "Could not access MC info for first daughter of V0, continuing");
2253  isMCokV0 = kFALSE;
2254  }
2255  if (!daughv02 && labelsv0daugh[1] > 0){
2256  AliDebug(2, "Could not access MC info for second daughter of V0, continuing");
2257  isMCokV0 = kFALSE;
2258  }
2259  if (isMCokV0){
2260  if( daughv01->GetMother() == daughv02->GetMother() && daughv01->GetMother()>=0 ){
2261  AliDebug(3, Form("The mother has label %d", daughv01->GetMother()));
2262  motherV0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughv01->GetMother()));
2263  if( motherV0 && TMath::Abs(motherV0->GetPdgCode()) != 21 ){ // These are all V0 that are truly V0, not only K0S, but no gluons
2264  if( motherV0->GetNDaughters() == 2 ){
2265  fHistoMassV0True->Fill(massV0);
2266  fHistoDecayLengthV0True->Fill(decayLengthV0);
2267  fHistoLifeTimeV0True->Fill(lifeTimeV0);
2268  fHistoMassV0TrueFromAOD->Fill(v0part->MassK0Short());
2269  if (TMath::Abs(motherV0->GetPdgCode()) == 310){ // These are true V0 that are also K0S
2270  fHistoMassV0TrueK0S->Fill(massV0);
2271  fHistoDecayLengthV0TrueK0S->Fill(decayLengthV0);
2272  fHistoLifeTimeV0TrueK0S->Fill(lifeTimeV0);
2273  fHistoMassV0TrueK0SFromAOD->Fill(v0part->MassK0Short());
2274  }
2275  }
2276  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()));
2277  }
2278  else if (!motherV0){
2279  AliDebug(3, "could not access MC info for mother, continuing");
2280  }
2281  else {
2282  AliDebug(3, "MC mother is a gluon, continuing");
2283  }
2284  }
2285  else {
2286  AliDebug(3, "Background V0!");
2287  isBkgV0 = kTRUE;
2288  }
2289  }
2290  }
2291 
2292  AliDebug(2, Form("isMCokV0 = %d, isBkgV0 = %d", (Int_t)isMCokV0, (Int_t)isBkgV0));
2293 
2294  // Going back to Lc
2295 
2296  //* Get Lc invariant mass
2297  Double_t massLc = 999999, sigmaMassLc= 999999;
2298  Int_t retMLc = Lc.GetMass( massLc, sigmaMassLc );
2299  if( retMLc ) {
2300  AliDebug(3, Form("----> Could not get mass (%e), and sigma(%e) for Lc, continuing", massLc, sigmaMassLc));
2301  if (massLc < 0) {
2302  codeKFLc = 1; // Mass not ok
2303  if (sigmaMassLc > 1e19) codeKFLc = 5; // Mass and SigmaMass not ok
2304  }
2305  else if (sigmaMassLc > 1e19) codeKFLc = 2; // SigmaMass not ok
2306  }
2307  fHistoMassKFLc->Fill(massLc, sigmaMassLc);
2308 
2309  //* Get Lc Decay length
2310  Double_t decayLengthLc = 999999, sigmaDecayLengthLc = 999999;
2311  Int_t retDLLc = Lc.GetDecayLength( decayLengthLc, sigmaDecayLengthLc );
2312  if( retDLLc ) {
2313  AliDebug(3, "----> Lc: Could not get decay length, and sigma");
2314  if (sigmaDecayLengthLc > 1e19) {
2315  codeKFLc = 3; // DecayLength not ok
2316  if (massLc < 0) {
2317  codeKFLc = 6; // DecayLength and Mass not ok
2318  if (sigmaMassLc > 1e19) codeKFLc = 11; // DecayLength and Mass and SigmaMass not ok
2319  }
2320  else if (sigmaMassLc > 1e19) codeKFLc = 8; // DecayLength and SigmaMass not ok
2321  }
2322  }
2323  AliDebug(3, Form("retDLLc = %d, with decayLength = %f and error = %e", retDLLc, decayLengthLc, sigmaDecayLengthLc));
2324 
2325  fHistoDecayLengthKFLc->Fill(decayLengthLc, sigmaDecayLengthLc);
2326 
2327  //* Get Lc life time
2328  Double_t lifeTimeLc = 999999, sigmaLifeTimeLc = 999999;
2329  Int_t retTLLc = Lc.GetLifeTime( lifeTimeLc, sigmaLifeTimeLc );
2330  if( retTLLc ) {
2331  AliDebug(3, "----> Lc: Could not get lifeTime, and sigma");
2332  if (sigmaLifeTimeLc > 1e19) {
2333  codeKFLc = 4; // LifeTime not ok
2334  if (sigmaDecayLengthLc > 1e19) {
2335  codeKFLc = 9; // LifeTime and DecayLength not ok
2336  if (massLc < 0) {
2337  codeKFLc = 14; // LifeTime and DecayLength and Mass not ok
2338  if (sigmaMassLc > 1e19) codeKFLc = 15; // LifeTime and DecayLength and Mass and SigmaMass not ok
2339  }
2340  else if (sigmaMassLc > 1e19) codeKFLc = 13; // LifeTime and DecayLength and SigmaMass not ok
2341  }
2342  else if (massLc < 0) { // LifeTime and Mass and SigmaMass not ok
2343  codeKFLc = 7; // LifeTime and Mass not ok
2344  if (sigmaMassLc > 1e19) codeKFLc = 12; // LifeTime and Mass and SigmaMass not ok
2345  }
2346  else if (sigmaMassLc > 1e19) codeKFLc = 10; // LifeTime and SigmaMass not ok
2347  }
2348  }
2349 
2350  fHistoLifeTimeKFLc->Fill(lifeTimeLc, sigmaLifeTimeLc);
2351 
2352  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));
2353 
2354  if (codeKFLc == -1) codeKFLc = 0;
2355  fHistoKFLc->Fill(codeKFLc);
2356 
2357  fHistoKF->Fill(codeKFV0, codeKFLc);
2358 
2359  // here we fill the histgrams for all the reconstructed KF vertices for the cascade
2360  fHistoMassLcAll->Fill(massLc);
2361  fHistoDecayLengthLcAll->Fill(decayLengthLc);
2362  fHistoLifeTimeLcAll->Fill(lifeTimeLc);
2363 
2364  fHistoMassV0fromLcAll->Fill(massV0);
2365  fHistoDecayLengthV0fromLcAll->Fill(decayLengthV0);
2366  fHistoLifeTimeV0fromLcAll->Fill(lifeTimeV0);
2367 
2368  Double_t xV0 = V0.GetX();
2369  Double_t yV0 = V0.GetY();
2370  Double_t zV0 = V0.GetZ();
2371 
2372  Double_t xLc = Lc.GetX();
2373  Double_t yLc = Lc.GetY();
2374  Double_t zLc = Lc.GetZ();
2375 
2376  Double_t xPrimVtx = primVtxCopy.GetX();
2377  Double_t yPrimVtx = primVtxCopy.GetY();
2378  Double_t zPrimVtx = primVtxCopy.GetZ();
2379 
2380  Double_t distanceLcToPrimVtx = TMath::Sqrt((xPrimVtx - xLc) * (xPrimVtx - xLc) +
2381  (yPrimVtx - yLc) * (yPrimVtx - yLc) +
2382  (zPrimVtx - zLc) * (zPrimVtx - zLc));
2383 
2384  Double_t distanceV0ToPrimVtx = TMath::Sqrt((xPrimVtx - xV0) * (xPrimVtx - xV0) +
2385  (yPrimVtx - yV0) * (yPrimVtx - yV0) +
2386  (zPrimVtx - zV0) * (zPrimVtx - zV0));
2387 
2388  Double_t distanceV0ToLc = TMath::Sqrt((xLc - xV0)*(xLc - xV0) +
2389  (yLc - yV0)*(yLc - yV0) +
2390  (zLc - zV0)*(zLc - zV0));
2391 
2392  //Printf("distanceLcToPrimVtx = %e, distanceV0ToPrimVtx= %f, distanceV0ToLc = %f", distanceLcToPrimVtx, distanceV0ToPrimVtx, distanceV0ToLc);
2393 
2394  fHistoDistanceLcToPrimVtx->Fill(distanceLcToPrimVtx);
2395  fHistoDistanceV0ToPrimVtx->Fill(distanceV0ToPrimVtx);
2396  fHistoDistanceV0ToLc->Fill(distanceV0ToLc);
2397 
2398  distances[0] = distanceLcToPrimVtx;
2399  distances[1] = distanceV0ToPrimVtx;
2400  distances[2] = distanceV0ToLc;
2401 
2402  if (fUseMCInfo) {
2403 
2404  AliAODMCParticle *daughv01Lc = 0x0;
2405  AliAODMCParticle *K0S = 0x0;
2406  AliAODMCParticle *daughv02Lc = 0x0;
2407 
2408  if (labelsLcdaugh[0] >= 0) {
2409  // Printf("Getting Bachelor from label %d", labelsLcdaugh[1]);
2410  daughv01Lc = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelsLcdaugh[0]));
2411  if (!daughv01Lc){
2412  AliDebug(3, "Could not access MC info for first daughter of Lc");
2413  isMCokLc = kFALSE;
2414  }
2415  else {
2416  AliDebug(2, Form("The bachelor has label = %d", daughv01Lc->GetLabel()));
2417  if (TMath::Abs(daughv01Lc->GetPdgCode()) != 2212) isBkgLc = kTRUE;
2418  }
2419  }
2420  else { // The bachelor is a fake
2421  isBkgLc = kTRUE;
2422  }
2423 
2424  if (labelsLcdaugh[1] >= 0) {
2425  //Printf("Getting K0S");
2426  K0S = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelsLcdaugh[1]));
2427  if (!K0S) {
2428  AliDebug(3, "Could not access MC info for second daughter of Lc");
2429  isMCokLc = kFALSE;
2430  }
2431  else{
2432  if (TMath::Abs(K0S->GetPdgCode()) != 310) isBkgLc = kTRUE;
2433  }
2434  }
2435  else{
2436  AliDebug(2, "The K0S is not true --> it does not have a label, continuing...");
2437  isBkgLc = kTRUE;
2438  }
2439 
2440  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
2441  if (isMCokLc) { // We can then access its MC info, and it might then be that also the Lc is a true Lc
2442  Int_t iK0 = K0S->GetMother();
2443  if (iK0 < 0) {
2444  Printf("The K0S has no mother... IMPOSSIBLE"); // the K0S MUST have a mother!
2445  }
2446  else { // The K0S has a mother
2447  daughv02Lc = dynamic_cast<AliAODMCParticle*>(mcArray->At(iK0));
2448  if (!daughv02Lc){
2449  AliDebug(3, "Could not access MC info for second daughter of Lc");
2450  }
2451  else { // we can access safely the K0S mother in the MC
2452  if( daughv01Lc && (daughv01Lc->GetMother() == daughv02Lc->GetMother()) && (daughv01Lc->GetMother()>=0) ){ // This is a true cascade! bachelor and V0 come from the same mother
2453  //Printf("Lc: The mother has label %d", daughv01Lc->GetMother());
2454  AliAODMCParticle *motherLc = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughv01Lc->GetMother()));
2455  Int_t pdgMum = 0, pdgBach = 0, pdgV0 = 0;
2456  if (motherLc) pdgMum = motherLc->GetPdgCode();
2457  if (daughv01Lc) pdgBach = daughv01Lc->GetPdgCode();
2458  if (daughv02Lc) pdgV0 = daughv02Lc->GetPdgCode();
2459  AliDebug(2, Form("pdgMum = %d, pdgBach = %d, pdgV0 = %d", pdgMum, pdgBach, pdgV0));
2460 
2461  if( motherLc && TMath::Abs(motherLc->GetPdgCode()) != 21 ){ // These are true cascades, we don't know yet if they are Lc
2462  fHistoMassLcTrue->Fill(massLc);
2463  fHistoDecayLengthLcTrue->Fill(decayLengthLc);
2464  fHistoLifeTimeLcTrue->Fill(lifeTimeLc);
2465  fHistoMassLcTrueFromAOD->Fill(cascade->InvMassLctoK0sP());
2466 
2467  fHistoMassV0fromLcTrue->Fill(massV0);
2468  fHistoDecayLengthV0fromLcTrue->Fill(decayLengthV0);
2469  fHistoLifeTimeV0fromLcTrue->Fill(lifeTimeV0);
2470 
2471  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...
2472  AliDebug(2, Form("IT IS SIGNAL!!! with label = %d", motherLc->GetLabel()));
2473 
2474  fHistoArmenterosPodolanskiV0KFSgn->Fill(qtAlphaV0[1], qtAlphaV0[0]);
2475  fHistoArmenterosPodolanskiV0AODSgn->Fill(v0part->AlphaV0(), v0part->PtArmV0());
2476 
2477  fHistoDistanceLcToPrimVtxSgn->Fill(distanceLcToPrimVtx);
2478  fHistoDistanceV0ToPrimVtxSgn->Fill(distanceV0ToPrimVtx);
2479  fHistoDistanceV0ToLcSgn->Fill(distanceV0ToLc);
2480 
2481  fHistoMassLcSgn->Fill(massLc);
2482  fHistoMassLcSgnFromAOD->Fill(cascade->InvMassLctoK0sP());
2483 
2484  fHistoDecayLengthLcSgn->Fill(decayLengthLc);
2485  fHistoLifeTimeLcSgn->Fill(lifeTimeLc);
2486 
2487  fHistoMassV0fromLcSgn->Fill(massV0);
2488  fHistoDecayLengthV0fromLcSgn->Fill(decayLengthV0);
2489  fHistoLifeTimeV0fromLcSgn->Fill(lifeTimeV0);
2490  }
2491  else {
2492  isBkgLc = kTRUE; // it is not a Lc, since the pdg != 4122
2493  }
2494 
2495  // if we got here, we can compare with MC information; this part is done also for cases where the cascade is not a Lc!
2496  // First, we compare the vtx of the Lc
2497  Double_t xLcMC = motherLc->Xv();
2498  Double_t yLcMC = motherLc->Yv();
2499  Double_t zLcMC = motherLc->Zv();
2500  //Printf("Got MC vtx for Lc");
2501  Double_t xProtonMC = daughv01Lc->Xv();
2502  Double_t yProtonMC = daughv01Lc->Yv();
2503  Double_t zProtonMC = daughv01Lc->Zv();
2504  //Printf("Got MC vtx for Proton");
2505 
2506  Double_t vtxLcResidualToPrimVtx = TMath::Sqrt((xLcMC - xProtonMC) * (xLcMC - xProtonMC) +
2507  (yLcMC - yProtonMC) * (yLcMC - yProtonMC) +
2508  (zLcMC - zProtonMC) * (zLcMC - zProtonMC)) - distanceLcToPrimVtx;
2509 
2510  // Then, we compare the vtx of the K0s
2511  Double_t xV0MC = motherV0->Xv(); //Production vertex of K0S --> Where Lc decays
2512  Double_t yV0MC = motherV0->Yv();
2513  Double_t zV0MC = motherV0->Zv();
2514 
2515  //Printf("Got MC vtx for V0");
2516  Double_t xPionMC = daughv01->Xv(); //Production vertex of Pion --> Where K0S decays
2517  Double_t yPionMC = daughv01->Yv();
2518  Double_t zPionMC = daughv01->Zv();
2519  //Printf("Got MC vtx for Pion");
2520  Printf("Vertices: MC: x = %f, y = %f, z = %f", xPionMC, yPionMC, zPionMC);
2521 
2522  Double_t vtxV0ResidualToLc = TMath::Sqrt((xV0MC - xPionMC) * (xV0MC - xPionMC) +
2523  (yV0MC - yPionMC) * (yV0MC - yPionMC) +
2524  (zV0MC - zPionMC) * (zV0MC - zPionMC)) - distanceV0ToLc;
2525 
2526  // Then, the K0S vertex wrt primary vertex
2527 
2528  Double_t vtxV0ResidualToPrimVtx = TMath::Sqrt((xPionMC - xLcMC) * (xPionMC - xLcMC) +
2529  (yPionMC - yLcMC) * (yPionMC - yLcMC) +
2530  (zPionMC - zLcMC) * (zPionMC - zLcMC)) - distanceV0ToPrimVtx;
2531 
2532  fHistoVtxLcResidualToPrimVtx->Fill(vtxLcResidualToPrimVtx);
2533  fHistoVtxV0ResidualToLc->Fill(vtxV0ResidualToLc);
2534  fHistoVtxV0ResidualToPrimVtx->Fill(vtxV0ResidualToPrimVtx);
2535 
2536  } //closing: if( motherLc && TMath::Abs(motherLc->GetPdgCode()) != 21 )
2537  else if (!motherLc){
2538  AliDebug(2, "We could not access MC info for Lc mother, so we did nothing");
2539  }
2540  else {
2541  AliDebug(2, "MC Lc mother is a gluon, so we did nothing");
2542  }
2543  } //closing: if( daughv01Lc->GetMother() == daughv02Lc->GetMother() && daughv01Lc->GetMother()>=0 )
2544  else {
2545  isBkgLc = kTRUE; // it cannot be a Lc, since the daughters do not have the same mother
2546  }
2547  } // closing: else { // we can access safely the K0S mother in the MC
2548  } // closing: else { // The K0S has a mother
2549  } // closing isMCLcok
2550  } // closing !isBkgLc
2551  } // closing fUseMCInfo
2552 
2553  //Printf("retMV0 = %d, retMLc = %d", retMV0, retMLc);
2554  if ( retMV0 == 0 && retMLc == 0){
2555  V0KF[0] = massV0;
2556  errV0KF[0] = sigmaMassV0;
2557  V0KF[1] = decayLengthV0;
2558  errV0KF[1] = sigmaDecayLengthV0;
2559  V0KF[2] = lifeTimeV0;
2560  errV0KF[2] = sigmaLifeTimeV0;
2561  LcKF[0] = massLc;
2562  errLcKF[0] = sigmaMassLc;
2563  LcKF[1] = decayLengthLc;
2564  errLcKF[1] = sigmaDecayLengthLc;
2565  LcKF[2] = lifeTimeLc;
2566  errLcKF[2] = sigmaLifeTimeLc;
2567  }
2568 
2569  return 1;
2570 
2571 }
2572 //________________________________________________________________________
2574  AliAODTrack* bachelor,
2575  TClonesArray *mcArray ){
2576 
2577  //Printf("In CheckBachelor");
2578 
2581 
2582  Int_t label = bachelor->GetLabel();
2583  if (label == -1) {
2584  return kBachFake;
2585  }
2586 
2587  AliAODMCParticle *mcpart = dynamic_cast<AliAODMCParticle*>(mcArray->At(TMath::Abs(label)));
2588  if (!mcpart) return kBachInvalid;
2589  Int_t pdg = mcpart->PdgCode();
2590  if (TMath::Abs(pdg) != 2212) {
2591  AliDebug(2, Form("Bachelor is not a p, but a particle with pdg code = %d", pdg));
2592  return kBachNoProton;
2593  }
2594  else { // it is a proton
2595  //Int_t labelLc = part->GetLabel();
2596  Int_t labelLc = FindLcLabel(part, mcArray);
2597  //Printf(">>>>>>>>>>>>> label for Lc = %d", labelLc);
2598  Int_t bachelorMotherLabel = mcpart->GetMother();
2599  //Printf(">>>>>>>>>>>>> label for bachelorMother = %d", bachelorMotherLabel);
2600  if (bachelorMotherLabel == -1) {
2601  return kBachPrimary;
2602  }
2603  AliAODMCParticle *bachelorMother = dynamic_cast<AliAODMCParticle*>(mcArray->At(bachelorMotherLabel));
2604  if (!bachelorMother) return kBachInvalid;
2605  Int_t pdgMother = bachelorMother->PdgCode();
2606  if (TMath::Abs(pdgMother) != 4122) {
2607  AliDebug(2, Form("The proton does not come from a Lc, but from a particle with pdgcode = %d", pdgMother));
2608  return kBachNoLambdaMother;
2609  }
2610  else { // it comes from Lc
2611  if (labelLc != bachelorMotherLabel){
2612  //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));
2613  AliDebug(2, Form("The proton comes from a Lc, but it is not the candidate we are analyzing"));
2615  }
2616  else { // it comes from the correct Lc
2617  AliDebug(2, Form("The proton comes from a Lc, and it is EXACTLY the candidate we are analyzing"));
2618  return kBachCorrectLambdaMother;
2619  }
2620  }
2621  }
2622 
2623  return kBachInvalid;
2624 
2625 }
2626 
2627 //________________________________________________________________________
2629  AliAODv0* v0part,
2630  //AliAODTrack* v0part,
2631  TClonesArray *mcArray ){
2632 
2635 
2636  //Printf(" CheckK0S");
2637 
2638  //Int_t labelMatchToMC = v0part->MatchToMC(310, mcArray);
2639  //Int_t label = v0part->GetLabel();
2640  Int_t labelFind = FindV0Label(v0part, mcArray);
2641  //Printf("\n >>>>>>>>>>>>> label for V0 = %d, from MatchToMC = %d, from Find = %d", label, labelMatchToMC, labelFind);
2642  if (labelFind == -1) {
2643  return kK0SFake;
2644  }
2645 
2646  AliAODMCParticle *mcpart = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelFind));
2647  if (!mcpart) return kK0SInvalid;
2648  Int_t pdg = mcpart->PdgCode();
2649  if (TMath::Abs(pdg) != 310) {
2650  AliDebug(2, Form("K0Spart is not a K0S, but a particle with pdg code = %d", pdg));
2651  //AliInfo(Form("K0Spart is not a K0S, but a particle with pdg code = %d", pdg));
2652  return kK0SNoK0S;
2653  }
2654  else { // it is a K0S
2655  //Int_t labelLc = part->GetLabel();
2656  Int_t labelLc = FindLcLabel(part, mcArray);
2657  Int_t K0SpartMotherLabel = mcpart->GetMother();
2658  if (K0SpartMotherLabel == -1) {
2659  return kK0SWithoutMother;
2660  }
2661  AliAODMCParticle *K0SpartMother = dynamic_cast<AliAODMCParticle*>(mcArray->At(K0SpartMotherLabel)); // should be a K0 in case of signal
2662  if (!K0SpartMother) return kK0SInvalid;
2663  Int_t pdgMotherK0S = K0SpartMother->PdgCode();
2664  if (TMath::Abs(pdgMotherK0S) != 311) {
2665  AliDebug(2, Form("The K0S does not come from a K0, but from a particle with pdgcode = %d", pdgMotherK0S));
2666  //AliInfo(Form("The K0S does not come from a K0, but from a particle with pdgcode = %d", pdgMotherK0S));
2667  return kK0SNotFromK0;
2668  }
2669  else { // the K0S comes from a K0
2670  Int_t K0MotherLabel = K0SpartMother->GetMother(); // mother of K0 --> Lc in case of signal
2671  if (K0MotherLabel == -1) {
2672  return kK0Primary;
2673  }
2674  AliAODMCParticle *K0Mother = dynamic_cast<AliAODMCParticle*>(mcArray->At(K0MotherLabel));
2675  if (!K0Mother) return kK0SInvalid;
2676  Int_t pdgK0Mother = K0Mother->PdgCode();
2677  if (TMath::Abs(pdgK0Mother) != 4122) { // the K0 does not come from a Lc
2678  AliDebug(2, Form("The K0 does not come from a Lc, but from a particle with pdgcode = %d", pdgK0Mother));
2679  //AliInfo(Form("The K0 does not come from a Lc, but from a particle with pdgcode = %d", pdgK0Mother));
2680  return kK0NoLambdaMother;
2681  }
2682  else { // the K0 comes from Lc
2683  if (labelLc != K0MotherLabel){ // The K0 comes from a different Lc
2684  AliDebug(2, Form("The K0S comes from a Lc, but it is not the candidate we are analyzing"));
2685  //AliInfo(Form("The K0S comes from a Lc, but it is not the candidate we are analyzing"));
2686  return kK0DifferentLambdaMother;
2687  }
2688  else { // it comes from the correct Lc
2689  AliDebug(2, Form("The K0S comes from a Lc, and it is EXACTLY the candidate we are analyzing"));
2690  //AliInfo(Form("The K0S comes from a Lc, and it is EXACTLY the candidate we are analyzing"));
2691  return kK0CorrectLambdaMother;
2692  }
2693  }
2694  }
2695  }
2696 
2697  return kK0SInvalid;
2698 
2699 }
2700 
2701 //----------------------------------------------------------------------------
2702 Int_t AliAnalysisTaskSELc2V0bachelorTMVAApp::FindV0Label(AliAODRecoDecay* v0part, TClonesArray *mcArray) const
2703 {
2704 
2705  //Printf(" FindV0Label");
2706 
2708 
2709  Int_t labMother[2] = {-1, -1};
2710  AliAODMCParticle *part = 0;
2711  AliAODMCParticle *mother = 0;
2712  Int_t dgLabels = -1;
2713 
2714  Int_t ndg = v0part->GetNDaughters();
2715  if (ndg != 2) {
2716  AliFatal(Form("IMPOSSIBLE!! There are %d daughters, but they should be 2!", ndg));
2717  }
2718 
2719  for(Int_t i = 0; i < 2; i++) {
2720  AliAODTrack *trk = (AliAODTrack*)v0part->GetDaughter(i);
2721  dgLabels = trk->GetLabel();
2722  if (dgLabels == -1) {
2723  //printf("daughter with negative label %d\n", dgLabels);
2724  return -1;
2725  }
2726  part = (AliAODMCParticle*)mcArray->At(TMath::Abs(dgLabels));
2727  if (!part) {
2728  //printf("no MC particle\n");
2729  return -1;
2730  }
2731  labMother[i] = part->GetMother();
2732  if (labMother[i] != -1){
2733  mother = (AliAODMCParticle*)mcArray->At(labMother[i]);
2734  if(!mother) {
2735  //printf("no MC mother particle\n");
2736  return -1;
2737  }
2738  }
2739  else {
2740  return -1;
2741  }
2742  }
2743 
2744  if (labMother[0] == labMother[1]) return labMother[0];
2745 
2746  else return -1;
2747 
2748 }
2749 
2750 //----------------------------------------------------------------------------
2752 {
2753 
2755 
2756  //Printf(" FindLcLabel");
2757 
2758  AliAODMCParticle *part = 0;
2759  AliAODMCParticle *mother = 0;
2760  AliAODMCParticle *grandmother = 0;
2761  Int_t dgLabels = -1;
2762 
2763  Int_t ndg = cascade->GetNDaughters();
2764  if (ndg != 2) {
2765  AliFatal(Form("IMPOSSIBLE!! There are %d daughters, but they should be 2!", ndg));
2766  }
2767 
2768  // daughter 0 --> bachelor, daughter 1 --> V0
2769  AliAODTrack *trk = (AliAODTrack*)cascade->GetDaughter(0); // bachelor
2770  dgLabels = trk->GetLabel();
2771  if (dgLabels == -1 ) {
2772  //printf("daughter with negative label %d\n", dgLabels);
2773  return -1;
2774  }
2775  part = (AliAODMCParticle*)mcArray->At(TMath::Abs(dgLabels));
2776  if (!part) {
2777  //printf("no MC particle\n");
2778  return -1;
2779  }
2780  Int_t labMotherBach = part->GetMother();
2781  if (labMotherBach == -1){
2782  return -1;
2783  }
2784  mother = (AliAODMCParticle*)mcArray->At(labMotherBach);
2785  if(!mother) {
2786  //printf("no MC mother particle\n");
2787  return -1;
2788  }
2789 
2790  AliAODv0 *v0 = (AliAODv0*)cascade->GetDaughter(1); // V0
2791  dgLabels = FindV0Label(v0, mcArray);
2792  if (dgLabels == -1 ) {
2793  //printf("daughter with negative label (v0 was a fake) %d\n", dgLabels);
2794  return -1;
2795  }
2796  part = (AliAODMCParticle*)mcArray->At(TMath::Abs(dgLabels));
2797  if (!part) {
2798  //printf("no MC particle\n");
2799  return -1;
2800  }
2801  Int_t labMotherv0 = part->GetMother();
2802  if (labMotherv0 == -1){
2803  return -1;
2804  }
2805  mother = (AliAODMCParticle*)mcArray->At(labMotherv0);
2806  if(!mother) {
2807  //printf("no MC mother particle\n");
2808  return -1;
2809  }
2810  Int_t labGrandMotherv0 = mother->GetMother();
2811  if (labGrandMotherv0 == -1){
2812  return -1;
2813  }
2814  grandmother = (AliAODMCParticle*)mcArray->At(labGrandMotherv0);
2815  if(!grandmother) {
2816  //printf("no MC mother particle\n");
2817  return -1;
2818  }
2819 
2820  //Printf("labMotherBach = %d, labMotherv0 = %d, labGrandMotherv0 = %d", labMotherBach, labMotherv0, labGrandMotherv0);
2821  if (labGrandMotherv0 == labMotherBach) return labMotherBach;
2822 
2823  else return -1;
2824 
2825 }
2826 
2827 
2828 // TProfile* AliAnalysisTaskSELc2V0bachelorTMVAApp::GetEstimatorHistogram(const AliVEvent* event){
2829 // /// Get estimator histogram from period based on event->GetRunNumber();
2830 
2831 // Int_t runNo = event->GetRunNumber();
2832 // Int_t period = -1;
2833 // switch (fAnalysisType) { // flag to set which system and year is being used
2834 // case kpPb2013: //0 = LHC13b, 1 = LHC13c
2835 // if (runNo > 195343 && runNo < 195484) period = 0;
2836 // if (runNo > 195528 && runNo < 195678) period = 1;
2837 // if (period < 0 || period > 1) { AliInfo(Form("Run number %d not found for LHC13!",runNo)); return 0;}
2838 // break;
2839 // case kpPb2016: //0 = LHC16q, 265499 -- 265525 || 265309 -- 265387, 1 = LHC16q, 265435, 2 = LHC16q, 265388 -- 265427, 3 = LHC16t, 267163 -- 267166
2840 // if ((runNo >=265499 && runNo <=265525) || (runNo >= 265309 && runNo <= 265387)) period = 0;
2841 // else if (runNo == 265435) period = 1;
2842 // else if (runNo >= 265388 && runNo <= 265427) period = 2;
2843 // else if (runNo >=267163 && runNo <=276166) period = 3;
2844 // if (period < 0 || period > 3) { AliInfo(Form("Run number %d not found for LHC16 pPb!",runNo)); return 0;}
2845 // break;
2846 // case kpp2016: //0 = LHC16j, 1 = LHC16k, 2 = LHC16l
2847 // if (runNo >= 256219 && runNo <= 256418) period = 0;
2848 // else if (runNo >= 256504 && runNo <= 258537) period = 1;
2849 // else if (runNo >= 258883 && runNo <= 260187) period = 2;
2850 // if (period < 0 || period > 2) {AliInfo(Form("Run number %d not found for LHC16 pp!",runNo)); return 0;}
2851 // break;
2852 // case kpp2010: // 0 = LHC10b, 1 = LHC10c, 2 = LHC10d, 3 = LHC10e
2853 // if (runNo > 114930 && runNo < 117223) period = 0;
2854 // if (runNo > 119158 && runNo < 120830) period = 1;
2855 // if (runNo > 122373 && runNo < 126438) period = 2;
2856 // if (runNo > 127711 && runNo < 130851) period = 3;
2857 // if (period < 0 || period > 3) {AliInfo(Form("Run number %d not found for LHC10 pp!",runNo)); return 0;}
2858 // break;
2859 // default: //no valid switch
2860 // return 0;
2861 // break;
2862 // }
2863 
2864 // return fMultEstimatorAvg[period];
2865 // }
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 pdg
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]
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
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
char Char_t
Definition: External.C:18
Double_t InvMassLctoLambdaPi() const
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
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
TH1D * fHistoMassV0All
! KF: mass for all V0 reconstructed with KF
TH1F * fCEvents
! Histogram to check selected events
TH2D * fHistoMassKFV0
! KF: mass vs mass error for V0 from KF
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:271
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
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
TH2D * fBDTHistoVsCosThetaStar
! BDT classifier vs proton emission angle in pK0s pair rest frame
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]
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
Float_t GetCentrality(AliAODEvent *aodEvent)
Definition: AliRDHFCuts.h:272
TH1D * fHistoMassV0TrueK0SFromAOD
! KF: AOD mass for true V0 which are really K0S reconstructed with KF
TH1D * fHistoMassV0fromLcSgn
! KF: mass of V0 for signal Lc reconstructed with KF
Bool_t HasCascadeCandidateAnyDaughInjected(AliAODRecoCascadeHF *cand, AliAODMCHeader *header, TClonesArray *arrayMC)
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.
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
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
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.)
AliVertexingHFUtils * fUtils
flag to fill bkg with only candidates that have daughters generated by HIJING (to be used for enriche...
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
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
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
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:207
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
Double_t DecayLength() const
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)
void SetTriggerMask(ULong64_t mask=0)
Definition: AliRDHFCuts.h:70
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
flag to decide whether to fill the sgn and bkg trees
TH1D * fHistoDecayLengthLcAll
! KF: decay length for all Lc reconstructed with KF
Class with functions useful for different D2H analyses //.