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