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