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