AliPhysics  648edd6 (648edd6)
 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 //__________________________________________________________________________
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 //___________________________________________________________________________
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  }else if(fDecayChannel==22){
1124  if(!((dynamic_cast<AliAODRecoCascadeHF*>(charmCandidate))->CheckCascadeFlags())) isBitSelected = kFALSE; // select only Lc among cascade candidates
1125  }
1126  if(!isBitSelected){
1127  if(unsetvtx) charmCandidate->UnsetOwnPrimaryVtx();
1128  continue;
1129  }
1130 
1131 
1132  if (recoStep && recoContFilled && vtxCheck){
1133  fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed, fWeight) ;
1134  icountReco++;
1135  AliDebug(3,"Reco step passed and container filled\n");
1136 
1137  //Reco in the acceptance -- take care of UNFOLDING!!!!
1138  Bool_t recoAcceptanceStep = cfVtxHF->RecoAcceptStep(&trackCuts[0]);
1139  if (recoAcceptanceStep) {
1140  fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoAcceptance, fWeight) ;
1141  icountRecoAcc++;
1142  AliDebug(3,"Reco acceptance cut passed and container filled\n");
1143 
1144  if(fAcceptanceUnf){
1145  Double_t fill[4]; //fill response matrix
1146  Bool_t bUnfolding = cfVtxHF -> FillUnfoldingMatrix(fPDGcode,fill);
1147  if (bUnfolding) fCorrelation->Fill(fill);
1148  }
1149 
1150  //Number of ITS cluster requirements
1151  Int_t recoITSnCluster = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kTracks);
1152  if (recoITSnCluster){
1153  fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoITSClusters, fWeight) ;
1154  icountRecoITSClusters++;
1155  AliDebug(3,"Reco n ITS cluster cut passed and container filled\n");
1156 
1157  Bool_t iscutsusingpid = fCuts->GetIsUsePID();
1158  Int_t recoAnalysisCuts = -1, recoPidSelection = -1;
1159  fCuts->SetUsePID(kFALSE);
1160  recoAnalysisCuts = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kCandidate, aodEvent);
1161 
1162  if (fDecayChannel==33){ // Ds case, where more possibilities are considered
1163  Bool_t keepDs=ProcessDs(recoAnalysisCuts);
1164  if(keepDs) recoAnalysisCuts=3;
1165  }
1166  else if (fDecayChannel==22){ // Lc->V0+bachelor case, where more possibilities are considered
1167  Bool_t keepLctoV0bachelor = ProcessLctoV0Bachelor(recoAnalysisCuts);
1168  if (keepLctoV0bachelor) recoAnalysisCuts=3;
1169  AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
1170  AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
1171  AliPIDResponse* pidResponse = inputHandler->GetPIDResponse();
1172  if (fUseAdditionalCuts){
1173  if (!((AliCFVertexingHFCascade*)cfVtxHF)->CheckAdditionalCuts(pidResponse)) recoAnalysisCuts = -1;
1174  }
1175  }
1176 
1177  fCuts->SetUsePID(iscutsusingpid); //restoring usage of the PID from the cuts object
1178  Bool_t tempAn=(recoAnalysisCuts == 3 || recoAnalysisCuts == isPartOrAntipart);
1179  if (fDecayChannel == 22) tempAn = (recoAnalysisCuts == 3);
1180  if (fDecayChannel == 32) tempAn=(recoAnalysisCuts >0 || recoAnalysisCuts == isPartOrAntipart);
1181 
1182  if (tempAn){
1183  fCFManager->GetParticleContainer()->Fill(containerInput, kStepRecoPPR, fWeight);
1184  icountRecoPPR++;
1185  AliDebug(3,"Reco Analysis cuts passed and container filled \n");
1186  //pid selection
1187  //recoPidSelection = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kPID);
1188  //if((fCuts->CombineSelectionLevels(3,recoAnalysisCuts,recoPidSelection)==isPartOrAntipart)||(fCuts->CombineSelectionLevels(3,recoAnalysisCuts,recoPidSelection)==3)){
1189  recoPidSelection = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kCandidate, aodEvent);
1190 
1191  if (fDecayChannel==33){ // Ds case, where more possibilities are considered
1192  Bool_t keepDs=ProcessDs(recoPidSelection);
1193  if(keepDs) recoPidSelection=3;
1194  } else if (fDecayChannel==22){ // Lc->V0+bachelor case, where more possibilities are considered
1195  Bool_t keepLctoV0bachelor=ProcessLctoV0Bachelor(recoPidSelection);
1196  if (keepLctoV0bachelor) recoPidSelection=3;
1197  }
1198 
1199  Bool_t tempPid=(recoPidSelection == 3 || recoPidSelection == isPartOrAntipart);
1200  if (fDecayChannel == 22) tempPid = (recoPidSelection == 3);
1201  if (fDecayChannel == 32) tempPid=(recoPidSelection >0 || recoPidSelection == isPartOrAntipart);
1202 
1203  if (tempPid){
1204  Double_t weigPID = 1.;
1205  if (fDecayChannel == 2){ // D0 with Bayesian PID using weights
1206  if(((AliRDHFCutsD0toKpi*)fCuts)->GetCombPID() && (((AliRDHFCutsD0toKpi*)fCuts)->GetBayesianStrategy() == AliRDHFCutsD0toKpi::kBayesWeight || ((AliRDHFCutsD0toKpi*)fCuts)->GetBayesianStrategy() == AliRDHFCutsD0toKpi::kBayesWeightNoFilter)){
1207  if (isPartOrAntipart == 1){
1208  weigPID = ((AliRDHFCutsD0toKpi*)fCuts)->GetWeightsNegative()[AliPID::kKaon] * ((AliRDHFCutsD0toKpi*)fCuts)->GetWeightsPositive()[AliPID::kPion];
1209  }else if (isPartOrAntipart == 2){
1210  weigPID = ((AliRDHFCutsD0toKpi*)fCuts)->GetWeightsPositive()[AliPID::kKaon] * ((AliRDHFCutsD0toKpi*)fCuts)->GetWeightsNegative()[AliPID::kPion];
1211  }
1212  if ((weigPID < 0) || (weigPID > 1)) weigPID = 0.;
1213  }
1214  }else if (fDecayChannel == 33){ // Ds with Bayesian PID using weights
1216  Int_t labDau0=((AliAODTrack*)charmCandidate->GetDaughter(0))->GetLabel();
1217  AliAODMCParticle* firstDau=(AliAODMCParticle*)mcArray->UncheckedAt(TMath::Abs(labDau0));
1218  if(firstDau){
1219  Int_t pdgCode0=TMath::Abs(firstDau->GetPdgCode());
1220  if(pdgCode0==321){
1221  weigPID=((AliRDHFCutsDstoKKpi*)fCuts)->GetWeightForKKpi();
1222  }else if(pdgCode0==211){
1223  weigPID=((AliRDHFCutsDstoKKpi*)fCuts)->GetWeightForpiKK();
1224  }
1225  if ((weigPID < 0) || (weigPID > 1)) weigPID = 0.;
1226  }else{
1227  weigPID=0.;
1228  }
1229  }
1230  }
1231  fCFManager->GetParticleContainer()->Fill(containerInput, kStepRecoPID, fWeight*weigPID);
1232 
1233  AliAODRecoDecayHF2Prong *d0toKpi = (AliAODRecoDecayHF2Prong*)charmCandidate;
1234  fObjSpr->SetFillMC(kTRUE);
1235  fObjSpr->FillSparses(d0toKpi, recoAnalysisCuts, d0toKpi->Pt(),d0toKpi->InvMassD0(), d0toKpi->InvMassD0bar(), fWeight*weigPID, mcArray, aodEvent, mcHeader);
1236 
1237  icountRecoPID++;
1238  AliDebug(3,"Reco PID cuts passed and container filled \n");
1239  if(!fAcceptanceUnf){
1240  Double_t fill[4]; //fill response matrix
1241  Bool_t bUnfolding = cfVtxHF -> FillUnfoldingMatrix(fPDGcode,fill);
1242  if (bUnfolding) fCorrelation->Fill(fill);
1243  }
1244  }
1245  else {
1246  AliDebug(3, "Analysis Cuts step not passed \n");
1247  if(unsetvtx) charmCandidate->UnsetOwnPrimaryVtx();
1248  continue;
1249  }
1250  }
1251  else {
1252  AliDebug(3, "PID selection not passed \n");
1253  if(unsetvtx) charmCandidate->UnsetOwnPrimaryVtx();
1254  continue;
1255  }
1256  }
1257  else{
1258  AliDebug(3, "Number of ITS cluster step not passed\n");
1259  if(unsetvtx) charmCandidate->UnsetOwnPrimaryVtx();
1260  continue;
1261  }
1262  }
1263  else{
1264  AliDebug(3, "Reco acceptance step not passed\n");
1265  if(unsetvtx) charmCandidate->UnsetOwnPrimaryVtx();
1266  continue;
1267  }
1268  }
1269  else {
1270  AliDebug(3, "Reco step not passed\n");
1271  if(unsetvtx) charmCandidate->UnsetOwnPrimaryVtx();
1272  continue;
1273  }
1274  }
1275 
1276  if(unsetvtx) charmCandidate->UnsetOwnPrimaryVtx();
1277  } // end loop on candidate
1278 
1279  fCountReco+= icountReco;
1280  fCountRecoAcc+= icountRecoAcc;
1281  fCountRecoITSClusters+= icountRecoITSClusters;
1282  fCountRecoPPR+= icountRecoPPR;
1283  fCountRecoPID+= icountRecoPID;
1284 
1285  fHistEventsProcessed->Fill(1.5);
1286 
1287  delete[] containerInput;
1288  delete[] containerInputMC;
1289  delete cfVtxHF;
1290  if (trackCuts){
1291  // for (Int_t i=0; i<cfVtxHF->GetNProngs(); i++){
1292  // delete [] trackCuts[i];
1293  // }
1294  delete [] trackCuts;
1295  }
1296 
1297 
1298 }
1299 
1300 //___________________________________________________________________________
1302 {
1303  // The Terminate() function is the last function to be called during
1304  // a query. It always runs on the client, it can be used to present
1305  // the results graphically or save the results to file.
1306 
1307  AliAnalysisTaskSE::Terminate();
1308 
1309  AliInfo(Form("Found %i MC particles that are %s in MC, in %d events",fCountMC,fPartName.Data(),fEvents));
1310  AliInfo(Form("Found %i MC particles that are %s in MC and satisfy Acc cuts, in %d events",fCountAcc,fPartName.Data(),fEvents));
1311  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));
1312  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));
1313  AliInfo(Form("Found %i reco %s that are decaying in %s, in %d events",fCountReco,fPartName.Data(),fDauNames.Data(),fEvents));
1314  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));
1315  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));
1316  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));
1317  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));
1318 
1319  // draw some example plots....
1320  AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
1321  if(!cont) {
1322  printf("CONTAINER NOT FOUND\n");
1323  return;
1324  }
1325  // projecting the containers to obtain histograms
1326  // first argument = variable, second argument = step
1327 
1328  TH1D** h = new TH1D*[3];
1329  Int_t nvarToPlot = 0;
1330  if (fConfiguration == kSnail){
1331  //h = new TH1D[3][12];
1332  for (Int_t ih = 0; ih<3; ih++){
1333  if(fDecayChannel==22){
1334  nvarToPlot = 16;
1335  } else {
1336  nvarToPlot = 12;
1337  }
1338  h[ih] = new TH1D[nvarToPlot];
1339  }
1340  for(Int_t iC=1;iC<nvarToPlot; iC++){
1341  // MC-level
1342  h[0][iC] = *(cont->ShowProjection(iC,0));
1343  // MC-Acceptance level
1344  h[1][iC] = *(cont->ShowProjection(iC,1));
1345  // Reco-level
1346  h[2][iC] = *(cont->ShowProjection(iC,4));
1347  }
1348  }
1349  else {
1350  //h = new TH1D[3][12];
1351  nvarToPlot = 8;
1352  for (Int_t ih = 0; ih<3; ih++){
1353  h[ih] = new TH1D[nvarToPlot];
1354  }
1355  for(Int_t iC=0;iC<nvarToPlot; iC++){
1356  // MC-level
1357  h[0][iC] = *(cont->ShowProjection(iC,0));
1358  // MC-Acceptance level
1359  h[1][iC] = *(cont->ShowProjection(iC,1));
1360  // Reco-level
1361  h[2][iC] = *(cont->ShowProjection(iC,4));
1362  }
1363  }
1364  TString* titles;
1365  //Int_t nvarToPlot = 0;
1366  if (fConfiguration == kSnail){
1367  if(fDecayChannel==31){
1368  //nvarToPlot = 12;
1369  titles = new TString[nvarToPlot];
1370  titles[0]="pT_Dplus (GeV/c)";
1371  titles[1]="rapidity";
1372  titles[2]="phi (rad)";
1373  titles[3]="cT (#mum)";
1374  titles[4]="cosPointingAngle";
1375  titles[5]="pT_1 (GeV/c)";
1376  titles[6]="pT_2 (GeV/c)";
1377  titles[7]="pT_3 (GeV/c)";
1378  titles[8]="d0_1 (#mum)";
1379  titles[9]="d0_2 (#mum)";
1380  titles[10]="d0_3 (#mum)";
1381  titles[11]="zVertex (cm)";
1382  } else if (fDecayChannel==22) {
1383  //nvarToPlot = 16;
1384  titles = new TString[nvarToPlot];
1385  titles[0]="p_{T}(#Lambda_{c}) [GeV/c]";
1386  titles[1]="y(#Lambda_{c})";
1387  titles[2]="#varphi(#Lambda_{c}) [rad]";
1388  titles[3]="onTheFlyStatusV0";
1389  titles[4]="z_{vtx} [cm]";
1390  titles[5]="centrality";
1391  titles[6]="fake";
1392  titles[7]="multiplicity";
1393  //titles[8]="pT(bachelor) [GeV/c]";
1394  titles[8]="p(bachelor) [GeV/c]";
1395  titles[9]="p_{T}(V0) [GeV/c]";
1396  titles[10]="y(V0)";
1397  titles[11]="#varphi(V0) [rad]";
1398  titles[12]="m_{inv}(#pi^{+}#pi^{+}) [GeV/c^{2}]";
1399  titles[13]="dcaV0 (nSigma)";
1400  titles[14]="cosine pointing angle (V0)";
1401  titles[15]="cosine pointing angle (#Lambda_{c})";
1402  //titles[16]="c#tauV0 (#mum)";
1403  //titles[17]="c#tau (#mum)";
1404  } else {
1405  //nvarToPlot = 12;
1406  titles = new TString[nvarToPlot];
1407  titles[0]="pT_D0 (GeV/c)";
1408  titles[1]="rapidity";
1409  titles[2]="cosThetaStar";
1410  titles[3]="pT_pi (GeV/c)";
1411  titles[4]="pT_K (Gev/c)";
1412  titles[5]="cT (#mum)";
1413  titles[6]="dca (#mum)";
1414  titles[7]="d0_pi (#mum)";
1415  titles[8]="d0_K (#mum)";
1416  titles[9]="d0xd0 (#mum^2)";
1417  titles[10]="cosPointingAngle";
1418  titles[11]="phi (rad)";
1419  }
1420  }
1421  else {
1422  //nvarToPlot = 8;
1423  titles = new TString[nvarToPlot];
1424  if (fDecayChannel==22) {
1425  titles[0]="p_{T}(#Lambda_{c}) [GeV/c]";
1426  titles[1]="y(#Lambda_{c})";
1427  titles[2]="#varphi(#Lambda_{c}) [rad]";
1428  titles[3]="onTheFlyStatusV0";
1429  titles[4]="z_{vtx} [cm]";
1430  titles[5]="centrality";
1431  titles[6]="fake";
1432  titles[7]="multiplicity";
1433  } else {
1434  titles[0]="pT_candidate (GeV/c)";
1435  titles[1]="rapidity";
1436  titles[2]="cT (#mum)";
1437  titles[3]="phi";
1438  titles[4]="z_{vtx}";
1439  titles[5]="centrality";
1440  titles[6]="fake";
1441  titles[7]="multiplicity";
1442  }
1443  }
1444 
1445  Int_t markers[16]={20,24,21,25,27,28,
1446  20,24,21,25,27,28,
1447  20,24,21,25};
1448  Int_t colors[3]={2,8,4};
1449  for(Int_t iC=0;iC<nvarToPlot; iC++){
1450  for(Int_t iStep=0;iStep<3;iStep++){
1451  h[iStep][iC].SetTitle(titles[iC].Data());
1452  h[iStep][iC].GetXaxis()->SetTitle(titles[iC].Data());
1453  Double_t maxh=h[iStep][iC].GetMaximum();
1454  h[iStep][iC].GetYaxis()->SetRangeUser(0,maxh*1.2);
1455  h[iStep][iC].SetMarkerStyle(markers[iC]);
1456  h[iStep][iC].SetMarkerColor(colors[iStep]);
1457  }
1458  }
1459 
1460  gStyle->SetCanvasColor(0);
1461  gStyle->SetFrameFillColor(0);
1462  gStyle->SetTitleFillColor(0);
1463  gStyle->SetStatColor(0);
1464 
1465  // drawing in 2 separate canvas for a matter of clearity
1466  TCanvas * c1 =new TCanvas(Form("c1New_%d",fDecayChannel),"Vars 0, 1, 2, 3",1100,1200);
1467  c1->Divide(3,4);
1468  Int_t iPad=1;
1469  for(Int_t iVar=0; iVar<4; iVar++){
1470  c1->cd(iPad++);
1471  h[0][iVar].DrawCopy("p");
1472  c1->cd(iPad++);
1473  h[1][iVar].DrawCopy("p");
1474  c1->cd(iPad++);
1475  h[2][iVar].DrawCopy("p");
1476  }
1477 
1478  TCanvas * c2 =new TCanvas(Form("c2New_%d",fDecayChannel),"Vars 4, 5, 6, 7",1100,1200);
1479  c2->Divide(3,4);
1480  iPad=1;
1481  for(Int_t iVar=4; iVar<8; iVar++){
1482  c2->cd(iPad++);
1483  h[0][iVar].DrawCopy("p");
1484  c2->cd(iPad++);
1485  h[1][iVar].DrawCopy("p");
1486  c2->cd(iPad++);
1487  h[2][iVar].DrawCopy("p");
1488  }
1489 
1490  if (fConfiguration == kSnail){
1491  TCanvas * c3 =new TCanvas(Form("c3New_%d",fDecayChannel),"Vars 8, 9, 10, 11",1100,1200);
1492  c3->Divide(3,4);
1493  iPad=1;
1494  for(Int_t iVar=8; iVar<12; iVar++){
1495  c3->cd(iPad++);
1496  h[0][iVar].DrawCopy("p");
1497  c3->cd(iPad++);
1498  h[1][iVar].DrawCopy("p");
1499  c3->cd(iPad++);
1500  h[2][iVar].DrawCopy("p");
1501  }
1502  if (fDecayChannel==22) {
1503  TCanvas * c4 =new TCanvas(Form("c4New_%d",fDecayChannel),"Vars 12, 13, 14, 15",1100,1200);
1504  c4->Divide(3,4);
1505  iPad=1;
1506  for(Int_t iVar=12; iVar<16; iVar++){
1507  c4->cd(iPad++);
1508  h[0][iVar].DrawCopy("p");
1509  c4->cd(iPad++);
1510  h[1][iVar].DrawCopy("p");
1511  c4->cd(iPad++);
1512  h[2][iVar].DrawCopy("p");
1513  }
1514  }
1515  }
1516 
1517  /*
1518  THnSparseD* hcorr = dynamic_cast<THnSparseD*> (GetOutputData(3));
1519 
1520  TH2D* corr1 = hcorr->Projection(0,2);
1521  TH2D* corr2 = hcorr->Projection(1,3);
1522 
1523  TCanvas * c7 =new TCanvas(Form("c7New_%d",fDecayChannel),"",800,400);
1524  c7->Divide(2,1);
1525  c7->cd(1);
1526  corr1->DrawCopy("text");
1527  c7->cd(2);
1528  corr2->DrawCopy("text");
1529  */
1530  TFile* file_projection = new TFile("CFtaskHFprojectionNew.root","RECREATE");
1531 
1532  // corr1->Write();
1533  // corr2->Write();
1534 
1535  for(Int_t iC=0;iC<nvarToPlot; iC++){
1536  for(Int_t iStep=0;iStep<3;iStep++){
1537  h[iStep][iC].Write(Form("Step%d_%s",iStep,titles[iC].Data()));
1538  }
1539  }
1540  file_projection->Close();
1541  for (Int_t ih = 0; ih<3; ih++) delete [] h[ih];
1542  delete [] h;
1543  delete [] titles;
1544 
1545 }
1546 
1547 //___________________________________________________________________________
1549 {
1550  //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
1551  //TO BE SET BEFORE THE EXECUTION OF THE TASK
1552  //
1553  Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
1554 
1555  //slot #1
1556  OpenFile(1);
1557  const char* nameoutput=GetOutputSlot(1)->GetContainer()->GetName();
1558  fHistEventsProcessed = new TH1I(nameoutput,"",2,0,2) ;
1559  fHistEventsProcessed->GetXaxis()->SetBinLabel(1,"Events processed (all)");
1560  fHistEventsProcessed->GetXaxis()->SetBinLabel(2,"Events analyzed (after selection)");
1561 
1562  if(!fObjSpr){
1563  fObjSpr=new AliHFsubtractBFDcuts("feedDownCuts","feedDownCuts");
1564  fObjSpr->InitHistos();
1569  }
1570 
1571  fhBptCutVar = new TH1F("hBptCutVar", "B meson #it{p}_{T} spectrum;#it{p}_{T};Counts (a.u.)",201,0.,50.25);
1572  fhBquarkPt = new TH1F("hBquarkPt", "B quark #it{p}_{T} spectrum;#it{p}_{T};Counts (a.u.)",201,0.,50.25);
1573  fhCquarkPt = new TH1F("hCquarkPt", "C quark #it{p}_{T} spectrum;#it{p}_{T};Counts (a.u.)",201,0.,50.25);
1574 
1575  PostData( 1,fHistEventsProcessed) ;
1576  PostData( 2,fCFManager->GetParticleContainer()) ;
1577  PostData( 3,fCorrelation) ;
1578  PostData( 6,fTHnAnalysis);
1579  PostData( 7,fTHnGenerator);
1580  PostData( 8,fhBptCutVar);
1581  PostData( 9,fhBquarkPt);
1582  PostData(10,fhCquarkPt);
1583  PostData(11,fListBdecays);
1584  PostData(12,fQAHists);
1585 }
1586 
1587 
1588 //_________________________________________________________________________
1590  // ad-hoc weight function from ratio of
1591  // pt spectra from FONLL 2.76 TeV and
1592  // pt spectra from MC production LHC12a17a (PYTHIA Perugia0 with pthard bins)
1593  if(fFuncWeight) delete fFuncWeight;
1594  fFuncWeight=new TF1("funcWeight","[0]+[1]*TMath::Exp(-[2]*x)",0.,50.);
1595  fFuncWeight->SetParameter(0,4.63891e-02);
1596  fFuncWeight->SetParameter(1,1.51674e+01);
1597  fFuncWeight->SetParameter(2,4.09941e-01);
1598  fUseWeight=kTRUE;
1599 }
1600 //_________________________________________________________________________
1602  // ad-hoc weight function from ratio of
1603  // D0 pt spectra in PbPb 2011 0-10% centrality and
1604  // pt spectra from MC production LHC12a17a (PYTHIA Perugia0 with pthard bins)
1605  if(fFuncWeight) delete fFuncWeight;
1606  fFuncWeight=new TF1("funcWeight","[0]+[1]/TMath::Power(x,[2])",0.05,50.);
1607  fFuncWeight->SetParameter(0,1.43116e-02);
1608  fFuncWeight->SetParameter(1,4.37758e+02);
1609  fFuncWeight->SetParameter(2,3.08583);
1610  fUseWeight=kTRUE;
1611 }
1612 
1613 //_________________________________________________________________________
1615  // weight function from the ratio of the LHC12a17b MC
1616  // and FONLL calculations for pp data
1617  if(fFuncWeight) delete fFuncWeight;
1618  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.);
1619  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);
1620  fUseWeight=kTRUE;
1621 }
1622 
1623 //_________________________________________________________________________
1625  // weight function from the ratio of the LHC12a17b MC
1626  // and FONLL calculations for pp data
1627  // corrected by the BAMPS Raa calculation for 30-50% CC
1628  if(fFuncWeight) delete fFuncWeight;
1629  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.);
1630  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);
1631  fUseWeight=kTRUE;
1632 }
1633 
1634 //_________________________________________________________________________
1636  // weight function from the ratio of the LHC13d3 MC
1637  // and FONLL calculations for pp data
1638  if(fFuncWeight) delete fFuncWeight;
1639  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.);
1640  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);
1641  fUseWeight=kTRUE;
1642 }
1643 
1644 //_________________________________________________________________________
1649 
1650  const Int_t nBins=48;
1651 
1652  Double_t weights[nBins] = { 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 };
1653 
1654  Double_t xMin=0.;
1655  Double_t xMax=24.;
1656  fHistoMotherPtWeight = new TH1F("hMotherPtWeight",";Mother #it{p}_T;Weight",nBins,xMin,xMax);
1657  for (Int_t ix=1;ix<nBins+1;++ix) fHistoMotherPtWeight->SetBinContent(ix, weights[ix-1]);
1658  fUseMotherPtWeight=kTRUE;
1659 }
1660 
1661 //_________________________________________________________________________
1663  // weight function from the ratio of the LHC10f6a MC
1664  // and FONLL calculations for pp data
1665  if(fFuncWeight) delete fFuncWeight;
1666  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.);
1667  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);
1668  fUseWeight=kTRUE;
1669 }
1670 
1671 //_________________________________________________________________________
1673  // weight function from the ratio of the LHC12a12 MC
1674  // and FONLL calculations for pp data
1675  if(fFuncWeight) delete fFuncWeight;
1676  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.);
1677  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);
1678  fUseWeight=kTRUE;
1679 }
1680 
1681 //_________________________________________________________________________
1683  // weight function from the ratio of the LHC12a12bis MC
1684  // and FONLL calculations for pp data
1685  if(fFuncWeight) delete fFuncWeight;
1686  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.);
1687  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);
1688  fUseWeight=kTRUE;
1689 }
1690 
1691 //_________________________________________________________________________
1693  // weight function from the ratio of the LHC13e2fix MC
1694  // and FONLL calculations for pp data
1695  if(fFuncWeight) delete fFuncWeight;
1696  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.);
1697  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);
1698  fUseWeight=kTRUE;
1699 }
1700 
1701 //________________________________________________________________
1703  // weight function from the ratio of the LHC10f6a MC
1704  // and FONLL calculations for pp data
1705  if(fFuncWeight) delete fFuncWeight;
1706  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.);
1707  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);
1708  fUseWeight=kTRUE;
1709 }
1710 
1711 //________________________________________________________________
1713  // weight function from the ratio of the LHC10f6a MC
1714  // and FONLL calculations for pp data
1715  if(fFuncWeight) delete fFuncWeight;
1716  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.);
1717  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);
1718  fUseWeight=kTRUE;
1719 }
1720 
1721 //_________________________________________________________________________
1723 {
1724  //
1725  // calculating the weight to fill the container
1726  //
1727 
1728  // FNOLL central:
1729  // p0 = 1.63297e-01 --> 0.322643
1730  // p1 = 2.96275e+00
1731  // p2 = 2.30301e+00
1732  // p3 = 2.50000e+00
1733 
1734  // PYTHIA
1735  // p0 = 1.85906e-01 --> 0.36609
1736  // p1 = 1.94635e+00
1737  // p2 = 1.40463e+00
1738  // p3 = 2.50000e+00
1739 
1740  Double_t func1[4] = {0.322643,2.96275,2.30301,2.5};
1741  Double_t func2[4] = {0.36609,1.94635,1.40463,2.5};
1742 
1743  Double_t dndpt_func1 = dNdptFit(pt,func1);
1744  if(fUseFlatPtWeight) dndpt_func1 = 1./30.;
1745  Double_t dndpt_func2 = dNdptFit(pt,func2);
1746  AliDebug(2,Form("pt = %f, FONLL = %f, Pythia = %f, ratio = %f",pt,dndpt_func1,dndpt_func2,dndpt_func1/dndpt_func2));
1747  return dndpt_func1/dndpt_func2;
1748 }
1749 
1750 //__________________________________________________________________________________________________
1752 {
1753  //
1754  // calculating dNdpt
1755  //
1756 
1757  Double_t denom = TMath::Power((pt/par[1]), par[3] );
1758  Double_t dNdpt = par[0]*pt/TMath::Power(1.+denom, par[2]);
1759 
1760  return dNdpt;
1761 }
1762 
1763 //_________________________________________________________________________
1765 {
1766  //
1767  // Using an histogram as weight function
1768  // weight = 0 in the range outside the function
1769  //
1770  Double_t weight = 0.0;
1771  Int_t histoNbins = fHistoPtWeight->GetNbinsX();
1772  Int_t histobin = fHistoPtWeight->FindBin(pt);
1773  if( (histobin>0) && (histobin<=histoNbins) ) {
1774  weight = fHistoPtWeight->GetBinContent(histobin);
1775  }
1776 
1777  return weight;
1778 }
1779 
1780 //_________________________________________________________________________
1782 {
1783  //
1786  //
1787  Double_t weight = 0.0;
1788  Int_t histoNbins = fHistoMotherPtWeight->GetNbinsX();
1789  Int_t histobin = fHistoMotherPtWeight->FindBin(motherPt);
1790  if( (histobin>0) && (histobin<=histoNbins) ) {
1791  weight = fHistoMotherPtWeight->GetBinContent(histobin);
1792  }
1793 
1794  return weight;
1795 }
1796 
1797 //_________________________________________________________________________
1799 {
1800  //
1802  //
1803  Double_t motherPt = -1.;
1804  Int_t labDau0=((AliAODTrack*)charmCandidate->GetDaughter(0))->GetLabel();
1805  if (labDau0<0) labDau0*=-1;
1806  AliAODMCParticle* firstDau=(AliAODMCParticle*)mcArray->UncheckedAt(labDau0);
1807  Int_t labMother=firstDau->GetMother();
1808  if (labMother<0) labMother*=-1;
1809  AliAODMCParticle* mother=(AliAODMCParticle*)mcArray->UncheckedAt(labMother);
1810  Int_t pdgMother=mother->GetPdgCode();
1811  while ((pdgMother%1000)/100==4 || (pdgMother%10000)/1000==4) {
1813  labMother=mother->GetMother();
1814  if (labMother<0) labMother*=-1;
1815  mother=(AliAODMCParticle*)mcArray->UncheckedAt(labMother);
1816  pdgMother=mother->GetPdgCode();
1817  }
1818  if ((pdgMother%1000)/100!=5 && (pdgMother%10000)/1000!=5) {
1819  AliDebug(3, "Found strange decay, expected the mother to be a beauty hadron!");
1820  }
1821  else {
1822  motherPt=mother->Pt();
1823  }
1824  return motherPt;
1825 }
1826 
1827 
1828 //__________________________________________________________________________________________________
1830  //
1831  // calculating the z-vtx weight for the given run range
1832  //
1833 
1834  if(runnumber>146824 || runnumber<146803) return 1.0;
1835 
1836  Double_t func1[3] = {1.0, -0.5, 6.5 };
1837  Double_t func2[3] = {1.0, -0.5, 5.5 };
1838 
1839  Double_t dzFunc1 = DodzFit(z,func1);
1840  Double_t dzFunc2 = DodzFit(z,func2);
1841 
1842  return dzFunc1/dzFunc2;
1843 }
1844 
1845 //__________________________________________________________________________________________________
1847 
1848  //
1849  // Gaussian z-vtx shape
1850  //
1851  //gaussian = [0]/TMath::Sqrt(2.*TMath::Pi())/[2]*exp[-(x-[1])*(x-[1])/(2*[2]*[2])]
1852 
1853  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]);
1854 
1855  return value;
1856 }
1857 //__________________________________________________________________________________________________
1859  //
1860  // calculates the Nch weight using the measured and generateed Nch distributions
1861  //
1862  if(nch<=0) return 0.;
1863  if(!fHistoMeasNch || !fHistoMCNch) { AliError("Input histos to evaluate Nch weights missing"); return 0.; }
1864  Double_t pMeas=fHistoMeasNch->GetBinContent(fHistoMeasNch->FindBin(nch));
1865  Double_t pMC=fHistoMCNch->GetBinContent(fHistoMCNch->FindBin(nch));
1866  Double_t weight = pMC>0 ? pMeas/pMC : 0.;
1867  if(fUseMultRatioAsWeight) weight = pMC;
1868  return weight;
1869 }
1870 //__________________________________________________________________________________________________
1872  // creates historgam with measured multiplcity distribution in pp 7 TeV collisions (from Eur. Phys. J. C (2010) 68: 345–354)
1873  //
1874  // for Nch > 70 the points were obtained with a double NBD distribution fit
1875  // 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
1876  // fit1->SetParameter(1,1.63); // k parameter
1877  // fit1->SetParameter(2,12.8); // mean multiplicity
1878  //
1879  Double_t nchbins[82]={0.50,1.50,2.50,3.50,4.50,5.50,6.50,7.50,8.50,9.50,
1880  10.50,11.50,12.50,13.50,14.50,15.50,16.50,17.50,18.50,19.50,
1881  20.50,21.50,22.50,23.50,24.50,25.50,26.50,27.50,28.50,29.50,
1882  30.50,31.50,32.50,33.50,34.50,35.50,36.50,37.50,38.50,39.50,
1883  40.50,41.50,42.50,43.50,44.50,45.50,46.50,47.50,48.50,49.50,
1884  50.50,51.50,52.50,53.50,54.50,55.50,56.50,57.50,58.50,59.50,
1885  60.50,62.50,64.50,66.50,68.50,70.50,72.50,74.50,76.50,78.50,
1886  80.50,82.50,84.50,86.50,88.50,90.50,92.50,94.50,96.50,98.50,
1887  100.50,102.50};
1888  Double_t pch[81]={0.062011,0.072943,0.070771,0.067245,0.062834,0.057383,0.051499,0.04591,0.041109,0.036954,
1889  0.03359,0.030729,0.028539,0.026575,0.024653,0.0229,0.021325,0.019768,0.018561,0.017187,
1890  0.01604,0.014836,0.013726,0.012576,0.011481,0.010393,0.009502,0.008776,0.008024,0.007452,
1891  0.006851,0.006428,0.00594,0.005515,0.005102,0.00469,0.004162,0.003811,0.003389,0.003071,
1892  0.002708,0.002422,0.002184,0.001968,0.00186,0.00165,0.001577,0.001387,0.001254,0.001118,
1893  0.001037,0.000942,0.000823,0.000736,0.000654,0.000579,0.000512,0.00049,0.00045,0.000355,
1894  0.000296,0.000265,0.000193,0.00016,0.000126,0.0000851, 0.0000676,0.0000537,0.0000426, 0.0000338,
1895  0.0000268,0.0000213,0.0000166,0.0000133,0.0000106,0.00000837,0.00000662, 0.00000524,0.00000414, 0.00000327,
1896  0.00000258};
1897 
1898  if(fHistoMeasNch) delete fHistoMeasNch;
1899  fHistoMeasNch=new TH1F("hMeaseNch","",81,nchbins);
1900  for(Int_t i=0; i<81; i++){
1901  fHistoMeasNch->SetBinContent(i+1,pch[i]);
1902  fHistoMeasNch->SetBinError(i+1,0.);
1903  }
1904 }
1905 
1906 //__________________________________________________________________________________________________
1908  // processes output of Ds is selected
1909  Bool_t keep=kFALSE;
1910  if(recoAnalysisCuts > 0){
1911  Int_t isKKpi=recoAnalysisCuts&1;
1912  Int_t ispiKK=recoAnalysisCuts&2;
1913  Int_t isPhiKKpi=recoAnalysisCuts&4;
1914  Int_t isPhipiKK=recoAnalysisCuts&8;
1915  Int_t isK0starKKpi=recoAnalysisCuts&16;
1916  Int_t isK0starpiKK=recoAnalysisCuts&32;
1917  if(fDsOption==1){
1918  if(isKKpi && isPhiKKpi) keep=kTRUE;
1919  if(ispiKK && isPhipiKK) keep=kTRUE;
1920  }
1921  else if(fDsOption==2){
1922  if(isKKpi && isK0starKKpi) keep=kTRUE;
1923  if(ispiKK && isK0starpiKK) keep=kTRUE;
1924  }
1925  else if(fDsOption==3)keep=kTRUE;
1926  }
1927  return keep;
1928 }
1929 //__________________________________________________________________________________________________
1931  // processes output of Lc->V0+bnachelor is selected
1932  Bool_t keep=kFALSE;
1933  if(recoAnalysisCuts > 0){
1934 
1935  Int_t isK0Sp = recoAnalysisCuts&1;
1936  Int_t isLambdaBarpi = recoAnalysisCuts&2;
1937  Int_t isLambdapi = recoAnalysisCuts&4;
1938 
1939  if(fLctoV0bachelorOption==1){
1940  if(isK0Sp) keep=kTRUE;
1941  }
1942  else if(fLctoV0bachelorOption==2){
1943  if(isLambdaBarpi) keep=kTRUE;
1944  }
1945  else if(fLctoV0bachelorOption==4){
1946  if(isLambdapi) keep=kTRUE;
1947  }
1948  else if(fLctoV0bachelorOption==7) {
1949  if (isK0Sp || isLambdaBarpi || isLambdapi) keep=kTRUE;
1950  }
1951  }
1952  return keep;
1953 }
1954 
1955 //____________________________________________________________________________
1957  // Get Estimator Histogram from period event->GetRunNumber();
1958  //
1959  // If you select SPD tracklets in |eta|<1 you should use type == 1
1960  //
1961 
1962  Int_t runNo = event->GetRunNumber();
1963  Int_t period = -1; // pp: 0-LHC10b, 1-LHC10c, 2-LHC10d, 3-LHC10e
1964  // pPb: 0-LHC13b, 1-LHC13c
1965 
1966  if (fIsPPbData) { // setting run numbers for LHC13 if pPb
1967  if (runNo>195343 && runNo<195484) period = 0;
1968  if (runNo>195528 && runNo<195678) period = 1;
1969  if (period<0 || period>1) return 0;
1970  } else { //else assume pp 2010
1971  if(runNo>114930 && runNo<117223) period = 0;
1972  if(runNo>119158 && runNo<120830) period = 1;
1973  if(runNo>122373 && runNo<126438) period = 2;
1974  if(runNo>127711 && runNo<130841) period = 3;
1975  if(period<0 || period>3) return 0;
1976  }
1977 
1978  return fMultEstimatorAvg[period];
1979 }
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:254
void SetFillFromGenerated(Bool_t flag)
Bool_t fRejectIfNoQuark
selection flag for fakes tracks
Int_t fCountRecoAcc
Reco particle found that satisfy cuts.
double Double_t
Definition: External.C:58
virtual AliESDtrackCuts * GetTrackCutsSoftPi() const
Definition: AliRDHFCuts.h:253
void SetUseMCVertex()
Definition: AliRDHFCuts.h:353
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)
char Char_t
Definition: External.C:18
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)
TCanvas * c
Definition: TestFitELoss.C:172
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:247
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.
int Int_t
Definition: External.C:63
Bool_t ProcessDs(Int_t returnCodeDs) const
Definition: External.C:204
void SetMinCentrality(Float_t minCentrality=0.)
Definition: AliRDHFCuts.h:51
THnSparseF * GetSparseMCgen() const
float Float_t
Definition: External.C:68
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)
Definition: External.C:212
AliAODVertex * GetOwnPrimaryVtx() const
AliESDtrackCuts * GetTrackCuts() const
Definition: AliRDHFCuts.h:252
Float_t GetCentrality(AliAODEvent *aodEvent)
Definition: AliRDHFCuts.h:248
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
Int_t colors[nPtBins]
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)
ClassImp(AliAnalysisTaskDeltaPt) AliAnalysisTaskDeltaPt
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:206
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 ...
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:281
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:260
const char Option_t
Definition: External.C:48
UInt_t fResonantDecay
histogram with Nch distribution from MC production
Float_t fCutOnMomConservation
flag to define which task to use for Lc –> K0S+p
bool Bool_t
Definition: External.C:53
AliRDHFCuts * fCuts
flag for unfolding before or after cuts.
void SetTriggerClass(TString trclass0, TString trclass1="")
Definition: AliRDHFCuts.h:194
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:301
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.
TList * OpenFile(const char *fname)
Definition: DrawAnaELoss.C:65
Int_t fCountVertex
MC particle found that satisfy acceptance cuts.