AliPhysics  9b6b435 (9b6b435)
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 
572  fCandidateVariables = new Float_t [nVar];
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 
1538  }
1539 
1540  return;
1541 
1542 }
1543 //________________________________________________________________________
1545  Int_t isLc,
1546  Int_t &nSelectedAnal,
1547  AliRDHFCutsLctoV0 *cutsAnal,
1548  TClonesArray *mcArray, Int_t iLctopK0s){
1549  //
1551  //
1552 
1553  /*
1554  if (!part->GetOwnPrimaryVtx()) {
1555  //Printf("No primary vertex for Lc found!!");
1556  part->SetOwnPrimaryVtx(fVtx1);
1557  }
1558  else {
1559  //Printf("Yu-huuuu!!! primary vertex for Lc found!!");
1560  }
1561  */
1562  Double_t invmassLc = part->InvMassLctoK0sP();
1563 
1564  AliAODv0 * v0part = part->Getv0();
1565  Bool_t onFlyV0 = v0part->GetOnFlyStatus(); // on-the-flight V0s
1566  if (onFlyV0){ // on-the-fly V0
1567  if (isLc) { // Lc
1568  fHistoLcOnTheFly->Fill(2.);
1569  }
1570  else { // not Lc
1571  fHistoLcOnTheFly->Fill(0.);
1572  }
1573  }
1574  else { // offline V0
1575  if (isLc) { // Lc
1576  fHistoLcOnTheFly->Fill(3.);
1577  }
1578  else { // not Lc
1579  fHistoLcOnTheFly->Fill(1.);
1580  }
1581  }
1582 
1583  Double_t dcaV0 = v0part->GetDCA();
1584  Double_t invmassK0s = v0part->MassK0Short();
1585 
1586  if ( (cutsAnal->IsInFiducialAcceptance(part->Pt(),part->Y(4122))) ) {
1587  if (isLc) {
1588  fHistoFiducialAcceptance->Fill(3.);
1589  }
1590  else {
1591  fHistoFiducialAcceptance->Fill(1.);
1592  }
1593  }
1594  else {
1595  if (isLc) {
1596  fHistoFiducialAcceptance->Fill(2.);
1597  }
1598  else {
1599  fHistoFiducialAcceptance->Fill(0.);
1600  }
1601  }
1602 
1603  Int_t isInV0window = (((cutsAnal->IsSelectedSingleCut(part, AliRDHFCuts::kCandidate, 2)) & (AliRDHFCutsLctoV0::kLcToK0Spr)) == (AliRDHFCutsLctoV0::kLcToK0Spr)); // cut on V0 invMass
1604 
1605  if (isInV0window == 0) {
1606  AliDebug(2, "No: The candidate has NOT passed the V0 window cuts!");
1607  if (isLc) Printf("SIGNAL candidate rejected: V0 window cuts");
1608  return;
1609  }
1610  else AliDebug(2, "Yes: The candidate has passed the mass cuts!");
1611 
1612  Bool_t isInCascadeWindow = (((cutsAnal->IsSelectedSingleCut(part, AliRDHFCuts::kCandidate, 0)) & (AliRDHFCutsLctoV0::kLcToK0Spr)) == (AliRDHFCutsLctoV0::kLcToK0Spr)); // cut on Lc->p+K0S invMass
1613 
1614  if (!isInCascadeWindow) {
1615  AliDebug(2, "No: The candidate has NOT passed the cascade window cuts!");
1616  if (isLc) Printf("SIGNAL candidate rejected: cascade window cuts");
1617  return;
1618  }
1619  else AliDebug(2, "Yes: The candidate has passed the cascade window cuts!");
1620 
1621  Bool_t isCandidateSelectedCuts = (((cutsAnal->IsSelected(part, AliRDHFCuts::kCandidate)) & (AliRDHFCutsLctoV0::kLcToK0Spr)) == (AliRDHFCutsLctoV0::kLcToK0Spr)); // kinematic/topological cuts
1622  AliDebug(2, Form("recoAnalysisCuts = %d", cutsAnal->IsSelected(part, AliRDHFCuts::kCandidate) & (AliRDHFCutsLctoV0::kLcToK0Spr)));
1623  if (!isCandidateSelectedCuts){
1624  AliDebug(2, "No: Analysis cuts kCandidate level NOT passed");
1625  if (isLc) Printf("SIGNAL candidate rejected");
1626  return;
1627  }
1628  else {
1629  AliDebug(2, "Yes: Analysis cuts kCandidate level passed");
1630  }
1631 
1632  AliAODTrack *bachelor = (AliAODTrack*)part->GetBachelor();
1633  if (!bachelor) {
1634  AliDebug(2, Form("Very weird, the bachelor is not there... returning for this candidate"));
1635  return;
1636  }
1637 
1638  //Bool_t isBachelorID = (((cutsAnal->IsSelected(part,AliRDHFCuts::kPID))&(AliRDHFCutsLctoV0::kLcToK0Spr))==(AliRDHFCutsLctoV0::kLcToK0Spr)); // ID x bachelor
1639  Double_t probTPCTOF[AliPID::kSPECIES]={-1.};
1640 
1641  UInt_t detUsed = fPIDCombined->ComputeProbabilities(bachelor, fPIDResponse, probTPCTOF);
1642  AliDebug(2, Form("detUsed (TPCTOF case) = %d", detUsed));
1643  Double_t probProton = -1.;
1644  // Double_t probPion = -1.;
1645  // Double_t probKaon = -1.;
1646  if (detUsed == (UInt_t)fPIDCombined->GetDetectorMask() ) {
1647  AliDebug(2, Form("We have found the detector mask for TOF + TPC: probProton will be set to %f", probTPCTOF[AliPID::kProton]));
1648  probProton = probTPCTOF[AliPID::kProton];
1649  // probPion = probTPCTOF[AliPID::kPion];
1650  // probKaon = probTPCTOF[AliPID::kKaon];
1651  }
1652  else { // if you don't have both TOF and TPC, try only TPC
1653  fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
1654  AliDebug(2, "We did not find the detector mask for TOF + TPC, let's see only TPC");
1655  detUsed = fPIDCombined->ComputeProbabilities(bachelor, fPIDResponse, probTPCTOF);
1656  AliDebug(2,Form(" detUsed (TPC case) = %d", detUsed));
1657  if (detUsed == (UInt_t)fPIDCombined->GetDetectorMask()) {
1658  probProton = probTPCTOF[AliPID::kProton];
1659  // probPion = probTPCTOF[AliPID::kPion];
1660  // probKaon = probTPCTOF[AliPID::kKaon];
1661  AliDebug(2, Form("TPC only worked: probProton will be set to %f", probTPCTOF[AliPID::kProton]));
1662  }
1663  else {
1664  AliDebug(2, "Only TPC did not work...");
1665  }
1666  // resetting mask to ask for both TPC+TOF
1667  fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC+AliPIDResponse::kDetTOF);
1668  }
1669  AliDebug(2, Form("probProton = %f", probProton));
1670 
1671  // now we get the TPC and TOF single PID probabilities (only for Proton, or the tree will explode :) )
1672  Double_t probProtonTPC = -1.;
1673  Double_t probProtonTOF = -1.;
1674  Double_t pidTPC[AliPID::kSPECIES]={-1.};
1675  Double_t pidTOF[AliPID::kSPECIES]={-1.};
1676  Int_t respTPC = fPIDResponse->ComputePIDProbability(AliPIDResponse::kDetTPC, bachelor, AliPID::kSPECIES, pidTPC);
1677  Int_t respTOF = fPIDResponse->ComputePIDProbability(AliPIDResponse::kDetTOF, bachelor, AliPID::kSPECIES, pidTOF);
1678  if (respTPC == AliPIDResponse::kDetPidOk) probProtonTPC = pidTPC[AliPID::kProton];
1679  if (respTOF == AliPIDResponse::kDetPidOk) probProtonTOF = pidTOF[AliPID::kProton];
1680 
1681  // checking V0 status (on-the-fly vs offline)
1682  if ( !( !onFlyV0 || (onFlyV0 && fUseOnTheFlyV0) ) ) {
1683  AliDebug(2, "On-the-fly discarded");
1684  return;
1685  }
1686 
1688  nSelectedAnal++;
1689  }
1690 
1691  if ( !(cutsAnal->IsInFiducialAcceptance(part->Pt(),part->Y(4122))) ) return;
1692 
1693  if ( !( ( (cutsAnal->IsSelected(part, AliRDHFCuts::kTracks)) & (AliRDHFCutsLctoV0::kLcToK0Spr)) == (AliRDHFCutsLctoV0::kLcToK0Spr) ) ) { // esd track cuts
1694  if (isLc) Printf("SIGNAL candidate rejected");
1695  AliDebug(2, "No: Analysis cuts kTracks level NOT passed");
1696  return;
1697  }
1698  else {
1699  AliDebug(2, "Yes: Analysis cuts kTracks level passed");
1700  }
1701 
1702  Int_t pdgCand = 4122;
1703  Int_t pdgDgLctoV0bachelor[2]={211, 3122}; // case of wrong decay! Lc --> L + pi
1704  Int_t pdgDgV0toDaughters[2]={2212, 211}; // case of wrong decay! Lc --> L + pi
1705  Int_t isLc2LBarpi=0, isLc2Lpi=0;
1706  Int_t currentLabel = part->GetLabel();
1707  Int_t mcLabel = 0;
1708  if (fUseMCInfo) {
1709  mcLabel = part->MatchToMC(pdgCand, pdgDgLctoV0bachelor[1], pdgDgLctoV0bachelor, pdgDgV0toDaughters, mcArray, kTRUE);
1710  if (mcLabel>=0) {
1711  if (bachelor->Charge()==-1) isLc2LBarpi=1;
1712  if (bachelor->Charge()==+1) isLc2Lpi=1;
1713  }
1714  }
1715 
1716  Int_t pdgDg2prong[2] = {211, 211};
1717  Int_t labelK0S = 0;
1718  Int_t isK0S = 0;
1719  if (fUseMCInfo) {
1720  labelK0S = v0part->MatchToMC(310, mcArray, 2, pdgDg2prong);
1721  if (labelK0S>=0) isK0S = 1;
1722  }
1723 
1724  pdgDg2prong[0] = 211;
1725  pdgDg2prong[1] = 2212;
1726  Int_t isLambda = 0;
1727  Int_t isLambdaBar = 0;
1728  Int_t lambdaLabel = 0;
1729  if (fUseMCInfo) {
1730  lambdaLabel = v0part->MatchToMC(3122, mcArray, 2, pdgDg2prong);
1731  if (lambdaLabel>=0) {
1732  AliAODMCParticle *lambdaTrack = (AliAODMCParticle*)mcArray->At(lambdaLabel);
1733  if (lambdaTrack->GetPdgCode()==3122) isLambda = 1;
1734  else if (lambdaTrack->GetPdgCode()==-3122) isLambdaBar = 1;
1735  }
1736  }
1737 
1738  pdgDg2prong[0] = 11;
1739  pdgDg2prong[1] = 11;
1740  Int_t isGamma = 0;
1741  Int_t gammaLabel = 0;
1742  if (fUseMCInfo) {
1743  gammaLabel = v0part->MatchToMC(22, mcArray, 2, pdgDg2prong);
1744  if (gammaLabel>=0) {
1745  AliAODMCParticle *gammaTrack = (AliAODMCParticle*)mcArray->At(gammaLabel);
1746  if (gammaTrack->GetPdgCode()==22) isGamma = 1;
1747  }
1748  }
1749 
1750  Int_t pdgTemp = -1;
1751  if (currentLabel != -1){
1752  AliAODMCParticle *tempPart = (AliAODMCParticle*)mcArray->At(currentLabel);
1753  pdgTemp = tempPart->GetPdgCode();
1754  }
1755  if (isLc) AliDebug(2, Form("Signal: Candidate is a Lc in K0s+p"));
1756  else if (isLc2LBarpi) AliDebug(2, Form("Background: Candidate is a Lc in Lbar + pi"));
1757  else if (isLc2Lpi) AliDebug(2, Form("Background: Candidate is a Lc in L + pi"));
1758  else AliDebug(2, Form("Pure bkg: Candidate has pdg = %d", pdgTemp));
1759  if (isK0S) AliDebug(2, Form("V0 is a K0S"));
1760  else if (isLambda) AliDebug(2, Form("V0 is a Lambda"));
1761  else if (isLambdaBar) AliDebug(2, Form("V0 is a LambdaBar"));
1762  else if (isGamma) AliDebug(2, Form("V0 is a Gamma"));
1763  //else AliDebug(2, Form("V0 is something else!!"));
1764 
1765  Double_t invmassLc2Lpi = part->InvMassLctoLambdaPi();
1766  Double_t invmassLambda = v0part->MassLambda();
1767  Double_t invmassLambdaBar = v0part->MassAntiLambda();
1768 
1769  //Double_t nSigmaITSpr=-999.;
1770  Double_t nSigmaTPCpr=-999.;
1771  Double_t nSigmaTOFpr=-999.;
1772 
1773  //Double_t nSigmaITSpi=-999.;
1774  Double_t nSigmaTPCpi=-999.;
1775  Double_t nSigmaTOFpi=-999.;
1776 
1777  //Double_t nSigmaITSka=-999.;
1778  Double_t nSigmaTPCka=-999.;
1779  Double_t nSigmaTOFka=-999.;
1780 
1781  /*
1782  cutsAnal->GetPidHF()->GetnSigmaITS(bachelor,4,nSigmaITSpr);
1783  cutsAnal->GetPidHF()->GetnSigmaTPC(bachelor,4,nSigmaTPCpr);
1784  cutsAnal->GetPidHF()->GetnSigmaTOF(bachelor,4,nSigmaTOFpr);
1785  cutsAnal->GetPidHF()->GetnSigmaITS(bachelor,2,nSigmaITSpi);
1786  cutsAnal->GetPidHF()->GetnSigmaTPC(bachelor,2,nSigmaTPCpi);
1787  cutsAnal->GetPidHF()->GetnSigmaTOF(bachelor,2,nSigmaTOFpi);
1788  cutsAnal->GetPidHF()->GetnSigmaITS(bachelor,3,nSigmaITSka);
1789  cutsAnal->GetPidHF()->GetnSigmaTPC(bachelor,3,nSigmaTPCka);
1790  cutsAnal->GetPidHF()->GetnSigmaTOF(bachelor,3,nSigmaTOFka);
1791  */
1792 
1793  nSigmaTPCpi = fPIDResponse->NumberOfSigmasTPC(bachelor,(AliPID::kPion));
1794  nSigmaTPCka = fPIDResponse->NumberOfSigmasTPC(bachelor,(AliPID::kKaon));
1795  nSigmaTPCpr = fPIDResponse->NumberOfSigmasTPC(bachelor,(AliPID::kProton));
1796  nSigmaTOFpi = fPIDResponse->NumberOfSigmasTOF(bachelor,(AliPID::kPion));
1797  nSigmaTOFka = fPIDResponse->NumberOfSigmasTOF(bachelor,(AliPID::kKaon));
1798  nSigmaTOFpr = fPIDResponse->NumberOfSigmasTOF(bachelor,(AliPID::kProton));
1799 
1800  Double_t ptLcMC = -1;
1801  Double_t weightPythia = -1, weight5LHC13d3 = -1, weight5LHC13d3Lc = -1;
1802 
1803  if (fUseMCInfo) {
1804  if (iLctopK0s >= 0) {
1805  AliAODMCParticle *partLcMC = (AliAODMCParticle*)mcArray->At(iLctopK0s);
1806  ptLcMC = partLcMC->Pt();
1807  //Printf("--------------------- Reco pt = %f, MC particle pt = %f", part->Pt(), ptLcMC);
1808  weightPythia = fFuncWeightPythia->Eval(ptLcMC);
1809  weight5LHC13d3 = fFuncWeightFONLL5overLHC13d3->Eval(ptLcMC);
1810  weight5LHC13d3Lc = fFuncWeightFONLL5overLHC13d3Lc->Eval(ptLcMC);
1811  }
1812  }
1813 
1814  Double_t weightNch = 1;
1815  if (fUseMCInfo) {
1816  //Int_t nChargedMCPhysicalPrimary=AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(mcArray,-1.0,1.0);
1817  // if(nChargedMCPhysicalPrimary > 0)
1818 
1819  if(fNTracklets > 0){
1820  if(!fHistoMCNch) AliInfo("Input histos to evaluate Nch weights missing");
1821  if(fHistoMCNch) weightNch *= fHistoMCNch->GetBinContent(fHistoMCNch->FindBin(fNTracklets));
1822  }
1823  }
1824 
1825 
1826 
1827 
1828  // Fill candidate variable Tree (track selection, V0 invMass selection)
1829  if (!onFlyV0 && isInV0window && isInCascadeWindow && part->CosV0PointingAngle()>0.99 && TMath::Abs(nSigmaTPCpr) <= 3 && v0part->Getd0Prong(0) < 20 && v0part->Getd0Prong(1) < 20) {
1830 
1831  EBachelor bachCode = kBachInvalid;
1832  EK0S k0SCode = kK0SInvalid;
1833  if (fUseMCInfo) {
1834  bachCode = CheckBachelor(part, bachelor, mcArray);
1835  k0SCode = CheckK0S(part, v0part, mcArray);
1836  }
1837 
1838  Double_t V0KF[3] = {-999999, -999999, -999999};
1839  Double_t errV0KF[3] = {-999999, -999999, -999999};
1840  Double_t LcKF[3] = {-999999, -999999, -999999};
1841  Double_t errLcKF[3] = {-999999, -999999, -999999};
1842  Double_t distances[3] = {-999999, -999999, -999999};
1843  Double_t armPolKF[2] = {-999999, -999999};
1844 
1845  AliAODTrack *v0pos = (AliAODTrack*)part->Getv0PositiveTrack();
1846  AliAODTrack *v0neg = (AliAODTrack*)part->Getv0NegativeTrack();
1847  Int_t kfResult;
1848  TVector3 mom1(bachelor->Px(), bachelor->Py(), bachelor->Pz());
1849  TVector3 mom2(v0part->Px(), v0part->Py(), v0part->Pz());
1850  TVector3 momTot(part->Px(), part->Py(), part->Pz());
1851 
1852  Double_t Ql1 = mom1.Dot(momTot)/momTot.Mag();
1853  Double_t Ql2 = mom2.Dot(momTot)/momTot.Mag();
1854 
1855  Double_t alphaArmLc = (Ql1 - Ql2)/(Ql1 + Ql2);
1856  Double_t alphaArmLcCharge = ( bachelor->Charge() > 0 ? (Ql1 - Ql2)/(Ql1 + Ql2) : (Ql2 - Ql1)/(Ql1 + Ql2) );
1857  Double_t ptArmLc = mom1.Perp(momTot);
1858 
1859  Double_t massK0SPDG = TDatabasePDG::Instance()->GetParticle(310)->Mass(); // mass K0S PDG
1860  Double_t massPrPDG = TDatabasePDG::Instance()->GetParticle(2212)->Mass(); // mass Proton PDG
1861  Double_t massLcPDG = TDatabasePDG::Instance()->GetParticle(4122)->Mass(); // mass Lc PDG
1862 
1863  Double_t pStar = TMath::Sqrt((massLcPDG*massLcPDG-massPrPDG*massPrPDG-massK0SPDG*massK0SPDG)*(massLcPDG*massLcPDG-massPrPDG*massPrPDG-massK0SPDG*massK0SPDG)-4.*massPrPDG*massPrPDG*massK0SPDG*massK0SPDG)/(2.*massLcPDG);
1864  Double_t e = part->E(4122);
1865  Double_t beta = part->P()/e;
1866  Double_t gamma = e/massLcPDG;
1867 
1868  if (fCallKFVertexing){
1869  kfResult = CallKFVertexing(part, v0part, bachelor, mcArray, &V0KF[0], &errV0KF[0], &LcKF[0], &errLcKF[0], &distances[0], &armPolKF[0]);
1870  AliDebug(2, Form("Result from KF = %d", kfResult));
1871  }
1872 
1873  /*
1874  for (Int_t i = 0; i< 3; i++){
1875  Printf("i = %d, V0KF = %f, errV0KF = %f, LcKF = %f, errLcKF = %f", V0KF[i], errV0KF[i], LcKF[i], errLcKF[i]);
1876  }
1877  */
1878  Double_t cts = (Ql1/gamma-beta*TMath::Sqrt(pStar*pStar+massPrPDG*massPrPDG))/pStar;
1879 
1880  Double_t countTreta1corr = fNTracklets;
1881 
1882  AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
1883  TProfile* estimatorAvg = GetEstimatorHistogram(aodEvent);
1884  if(estimatorAvg) {
1885  countTreta1corr = static_cast<Int_t>(AliVertexingHFUtils::GetCorrectedNtracklets(estimatorAvg,fNTracklets,fVtx1->GetZ(),fRefMult));
1886  }
1887 
1888 
1889 
1890 
1891  if (fUseMCInfo) { // save full tree if on MC
1892  fCandidateVariables[0] = invmassLc;
1893  fCandidateVariables[1] = invmassLc2Lpi;
1894  fCandidateVariables[2] = invmassK0s;
1895  fCandidateVariables[3] = invmassLambda;
1896  fCandidateVariables[4] = invmassLambdaBar;
1898  fCandidateVariables[6] = dcaV0;
1899  fCandidateVariables[7] = part->Getd0Prong(0);
1900  fCandidateVariables[8] = part->Getd0Prong(1);
1901  fCandidateVariables[9] = nSigmaTPCpr;
1902 // fCandidateVariables[10] = nSigmaTPCpi;
1903 // fCandidateVariables[11] = nSigmaTPCka;
1904  fCandidateVariables[10] = nSigmaTOFpr;
1905 // fCandidateVariables[14] = nSigmaTOFka;
1906 // fCandidateVariables[13] = nSigmaTOFpi;
1907  fCandidateVariables[11] = bachelor->Pt();
1908  fCandidateVariables[12] = v0pos->Pt();
1909  fCandidateVariables[13] = v0neg->Pt();
1910  fCandidateVariables[14] = v0part->Getd0Prong(0);
1911  fCandidateVariables[15] = v0part->Getd0Prong(1);
1912  fCandidateVariables[16] = v0part->Pt();
1913  fCandidateVariables[17] = v0part->InvMass2Prongs(0,1,11,11);
1914  fCandidateVariables[18] = part->Pt();
1915  fCandidateVariables[19] = probProton;
1916  fCandidateVariables[20] = part->Eta();
1917  fCandidateVariables[21] = v0pos->Eta();
1918  fCandidateVariables[22] = v0neg->Eta();
1919  fCandidateVariables[23] = probProtonTPC;
1920  fCandidateVariables[24] = probProtonTOF;
1921  fCandidateVariables[25] = bachelor->Eta();
1922 
1923  fCandidateVariables[26] = part->P();
1924  fCandidateVariables[27] = bachelor->P();
1925  fCandidateVariables[28] = v0part->P();
1926  fCandidateVariables[29] = v0pos->P();
1927  fCandidateVariables[30] = v0neg->P();
1928 
1929 // fCandidateVariables[35] = part->Y(4122);
1930 // fCandidateVariables[36] = bachelor->Y(2212);
1931 // fCandidateVariables[37] = v0part->Y(310);
1932 // fCandidateVariables[38] = v0pos->Y(211);
1933 // fCandidateVariables[39] = v0neg->Y(211);
1934 
1935  fCandidateVariables[31] = v0part->Eta();
1936 
1937  fCandidateVariables[32] = part->DecayLength();
1938  fCandidateVariables[33] = part->DecayLengthV0();
1939  fCandidateVariables[34] = part->Ct(4122);
1940  fCandidateVariables[35] = v0part->Ct(310, v0part->GetSecondaryVtx());
1941 
1942 
1943 
1944  fCandidateVariables[36] = bachCode;
1945  fCandidateVariables[37] = k0SCode;
1946 // fCandidateVariables[47] = V0KF[0];
1947 // fCandidateVariables[48] = V0KF[1];
1948 // fCandidateVariables[49] = V0KF[2];
1949 //
1950 // fCandidateVariables[50] = errV0KF[0];
1951 // fCandidateVariables[51] = errV0KF[1];
1952 // fCandidateVariables[52] = errV0KF[2];
1953 //
1954 // fCandidateVariables[53] = LcKF[0];
1955 // fCandidateVariables[54] = LcKF[1];
1956 // fCandidateVariables[55] = LcKF[2];
1957 //
1958 // fCandidateVariables[56] = errLcKF[0];
1959 // fCandidateVariables[57] = errLcKF[1];
1960 // fCandidateVariables[58] = errLcKF[2];
1961 //
1962 // fCandidateVariables[59] = distances[0];
1963 // fCandidateVariables[60] = distances[1];
1964 // fCandidateVariables[61] = distances[2];
1965 // fCandidateVariables[62] = armPolKF[0];
1966 // fCandidateVariables[63] = armPolKF[1];
1967  fCandidateVariables[38] = v0part->AlphaV0();
1968  fCandidateVariables[39] = v0part->PtArmV0();
1969 
1970  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)));
1971 // fCandidateVariables[66] = v0pos->GetStatus() & AliESDtrack::kITSrefit;
1972 // fCandidateVariables[67] = v0neg->GetStatus() & AliESDtrack::kITSrefit;
1973 // fCandidateVariables[68] = v0pos->GetTPCClusterInfo(2, 1);
1974 // fCandidateVariables[69] = v0neg->GetTPCClusterInfo(2, 1);
1975 //
1976 // fCandidateVariables[70] = v0part->Xv();
1977 // fCandidateVariables[71] = v0part->Yv();
1978 // fCandidateVariables[72] = v0part->Zv();
1979 //
1980 // fCandidateVariables[73] = fVtx1->GetX();
1981 // fCandidateVariables[74] = fVtx1->GetY();
1982 // fCandidateVariables[75] = fVtx1->GetZ();
1983 //
1984 // fCandidateVariables[76] = bachelor->GetITSNcls();
1985 // fCandidateVariables[77] = bachelor->HasPointOnITSLayer(0) + bachelor->HasPointOnITSLayer(1);
1986 //
1987 // fCandidateVariables[78] = v0pos->GetITSNcls();
1988 // fCandidateVariables[79] = v0pos->HasPointOnITSLayer(0) + v0pos->HasPointOnITSLayer(1);
1989 //
1990 // fCandidateVariables[80] = v0neg->GetITSNcls();
1991 // fCandidateVariables[81] = v0neg->HasPointOnITSLayer(0) + v0neg->HasPointOnITSLayer(1);
1992 //
1993 //
1994 // fCandidateVariables[82] = alphaArmLc;
1995 // fCandidateVariables[83] = alphaArmLcCharge;
1996 // fCandidateVariables[84] = ptArmLc;
1997 
1998 
1999  fCandidateVariables[40] = cts;
2000 
2001  fCandidateVariables[41] = weightPythia;
2002  fCandidateVariables[42] = weight5LHC13d3;
2003  fCandidateVariables[43] = weight5LHC13d3Lc;
2004  fCandidateVariables[44] = weightNch;
2006 
2007  fCandidateVariables[46] = countTreta1corr;
2008  }
2009 
2010  else { //remove MC-only variables from tree if data
2011  fCandidateVariables[0] = invmassLc;
2012  fCandidateVariables[1] = invmassLc2Lpi;
2013  fCandidateVariables[2] = invmassK0s;
2014  fCandidateVariables[3] = invmassLambda;
2015  fCandidateVariables[4] = invmassLambdaBar;
2017  fCandidateVariables[6] = dcaV0;
2018  fCandidateVariables[7] = part->Getd0Prong(0);
2019  fCandidateVariables[8] = part->Getd0Prong(1);
2020  fCandidateVariables[9] = nSigmaTPCpr;
2021  fCandidateVariables[10] = nSigmaTOFpr;
2022  fCandidateVariables[11] = bachelor->Pt();
2023  fCandidateVariables[12] = v0pos->Pt();
2024  fCandidateVariables[13] = v0neg->Pt();
2025  fCandidateVariables[14] = v0part->Getd0Prong(0);
2026  fCandidateVariables[15] = v0part->Getd0Prong(1);
2027  fCandidateVariables[16] = v0part->Pt();
2028  fCandidateVariables[17] = v0part->InvMass2Prongs(0,1,11,11);
2029  fCandidateVariables[18] = part->Pt();
2030  fCandidateVariables[19] = probProton;
2031  fCandidateVariables[20] = v0pos->Eta();
2032  fCandidateVariables[21] = v0neg->Eta();
2033  fCandidateVariables[22] = bachelor->Eta();
2034  fCandidateVariables[23] = v0part->P();
2035  fCandidateVariables[24] = part->DecayLengthV0();
2036  fCandidateVariables[25] = v0part->Ct(310, v0part->GetSecondaryVtx());
2037  fCandidateVariables[26] = v0part->AlphaV0();
2038  fCandidateVariables[27] = v0part->PtArmV0();
2040  fCandidateVariables[29] = countTreta1corr;
2041 
2042  }
2043 
2044  // fill multiplicity histograms for events with a candidate
2045  fHistNtrUnCorrEvWithCand->Fill(fNTracklets,weightNch);
2046  fHistNtrCorrEvWithCand->Fill(countTreta1corr,weightNch);
2047  if (fUseMCInfo) {
2048  if (isLc){
2049  AliDebug(2, Form("Reco particle %d --> Filling Sgn", iLctopK0s));
2050  fVariablesTreeSgn->Fill();
2051  fHistoCodesSgn->Fill(bachCode, k0SCode);
2052  }
2053  else {
2054  if (fFillOnlySgn == kFALSE){
2055  AliDebug(2, "Filling Bkg");
2056  fVariablesTreeBkg->Fill();
2057  fHistoCodesBkg->Fill(bachCode, k0SCode);
2058  }
2059  }
2060  }
2061  else {
2062  fVariablesTreeSgn->Fill();
2063  }
2064  }
2065 
2066  return;
2067 
2068 
2069  }
2070 
2071 //________________________________________________________________________
2072 Int_t AliAnalysisTaskSELc2V0bachelorTMVA::CallKFVertexing(AliAODRecoCascadeHF *cascade, AliAODv0* v0part, AliAODTrack* bach, TClonesArray *mcArray,
2073  Double_t* V0KF, Double_t* errV0KF, Double_t* LcKF, Double_t* errLcKF,
2074  Double_t* distances, Double_t* armPolKF) {
2075 
2076  //
2079  //
2080 
2081  Int_t codeKFV0 = -1, codeKFLc = -1;
2082 
2083  AliKFVertex primVtxCopy;
2084  Int_t nt = 0, ntcheck = 0;
2085  Double_t pos[3] = {0., 0., 0.};
2086 
2087  fVtx1->GetXYZ(pos);
2088  Int_t contr = fVtx1->GetNContributors();
2089  Double_t covmatrix[6] = {0.};
2090  fVtx1->GetCovarianceMatrix(covmatrix);
2091  Double_t chi2 = fVtx1->GetChi2();
2092  AliESDVertex primaryESDVtxCopy(pos, covmatrix, chi2, contr, "Vertex");
2093 
2094  // topological constraint
2095  primVtxCopy = AliKFVertex(primaryESDVtxCopy);
2096  nt = primaryESDVtxCopy.GetNContributors();
2097  ntcheck = nt;
2098 
2099  Int_t pdg[2] = {211, -211};
2100  Int_t pdgLc[2] = {2212, 310};
2101 
2102  Int_t pdgDgV0toDaughters[2] = {211, 211};
2103 
2104  Int_t mcLabelV0 = v0part->MatchToMC(310, mcArray, 2, pdgDgV0toDaughters);
2105 
2106  // the KF vertex for the V0 has to be built with the prongs of the V0!
2107  Bool_t isMCokV0 = kTRUE, isBkgV0 = kFALSE;
2108  AliKFParticle V0, positiveV0KF, negativeV0KF;
2109  Int_t labelsv0daugh[2] = {-1, -1};
2110  Int_t idv0daugh[2] = {-1, -1};
2111  AliExternalTrackParam* esdv0Daugh1 = 0x0;
2112  AliExternalTrackParam* esdv0Daugh2 = 0x0;
2113  for(Int_t ipr= 0; ipr < 2; ipr++){ // 0 is positive, 1 is negative
2114  AliAODTrack *aodTrack = (AliAODTrack*)v0part->GetDaughter(ipr);
2115  if(!aodTrack) {
2116  AliDebug(2, "No V0 daughters available");
2117  return -1;
2118  }
2119  Double_t xyz[3], pxpypz[3], cv[21];
2120  Short_t sign;
2121  aodTrack->GetXYZ(xyz);
2122  aodTrack->PxPyPz(pxpypz);
2123  aodTrack->GetCovarianceXYZPxPyPz(cv);
2124  sign = aodTrack->Charge();
2125  AliExternalTrackParam tmp1( xyz, pxpypz, cv, sign);
2126 
2127  if (ipr == 0) esdv0Daugh1 = new AliExternalTrackParam( xyz, pxpypz, cv, sign);
2128  else esdv0Daugh2 = new AliExternalTrackParam( xyz, pxpypz, cv, sign);
2129  labelsv0daugh[ipr] = TMath::Abs(aodTrack->GetLabel());
2130  idv0daugh[ipr] = aodTrack->GetID();
2131  if (labelsv0daugh[ipr] == -1) isBkgV0 = kTRUE;
2132 
2133  //Printf("v0 daughter %d has label %d", ipr, labelsv0daugh[ipr]);
2134 
2135  AliKFParticle daughterKF(*aodTrack, pdg[ipr]); // we assume that the PDG is correct
2136  if (aodTrack->Charge() > 0) { // assigning positive and negative track to KF V0 for Armenteros-Podolanski plot
2137  positiveV0KF = daughterKF;
2138  }
2139  else {
2140  negativeV0KF = daughterKF;
2141  }
2142  }
2143 
2144  Double_t xn=0., xp=0.;//, dca;
2145  AliDebug(2, Form("bField = %f, esdv0Daugh1 = %p, esdv0Daugh2 = %p", fBField, esdv0Daugh1, esdv0Daugh2));
2146  // dca = esdv0Daugh1->GetDCA(esdv0Daugh2, fBField, xn, xp);
2147 
2148  AliExternalTrackParam tr1(*esdv0Daugh1);
2149  AliExternalTrackParam tr2(*esdv0Daugh2);
2150  tr1.PropagateTo(xn, fBField);
2151  tr2.PropagateTo(xp, fBField);
2152 
2153  AliKFParticle daughterKF1(tr1, 211);
2154  AliKFParticle daughterKF2(tr2, 211);
2155  V0.AddDaughter(positiveV0KF);
2156  V0.AddDaughter(negativeV0KF);
2157  //V0.AddDaughter(daughterKF1);
2158  //V0.AddDaughter(daughterKF2);
2159 
2160  delete esdv0Daugh1;
2161  delete esdv0Daugh2;
2162  esdv0Daugh1=0;
2163  esdv0Daugh2=0;
2164  // Checking the quality of the KF V0 vertex
2165  if( V0.GetNDF() < 1 ) {
2166  //Printf("Number of degrees of freedom < 1, continuing");
2167  return -1;
2168  }
2169  if( TMath::Sqrt(TMath::Abs(V0.GetChi2()/V0.GetNDF())) > fCutKFChi2NDF ) {
2170  //Printf("Chi2 per DOF too big, continuing");
2171  return -1;
2172  }
2173 
2174  if(ftopoConstraint && nt > 0){
2175  for(Int_t ipr = 0; ipr < 2; ipr++){ // 0 is positive, 1 is negative
2176  AliAODTrack *aodTrack = (AliAODTrack*)v0part->GetDaughter(ipr);
2177  //* subtruct daughters from primary vertex
2178  if(!aodTrack->GetUsedForPrimVtxFit()) {
2179  //Printf("Track %d was not used for primary vertex, continuing", i);
2180  continue;
2181  }
2182  AliKFParticle daughterKF(*aodTrack, pdg[ipr]); // we assume that the PDG is correct
2183  primVtxCopy -= daughterKF;
2184  ntcheck--;
2185  }
2186  }
2187 
2188  // Check V0 Chi^2 deviation from primary vertex // not needed for V0 for Lc decay!!
2189  /*
2190  if( V0.GetDeviationFromVertex( primVtxCopy ) < fCutKFDeviationFromVtxV0) {
2191  //Printf("Deviation from vertex too big, continuing");
2192  return -1;
2193  }
2194  */
2195 
2196  //* Get V0 invariant mass
2197  Double_t massV0 = 999999, sigmaMassV0 = 999999;
2198  Int_t retMV0 = V0.GetMass( massV0, sigmaMassV0 );
2199  if( retMV0 ) {
2200  if (massV0 < 0) {
2201  codeKFV0 = 1; // Mass not ok
2202  if (sigmaMassV0 > 1e19) codeKFV0 = 5; // Mass and SigmaMass not ok
2203  }
2204  else if (sigmaMassV0 > 1e19) codeKFV0 = 2; // SigmaMass not ok
2205  }
2206  fHistoMassKFV0->Fill(massV0, sigmaMassV0);
2207 
2208  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]);
2209  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]);
2210 
2211  Printf("Vertices: KF: x = %f, y = %f, z = %f", V0.GetX(), V0.GetY(), V0.GetZ());
2212  Printf("Vertices: AOD: x = %f, y = %f, z = %f", v0part->Xv(), v0part->Yv(), v0part->Zv());
2213 
2214  //Printf("Got MC vtx for V0");
2215  if (fUseMCInfo && TMath::Abs(labelsv0daugh[0] - labelsv0daugh[1]) == 1) {
2216  AliAODMCParticle* tmpdaughv01 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelsv0daugh[0]));
2217  AliAODMCParticle* tmpdaughv02 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelsv0daugh[1]));
2218  if (!tmpdaughv01 && labelsv0daugh[0] > 0){
2219  AliDebug(2, "Could not access MC info for first daughter of V0, continuing");
2220  }
2221  if (!tmpdaughv02 && labelsv0daugh[1] > 0){
2222  AliDebug(2, "Could not access MC info for second daughter of V0, continuing");
2223  }
2224  if(tmpdaughv01){
2225  Double_t xPionMC = tmpdaughv01->Xv(); //Production vertex of Pion --> Where K0S decays
2226  Double_t yPionMC = tmpdaughv01->Yv();
2227  Double_t zPionMC = tmpdaughv01->Zv();
2228  //Printf("Got MC vtx for Pion");
2229  Printf("Vertices: MC: x = %f, y = %f, z = %f", xPionMC, yPionMC, zPionMC);
2230  }
2231  }
2232  else {
2233  Printf("Not a true V0");
2234  }
2235  //massV0=-1;//return -1;// !!!!
2236 
2237  // now use what we just try with the bachelor, to build the Lc
2238 
2239  // topological constraint
2240  nt = primVtxCopy.GetNContributors();
2241  ntcheck = nt;
2242 
2243  Bool_t isMCokLc = kTRUE, isBkgLc = kFALSE;
2244  AliKFParticle Lc;
2245  Int_t labelsLcdaugh[2] = {-1, -1};
2246  labelsLcdaugh[0] = TMath::Abs(bach->GetLabel());
2247  labelsLcdaugh[1] = mcLabelV0;
2248 
2249  if (bach->Charge() < 0) pdgLc[0] = -pdgLc[0];
2250  AliKFParticle daughterKFLc(*bach, pdgLc[0]);
2251  Lc.AddDaughter(daughterKFLc);
2252  TParticlePDG* particlePDG = TDatabasePDG::Instance()->GetParticle(310);
2253  Double_t massPDGK0S = particlePDG->Mass();
2254  V0.SetMassConstraint(massPDGK0S);
2255  Lc.AddDaughter(V0);
2256  if( Lc.GetNDF() < 1 ) {
2257  AliDebug(2, Form("Lc: Number of degrees of freedom < 1 (%d), continuing", Lc.GetNDF()));
2258  return -1;
2259  }
2260  if( TMath::Sqrt(TMath::Abs(Lc.GetChi2()/Lc.GetNDF())) > fCutKFChi2NDF) {
2261  AliDebug(2, Form("Lc: Chi2 per DOF too big, continuing (%f)", TMath::Sqrt(TMath::Abs(Lc.GetChi2()/Lc.GetNDF()))));
2262  return -1;
2263  }
2264 
2265  if(ftopoConstraint && nt > 0){
2266  //* subtruct daughters from primary vertex
2267  if(!bach->GetUsedForPrimVtxFit()) {
2268  AliDebug(3, "Lc: Bachelor was not used for primary vertex, not subtracting it from primary vertex");
2269  }
2270  else{
2271  primVtxCopy -= daughterKFLc;
2272  ntcheck--;
2273  }
2274  /* the V0 was added above, so it is ok to remove it without checking
2275  if(!V0->GetUsedForPrimVtxFit()) {
2276  Printf("Lc: V0 was not used for primary vertex, continuing");
2277  continue;
2278  }
2279  */
2280  //primVtxCopy -= V0;
2281  //ntcheck--;
2282  }
2283 
2284  // Check Lc Chi^2 deviation from primary vertex
2285  /*
2286  if( Lc.GetDeviationFromVertex( primVtxCopy ) > fCutKFDeviationFromVtx) {
2287  AliDebug(2, Form("Lc: Deviation from vertex too big, continuing (%f)", Lc.GetDeviationFromVertex( primVtxCopy )));
2288  return -1;
2289  }
2290 
2291  if(ftopoConstraint){
2292  if(ntcheck>0) {
2293  // Add Lc to primary vertex to improve the primary vertex resolution
2294  primVtxCopy += Lc;
2295  Lc.SetProductionVertex(primVtxCopy);
2296  }
2297  }
2298  */
2299  //* Check chi^2
2300  if( TMath::Sqrt( TMath::Abs(Lc.GetChi2()/Lc.GetNDF())) > fCutKFChi2NDF) {
2301  AliDebug(2, Form("Lc: Final Chi2 per DOF too big, continuing (%f)", TMath::Sqrt( TMath::Abs(Lc.GetChi2()/Lc.GetNDF()))));
2302  return -1;
2303  }
2304 
2305  if(ftopoConstraint){
2306  V0.SetProductionVertex(Lc);
2307  }
2308 
2309  // After setting the vertex of the V0, getting/filling some info
2310 
2311  //* Get V0 decayLength
2312  Double_t decayLengthV0 = 999999, sigmaDecayLengthV0 = 999999;
2313  Int_t retDLV0 = V0.GetDecayLength( decayLengthV0, sigmaDecayLengthV0 );
2314  if( retDLV0 ) {
2315  if (sigmaDecayLengthV0 > 1e19) {
2316  codeKFV0 = 3; // DecayLength not ok
2317  if (massV0 < 0) {
2318  codeKFV0 = 6; // DecayLength and Mass not ok
2319  if (sigmaMassV0 > 1e19) codeKFV0 = 11; // DecayLength and Mass and SigmaMass not ok
2320  }
2321  else if (sigmaMassV0 > 1e19) codeKFV0 = 8; // DecayLength and SigmaMass not ok
2322  }
2323  }
2324  fHistoDecayLengthKFV0->Fill(decayLengthV0, sigmaDecayLengthV0);
2325 
2326  //* Get V0 life time
2327  Double_t lifeTimeV0 = 999999, sigmaLifeTimeV0 = 999999;
2328  Int_t retTLV0 = V0.GetLifeTime( lifeTimeV0, sigmaLifeTimeV0 );
2329  if( retTLV0 ) {
2330  if (sigmaLifeTimeV0 > 1e19) {
2331  codeKFV0 = 4; // LifeTime not ok
2332  if (sigmaDecayLengthV0 > 1e19) {
2333  codeKFV0 = 9; // LifeTime and DecayLength not ok
2334  if (massV0 < 0) {
2335  codeKFV0 = 14; // LifeTime and DecayLength and Mass not ok
2336  if (sigmaMassV0 > 1e19) codeKFV0 = 15; // LifeTime and DecayLength and Mass and SigmaMass not ok
2337  }
2338  else if (sigmaMassV0 > 1e19) codeKFV0 = 13; // LifeTime and DecayLength and SigmaMass not ok
2339  }
2340  else if (massV0 < 0) { // LifeTime and Mass and SigmaMass not ok
2341  codeKFV0 = 7; // LifeTime and Mass not ok
2342  if (sigmaMassV0 > 1e19) codeKFV0 = 12; // LifeTime and Mass and SigmaMass not ok
2343  }
2344  else if (sigmaMassV0 > 1e19) codeKFV0 = 10; // LifeTime and SigmaMass not ok
2345  }
2346  }
2347  fHistoLifeTimeKFV0->Fill(lifeTimeV0, sigmaLifeTimeV0);
2348 
2349  if (codeKFV0 == -1) codeKFV0 = 0;
2350  fHistoKFV0->Fill(codeKFV0);
2351 
2352  AliDebug(2, Form("V0: mass = %f, decay length = %f, life time = %f", massV0, decayLengthV0, lifeTimeV0 ));
2353 
2354  fHistoMassV0All->Fill(massV0);
2355  fHistoDecayLengthV0All->Fill(decayLengthV0);
2356  fHistoLifeTimeV0All->Fill(lifeTimeV0);
2357 
2358  Double_t qtAlphaV0[2] = {0., 0.};
2359  Double_t vtxV0KF[3] = {V0.GetX(), V0.GetY(), V0.GetZ()};
2360  positiveV0KF.TransportToPoint(vtxV0KF);
2361  negativeV0KF.TransportToPoint(vtxV0KF);
2362  V0.GetArmenterosPodolanski(positiveV0KF, negativeV0KF, qtAlphaV0);
2363  AliDebug(2, Form("Armenteros-Podolanski variables: alpha = %f, qt = %f", qtAlphaV0[1], qtAlphaV0[0]));
2364  fHistoArmenterosPodolanskiV0KF->Fill(qtAlphaV0[1], qtAlphaV0[0]);
2365  fHistoArmenterosPodolanskiV0AOD->Fill(v0part->AlphaV0(), v0part->PtArmV0());
2366  armPolKF[0] = qtAlphaV0[1];
2367  armPolKF[1] = qtAlphaV0[0];
2368 
2369  // Checking MC info for V0
2370 
2371  AliAODMCParticle *motherV0 = 0x0;
2372  AliAODMCParticle *daughv01 = 0x0;
2373  AliAODMCParticle *daughv02 = 0x0;
2374 
2375  if (fUseMCInfo) {
2376  daughv01 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelsv0daugh[0]));
2377  daughv02 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelsv0daugh[1]));
2378  if (!daughv01 && labelsv0daugh[0] > 0){
2379  AliDebug(2, "Could not access MC info for first daughter of V0, continuing");
2380  isMCokV0 = kFALSE;
2381  }
2382  if (!daughv02 && labelsv0daugh[1] > 0){
2383  AliDebug(2, "Could not access MC info for second daughter of V0, continuing");
2384  isMCokV0 = kFALSE;
2385  }
2386  if (isMCokV0){
2387  if( daughv01->GetMother() == daughv02->GetMother() && daughv01->GetMother()>=0 ){
2388  AliDebug(3, Form("The mother has label %d", daughv01->GetMother()));
2389  motherV0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughv01->GetMother()));
2390  if( motherV0 && TMath::Abs(motherV0->GetPdgCode()) != 21 ){ // These are all V0 that are truly V0, not only K0S, but no gluons
2391  if( motherV0->GetNDaughters() == 2 ){
2392  fHistoMassV0True->Fill(massV0);
2393  fHistoDecayLengthV0True->Fill(decayLengthV0);
2394  fHistoLifeTimeV0True->Fill(lifeTimeV0);
2395  fHistoMassV0TrueFromAOD->Fill(v0part->MassK0Short());
2396  if (TMath::Abs(motherV0->GetPdgCode()) == 310){ // These are true V0 that are also K0S
2397  fHistoMassV0TrueK0S->Fill(massV0);
2398  fHistoDecayLengthV0TrueK0S->Fill(decayLengthV0);
2399  fHistoLifeTimeV0TrueK0S->Fill(lifeTimeV0);
2400  fHistoMassV0TrueK0SFromAOD->Fill(v0part->MassK0Short());
2401  }
2402  }
2403  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()));
2404  }
2405  else if (!motherV0){
2406  AliDebug(3, "could not access MC info for mother, continuing");
2407  }
2408  else {
2409  AliDebug(3, "MC mother is a gluon, continuing");
2410  }
2411  }
2412  else {
2413  AliDebug(3, "Background V0!");
2414  isBkgV0 = kTRUE;
2415  }
2416  }
2417  }
2418 
2419  AliDebug(2, Form("isMCokV0 = %d, isBkgV0 = %d", (Int_t)isMCokV0, (Int_t)isBkgV0));
2420 
2421  // Going back to Lc
2422 
2423  //* Get Lc invariant mass
2424  Double_t massLc = 999999, sigmaMassLc= 999999;
2425  Int_t retMLc = Lc.GetMass( massLc, sigmaMassLc );
2426  if( retMLc ) {
2427  AliDebug(3, Form("----> Could not get mass (%e), and sigma(%e) for Lc, continuing", massLc, sigmaMassLc));
2428  if (massLc < 0) {
2429  codeKFLc = 1; // Mass not ok
2430  if (sigmaMassLc > 1e19) codeKFLc = 5; // Mass and SigmaMass not ok
2431  }
2432  else if (sigmaMassLc > 1e19) codeKFLc = 2; // SigmaMass not ok
2433  }
2434  fHistoMassKFLc->Fill(massLc, sigmaMassLc);
2435 
2436  //* Get Lc Decay length
2437  Double_t decayLengthLc = 999999, sigmaDecayLengthLc = 999999;
2438  Int_t retDLLc = Lc.GetDecayLength( decayLengthLc, sigmaDecayLengthLc );
2439  if( retDLLc ) {
2440  AliDebug(3, "----> Lc: Could not get decay length, and sigma");
2441  if (sigmaDecayLengthLc > 1e19) {
2442  codeKFLc = 3; // DecayLength not ok
2443  if (massLc < 0) {
2444  codeKFLc = 6; // DecayLength and Mass not ok
2445  if (sigmaMassLc > 1e19) codeKFLc = 11; // DecayLength and Mass and SigmaMass not ok
2446  }
2447  else if (sigmaMassLc > 1e19) codeKFLc = 8; // DecayLength and SigmaMass not ok
2448  }
2449  }
2450  AliDebug(3, Form("retDLLc = %d, with decayLength = %f and error = %e", retDLLc, decayLengthLc, sigmaDecayLengthLc));
2451 
2452  fHistoDecayLengthKFLc->Fill(decayLengthLc, sigmaDecayLengthLc);
2453 
2454  //* Get Lc life time
2455  Double_t lifeTimeLc = 999999, sigmaLifeTimeLc = 999999;
2456  Int_t retTLLc = Lc.GetLifeTime( lifeTimeLc, sigmaLifeTimeLc );
2457  if( retTLLc ) {
2458  AliDebug(3, "----> Lc: Could not get lifeTime, and sigma");
2459  if (sigmaLifeTimeLc > 1e19) {
2460  codeKFLc = 4; // LifeTime not ok
2461  if (sigmaDecayLengthLc > 1e19) {
2462  codeKFLc = 9; // LifeTime and DecayLength not ok
2463  if (massLc < 0) {
2464  codeKFLc = 14; // LifeTime and DecayLength and Mass not ok
2465  if (sigmaMassLc > 1e19) codeKFLc = 15; // LifeTime and DecayLength and Mass and SigmaMass not ok
2466  }
2467  else if (sigmaMassLc > 1e19) codeKFLc = 13; // LifeTime and DecayLength and SigmaMass not ok
2468  }
2469  else if (massLc < 0) { // LifeTime and Mass and SigmaMass not ok
2470  codeKFLc = 7; // LifeTime and Mass not ok
2471  if (sigmaMassLc > 1e19) codeKFLc = 12; // LifeTime and Mass and SigmaMass not ok
2472  }
2473  else if (sigmaMassLc > 1e19) codeKFLc = 10; // LifeTime and SigmaMass not ok
2474  }
2475  }
2476 
2477  fHistoLifeTimeKFLc->Fill(lifeTimeLc, sigmaLifeTimeLc);
2478 
2479  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));
2480 
2481  if (codeKFLc == -1) codeKFLc = 0;
2482  fHistoKFLc->Fill(codeKFLc);
2483 
2484  fHistoKF->Fill(codeKFV0, codeKFLc);
2485 
2486  // here we fill the histgrams for all the reconstructed KF vertices for the cascade
2487  fHistoMassLcAll->Fill(massLc);
2488  fHistoDecayLengthLcAll->Fill(decayLengthLc);
2489  fHistoLifeTimeLcAll->Fill(lifeTimeLc);
2490 
2491  fHistoMassV0fromLcAll->Fill(massV0);
2492  fHistoDecayLengthV0fromLcAll->Fill(decayLengthV0);
2493  fHistoLifeTimeV0fromLcAll->Fill(lifeTimeV0);
2494 
2495  Double_t xV0 = V0.GetX();
2496  Double_t yV0 = V0.GetY();
2497  Double_t zV0 = V0.GetZ();
2498 
2499  Double_t xLc = Lc.GetX();
2500  Double_t yLc = Lc.GetY();
2501  Double_t zLc = Lc.GetZ();
2502 
2503  Double_t xPrimVtx = primVtxCopy.GetX();
2504  Double_t yPrimVtx = primVtxCopy.GetY();
2505  Double_t zPrimVtx = primVtxCopy.GetZ();
2506 
2507  Double_t distanceLcToPrimVtx = TMath::Sqrt((xPrimVtx - xLc) * (xPrimVtx - xLc) +
2508  (yPrimVtx - yLc) * (yPrimVtx - yLc) +
2509  (zPrimVtx - zLc) * (zPrimVtx - zLc));
2510 
2511  Double_t distanceV0ToPrimVtx = TMath::Sqrt((xPrimVtx - xV0) * (xPrimVtx - xV0) +
2512  (yPrimVtx - yV0) * (yPrimVtx - yV0) +
2513  (zPrimVtx - zV0) * (zPrimVtx - zV0));
2514 
2515  Double_t distanceV0ToLc = TMath::Sqrt((xLc - xV0)*(xLc - xV0) +
2516  (yLc - yV0)*(yLc - yV0) +
2517  (zLc - zV0)*(zLc - zV0));
2518 
2519  //Printf("distanceLcToPrimVtx = %e, distanceV0ToPrimVtx= %f, distanceV0ToLc = %f", distanceLcToPrimVtx, distanceV0ToPrimVtx, distanceV0ToLc);
2520 
2521  fHistoDistanceLcToPrimVtx->Fill(distanceLcToPrimVtx);
2522  fHistoDistanceV0ToPrimVtx->Fill(distanceV0ToPrimVtx);
2523  fHistoDistanceV0ToLc->Fill(distanceV0ToLc);
2524 
2525  distances[0] = distanceLcToPrimVtx;
2526  distances[1] = distanceV0ToPrimVtx;
2527  distances[2] = distanceV0ToLc;
2528 
2529  if (fUseMCInfo) {
2530 
2531  AliAODMCParticle *daughv01Lc = 0x0;
2532  AliAODMCParticle *K0S = 0x0;
2533  AliAODMCParticle *daughv02Lc = 0x0;
2534 
2535  if (labelsLcdaugh[0] >= 0) {
2536  // Printf("Getting Bachelor from label %d", labelsLcdaugh[1]);
2537  daughv01Lc = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelsLcdaugh[0]));
2538  if (!daughv01Lc){
2539  AliDebug(3, "Could not access MC info for first daughter of Lc");
2540  isMCokLc = kFALSE;
2541  }
2542  else {
2543  AliDebug(2, Form("The bachelor has label = %d", daughv01Lc->GetLabel()));
2544  if (TMath::Abs(daughv01Lc->GetPdgCode()) != 2212) isBkgLc = kTRUE;
2545  }
2546  }
2547  else { // The bachelor is a fake
2548  isBkgLc = kTRUE;
2549  }
2550 
2551  if (labelsLcdaugh[1] >= 0) {
2552  //Printf("Getting K0S");
2553  K0S = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelsLcdaugh[1]));
2554  if (!K0S) {
2555  AliDebug(3, "Could not access MC info for second daughter of Lc");
2556  isMCokLc = kFALSE;
2557  }
2558  else{
2559  if (TMath::Abs(K0S->GetPdgCode()) != 310) isBkgLc = kTRUE;
2560  }
2561  }
2562  else{
2563  AliDebug(2, "The K0S is not true --> it does not have a label, continuing...");
2564  isBkgLc = kTRUE;
2565  }
2566 
2567  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
2568  if (isMCokLc) { // We can then access its MC info, and it might then be that also the Lc is a true Lc
2569  Int_t iK0 = K0S->GetMother();
2570  if (iK0 < 0) {
2571  Printf("The K0S has no mother... IMPOSSIBLE"); // the K0S MUST have a mother!
2572  }
2573  else { // The K0S has a mother
2574  daughv02Lc = dynamic_cast<AliAODMCParticle*>(mcArray->At(iK0));
2575  if (!daughv02Lc){
2576  AliDebug(3, "Could not access MC info for second daughter of Lc");
2577  }
2578  else { // we can access safely the K0S mother in the MC
2579  if( daughv01Lc && (daughv01Lc->GetMother() == daughv02Lc->GetMother()) && (daughv01Lc->GetMother()>=0) ){ // This is a true cascade! bachelor and V0 come from the same mother
2580  //Printf("Lc: The mother has label %d", daughv01Lc->GetMother());
2581  AliAODMCParticle *motherLc = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughv01Lc->GetMother()));
2582  Int_t pdgMum = 0, pdgBach = 0, pdgV0 = 0;
2583  if (motherLc) pdgMum = motherLc->GetPdgCode();
2584  if (daughv01Lc) pdgBach = daughv01Lc->GetPdgCode();
2585  if (daughv02Lc) pdgV0 = daughv02Lc->GetPdgCode();
2586  AliDebug(2, Form("pdgMum = %d, pdgBach = %d, pdgV0 = %d", pdgMum, pdgBach, pdgV0));
2587 
2588  if( motherLc && TMath::Abs(motherLc->GetPdgCode()) != 21 ){ // These are true cascades, we don't know yet if they are Lc
2589  fHistoMassLcTrue->Fill(massLc);
2590  fHistoDecayLengthLcTrue->Fill(decayLengthLc);
2591  fHistoLifeTimeLcTrue->Fill(lifeTimeLc);
2592  fHistoMassLcTrueFromAOD->Fill(cascade->InvMassLctoK0sP());
2593 
2594  fHistoMassV0fromLcTrue->Fill(massV0);
2595  fHistoDecayLengthV0fromLcTrue->Fill(decayLengthV0);
2596  fHistoLifeTimeV0fromLcTrue->Fill(lifeTimeV0);
2597 
2598  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...
2599  AliDebug(2, Form("IT IS SIGNAL!!! with label = %d", motherLc->GetLabel()));
2600 
2601  fHistoArmenterosPodolanskiV0KFSgn->Fill(qtAlphaV0[1], qtAlphaV0[0]);
2602  fHistoArmenterosPodolanskiV0AODSgn->Fill(v0part->AlphaV0(), v0part->PtArmV0());
2603 
2604  fHistoDistanceLcToPrimVtxSgn->Fill(distanceLcToPrimVtx);
2605  fHistoDistanceV0ToPrimVtxSgn->Fill(distanceV0ToPrimVtx);
2606  fHistoDistanceV0ToLcSgn->Fill(distanceV0ToLc);
2607 
2608  fHistoMassLcSgn->Fill(massLc);
2609  fHistoMassLcSgnFromAOD->Fill(cascade->InvMassLctoK0sP());
2610 
2611  fHistoDecayLengthLcSgn->Fill(decayLengthLc);
2612  fHistoLifeTimeLcSgn->Fill(lifeTimeLc);
2613 
2614  fHistoMassV0fromLcSgn->Fill(massV0);
2615  fHistoDecayLengthV0fromLcSgn->Fill(decayLengthV0);
2616  fHistoLifeTimeV0fromLcSgn->Fill(lifeTimeV0);
2617  }
2618  else {
2619  isBkgLc = kTRUE; // it is not a Lc, since the pdg != 4122
2620  }
2621 
2622  // if we got here, we can compare with MC information; this part is done also for cases where the cascade is not a Lc!
2623  // First, we compare the vtx of the Lc
2624  Double_t xLcMC = motherLc->Xv();
2625  Double_t yLcMC = motherLc->Yv();
2626  Double_t zLcMC = motherLc->Zv();
2627  //Printf("Got MC vtx for Lc");
2628  Double_t xProtonMC = daughv01Lc->Xv();
2629  Double_t yProtonMC = daughv01Lc->Yv();
2630  Double_t zProtonMC = daughv01Lc->Zv();
2631  //Printf("Got MC vtx for Proton");
2632 
2633  Double_t vtxLcResidualToPrimVtx = TMath::Sqrt((xLcMC - xProtonMC) * (xLcMC - xProtonMC) +
2634  (yLcMC - yProtonMC) * (yLcMC - yProtonMC) +
2635  (zLcMC - zProtonMC) * (zLcMC - zProtonMC)) - distanceLcToPrimVtx;
2636 
2637  // Then, we compare the vtx of the K0s
2638  Double_t xV0MC = motherV0->Xv(); //Production vertex of K0S --> Where Lc decays
2639  Double_t yV0MC = motherV0->Yv();
2640  Double_t zV0MC = motherV0->Zv();
2641 
2642  //Printf("Got MC vtx for V0");
2643  Double_t xPionMC = daughv01->Xv(); //Production vertex of Pion --> Where K0S decays
2644  Double_t yPionMC = daughv01->Yv();
2645  Double_t zPionMC = daughv01->Zv();
2646  //Printf("Got MC vtx for Pion");
2647  Printf("Vertices: MC: x = %f, y = %f, z = %f", xPionMC, yPionMC, zPionMC);
2648 
2649  Double_t vtxV0ResidualToLc = TMath::Sqrt((xV0MC - xPionMC) * (xV0MC - xPionMC) +
2650  (yV0MC - yPionMC) * (yV0MC - yPionMC) +
2651  (zV0MC - zPionMC) * (zV0MC - zPionMC)) - distanceV0ToLc;
2652 
2653  // Then, the K0S vertex wrt primary vertex
2654 
2655  Double_t vtxV0ResidualToPrimVtx = TMath::Sqrt((xPionMC - xLcMC) * (xPionMC - xLcMC) +
2656  (yPionMC - yLcMC) * (yPionMC - yLcMC) +
2657  (zPionMC - zLcMC) * (zPionMC - zLcMC)) - distanceV0ToPrimVtx;
2658 
2659  fHistoVtxLcResidualToPrimVtx->Fill(vtxLcResidualToPrimVtx);
2660  fHistoVtxV0ResidualToLc->Fill(vtxV0ResidualToLc);
2661  fHistoVtxV0ResidualToPrimVtx->Fill(vtxV0ResidualToPrimVtx);
2662 
2663  } //closing: if( motherLc && TMath::Abs(motherLc->GetPdgCode()) != 21 )
2664  else if (!motherLc){
2665  AliDebug(2, "We could not access MC info for Lc mother, so we did nothing");
2666  }
2667  else {
2668  AliDebug(2, "MC Lc mother is a gluon, so we did nothing");
2669  }
2670  } //closing: if( daughv01Lc->GetMother() == daughv02Lc->GetMother() && daughv01Lc->GetMother()>=0 )
2671  else {
2672  isBkgLc = kTRUE; // it cannot be a Lc, since the daughters do not have the same mother
2673  }
2674  } // closing: else { // we can access safely the K0S mother in the MC
2675  } // closing: else { // The K0S has a mother
2676  } // closing isMCLcok
2677  } // closing !isBkgLc
2678  } // closing fUseMCInfo
2679 
2680  //Printf("retMV0 = %d, retMLc = %d", retMV0, retMLc);
2681  if ( retMV0 == 0 && retMLc == 0){
2682  V0KF[0] = massV0;
2683  errV0KF[0] = sigmaMassV0;
2684  V0KF[1] = decayLengthV0;
2685  errV0KF[1] = sigmaDecayLengthV0;
2686  V0KF[2] = lifeTimeV0;
2687  errV0KF[2] = sigmaLifeTimeV0;
2688  LcKF[0] = massLc;
2689  errLcKF[0] = sigmaMassLc;
2690  LcKF[1] = decayLengthLc;
2691  errLcKF[1] = sigmaDecayLengthLc;
2692  LcKF[2] = lifeTimeLc;
2693  errLcKF[2] = sigmaLifeTimeLc;
2694  }
2695 
2696  return 1;
2697 
2698 }
2699 //________________________________________________________________________
2701  AliAODTrack* bachelor,
2702  TClonesArray *mcArray ){
2703 
2704  //Printf("In CheckBachelor");
2705 
2708 
2709  Int_t label = bachelor->GetLabel();
2710  if (label == -1) {
2711  return kBachFake;
2712  }
2713 
2714  AliAODMCParticle *mcpart = dynamic_cast<AliAODMCParticle*>(mcArray->At(TMath::Abs(label)));
2715  if (!mcpart) return kBachInvalid;
2716  Int_t pdg = mcpart->PdgCode();
2717  if (TMath::Abs(pdg) != 2212) {
2718  AliDebug(2, Form("Bachelor is not a p, but a particle with pdg code = %d", pdg));
2719  return kBachNoProton;
2720  }
2721  else { // it is a proton
2722  //Int_t labelLc = part->GetLabel();
2723  Int_t labelLc = FindLcLabel(part, mcArray);
2724  //Printf(">>>>>>>>>>>>> label for Lc = %d", labelLc);
2725  Int_t bachelorMotherLabel = mcpart->GetMother();
2726  //Printf(">>>>>>>>>>>>> label for bachelorMother = %d", bachelorMotherLabel);
2727  if (bachelorMotherLabel == -1) {
2728  return kBachPrimary;
2729  }
2730  AliAODMCParticle *bachelorMother = dynamic_cast<AliAODMCParticle*>(mcArray->At(bachelorMotherLabel));
2731  if (!bachelorMother) return kBachInvalid;
2732  Int_t pdgMother = bachelorMother->PdgCode();
2733  if (TMath::Abs(pdgMother) != 4122) {
2734  AliDebug(2, Form("The proton does not come from a Lc, but from a particle with pdgcode = %d", pdgMother));
2735  return kBachNoLambdaMother;
2736  }
2737  else { // it comes from Lc
2738  if (labelLc != bachelorMotherLabel){
2739  //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));
2740  AliDebug(2, Form("The proton comes from a Lc, but it is not the candidate we are analyzing"));
2742  }
2743  else { // it comes from the correct Lc
2744  AliDebug(2, Form("The proton comes from a Lc, and it is EXACTLY the candidate we are analyzing"));
2745  return kBachCorrectLambdaMother;
2746  }
2747  }
2748  }
2749 
2750  return kBachInvalid;
2751 
2752 }
2753 
2754 //________________________________________________________________________
2756  AliAODv0* v0part,
2757  //AliAODTrack* v0part,
2758  TClonesArray *mcArray ){
2759 
2762 
2763  //Printf(" CheckK0S");
2764 
2765  //Int_t labelMatchToMC = v0part->MatchToMC(310, mcArray);
2766  //Int_t label = v0part->GetLabel();
2767  Int_t labelFind = FindV0Label(v0part, mcArray);
2768  //Printf("\n >>>>>>>>>>>>> label for V0 = %d, from MatchToMC = %d, from Find = %d", label, labelMatchToMC, labelFind);
2769  if (labelFind == -1) {
2770  return kK0SFake;
2771  }
2772 
2773  AliAODMCParticle *mcpart = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelFind));
2774  if (!mcpart) return kK0SInvalid;
2775  Int_t pdg = mcpart->PdgCode();
2776  if (TMath::Abs(pdg) != 310) {
2777  AliDebug(2, Form("K0Spart is not a K0S, but a particle with pdg code = %d", pdg));
2778  //AliInfo(Form("K0Spart is not a K0S, but a particle with pdg code = %d", pdg));
2779  return kK0SNoK0S;
2780  }
2781  else { // it is a K0S
2782  //Int_t labelLc = part->GetLabel();
2783  Int_t labelLc = FindLcLabel(part, mcArray);
2784  Int_t K0SpartMotherLabel = mcpart->GetMother();
2785  if (K0SpartMotherLabel == -1) {
2786  return kK0SWithoutMother;
2787  }
2788  AliAODMCParticle *K0SpartMother = dynamic_cast<AliAODMCParticle*>(mcArray->At(K0SpartMotherLabel)); // should be a K0 in case of signal
2789  if (!K0SpartMother) return kK0SInvalid;
2790  Int_t pdgMotherK0S = K0SpartMother->PdgCode();
2791  if (TMath::Abs(pdgMotherK0S) != 311) {
2792  AliDebug(2, Form("The K0S does not come from a K0, but from a particle with pdgcode = %d", pdgMotherK0S));
2793  //AliInfo(Form("The K0S does not come from a K0, but from a particle with pdgcode = %d", pdgMotherK0S));
2794  return kK0SNotFromK0;
2795  }
2796  else { // the K0S comes from a K0
2797  Int_t K0MotherLabel = K0SpartMother->GetMother(); // mother of K0 --> Lc in case of signal
2798  if (K0MotherLabel == -1) {
2799  return kK0Primary;
2800  }
2801  AliAODMCParticle *K0Mother = dynamic_cast<AliAODMCParticle*>(mcArray->At(K0MotherLabel));
2802  if (!K0Mother) return kK0SInvalid;
2803  Int_t pdgK0Mother = K0Mother->PdgCode();
2804  if (TMath::Abs(pdgK0Mother) != 4122) { // the K0 does not come from a Lc
2805  AliDebug(2, Form("The K0 does not come from a Lc, but from a particle with pdgcode = %d", pdgK0Mother));
2806  //AliInfo(Form("The K0 does not come from a Lc, but from a particle with pdgcode = %d", pdgK0Mother));
2807  return kK0NoLambdaMother;
2808  }
2809  else { // the K0 comes from Lc
2810  if (labelLc != K0MotherLabel){ // The K0 comes from a different Lc
2811  AliDebug(2, Form("The K0S comes from a Lc, but it is not the candidate we are analyzing"));
2812  //AliInfo(Form("The K0S comes from a Lc, but it is not the candidate we are analyzing"));
2813  return kK0DifferentLambdaMother;
2814  }
2815  else { // it comes from the correct Lc
2816  AliDebug(2, Form("The K0S comes from a Lc, and it is EXACTLY the candidate we are analyzing"));
2817  //AliInfo(Form("The K0S comes from a Lc, and it is EXACTLY the candidate we are analyzing"));
2818  return kK0CorrectLambdaMother;
2819  }
2820  }
2821  }
2822  }
2823 
2824  return kK0SInvalid;
2825 
2826 }
2827 
2828 //----------------------------------------------------------------------------
2829 Int_t AliAnalysisTaskSELc2V0bachelorTMVA::FindV0Label(AliAODRecoDecay* v0part, TClonesArray *mcArray) const
2830 {
2831 
2832  //Printf(" FindV0Label");
2833 
2835 
2836  Int_t labMother[2]={-1, -1};
2837  AliAODMCParticle *part=0;
2838  AliAODMCParticle *mother=0;
2839  Int_t dgLabels = -1;
2840 
2841  Int_t ndg = v0part->GetNDaughters();
2842  if (ndg != 2) {
2843  AliFatal(Form("IMPOSSIBLE!! There are %d daughters, but they should be 2!", ndg));
2844  }
2845 
2846  for(Int_t i = 0; i < 2; i++) {
2847  AliAODTrack *trk = (AliAODTrack*)v0part->GetDaughter(i);
2848  dgLabels = trk->GetLabel();
2849  if (dgLabels == -1) {
2850  //printf("daughter with negative label %d\n", dgLabels);
2851  return -1;
2852  }
2853  part = (AliAODMCParticle*)mcArray->At(TMath::Abs(dgLabels));
2854  if (!part) {
2855  //printf("no MC particle\n");
2856  return -1;
2857  }
2858  labMother[i] = part->GetMother();
2859  if (labMother[i] != -1){
2860  mother = (AliAODMCParticle*)mcArray->At(labMother[i]);
2861  if(!mother) {
2862  //printf("no MC mother particle\n");
2863  return -1;
2864  }
2865  }
2866  else {
2867  return -1;
2868  }
2869  }
2870 
2871  if (labMother[0] == labMother[1]) return labMother[0];
2872 
2873  else return -1;
2874 
2875 }
2876 
2877 //----------------------------------------------------------------------------
2879 {
2880 
2882 
2883  //Printf(" FindLcLabel");
2884 
2885  AliAODMCParticle *part=0;
2886  AliAODMCParticle *mother=0;
2887  AliAODMCParticle *grandmother=0;
2888  Int_t dgLabels = -1;
2889 
2890  Int_t ndg = cascade->GetNDaughters();
2891  if (ndg != 2) {
2892  AliFatal(Form("IMPOSSIBLE!! There are %d daughters, but they should be 2!", ndg));
2893  }
2894 
2895  // daughter 0 --> bachelor, daughter 1 --> V0
2896  AliAODTrack *trk = (AliAODTrack*)cascade->GetDaughter(0); // bachelor
2897  dgLabels = trk->GetLabel();
2898  if (dgLabels == -1 ) {
2899  //printf("daughter with negative label %d\n", dgLabels);
2900  return -1;
2901  }
2902  part = (AliAODMCParticle*)mcArray->At(TMath::Abs(dgLabels));
2903  if (!part) {
2904  //printf("no MC particle\n");
2905  return -1;
2906  }
2907  Int_t labMotherBach = part->GetMother();
2908  if (labMotherBach == -1){
2909  return -1;
2910  }
2911  mother = (AliAODMCParticle*)mcArray->At(labMotherBach);
2912  if(!mother) {
2913  //printf("no MC mother particle\n");
2914  return -1;
2915  }
2916 
2917  AliAODv0 *v0 = (AliAODv0*)cascade->GetDaughter(1); // V0
2918  dgLabels = FindV0Label(v0, mcArray);
2919  if (dgLabels == -1 ) {
2920  //printf("daughter with negative label (v0 was a fake) %d\n", dgLabels);
2921  return -1;
2922  }
2923  part = (AliAODMCParticle*)mcArray->At(TMath::Abs(dgLabels));
2924  if (!part) {
2925  //printf("no MC particle\n");
2926  return -1;
2927  }
2928  Int_t labMotherv0 = part->GetMother();
2929  if (labMotherv0 == -1){
2930  return -1;
2931  }
2932  mother = (AliAODMCParticle*)mcArray->At(labMotherv0);
2933  if(!mother) {
2934  //printf("no MC mother particle\n");
2935  return -1;
2936  }
2937  Int_t labGrandMotherv0 = mother->GetMother();
2938  if (labGrandMotherv0 == -1){
2939  return -1;
2940  }
2941  grandmother = (AliAODMCParticle*)mcArray->At(labGrandMotherv0);
2942  if(!grandmother) {
2943  //printf("no MC mother particle\n");
2944  return -1;
2945  }
2946 
2947  //Printf("labMotherBach = %d, labMotherv0 = %d, labGrandMotherv0 = %d", labMotherBach, labMotherv0, labGrandMotherv0);
2948  if (labGrandMotherv0 == labMotherBach) return labMotherBach;
2949 
2950  else return -1;
2951 
2952 }
2953 
2954 
2957 
2958  Int_t runNo = event->GetRunNumber();
2959  Int_t period = -1;
2960  switch (fAnalysisType) { // flag to set which system and year is being used
2961  case kpPb2013: //0 = LHC13b, 1 = LHC13c
2962  if (runNo > 195343 && runNo < 195484) period = 0;
2963  if (runNo > 195528 && runNo < 195678) period = 1;
2964  if (period < 0 || period > 1) { AliInfo(Form("Run number %d not found for LHC13!",runNo)); return 0;}
2965  break;
2966  case kpPb2016: //0 = LHC16q, 265499 -- 265525 || 265309 -- 265387, 1 = LHC16q, 265435, 2 = LHC16q, 265388 -- 265427, 3 = LHC16t, 267163 -- 267166
2967  if ((runNo >=265499 && runNo <=265525) || (runNo >= 265309 && runNo <= 265387)) period = 0;
2968  else if (runNo == 265435) period = 1;
2969  else if (runNo >= 265388 && runNo <= 265427) period = 2;
2970  else if (runNo >=267163 && runNo <=276166) period = 3;
2971  if (period < 0 || period > 3) { AliInfo(Form("Run number %d not found for LHC16 pPb!",runNo)); return 0;}
2972  break;
2973  case kpp2016: //0 = LHC16j, 1 = LHC16k, 2 = LHC16l
2974  if (runNo >= 256219 && runNo <= 256418) period = 0;
2975  else if (runNo >= 256504 && runNo <= 258537) period = 1;
2976  else if (runNo >= 258883 && runNo <= 260187) period = 2;
2977  if (period < 0 || period > 2) {AliInfo(Form("Run number %d not found for LHC16 pp!",runNo)); return 0;}
2978  break;
2979  case kpp2010: // 0 = LHC10b, 1 = LHC10c, 2 = LHC10d, 3 = LHC10e
2980  if (runNo > 114930 && runNo < 117223) period = 0;
2981  if (runNo > 119158 && runNo < 120830) period = 1;
2982  if (runNo > 122373 && runNo < 126438) period = 2;
2983  if (runNo > 127711 && runNo < 130851) period = 3;
2984  if (period < 0 || period > 3) {AliInfo(Form("Run number %d not found for LHC10 pp!",runNo)); return 0;}
2985  break;
2986  default: //no valid switch
2987  return 0;
2988  break;
2989  }
2990 
2991  return fMultEstimatorAvg[period];
2992 }
2993 
2994 
2995 
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:258
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:210
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
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:271
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:198
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:68
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.