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