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