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