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