AliPhysics  9fe175b (9fe175b)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliCFTaskVertexingHFCutVarFDSub.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  * appear 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 //-----------------------------------------------------------------------
17 // Class for HF corrections as a function of many variables
18 // 6 Steps introduced: MC, MC Acc, Reco, Reco Acc, Reco Acc + ITS Cl,
19 // Reco Acc + ITS Cl + PPR cuts
20 // 12 variables used: pt, y, cosThetaStar, ptPi, ptK, ct,
21 // dca, d0Pi, d0K, d0Pixd0K, cosPointingAngle, phi
22 //
23 //-----------------------------------------------------------------------
24 // Author : C. Zampolli, CERN
25 // D. Caffarri, Univ & INFN Padova caffarri@pd.infn.it
26 //-----------------------------------------------------------------------
27 //-----------------------------------------------------------------------
28 // Base class for HF Unfolding (pt and eta)
29 // correlation matrix filled at Acceptance and PPR level
30 // Author: A.Grelli , Utrecht - agrelli@uu.nl
31 //-----------------------------------------------------------------------
32 //-----------------------------------------------------------------------
33 // Modified by for analysis of D0 candidates using the cut-variation
34 // feed-down subtraction method
35 // by
36 // Andrea Rossi <Andrea.Rossi@cern.ch>
37 // Felix Reidt <Felix.Reidt@cern.ch>
38 //-----------------------------------------------------------------------
40 
41 #include "TCanvas.h"
42 #include "TProfile.h"
43 #include "TH1I.h"
44 #include "TH3F.h"
45 #include "TH1F.h"
46 #include "TStyle.h"
47 #include "TFile.h"
48 #include "TF1.h"
49 
50 #include "AliCFManager.h"
51 #include "AliCFContainer.h"
52 #include "AliLog.h"
53 #include "AliInputEventHandler.h"
54 #include "AliAnalysisManager.h"
55 #include "AliAODHandler.h"
56 #include "AliAODEvent.h"
57 #include "AliAODRecoDecayHF.h"
61 #include "AliAODRecoCascadeHF.h"
62 #include "AliAODMCParticle.h"
63 #include "AliAODMCHeader.h"
64 #include "AliESDtrack.h"
65 
66 #include "THnSparse.h"
67 #include "AliESDtrackCuts.h"
68 #include "AliRDHFCuts.h"
69 #include "AliRDHFCutsD0toKpi.h"
72 #include "AliRDHFCutsDstoKKpi.h"
73 #include "AliRDHFCutsLctopKpi.h"
74 #include "AliRDHFCutsD0toKpipipi.h"
75 #include "AliRDHFCutsLctoV0.h"
76 #include "AliCFVertexingHF2Prong.h"
77 #include "AliCFVertexingHF3Prong.h"
80 #include "AliCFVertexingHF.h"
81 #include "AliVertexingHFUtils.h"
82 #include "AliAnalysisDataSlot.h"
83 #include "AliAnalysisDataContainer.h"
84 #include "AliAnalysisTaskSE.h"
85 #include "AliPIDResponse.h"
86 #include "AliHFsubtractBFDcuts.h"
87 
91 
92 //__________________________________________________________________________
94  AliAnalysisTaskSE(),
95  fCFManager(0x0),
96  fHistEventsProcessed(0x0),
97  fCorrelation(0x0),
98  fListProfiles(0),
99  fCountMC(0),
100  fCountAcc(0),
101  fCountVertex(0),
102  fCountRefit(0),
103  fCountReco(0),
104  fCountRecoAcc(0),
105  fCountRecoITSClusters(0),
106  fCountRecoPPR(0),
107  fCountRecoPID(0),
108  fEvents(0),
109  fDecayChannel(0),
110  fFillFromGenerated(kFALSE),
111  fOriginDselection(0),
112  fAcceptanceUnf(kTRUE),
113  fCuts(0),
114  fUseWeight(kFALSE),
115  fWeight(1.),
116  fUseFlatPtWeight(kFALSE),
117  fUseZWeight(kFALSE),
118  fUseNchWeight(kFALSE),
119  fUseTrackletsWeight(kFALSE),
120  fUseMultRatioAsWeight(kFALSE),
121  fNvar(0),
122  fPartName(""),
123  fDauNames(""),
124  fSign(2),
125  fCentralitySelection(kTRUE),
126  fFakeSelection(0),
127  fRejectIfNoQuark(kTRUE),
128  fUseMCVertex(kFALSE),
129  fDsOption(1),
130  fGenDsOption(3),
131  fConfiguration(kCheetah), // by default, setting the fast configuration
132  fFuncWeight(0x0),
133  fHistoPtWeight(0x0),
134  fHistoMotherPtWeight(0x0),
135  fHistoMeasNch(0x0),
136  fHistoMCNch(0x0),
137  fResonantDecay(0),
138  fLctoV0bachelorOption(1),
139  fGenLctoV0bachelorOption(0),
140  fUseSelectionBit(kTRUE),
141  fPDGcode(0),
142  fMultiplicityEstimator(kNtrk10),
143  fRefMult(9.26),
144  fZvtxCorrectedNtrkEstimator(kFALSE),
145  fIsPPData(kFALSE),
146  fIsPPbData(kFALSE),
147  fUseAdditionalCuts(kFALSE),
148  fUseCutsForTMVA(kFALSE),
149  fUseCascadeTaskForLctoV0bachelor(kFALSE),
150  fCutOnMomConservation(0.00001),
151  fObjSpr(0x0),
152  fTHnAnalysis(0x0),
153  fTHnGenerator(0x0),
154  fhBptCutVar(0x0),
155  fhBquarkPt(0x0),
156  fhCquarkPt(0x0),
157  fListBdecays(0x0),
158  fQAHists(0x0)
159 {
160  //
161  //Default ctor
162  //
163  for(Int_t i=0; i<4; i++) fMultEstimatorAvg[i]=0;
164 }
165 //___________________________________________________________________________
167  AliAnalysisTaskSE(name),
168  fCFManager(0x0),
169  fHistEventsProcessed(0x0),
170  fCorrelation(0x0),
171  fListProfiles(0),
172  fCountMC(0),
173  fCountAcc(0),
174  fCountVertex(0),
175  fCountRefit(0),
176  fCountReco(0),
177  fCountRecoAcc(0),
178  fCountRecoITSClusters(0),
179  fCountRecoPPR(0),
180  fCountRecoPID(0),
181  fEvents(0),
182  fDecayChannel(0),
183  fFillFromGenerated(kFALSE),
184  fOriginDselection(0),
185  fAcceptanceUnf(kTRUE),
186  fCuts(cuts),
187  fUseWeight(kFALSE),
188  fWeight(1.),
189  fUseFlatPtWeight(kFALSE),
190  fUseZWeight(kFALSE),
191  fUseNchWeight(kFALSE),
192  fUseTrackletsWeight(kFALSE),
193  fUseMultRatioAsWeight(kFALSE),
194  fNvar(0),
195  fPartName(""),
196  fDauNames(""),
197  fSign(2),
198  fCentralitySelection(kTRUE),
199  fFakeSelection(0),
200  fRejectIfNoQuark(kTRUE),
201  fUseMCVertex(kFALSE),
202  fDsOption(1),
203  fGenDsOption(3),
204  fConfiguration(kCheetah), // by default, setting the fast configuration
205  fFuncWeight(func),
206  fHistoPtWeight(0x0),
207  fHistoMotherPtWeight(0x0),
208  fHistoMeasNch(0x0),
209  fHistoMCNch(0x0),
210  fResonantDecay(0),
211  fLctoV0bachelorOption(1),
212  fGenLctoV0bachelorOption(0),
213  fUseSelectionBit(kTRUE),
214  fPDGcode(0),
215  fMultiplicityEstimator(kNtrk10),
216  fRefMult(9.26),
217  fZvtxCorrectedNtrkEstimator(kFALSE),
218  fIsPPData(kFALSE),
219  fIsPPbData(kFALSE),
220  fUseAdditionalCuts(kFALSE),
221  fUseCutsForTMVA(kFALSE),
222  fUseCascadeTaskForLctoV0bachelor(kFALSE),
223  fCutOnMomConservation(0.00001),
224  fObjSpr(0x0),
225  fTHnAnalysis(0x0),
226  fTHnGenerator(0x0),
227  fhBptCutVar(0x0),
228  fhBquarkPt(0x0),
229  fhCquarkPt(0x0),
230  fListBdecays(0x0),
231  fQAHists(0x0)
232 {
233  //
234  // Constructor. Initialization of Inputs and Outputs
235  //
236  /*
237  DefineInput(0) and DefineOutput(0)
238  are taken care of by AliAnalysisTaskSE constructor
239  */
240  DefineOutput(1,TH1I::Class());
241  DefineOutput(2,AliCFContainer::Class());
242  DefineOutput(3,THnSparseD::Class());
243  DefineOutput(4,AliRDHFCuts::Class());
244  for(Int_t i=0; i<4; i++) fMultEstimatorAvg[i]=0;
245  DefineOutput(5,TList::Class()); // slot #5 keeps the zvtx Ntrakclets correction profiles
246  DefineOutput(6,THnSparseF::Class());
247  DefineOutput(7,TH3F::Class());
248  DefineOutput(8,TH1F::Class());
249  DefineOutput(9,TH1F::Class());
250  DefineOutput(10,TH1F::Class());
251  DefineOutput(11,TList::Class());
252  DefineOutput(12,TList::Class());
253 
254  fCuts->PrintAll();
255 }
256 
257 //___________________________________________________________________________
259 {
260  //
261  // Assignment operator
262  //
263  if (this!=&c) {
264  AliAnalysisTaskSE::operator=(c) ;
267  fCuts = c.fCuts;
273  for(Int_t i=0; i<4; i++) fMultEstimatorAvg[i]=c.fMultEstimatorAvg[i];
274  fObjSpr=c.fObjSpr;
280  fListBdecays=c.fListBdecays; // FIXME: TList copy contructor not implemented
281  fQAHists=c.fQAHists; // FIXME: TList copy contructor not implemented
282  }
283  return *this;
284 }
285 
286 //___________________________________________________________________________
288  AliAnalysisTaskSE(c),
289  fCFManager(c.fCFManager),
290  fHistEventsProcessed(c.fHistEventsProcessed),
291  fCorrelation(c.fCorrelation),
292  fListProfiles(c.fListProfiles),
293  fCountMC(c.fCountMC),
294  fCountAcc(c.fCountAcc),
295  fCountVertex(c.fCountVertex),
296  fCountRefit(c.fCountRefit),
297  fCountReco(c.fCountReco),
298  fCountRecoAcc(c.fCountRecoAcc),
299  fCountRecoITSClusters(c.fCountRecoITSClusters),
300  fCountRecoPPR(c.fCountRecoPPR),
301  fCountRecoPID(c.fCountRecoPID),
302  fEvents(c.fEvents),
303  fDecayChannel(c.fDecayChannel),
304  fFillFromGenerated(c.fFillFromGenerated),
305  fOriginDselection(c.fOriginDselection),
306  fAcceptanceUnf(c.fAcceptanceUnf),
307  fCuts(c.fCuts),
308  fUseWeight(c.fUseWeight),
309  fWeight(c.fWeight),
310  fUseFlatPtWeight(c.fUseFlatPtWeight),
311  fUseZWeight(c.fUseZWeight),
312  fUseNchWeight(c.fUseNchWeight),
313  fUseTrackletsWeight(c.fUseTrackletsWeight),
314  fUseMultRatioAsWeight(c.fUseMultRatioAsWeight),
315  fNvar(c.fNvar),
316  fPartName(c.fPartName),
317  fDauNames(c.fDauNames),
318  fSign(c.fSign),
319  fCentralitySelection(c.fCentralitySelection),
320  fFakeSelection(c.fFakeSelection),
321  fRejectIfNoQuark(c.fRejectIfNoQuark),
322  fUseMCVertex(c.fUseMCVertex),
323  fDsOption(c.fDsOption),
324  fGenDsOption(c.fGenDsOption),
325  fConfiguration(c.fConfiguration),
326  fFuncWeight(c.fFuncWeight),
327  fHistoPtWeight(c.fHistoPtWeight),
328  fHistoMotherPtWeight(c.fHistoMotherPtWeight),
329  fHistoMeasNch(c.fHistoMeasNch),
330  fHistoMCNch(c.fHistoMCNch),
331  fResonantDecay(c.fResonantDecay),
332  fLctoV0bachelorOption(c.fLctoV0bachelorOption),
333  fGenLctoV0bachelorOption(c.fGenLctoV0bachelorOption),
334  fUseSelectionBit(c.fUseSelectionBit),
335  fPDGcode(c.fPDGcode),
336  fMultiplicityEstimator(c.fMultiplicityEstimator),
337  fRefMult(c.fRefMult),
338  fZvtxCorrectedNtrkEstimator(c.fZvtxCorrectedNtrkEstimator),
339  fIsPPData(c.fIsPPData),
340  fIsPPbData(c.fIsPPbData),
341  fUseAdditionalCuts(c.fUseAdditionalCuts),
342  fUseCutsForTMVA(c.fUseCutsForTMVA),
343  fUseCascadeTaskForLctoV0bachelor(c.fUseCascadeTaskForLctoV0bachelor),
344  fCutOnMomConservation(c.fCutOnMomConservation),
345  fObjSpr(c.fObjSpr),
346  fTHnAnalysis(c.fTHnAnalysis),
347  fTHnGenerator(c.fTHnGenerator),
348  fhBptCutVar(c.fhBptCutVar),
349  fhBquarkPt(c.fhBquarkPt),
350  fhCquarkPt(c.fhCquarkPt),
351  fListBdecays(c.fListBdecays), // FIXME: TList copy contructor not implemented
352  fQAHists(c.fQAHists) // FIXME: TList copy contructor not implemented
353 {
354  //
355  // Copy Constructor
356  //
357  for(Int_t i=0; i<4; i++) fMultEstimatorAvg[i]=c.fMultEstimatorAvg[i];
358 }
359 
360 //___________________________________________________________________________
362 {
363  //
364  //destructor
365  //
366  if (fCFManager) { delete fCFManager; fCFManager=0x0; }
368  if (fCorrelation) { delete fCorrelation; fCorrelation=0x0; }
369  if (fListProfiles) { delete fListProfiles; fListProfiles=0x0; }
370  if (fCuts) { delete fCuts; fCuts=0x0; }
371  if (fFuncWeight) { delete fFuncWeight; fFuncWeight=0x0; }
372  if (fHistoPtWeight) { delete fHistoPtWeight; fHistoPtWeight=0x0; }
374  if (fHistoMeasNch) { delete fHistoMeasNch; fHistoMeasNch=0x0; }
375  if (fHistoMCNch) { delete fHistoMCNch; fHistoMCNch=0x0; }
376  if (fObjSpr) { delete fObjSpr; fObjSpr=0x0; }
377  if (fTHnAnalysis) { delete fTHnAnalysis; fTHnAnalysis=0x0; }
378  if (fTHnGenerator) { delete fTHnGenerator; fTHnGenerator=0x0; }
379  if (fhBptCutVar) { delete fhBptCutVar; fhBptCutVar=0x0; }
380  if (fhBquarkPt) { delete fhBquarkPt; fhBquarkPt=0x0; }
381  if (fhCquarkPt) { delete fhCquarkPt; fhCquarkPt=0x0; }
382  if (fListBdecays) { delete fListBdecays; fListBdecays=0x0; }
383  if (fQAHists) { delete fQAHists; fQAHists=0x0; }
384  for(Int_t i=0; i<4; i++) {
385  if(fMultEstimatorAvg[i]) { delete fMultEstimatorAvg[i]; fMultEstimatorAvg[i]=0x0; }
386  }
387 }
388 
389 //__________________________________________________________________________
391 {
392  //
393  // Initialization
394  //
395 
396  if (fDebug>1) printf("AliCFTaskVertexingHFCutVarFDSub::Init()");
397  if(fUseWeight && fUseZWeight) { AliFatal("Can not use at the same time pt and z-vtx weights, please choose"); return; }
398  if(fUseWeight && fUseNchWeight) { AliInfo("Beware, using at the same time pt and Nch weights, please check"); }
399  if(fUseNchWeight && !fHistoMCNch) { AliFatal("Need to pass the MC Nch distribution to use Nch weights"); return; }
401 
402  AliRDHFCuts *copyfCuts = 0x0;
403  if (!fCuts){
404  AliFatal("No cuts defined - Exiting...");
405  return;
406  }
407 
408  switch (fDecayChannel){
409  case 2:{
410  fPDGcode = 421;
411  copyfCuts = new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fCuts)));
412  switch (fConfiguration) {
413  case kSnail: // slow configuration: all variables in
414  fNvar = 16;
415  break;
416  case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled
417  fNvar = 8;
418  break;
419  }
420  fPartName="D0";
421  fDauNames="K+pi";
422  break;
423  }
424  case 21:{
425  fPDGcode = 413;
426  copyfCuts = new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fCuts)));
427  switch (fConfiguration) {
428  case kSnail: // slow configuration: all variables in
429  fNvar = 16;
430  break;
431  case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled
432  fNvar = 8;
433  break;
434  }
435  fPartName="Dstar";
436  fDauNames="K+pi+pi";
437  break;
438  }
439  case 22:{
440  fPDGcode = 4122;
441  copyfCuts = new AliRDHFCutsLctoV0(*(static_cast<AliRDHFCutsLctoV0*>(fCuts)));
442  switch (fConfiguration) {
443  case kSnail: // slow configuration: all variables in
444  fNvar = 16;
445  break;
446  case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled
447  fNvar = 8;
448  break;
449  }
450  fPartName="Lambdac";
451  fDauNames="V0+bachelor";
452  break;
453  }
454  case 31:{
455  fPDGcode = 411;
456  copyfCuts = new AliRDHFCutsDplustoKpipi(*(static_cast<AliRDHFCutsDplustoKpipi*>(fCuts)));
457  switch (fConfiguration) {
458  case kSnail: // slow configuration: all variables in
459  fNvar = 14;
460  break;
461  case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled
462  fNvar = 8;
463  break;
464  }
465  fPartName="Dplus";
466  fDauNames="K+pi+pi";
467  break;
468  }
469  case 32:{
470  fPDGcode = 4122;
471  copyfCuts = new AliRDHFCutsLctopKpi(*(static_cast<AliRDHFCutsLctopKpi*>(fCuts)));
472  switch (fConfiguration) {
473  case kSnail: // slow configuration: all variables in
474  fNvar = 18;
475  break;
476  case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled
477  fNvar = 8;
478  break;
479  }
480  fPartName="Lambdac";
481  fDauNames="p+K+pi";
482  break;
483  }
484  case 33:{
485  fPDGcode = 431;
486  copyfCuts = new AliRDHFCutsDstoKKpi(*(static_cast<AliRDHFCutsDstoKKpi*>(fCuts)));
487  switch (fConfiguration) {
488  case kSnail: // slow configuration: all variables in
489  fNvar = 14;
490  break;
491  case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled
492  fNvar = 8;
493  break;
494  }
495  fPartName="Ds";
496  fDauNames="K+K+pi";
497  break;
498  }
499  case 4:{
500  fPDGcode = 421;
501  copyfCuts = new AliRDHFCutsD0toKpipipi(*(static_cast<AliRDHFCutsD0toKpipipi*>(fCuts)));
502  switch (fConfiguration) {
503  case kSnail: // slow configuration: all variables in
504  fNvar = 16;
505  break;
506  case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled
507  fNvar = 8;
508  break;
509  }
510  fPartName="D0";
511  fDauNames="K+pi+pi+pi";
512  break;
513  }
514  default:
515  AliFatal("The decay channel MUST be defined according to AliCFVertexing::DecayChannel - Exiting...");
516  break;
517  }
518 
519  const char* nameoutput=GetOutputSlot(4)->GetContainer()->GetName();
520  if (copyfCuts){
521  copyfCuts->SetName(nameoutput);
522 
523  //Post the data
524  PostData(4, copyfCuts);
525  }
526  else{
527  AliFatal("Failing initializing AliRDHFCuts object - Exiting...");
528  }
529 
530  fListProfiles = new TList();
531  fListProfiles->SetOwner();
532  TString period[4];
533  Int_t nProfiles=4;
534 
535  if (fIsPPbData) { //if pPb, use only two estimator histos
536  period[0] = "LHC13b"; period[1] = "LHC13c";
537  nProfiles = 2;
538  } else { // else assume pp (four histos for LHC10)
539  period[0] = "LHC10b"; period[1] = "LHC10c"; period[2] = "LHC10d"; period[3] = "LHC10e";
540  nProfiles = 4;
541  }
542 
543  for(Int_t i=0; i<nProfiles; i++){
544  if(fMultEstimatorAvg[i]){
545  TProfile* hprof=new TProfile(*fMultEstimatorAvg[i]);
546  hprof->SetName(Form("ProfileTrkVsZvtx%s\n",period[i].Data()));
547  fListProfiles->Add(hprof);
548  }
549  }
550 
557 
558  PostData(5,fListProfiles);
559 
560  return;
561 }
562 
563 //_________________________________________________
565 {
566  //
567  // Main loop function
568  //
569 
570  PostData(1,fHistEventsProcessed) ;
571  PostData(2,fCFManager->GetParticleContainer()) ;
572  PostData(3,fCorrelation);
573 
574  AliDebug(3,Form("*** Processing event %d\n", fEvents));
575 
576  if (fFillFromGenerated){
577  AliWarning("Flag to fill container with generated value ON ---> dca, d0pi, d0K, d0xd0, cosPointingAngle will be set as dummy!");
578  }
579 
580  if (!fInputEvent) {
581  Error("UserExec","NO EVENT FOUND!");
582  return;
583  }
584 
585  AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
586 
587  TClonesArray *arrayBranch=0;
588 
589  if(!aodEvent && AODEvent() && IsStandardAOD()) {
590  // In case there is an AOD handler writing a standard AOD, use the AOD
591  // event in memory rather than the input (ESD) event.
592  aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
593  // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
594  // have to taken from the AOD event hold by the AliAODExtension
595  AliAODHandler* aodHandler = (AliAODHandler*)
596  ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
597  if(aodHandler->GetExtensions()) {
598  AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
599  AliAODEvent *aodFromExt = ext->GetAOD();
600 
601  switch (fDecayChannel){
602  case 2:{
603  arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
604  break;
605  }
606  case 21:{
607  arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
608  break;
609  }
610  case 22:{
611  arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("CascadesHF");
612  break;
613  }
614  case 31:
615  case 32:
616  case 33:{
617  arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
618  break;
619  }
620  case 4:{
621  arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm4Prong");
622  break;
623  }
624  default:
625  break;
626  }
627  }
628  }
629  else {
630  switch (fDecayChannel){
631  case 2:{
632  arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("D0toKpi");
633  break;
634  }
635  case 21:{
636  arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("Dstar");
637  break;
638  }
639  case 22:{
640  arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("CascadesHF");
641  break;
642  }
643  case 31:
644  case 32:
645  case 33:{
646  arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("Charm3Prong");
647  break;
648  }
649  case 4:{
650  arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("Charm4Prong");
651  break;
652  }
653  default:
654  break;
655  }
656  }
657 
658  AliAODVertex *aodVtx = (AliAODVertex*)aodEvent->GetPrimaryVertex();
659  if (!aodVtx) {
660  AliDebug(3, "The event was skipped due to missing vertex");
661  return;
662  }
663 
664  if (!arrayBranch) {
665  AliError("Could not find array of HF vertices");
666  return;
667  }
668 
669  fEvents++;
670 
671  fCFManager->SetRecEventInfo(aodEvent);
672  fCFManager->SetMCEventInfo(aodEvent);
673 
674  //******** DEFINE number of variables of the container***** for now set at 13, in the future in the config macro.
675 
676  //loop on the MC event
677 
678  TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
679  if (!mcArray) {
680  AliError("Could not find Monte-Carlo in AOD");
681  return;
682  }
683  Int_t icountMC = 0;
684  Int_t icountAcc = 0;
685  Int_t icountReco = 0;
686  Int_t icountVertex = 0;
687  Int_t icountRefit = 0;
688  Int_t icountRecoAcc = 0;
689  Int_t icountRecoITSClusters = 0;
690  Int_t icountRecoPPR = 0;
691  Int_t icountRecoPID = 0;
692  Int_t cquarks = 0;
693 
694  AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
695  if (!mcHeader) {
696  AliError("Could not find MC Header in AOD");
697  return;
698  }
699 
700  fHistEventsProcessed->Fill(0.5);
701 
702  Double_t* containerInput = new Double_t[fNvar];
703  Double_t* containerInputMC = new Double_t[fNvar];
704 
705 
706  AliCFVertexingHF* cfVtxHF=0x0;
707  switch (fDecayChannel){
708  case 2:{
709  cfVtxHF = new AliCFVertexingHF2Prong(mcArray, fOriginDselection);
710  break;
711  }
712  case 21:{
713  cfVtxHF = new AliCFVertexingHFCascade(mcArray, fOriginDselection);
714  ((AliCFVertexingHFCascade*)cfVtxHF)->SetPDGcascade(413);
715  ((AliCFVertexingHFCascade*)cfVtxHF)->SetPDGbachelor(211);
716  ((AliCFVertexingHFCascade*)cfVtxHF)->SetPDGneutrDaugh(421);
717  ((AliCFVertexingHFCascade*)cfVtxHF)->SetPDGneutrDaughForMC(421);
718  ((AliCFVertexingHFCascade*)cfVtxHF)->SetPDGneutrDaughPositive(211);
719  ((AliCFVertexingHFCascade*)cfVtxHF)->SetPDGneutrDaughNegative(321);
720  ((AliCFVertexingHFCascade*)cfVtxHF)->SetPrimaryVertex(aodVtx);
721  break;
722  }
723  case 22:{
724  // Lc -> K0S+proton
726  cfVtxHF = new AliCFVertexingHFCascade(mcArray, fOriginDselection);
727  ((AliCFVertexingHFCascade*)cfVtxHF)->SetPDGcascade(4122);
728  ((AliCFVertexingHFCascade*)cfVtxHF)->SetPDGbachelor(2212);
729  ((AliCFVertexingHFCascade*)cfVtxHF)->SetPDGneutrDaugh(310);
730  ((AliCFVertexingHFCascade*)cfVtxHF)->SetPDGneutrDaughForMC(311);
731  ((AliCFVertexingHFCascade*)cfVtxHF)->SetPDGneutrDaughPositive(211);
732  ((AliCFVertexingHFCascade*)cfVtxHF)->SetPDGneutrDaughNegative(211);
733  ((AliCFVertexingHFCascade*)cfVtxHF)->SetPrimaryVertex(aodVtx);
734  ((AliCFVertexingHFCascade*)cfVtxHF)->SetCutOnMomConservation(fCutOnMomConservation);
735  if (fUseAdditionalCuts) ((AliCFVertexingHFCascade*)cfVtxHF)->SetUseCutsForTMVA(fUseCutsForTMVA);
736  }
737  else {
739  }
740  break;
741  }
742  case 31:
743  // case 32:
744  case 33:{
745  cfVtxHF = new AliCFVertexingHF3Prong(mcArray, fOriginDselection, fDecayChannel);
746  if(fDecayChannel==33){
747  ((AliCFVertexingHF3Prong*)cfVtxHF)->SetGeneratedDsOption(fGenDsOption);
748  }
749  break;
750  }
751  case 32:{
753  }
754  case 4:{
755  //cfVtxHF = new AliCFVertexingHF4Prong(mcArray, originDselection); // not there yet
756  break;
757  }
758  default:
759  break;
760  }
761  if (!cfVtxHF){
762  AliError("No AliCFVertexingHF initialized");
763  delete[] containerInput;
764  delete[] containerInputMC;
765  return;
766  }
767 
768  Double_t zPrimVertex = aodVtx ->GetZ();
769  Double_t zMCVertex = mcHeader->GetVtxZ();
770  Int_t runnumber = aodEvent->GetRunNumber();
771 
772  // Multiplicity definition with tracklets
773  Double_t nTracklets = 0;
774  Int_t nTrackletsEta10 = static_cast<Int_t>(AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aodEvent,-1.,1.));
775  Int_t nTrackletsEta16 = static_cast<Int_t>(AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aodEvent,-1.6,1.6));
776  nTracklets = (Double_t)nTrackletsEta10;
777  if(fMultiplicityEstimator==kNtrk10to16) { nTracklets = (Double_t)(nTrackletsEta16 - nTrackletsEta10); }
778 
779  // Apply the Ntracklets z-vtx data driven correction
781  TProfile* estimatorAvg = GetEstimatorHistogram(aodEvent);
782  if(estimatorAvg) {
783  Int_t nTrackletsEta10Corr = static_cast<Int_t>(AliVertexingHFUtils::GetCorrectedNtracklets(estimatorAvg,nTrackletsEta10,zPrimVertex,fRefMult));
784  Int_t nTrackletsEta16Corr = static_cast<Int_t>(AliVertexingHFUtils::GetCorrectedNtracklets(estimatorAvg,nTrackletsEta16,zPrimVertex,fRefMult));
785  nTracklets = (Double_t)nTrackletsEta10Corr;
786  if(fMultiplicityEstimator==kNtrk10to16) { nTracklets = (Double_t)(nTrackletsEta16Corr - nTrackletsEta10Corr); }
787  }
788  }
789 
790 
791  fWeight=1.;
792  if(fUseZWeight) fWeight *= GetZWeight(zMCVertex,runnumber);
793  if(fUseNchWeight){
794  Int_t nChargedMCPhysicalPrimary=AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(mcArray,-1.0,1.0);
795  if(!fUseTrackletsWeight) fWeight *= GetNchWeight(nChargedMCPhysicalPrimary);
796  else fWeight *= GetNchWeight(static_cast<Int_t>(nTracklets));
797  AliDebug(2,Form("Using Nch weights, Mult=%d Weight=%f\n",nChargedMCPhysicalPrimary,fWeight));
798  }
799  Double_t eventWeight=fWeight;
800 
801  if (TMath::Abs(zMCVertex) > fCuts->GetMaxVtxZ()){
802  AliDebug(3,Form("z coordinate of MC vertex = %f, it was required to be within [-%f, +%f], skipping event", zMCVertex, fCuts->GetMaxVtxZ(), fCuts->GetMaxVtxZ()));
803  delete[] containerInput;
804  delete[] containerInputMC;
805  delete cfVtxHF;
806  return;
807  }
808 
809  if(aodEvent->GetTriggerMask()==0 &&
810  (runnumber>=195344 && runnumber<=195677)){
811  AliDebug(3,"Event rejected because of null trigger mask");
812  delete[] containerInput;
813  delete[] containerInputMC;
814  delete cfVtxHF;
815  return;
816  }
817 
818  AliESDtrackCuts** trackCuts = new AliESDtrackCuts*[cfVtxHF->GetNProngs()];
819  if (fDecayChannel == 21){
820  // for the D*, setting the third element of the array of the track cuts to those for the soft pion
821  for (Int_t iProng = 0; iProng<cfVtxHF->GetNProngs()-1; iProng++){
822  trackCuts[iProng]=fCuts->GetTrackCuts();
823  }
824  trackCuts[2] = fCuts->GetTrackCutsSoftPi();
825  }
826  else if (fDecayChannel == 22) {
827  // for the Lc->V0+bachelor, setting the second and third elements of the array of the track cuts to those for the V0 daughters
828  trackCuts[0]=fCuts->GetTrackCuts();
829  trackCuts[1]=fCuts->GetTrackCutsV0daughters();
830  trackCuts[2]=fCuts->GetTrackCutsV0daughters();
831  }
832  else {
833  for (Int_t iProng = 0; iProng<cfVtxHF->GetNProngs(); iProng++){
834  trackCuts[iProng]=fCuts->GetTrackCuts();
835  }
836  }
837 
838  //General settings: vertex, feed down and fill reco container with generated values.
839  cfVtxHF->SetRecoPrimVertex(zPrimVertex);
840  cfVtxHF->SetMCPrimaryVertex(zMCVertex);
842  cfVtxHF->SetNVar(fNvar);
846 
847  // switch-off the trigger class selection (doesn't work for MC)
848  fCuts->SetTriggerClass("");
849 
850  // MC vertex, to be used, in case, for pp
852 
853  if (fCentralitySelection){ // keep only the requested centrality
854 
855  if(fCuts->IsEventSelectedInCentrality(aodEvent)!=0) {
856  delete[] containerInput;
857  delete[] containerInputMC;
858  delete [] trackCuts;
859  delete cfVtxHF;
860  return;
861  }
862  } else { // keep all centralities
863  fCuts->SetMinCentrality(0.);
864  fCuts->SetMaxCentrality(100.);
865  }
866 
867  Float_t centValue = 0.;
868  if(!fIsPPData) centValue = fCuts->GetCentrality(aodEvent);
869  cfVtxHF->SetCentralityValue(centValue);
870 
871  // multiplicity estimator with VZERO
872  Double_t vzeroMult=0;
873  AliAODVZERO *vzeroAOD = (AliAODVZERO*)aodEvent->GetVZEROData();
874  if(vzeroAOD) vzeroMult = vzeroAOD->GetMTotV0A() + vzeroAOD->GetMTotV0C();
875 
876  Double_t multiplicity = nTracklets; // set to the Ntracklet estimator
877  if(fMultiplicityEstimator==kVZERO) { multiplicity = vzeroMult; }
878 
879  cfVtxHF->SetMultiplicity(multiplicity);
880 
881  // printf("Multiplicity estimator %d, value %2.2f\n",fMultiplicityEstimator,multiplicity);
882 
883  for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) {
884  AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart));
885  if (!mcPart){
886  AliError("Failed casting particle from MC array!, Skipping particle");
887  continue;
888  }
889 
890  // Obtain B0 pt spectrum
891  if ((((mcPart->GetPdgCode()%1000)/100==5||(mcPart->GetPdgCode()%1000)/100==-5)) && //||(mcPart->GetPdgCode()%10000)/1000==-5)||(mcPart->GetPdgCode()%10000)/1000==-5)) &&
892  (mcPart->Y()<1.)&&(mcPart->Y()>-1.)) {
893  fhBptCutVar->Fill(mcPart->Pt());
894  }
895 
896  // Obtain the pt spectrum of the injected quarks
897  if (TMath::Abs(mcPart->GetPdgCode())==5&&TMath::Abs(mcPart->Y())<1.) fhBquarkPt->Fill(mcPart->Pt());
898  if (TMath::Abs(mcPart->GetPdgCode())==4&&TMath::Abs(mcPart->Y())<1.) fhCquarkPt->Fill(mcPart->Pt());
899 
900  //counting c quarks
901  cquarks += cfVtxHF->MCcquarkCounting(mcPart);
902 
903  // check the MC-level cuts, must be the desidered particle
904  if (!fCFManager->CheckParticleCuts(0, mcPart)) {
905  AliDebug(2,"Check the MC-level cuts - not desidered particle");
906  continue; // 0 stands for MC level
907  }
908  else {
909  AliDebug(3, Form("\n\n---> COOL! we found a particle (particle %d)!!! with PDG code = %d \n\n", iPart, mcPart->GetPdgCode()));
910  }
911  cfVtxHF->SetMCCandidateParam(iPart);
912 
913 
914  if (!(cfVtxHF->SetLabelArray())){
915  AliDebug(2,Form("Impossible to set the label array for particle %d (decaychannel = %d)", iPart, fDecayChannel));
916  continue;
917  }
918 
919  //check the candiate family at MC level
920  if (!(cfVtxHF->CheckMCPartFamily(mcPart, mcArray))) {
921  AliDebug(2,Form("Check on the family wrong for particle %d!!! (decaychannel = %d)", iPart, fDecayChannel));
922  continue;
923  }
924  else{
925  AliDebug(2,Form("Check on the family OK for particle %d!!! (decaychannel = %d)", iPart, fDecayChannel));
926  }
927 
928  //Fill the MC container
929  Bool_t mcContainerFilled = cfVtxHF -> FillMCContainer(containerInputMC);
930  AliDebug(2, Form("particle = %d mcContainerFilled = %d", iPart, mcContainerFilled));
931  if (mcContainerFilled) {
932  if (fUseWeight){
933  if (fHistoPtWeight) { // using an histogram as weight function
934  AliDebug(2,"Using Histogram as Pt weight function");
935  fWeight = eventWeight*GetPtWeightFromHistogram(containerInputMC[0]);
936  }
937  else if (fFuncWeight){ // using user-defined function
938  AliDebug(2,"Using function");
939  fWeight = eventWeight*fFuncWeight->Eval(containerInputMC[0]);
940  }
941  else{ // using FONLL
942  AliDebug(2,"Using FONLL");
943  fWeight = eventWeight*GetWeight(containerInputMC[0]);
944  }
945  AliDebug(2,Form("pt = %f, weight = %f",containerInputMC[0], fWeight));
946  }
947  if (!fCuts->IsInFiducialAcceptance(containerInputMC[0],containerInputMC[1])) {
948  AliDebug(3, Form("Not in limited acceptance, containerInputMC[0] = %f, containerInputMC[1] = %f", containerInputMC[0], containerInputMC[1]));
949  continue;
950  }
951  else{
952  AliDebug(3, Form("YES!! in limited acceptance, containerInputMC[0] = %f, containerInputMC[1] = %f", containerInputMC[0],containerInputMC[1]));
953  }
954 
955  //MC Limited Acceptance
956  if (TMath::Abs(containerInputMC[1]) < 0.5) {
957  fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepGeneratedLimAcc, fWeight);
958  AliDebug(3,"MC Lim Acc container filled\n");
959  }
960 
961  //MC
962  fCFManager->GetParticleContainer()->Fill(containerInputMC, kStepGenerated, fWeight);
963  icountMC++;
964  AliDebug(3,"MC container filled \n");
965 
966  // MC in acceptance
967  // check the MC-Acceptance level cuts
968  // since standard CF functions are not applicable, using Kine Cuts on daughters
969  Bool_t mcAccepStep = cfVtxHF-> MCAcceptanceStep();
970  if (mcAccepStep){
971  fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepAcceptance, fWeight);
972  AliDebug(3,"MC acceptance cut passed\n");
973  icountAcc++;
974  if(fTHnGenerator){fObjSpr->FillGenStep(mcPart, mcPart->Pt(), fWeight, mcArray, mcHeader);}
975 
976  //MC Vertex step
977  if (fCuts->IsEventSelected(aodEvent)){
978  // filling the container if the vertex is ok
979  fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepVertex, fWeight) ;
980  AliDebug(3,"Vertex cut passed and container filled\n");
981  icountVertex++;
982 
983  //mc Refit requirement
984  Bool_t mcRefitStep = cfVtxHF->MCRefitStep(aodEvent, &trackCuts[0]);
985  if (mcRefitStep){
986  fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepRefit, fWeight);
987  AliDebug(3,"MC Refit cut passed and container filled\n");
988  icountRefit++;
989  }
990  else{
991  AliDebug(3,"MC Refit cut not passed\n");
992  continue;
993  }
994  }
995  else{
996  AliDebug (3, "MC vertex step not passed\n");
997  continue;
998  }
999  }
1000  else{
1001  AliDebug (3, "MC in acceptance step not passed\n");
1002  continue;
1003  }
1004  }
1005  else {
1006  AliDebug (3, "MC container not filled\n");
1007  }
1008  }
1009 
1010  if (cquarks<2) AliDebug(2,Form("Event with %d c-quarks", cquarks));
1011  AliDebug(2,Form("Found %i MC particles that are %s!!",icountMC,fPartName.Data()));
1012  AliDebug(2,Form("Found %i MC particles that are %s and satisfy Acc cuts!!",icountAcc,fPartName.Data()));
1013  AliDebug(2,Form("Found %i MC particles that are %s and satisfy Vertex cuts!!",icountVertex,fPartName.Data()));
1014  AliDebug(2,Form("Found %i MC particles that are %s and satisfy Refit cuts!!",icountRefit,fPartName.Data()));
1015 
1016  // Now go to rec level
1017  fCountMC += icountMC;
1018  fCountAcc += icountAcc;
1019  fCountVertex+= icountVertex;
1020  fCountRefit+= icountRefit;
1021 
1022  AliDebug(2,Form("Found %d vertices for decay channel %d",arrayBranch->GetEntriesFast(),fDecayChannel));
1023 
1024  for(Int_t iCandid = 0; iCandid<arrayBranch->GetEntriesFast();iCandid++){
1025  AliAODRecoDecayHF* charmCandidate=0x0;
1026  switch (fDecayChannel){
1027  case 2:{
1028  charmCandidate = (AliAODRecoDecayHF2Prong*)arrayBranch->At(iCandid);
1029  break;
1030  }
1031  case 21:
1032  case 22:{
1033  charmCandidate = (AliAODRecoCascadeHF*)arrayBranch->At(iCandid);
1034  break;
1035  }
1036  case 31:
1037  case 32:
1038  case 33:{
1039  charmCandidate = (AliAODRecoDecayHF3Prong*)arrayBranch->At(iCandid);
1040  break;
1041  }
1042  case 4:{
1043  charmCandidate = (AliAODRecoDecayHF4Prong*)arrayBranch->At(iCandid);
1044  break;
1045  }
1046  default:
1047  break;
1048  }
1049 
1050  Double_t motherPtWeight=1.;
1052  if (fHistoMotherPtWeight) {
1053  AliDebug(2,"Using Histogram as mother pt weight function");
1054  Double_t motherPt = GetMotherPtFromRecoDecay(charmCandidate,mcArray);
1055  motherPtWeight = GetMotherPtWeightFromHistogram(motherPt);
1056  AliDebug(2,Form("mother pt = %f, weight = %f", motherPt, fWeight));
1057  }
1058  else AliDebug(2, "Couldn't find mother pt dependend weights");
1059  }
1060 
1061  Bool_t unsetvtx=kFALSE;
1062  if(!charmCandidate->GetOwnPrimaryVtx()) {
1063  charmCandidate->SetOwnPrimaryVtx(aodVtx); // needed to compute all variables
1064  unsetvtx=kTRUE;
1065  }
1066 
1067  Bool_t signAssociation = cfVtxHF->SetRecoCandidateParam((AliAODRecoDecayHF*)charmCandidate);
1068  if (!signAssociation){
1069  if(unsetvtx) charmCandidate->UnsetOwnPrimaryVtx();
1070  continue;
1071  }
1072 
1073  Int_t isPartOrAntipart = cfVtxHF->CheckReflexion(fSign);
1074  if (isPartOrAntipart == 0){
1075  AliDebug(2, Form("The candidate pdg code doesn't match the requirement set in the task (fSign = %d)",fSign));
1076  if(unsetvtx) charmCandidate->UnsetOwnPrimaryVtx();
1077  continue;
1078  }
1079 
1080  AliDebug(3,Form("iCandid=%d - signAssociation=%d, isPartOrAntipart=%d",iCandid, signAssociation, isPartOrAntipart));
1081 
1082  Bool_t recoContFilled = cfVtxHF->FillRecoContainer(containerInput);
1083  AliDebug(3, Form("CF task: RecoContFilled for candidate %d is %d", iCandid, (Int_t)recoContFilled));
1084  if (recoContFilled){
1085 
1086  // weight according to pt
1087  if (fUseWeight){
1088  if (fHistoPtWeight) {
1089  AliDebug(2,"Using Histogram as Pt weight function");
1090  fWeight = eventWeight*GetPtWeightFromHistogram(containerInput[0]);
1091  }
1092  else if (fFuncWeight){ // using user-defined function
1093  AliDebug(2, "Using function");
1094  fWeight = eventWeight*fFuncWeight->Eval(containerInput[0]);
1095  }
1096  else{ // using FONLL
1097  AliDebug(2, "Using FONLL");
1098  fWeight = eventWeight*GetWeight(containerInput[0]);
1099  }
1100  AliDebug(2, Form("pt = %f, weight = %f",containerInput[0], fWeight));
1101  }
1102 
1103  if (!fCuts->IsInFiducialAcceptance(containerInput[0],containerInput[1])){
1104  if(unsetvtx) charmCandidate->UnsetOwnPrimaryVtx();
1105  continue;
1106  }
1107  //Reco Step
1108  Bool_t recoStep = cfVtxHF->RecoStep();
1109  if (recoStep) AliDebug(2, Form("particle = %d --> CF task: Reco step for candidate %d is %d", iCandid, iCandid, (Int_t)recoStep));
1110  Bool_t vtxCheck = fCuts->IsEventSelected(aodEvent);
1111 
1112 
1113  // Selection on the filtering bit
1114  Bool_t isBitSelected = kTRUE;
1115  if(fDecayChannel==2) {
1116  if(fUseSelectionBit && !charmCandidate->HasSelectionBit(AliRDHFCuts::kD0toKpiCuts)) isBitSelected = kFALSE;
1117  }else if(fDecayChannel==31){
1118  if(fUseSelectionBit && !charmCandidate->HasSelectionBit(AliRDHFCuts::kDplusCuts)) isBitSelected = kFALSE;
1119  }else if(fDecayChannel==33){
1120  if(fUseSelectionBit && !charmCandidate->HasSelectionBit(AliRDHFCuts::kDsCuts)) isBitSelected = kFALSE;
1121  }else if(fDecayChannel==32){
1122  if(fUseSelectionBit && !charmCandidate->HasSelectionBit(AliRDHFCuts::kLcCuts)) isBitSelected = kFALSE;
1123  }
1124  if(!isBitSelected){
1125  if(unsetvtx) charmCandidate->UnsetOwnPrimaryVtx();
1126  continue;
1127  }
1128 
1129 
1130  if (recoStep && recoContFilled && vtxCheck){
1131  fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed, fWeight) ;
1132  icountReco++;
1133  AliDebug(3,"Reco step passed and container filled\n");
1134 
1135  //Reco in the acceptance -- take care of UNFOLDING!!!!
1136  Bool_t recoAcceptanceStep = cfVtxHF->RecoAcceptStep(&trackCuts[0]);
1137  if (recoAcceptanceStep) {
1138  fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoAcceptance, fWeight) ;
1139  icountRecoAcc++;
1140  AliDebug(3,"Reco acceptance cut passed and container filled\n");
1141 
1142  if(fAcceptanceUnf){
1143  Double_t fill[4]; //fill response matrix
1144  Bool_t bUnfolding = cfVtxHF -> FillUnfoldingMatrix(fPDGcode,fill);
1145  if (bUnfolding) fCorrelation->Fill(fill);
1146  }
1147 
1148  //Number of ITS cluster requirements
1149  Int_t recoITSnCluster = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kTracks);
1150  if (recoITSnCluster){
1151  fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoITSClusters, fWeight) ;
1152  icountRecoITSClusters++;
1153  AliDebug(3,"Reco n ITS cluster cut passed and container filled\n");
1154 
1155  Bool_t iscutsusingpid = fCuts->GetIsUsePID();
1156  Int_t recoAnalysisCuts = -1, recoPidSelection = -1;
1157  fCuts->SetUsePID(kFALSE);
1158  recoAnalysisCuts = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kCandidate, aodEvent);
1159 
1160  if (fDecayChannel==33){ // Ds case, where more possibilities are considered
1161  Bool_t keepDs=ProcessDs(recoAnalysisCuts);
1162  if(keepDs) recoAnalysisCuts=3;
1163  }
1164  else if (fDecayChannel==22){ // Lc->V0+bachelor case, where more possibilities are considered
1165  Bool_t keepLctoV0bachelor = ProcessLctoV0Bachelor(recoAnalysisCuts);
1166  if (keepLctoV0bachelor) recoAnalysisCuts=3;
1167  AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
1168  AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
1169  AliPIDResponse* pidResponse = inputHandler->GetPIDResponse();
1170  if (fUseAdditionalCuts){
1171  if (!((AliCFVertexingHFCascade*)cfVtxHF)->CheckAdditionalCuts(pidResponse)) recoAnalysisCuts = -1;
1172  }
1173  }
1174 
1175  fCuts->SetUsePID(iscutsusingpid); //restoring usage of the PID from the cuts object
1176  Bool_t tempAn=(recoAnalysisCuts == 3 || recoAnalysisCuts == isPartOrAntipart);
1177  if (fDecayChannel == 22) tempAn = (recoAnalysisCuts == 3);
1178  if (fDecayChannel == 32) tempAn=(recoAnalysisCuts >0 || recoAnalysisCuts == isPartOrAntipart);
1179 
1180  if (tempAn){
1181  fCFManager->GetParticleContainer()->Fill(containerInput, kStepRecoPPR, fWeight);
1182  icountRecoPPR++;
1183  AliDebug(3,"Reco Analysis cuts passed and container filled \n");
1184  //pid selection
1185  //recoPidSelection = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kPID);
1186  //if((fCuts->CombineSelectionLevels(3,recoAnalysisCuts,recoPidSelection)==isPartOrAntipart)||(fCuts->CombineSelectionLevels(3,recoAnalysisCuts,recoPidSelection)==3)){
1187  recoPidSelection = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kCandidate, aodEvent);
1188 
1189  if (fDecayChannel==33){ // Ds case, where more possibilities are considered
1190  Bool_t keepDs=ProcessDs(recoPidSelection);
1191  if(keepDs) recoPidSelection=3;
1192  } else if (fDecayChannel==22){ // Lc->V0+bachelor case, where more possibilities are considered
1193  Bool_t keepLctoV0bachelor=ProcessLctoV0Bachelor(recoPidSelection);
1194  if (keepLctoV0bachelor) recoPidSelection=3;
1195  }
1196 
1197  Bool_t tempPid=(recoPidSelection == 3 || recoPidSelection == isPartOrAntipart);
1198  if (fDecayChannel == 22) tempPid = (recoPidSelection == 3);
1199  if (fDecayChannel == 32) tempPid=(recoPidSelection >0 || recoPidSelection == isPartOrAntipart);
1200 
1201  if (tempPid){
1202  Double_t weigPID = 1.;
1203  if (fDecayChannel == 2){ // D0 with Bayesian PID using weights
1204  if(((AliRDHFCutsD0toKpi*)fCuts)->GetCombPID() && (((AliRDHFCutsD0toKpi*)fCuts)->GetBayesianStrategy() == AliRDHFCutsD0toKpi::kBayesWeight || ((AliRDHFCutsD0toKpi*)fCuts)->GetBayesianStrategy() == AliRDHFCutsD0toKpi::kBayesWeightNoFilter)){
1205  if (isPartOrAntipart == 1){
1206  weigPID = ((AliRDHFCutsD0toKpi*)fCuts)->GetWeightsNegative()[AliPID::kKaon] * ((AliRDHFCutsD0toKpi*)fCuts)->GetWeightsPositive()[AliPID::kPion];
1207  }else if (isPartOrAntipart == 2){
1208  weigPID = ((AliRDHFCutsD0toKpi*)fCuts)->GetWeightsPositive()[AliPID::kKaon] * ((AliRDHFCutsD0toKpi*)fCuts)->GetWeightsNegative()[AliPID::kPion];
1209  }
1210  if ((weigPID < 0) || (weigPID > 1)) weigPID = 0.;
1211  }
1212  }else if (fDecayChannel == 33){ // Ds with Bayesian PID using weights
1214  Int_t labDau0=((AliAODTrack*)charmCandidate->GetDaughter(0))->GetLabel();
1215  AliAODMCParticle* firstDau=(AliAODMCParticle*)mcArray->UncheckedAt(TMath::Abs(labDau0));
1216  if(firstDau){
1217  Int_t pdgCode0=TMath::Abs(firstDau->GetPdgCode());
1218  if(pdgCode0==321){
1219  weigPID=((AliRDHFCutsDstoKKpi*)fCuts)->GetWeightForKKpi();
1220  }else if(pdgCode0==211){
1221  weigPID=((AliRDHFCutsDstoKKpi*)fCuts)->GetWeightForpiKK();
1222  }
1223  if ((weigPID < 0) || (weigPID > 1)) weigPID = 0.;
1224  }else{
1225  weigPID=0.;
1226  }
1227  }
1228  }
1229  fCFManager->GetParticleContainer()->Fill(containerInput, kStepRecoPID, fWeight*weigPID);
1230 
1231  AliAODRecoDecayHF2Prong *d0toKpi = (AliAODRecoDecayHF2Prong*)charmCandidate;
1232  fObjSpr->SetFillMC(kTRUE);
1233  fObjSpr->FillSparses(d0toKpi, recoAnalysisCuts, d0toKpi->Pt(),d0toKpi->InvMassD0(), d0toKpi->InvMassD0bar(), fWeight*weigPID, mcArray, aodEvent, mcHeader);
1234 
1235  icountRecoPID++;
1236  AliDebug(3,"Reco PID cuts passed and container filled \n");
1237  if(!fAcceptanceUnf){
1238  Double_t fill[4]; //fill response matrix
1239  Bool_t bUnfolding = cfVtxHF -> FillUnfoldingMatrix(fPDGcode,fill);
1240  if (bUnfolding) fCorrelation->Fill(fill);
1241  }
1242  }
1243  else {
1244  AliDebug(3, "Analysis Cuts step not passed \n");
1245  if(unsetvtx) charmCandidate->UnsetOwnPrimaryVtx();
1246  continue;
1247  }
1248  }
1249  else {
1250  AliDebug(3, "PID selection not passed \n");
1251  if(unsetvtx) charmCandidate->UnsetOwnPrimaryVtx();
1252  continue;
1253  }
1254  }
1255  else{
1256  AliDebug(3, "Number of ITS cluster step not passed\n");
1257  if(unsetvtx) charmCandidate->UnsetOwnPrimaryVtx();
1258  continue;
1259  }
1260  }
1261  else{
1262  AliDebug(3, "Reco acceptance step not passed\n");
1263  if(unsetvtx) charmCandidate->UnsetOwnPrimaryVtx();
1264  continue;
1265  }
1266  }
1267  else {
1268  AliDebug(3, "Reco step not passed\n");
1269  if(unsetvtx) charmCandidate->UnsetOwnPrimaryVtx();
1270  continue;
1271  }
1272  }
1273 
1274  if(unsetvtx) charmCandidate->UnsetOwnPrimaryVtx();
1275  } // end loop on candidate
1276 
1277  fCountReco+= icountReco;
1278  fCountRecoAcc+= icountRecoAcc;
1279  fCountRecoITSClusters+= icountRecoITSClusters;
1280  fCountRecoPPR+= icountRecoPPR;
1281  fCountRecoPID+= icountRecoPID;
1282 
1283  fHistEventsProcessed->Fill(1.5);
1284 
1285  delete[] containerInput;
1286  delete[] containerInputMC;
1287  delete cfVtxHF;
1288  if (trackCuts){
1289  // for (Int_t i=0; i<cfVtxHF->GetNProngs(); i++){
1290  // delete [] trackCuts[i];
1291  // }
1292  delete [] trackCuts;
1293  }
1294 
1295 
1296 }
1297 
1298 //___________________________________________________________________________
1300 {
1301  // The Terminate() function is the last function to be called during
1302  // a query. It always runs on the client, it can be used to present
1303  // the results graphically or save the results to file.
1304 
1305  AliAnalysisTaskSE::Terminate();
1306 
1307  AliInfo(Form("Found %i MC particles that are %s in MC, in %d events",fCountMC,fPartName.Data(),fEvents));
1308  AliInfo(Form("Found %i MC particles that are %s in MC and satisfy Acc cuts, in %d events",fCountAcc,fPartName.Data(),fEvents));
1309  AliInfo(Form("Found %i MC particles that are %s in MC and satisfy Acc cuts, and satisfy Vertex requirement in %d events",fCountVertex,fPartName.Data(),fEvents));
1310  AliInfo(Form("Found %i MC particles that are %s in MC and satisfy Acc cuts, and satisfy ITS+TPC refit requirementin %d events",fCountRefit,fPartName.Data(),fEvents));
1311  AliInfo(Form("Found %i reco %s that are decaying in %s, in %d events",fCountReco,fPartName.Data(),fDauNames.Data(),fEvents));
1312  AliInfo(Form("Among the above, found %i reco %s that are decaying in %s and are in the requested acceptance, in %d events",fCountRecoAcc,fPartName.Data(),fDauNames.Data(),fEvents));
1313  AliInfo(Form("Among the above, found %i reco %s that are decaying in %s and have at least 5 clusters in ITS, in %d events",fCountRecoITSClusters,fPartName.Data(),fDauNames.Data(),fEvents));
1314  AliInfo(Form("Among the above, found %i reco %s that are decaying in %s and satisfy PPR cuts, in %d events",fCountRecoPPR,fPartName.Data(),fDauNames.Data(),fEvents));
1315  AliInfo(Form("Among the above, found %i reco %s that are decaying in %s and satisfy PPR+PID cuts, in %d events",fCountRecoPID,fPartName.Data(),fDauNames.Data(),fEvents));
1316 
1317  // draw some example plots....
1318  AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
1319  if(!cont) {
1320  printf("CONTAINER NOT FOUND\n");
1321  return;
1322  }
1323  // projecting the containers to obtain histograms
1324  // first argument = variable, second argument = step
1325 
1326  TH1D** h = new TH1D*[3];
1327  Int_t nvarToPlot = 0;
1328  if (fConfiguration == kSnail){
1329  //h = new TH1D[3][12];
1330  for (Int_t ih = 0; ih<3; ih++){
1331  if(fDecayChannel==22){
1332  nvarToPlot = 16;
1333  } else {
1334  nvarToPlot = 12;
1335  }
1336  h[ih] = new TH1D[nvarToPlot];
1337  }
1338  for(Int_t iC=1;iC<nvarToPlot; iC++){
1339  // MC-level
1340  h[0][iC] = *(cont->ShowProjection(iC,0));
1341  // MC-Acceptance level
1342  h[1][iC] = *(cont->ShowProjection(iC,1));
1343  // Reco-level
1344  h[2][iC] = *(cont->ShowProjection(iC,4));
1345  }
1346  }
1347  else {
1348  //h = new TH1D[3][12];
1349  nvarToPlot = 8;
1350  for (Int_t ih = 0; ih<3; ih++){
1351  h[ih] = new TH1D[nvarToPlot];
1352  }
1353  for(Int_t iC=0;iC<nvarToPlot; iC++){
1354  // MC-level
1355  h[0][iC] = *(cont->ShowProjection(iC,0));
1356  // MC-Acceptance level
1357  h[1][iC] = *(cont->ShowProjection(iC,1));
1358  // Reco-level
1359  h[2][iC] = *(cont->ShowProjection(iC,4));
1360  }
1361  }
1362  TString* titles;
1363  //Int_t nvarToPlot = 0;
1364  if (fConfiguration == kSnail){
1365  if(fDecayChannel==31){
1366  //nvarToPlot = 12;
1367  titles = new TString[nvarToPlot];
1368  titles[0]="pT_Dplus (GeV/c)";
1369  titles[1]="rapidity";
1370  titles[2]="phi (rad)";
1371  titles[3]="cT (#mum)";
1372  titles[4]="cosPointingAngle";
1373  titles[5]="pT_1 (GeV/c)";
1374  titles[6]="pT_2 (GeV/c)";
1375  titles[7]="pT_3 (GeV/c)";
1376  titles[8]="d0_1 (#mum)";
1377  titles[9]="d0_2 (#mum)";
1378  titles[10]="d0_3 (#mum)";
1379  titles[11]="zVertex (cm)";
1380  } else if (fDecayChannel==22) {
1381  //nvarToPlot = 16;
1382  titles = new TString[nvarToPlot];
1383  titles[0]="p_{T}(#Lambda_{c}) [GeV/c]";
1384  titles[1]="y(#Lambda_{c})";
1385  titles[2]="#varphi(#Lambda_{c}) [rad]";
1386  titles[3]="onTheFlyStatusV0";
1387  titles[4]="z_{vtx} [cm]";
1388  titles[5]="centrality";
1389  titles[6]="fake";
1390  titles[7]="multiplicity";
1391  //titles[8]="pT(bachelor) [GeV/c]";
1392  titles[8]="p(bachelor) [GeV/c]";
1393  titles[9]="p_{T}(V0) [GeV/c]";
1394  titles[10]="y(V0)";
1395  titles[11]="#varphi(V0) [rad]";
1396  titles[12]="m_{inv}(#pi^{+}#pi^{+}) [GeV/c^{2}]";
1397  titles[13]="dcaV0 (nSigma)";
1398  titles[14]="cosine pointing angle (V0)";
1399  titles[15]="cosine pointing angle (#Lambda_{c})";
1400  //titles[16]="c#tauV0 (#mum)";
1401  //titles[17]="c#tau (#mum)";
1402  } else {
1403  //nvarToPlot = 12;
1404  titles = new TString[nvarToPlot];
1405  titles[0]="pT_D0 (GeV/c)";
1406  titles[1]="rapidity";
1407  titles[2]="cosThetaStar";
1408  titles[3]="pT_pi (GeV/c)";
1409  titles[4]="pT_K (Gev/c)";
1410  titles[5]="cT (#mum)";
1411  titles[6]="dca (#mum)";
1412  titles[7]="d0_pi (#mum)";
1413  titles[8]="d0_K (#mum)";
1414  titles[9]="d0xd0 (#mum^2)";
1415  titles[10]="cosPointingAngle";
1416  titles[11]="phi (rad)";
1417  }
1418  }
1419  else {
1420  //nvarToPlot = 8;
1421  titles = new TString[nvarToPlot];
1422  if (fDecayChannel==22) {
1423  titles[0]="p_{T}(#Lambda_{c}) [GeV/c]";
1424  titles[1]="y(#Lambda_{c})";
1425  titles[2]="#varphi(#Lambda_{c}) [rad]";
1426  titles[3]="onTheFlyStatusV0";
1427  titles[4]="z_{vtx} [cm]";
1428  titles[5]="centrality";
1429  titles[6]="fake";
1430  titles[7]="multiplicity";
1431  } else {
1432  titles[0]="pT_candidate (GeV/c)";
1433  titles[1]="rapidity";
1434  titles[2]="cT (#mum)";
1435  titles[3]="phi";
1436  titles[4]="z_{vtx}";
1437  titles[5]="centrality";
1438  titles[6]="fake";
1439  titles[7]="multiplicity";
1440  }
1441  }
1442 
1443  Int_t markers[16]={20,24,21,25,27,28,
1444  20,24,21,25,27,28,
1445  20,24,21,25};
1446  Int_t colors[3]={2,8,4};
1447  for(Int_t iC=0;iC<nvarToPlot; iC++){
1448  for(Int_t iStep=0;iStep<3;iStep++){
1449  h[iStep][iC].SetTitle(titles[iC].Data());
1450  h[iStep][iC].GetXaxis()->SetTitle(titles[iC].Data());
1451  Double_t maxh=h[iStep][iC].GetMaximum();
1452  h[iStep][iC].GetYaxis()->SetRangeUser(0,maxh*1.2);
1453  h[iStep][iC].SetMarkerStyle(markers[iC]);
1454  h[iStep][iC].SetMarkerColor(colors[iStep]);
1455  }
1456  }
1457 
1458  gStyle->SetCanvasColor(0);
1459  gStyle->SetFrameFillColor(0);
1460  gStyle->SetTitleFillColor(0);
1461  gStyle->SetStatColor(0);
1462 
1463  // drawing in 2 separate canvas for a matter of clearity
1464  TCanvas * c1 =new TCanvas(Form("c1New_%d",fDecayChannel),"Vars 0, 1, 2, 3",1100,1200);
1465  c1->Divide(3,4);
1466  Int_t iPad=1;
1467  for(Int_t iVar=0; iVar<4; iVar++){
1468  c1->cd(iPad++);
1469  h[0][iVar].DrawCopy("p");
1470  c1->cd(iPad++);
1471  h[1][iVar].DrawCopy("p");
1472  c1->cd(iPad++);
1473  h[2][iVar].DrawCopy("p");
1474  }
1475 
1476  TCanvas * c2 =new TCanvas(Form("c2New_%d",fDecayChannel),"Vars 4, 5, 6, 7",1100,1200);
1477  c2->Divide(3,4);
1478  iPad=1;
1479  for(Int_t iVar=4; iVar<8; iVar++){
1480  c2->cd(iPad++);
1481  h[0][iVar].DrawCopy("p");
1482  c2->cd(iPad++);
1483  h[1][iVar].DrawCopy("p");
1484  c2->cd(iPad++);
1485  h[2][iVar].DrawCopy("p");
1486  }
1487 
1488  if (fConfiguration == kSnail){
1489  TCanvas * c3 =new TCanvas(Form("c3New_%d",fDecayChannel),"Vars 8, 9, 10, 11",1100,1200);
1490  c3->Divide(3,4);
1491  iPad=1;
1492  for(Int_t iVar=8; iVar<12; iVar++){
1493  c3->cd(iPad++);
1494  h[0][iVar].DrawCopy("p");
1495  c3->cd(iPad++);
1496  h[1][iVar].DrawCopy("p");
1497  c3->cd(iPad++);
1498  h[2][iVar].DrawCopy("p");
1499  }
1500  if (fDecayChannel==22) {
1501  TCanvas * c4 =new TCanvas(Form("c4New_%d",fDecayChannel),"Vars 12, 13, 14, 15",1100,1200);
1502  c4->Divide(3,4);
1503  iPad=1;
1504  for(Int_t iVar=12; iVar<16; iVar++){
1505  c4->cd(iPad++);
1506  h[0][iVar].DrawCopy("p");
1507  c4->cd(iPad++);
1508  h[1][iVar].DrawCopy("p");
1509  c4->cd(iPad++);
1510  h[2][iVar].DrawCopy("p");
1511  }
1512  }
1513  }
1514 
1515  /*
1516  THnSparseD* hcorr = dynamic_cast<THnSparseD*> (GetOutputData(3));
1517 
1518  TH2D* corr1 = hcorr->Projection(0,2);
1519  TH2D* corr2 = hcorr->Projection(1,3);
1520 
1521  TCanvas * c7 =new TCanvas(Form("c7New_%d",fDecayChannel),"",800,400);
1522  c7->Divide(2,1);
1523  c7->cd(1);
1524  corr1->DrawCopy("text");
1525  c7->cd(2);
1526  corr2->DrawCopy("text");
1527  */
1528  TFile* file_projection = new TFile("CFtaskHFprojectionNew.root","RECREATE");
1529 
1530  // corr1->Write();
1531  // corr2->Write();
1532 
1533  for(Int_t iC=0;iC<nvarToPlot; iC++){
1534  for(Int_t iStep=0;iStep<3;iStep++){
1535  h[iStep][iC].Write(Form("Step%d_%s",iStep,titles[iC].Data()));
1536  }
1537  }
1538  file_projection->Close();
1539  for (Int_t ih = 0; ih<3; ih++) delete [] h[ih];
1540  delete [] h;
1541  delete [] titles;
1542 
1543 }
1544 
1545 //___________________________________________________________________________
1547 {
1548  //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
1549  //TO BE SET BEFORE THE EXECUTION OF THE TASK
1550  //
1551  Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
1552 
1553  //slot #1
1554  OpenFile(1);
1555  const char* nameoutput=GetOutputSlot(1)->GetContainer()->GetName();
1556  fHistEventsProcessed = new TH1I(nameoutput,"",2,0,2) ;
1557  fHistEventsProcessed->GetXaxis()->SetBinLabel(1,"Events processed (all)");
1558  fHistEventsProcessed->GetXaxis()->SetBinLabel(2,"Events analyzed (after selection)");
1559 
1560  if(!fObjSpr){
1561  fObjSpr=new AliHFsubtractBFDcuts("feedDownCuts","feedDownCuts");
1562  fObjSpr->InitHistos();
1567  }
1568 
1569  fhBptCutVar = new TH1F("hBptCutVar", "B meson #it{p}_{T} spectrum;#it{p}_{T};Counts (a.u.)",201,0.,50.25);
1570  fhBquarkPt = new TH1F("hBquarkPt", "B quark #it{p}_{T} spectrum;#it{p}_{T};Counts (a.u.)",201,0.,50.25);
1571  fhCquarkPt = new TH1F("hCquarkPt", "C quark #it{p}_{T} spectrum;#it{p}_{T};Counts (a.u.)",201,0.,50.25);
1572 
1573  PostData( 1,fHistEventsProcessed) ;
1574  PostData( 2,fCFManager->GetParticleContainer()) ;
1575  PostData( 3,fCorrelation) ;
1576  PostData( 6,fTHnAnalysis);
1577  PostData( 7,fTHnGenerator);
1578  PostData( 8,fhBptCutVar);
1579  PostData( 9,fhBquarkPt);
1580  PostData(10,fhCquarkPt);
1581  PostData(11,fListBdecays);
1582  PostData(12,fQAHists);
1583 }
1584 
1585 
1586 //_________________________________________________________________________
1588  // ad-hoc weight function from ratio of
1589  // pt spectra from FONLL 2.76 TeV and
1590  // pt spectra from MC production LHC12a17a (PYTHIA Perugia0 with pthard bins)
1591  if(fFuncWeight) delete fFuncWeight;
1592  fFuncWeight=new TF1("funcWeight","[0]+[1]*TMath::Exp(-[2]*x)",0.,50.);
1593  fFuncWeight->SetParameter(0,4.63891e-02);
1594  fFuncWeight->SetParameter(1,1.51674e+01);
1595  fFuncWeight->SetParameter(2,4.09941e-01);
1596  fUseWeight=kTRUE;
1597 }
1598 //_________________________________________________________________________
1600  // ad-hoc weight function from ratio of
1601  // D0 pt spectra in PbPb 2011 0-10% centrality and
1602  // pt spectra from MC production LHC12a17a (PYTHIA Perugia0 with pthard bins)
1603  if(fFuncWeight) delete fFuncWeight;
1604  fFuncWeight=new TF1("funcWeight","[0]+[1]/TMath::Power(x,[2])",0.05,50.);
1605  fFuncWeight->SetParameter(0,1.43116e-02);
1606  fFuncWeight->SetParameter(1,4.37758e+02);
1607  fFuncWeight->SetParameter(2,3.08583);
1608  fUseWeight=kTRUE;
1609 }
1610 
1611 //_________________________________________________________________________
1613  // weight function from the ratio of the LHC12a17b MC
1614  // and FONLL calculations for pp data
1615  if(fFuncWeight) delete fFuncWeight;
1616  fFuncWeight=new TF1("funcWeight","([0]*x)/TMath::Power([2],(1+TMath::Power([3],x/[1])))+[4]*TMath::Exp([5]+[6]*x)+[7]*TMath::Exp([8]*x)",0.15,50.);
1617  fFuncWeight->SetParameters(1.92381e+01, 5.05055e+00, 1.05314e+01, 2.5, 1.88214e-03, 3.44871e+00, -9.74325e-02, 1.97671e+00, -3.21278e-01);
1618  fUseWeight=kTRUE;
1619 }
1620 
1621 //_________________________________________________________________________
1623  // weight function from the ratio of the LHC12a17b MC
1624  // and FONLL calculations for pp data
1625  // corrected by the BAMPS Raa calculation for 30-50% CC
1626  if(fFuncWeight) delete fFuncWeight;
1627  fFuncWeight=new TF1("funcWeight","([0]*x)/TMath::Power([2],(1+TMath::Power([3],x/[1])))+[4]*TMath::Exp([5]+[6]*x)+[7]*TMath::Exp([8]*x)",0.15,50.);
1628  fFuncWeight->SetParameters(6.10443e+00, 1.53487e+00, 1.99474e+00, 2.5, 5.51172e-03, 5.86590e+00, -5.46963e-01, 9.41201e-02, -1.64323e-01);
1629  fUseWeight=kTRUE;
1630 }
1631 
1632 //_________________________________________________________________________
1634  // weight function from the ratio of the LHC13d3 MC
1635  // and FONLL calculations for pp data
1636  if(fFuncWeight) delete fFuncWeight;
1637  fFuncWeight=new TF1("funcWeight","([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.);
1638  fFuncWeight->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);
1639  fUseWeight=kTRUE;
1640 }
1641 
1642 //_________________________________________________________________________
1647 
1648  Double_t weights[] = { 0.958796, 0.94262, 1.00746, 1.00264, 1.0103, 1.00033, 0.970773, 0.925759, 0.873651, 0.827537, 0.793838, 0.770711, 0.753224, 0.740698, 0.729354, 0.720129, 0.71487, 0.711469, 0.706371, 0.704551, 0.705003, 0.702572, 0.704875, 0.706797, 0.70574, 0.710273, 0.714451, 0.716691, 0.720462, 0.719736, 0.722672, 0.724752, 0.731169, 0.734225, 0.741758, 0.749343, 0.741116, 0.747457, 0.746682, 0.754226, 0.763797, 0.765175, 0.762687, 0.772638, 0.767542, 0.779808, 0.782834, 0.781626 };
1649  Int_t nBins=48;
1650  Double_t xMin=0.;
1651  Double_t xMax=24.;
1652  fHistoMotherPtWeight = new TH1F("hMotherPtWeight",";Mother #it{p}_T;Weight",nBins,xMin,xMax);
1653  for (Int_t ix=1;ix<nBins+1;++ix) fHistoMotherPtWeight->SetBinContent(ix, weights[ix]);
1654  fUseMotherPtWeight=kTRUE;
1655 }
1656 
1657 //_________________________________________________________________________
1659  // weight function from the ratio of the LHC10f6a MC
1660  // and FONLL calculations for pp data
1661  if(fFuncWeight) delete fFuncWeight;
1662  fFuncWeight=new TF1("funcWeight","([0]*x)/TMath::Power([2],(1+TMath::Power([3],x/[1])))+[4]*TMath::Exp([5]+[6]*x)+[7]*TMath::Exp([8]*x)",0.15,40.);
1663  fFuncWeight->SetParameters(2.41522e+01,4.92146e+00,6.72495e+00,2.5,6.15361e-03,4.78995e+00,-4.29135e-01,3.99421e-01,-1.57220e-02);
1664  fUseWeight=kTRUE;
1665 }
1666 
1667 //_________________________________________________________________________
1669  // weight function from the ratio of the LHC12a12 MC
1670  // and FONLL calculations for pp data
1671  if(fFuncWeight) delete fFuncWeight;
1672  fFuncWeight=new TF1("funcWeight","([0]*x)/TMath::Power([2],(1+TMath::Power([3],x/[1])))+[4]*TMath::Exp([5]+[6]*x)+[7]*TMath::Exp([8]*x+[9])",0.15,50.);
1673  fFuncWeight->SetParameters(1.31497e+01,3.30503e+00,3.45594e+00,2.5,2.28642e-02,1.42372e+00,2.32892e-04,5.21986e-02,-2.14236e-01,3.86200e+00);
1674  fUseWeight=kTRUE;
1675 }
1676 
1677 //_________________________________________________________________________
1679  // weight function from the ratio of the LHC12a12bis MC
1680  // and FONLL calculations for pp data
1681  if(fFuncWeight) delete fFuncWeight;
1682  fFuncWeight=new TF1("funcWeight","([0]*x)/TMath::Power([2],(1+TMath::Power([3],x/[1])))+[4]*TMath::Exp([5]+[6]*x)+[7]*TMath::Exp([8]*x+[9])",0.15,50.);
1683  fFuncWeight->SetParameters(6.54519e+00,2.74007e+00,2.48325e+00,2.5,1.61113e-01,-5.32546e-01,-3.75916e-04,2.38189e-01,-2.17561e-01,2.35975e+00);
1684  fUseWeight=kTRUE;
1685 }
1686 
1687 //_________________________________________________________________________
1689  // weight function from the ratio of the LHC13e2fix MC
1690  // and FONLL calculations for pp data
1691  if(fFuncWeight) delete fFuncWeight;
1692  fFuncWeight=new TF1("funcWeight","([0]*x)/TMath::Power([2],(1+TMath::Power([3],x/[1])))+[4]*TMath::Exp([5]+[6]*x)+[7]*TMath::Exp([8]*x+[9])",0.15,50.);
1693  fFuncWeight->SetParameters(1.85862e+01,2.48171e+00,3.39356e+00,2.5,1.70426e-02,2.28453e+00,-4.57749e-02,5.84585e-02,-3.19719e-01,4.16789e+00);
1694  fUseWeight=kTRUE;
1695 }
1696 
1697 //________________________________________________________________
1699  // weight function from the ratio of the LHC10f6a MC
1700  // and FONLL calculations for pp data
1701  if(fFuncWeight) delete fFuncWeight;
1702  fFuncWeight=new TF1("funcWeight","([0]*x)/TMath::Power([2],(1+TMath::Power([3],x/[1])))+[4]*TMath::Exp([5]+[6]*x)+[7]*TMath::Exp([8]*x)",0.15,40.);
1703  fFuncWeight->SetParameters(2.77730e+01,4.78942e+00,7.45378e+00,2.5,9.86255e-02,2.30120e+00,-4.16435e-01,3.43770e-01,-2.29380e-02);
1704  fUseWeight=kTRUE;
1705 }
1706 
1707 //________________________________________________________________
1709  // weight function from the ratio of the LHC10f6a MC
1710  // and FONLL calculations for pp data
1711  if(fFuncWeight) delete fFuncWeight;
1712  fFuncWeight=new TF1("funcWeight","([0]*x)/TMath::Power([2],(1+TMath::Power([3],x/[1])))+[4]*TMath::Exp([5]+[6]*x)+[7]*TMath::Exp([8]*x+[9])",0.15,40.);
1713  fFuncWeight->SetParameters(1.34412e+01,3.20068e+00,5.14481e+00,2.5,7.59405e-04,7.51821e+00,-3.93811e-01,2.16849e-02,-3.37768e-02,2.40308e+00);
1714  fUseWeight=kTRUE;
1715 }
1716 
1717 //_________________________________________________________________________
1719 {
1720  //
1721  // calculating the weight to fill the container
1722  //
1723 
1724  // FNOLL central:
1725  // p0 = 1.63297e-01 --> 0.322643
1726  // p1 = 2.96275e+00
1727  // p2 = 2.30301e+00
1728  // p3 = 2.50000e+00
1729 
1730  // PYTHIA
1731  // p0 = 1.85906e-01 --> 0.36609
1732  // p1 = 1.94635e+00
1733  // p2 = 1.40463e+00
1734  // p3 = 2.50000e+00
1735 
1736  Double_t func1[4] = {0.322643,2.96275,2.30301,2.5};
1737  Double_t func2[4] = {0.36609,1.94635,1.40463,2.5};
1738 
1739  Double_t dndpt_func1 = dNdptFit(pt,func1);
1740  if(fUseFlatPtWeight) dndpt_func1 = 1./30.;
1741  Double_t dndpt_func2 = dNdptFit(pt,func2);
1742  AliDebug(2,Form("pt = %f, FONLL = %f, Pythia = %f, ratio = %f",pt,dndpt_func1,dndpt_func2,dndpt_func1/dndpt_func2));
1743  return dndpt_func1/dndpt_func2;
1744 }
1745 
1746 //__________________________________________________________________________________________________
1747 Double_t AliCFTaskVertexingHFCutVarFDSub::dNdptFit(Float_t pt, Double_t* par)
1748 {
1749  //
1750  // calculating dNdpt
1751  //
1752 
1753  Double_t denom = TMath::Power((pt/par[1]), par[3] );
1754  Double_t dNdpt = par[0]*pt/TMath::Power(1.+denom, par[2]);
1755 
1756  return dNdpt;
1757 }
1758 
1759 //_________________________________________________________________________
1761 {
1762  //
1763  // Using an histogram as weight function
1764  // weight = 0 in the range outside the function
1765  //
1766  Double_t weight = 0.0;
1767  Int_t histoNbins = fHistoPtWeight->GetNbinsX();
1768  Int_t histobin = fHistoPtWeight->FindBin(pt);
1769  if( (histobin>0) && (histobin<=histoNbins) ) {
1770  weight = fHistoPtWeight->GetBinContent(histobin);
1771  }
1772 
1773  return weight;
1774 }
1775 
1776 //_________________________________________________________________________
1778 {
1779  //
1782  //
1783  Double_t weight = 0.0;
1784  Int_t histoNbins = fHistoMotherPtWeight->GetNbinsX();
1785  Int_t histobin = fHistoMotherPtWeight->FindBin(motherPt);
1786  if( (histobin>0) && (histobin<=histoNbins) ) {
1787  weight = fHistoMotherPtWeight->GetBinContent(histobin);
1788  }
1789 
1790  return weight;
1791 }
1792 
1793 //_________________________________________________________________________
1795 {
1796  //
1798  //
1799  Double_t motherPt = -1.;
1800  Int_t labDau0=((AliAODTrack*)charmCandidate->GetDaughter(0))->GetLabel();
1801  if (labDau0<0) labDau0*=-1;
1802  AliAODMCParticle* firstDau=(AliAODMCParticle*)mcArray->UncheckedAt(labDau0);
1803  Int_t labMother=firstDau->GetMother();
1804  if (labMother<0) labMother*=-1;
1805  AliAODMCParticle* mother=(AliAODMCParticle*)mcArray->UncheckedAt(labMother);
1806  Int_t pdgMother=mother->GetPdgCode();
1807  while ((pdgMother%1000)/100==4 || (pdgMother%10000)/1000==4) {
1809  labMother=mother->GetMother();
1810  if (labMother<0) labMother*=-1;
1811  mother=(AliAODMCParticle*)mcArray->UncheckedAt(labMother);
1812  pdgMother=mother->GetPdgCode();
1813  }
1814  if ((pdgMother%1000)/100!=5 && (pdgMother%10000)/1000!=5) {
1815  AliDebug(3, "Found strange decay, expected the mother to be a beauty hadron!");
1816  }
1817  else {
1818  motherPt=mother->Pt();
1819  }
1820  return motherPt;
1821 }
1822 
1823 
1824 //__________________________________________________________________________________________________
1825 Double_t AliCFTaskVertexingHFCutVarFDSub::GetZWeight(Float_t z, Int_t runnumber){
1826  //
1827  // calculating the z-vtx weight for the given run range
1828  //
1829 
1830  if(runnumber>146824 || runnumber<146803) return 1.0;
1831 
1832  Double_t func1[3] = {1.0, -0.5, 6.5 };
1833  Double_t func2[3] = {1.0, -0.5, 5.5 };
1834 
1835  Double_t dzFunc1 = DodzFit(z,func1);
1836  Double_t dzFunc2 = DodzFit(z,func2);
1837 
1838  return dzFunc1/dzFunc2;
1839 }
1840 
1841 //__________________________________________________________________________________________________
1842 Double_t AliCFTaskVertexingHFCutVarFDSub::DodzFit(Float_t z, Double_t* par) {
1843 
1844  //
1845  // Gaussian z-vtx shape
1846  //
1847  //gaussian = [0]/TMath::Sqrt(2.*TMath::Pi())/[2]*exp[-(x-[1])*(x-[1])/(2*[2]*[2])]
1848 
1849  Double_t value = par[0]/TMath::Sqrt(2.*TMath::Pi())/par[2]*TMath::Exp(-(z-par[1])*(z-par[1])/2./par[2]/par[2]);
1850 
1851  return value;
1852 }
1853 //__________________________________________________________________________________________________
1855  //
1856  // calculates the Nch weight using the measured and generateed Nch distributions
1857  //
1858  if(nch<=0) return 0.;
1859  if(!fHistoMeasNch || !fHistoMCNch) { AliError("Input histos to evaluate Nch weights missing"); return 0.; }
1860  Double_t pMeas=fHistoMeasNch->GetBinContent(fHistoMeasNch->FindBin(nch));
1861  Double_t pMC=fHistoMCNch->GetBinContent(fHistoMCNch->FindBin(nch));
1862  Double_t weight = pMC>0 ? pMeas/pMC : 0.;
1863  if(fUseMultRatioAsWeight) weight = pMC;
1864  return weight;
1865 }
1866 //__________________________________________________________________________________________________
1868  // creates historgam with measured multiplcity distribution in pp 7 TeV collisions (from Eur. Phys. J. C (2010) 68: 345–354)
1869  //
1870  // for Nch > 70 the points were obtained with a double NBD distribution fit
1871  // TF1 *fit1 = new TF1("fit1","[0]*(TMath::Gamma(x+[1])/(TMath::Gamma(x+1)*TMath::Gamma([1])))*(TMath::Power(([2]/[1]),x))*(TMath::Power((1+([2]/[1])),-x-[1]))"); fit1->SetParameter(0,1.);// normalization constant
1872  // fit1->SetParameter(1,1.63); // k parameter
1873  // fit1->SetParameter(2,12.8); // mean multiplicity
1874  //
1875  Double_t nchbins[82]={0.50,1.50,2.50,3.50,4.50,5.50,6.50,7.50,8.50,9.50,
1876  10.50,11.50,12.50,13.50,14.50,15.50,16.50,17.50,18.50,19.50,
1877  20.50,21.50,22.50,23.50,24.50,25.50,26.50,27.50,28.50,29.50,
1878  30.50,31.50,32.50,33.50,34.50,35.50,36.50,37.50,38.50,39.50,
1879  40.50,41.50,42.50,43.50,44.50,45.50,46.50,47.50,48.50,49.50,
1880  50.50,51.50,52.50,53.50,54.50,55.50,56.50,57.50,58.50,59.50,
1881  60.50,62.50,64.50,66.50,68.50,70.50,72.50,74.50,76.50,78.50,
1882  80.50,82.50,84.50,86.50,88.50,90.50,92.50,94.50,96.50,98.50,
1883  100.50,102.50};
1884  Double_t pch[81]={0.062011,0.072943,0.070771,0.067245,0.062834,0.057383,0.051499,0.04591,0.041109,0.036954,
1885  0.03359,0.030729,0.028539,0.026575,0.024653,0.0229,0.021325,0.019768,0.018561,0.017187,
1886  0.01604,0.014836,0.013726,0.012576,0.011481,0.010393,0.009502,0.008776,0.008024,0.007452,
1887  0.006851,0.006428,0.00594,0.005515,0.005102,0.00469,0.004162,0.003811,0.003389,0.003071,
1888  0.002708,0.002422,0.002184,0.001968,0.00186,0.00165,0.001577,0.001387,0.001254,0.001118,
1889  0.001037,0.000942,0.000823,0.000736,0.000654,0.000579,0.000512,0.00049,0.00045,0.000355,
1890  0.000296,0.000265,0.000193,0.00016,0.000126,0.0000851, 0.0000676,0.0000537,0.0000426, 0.0000338,
1891  0.0000268,0.0000213,0.0000166,0.0000133,0.0000106,0.00000837,0.00000662, 0.00000524,0.00000414, 0.00000327,
1892  0.00000258};
1893 
1894  if(fHistoMeasNch) delete fHistoMeasNch;
1895  fHistoMeasNch=new TH1F("hMeaseNch","",81,nchbins);
1896  for(Int_t i=0; i<81; i++){
1897  fHistoMeasNch->SetBinContent(i+1,pch[i]);
1898  fHistoMeasNch->SetBinError(i+1,0.);
1899  }
1900 }
1901 
1902 //__________________________________________________________________________________________________
1903 Bool_t AliCFTaskVertexingHFCutVarFDSub::ProcessDs(Int_t recoAnalysisCuts) const{
1904  // processes output of Ds is selected
1905  Bool_t keep=kFALSE;
1906  if(recoAnalysisCuts > 0){
1907  Int_t isKKpi=recoAnalysisCuts&1;
1908  Int_t ispiKK=recoAnalysisCuts&2;
1909  Int_t isPhiKKpi=recoAnalysisCuts&4;
1910  Int_t isPhipiKK=recoAnalysisCuts&8;
1911  Int_t isK0starKKpi=recoAnalysisCuts&16;
1912  Int_t isK0starpiKK=recoAnalysisCuts&32;
1913  if(fDsOption==1){
1914  if(isKKpi && isPhiKKpi) keep=kTRUE;
1915  if(ispiKK && isPhipiKK) keep=kTRUE;
1916  }
1917  else if(fDsOption==2){
1918  if(isKKpi && isK0starKKpi) keep=kTRUE;
1919  if(ispiKK && isK0starpiKK) keep=kTRUE;
1920  }
1921  else if(fDsOption==3)keep=kTRUE;
1922  }
1923  return keep;
1924 }
1925 //__________________________________________________________________________________________________
1926 Bool_t AliCFTaskVertexingHFCutVarFDSub::ProcessLctoV0Bachelor(Int_t recoAnalysisCuts) const{
1927  // processes output of Lc->V0+bnachelor is selected
1928  Bool_t keep=kFALSE;
1929  if(recoAnalysisCuts > 0){
1930 
1931  Int_t isK0Sp = recoAnalysisCuts&1;
1932  Int_t isLambdaBarpi = recoAnalysisCuts&2;
1933  Int_t isLambdapi = recoAnalysisCuts&4;
1934 
1935  if(fLctoV0bachelorOption==1){
1936  if(isK0Sp) keep=kTRUE;
1937  }
1938  else if(fLctoV0bachelorOption==2){
1939  if(isLambdaBarpi) keep=kTRUE;
1940  }
1941  else if(fLctoV0bachelorOption==4){
1942  if(isLambdapi) keep=kTRUE;
1943  }
1944  else if(fLctoV0bachelorOption==7) {
1945  if (isK0Sp || isLambdaBarpi || isLambdapi) keep=kTRUE;
1946  }
1947  }
1948  return keep;
1949 }
1950 
1951 //____________________________________________________________________________
1953  // Get Estimator Histogram from period event->GetRunNumber();
1954  //
1955  // If you select SPD tracklets in |eta|<1 you should use type == 1
1956  //
1957 
1958  Int_t runNo = event->GetRunNumber();
1959  Int_t period = -1; // pp: 0-LHC10b, 1-LHC10c, 2-LHC10d, 3-LHC10e
1960  // pPb: 0-LHC13b, 1-LHC13c
1961 
1962  if (fIsPPbData) { // setting run numbers for LHC13 if pPb
1963  if (runNo>195343 && runNo<195484) period = 0;
1964  if (runNo>195528 && runNo<195678) period = 1;
1965  if (period<0 || period>1) return 0;
1966  } else { //else assume pp 2010
1967  if(runNo>114930 && runNo<117223) period = 0;
1968  if(runNo>119158 && runNo<120830) period = 1;
1969  if(runNo>122373 && runNo<126438) period = 2;
1970  if(runNo>127711 && runNo<130841) period = 3;
1971  if(period<0 || period>3) return 0;
1972  }
1973 
1974  return fMultEstimatorAvg[period];
1975 }
Int_t fNvar
flag to use directly the ratio of the distributions (fHistoMCNch) instead of computing it ...
void SetCentralityValue(Float_t centValue)
Bool_t fAcceptanceUnf
flag to select D0 origins. 0 Only from charm 1 only from beauty 2 both from charm and beauty ...
virtual AliESDtrackCuts * GetTrackCutsV0daughters() const
Definition: AliRDHFCuts.h:248
void SetFillFromGenerated(Bool_t flag)
Bool_t fRejectIfNoQuark
selection flag for fakes tracks
Int_t fCountRecoAcc
Reco particle found that satisfy cuts.
virtual AliESDtrackCuts * GetTrackCutsSoftPi() const
Definition: AliRDHFCuts.h:247
void SetUseMCVertex()
Definition: AliRDHFCuts.h:337
Bool_t fFillFromGenerated
decay channel to configure the task
Int_t MCcquarkCounting(AliAODMCParticle *mcPart) const
UShort_t fOriginDselection
flag to indicate whether data container should be filled with generated values also for reconstructed...
Bool_t fUseFlatPtWeight
weight used to fill the container
Class for HF corrections as a function of many variables and steps For D* and other cascades...
Int_t IsEventSelectedInCentrality(AliVEvent *event)
Bool_t HasSelectionBit(Int_t i) const
Class for HF corrections as a function of many variables and step.
Bool_t fIsPPData
flag to use the z-vtx corrected (if not use uncorrected) multiplicity estimator
AliHFsubtractBFDcuts * fObjSpr
cut on momentum conservation
Bool_t FillRecoContainer(Double_t *containerInput)
Double_t fWeight
flag to decide whether to use pt-weights != 1 when filling the container or not
Double_t GetMotherPtFromRecoDecay(AliAODRecoDecayHF *charmCandidate, TClonesArray *mcArray)
Int_t fConfiguration
Ds decay option (generation level)
THnSparseF * GetSparseMC() const
void SetNVar(Int_t nVar)
Int_t fCountRecoITSClusters
Reco particle found that satisfy cuts in requested acceptance.
Bool_t fIsPPbData
flag for pp data (not checking centrality)
Double_t dNdptFit(Float_t pt, Double_t *par)
Class for storing and handling D0 meson candidates properties // for estimating the feed-down fractio...
void SetRejectCandidateIfNotFromQuark(Bool_t opt)
virtual Bool_t SetLabelArray()
Bool_t fCentralitySelection
flag to decide wheter to keep D0 only (0), D0bar only (1), or both D0 and D0bar (2) ...
UInt_t fPDGcode
flag to use selection bit
Bool_t fUseMCVertex
flag to remove events not geenrated with PYTHIA
Int_t GetNProngs() const
TH1F * fHistoMCNch
histogram with measured Nch distribution (pp 7 TeV)
Double_t GetMaxVtxZ() const
Definition: AliRDHFCuts.h:241
Double_t GetMotherPtWeightFromHistogram(Float_t motherPt)
void SetFillMC(Bool_t fillMC=kTRUE)
Int_t fCountMC
list of profile histos for z-vtx correction
AliCFTaskVertexingHFCutVarFDSub()
multiplicity estimators
TH1F * fHistoMeasNch
user-defined histogram to calculate the Mother Pt weights
TH1F * fHistoPtWeight
user-defined function to be used to calculate weights
Class for cuts on AOD reconstructed D+->Kpipi.
TF1 * fFuncWeight
configuration (slow / fast) of the CF –> different variables will be allocated (all / reduced number)...
TH1F * fhBptCutVar
Generator level histogram D0 pt, Nprongs of the decay, Mother pt.
static Int_t GetNumberOfTrackletsInEtaRange(AliAODEvent *ev, Double_t mineta, Double_t maxeta)
Int_t fCountReco
Reco particle found that satisfy kTPCrefit and kITSrefit.
Bool_t ProcessDs(Int_t returnCodeDs) const
void SetMinCentrality(Float_t minCentrality=0.)
Definition: AliRDHFCuts.h:51
THnSparseF * GetSparseMCgen() const
Bool_t fUseAdditionalCuts
flag for pPb data (used for multiplicity corrections)
Bool_t fUseMultRatioAsWeight
flag to decide whether to use Ncharged weights != 1 when filling the container or not ...
Int_t fGenLctoV0bachelorOption
Lc->V0+bachelor decay option (selection level)
TH1F * fHistoMotherPtWeight
user-defined histogram to calculate the Pt weights
Bool_t fUseTrackletsWeight
flag to decide whether to use Ncharged weights != 1 when filling the container or not ...
Double_t DodzFit(Float_t z, Double_t *par)
AliAODVertex * GetOwnPrimaryVtx() const
AliESDtrackCuts * GetTrackCuts() const
Definition: AliRDHFCuts.h:246
Float_t GetCentrality(AliAODEvent *aodEvent)
Definition: AliRDHFCuts.h:242
Int_t fLctoV0bachelorOption
resonant deacy channel to be used if the CF should be run on resonant channels only ...
void SetFakeSelection(Int_t fakeSel)
Bool_t fUseSelectionBit
Lc->V0+bachelor decay option (generation level)
Bool_t fZvtxCorrectedNtrkEstimator
refrence multiplcity (period b)
Int_t fCountRefit
Reco particle found that satisfy vertex constrained.
Double_t fRefMult
TProfile with mult vas. Z per period.
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)
void SetMaxCentrality(Float_t maxCentrality=100.)
Definition: AliRDHFCuts.h:52
void SetOwnPrimaryVtx(const AliAODVertex *vtx)
void SetRecoPrimVertex(Double_t zPrimVertex)
Double_t GetZWeight(Float_t z, Int_t runnumber)
TList * fQAHists
List with B hadron decays (unique entries)
void SetMCCandidateParam(Int_t label)
TH1I * fHistEventsProcessed
pointer to the CF manager
Int_t fCountRecoPPR
Reco particle found that satisfy cuts in n. of ITS clusters.
Bool_t CheckMCPartFamily(AliAODMCParticle *, TClonesArray *) const
Bool_t fUseCutsForTMVA
flag to use additional cuts needed for Lc –> K0S + p, TMVA
virtual Bool_t SetRecoCandidateParam(AliAODRecoDecayHF *)
static Double_t GetCorrectedNtracklets(TProfile *estimatorAvg, Double_t uncorrectedNacc, Double_t vtxZ, Double_t refMult)
Int_t CheckReflexion(Char_t isSign)
Bool_t IsEventSelected(AliVEvent *event)
Bool_t fUseMotherPtWeight
flag to decide whether to use pt-weights != 1 when filling the container or not
Bool_t fUseZWeight
flag to decide to use a flat pt shape
void SetUsePID(Bool_t flag=kTRUE)
Definition: AliRDHFCuts.h:204
virtual void PrintAll() const
void FillSparses(AliAODRecoDecayHF2Prong *dzeroPart, Int_t isSelected, Double_t pt=-1, Double_t massD0=-1, Double_t massD0bar=-1, Double_t weight=1., TClonesArray *mcArray=0x0, AliAODEvent *aodEvent=0x0, AliAODMCHeader *mcHeader=0x0)
AliCFTaskVertexingHFCutVarFDSub & operator=(const AliCFTaskVertexingHFCutVarFDSub &c)
Bool_t fUseNchWeight
flag to decide whether to use z-vtx weights != 1 when filling the container or not ...
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
Bool_t RecoAcceptStep(AliESDtrackCuts **trackCuts) const
void FillGenStep(AliAODMCParticle *dzeroMC, Double_t pt=-1, Double_t weight=1., TClonesArray *mcArray=0x0, AliAODMCHeader *mcHeader=0x0)
TList * fListProfiles
response matrix for unfolding
THnSparseF * fTHnAnalysis
object for cut variation study
static Int_t GetGeneratedPhysicalPrimariesInEtaRange(TClonesArray *arrayMC, Double_t mineta, Double_t maxeta)
Bool_t IsSelected(TObject *obj)
Definition: AliRDHFCuts.h:274
TProfile * fMultEstimatorAvg[4]
Definition of the multiplicity estimator: kNtrk10=0, kNtrk10to16=1, kVZERO=2.
Bool_t MCRefitStep(AliAODEvent *aodEvent, AliESDtrackCuts **trackCuts) const
Int_t fFakeSelection
flag to switch off the centrality selection
Bool_t GetIsUsePID() const
Definition: AliRDHFCuts.h:254
UInt_t fResonantDecay
histogram with Nch distribution from MC production
Float_t fCutOnMomConservation
flag to define which task to use for Lc –> K0S+p
AliRDHFCuts * fCuts
flag for unfolding before or after cuts.
void SetTriggerClass(TString trclass0, TString trclass1="")
Definition: AliRDHFCuts.h:192
Class to compute variables for correction framework // for 3-body decays of D mesons (D+...
virtual Bool_t IsInFiducialAcceptance(Double_t, Double_t) const
Definition: AliRDHFCuts.h:291
TProfile * GetEstimatorHistogram(const AliVEvent *event)
void SetConfiguration(Int_t configuration)
void SetMultiplicity(Double_t multiplicity)
Int_t fGenDsOption
Ds decay option (selection level)
void SetMCPrimaryVertex(Double_t zMCVertex)
Int_t fCountRecoPID
Reco particle found that satisfy cuts in PPR.
Int_t fDsOption
flag to use MC vertex (useful when runnign in pp)
void UserCreateOutputObjects()
ANALYSIS FRAMEWORK STUFF to loop on data and fill output objects.
Bool_t ProcessLctoV0Bachelor(Int_t returnCodeDs) const
TString fPartName
number of variables for the container
Class for HF corrections as a function of many variables and step.
Int_t fCountVertex
MC particle found that satisfy acceptance cuts.