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