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