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