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