Changes for PbPb (Zaida)
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliHFPtSpectrum.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2010, 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 /* $Id$ */
17
18 //***********************************************************************
19 // Class AliHFPtSpectrum
20 // Base class for feed-down corrections on heavy-flavour decays
21 // computes the cross-section via one of the three implemented methods:
22 //   0) Consider no feed-down prediction 
23 //   1) Subtract the feed-down with the "fc" method 
24 //       Yield = Reco * fc;  where fc = 1 / ( 1 + (eff_b/eff_c)*(N_b/N_c) ) ;
25 //   2) Subtract the feed-down with the "Nb" method
26 //       Yield = Reco - Feed-down (exact formula on the function implementation)
27 //
28 //  (the corrected yields per bin are divided by the bin-width)
29 //
30 // Author: Z.Conesa, zconesa@in2p3.fr
31 //***********************************************************************
32
33 #include <Riostream.h>
34
35 #include "TMath.h"
36 #include "TH1.h"
37 #include "TH1D.h"
38 #include "TGraphAsymmErrors.h"
39 #include "TNamed.h"
40 #include "TCanvas.h"
41 #include "TLegend.h"
42
43 #include "AliLog.h"
44 #include "AliHFSystErr.h"
45 #include "AliHFPtSpectrum.h"
46
47 ClassImp(AliHFPtSpectrum)
48
49 //_________________________________________________________________________________________________________
50 AliHFPtSpectrum::AliHFPtSpectrum(const char* name, const char* title, Int_t option):
51   TNamed(name,title),
52   fhDirectMCpt(NULL),
53   fhFeedDownMCpt(NULL),
54   fhDirectMCptMax(NULL),
55   fhDirectMCptMin(NULL),
56   fhFeedDownMCptMax(NULL),
57   fhFeedDownMCptMin(NULL),
58   fhDirectEffpt(NULL),
59   fhFeedDownEffpt(NULL),
60   fhRECpt(NULL),
61   fgRECSystematics(NULL),
62   fLuminosity(),
63   fTrigEfficiency(),
64   fGlobalEfficiencyUncertainties(),
65   fhFc(NULL),
66   fhFcMax(NULL),
67   fhFcMin(NULL),
68   fgFcExtreme(NULL),
69   fgFcConservative(NULL),
70   fhYieldCorr(NULL),
71   fhYieldCorrMax(NULL),
72   fhYieldCorrMin(NULL),
73   fgYieldCorr(NULL),
74   fgYieldCorrExtreme(NULL),
75   fgYieldCorrConservative(NULL),
76   fhSigmaCorr(NULL),
77   fhSigmaCorrMax(NULL),
78   fhSigmaCorrMin(NULL),
79   fgSigmaCorr(NULL),
80   fgSigmaCorrExtreme(NULL),
81   fgSigmaCorrConservative(NULL),
82   fFeedDownOption(option),
83   fAsymUncertainties(kTRUE)
84 {
85   //
86   // Default constructor
87   //
88
89   fLuminosity[0]=1.;  fLuminosity[1]=0.;  
90   fTrigEfficiency[0]=1.; fTrigEfficiency[1]=0.; 
91   fGlobalEfficiencyUncertainties[0]=0.; fGlobalEfficiencyUncertainties[1]=0.;
92
93 }
94
95 //_________________________________________________________________________________________________________
96 AliHFPtSpectrum::AliHFPtSpectrum(const AliHFPtSpectrum &rhs):
97   TNamed(rhs),
98   fhDirectMCpt(rhs.fhDirectMCpt),
99   fhFeedDownMCpt(rhs.fhFeedDownMCpt),
100   fhDirectMCptMax(rhs.fhDirectMCptMax),
101   fhDirectMCptMin(rhs.fhDirectMCptMin),
102   fhFeedDownMCptMax(rhs.fhFeedDownMCptMax),
103   fhFeedDownMCptMin(rhs.fhFeedDownMCptMin),
104   fhDirectEffpt(rhs.fhDirectEffpt),
105   fhFeedDownEffpt(rhs.fhFeedDownEffpt),
106   fhRECpt(rhs.fhRECpt),
107   fgRECSystematics(rhs.fgRECSystematics),
108   fLuminosity(),
109   fTrigEfficiency(),
110   fGlobalEfficiencyUncertainties(),
111   fhFc(rhs.fhFc),
112   fhFcMax(rhs.fhFcMax),
113   fhFcMin(rhs.fhFcMin),
114   fgFcExtreme(rhs.fgFcExtreme),
115   fgFcConservative(rhs.fgFcConservative),
116   fhYieldCorr(rhs.fhYieldCorr),
117   fhYieldCorrMax(rhs.fhYieldCorrMax),
118   fhYieldCorrMin(rhs.fhYieldCorrMin),
119   fgYieldCorr(rhs.fgYieldCorr),
120   fgYieldCorrExtreme(rhs.fgYieldCorrExtreme),
121   fgYieldCorrConservative(rhs.fgYieldCorrConservative),
122   fhSigmaCorr(rhs.fhSigmaCorr),
123   fhSigmaCorrMax(rhs.fhSigmaCorrMax),
124   fhSigmaCorrMin(rhs.fhSigmaCorrMin),
125   fgSigmaCorr(rhs.fgSigmaCorr),
126   fgSigmaCorrExtreme(rhs.fgSigmaCorrExtreme),
127   fgSigmaCorrConservative(rhs.fgSigmaCorrConservative),
128   fFeedDownOption(rhs.fFeedDownOption),
129   fAsymUncertainties(rhs.fAsymUncertainties)
130 {
131   //
132   // Copy constructor
133   //
134
135   for(Int_t i=0; i<2; i++){
136     fLuminosity[i] = rhs.fLuminosity[i];
137     fTrigEfficiency[i] = rhs.fTrigEfficiency[i];
138     fGlobalEfficiencyUncertainties[i] = rhs.fGlobalEfficiencyUncertainties[i];
139   }
140
141 }
142
143 //_________________________________________________________________________________________________________
144 AliHFPtSpectrum &AliHFPtSpectrum::operator=(const AliHFPtSpectrum &source){
145   //
146   // Assignment operator
147   //
148
149   if (&source == this) return *this;
150   
151   fhDirectMCpt = source.fhDirectMCpt;
152   fhFeedDownMCpt = source.fhFeedDownMCpt;
153   fhDirectMCptMax = source.fhDirectMCptMax;
154   fhDirectMCptMin = source.fhDirectMCptMin;
155   fhFeedDownMCptMax = source.fhFeedDownMCptMax;
156   fhFeedDownMCptMin = source.fhFeedDownMCptMin;
157   fhDirectEffpt = source.fhDirectEffpt;
158   fhFeedDownEffpt = source.fhFeedDownEffpt;
159   fhRECpt = source.fhRECpt;
160   fgRECSystematics = source.fgRECSystematics;
161   fhFc = source.fhFc;
162   fhFcMax = source.fhFcMax;
163   fhFcMin = source.fhFcMin;
164   fgFcExtreme = source.fgFcExtreme;
165   fgFcConservative = source.fgFcConservative;
166   fhYieldCorr = source.fhYieldCorr;
167   fhYieldCorrMax = source.fhYieldCorrMax;
168   fhYieldCorrMin = source.fhYieldCorrMin;
169   fgYieldCorr = source.fgYieldCorr;
170   fgYieldCorrExtreme = source.fgYieldCorrExtreme;
171   fgYieldCorrConservative = source.fgYieldCorrConservative;
172   fhSigmaCorr = source.fhSigmaCorr;
173   fhSigmaCorrMax = source.fhSigmaCorrMax;
174   fhSigmaCorrMin = source.fhSigmaCorrMin;
175   fgSigmaCorr = source.fgSigmaCorr;
176   fgSigmaCorrExtreme = source.fgSigmaCorrExtreme;
177   fgSigmaCorrConservative = source.fgSigmaCorrConservative;
178   fFeedDownOption = source.fFeedDownOption;
179   fAsymUncertainties = source.fAsymUncertainties;
180   
181   for(Int_t i=0; i<2; i++){
182     fLuminosity[i] = source.fLuminosity[i];
183     fTrigEfficiency[i] = source.fTrigEfficiency[i];
184     fGlobalEfficiencyUncertainties[i] = source.fGlobalEfficiencyUncertainties[i];
185   }
186
187   return *this;
188 }
189
190 //_________________________________________________________________________________________________________
191 AliHFPtSpectrum::~AliHFPtSpectrum(){
192   //
193   // Destructor
194   //
195   if (fhDirectMCpt) delete fhDirectMCpt;    
196   if (fhFeedDownMCpt) delete fhFeedDownMCpt;  
197   if (fhDirectMCptMax) delete fhDirectMCptMax; 
198   if (fhDirectMCptMin) delete fhDirectMCptMin; 
199   if (fhFeedDownMCptMax) delete fhFeedDownMCptMax;
200   if (fhFeedDownMCptMin) delete fhFeedDownMCptMin;
201   if (fhDirectEffpt) delete fhDirectEffpt;    
202   if (fhFeedDownEffpt) delete fhFeedDownEffpt;  
203   if (fhRECpt) delete fhRECpt;    
204   if (fgRECSystematics) delete fgRECSystematics;
205   if (fhFc) delete fhFc;
206   if (fhFcMax) delete fhFcMax;
207   if (fhFcMin) delete fhFcMin;
208   if (fgFcExtreme) delete fgFcExtreme;  
209   if (fgFcConservative) delete fgFcConservative;
210   if (fhYieldCorr) delete fhYieldCorr;                
211   if (fhYieldCorrMax) delete fhYieldCorrMax;             
212   if (fhYieldCorrMin) delete fhYieldCorrMin;             
213   if (fgYieldCorr) delete fgYieldCorr;  
214   if (fgYieldCorrExtreme) delete fgYieldCorrExtreme;
215   if (fgYieldCorrConservative) delete fgYieldCorrConservative;
216   if (fhSigmaCorr) delete fhSigmaCorr;                  
217   if (fhSigmaCorrMax) delete fhSigmaCorrMax;               
218   if (fhSigmaCorrMin) delete fhSigmaCorrMin;               
219   if (fgSigmaCorr) delete fgSigmaCorr;    
220   if (fgSigmaCorrExtreme) delete fgSigmaCorrExtreme;
221   if (fgSigmaCorrConservative) delete fgSigmaCorrConservative;
222 }
223   
224
225 //_________________________________________________________________________________________________________
226 TH1D * AliHFPtSpectrum::RebinTheoreticalSpectra(TH1D *hTheory, const char *name) {
227   //
228   // Function to rebin the theoretical spectrum 
229   //  with respect to the real-data reconstructed spectrum binning 
230   //
231   
232   if (!hTheory || !fhRECpt) {
233     AliError("Feed-down or reconstructed spectra don't exist");
234     return NULL;
235   }
236
237   //
238   // Get the reconstructed spectra bins & limits
239   Int_t nbins = fhRECpt->GetNbinsX();
240   Int_t nbinsMC = hTheory->GetNbinsX();
241   Double_t *limits = new Double_t[nbins+1];
242   Double_t xlow=0., binwidth=0.;
243   for (Int_t i=1; i<=nbins; i++) {
244     binwidth = fhRECpt->GetBinWidth(i);
245     xlow = fhRECpt->GetBinLowEdge(i);
246     limits[i-1] = xlow;
247   }
248   limits[nbins] = xlow + binwidth;
249
250   // Check that the reconstructed spectra binning 
251   // is larger than the theoretical one
252   Double_t thbinwidth = hTheory->GetBinWidth(1);
253   for (Int_t i=1; i<=nbins; i++) {
254     binwidth = fhRECpt->GetBinWidth(i);
255     if ( thbinwidth > binwidth ) {
256       AliInfo(" Beware it seems that the reconstructed spectra has a smaller binning than the theoretical predictions !! ");
257     }
258   }
259
260   //
261   // Define a new histogram with the real-data reconstructed spectrum binning 
262   TH1D * hTheoryRebin = new TH1D(name," theoretical rebinned prediction",nbins,limits);
263
264   Double_t sum[nbins], items[nbins];
265   for (Int_t ibin=0; ibin<nbins; ibin++) {
266     sum[ibin]=0.; items[ibin]=0.;
267   }
268   for (Int_t ibin=0; ibin<=nbinsMC; ibin++){
269     
270     for (Int_t ibinrec=0; ibinrec<nbins; ibinrec++){
271       if (hTheory->GetBinCenter(ibin)>limits[ibinrec] && 
272           hTheory->GetBinCenter(ibin)<limits[ibinrec+1]){
273         sum[ibinrec]+=hTheory->GetBinContent(ibin);
274         items[ibinrec]+=1.;
275       }
276     }
277     
278   }
279
280   // set the theoretical rebinned spectra to ( sum-bins / n-bins ) per new bin
281   for (Int_t ibinrec=0; ibinrec<nbins; ibinrec++) {
282     hTheoryRebin->SetBinContent(ibinrec+1,sum[ibinrec]/items[ibinrec]);
283   }
284   
285   return (TH1D*)hTheoryRebin;
286 }
287
288 //_________________________________________________________________________________________________________
289 void AliHFPtSpectrum::SetMCptSpectra(TH1D *hDirect, TH1D *hFeedDown){
290   //
291   // Set the MonteCarlo or Theoretical spectra
292   //  both for direct and feed-down contributions
293   //
294   
295   if (!hDirect || !hFeedDown || !fhRECpt) {
296     AliError("One or both (direct, feed-down) spectra or the reconstructed spectra don't exist");
297     return;
298   }
299
300   Bool_t areconsistent = kTRUE;
301   areconsistent = CheckHistosConsistency(hDirect,hFeedDown);
302   if (!areconsistent) {
303     AliInfo("Histograms are not consistent (bin width, bounds)"); 
304     return;
305   }
306
307   //
308   // Rebin the theoretical predictions to the reconstructed spectra binning 
309   //
310   fhDirectMCpt = RebinTheoreticalSpectra(hDirect,"fhDirectMCpt");
311   fhDirectMCpt->SetNameTitle("fhDirectMCpt"," direct theoretical prediction");
312   fhFeedDownMCpt = RebinTheoreticalSpectra(hFeedDown,"fhFeedDownMCpt");
313   fhFeedDownMCpt->SetNameTitle("fhFeedDownMCpt"," feed-down theoretical prediction");
314
315 }
316
317 //_________________________________________________________________________________________________________
318 void AliHFPtSpectrum::SetFeedDownMCptSpectra(TH1D *hFeedDown){
319   //
320   // Set the MonteCarlo or Theoretical spectra
321   //  for feed-down contribution
322   //
323   
324   if (!hFeedDown || !fhRECpt) {
325     AliError("Feed-down or reconstructed spectra don't exist");
326     return;
327   }
328
329   //
330   // Rebin the theoretical predictions to the reconstructed spectra binning 
331   //
332   fhFeedDownMCpt = RebinTheoreticalSpectra(hFeedDown,"fhFeedDownMCpt");
333   fhFeedDownMCpt->SetNameTitle("fhFeedDownMCpt"," feed-down theoretical prediction");
334
335 }
336
337 //_________________________________________________________________________________________________________
338 void AliHFPtSpectrum::SetMCptDistributionsBounds(TH1D *hDirectMax, TH1D *hDirectMin, TH1D *hFeedDownMax, TH1D *hFeedDownMin){
339   //
340   // Set the maximum and minimum MonteCarlo or Theoretical spectra
341   //  both for direct and feed-down contributions
342   // used in case uncertainties are asymmetric and ca not be on the "basic histograms"
343   //
344
345   if (!hDirectMax || !hDirectMin || !hFeedDownMax|| !hFeedDownMin || !fhRECpt) {
346     AliError("One or all of the max/min direct/feed-down or the reconstructed spectra don't exist");
347     return;
348   }
349
350   Bool_t areconsistent = kTRUE; 
351   areconsistent &= CheckHistosConsistency(hDirectMax,hDirectMin);
352   areconsistent &= CheckHistosConsistency(hFeedDownMax,hFeedDownMin);
353   areconsistent &= CheckHistosConsistency(hDirectMax,hFeedDownMax);
354   if (!areconsistent) {
355     AliInfo("Histograms are not consistent (bin width, bounds)"); 
356     return;
357   }
358
359
360   //
361   // Rebin the theoretical predictions to the reconstructed spectra binning 
362   //
363   fhDirectMCptMax = RebinTheoreticalSpectra(hDirectMax,"fhDirectMCptMax");
364   fhDirectMCptMax->SetNameTitle("fhDirectMCptMax"," maximum direct theoretical prediction");
365   fhDirectMCptMin = RebinTheoreticalSpectra(hDirectMin,"fhDirectMCptMin");
366   fhDirectMCptMin->SetNameTitle("fhDirectMCptMin"," minimum direct theoretical prediction");
367   fhFeedDownMCptMax = RebinTheoreticalSpectra(hFeedDownMax,"fhFeedDownMCptMax");
368   fhFeedDownMCptMax->SetNameTitle("fhFeedDownMCptMax"," maximum feed-down theoretical prediction");
369   fhFeedDownMCptMin = RebinTheoreticalSpectra(hFeedDownMin,"fhFeedDownMCptMin");
370   fhFeedDownMCptMin->SetNameTitle("fhFeedDownMCptMin"," minimum feed-down theoretical prediction");
371
372 }
373
374 //_________________________________________________________________________________________________________
375 void AliHFPtSpectrum::SetFeedDownMCptDistributionsBounds(TH1D *hFeedDownMax, TH1D *hFeedDownMin){
376   //
377   // Set the maximum and minimum MonteCarlo or Theoretical spectra
378   //   for feed-down contributions
379   // used in case uncertainties are asymmetric and can not be on the "basic histogram"
380   //
381
382   if (!hFeedDownMax || !hFeedDownMin || !fhRECpt) {
383     AliError("One or all of the max/min direct/feed-down spectra don't exist");
384     return;
385   }
386
387   Bool_t areconsistent = kTRUE; 
388   areconsistent &= CheckHistosConsistency(hFeedDownMax,hFeedDownMin);
389   if (!areconsistent) {
390     AliInfo("Histograms are not consistent (bin width, bounds)"); 
391     return;
392   }
393
394
395   //
396   // Rebin the theoretical predictions to the reconstructed spectra binning 
397   //
398   fhFeedDownMCptMax = RebinTheoreticalSpectra(hFeedDownMax,"fhFeedDownMCptMax");
399   fhFeedDownMCptMax->SetNameTitle("fhFeedDownMCptMax"," maximum feed-down theoretical prediction");
400   fhFeedDownMCptMin = RebinTheoreticalSpectra(hFeedDownMin,"fhFeedDownMCptMin");
401   fhFeedDownMCptMin->SetNameTitle("fhFeedDownMCptMin"," minimum feed-down theoretical prediction");
402
403 }
404
405 //_________________________________________________________________________________________________________
406 void AliHFPtSpectrum::SetDirectAccEffCorrection(TH1D *hDirectEff){
407   //
408   // Set the Acceptance and Efficiency corrections 
409   //   for the direct contribution
410   //
411   
412   if (!hDirectEff) {
413     AliError("The direct acceptance and efficiency corrections doesn't exist");
414     return;
415   }
416
417   fhDirectEffpt = (TH1D*)hDirectEff->Clone();
418   fhDirectEffpt->SetNameTitle("fhDirectEffpt"," direct acceptance x efficiency correction");
419 }
420
421 //_________________________________________________________________________________________________________
422 void AliHFPtSpectrum::SetAccEffCorrection(TH1D *hDirectEff, TH1D *hFeedDownEff){
423   //
424   // Set the Acceptance and Efficiency corrections 
425   //  both for direct and feed-down contributions
426   //
427   
428   if (!hDirectEff || !hFeedDownEff) {
429     AliError("One or both (direct, feed-down) acceptance and efficiency corrections don't exist");
430     return;
431   }
432
433   Bool_t areconsistent=kTRUE;
434   areconsistent = CheckHistosConsistency(hDirectEff,hFeedDownEff);
435   if (!areconsistent) {
436     AliInfo("Histograms are not consistent (bin width, bounds)"); 
437     return;
438   }
439
440   fhDirectEffpt = (TH1D*)hDirectEff->Clone();
441   fhFeedDownEffpt = (TH1D*)hFeedDownEff->Clone();
442   fhDirectEffpt->SetNameTitle("fhDirectEffpt"," direct acceptance x efficiency correction");
443   fhFeedDownEffpt->SetNameTitle("fhFeedDownEffpt"," feed-down acceptance x efficiency correction");
444 }
445
446 //_________________________________________________________________________________________________________
447 void AliHFPtSpectrum::SetReconstructedSpectrum(TH1D *hRec) {
448   //
449   // Set the reconstructed spectrum
450   //
451   
452   if (!hRec) {
453     AliError("The reconstructed spectrum doesn't exist");
454     return;
455   }
456
457   fhRECpt = (TH1D*)hRec->Clone();
458   fhRECpt->SetNameTitle("fhRECpt"," reconstructed spectrum");
459 }
460
461 //_________________________________________________________________________________________________________
462 void AliHFPtSpectrum::SetReconstructedSpectrumSystematics(TGraphAsymmErrors *gRec) {
463   //
464   // Set the reconstructed spectrum (uncorrected yield) systematic uncertainties
465   // 
466
467   // Check the compatibility with the reconstructed spectrum
468   Double_t gbinwidth = gRec->GetErrorXlow(1) + gRec->GetErrorXhigh(1) ;
469   Double_t hbinwidth = fhRECpt->GetBinWidth(1);
470   Double_t gxbincenter=0., gybincenter=0.;
471   gRec->GetPoint(1,gxbincenter,gybincenter);
472   Double_t hbincenter = fhRECpt->GetBinCenter(1);
473   if ( (gbinwidth != hbinwidth) || (gxbincenter!=hbincenter) ) {
474     AliError(" The reconstructed spectrum and its systematics don't seem compatible");
475     return;
476   }
477   
478   fgRECSystematics = gRec; 
479 }
480
481 //_________________________________________________________________________________________________________
482 void AliHFPtSpectrum::ComputeHFPtSpectrum(Double_t deltaY, Double_t branchingRatioC, Double_t branchingRatioBintoFinalDecay) {
483   //
484   // Main function to compute the corrected cross-section:
485   // variables : analysed delta_y, BR for the final correction,
486   //             BR b --> D --> decay (relative to the input theoretical prediction)
487   //
488   //   Sigma = ( 1. / (lumi * delta_y * BR_c * eff_trig * eff_c ) ) * spectra (corrected for feed-down)
489   //
490   // Uncertainties: (stat) delta_sigma = sigma * sqrt ( (delta_spectra/spectra)^2 )
491   //  (syst but feed-down) delta_sigma = sigma * sqrt ( (delta_spectra_syst/spectra)^2 + (delta_lumi/lumi)^2 + (delta_eff_trig/eff_trig)^2 + (delta_eff/eff)^2 )
492   //      (feed-down syst) delta_sigma = sigma * sqrt ( (delta_spectra_fd/spectra_fd)^2 )
493
494   //
495   // First: Initialization
496   //
497   Bool_t areHistosOk = Initialize();
498   if (!areHistosOk) {
499     AliInfo(" Histos not properly initialized. Check : inconsistent binning ? missing histos ?");
500     return;
501   }
502
503   //
504   // Second: Correct for feed-down
505   //
506   if (fFeedDownOption==1) {
507     // Compute the feed-down correction via fc-method
508     CalculateFeedDownCorrectionFc(); 
509     // Correct the yield for feed-down correction via fc-method
510     CalculateFeedDownCorrectedSpectrumFc(); 
511   }
512   else if (fFeedDownOption==2) {
513     // Correct the yield for feed-down correction via Nb-method
514     CalculateFeedDownCorrectedSpectrumNb(deltaY,branchingRatioBintoFinalDecay); 
515   }
516   else if (fFeedDownOption==0) { 
517     // If there is no need for feed-down correction,
518     //    the "corrected" yield is equal to the raw yield
519     fhYieldCorr = (TH1D*)fhRECpt->Clone();
520     fhYieldCorr->SetNameTitle("fhYieldCorr","un-corrected yield");
521     fhYieldCorrMax = (TH1D*)fhRECpt->Clone();
522     fhYieldCorrMin = (TH1D*)fhRECpt->Clone();
523     fhYieldCorrMax->SetNameTitle("fhYieldCorrMax","un-corrected yield");
524     fhYieldCorrMin->SetNameTitle("fhYieldCorrMin","un-corrected yield");
525     fAsymUncertainties=kFALSE;
526   }
527   else { 
528     AliInfo(" Are you sure the feed-down correction option is right ?"); 
529   }
530
531   // Print out information
532   //  printf("\n\n     Correcting the spectra with : \n   luminosity = %2.2e +- %2.2e, trigger efficiency = %2.2e +- %2.2e, \n    delta_y = %2.2f, BR_c = %2.2e, BR_b_decay = %2.2e \n\n",fLuminosity[0],fLuminosity[1],fTrigEfficiency[0],fTrigEfficiency[1],deltaY,branchingRatioC,branchingRatioBintoFinalDecay);
533   printf("\n\n     Correcting the spectra with : \n   luminosity = %2.2e +- %2.2e, trigger efficiency = %2.2e +- %2.2e, \n    delta_y = %2.2f, BR_c = %2.2e, BR_b_decay = %2.2e \n    %2.2f percent uncertainty on the efficiencies, and %2.2f percent uncertainty on the b/c efficiencies ratio \n\n",fLuminosity[0],fLuminosity[1],fTrigEfficiency[0],fTrigEfficiency[1],deltaY,branchingRatioC,branchingRatioBintoFinalDecay,fGlobalEfficiencyUncertainties[0],fGlobalEfficiencyUncertainties[1]);
534
535   //
536   // Finally: Correct from yields to cross-section
537   //
538   Int_t nbins = fhRECpt->GetNbinsX();
539   Double_t binwidth = fhRECpt->GetBinWidth(1);
540   Double_t *limits = new Double_t[nbins+1]; 
541   Double_t *binwidths = new Double_t[nbins]; 
542   Double_t xlow=0.;
543   for (Int_t i=1; i<=nbins; i++) {
544     binwidth = fhRECpt->GetBinWidth(i);
545     xlow = fhRECpt->GetBinLowEdge(i);
546     limits[i-1] = xlow;
547     binwidths[i-1] = binwidth;
548   }
549   limits[nbins] = xlow + binwidth;
550
551   
552   // declare the output histograms
553   fhSigmaCorr = new TH1D("fhSigmaCorr","corrected sigma",nbins,limits);
554   fhSigmaCorrMax = new TH1D("fhSigmaCorrMax","max corrected sigma",nbins,limits);
555   fhSigmaCorrMin = new TH1D("fhSigmaCorrMin","min corrected sigma",nbins,limits);
556   // and the output TGraphAsymmErrors
557   if (fAsymUncertainties){
558     fgSigmaCorr = new TGraphAsymmErrors(nbins+1);
559     fgSigmaCorrExtreme = new TGraphAsymmErrors(nbins+1);
560     fgSigmaCorrConservative = new TGraphAsymmErrors(nbins+1);
561   }
562
563   // protect against null denominator
564   if (deltaY==0. || fLuminosity[0]==0. || fTrigEfficiency[0]==0. || branchingRatioC==0.) {
565     AliError(" Hey you ! Why luminosity or trigger-efficiency or the c-BR or delta_y are set to zero ?! ");
566     delete [] limits;
567     delete [] binwidths;
568     return ;
569   }
570
571   Double_t value=0, errValue=0, errvalueMax=0., errvalueMin=0.;
572   Double_t errvalueExtremeMax=0., errvalueExtremeMin=0.;
573   Double_t errvalueConservativeMax=0., errvalueConservativeMin=0.;
574   for(Int_t ibin=1; ibin<=nbins; ibin++){
575
576     // Sigma calculation
577     //   Sigma = ( 1. / (lumi * delta_y * BR_c * eff_trig * eff_c ) ) * spectra (corrected for feed-down)
578     value = (fhDirectEffpt->GetBinContent(ibin) && fhDirectEffpt->GetBinContent(ibin)!=0.) ? 
579       ( fhYieldCorr->GetBinContent(ibin) / ( deltaY * branchingRatioC * fLuminosity[0] * fTrigEfficiency[0] * fhDirectEffpt->GetBinContent(ibin) ) )
580       : 0. ;
581
582     // Sigma statistical uncertainty:
583     //   delta_sigma = sigma * sqrt ( (delta_spectra/spectra)^2 )
584     errValue = (value!=0.) ?  value * (fhYieldCorr->GetBinError(ibin)/fhYieldCorr->GetBinContent(ibin))  : 0. ;
585
586
587     //
588     // Sigma systematic uncertainties
589     //
590     if (fAsymUncertainties) {
591       
592       //  (syst but feed-down) delta_sigma = sigma * sqrt ( (delta_spectra_syst/spectra)^2 +
593       //                                     (delta_lumi/lumi)^2 + (delta_eff_trig/eff_trig)^2 + (delta_eff/eff)^2  + (global_eff)^2 )
594       errvalueMax = value * TMath::Sqrt( (fgYieldCorr->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin))*(fgYieldCorr->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin)) + 
595                                          (fLuminosity[1]/fLuminosity[0])*(fLuminosity[1]/fLuminosity[0]) + 
596                                          (fTrigEfficiency[1]/fTrigEfficiency[0])*(fTrigEfficiency[1]/fTrigEfficiency[0])  +
597                                          (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin)) +
598                                          fGlobalEfficiencyUncertainties[0]*fGlobalEfficiencyUncertainties[0] );
599       errvalueMin = value * TMath::Sqrt( (fgYieldCorr->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin))*(fgYieldCorr->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin)) + 
600                                          (fLuminosity[1]/fLuminosity[0])*(fLuminosity[1]/fLuminosity[0]) + 
601                                          (fTrigEfficiency[1]/fTrigEfficiency[0])*(fTrigEfficiency[1]/fTrigEfficiency[0])  +
602                                          (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))  +
603                                          fGlobalEfficiencyUncertainties[0]*fGlobalEfficiencyUncertainties[0] );
604
605       // Uncertainties from feed-down
606       //      (feed-down syst) delta_sigma = sigma * sqrt ( (delta_spectra_fd/spectra_fd)^2 )
607       //   extreme case
608       errvalueExtremeMax = value * (fgYieldCorrExtreme->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin));
609       errvalueExtremeMin =  value * (fgYieldCorrExtreme->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin));
610       //
611       //   conservative case
612       errvalueConservativeMax = value * (fgYieldCorrConservative->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin));
613       errvalueConservativeMin =  value * (fgYieldCorrConservative->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin));
614
615     }
616     else {
617       // protect against null denominator
618       errvalueMax = (value!=0.) ?
619         value * TMath::Sqrt( (fLuminosity[1]/fLuminosity[0])*(fLuminosity[1]/fLuminosity[0]) + 
620                              (fTrigEfficiency[1]/fTrigEfficiency[0])*(fTrigEfficiency[1]/fTrigEfficiency[0])  +
621                              (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))  +
622                              fGlobalEfficiencyUncertainties[0]*fGlobalEfficiencyUncertainties[0] )
623         : 0. ;
624       errvalueMin = errvalueMax;
625     }
626     
627     // Fill the histograms
628     fhSigmaCorr->SetBinContent(ibin,value);
629     fhSigmaCorr->SetBinError(ibin,errValue);
630     // Fill the TGraphAsymmErrors
631     if (fAsymUncertainties) {
632       Double_t x = fhYieldCorr->GetBinCenter(ibin);
633       fgSigmaCorr->SetPoint(ibin,x,value); // i,x,y
634       fgSigmaCorr->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueMin,errvalueMax); // i,xl,xh,yl,yh
635       fhSigmaCorrMax->SetBinContent(ibin,value+errvalueMax);
636       fhSigmaCorrMin->SetBinContent(ibin,value-errvalueMin);
637       fgSigmaCorrExtreme->SetPoint(ibin,x,value); // i,x,y
638       fgSigmaCorrExtreme->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueExtremeMin,errvalueExtremeMax); // i,xl,xh,yl,yh
639       fgSigmaCorrConservative->SetPoint(ibin,x,value); // i,x,y
640       fgSigmaCorrConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueConservativeMin,errvalueConservativeMax); // i,xl,xh,yl,yh
641
642     }
643     
644   }
645   delete [] binwidths;
646   delete [] limits;
647
648 }
649
650 //_________________________________________________________________________________________________________
651 TH1D * AliHFPtSpectrum::EstimateEfficiencyRecoBin(TH1D *hSimu, TH1D *hReco, const char *name) {
652   //
653   // Function that computes the acceptance and efficiency correction
654   //  based on the simulated and reconstructed spectra
655   //  and using the reconstructed spectra bin width
656   //
657   //  eff = reco/sim ; err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
658   // 
659
660   if(!fhRECpt){
661     AliInfo("Hey, the reconstructed histogram was not set yet !"); 
662     return NULL;
663   }
664
665   Int_t nbins = fhRECpt->GetNbinsX();
666   Double_t *limits = new Double_t[nbins+1];
667   Double_t xlow=0.,binwidth=0.;
668   for (Int_t i=1; i<=nbins; i++) {
669     binwidth = fhRECpt->GetBinWidth(i);
670     xlow = fhRECpt->GetBinLowEdge(i);
671     limits[i-1] = xlow;
672   }
673   limits[nbins] = xlow + binwidth;
674
675   TH1D * hEfficiency = new TH1D(name," acceptance #times efficiency",nbins,limits);
676   
677   Double_t *sumSimu=new Double_t[nbins];
678   Double_t *sumReco=new Double_t[nbins];
679   for (Int_t ibin=0; ibin<nbins; ibin++){
680     sumSimu[ibin]=0.;  sumReco[ibin]=0.;
681   }
682   for (Int_t ibin=0; ibin<=hSimu->GetNbinsX(); ibin++){
683
684     for (Int_t ibinrec=0; ibinrec<nbins; ibinrec++){
685       if ( hSimu->GetBinCenter(ibin)>limits[ibinrec] && 
686            hSimu->GetBinCenter(ibin)<limits[ibinrec+1] ) {
687         sumSimu[ibinrec]+=hSimu->GetBinContent(ibin);
688       }
689       if ( hReco->GetBinCenter(ibin)>limits[ibinrec] && 
690            hReco->GetBinCenter(ibin)<limits[ibinrec+1] ) {
691         sumReco[ibinrec]+=hReco->GetBinContent(ibin);
692       }
693     }
694     
695   }
696
697
698   // the efficiency is computed as reco/sim (in each bin)
699   //  its uncertainty is err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
700   Double_t eff=0., erreff=0.;
701   for (Int_t ibinrec=0; ibinrec<nbins; ibinrec++) {
702     if (sumSimu[ibinrec]!= 0. && sumReco[ibinrec]!=0.) {
703       eff = sumReco[ibinrec] / sumSimu[ibinrec] ;
704       // protection in case eff > 1.0
705       // test calculation (make the argument of the sqrt positive)
706       erreff = TMath::Sqrt( eff * TMath::Abs(1.0 - eff) ) / TMath::Sqrt( sumSimu[ibinrec] );
707     }
708     else { eff=0.0; erreff=0.; }
709     hEfficiency->SetBinContent(ibinrec+1,eff);
710     hEfficiency->SetBinError(ibinrec+1,erreff);
711   }
712
713   delete [] sumSimu;
714   delete [] sumReco;
715
716   return (TH1D*)hEfficiency;
717 }
718
719 //_________________________________________________________________________________________________________
720 void AliHFPtSpectrum::EstimateAndSetDirectEfficiencyRecoBin(TH1D *hSimu, TH1D *hReco) {
721   //
722   // Function that computes the Direct  acceptance and efficiency correction
723   //  based on the simulated and reconstructed spectra
724   //  and using the reconstructed spectra bin width
725   //
726   //  eff = reco/sim ; err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
727   // 
728
729   if(!fhRECpt || !hSimu || !hReco){
730     AliError("Hey, the reconstructed histogram was not set yet !"); 
731     return;
732   }
733
734   fhDirectEffpt = EstimateEfficiencyRecoBin(hSimu,hReco,"fhDirectEffpt");
735   fhDirectEffpt->SetNameTitle("fhDirectEffpt"," direct acceptance #times efficiency");
736
737 }
738
739 //_________________________________________________________________________________________________________
740 void AliHFPtSpectrum::EstimateAndSetFeedDownEfficiencyRecoBin(TH1D *hSimu, TH1D *hReco) {
741   //
742   // Function that computes the Feed-Down acceptance and efficiency correction
743   //  based on the simulated and reconstructed spectra
744   //  and using the reconstructed spectra bin width
745   //
746   //  eff = reco/sim ; err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
747   // 
748   
749   if(!fhRECpt || !hSimu || !hReco){
750     AliError("Hey, the reconstructed histogram was not set yet !"); 
751     return;
752   }
753
754   fhFeedDownEffpt = EstimateEfficiencyRecoBin(hSimu,hReco,"fhFeedDownEffpt");
755   fhFeedDownEffpt->SetNameTitle("fhFeedDownEffpt"," feed-down acceptance #times efficiency");
756
757 }
758
759 //_________________________________________________________________________________________________________
760 Bool_t AliHFPtSpectrum::Initialize(){
761   //
762   // Initialization of the variables (histograms)
763   //
764
765   if (fFeedDownOption==0) { 
766     AliInfo("Getting ready for the corrections without feed-down consideration");
767   } else if (fFeedDownOption==1) { 
768     AliInfo("Getting ready for the fc feed-down correction calculation");
769   } else if (fFeedDownOption==2) {
770     AliInfo("Getting ready for the Nb feed-down correction calculation");
771   } else { AliError("The calculation option must be <=2"); return kFALSE; }
772
773   // Start checking the input histograms consistency
774   Bool_t areconsistent=kTRUE;
775
776   // General checks 
777   if (!fhDirectEffpt || !fhRECpt) {
778     AliError(" Reconstructed spectra and/or the Nc efficiency distributions are not defined");
779     return kFALSE;
780   }
781   areconsistent &= CheckHistosConsistency(fhRECpt,fhDirectEffpt);
782   if (!areconsistent) {
783     AliInfo("Histograms required for Nb correction are not consistent (bin width, bounds)"); 
784     return kFALSE;
785   }
786   if (fFeedDownOption==0) return kTRUE;
787
788   //
789   // Common checks for options 1 (fc) & 2(Nb)
790   if (!fhFeedDownMCpt || !fhFeedDownEffpt) {
791     AliError(" Theoretical Nb and/or the Nb efficiency distributions are not defined");
792     return kFALSE;
793   }
794   areconsistent &= CheckHistosConsistency(fhRECpt,fhFeedDownMCpt);
795   areconsistent &= CheckHistosConsistency(fhFeedDownMCpt,fhFeedDownEffpt);
796   if (fAsymUncertainties) {
797     if (!fhFeedDownMCptMax || !fhFeedDownMCptMin) {
798       AliError(" Max/Min theoretical Nb distributions are not defined");
799       return kFALSE;
800     }
801     areconsistent &= CheckHistosConsistency(fhFeedDownMCpt,fhFeedDownMCptMax);
802   }
803   if (!areconsistent) {
804     AliInfo("Histograms required for Nb correction are not consistent (bin width, bounds)"); 
805     return kFALSE;
806   }
807   if (fFeedDownOption>1) return kTRUE;  
808
809   //
810   // Now checks for option 1 (fc correction) 
811   if (!fhDirectMCpt) {
812     AliError("Theoretical Nc distributions is not defined");
813     return kFALSE;
814   }
815   areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhFeedDownMCpt);
816   areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhDirectEffpt);
817   if (fAsymUncertainties) {
818     if (!fhDirectMCptMax || !fhDirectMCptMin) {
819       AliError(" Max/Min theoretical Nc distributions are not defined");
820       return kFALSE;
821     }
822     areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhDirectMCptMax);
823   }
824   if (!areconsistent) {
825     AliInfo("Histograms required for fc correction are not consistent (bin width, bounds)"); 
826     return kFALSE;
827   }
828
829   return kTRUE;
830 }
831
832 //_________________________________________________________________________________________________________
833 Bool_t AliHFPtSpectrum::CheckHistosConsistency(TH1D *h1, TH1D *h2){
834   //
835   // Check the histograms consistency (bins, limits)
836   //
837
838   if (!h1 || !h2) {
839     AliError("One or both histograms don't exist");
840     return kFALSE;
841   }
842
843   Double_t binwidth1 = h1->GetBinWidth(1);
844   Double_t binwidth2 = h2->GetBinWidth(1);
845   Double_t min1 = h1->GetBinCenter(1) - (binwidth1/2.) ; 
846 //   Double_t max1 = h1->GetBinCenter(nbins1) + (binwidth1/2.) ;
847   Double_t min2 = h2->GetBinCenter(1) - (binwidth2/2.) ;
848 //   Double_t max2 = h2->GetBinCenter(nbins2) + (binwidth2/2.) ;
849
850   if (binwidth1!=binwidth2) {
851     AliInfo(" histograms with different bin width");
852     return kFALSE;
853   }
854   if (min1!=min2) {
855     AliInfo(" histograms with different minimum");
856     return kFALSE;
857   }
858 //   if (max1!=max2) {
859 //     AliInfo(" histograms with different maximum");
860 //     return kFALSE;
861 //   }
862
863   return kTRUE;
864 }
865
866 //_________________________________________________________________________________________________________
867 void AliHFPtSpectrum::CalculateFeedDownCorrectionFc(){ 
868   //
869   // Compute fc factor and its uncertainties bin by bin
870   //   fc = 1 / ( 1 + (eff_b/eff_c)*(N_b/N_c) ) 
871   //
872   // uncertainties: (conservative) combine the upper/lower N_b & N_c predictions together
873   //                (extreme) combine the upper N_b predictions with the lower N_c predictions & viceversa
874   //                systematic uncertainty on the acceptance x efficiency b/c ratio are included 
875   //
876   
877   // define the variables
878   Int_t nbins = fhRECpt->GetNbinsX();
879   Double_t binwidth = fhRECpt->GetBinWidth(1);
880   Double_t *limits = new Double_t[nbins+1];
881   Double_t *binwidths = new Double_t[nbins];
882   Double_t xlow=0.;
883   for (Int_t i=1; i<=nbins; i++) {
884     binwidth = fhRECpt->GetBinWidth(i);
885     xlow = fhRECpt->GetBinLowEdge(i);
886     limits[i-1] = xlow;
887     binwidths[i-1] = binwidth;
888   }
889   limits[nbins] = xlow + binwidth;
890
891   Double_t correction=1.;
892   Double_t theoryRatio=1.;
893   Double_t effRatio=1.; 
894   Double_t correctionExtremeA=1., correctionExtremeB=1.;
895   Double_t theoryRatioExtremeA=1., theoryRatioExtremeB=1.;
896   Double_t correctionConservativeA=1., correctionConservativeB=1.;
897   Double_t theoryRatioConservativeA=1., theoryRatioConservativeB=1.;
898   Double_t correctionUnc=0.;
899   Double_t correctionExtremeAUnc=0., correctionExtremeBUnc=0.;
900   Double_t correctionConservativeAUnc=0., correctionConservativeBUnc=0.;
901
902   // declare the output histograms
903   fhFc = new TH1D("fhFc","fc correction factor",nbins,limits);
904   fhFcMax = new TH1D("fhFcMax","max fc correction factor",nbins,limits);
905   fhFcMin = new TH1D("fhFcMin","min fc correction factor",nbins,limits);
906   // two local control histograms
907   TH1D *hTheoryRatio = new TH1D("hTheoryRatio","Theoretical B-->D over c-->D (feed-down/direct) ratio",nbins,limits);
908   TH1D *hEffRatio = new TH1D("hEffRatio","Efficiency B-->D over c-->D (feed-down/direct) ratio",nbins,limits);
909   // and the output TGraphAsymmErrors
910   if (fAsymUncertainties) {
911     fgFcExtreme = new TGraphAsymmErrors(nbins+1);
912     fgFcExtreme->SetNameTitle("fgFcExtreme","fgFcExtreme");
913     fgFcConservative = new TGraphAsymmErrors(nbins+1);
914     fgFcConservative->SetNameTitle("fgFcConservative","fgFcConservative");
915   }
916
917
918   //
919   // Compute fc
920   //
921   for (Int_t ibin=1; ibin<=nbins; ibin++) {
922
923     //  theory_ratio = (N_b/N_c) 
924     theoryRatio = (fhDirectMCpt->GetBinContent(ibin)>0. && fhFeedDownMCpt->GetBinContent(ibin)>0.) ? 
925       fhFeedDownMCpt->GetBinContent(ibin) / fhDirectMCpt->GetBinContent(ibin) : 1.0 ;
926
927     //
928     // Calculate the uncertainty [ considering only the theoretical uncertainties on Nb & Nc for now !!! ]
929     //
930     // extreme A = direct-max, feed-down-min
931     theoryRatioExtremeA = (fhDirectMCptMax->GetBinContent(ibin)>0. && fhFeedDownMCptMin->GetBinContent(ibin)>0.) ? 
932       fhFeedDownMCptMin->GetBinContent(ibin) / fhDirectMCptMax->GetBinContent(ibin) : 1.0 ;
933     // extreme B = direct-min, feed-down-max
934     theoryRatioExtremeB = (fhDirectMCptMin->GetBinContent(ibin)>0. && fhDirectMCptMax->GetBinContent(ibin)>0.) ? 
935       fhFeedDownMCptMax->GetBinContent(ibin) / fhDirectMCptMin->GetBinContent(ibin) : 1.0 ;
936     // conservative A = direct-max, feed-down-max
937     theoryRatioConservativeA = (fhDirectMCptMax->GetBinContent(ibin)>0. && fhFeedDownMCptMin->GetBinContent(ibin)>0.) ? 
938       fhFeedDownMCptMax->GetBinContent(ibin) / fhDirectMCptMax->GetBinContent(ibin) : 1.0 ;
939     // conservative B = direct-min, feed-down-min
940     theoryRatioConservativeB = (fhDirectMCptMin->GetBinContent(ibin)>0. && fhDirectMCptMax->GetBinContent(ibin)>0.) ? 
941       fhFeedDownMCptMin->GetBinContent(ibin) / fhDirectMCptMin->GetBinContent(ibin) : 1.0 ;
942
943     //  eff_ratio = (eff_b/eff_c)
944     effRatio = (fhDirectEffpt->GetBinContent(ibin) && fhDirectEffpt->GetBinContent(ibin)!=0.) ? 
945       fhFeedDownEffpt->GetBinContent(ibin) / fhDirectEffpt->GetBinContent(ibin) : 1.0 ;
946
947     //   fc = 1 / ( 1 + (eff_b/eff_c)*(N_b/N_c) ) 
948     if( TMath::Abs(effRatio - 1.0)<0.01 || TMath::Abs(theoryRatio - 1.0)<0.01 ) {
949       correction = 1.0;
950       correctionExtremeA = 1.0;
951       correctionExtremeB = 1.0;
952       correctionConservativeA = 1.0; 
953       correctionConservativeB = 1.0;
954     }
955     else {
956       correction = ( 1. / ( 1 + ( effRatio * theoryRatio ) ) );
957       correctionExtremeA = ( 1. / ( 1 + ( effRatio * theoryRatioExtremeA ) ) );
958       correctionExtremeB = ( 1. / ( 1 + ( effRatio * theoryRatioExtremeB ) ) );
959       correctionConservativeA = ( 1. / ( 1 + ( effRatio * theoryRatioConservativeA ) ) );
960       correctionConservativeB = ( 1. / ( 1 + ( effRatio * theoryRatioConservativeB ) ) );
961     }
962
963
964     // fc uncertainty from (eff_b/eff_c) = fc^2 * (N_b/N_c) * delta(eff_b/eff_c)
965     //  delta(eff_b/eff_c) is a percentage = fGlobalEfficiencyUncertainties[1]*effRatio
966     //    correctionUnc = correction*correction * theoryRatio * fGlobalEfficiencyUncertainties[1]*effRatio;
967     //     correctionExtremeAUnc = correctionExtremeA*correctionExtremeA * theoryRatioExtremeA  * fGlobalEfficiencyUncertainties[1]*effRatio;
968     //     correctionExtremeBUnc = correctionExtremeB*correctionExtremeB * theoryRatioExtremeB  * fGlobalEfficiencyUncertainties[1]*effRatio;
969     //     correctionConservativeAUnc = correctionConservativeA*correctionConservativeA * theoryRatioConservativeA  * fGlobalEfficiencyUncertainties[1]*effRatio;
970     //     correctionConservativeBUnc = correctionConservativeB*correctionConservativeB * theoryRatioConservativeB  * fGlobalEfficiencyUncertainties[1]*effRatio;
971     correctionUnc = correction*correction * theoryRatio * effRatio *
972       TMath::Sqrt( fGlobalEfficiencyUncertainties[1]*fGlobalEfficiencyUncertainties[1] + 
973                    (fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) +
974                    (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin)) 
975                    );
976     correctionExtremeAUnc = correctionExtremeA*correctionExtremeA * theoryRatioExtremeA  * effRatio *
977       TMath::Sqrt( fGlobalEfficiencyUncertainties[1]*fGlobalEfficiencyUncertainties[1] + 
978                    (fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) +
979                    (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin)) 
980                    );
981     correctionExtremeBUnc = correctionExtremeB*correctionExtremeB * theoryRatioExtremeB  * effRatio *
982       TMath::Sqrt( fGlobalEfficiencyUncertainties[1]*fGlobalEfficiencyUncertainties[1] + 
983                    (fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) +
984                    (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin)) 
985                    );
986     correctionConservativeAUnc = correctionConservativeA*correctionConservativeA * theoryRatioConservativeA  *effRatio *
987       TMath::Sqrt( fGlobalEfficiencyUncertainties[1]*fGlobalEfficiencyUncertainties[1] + 
988                    (fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) +
989                    (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin)) 
990                    );
991     correctionConservativeBUnc = correctionConservativeB*correctionConservativeB * theoryRatioConservativeB  *effRatio *
992       TMath::Sqrt( fGlobalEfficiencyUncertainties[1]*fGlobalEfficiencyUncertainties[1] + 
993                    (fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) +
994                    (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin)) 
995                    );
996
997
998     // Fill in the histograms
999     hTheoryRatio->SetBinContent(ibin,theoryRatio);
1000     hEffRatio->SetBinContent(ibin,effRatio);
1001     fhFc->SetBinContent(ibin,correction);
1002     if (fAsymUncertainties) {
1003       Double_t x = fhDirectMCpt->GetBinCenter(ibin);
1004       Double_t val[4] = { correctionExtremeA + correctionExtremeAUnc, correctionExtremeA - correctionExtremeAUnc, 
1005                           correctionExtremeB + correctionExtremeBUnc, correctionExtremeB - correctionExtremeBUnc };
1006       Double_t uncExtremeMin = correction - TMath::MinElement(4,val);
1007       Double_t uncExtremeMax = TMath::MaxElement(4,val) - correction;
1008       fgFcExtreme->SetPoint(ibin,x,correction); // i,x,y
1009       fgFcExtreme->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),uncExtremeMin,uncExtremeMax); // i,xl,xh,yl,yh
1010       fhFcMax->SetBinContent(ibin,correction+uncExtremeMax);
1011       fhFcMin->SetBinContent(ibin,correction-uncExtremeMin);
1012       Double_t consval[4] = { correctionConservativeA - correctionConservativeAUnc, correctionConservativeA + correctionConservativeAUnc,
1013                               correctionConservativeB - correctionConservativeBUnc, correctionConservativeB + correctionConservativeBUnc};
1014       Double_t uncConservativeMin = correction - TMath::MinElement(4,consval);
1015       Double_t uncConservativeMax =  TMath::MaxElement(4,consval) - correction;
1016       fgFcConservative->SetPoint(ibin,x,correction); // i,x,y
1017       fgFcConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),uncConservativeMin,uncConservativeMax); // i,xl,xh,yl,yh
1018       if( !(correction>0.) ){
1019         fgFcExtreme->SetPoint(ibin,x,0.); // i,x,y
1020         fgFcExtreme->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),0.,0.); // i,xl,xh,yl,yh
1021         fgFcConservative->SetPoint(ibin,x,0.); // i,x,y
1022         fgFcConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),0.,0.); // i,xl,xh,yl,yh
1023       }
1024     }
1025
1026   }
1027   delete [] binwidths;
1028   delete [] limits;
1029
1030 }
1031
1032 //_________________________________________________________________________________________________________
1033 void AliHFPtSpectrum::CalculateFeedDownCorrectedSpectrumFc(){
1034   //
1035   // Compute the feed-down corrected spectrum if feed-down correction is done via fc factor (bin by bin)
1036   //    physics = reco * fc / bin-width
1037   //
1038   //    uncertainty:             (stat) delta_physics = physics * sqrt ( (delta_reco/reco)^2 )
1039   //               (syst but feed-down) delta_physics = physics * sqrt ( (delta_reco_syst/reco)^2 )
1040   //                   (feed-down syst) delta_physics = physics * sqrt ( (delta_fc/fc)^2 )
1041   //
1042   //    ( Calculation done bin by bin )
1043
1044   if (!fhFc || !fhRECpt) {
1045     AliError(" Reconstructed or fc distributions are not defined");
1046     return;
1047   }
1048
1049   Int_t nbins = fhRECpt->GetNbinsX();
1050   Double_t value = 0., errvalue = 0., errvalueMax= 0., errvalueMin= 0.;
1051   Double_t valueExtremeMax= 0., valueExtremeMin= 0.;
1052   Double_t valueConservativeMax= 0., valueConservativeMin= 0.;
1053   Double_t binwidth = fhRECpt->GetBinWidth(1);
1054   Double_t *limits = new Double_t[nbins+1];
1055   Double_t *binwidths = new Double_t[nbins];
1056   Double_t xlow=0.;
1057   for (Int_t i=1; i<=nbins; i++) {
1058     binwidth = fhRECpt->GetBinWidth(i);
1059     xlow = fhRECpt->GetBinLowEdge(i);
1060     limits[i-1] = xlow;
1061     binwidths[i-1] = binwidth;
1062   }
1063   limits[nbins] = xlow + binwidth;
1064   
1065   // declare the output histograms
1066   fhYieldCorr = new TH1D("fhYieldCorr","corrected yield (by fc)",nbins,limits);
1067   fhYieldCorrMax = new TH1D("fhYieldCorrMax","max corrected yield (by fc)",nbins,limits);
1068   fhYieldCorrMin = new TH1D("fhYieldCorrMin","min corrected yield (by fc)",nbins,limits);
1069   // and the output TGraphAsymmErrors
1070   if (fAsymUncertainties){
1071     fgYieldCorr = new TGraphAsymmErrors(nbins+1);
1072     fgYieldCorrExtreme = new TGraphAsymmErrors(nbins+1);
1073     fgYieldCorrConservative = new TGraphAsymmErrors(nbins+1);
1074   }
1075   
1076   //
1077   // Do the calculation
1078   // 
1079   for (Int_t ibin=1; ibin<=nbins; ibin++) {
1080
1081     // calculate the value 
1082     //    physics = reco * fc / bin-width
1083     value = (fhRECpt->GetBinContent(ibin) && fhFc->GetBinContent(ibin)) ? 
1084       fhRECpt->GetBinContent(ibin) * fhFc->GetBinContent(ibin) : 0. ;
1085     value /= fhRECpt->GetBinWidth(ibin) ;
1086     
1087     // Statistical uncertainty 
1088     //    (stat) delta_physics = physics * sqrt ( (delta_reco/reco)^2 )
1089     errvalue = (value!=0. && fhRECpt->GetBinContent(ibin) && fhRECpt->GetBinContent(ibin)!=0.) ?
1090       value * (fhRECpt->GetBinError(ibin)/fhRECpt->GetBinContent(ibin)) : 0. ; 
1091
1092     // Calculate the systematic uncertainties
1093     //    (syst but feed-down) delta_physics = physics * sqrt ( (delta_reco_syst/reco)^2 )
1094     //        (feed-down syst) delta_physics = physics * sqrt ( (delta_fc/fc)^2 )
1095     //
1096     // Protect against null denominator. If so, define uncertainty as null
1097     if (fhRECpt->GetBinContent(ibin) && fhRECpt->GetBinContent(ibin)!=0.) {
1098
1099       if (fAsymUncertainties) {
1100
1101         // Systematics but feed-down
1102         if (fgRECSystematics) { 
1103           errvalueMax = value * ( fgRECSystematics->GetErrorYhigh(ibin) / fhRECpt->GetBinContent(ibin) );
1104           errvalueMin = value * ( fgRECSystematics->GetErrorYlow(ibin) / fhRECpt->GetBinContent(ibin) );
1105         }
1106         else { errvalueMax = 0.; errvalueMin = 0.; }
1107         
1108         // Extreme feed-down systematics
1109         valueExtremeMax = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) + fgFcExtreme->GetErrorYhigh(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
1110         valueExtremeMin = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) - fgFcExtreme->GetErrorYlow(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
1111         
1112         // Conservative feed-down systematics
1113         valueConservativeMax = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) + fgFcConservative->GetErrorYhigh(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
1114         valueConservativeMin = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) - fgFcConservative->GetErrorYlow(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
1115
1116       }
1117
1118     }
1119     else { errvalueMax = 0.; errvalueMin = 0.; }
1120     
1121     // fill in the histograms
1122     fhYieldCorr->SetBinContent(ibin,value);
1123     fhYieldCorr->SetBinError(ibin,errvalue);
1124     if (fAsymUncertainties) {
1125       Double_t center = fhYieldCorr->GetBinCenter(ibin);
1126       fgYieldCorr->SetPoint(ibin,center,value); // i,x,y
1127       fgYieldCorr->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueMin,errvalueMax); // i,xl,xh,yl,yh
1128       fhYieldCorrMax->SetBinContent(ibin,value+errvalueMax); 
1129       fhYieldCorrMin->SetBinContent(ibin,value-errvalueMin);
1130       fgYieldCorrExtreme->SetPoint(ibin,center,value);
1131       fgYieldCorrExtreme->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),value-valueExtremeMin,valueExtremeMax-value);
1132       fgYieldCorrConservative->SetPoint(ibin,center,value);
1133       fgYieldCorrConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),value-valueConservativeMin,valueConservativeMax-value);
1134     }
1135
1136   }
1137   delete [] binwidths;
1138   delete [] limits;
1139   
1140 }
1141
1142 //_________________________________________________________________________________________________________
1143 void AliHFPtSpectrum::CalculateFeedDownCorrectedSpectrumNb(Double_t deltaY, Double_t branchingRatioBintoFinalDecay) {
1144   //
1145   // Compute the feed-down corrected spectrum if feed-down correction is done via Nb (bin by bin)
1146   //    physics =  [ reco  - (lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th) ] / bin-width
1147   //
1148   //    uncertainty:   (stat)  delta_physics = sqrt ( (delta_reco)^2 )  / bin-width
1149   //     (syst but feed-down)  delta_physics = sqrt ( (delta_reco_syst)^2 )  / bin-width
1150   //         (feed-down syst)  delta_physics = sqrt ( (k*delta_lumi/lumi)^2 + (k*delta_eff_trig/eff_trig)^2 
1151   //                                                   + (k*delta_Nb/Nb)^2 + (k*delta_eff/eff)^2  + (k*global_eff_ratio)^2 ) / bin-width
1152   //                    where k = lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th
1153   //
1154
1155   Int_t nbins = fhRECpt->GetNbinsX();
1156   Double_t binwidth = fhRECpt->GetBinWidth(1);
1157   Double_t value = 0., errvalue = 0., errvalueMax = 0., errvalueMin = 0., kfactor = 0.;
1158   Double_t errvalueExtremeMax = 0., errvalueExtremeMin = 0.;
1159   Double_t *limits = new Double_t[nbins+1];
1160   Double_t *binwidths = new Double_t[nbins];
1161   Double_t xlow=0.;
1162   for (Int_t i=1; i<=nbins; i++) {
1163     binwidth = fhRECpt->GetBinWidth(i);
1164     xlow = fhRECpt->GetBinLowEdge(i);
1165     limits[i-1] = xlow;
1166     binwidths[i-1] = binwidth;
1167   }
1168   limits[nbins] = xlow + binwidth;
1169   
1170   // declare the output histograms
1171   fhYieldCorr = new TH1D("fhYieldCorr","corrected yield (by Nb)",nbins,limits);
1172   fhYieldCorrMax = new TH1D("fhYieldCorrMax","max corrected yield (by Nb)",nbins,limits);
1173   fhYieldCorrMin = new TH1D("fhYieldCorrMin","min corrected yield (by Nb)",nbins,limits);
1174   // and the output TGraphAsymmErrors
1175   if (fAsymUncertainties){
1176     fgYieldCorr = new TGraphAsymmErrors(nbins+1);
1177     fgYieldCorrExtreme = new TGraphAsymmErrors(nbins+1);
1178     fgYieldCorrConservative = new TGraphAsymmErrors(nbins+1);
1179     // Define fc-conservative 
1180     fgFcConservative = new TGraphAsymmErrors(nbins+1);
1181     AliInfo(" Beware the conservative & extreme uncertainties are equal by definition !");
1182   }
1183
1184   // variables to define fc-conservative 
1185   double correction=0, correctionMax=0., correctionMin=0.;
1186
1187   //
1188   // Do the calculation
1189   // 
1190   for (Int_t ibin=1; ibin<=nbins; ibin++) {
1191     
1192     // Calculate the value
1193     //    physics =  [ reco  - (lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th) ] / bin-width
1194     value = ( fhRECpt->GetBinContent(ibin) && fhRECpt->GetBinContent(ibin)!=0. && 
1195               fhFeedDownMCpt->GetBinContent(ibin)>0. && fhFeedDownEffpt->GetBinContent(ibin)>0. ) ?
1196       fhRECpt->GetBinContent(ibin) - (deltaY*branchingRatioBintoFinalDecay*fLuminosity[0]*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) ) 
1197       : 0. ;
1198     value /= fhRECpt->GetBinWidth(ibin);
1199
1200     //  Statistical uncertainty:   delta_physics = sqrt ( (delta_reco)^2 )  / bin-width
1201     errvalue = (value!=0. && fhRECpt->GetBinError(ibin) && fhRECpt->GetBinError(ibin)!=0.)  ? 
1202       fhRECpt->GetBinError(ibin) : 0.;
1203     errvalue /= fhRECpt->GetBinWidth(ibin);
1204
1205     // Correction (fc) : Estimate of the relative amount feed-down subtracted
1206     // correction =  [ 1  - (lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th)/reco ] 
1207     correction = 1 - (deltaY*branchingRatioBintoFinalDecay*fLuminosity[0]*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) ) / fhRECpt->GetBinContent(ibin) ;
1208
1209     // Systematic uncertainties
1210     //     (syst but feed-down)  delta_physics = sqrt ( (delta_reco_syst)^2 )  / bin-width
1211     //         (feed-down syst)  delta_physics = sqrt ( (k*delta_lumi/lumi)^2 + (k*delta_eff_trig/eff_trig)^2 
1212     //                                                   + (k*delta_Nb/Nb)^2 + (k*delta_eff/eff)^2 + (k*global_eff_ratio)^2 ) / bin-width
1213     //                    where k = lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th
1214     kfactor = deltaY*branchingRatioBintoFinalDecay*fLuminosity[0]*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) ;
1215
1216     //
1217     if (fAsymUncertainties) {
1218       Double_t nb =  fhFeedDownMCpt->GetBinContent(ibin);
1219       Double_t nbDmax = fhFeedDownMCptMax->GetBinContent(ibin) - fhFeedDownMCpt->GetBinContent(ibin);
1220       Double_t nbDmin = fhFeedDownMCpt->GetBinContent(ibin) - fhFeedDownMCptMin->GetBinContent(ibin);
1221
1222       // Systematics but feed-down
1223       if (fgRECSystematics){
1224         errvalueMax = fgRECSystematics->GetErrorYhigh(ibin) / fhRECpt->GetBinWidth(ibin) ;
1225         errvalueMin = fgRECSystematics->GetErrorYlow(ibin) / fhRECpt->GetBinWidth(ibin);
1226       }
1227       else { errvalueMax = 0.; errvalueMin = 0.; }
1228   
1229       // Feed-down systematics
1230       // min value with the maximum Nb
1231       errvalueExtremeMin = TMath::Sqrt( ( (kfactor*fLuminosity[1]/fLuminosity[0])*(kfactor*fLuminosity[1]/fLuminosity[0]) ) +
1232                                         ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) +
1233                                         ( (kfactor*nbDmax/nb)*(kfactor*nbDmax/nb) )  +
1234                                         ( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) ) +
1235                                         ( (kfactor*fGlobalEfficiencyUncertainties[1])*(kfactor*fGlobalEfficiencyUncertainties[1]) ) 
1236                                         ) / fhRECpt->GetBinWidth(ibin);
1237       // max value with the minimum Nb
1238       errvalueExtremeMax =  TMath::Sqrt( ( (kfactor*fLuminosity[1]/fLuminosity[0])*(kfactor*fLuminosity[1]/fLuminosity[0]) ) +
1239                                          ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) +
1240                                          ( (kfactor*nbDmin/nb)*(kfactor*nbDmin/nb) )  +
1241                                          ( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))  ) +
1242                                          ( (kfactor*fGlobalEfficiencyUncertainties[1])*(kfactor*fGlobalEfficiencyUncertainties[1]) )
1243                                          ) / fhRECpt->GetBinWidth(ibin);
1244
1245       // Correction systematics (fc)
1246       // min value with the maximum Nb
1247       correctionMin = TMath::Sqrt( ( (kfactor*fLuminosity[1]/fLuminosity[0])*(kfactor*fLuminosity[1]/fLuminosity[0]) ) +
1248                                    ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) +
1249                                    ( (kfactor*nbDmax/nb)*(kfactor*nbDmax/nb) )  +
1250                                    ( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) ) +
1251                                    ( (kfactor*fGlobalEfficiencyUncertainties[1])*(kfactor*fGlobalEfficiencyUncertainties[1]) ) 
1252                                    ) / fhRECpt->GetBinContent(ibin) ;
1253       // max value with the minimum Nb
1254       correctionMax =  TMath::Sqrt( ( (kfactor*fLuminosity[1]/fLuminosity[0])*(kfactor*fLuminosity[1]/fLuminosity[0]) ) +
1255                                     ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) +
1256                                     ( (kfactor*nbDmin/nb)*(kfactor*nbDmin/nb) )  +
1257                                     ( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))       ) +
1258                                     ( (kfactor*fGlobalEfficiencyUncertainties[1])*(kfactor*fGlobalEfficiencyUncertainties[1]) )
1259                                     ) / fhRECpt->GetBinContent(ibin) ;
1260     }
1261     else{ // Don't consider Nb uncertainty in this case [ to be tested!!! ]
1262       errvalueExtremeMax =  TMath::Sqrt( ( (kfactor*fLuminosity[1]/fLuminosity[0])*(kfactor*fLuminosity[1]/fLuminosity[0]) ) +
1263                                          ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) )  +
1264                                          ( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))  )  +
1265                                          ( (kfactor*fGlobalEfficiencyUncertainties[1])*(kfactor*fGlobalEfficiencyUncertainties[1]) )
1266                                          ) / fhRECpt->GetBinWidth(ibin);
1267       errvalueExtremeMin =  errvalueExtremeMax ;
1268     }
1269     
1270
1271     // fill in histograms
1272     fhYieldCorr->SetBinContent(ibin,value);
1273     fhYieldCorr->SetBinError(ibin,errvalue);
1274     if (fAsymUncertainties) {
1275       Double_t x = fhYieldCorr->GetBinCenter(ibin);
1276       fgYieldCorr->SetPoint(ibin,x,value); // i,x,y
1277       fgYieldCorr->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueMin,errvalueMax); // i,xl,xh,yl,yh
1278       fhYieldCorrMax->SetBinContent(ibin,value+errvalueMax); 
1279       fhYieldCorrMin->SetBinContent(ibin,value-errvalueMin);
1280       fgYieldCorrExtreme->SetPoint(ibin,x,value); // i,x,y
1281       fgYieldCorrExtreme->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueExtremeMin,errvalueExtremeMax); // i,xl,xh,yl,yh
1282       fgYieldCorrConservative->SetPoint(ibin,x,value); // i,x,y
1283       fgYieldCorrConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueExtremeMin,errvalueExtremeMax); // i,xl,xh,yl,yh
1284       //      cout << " bin " << ibin << ", correction " << correction << ", min correction unc " << correctionMin << ", max correction unc " << correctionMax << endl;
1285       if(correction>0.){
1286         fgFcConservative->SetPoint(ibin,x,correction);
1287         fgFcConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),correctionMin,correctionMax);
1288       }
1289       else{
1290         fgFcConservative->SetPoint(ibin,x,0.);
1291         fgFcConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),0.,0.);
1292       }
1293     }
1294
1295   }
1296   delete [] binwidths;
1297   delete [] limits;
1298   
1299 }
1300
1301
1302 //_________________________________________________________________________________________________________
1303 void AliHFPtSpectrum::ComputeSystUncertainties(AliHFSystErr *systematics, Bool_t combineFeedDown) {
1304   //
1305   // Function that re-calculates the global systematic uncertainties
1306   //   by calling the class AliHFSystErr and combining those
1307   //   (in quadrature) with the feed-down subtraction uncertainties
1308   //
1309
1310
1311   // Estimate the feed-down uncertainty in percentage
1312   Int_t nentries = fgSigmaCorrConservative->GetN();
1313   TGraphAsymmErrors *grErrFeeddown = new TGraphAsymmErrors(nentries);
1314   Double_t x=0., y=0., errx=0., erryl=0., erryh=0;
1315   for(Int_t i=0; i<nentries; i++) {
1316     fgSigmaCorrConservative->GetPoint(i,x,y);
1317     errx = fgSigmaCorrConservative->GetErrorXlow(i) ;
1318     erryl = fgSigmaCorrConservative->GetErrorYlow(i) / y ;
1319     erryh = fgSigmaCorrConservative->GetErrorYhigh(i) / y ;
1320     grErrFeeddown->SetPoint(i,x,0.);
1321     grErrFeeddown->SetPointError(i,errx,errx,erryl,erryh); //i, xl, xh, yl, yh
1322   }
1323
1324   // Draw all the systematics independently
1325   systematics->DrawErrors(grErrFeeddown);
1326
1327   // Set the sigma systematic uncertainties
1328   // possibly combine with the feed-down uncertainties 
1329   Double_t errylcomb=0., erryhcomb=0;
1330   for(Int_t i=1; i<nentries; i++) {
1331     fgSigmaCorr->GetPoint(i,x,y);
1332     errx = grErrFeeddown->GetErrorXlow(i) ;
1333     erryl = grErrFeeddown->GetErrorYlow(i);
1334     erryh = grErrFeeddown->GetErrorYhigh(i);
1335     if (combineFeedDown) {
1336       errylcomb = systematics->GetTotalSystErr(x,erryl) * y ;
1337       erryhcomb = systematics->GetTotalSystErr(x,erryh) * y ;
1338     } else {
1339       errylcomb = systematics->GetTotalSystErr(x) * y ;
1340       erryhcomb = systematics->GetTotalSystErr(x) * y ;
1341     }
1342     fgSigmaCorr->SetPointError(i,errx,errx,errylcomb,erryhcomb);
1343   }
1344
1345 }
1346
1347
1348 //_________________________________________________________________________________________________________
1349 void AliHFPtSpectrum::DrawSpectrum(TGraphAsymmErrors *gPrediction) {
1350   //
1351   // Example method to draw the corrected spectrum & the theoretical prediction
1352   //
1353
1354   TCanvas *csigma = new TCanvas("csigma","Draw the corrected cross-section & the prediction");
1355   csigma->SetFillColor(0);
1356   gPrediction->GetXaxis()->SetTitleSize(0.05);
1357   gPrediction->GetXaxis()->SetTitleOffset(0.95);
1358   gPrediction->GetYaxis()->SetTitleSize(0.05);
1359   gPrediction->GetYaxis()->SetTitleOffset(0.95);
1360   gPrediction->GetXaxis()->SetTitle("p_{T}  [GeV]");
1361   gPrediction->GetYaxis()->SetTitle("BR #times #frac{d#sigma}{dp_{T}} |_{|y|<0.5}   [pb/GeV]");
1362   gPrediction->SetLineColor(kGreen+2);
1363   gPrediction->SetLineWidth(3);
1364   gPrediction->SetFillColor(kGreen+1);
1365   gPrediction->Draw("3CA");
1366   fgSigmaCorr->SetLineColor(kRed);
1367   fgSigmaCorr->SetLineWidth(1);
1368   fgSigmaCorr->SetFillColor(kRed);
1369   fgSigmaCorr->SetFillStyle(0);
1370   fgSigmaCorr->Draw("2");
1371   fhSigmaCorr->SetMarkerColor(kRed);
1372   fhSigmaCorr->Draw("esame");
1373   csigma->SetLogy();
1374   TLegend * leg = new TLegend(0.7,0.75,0.87,0.5);
1375   leg->SetBorderSize(0);
1376   leg->SetLineColor(0);
1377   leg->SetFillColor(0);
1378   leg->SetTextFont(42);
1379   leg->AddEntry(gPrediction,"FONLL ","fl");
1380   leg->AddEntry(fhSigmaCorr,"data stat. unc.","pl");
1381   leg->AddEntry(fgSigmaCorr,"data syst. unc.","f");
1382   leg->Draw();
1383   csigma->Draw();
1384
1385 }
1386
1387 //_________________________________________________________________________________________________________
1388 TH1D * AliHFPtSpectrum::ReweightHisto(TH1D *hToReweight, TH1D *hReference){
1389   //
1390   // Function to  reweight histograms for testing purposes: 
1391   // This function takes the histo hToReweight and reweights 
1392   //  it (its pt shape) with respect to hReference 
1393   // 
1394
1395   // check histograms consistency
1396   Bool_t areconsistent=kTRUE;
1397   areconsistent &= CheckHistosConsistency(hToReweight,hReference);
1398   if (!areconsistent) {
1399     AliInfo("the histograms to reweight are not consistent (bin width, bounds)"); 
1400     return NULL;
1401   }
1402
1403   // define a new empty histogram
1404   TH1D *hReweighted = (TH1D*)hToReweight->Clone("hReweighted");
1405   hReweighted->Reset();
1406   Double_t weight=1.0;
1407   Double_t yvalue=1.0; 
1408   Double_t integralRef = hReference->Integral();
1409   Double_t integralH = hToReweight->Integral();
1410
1411   // now reweight the spectra
1412   //
1413   // the weight is the relative probability of the given pt bin in the reference histo
1414   //  divided by its relative probability (to normalize it) on the histo to re-weight
1415   for (Int_t i=0; i<=hToReweight->GetNbinsX(); i++) {
1416     weight = (hReference->GetBinContent(i)/integralRef) / (hToReweight->GetBinContent(i)/integralH) ;
1417     yvalue = hToReweight->GetBinContent(i);
1418     hReweighted->SetBinContent(i,yvalue*weight);
1419   }
1420
1421   return (TH1D*)hReweighted;
1422 }
1423
1424 //_________________________________________________________________________________________________________
1425 TH1D * AliHFPtSpectrum::ReweightRecHisto(TH1D *hRecToReweight, TH1D *hMCToReweight, TH1D *hMCReference){
1426   //
1427   // Function to  reweight histograms for testing purposes: 
1428   // This function takes the histo hToReweight and reweights 
1429   //  it (its pt shape) with respect to hReference /hMCToReweight
1430   // 
1431
1432   // check histograms consistency
1433   Bool_t areconsistent=kTRUE;
1434   areconsistent &= CheckHistosConsistency(hMCToReweight,hMCReference);
1435   areconsistent &= CheckHistosConsistency(hRecToReweight,hMCReference);
1436   if (!areconsistent) {
1437     AliInfo("the histograms to reweight are not consistent (bin width, bounds)"); 
1438     return NULL;
1439   }
1440
1441   // define a new empty histogram
1442   TH1D *hReweighted = (TH1D*)hMCToReweight->Clone("hReweighted");
1443   hReweighted->Reset();
1444   TH1D *hRecReweighted = (TH1D*)hRecToReweight->Clone("hRecReweighted");
1445   hRecReweighted->Reset();
1446   Double_t weight=1.0;
1447   Double_t yvalue=1.0, yrecvalue=1.0; 
1448   Double_t integralRef = hMCReference->Integral();
1449   Double_t integralH = hMCToReweight->Integral();
1450
1451   // now reweight the spectra
1452   //
1453   // the weight is the relative probability of the given pt bin 
1454   //  that should be applied in the MC histo to get the reference histo shape
1455   //  Probabilities are properly normalized.
1456   for (Int_t i=0; i<=hMCToReweight->GetNbinsX(); i++) {
1457     weight = (hMCReference->GetBinContent(i)/integralRef) / (hMCToReweight->GetBinContent(i)/integralH) ;
1458     yvalue = hMCToReweight->GetBinContent(i);
1459     hReweighted->SetBinContent(i,yvalue*weight);
1460     yrecvalue = hRecToReweight->GetBinContent(i);
1461     hRecReweighted->SetBinContent(i,yrecvalue*weight);
1462   }
1463
1464   return (TH1D*)hRecReweighted;
1465 }
1466