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     // Variables initialization
 640     value=0.; errValue=0.; errvalueMax=0.; errvalueMin=0.;
 641     errvalueExtremeMax=0.; errvalueExtremeMin=0.;
 642     errvalueConservativeMax=0.; errvalueConservativeMin=0.;
 643     errvalueStatUncEffc=0.; errvalueStatUncEffb=0.;
 644
 645     // Sigma calculation
 646     //   Sigma = ( 1. / (lumi * delta_y * BR_c * eff_trig * eff_c ) ) * spectra (corrected for feed-down)
 647     value = (fhDirectEffpt->GetBinContent(ibin) && fhDirectEffpt->GetBinContent(ibin)!=0. && fhRECpt->GetBinContent(ibin)>0.) ?
 648       ( fhYieldCorr->GetBinContent(ibin) / ( deltaY * branchingRatioC * fLuminosity[0] * fTrigEfficiency[0] * fhDirectEffpt->GetBinContent(ibin) ) )
 649       : 0. ;
 650
 651     // Sigma statistical uncertainty:
 652     //   delta_sigma = sigma * sqrt ( (delta_spectra/spectra)^2 )
 653     errValue = (value!=0.) ?  value * (fhYieldCorr->GetBinError(ibin)/fhYieldCorr->GetBinContent(ibin))  : 0. ;
 654
 655     //    cout<< " x "<< fhRECpt->GetBinCenter(ibin) << " sigma " << value << " +- "<< errValue << " (stat)"<<endl;
 656
 657     //
 658     // Sigma systematic uncertainties
 659     //
 660     if (fAsymUncertainties && value>0.) {
 661
 662       //  (syst but feed-down) delta_sigma = sigma * sqrt ( (delta_spectra_syst/spectra)^2 +
 663       //                                     (delta_lumi/lumi)^2 + (delta_eff_trig/eff_trig)^2 + (delta_eff/eff)^2  + (global_eff)^2 )
 664       errvalueMax = value * TMath::Sqrt( (fgYieldCorr->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin))*(fgYieldCorr->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin)) +
 665                                          (fLuminosity[1]/fLuminosity[0])*(fLuminosity[1]/fLuminosity[0]) +
 666                                          (fTrigEfficiency[1]/fTrigEfficiency[0])*(fTrigEfficiency[1]/fTrigEfficiency[0])  +
 667                                          (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin)) +
 668                                          fGlobalEfficiencyUncertainties[0]*fGlobalEfficiencyUncertainties[0] );
 669       errvalueMin = value * TMath::Sqrt( (fgYieldCorr->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin))*(fgYieldCorr->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin)) +
 670                                          (fLuminosity[1]/fLuminosity[0])*(fLuminosity[1]/fLuminosity[0]) +
 671                                          (fTrigEfficiency[1]/fTrigEfficiency[0])*(fTrigEfficiency[1]/fTrigEfficiency[0])  +
 672                                          (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))  +
 673                                          fGlobalEfficiencyUncertainties[0]*fGlobalEfficiencyUncertainties[0] );
 674
 675       // Uncertainties from feed-down
 676       //      (feed-down syst) delta_sigma = sigma * sqrt ( (delta_spectra_fd/spectra_fd)^2 )
 677       //   extreme case
 678       errvalueExtremeMax = value * (fgYieldCorrExtreme->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin));
 679       errvalueExtremeMin =  value * (fgYieldCorrExtreme->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin));
 680       //
 681       //   conservative case
 682       errvalueConservativeMax = value * (fgYieldCorrConservative->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin));
 683       errvalueConservativeMin =  value * (fgYieldCorrConservative->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin));
 684
 685
 686       // stat unc of the efficiencies, separately
 687       errvalueStatUncEffc = value * (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin)) ;
 688       errvalueStatUncEffb = 0.;
 689
 690     }
 691     else {
 692       // protect against null denominator
 693       errvalueMax = (value!=0.) ?
 694         value * TMath::Sqrt( (fLuminosity[1]/fLuminosity[0])*(fLuminosity[1]/fLuminosity[0]) +
 695                              (fTrigEfficiency[1]/fTrigEfficiency[0])*(fTrigEfficiency[1]/fTrigEfficiency[0])  +
 696                              (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))  +
 697                              fGlobalEfficiencyUncertainties[0]*fGlobalEfficiencyUncertainties[0] )
 698         : 0. ;
 699       errvalueMin = errvalueMax;
 700     }
 701
 702     //
 703     // Fill the histograms
 704     //
 705     fhSigmaCorr->SetBinContent(ibin,value);
 706     fhSigmaCorr->SetBinError(ibin,errValue);
 707     //
 708     // Fill the histos and ntuple vs the Eloss hypothesis
 709     //
 710     if (fPbPbElossHypothesis) {
 711       // Loop over the Eloss hypothesis
 712       for (Float_t rval=0.0025; rval<4.0; rval+=0.005) {
 713         Int_t rbin = FindTH2YBin(fhYieldCorrRcb,rval);
 714         Double_t yieldRcbvalue = fhYieldCorrRcb->GetBinContent(ibin,rbin);
 715         // Sigma calculation
 716         //   Sigma = ( 1. / (lumi * delta_y * BR_c * eff_trig * eff_c ) ) * spectra (corrected for feed-down)
 717         Double_t sigmaRcbvalue = (fhDirectEffpt->GetBinContent(ibin) && fhDirectEffpt->GetBinContent(ibin)!=0.) ?
 718           ( yieldRcbvalue / ( deltaY * branchingRatioC * fLuminosity[0] * fTrigEfficiency[0] * fhDirectEffpt->GetBinContent(ibin) ) )
 719           : 0. ;
 720         fhSigmaCorrRcb->Fill( fhSigmaCorr->GetBinCenter(ibin) , rval, sigmaRcbvalue );
 721         //      if(ibin==3)
 722         //        cout << " pt "<< fhRECpt->GetBinCenter(ibin) <<" bin "<< ibin<<" rval="<<rval<<", rbin="<<rbin<<" fc-value="<< fhFcRcb->GetBinContent(ibin,rbin) <<", yield-fcRbvalue="<<yieldRcbvalue<<", sigma-fcRbvalue="<<sigmaRcbvalue<<endl;
 723         fnSigma->Fill(fhRECpt->GetBinCenter(ibin), fhRECpt->GetBinContent(ibin),
 724                       rval, fhFcRcb->GetBinContent(ibin,rbin),
 725                       yieldRcbvalue, sigmaRcbvalue );
 726       }
 727     }
 728     //
 729     // Fill the TGraphAsymmErrors
 730     if (fAsymUncertainties) {
 731       Double_t x = fhYieldCorr->GetBinCenter(ibin);
 732       fgSigmaCorr->SetPoint(ibin,x,value); // i,x,y
 733       fgSigmaCorr->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueMin,errvalueMax); // i,xl,xh,yl,yh
 734       fhSigmaCorrMax->SetBinContent(ibin,value+errvalueMax);
 735       fhSigmaCorrMin->SetBinContent(ibin,value-errvalueMin);
 736       fgSigmaCorrExtreme->SetPoint(ibin,x,value); // i,x,y
 737       fgSigmaCorrExtreme->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueExtremeMin,errvalueExtremeMax); // i,xl,xh,yl,yh
 738       fgSigmaCorrConservative->SetPoint(ibin,x,value); // i,x,y
 739       fgSigmaCorrConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueConservativeMin,errvalueConservativeMax); // i,xl,xh,yl,yh
 740
 741       fhStatUncEffcSigma->SetBinContent(ibin,0.);
 742       if(value>0.) fhStatUncEffcSigma->SetBinError(ibin,((errvalueStatUncEffc/value)*100.));
 743       fhStatUncEffbSigma->SetBinContent(ibin,0.); fhStatUncEffbSigma->SetBinError(ibin,0.);
 744       //      cout << " pt "<< fhRECpt->GetBinCenter(ibin) <<" bin "<< ibin<<" stat-unc-c-sigma "<< errvalueStatUncEffc/value << endl;
 745     }
 746
 747   }
 748   delete [] binwidths;
 749   delete [] limits;
 750
 751 }
 752
 753 //_________________________________________________________________________________________________________
 754 TH1D * AliHFPtSpectrum::EstimateEfficiencyRecoBin(TH1D *hSimu, TH1D *hReco, const char *name) {
 755   //
 756   // Function that computes the acceptance and efficiency correction
 757   //  based on the simulated and reconstructed spectra
 758   //  and using the reconstructed spectra bin width
 759   //
 760   //  eff = reco/sim ; err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
 761   //
 762
 763   if(!fhRECpt){
 764     AliInfo("Hey, the reconstructed histogram was not set yet !");
 765     return NULL;
 766   }
 767
 768   Int_t nbins = fhRECpt->GetNbinsX();
 769   Double_t *limits = new Double_t[nbins+1];
 770   Double_t xlow=0.,binwidth=0.;
 771   for (Int_t i=1; i<=nbins; i++) {
 772     binwidth = fhRECpt->GetBinWidth(i);
 773     xlow = fhRECpt->GetBinLowEdge(i);
 774     limits[i-1] = xlow;
 775   }
 776   limits[nbins] = xlow + binwidth;
 777
 778   TH1D * hEfficiency = new TH1D(name," acceptance #times efficiency",nbins,limits);
 779
 780   Double_t *sumSimu=new Double_t[nbins];
 781   Double_t *sumReco=new Double_t[nbins];
 782   for (Int_t ibin=0; ibin<nbins; ibin++){
 783     sumSimu[ibin]=0.;  sumReco[ibin]=0.;
 784   }
 785   for (Int_t ibin=0; ibin<=hSimu->GetNbinsX(); ibin++){
 786
 787     for (Int_t ibinrec=0; ibinrec<nbins; ibinrec++){
 788       if ( hSimu->GetBinCenter(ibin)>limits[ibinrec] &&
 789            hSimu->GetBinCenter(ibin)<limits[ibinrec+1] ) {
 790         sumSimu[ibinrec]+=hSimu->GetBinContent(ibin);
 791       }
 792       if ( hReco->GetBinCenter(ibin)>limits[ibinrec] &&
 793            hReco->GetBinCenter(ibin)<limits[ibinrec+1] ) {
 794         sumReco[ibinrec]+=hReco->GetBinContent(ibin);
 795       }
 796     }
 797
 798   }
 799
 800
 801   // the efficiency is computed as reco/sim (in each bin)
 802   //  its uncertainty is err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
 803   Double_t eff=0., erreff=0.;
 804   for (Int_t ibinrec=0; ibinrec<nbins; ibinrec++) {
 805     if (sumSimu[ibinrec]!= 0. && sumReco[ibinrec]!=0.) {
 806       eff = sumReco[ibinrec] / sumSimu[ibinrec] ;
 807       // protection in case eff > 1.0
 808       // test calculation (make the argument of the sqrt positive)
 809       erreff = TMath::Sqrt( eff * TMath::Abs(1.0 - eff) ) / TMath::Sqrt( sumSimu[ibinrec] );
 810     }
 811     else { eff=0.0; erreff=0.; }
 812     hEfficiency->SetBinContent(ibinrec+1,eff);
 813     hEfficiency->SetBinError(ibinrec+1,erreff);
 814   }
 815
 816   delete [] sumSimu;
 817   delete [] sumReco;
 818
 819   return (TH1D*)hEfficiency;
 820 }
 821
 822 //_________________________________________________________________________________________________________
 823 void AliHFPtSpectrum::EstimateAndSetDirectEfficiencyRecoBin(TH1D *hSimu, TH1D *hReco) {
 824   //
 825   // Function that computes the Direct  acceptance and efficiency correction
 826   //  based on the simulated and reconstructed spectra
 827   //  and using the reconstructed spectra bin width
 828   //
 829   //  eff = reco/sim ; err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
 830   //
 831
 832   if(!fhRECpt || !hSimu || !hReco){
 833     AliError("Hey, the reconstructed histogram was not set yet !");
 834     return;
 835   }
 836
 837   fhDirectEffpt = EstimateEfficiencyRecoBin(hSimu,hReco,"fhDirectEffpt");
 838   fhDirectEffpt->SetNameTitle("fhDirectEffpt"," direct acceptance #times efficiency");
 839
 840 }
 841
 842 //_________________________________________________________________________________________________________
 843 void AliHFPtSpectrum::EstimateAndSetFeedDownEfficiencyRecoBin(TH1D *hSimu, TH1D *hReco) {
 844   //
 845   // Function that computes the Feed-Down acceptance and efficiency correction
 846   //  based on the simulated and reconstructed spectra
 847   //  and using the reconstructed spectra bin width
 848   //
 849   //  eff = reco/sim ; err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
 850   //
 851
 852   if(!fhRECpt || !hSimu || !hReco){
 853     AliError("Hey, the reconstructed histogram was not set yet !");
 854     return;
 855   }
 856
 857   fhFeedDownEffpt = EstimateEfficiencyRecoBin(hSimu,hReco,"fhFeedDownEffpt");
 858   fhFeedDownEffpt->SetNameTitle("fhFeedDownEffpt"," feed-down acceptance #times efficiency");
 859
 860 }
 861
 862 //_________________________________________________________________________________________________________
 863 Bool_t AliHFPtSpectrum::Initialize(){
 864   //
 865   // Initialization of the variables (histograms)
 866   //
 867
 868   if (fFeedDownOption==0) {
 869     AliInfo("Getting ready for the corrections without feed-down consideration");
 870   } else if (fFeedDownOption==1) {
 871     AliInfo("Getting ready for the fc feed-down correction calculation");
 872   } else if (fFeedDownOption==2) {
 873     AliInfo("Getting ready for the Nb feed-down correction calculation");
 874   } else { AliError("The calculation option must be <=2"); return kFALSE; }
 875
 876   // Start checking the input histograms consistency
 877   Bool_t areconsistent=kTRUE;
 878
 879   // General checks
 880   if (!fhDirectEffpt || !fhRECpt) {
 881     AliError(" Reconstructed spectra and/or the Nc efficiency distributions are not defined");
 882     return kFALSE;
 883   }
 884   areconsistent &= CheckHistosConsistency(fhRECpt,fhDirectEffpt);
 885   if (!areconsistent) {
 886     AliInfo("Histograms required for Nb correction are not consistent (bin width, bounds)");
 887     return kFALSE;
 888   }
 889   if (fFeedDownOption==0) return kTRUE;
 890
 891   //
 892   // Common checks for options 1 (fc) & 2(Nb)
 893   if (!fhFeedDownMCpt || !fhFeedDownEffpt) {
 894     AliError(" Theoretical Nb and/or the Nb efficiency distributions are not defined");
 895     return kFALSE;
 896   }
 897   areconsistent &= CheckHistosConsistency(fhRECpt,fhFeedDownMCpt);
 898   areconsistent &= CheckHistosConsistency(fhFeedDownMCpt,fhFeedDownEffpt);
 899   if (fAsymUncertainties) {
 900     if (!fhFeedDownMCptMax || !fhFeedDownMCptMin) {
 901       AliError(" Max/Min theoretical Nb distributions are not defined");
 902       return kFALSE;
 903     }
 904     areconsistent &= CheckHistosConsistency(fhFeedDownMCpt,fhFeedDownMCptMax);
 905   }
 906   if (!areconsistent) {
 907     AliInfo("Histograms required for Nb correction are not consistent (bin width, bounds)");
 908     return kFALSE;
 909   }
 910   if (fFeedDownOption>1) return kTRUE;
 911
 912   //
 913   // Now checks for option 1 (fc correction)
 914   if (!fhDirectMCpt) {
 915     AliError("Theoretical Nc distributions is not defined");
 916     return kFALSE;
 917   }
 918   areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhFeedDownMCpt);
 919   areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhDirectEffpt);
 920   if (fAsymUncertainties) {
 921     if (!fhDirectMCptMax || !fhDirectMCptMin) {
 922       AliError(" Max/Min theoretical Nc distributions are not defined");
 923       return kFALSE;
 924     }
 925     areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhDirectMCptMax);
 926   }
 927   if (!areconsistent) {
 928     AliInfo("Histograms required for fc correction are not consistent (bin width, bounds)");
 929     return kFALSE;
 930   }
 931
 932   return kTRUE;
 933 }
 934
 935 //_________________________________________________________________________________________________________
 936 Bool_t AliHFPtSpectrum::CheckHistosConsistency(TH1D *h1, TH1D *h2){
 937   //
 938   // Check the histograms consistency (bins, limits)
 939   //
 940
 941   if (!h1 || !h2) {
 942     AliError("One or both histograms don't exist");
 943     return kFALSE;
 944   }
 945
 946   Double_t binwidth1 = h1->GetBinWidth(1);
 947   Double_t binwidth2 = h2->GetBinWidth(1);
 948   Double_t min1 = h1->GetBinCenter(1) - (binwidth1/2.) ;
 949 //   Double_t max1 = h1->GetBinCenter(nbins1) + (binwidth1/2.) ;
 950   Double_t min2 = h2->GetBinCenter(1) - (binwidth2/2.) ;
 951 //   Double_t max2 = h2->GetBinCenter(nbins2) + (binwidth2/2.) ;
 952
 953   if (binwidth1!=binwidth2) {
 954     AliInfo(" histograms with different bin width");
 955     return kFALSE;
 956   }
 957   if (min1!=min2) {
 958     AliInfo(" histograms with different minimum");
 959     return kFALSE;
 960   }
 961 //   if (max1!=max2) {
 962 //     AliInfo(" histograms with different maximum");
 963 //     return kFALSE;
 964 //   }
 965
 966   return kTRUE;
 967 }
 968
 969 //_________________________________________________________________________________________________________
 970 void AliHFPtSpectrum::CalculateFeedDownCorrectionFc(){
 971   //
 972   // Compute fc factor and its uncertainties bin by bin
 973   //   fc = 1 / ( 1 + (eff_b/eff_c)*(N_b/N_c) )
 974   //
 975   // uncertainties: (conservative) combine the upper/lower N_b & N_c predictions together
 976   //                (extreme) combine the upper N_b predictions with the lower N_c predictions & viceversa
 977   //                systematic uncertainty on the acceptance x efficiency b/c ratio are included
 978   //
 979   //  In addition, in HIC the feed-down correction varies with an energy loss hypothesis: Raa(c-->D) / Raa(b-->D) = Rcb
 980   //           fc (Rcb) = ( 1. / ( 1 + (eff_b/eff_c)*(N_b/N_c)* (1/Rcb) ) );
 981   //
 982   AliInfo("Calculating the feed-down correction factor (fc method)");
 983
 984   // define the variables
 985   Int_t nbins = fhRECpt->GetNbinsX();
 986   Double_t binwidth = fhRECpt->GetBinWidth(1);
 987   Double_t *limits = new Double_t[nbins+1];
 988   Double_t *binwidths = new Double_t[nbins];
 989   Double_t xlow=0.;
 990   for (Int_t i=1; i<=nbins; i++) {
 991     binwidth = fhRECpt->GetBinWidth(i);
 992     xlow = fhRECpt->GetBinLowEdge(i);
 993     limits[i-1] = xlow;
 994     binwidths[i-1] = binwidth;
 995   }
 996   limits[nbins] = xlow + binwidth;
 997
 998   Double_t correction=1.;
 999   Double_t theoryRatio=1.;
1000   Double_t effRatio=1.;
1001   Double_t correctionExtremeA=1., correctionExtremeB=1.;
1002   Double_t theoryRatioExtremeA=1., theoryRatioExtremeB=1.;
1003   Double_t correctionConservativeA=1., correctionConservativeB=1.;
1004   Double_t theoryRatioConservativeA=1., theoryRatioConservativeB=1.;
1005   Double_t correctionUnc=0.;
1006   Double_t correctionExtremeAUnc=0., correctionExtremeBUnc=0.;
1007   Double_t correctionConservativeAUnc=0., correctionConservativeBUnc=0.;
1008
1009   // declare the output histograms
1010   fhFc = new TH1D("fhFc","fc correction factor",nbins,limits);
1011   fhFcMax = new TH1D("fhFcMax","max fc correction factor",nbins,limits);
1012   fhFcMin = new TH1D("fhFcMin","min fc correction factor",nbins,limits);
1013   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.);
1014   // two local control histograms
1015   TH1D *hTheoryRatio = new TH1D("hTheoryRatio","Theoretical B-->D over c-->D (feed-down/direct) ratio",nbins,limits);
1016   TH1D *hEffRatio = new TH1D("hEffRatio","Efficiency B-->D over c-->D (feed-down/direct) ratio",nbins,limits);
1017   // and the output TGraphAsymmErrors
1018   if (fAsymUncertainties) {
1019     fgFcExtreme = new TGraphAsymmErrors(nbins+1);
1020     fgFcExtreme->SetNameTitle("fgFcExtreme","fgFcExtreme");
1021     fgFcConservative = new TGraphAsymmErrors(nbins+1);
1022     fgFcConservative->SetNameTitle("fgFcConservative","fgFcConservative");
1023   }
1024
1025   fhStatUncEffcFD = new TH1D("fhStatUncEffcFD","direct charm stat unc on the feed-down correction",nbins,limits);
1026   fhStatUncEffbFD = new TH1D("fhStatUncEffbFD","secondary charm stat unc on the feed-down correction",nbins,limits);
1027   Double_t correctionConservativeAUncStatEffc=0., correctionConservativeBUncStatEffc=0.;
1028   Double_t correctionConservativeAUncStatEffb=0., correctionConservativeBUncStatEffb=0.;
1029
1030   //
1031   // Compute fc
1032   //
1033   for (Int_t ibin=1; ibin<=nbins; ibin++) {
1034
1035     //  theory_ratio = (N_b/N_c)
1036     theoryRatio = (fhDirectMCpt->GetBinContent(ibin)>0. && fhFeedDownMCpt->GetBinContent(ibin)>0.) ?
1037       fhFeedDownMCpt->GetBinContent(ibin) / fhDirectMCpt->GetBinContent(ibin) : 1.0 ;
1038
1039     //
1040     // Calculate the uncertainty [ considering only the theoretical uncertainties on Nb & Nc for now !!! ]
1041     //
1042     // extreme A = direct-max, feed-down-min
1043     theoryRatioExtremeA = (fhDirectMCptMax->GetBinContent(ibin)>0. && fhFeedDownMCptMin->GetBinContent(ibin)>0.) ?
1044       fhFeedDownMCptMin->GetBinContent(ibin) / fhDirectMCptMax->GetBinContent(ibin) : 1.0 ;
1045     // extreme B = direct-min, feed-down-max
1046     theoryRatioExtremeB = (fhDirectMCptMin->GetBinContent(ibin)>0. && fhDirectMCptMax->GetBinContent(ibin)>0.) ?
1047       fhFeedDownMCptMax->GetBinContent(ibin) / fhDirectMCptMin->GetBinContent(ibin) : 1.0 ;
1048     // conservative A = direct-max, feed-down-max
1049     theoryRatioConservativeA = (fhDirectMCptMax->GetBinContent(ibin)>0. && fhFeedDownMCptMin->GetBinContent(ibin)>0.) ?
1050       fhFeedDownMCptMax->GetBinContent(ibin) / fhDirectMCptMax->GetBinContent(ibin) : 1.0 ;
1051     // conservative B = direct-min, feed-down-min
1052     theoryRatioConservativeB = (fhDirectMCptMin->GetBinContent(ibin)>0. && fhDirectMCptMax->GetBinContent(ibin)>0.) ?
1053       fhFeedDownMCptMin->GetBinContent(ibin) / fhDirectMCptMin->GetBinContent(ibin) : 1.0 ;
1054
1055     //  eff_ratio = (eff_b/eff_c)
1056     effRatio = (fhDirectEffpt->GetBinContent(ibin) && fhDirectEffpt->GetBinContent(ibin)!=0.) ?
1057       fhFeedDownEffpt->GetBinContent(ibin) / fhDirectEffpt->GetBinContent(ibin) : 1.0 ;
1058
1059     //   fc = 1 / ( 1 + (eff_b/eff_c)*(N_b/N_c) )
1060     if( TMath::Abs(effRatio - 1.0)<0.01 || TMath::Abs(theoryRatio - 1.0)<0.01 ) {
1061       correction = 1.0;
1062       correctionExtremeA = 1.0;
1063       correctionExtremeB = 1.0;
1064       correctionConservativeA = 1.0;
1065       correctionConservativeB = 1.0;
1066     }
1067     else {
1068       correction = ( 1. / ( 1 + ( effRatio * theoryRatio ) ) );
1069       correctionExtremeA = ( 1. / ( 1 + ( effRatio * theoryRatioExtremeA ) ) );
1070       correctionExtremeB = ( 1. / ( 1 + ( effRatio * theoryRatioExtremeB ) ) );
1071       correctionConservativeA = ( 1. / ( 1 + ( effRatio * theoryRatioConservativeA ) ) );
1072       correctionConservativeB = ( 1. / ( 1 + ( effRatio * theoryRatioConservativeB ) ) );
1073     }
1074
1075
1076     // fc uncertainty from (eff_b/eff_c) = fc^2 * (N_b/N_c) * delta(eff_b/eff_c)
1077     //  delta(eff_b/eff_c) is a percentage = effRatio * sqrt( fGlobalEfficiencyUncertainties[1]^2 + unc_eff_c ^2 + unc_eff_b ^2 )
1078     correctionUnc = correction*correction * theoryRatio * effRatio *
1079       TMath::Sqrt( fGlobalEfficiencyUncertainties[1]*fGlobalEfficiencyUncertainties[1] +
1080                    (fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) +
1081                    (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))
1082                    );
1083     correctionExtremeAUnc = correctionExtremeA*correctionExtremeA * theoryRatioExtremeA  * effRatio *
1084       TMath::Sqrt( fGlobalEfficiencyUncertainties[1]*fGlobalEfficiencyUncertainties[1] +
1085                    (fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) +
1086                    (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))
1087                    );
1088     correctionExtremeBUnc = correctionExtremeB*correctionExtremeB * theoryRatioExtremeB  * effRatio *
1089       TMath::Sqrt( fGlobalEfficiencyUncertainties[1]*fGlobalEfficiencyUncertainties[1] +
1090                    (fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) +
1091                    (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))
1092                    );
1093     correctionConservativeAUnc = correctionConservativeA*correctionConservativeA * theoryRatioConservativeA  *effRatio *
1094       TMath::Sqrt( fGlobalEfficiencyUncertainties[1]*fGlobalEfficiencyUncertainties[1] +
1095                    (fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) +
1096                    (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))
1097                    );
1098     //
1099     correctionConservativeAUncStatEffc = correctionConservativeA*correctionConservativeA * theoryRatioConservativeA  *effRatio *
1100       (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin));
1101     correctionConservativeAUncStatEffb = correctionConservativeA*correctionConservativeA * theoryRatioConservativeA  *effRatio *
1102       (fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin));
1103
1104     correctionConservativeBUnc = correctionConservativeB*correctionConservativeB * theoryRatioConservativeB  *effRatio *
1105       TMath::Sqrt( fGlobalEfficiencyUncertainties[1]*fGlobalEfficiencyUncertainties[1] +
1106                    (fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) +
1107                    (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))
1108                    );
1109     correctionConservativeBUncStatEffb = correctionConservativeB*correctionConservativeB * theoryRatioConservativeB  *effRatio *
1110       (fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin));
1111     correctionConservativeBUncStatEffc = correctionConservativeB*correctionConservativeB * theoryRatioConservativeB  *effRatio *
1112       (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin));
1113
1114
1115     // Fill in the histograms
1116     hTheoryRatio->SetBinContent(ibin,theoryRatio);
1117     hEffRatio->SetBinContent(ibin,effRatio);
1118     fhFc->SetBinContent(ibin,correction);
1119     //
1120     // Estimate how the result varies vs charm/beauty Eloss hypothesis
1121     //
1122     if ( TMath::Abs(correction-1.0)>0.01 && fPbPbElossHypothesis){
1123       // Loop over the Eloss hypothesis
1124       //      Int_t rbin=0;
1125       for (Float_t rval=0.0025; rval<4.0; rval+=0.005){
1126         Double_t correctionRcb = ( 1. / ( 1 + ( effRatio * theoryRatio * (1/rval) ) ) );
1127         fhFcRcb->Fill( fhFc->GetBinCenter(ibin) , rval, correctionRcb );
1128         //      if(ibin==3){
1129         //        cout << " pt "<< fhFc->GetBinCenter(ibin) <<" bin "<< ibin<<" rval="<<rval<<", rbin="<<rbin<<", fc-Rcb-value="<<correctionRcb<<endl;
1130         //        rbin++;
1131         //      }
1132       }
1133     }
1134     //
1135     // Fill the rest of (asymmetric) histograms
1136     //
1137     if (fAsymUncertainties) {
1138       Double_t x = fhDirectMCpt->GetBinCenter(ibin);
1139       Double_t val[4] = { correctionExtremeA + correctionExtremeAUnc, correctionExtremeA - correctionExtremeAUnc,
1140                           correctionExtremeB + correctionExtremeBUnc, correctionExtremeB - correctionExtremeBUnc };
1141       Double_t uncExtremeMin = correction - TMath::MinElement(4,val);
1142       Double_t uncExtremeMax = TMath::MaxElement(4,val) - correction;
1143       fgFcExtreme->SetPoint(ibin,x,correction); // i,x,y
1144       fgFcExtreme->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),uncExtremeMin,uncExtremeMax); // i,xl,xh,yl,yh
1145       fhFcMax->SetBinContent(ibin,correction+uncExtremeMax);
1146       fhFcMin->SetBinContent(ibin,correction-uncExtremeMin);
1147       Double_t consval[4] = { correctionConservativeA - correctionConservativeAUnc, correctionConservativeA + correctionConservativeAUnc,
1148                               correctionConservativeB - correctionConservativeBUnc, correctionConservativeB + correctionConservativeBUnc};
1149       Double_t uncConservativeMin = correction - TMath::MinElement(4,consval);
1150       Double_t uncConservativeMax =  TMath::MaxElement(4,consval) - correction;
1151       fgFcConservative->SetPoint(ibin,x,correction); // i,x,y
1152       fgFcConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),uncConservativeMin,uncConservativeMax); // i,xl,xh,yl,yh
1153       if( !(correction>0.) ){
1154         fgFcExtreme->SetPoint(ibin,x,0.); // i,x,y
1155         fgFcExtreme->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),0.,0.); // i,xl,xh,yl,yh
1156         fgFcConservative->SetPoint(ibin,x,0.); // i,x,y
1157         fgFcConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),0.,0.); // i,xl,xh,yl,yh
1158       }
1159
1160       Double_t valStatEffc[2] = { correctionConservativeAUncStatEffc/correctionConservativeA,
1161                                   correctionConservativeBUncStatEffc/correctionConservativeB };
1162       Double_t valStatEffb[2] = { correctionConservativeAUncStatEffb/correctionConservativeA,
1163                                   correctionConservativeBUncStatEffb/correctionConservativeB };
1164       Double_t uncConservativeStatEffc = TMath::MaxElement(2,valStatEffc);
1165       Double_t uncConservativeStatEffb = TMath::MaxElement(2,valStatEffb);
1166       fhStatUncEffcFD->SetBinContent(ibin,0.); fhStatUncEffcFD->SetBinError(ibin,uncConservativeStatEffc*100.);
1167       fhStatUncEffbFD->SetBinContent(ibin,0.); fhStatUncEffbFD->SetBinError(ibin,uncConservativeStatEffb*100.);
1168       //      cout << " pt "<< fhStatUncEffcFD->GetBinCenter(ibin) <<" bin "<< ibin<<" fc-stat-c ="<<uncConservativeStatEffc<<" fc-stat-b ="<<uncConservativeStatEffb<<endl;
1169
1170     }
1171
1172   }
1173   delete [] binwidths;
1174   delete [] limits;
1175
1176 }
1177
1178 //_________________________________________________________________________________________________________
1179 void AliHFPtSpectrum::CalculateFeedDownCorrectedSpectrumFc(){
1180   //
1181   // Compute the feed-down corrected spectrum if feed-down correction is done via fc factor (bin by bin)
1182   //    physics = reco * fc / bin-width
1183   //
1184   //    uncertainty:             (stat) delta_physics = physics * sqrt ( (delta_reco/reco)^2 )
1185   //               (syst but feed-down) delta_physics = physics * sqrt ( (delta_reco_syst/reco)^2 )
1186   //                   (feed-down syst) delta_physics = physics * sqrt ( (delta_fc/fc)^2 )
1187   //
1188   //    ( Calculation done bin by bin )
1189   //
1190   //  In addition, in HIC the feed-down correction varies with an energy loss hypothesis: Raa(c-->D) / Raa(b-->D) = Rcb
1191
1192   AliInfo(" Calculating the feed-down corrected spectrum (fc method)");
1193
1194   if (!fhFc || !fhRECpt) {
1195     AliError(" Reconstructed or fc distributions are not defined");
1196     return;
1197   }
1198
1199   Int_t nbins = fhRECpt->GetNbinsX();
1200   Double_t value = 0., errvalue = 0., errvalueMax= 0., errvalueMin= 0.;
1201   Double_t valueExtremeMax= 0., valueExtremeMin= 0.;
1202   Double_t valueConservativeMax= 0., valueConservativeMin= 0.;
1203   Double_t binwidth = fhRECpt->GetBinWidth(1);
1204   Double_t *limits = new Double_t[nbins+1];
1205   Double_t *binwidths = new Double_t[nbins];
1206   Double_t xlow=0.;
1207   for (Int_t i=1; i<=nbins; i++) {
1208     binwidth = fhRECpt->GetBinWidth(i);
1209     xlow = fhRECpt->GetBinLowEdge(i);
1210     limits[i-1] = xlow;
1211     binwidths[i-1] = binwidth;
1212   }
1213   limits[nbins] = xlow + binwidth;
1214
1215   // declare the output histograms
1216   fhYieldCorr = new TH1D("fhYieldCorr","corrected yield (by fc)",nbins,limits);
1217   fhYieldCorrMax = new TH1D("fhYieldCorrMax","max corrected yield (by fc)",nbins,limits);
1218   fhYieldCorrMin = new TH1D("fhYieldCorrMin","min corrected yield (by fc)",nbins,limits);
1219   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.);
1220   // and the output TGraphAsymmErrors
1221   if (fAsymUncertainties){
1222     fgYieldCorr = new TGraphAsymmErrors(nbins+1);
1223     fgYieldCorrExtreme = new TGraphAsymmErrors(nbins+1);
1224     fgYieldCorrConservative = new TGraphAsymmErrors(nbins+1);
1225   }
1226
1227   //
1228   // Do the calculation
1229   //
1230   for (Int_t ibin=1; ibin<=nbins; ibin++) {
1231
1232     // calculate the value
1233     //    physics = reco * fc / bin-width
1234     value = (fhRECpt->GetBinContent(ibin) && fhFc->GetBinContent(ibin)) ?
1235       fhRECpt->GetBinContent(ibin) * fhFc->GetBinContent(ibin) : 0. ;
1236     value /= fhRECpt->GetBinWidth(ibin) ;
1237
1238     // Statistical uncertainty
1239     //    (stat) delta_physics = physics * sqrt ( (delta_reco/reco)^2 )
1240     errvalue = (value!=0. && fhRECpt->GetBinContent(ibin) && fhRECpt->GetBinContent(ibin)!=0.) ?
1241       value * (fhRECpt->GetBinError(ibin)/fhRECpt->GetBinContent(ibin)) : 0. ;
1242
1243     // Calculate the systematic uncertainties
1244     //    (syst but feed-down) delta_physics = physics * sqrt ( (delta_reco_syst/reco)^2 )
1245     //        (feed-down syst) delta_physics = physics * sqrt ( (delta_fc/fc)^2 )
1246     //
1247     // Protect against null denominator. If so, define uncertainty as null
1248     if (fhRECpt->GetBinContent(ibin) && fhRECpt->GetBinContent(ibin)!=0.) {
1249
1250       if (fAsymUncertainties) {
1251
1252         // Systematics but feed-down
1253         if (fgRECSystematics) {
1254           errvalueMax = value * ( fgRECSystematics->GetErrorYhigh(ibin) / fhRECpt->GetBinContent(ibin) );
1255           errvalueMin = value * ( fgRECSystematics->GetErrorYlow(ibin) / fhRECpt->GetBinContent(ibin) );
1256         }
1257         else { errvalueMax = 0.; errvalueMin = 0.; }
1258
1259         // Extreme feed-down systematics
1260         valueExtremeMax = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) + fgFcExtreme->GetErrorYhigh(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
1261         valueExtremeMin = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) - fgFcExtreme->GetErrorYlow(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
1262
1263         // Conservative feed-down systematics
1264         valueConservativeMax = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) + fgFcConservative->GetErrorYhigh(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
1265         valueConservativeMin = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) - fgFcConservative->GetErrorYlow(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
1266
1267       }
1268
1269     }
1270     else { errvalueMax = 0.; errvalueMin = 0.; }
1271
1272     //
1273     // Fill in the histograms
1274     //
1275     fhYieldCorr->SetBinContent(ibin,value);
1276     fhYieldCorr->SetBinError(ibin,errvalue);
1277     //
1278     // Fill the histos and ntuple vs the Eloss hypothesis
1279     //
1280     if (fPbPbElossHypothesis) {
1281       // Loop over the Eloss hypothesis
1282       for (Float_t rval=0.0025; rval<4.0; rval+=0.005){
1283         Int_t rbin = FindTH2YBin(fhYieldCorrRcb,rval);
1284         Double_t fcRcbvalue = fhFcRcb->GetBinContent(ibin,rbin);
1285         //    physics = reco * fcRcb / bin-width
1286         Double_t Rcbvalue = (fhRECpt->GetBinContent(ibin) && fcRcbvalue) ?
1287           fhRECpt->GetBinContent(ibin) * fcRcbvalue : 0. ;
1288         Rcbvalue /= fhRECpt->GetBinWidth(ibin) ;
1289         fhYieldCorrRcb->Fill( fhYieldCorr->GetBinCenter(ibin) , rval, Rcbvalue );
1290         //        cout << " pt "<< fhRECpt->GetBinCenter(ibin) <<" bin "<< ibin<<" rval="<<rval<<", rbin="<<rbin<<" fc-fcRbvalue="<<fcRcbvalue<<", yield="<<Rcbvalue<<endl;
1291       }
1292     }
1293     if (fAsymUncertainties) {
1294       Double_t center = fhYieldCorr->GetBinCenter(ibin);
1295       fgYieldCorr->SetPoint(ibin,center,value); // i,x,y
1296       fgYieldCorr->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueMin,errvalueMax); // i,xl,xh,yl,yh
1297       fhYieldCorrMax->SetBinContent(ibin,value+errvalueMax);
1298       fhYieldCorrMin->SetBinContent(ibin,value-errvalueMin);
1299       fgYieldCorrExtreme->SetPoint(ibin,center,value);
1300       fgYieldCorrExtreme->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),value-valueExtremeMin,valueExtremeMax-value);
1301       fgYieldCorrConservative->SetPoint(ibin,center,value);
1302       fgYieldCorrConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),value-valueConservativeMin,valueConservativeMax-value);
1303     }
1304
1305   }
1306   delete [] binwidths;
1307   delete [] limits;
1308
1309 }
1310
1311 //_________________________________________________________________________________________________________
1312 void AliHFPtSpectrum::CalculateFeedDownCorrectedSpectrumNb(Double_t deltaY, Double_t branchingRatioBintoFinalDecay) {
1313   //
1314   // Compute the feed-down corrected spectrum if feed-down correction is done via Nb (bin by bin)
1315   //    physics =  [ reco  - (lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th) ] / bin-width
1316   //
1317   //    uncertainty:   (stat)  delta_physics = sqrt ( (delta_reco)^2 )  / bin-width
1318   //     (syst but feed-down)  delta_physics = sqrt ( (delta_reco_syst)^2 )  / bin-width
1319   //         (feed-down syst)  delta_physics = sqrt ( (k*delta_lumi/lumi)^2 + (k*delta_eff_trig/eff_trig)^2
1320   //                                                   + (k*delta_Nb/Nb)^2 + (k*delta_eff/eff)^2  + (k*global_eff_ratio)^2 ) / bin-width
1321   //                    where k = lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th
1322   //
1323   //  In addition, in HIC the feed-down correction varies with an energy loss hypothesis: Raa(b-->D) = Rb
1324   //    physics =  [ reco  - ( Tab * Nevt * delta_y * BR_b * eff_trig * eff_b * Nb_th * Rb ) ] / bin-width
1325   //
1326   AliInfo("Calculating the feed-down correction factor and spectrum (Nb method)");
1327
1328   Int_t nbins = fhRECpt->GetNbinsX();
1329   Double_t binwidth = fhRECpt->GetBinWidth(1);
1330   Double_t value = 0., errvalue = 0., errvalueMax = 0., errvalueMin = 0., kfactor = 0.;
1331   Double_t errvalueExtremeMax = 0., errvalueExtremeMin = 0.;
1332   Double_t *limits = new Double_t[nbins+1];
1333   Double_t *binwidths = new Double_t[nbins];
1334   Double_t xlow=0.;
1335   for (Int_t i=1; i<=nbins; i++) {
1336     binwidth = fhRECpt->GetBinWidth(i);
1337     xlow = fhRECpt->GetBinLowEdge(i);
1338     limits[i-1] = xlow;
1339     binwidths[i-1] = binwidth;
1340   }
1341   limits[nbins] = xlow + binwidth;
1342
1343   // declare the output histograms
1344   fhYieldCorr = new TH1D("fhYieldCorr","corrected yield (by Nb)",nbins,limits);
1345   fhYieldCorrMax = new TH1D("fhYieldCorrMax","max corrected yield (by Nb)",nbins,limits);
1346   fhYieldCorrMin = new TH1D("fhYieldCorrMin","min corrected yield (by Nb)",nbins,limits);
1347   if(fPbPbElossHypothesis) {
1348     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.);
1349     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.);
1350   }
1351   // and the output TGraphAsymmErrors
1352   if (fAsymUncertainties){
1353     fgYieldCorr = new TGraphAsymmErrors(nbins+1);
1354     fgYieldCorrExtreme = new TGraphAsymmErrors(nbins+1);
1355     fgYieldCorrConservative = new TGraphAsymmErrors(nbins+1);
1356     // Define fc-conservative
1357     fgFcConservative = new TGraphAsymmErrors(nbins+1);
1358     AliInfo(" Beware the conservative & extreme uncertainties are equal by definition !");
1359   }
1360
1361   // variables to define fc-conservative
1362   double correction=0, correctionMax=0., correctionMin=0.;
1363
1364   fhStatUncEffcFD = new TH1D("fhStatUncEffcFD","direct charm stat unc on the feed-down correction",nbins,limits);
1365   fhStatUncEffbFD = new TH1D("fhStatUncEffbFD","secondary charm stat unc on the feed-down correction",nbins,limits);
1366   Double_t correctionUncStatEffc=0.;
1367   Double_t correctionUncStatEffb=0.;
1368
1369
1370   //
1371   // Do the calculation
1372   //
1373   for (Int_t ibin=1; ibin<=nbins; ibin++) {
1374
1375     // Calculate the value
1376     //    physics =  [ reco  - (lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th) ] / bin-width
1377     //  In HIC :   physics =  [ reco  - ( Tab * Nevt * delta_y * BR_b * eff_trig * eff_b * Nb_th * Rb ) ] / bin-width
1378     //
1379     //
1380     Double_t frac = 1.0, errfrac =0.;
1381     if(fPbPbElossHypothesis) {
1382       frac = fTab[0]*fNevts;
1383       errfrac = frac * TMath::Sqrt( (fTab[1]/fTab[0])*(fTab[1]/fTab[0]) + (1/fNevts) );
1384     } else {
1385       frac = fLuminosity[0];
1386       errfrac = fLuminosity[1];
1387     }
1388
1389     value = ( fhRECpt->GetBinContent(ibin)>0. && fhRECpt->GetBinContent(ibin)!=0. &&
1390               fhFeedDownMCpt->GetBinContent(ibin)>0. && fhFeedDownEffpt->GetBinContent(ibin)>0. ) ?
1391       fhRECpt->GetBinContent(ibin) - frac*(deltaY*branchingRatioBintoFinalDecay*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) * fhRECpt->GetBinWidth(ibin) )
1392       : 0. ;
1393     value /= fhRECpt->GetBinWidth(ibin);
1394     if (value<0.) value =0.;
1395
1396     //  Statistical uncertainty:   delta_physics = sqrt ( (delta_reco)^2 )  / bin-width
1397     errvalue = (value!=0. && fhRECpt->GetBinError(ibin) && fhRECpt->GetBinError(ibin)!=0.)  ?
1398       fhRECpt->GetBinError(ibin) : 0.;
1399     errvalue /= fhRECpt->GetBinWidth(ibin);
1400
1401     // Correction (fc) : Estimate of the relative amount feed-down subtracted
1402     // correction =  [ 1  - (lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th)/reco ]
1403     // in HIC: correction =  [ 1  - ( Tab * Nevt * delta_y * BR_b * eff_trig * eff_b * Nb_th)/reco ]
1404     correction = (value>0.) ?
1405       1 - (frac*deltaY*branchingRatioBintoFinalDecay*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) * fhRECpt->GetBinWidth(ibin) ) / fhRECpt->GetBinContent(ibin)  : 0. ;
1406     if (correction<0.) correction = 0.;
1407
1408     // Systematic uncertainties
1409     //     (syst but feed-down)  delta_physics = sqrt ( (delta_reco_syst)^2 )  / bin-width
1410     //         (feed-down syst)  delta_physics = sqrt ( (k*delta_lumi/lumi)^2 + (k*delta_eff_trig/eff_trig)^2
1411     //                                                   + (k*delta_Nb/Nb)^2 + (k*delta_eff/eff)^2 + (k*global_eff_ratio)^2 ) / bin-width
1412     //                    where k = lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th * bin-width
1413     kfactor = frac*deltaY*branchingRatioBintoFinalDecay*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) * fhRECpt->GetBinWidth(ibin) ;
1414     //
1415     if (fAsymUncertainties && value>0.) {
1416       Double_t nb =  fhFeedDownMCpt->GetBinContent(ibin);
1417       Double_t nbDmax = fhFeedDownMCptMax->GetBinContent(ibin) - fhFeedDownMCpt->GetBinContent(ibin);
1418       Double_t nbDmin = fhFeedDownMCpt->GetBinContent(ibin) - fhFeedDownMCptMin->GetBinContent(ibin);
1419
1420       // Systematics but feed-down
1421       if (fgRECSystematics){
1422         errvalueMax = fgRECSystematics->GetErrorYhigh(ibin) / fhRECpt->GetBinWidth(ibin) ;
1423         errvalueMin = fgRECSystematics->GetErrorYlow(ibin) / fhRECpt->GetBinWidth(ibin);
1424       }
1425       else { errvalueMax = 0.; errvalueMin = 0.; }
1426
1427       // Feed-down systematics
1428       // min value with the maximum Nb
1429       errvalueExtremeMin = TMath::Sqrt( ( (kfactor*errfrac/frac)*(kfactor*errfrac/frac) ) +
1430                                         ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) +
1431                                         ( (kfactor*nbDmax/nb)*(kfactor*nbDmax/nb) )  +
1432                                         ( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) ) +
1433                                         ( (kfactor*fGlobalEfficiencyUncertainties[1])*(kfactor*fGlobalEfficiencyUncertainties[1]) )
1434                                         ) / fhRECpt->GetBinWidth(ibin);
1435       // max value with the minimum Nb
1436       errvalueExtremeMax =  TMath::Sqrt( ( (kfactor*errfrac/frac)*(kfactor*errfrac/frac) ) +
1437                                          ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) +
1438                                          ( (kfactor*nbDmin/nb)*(kfactor*nbDmin/nb) )  +
1439                                          ( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))  ) +
1440                                          ( (kfactor*fGlobalEfficiencyUncertainties[1])*(kfactor*fGlobalEfficiencyUncertainties[1]) )
1441                                          ) / fhRECpt->GetBinWidth(ibin);
1442
1443       // Correction systematics (fc)
1444       // min value with the maximum Nb
1445       correctionMin = TMath::Sqrt( ( (kfactor*errfrac/frac)*(kfactor*errfrac/frac) ) +
1446                                    ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) +
1447                                    ( (kfactor*nbDmax/nb)*(kfactor*nbDmax/nb) )  +
1448                                    ( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) ) +
1449                                    ( (kfactor*fGlobalEfficiencyUncertainties[1])*(kfactor*fGlobalEfficiencyUncertainties[1]) )
1450                                    ) / fhRECpt->GetBinContent(ibin) ;
1451       // max value with the minimum Nb
1452       correctionMax =  TMath::Sqrt( ( (kfactor*errfrac/frac)*(kfactor*errfrac/frac) ) +
1453                                     ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) +
1454                                     ( (kfactor*nbDmin/nb)*(kfactor*nbDmin/nb) )  +
1455                                     ( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))       ) +
1456                                     ( (kfactor*fGlobalEfficiencyUncertainties[1])*(kfactor*fGlobalEfficiencyUncertainties[1]) )
1457                                     ) / fhRECpt->GetBinContent(ibin) ;
1458       //
1459       correctionUncStatEffb = TMath::Sqrt(  ( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))       )
1460                                             ) / fhRECpt->GetBinContent(ibin) ;
1461       correctionUncStatEffc = 0.;
1462     }
1463     else{ // Don't consider Nb uncertainty in this case [ to be tested!!! ]
1464       errvalueExtremeMax =  TMath::Sqrt( ( (kfactor*errfrac/frac)*(kfactor*errfrac/frac) ) +
1465                                          ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) )  +
1466                                          ( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))  )  +
1467                                          ( (kfactor*fGlobalEfficiencyUncertainties[1])*(kfactor*fGlobalEfficiencyUncertainties[1]) )
1468                                          ) / fhRECpt->GetBinWidth(ibin);
1469       errvalueExtremeMin =  errvalueExtremeMax ;
1470     }
1471
1472
1473     // fill in histograms
1474     fhYieldCorr->SetBinContent(ibin,value);
1475     fhYieldCorr->SetBinError(ibin,errvalue);
1476     //
1477     // Estimate how the result varies vs charm/beauty Eloss hypothesis
1478     //
1479     if ( correction>0.0001 && fPbPbElossHypothesis){
1480       // Loop over the Eloss hypothesis
1481       //      Int_t rbin=0;
1482       for (Float_t rval=0.0025; rval<4.0; rval+=0.005){
1483         // correction =  [ 1  - (Tab *Nevt * delta_y * BR_b * eff_trig * eff_b * Nb_th *binwidth )* (rval) /reco ]
1484         Double_t fcRcbvalue = 1 - (fTab[0]*fNevts*deltaY*branchingRatioBintoFinalDecay*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin)*fhRECpt->GetBinWidth(ibin) * rval ) / fhRECpt->GetBinContent(ibin) ;
1485         if(fcRcbvalue<0.) fcRcbvalue=0.;
1486         fhFcRcb->Fill( fhRECpt->GetBinCenter(ibin) , rval, fcRcbvalue );
1487         //    physics = reco * fcRcb / bin-width
1488         Double_t Rcbvalue = (fhRECpt->GetBinContent(ibin) && fcRcbvalue) ?
1489           fhRECpt->GetBinContent(ibin) * fcRcbvalue : 0. ;
1490         Rcbvalue /= fhRECpt->GetBinWidth(ibin) ;
1491         fhYieldCorrRcb->Fill( fhYieldCorr->GetBinCenter(ibin) , rval, Rcbvalue );
1492         //      if(ibin==3){
1493         //        cout << " pt "<< fhFcRcb->GetBinCenter(ibin) <<" bin "<< ibin<<" rval="<<rval<<", rbin="<<rbin<<", fc-Rb-value="<< fcRcbvalue << ", yield-Rb-value="<< Rcbvalue <<endl;
1494         //      cout << " pt "<< fhFcRcb->GetBinCenter(ibin) <<" bin "<< ibin<<" rval="<<rval<<", fc-Rb-value="<< fcRcbvalue << ", yield-Rb-value="<< Rcbvalue <<endl;
1495         //        rbin++;
1496         //      }
1497       }
1498     }
1499     //
1500     // Fill the rest of (asymmetric) histograms
1501     //
1502     if (fAsymUncertainties) {
1503       Double_t x = fhYieldCorr->GetBinCenter(ibin);
1504       fgYieldCorr->SetPoint(ibin,x,value); // i,x,y
1505       fgYieldCorr->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueMin,errvalueMax); // i,xl,xh,yl,yh
1506       fhYieldCorrMax->SetBinContent(ibin,value+errvalueMax);
1507       fhYieldCorrMin->SetBinContent(ibin,value-errvalueMin);
1508       fgYieldCorrExtreme->SetPoint(ibin,x,value); // i,x,y
1509       fgYieldCorrExtreme->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueExtremeMin,errvalueExtremeMax); // i,xl,xh,yl,yh
1510       fgYieldCorrConservative->SetPoint(ibin,x,value); // i,x,y
1511       fgYieldCorrConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueExtremeMin,errvalueExtremeMax); // i,xl,xh,yl,yh
1512       //      cout << " bin " << ibin << ", correction " << correction << ", min correction unc " << correctionMin << ", max correction unc " << correctionMax << endl;
1513       if(correction>0.){
1514         fgFcConservative->SetPoint(ibin,x,correction);
1515         fgFcConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),correctionMin,correctionMax);
1516
1517         fhStatUncEffbFD->SetBinContent(ibin,0.); fhStatUncEffbFD->SetBinError(ibin,correctionUncStatEffb/correction*100.);
1518         fhStatUncEffcFD->SetBinContent(ibin,0.); fhStatUncEffcFD->SetBinError(ibin,correctionUncStatEffc/correction*100.);
1519         //      cout << " pt "<< fhStatUncEffcFD->GetBinCenter(ibin) <<" bin "<< ibin<<" fc-stat-c ="<< correctionUncStatEffc/correction <<" fc-stat-b ="<< correctionUncStatEffb/correction <<endl;
1520       }
1521       else{
1522         fgFcConservative->SetPoint(ibin,x,0.);
1523         fgFcConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),0.,0.);
1524       }
1525     }
1526
1527   }
1528   delete [] binwidths;
1529   delete [] limits;
1530
1531 }
1532
1533
1534 //_________________________________________________________________________________________________________
1535 void AliHFPtSpectrum::ComputeSystUncertainties(AliHFSystErr *systematics, Bool_t combineFeedDown) {
1536   //
1537   // Function that re-calculates the global systematic uncertainties
1538   //   by calling the class AliHFSystErr and combining those
1539   //   (in quadrature) with the feed-down subtraction uncertainties
1540   //
1541
1542   // Estimate the feed-down uncertainty in percentage
1543   Int_t nentries = fgSigmaCorrConservative->GetN();
1544   TGraphAsymmErrors *grErrFeeddown = new TGraphAsymmErrors(nentries);
1545   Double_t x=0., y=0., errx=0., erryl=0., erryh=0;
1546   for(Int_t i=0; i<nentries; i++) {
1547     x=0.; y=0.; errx=0.; erryl=0.; erryh=0.;
1548     fgSigmaCorrConservative->GetPoint(i,x,y);
1549     if(y>0.){
1550       errx = fgSigmaCorrConservative->GetErrorXlow(i) ;
1551       erryl = fgSigmaCorrConservative->GetErrorYlow(i) / y ;
1552       erryh = fgSigmaCorrConservative->GetErrorYhigh(i) / y ;
1553     }
1554     //    cout << " x "<< x << " +- "<<errx<<" , y "<<y<<" + "<<erryh<<" - "<<erryl<<endl;
1555     grErrFeeddown->SetPoint(i,x,0.);
1556     grErrFeeddown->SetPointError(i,errx,errx,erryl,erryh); //i, xl, xh, yl, yh
1557   }
1558
1559   // Draw all the systematics independently
1560   systematics->DrawErrors(grErrFeeddown);
1561
1562   // Set the sigma systematic uncertainties
1563   // possibly combine with the feed-down uncertainties
1564   Double_t errylcomb=0., erryhcomb=0;
1565   for(Int_t i=1; i<nentries; i++) {
1566     fgSigmaCorr->GetPoint(i,x,y);
1567     errx = grErrFeeddown->GetErrorXlow(i) ;
1568     erryl = grErrFeeddown->GetErrorYlow(i);
1569     erryh = grErrFeeddown->GetErrorYhigh(i);
1570     if (combineFeedDown) {
1571       errylcomb = systematics->GetTotalSystErr(x,erryl) * y ;
1572       erryhcomb = systematics->GetTotalSystErr(x,erryh) * y ;
1573     } else {
1574       errylcomb = systematics->GetTotalSystErr(x) * y ;
1575       erryhcomb = systematics->GetTotalSystErr(x) * y ;
1576     }
1577     fgSigmaCorr->SetPointError(i,errx,errx,errylcomb,erryhcomb);
1578     //
1579     fhSigmaCorrDataSyst->SetBinContent(i,y);
1580     erryl = systematics->GetTotalSystErr(x) * y ;
1581     fhSigmaCorrDataSyst->SetBinError(i,erryl);
1582   }
1583
1584 }
1585
1586
1587 //_________________________________________________________________________________________________________
1588 void AliHFPtSpectrum::DrawSpectrum(TGraphAsymmErrors *gPrediction) {
1589   //
1590   // Example method to draw the corrected spectrum & the theoretical prediction
1591   //
1592
1593   TCanvas *csigma = new TCanvas("csigma","Draw the corrected cross-section & the prediction");
1594   csigma->SetFillColor(0);
1595   gPrediction->GetXaxis()->SetTitleSize(0.05);
1596   gPrediction->GetXaxis()->SetTitleOffset(0.95);
1597   gPrediction->GetYaxis()->SetTitleSize(0.05);
1598   gPrediction->GetYaxis()->SetTitleOffset(0.95);
1599   gPrediction->GetXaxis()->SetTitle("p_{T}  [GeV]");
1600   gPrediction->GetYaxis()->SetTitle("BR #times #frac{d#sigma}{dp_{T}} |_{|y|<0.5}   [pb/GeV]");
1601   gPrediction->SetLineColor(kGreen+2);
1602   gPrediction->SetLineWidth(3);
1603   gPrediction->SetFillColor(kGreen+1);
1604   gPrediction->Draw("3CA");
1605   fgSigmaCorr->SetLineColor(kRed);
1606   fgSigmaCorr->SetLineWidth(1);
1607   fgSigmaCorr->SetFillColor(kRed);
1608   fgSigmaCorr->SetFillStyle(0);
1609   fgSigmaCorr->Draw("2");
1610   fhSigmaCorr->SetMarkerColor(kRed);
1611   fhSigmaCorr->Draw("esame");
1612   csigma->SetLogy();
1613   TLegend * leg = new TLegend(0.7,0.75,0.87,0.5);
1614   leg->SetBorderSize(0);
1615   leg->SetLineColor(0);
1616   leg->SetFillColor(0);
1617   leg->SetTextFont(42);
1618   leg->AddEntry(gPrediction,"FONLL ","fl");
1619   leg->AddEntry(fhSigmaCorr,"data stat. unc.","pl");
1620   leg->AddEntry(fgSigmaCorr,"data syst. unc.","f");
1621   leg->Draw();
1622   csigma->Draw();
1623
1624 }
1625
1626 //_________________________________________________________________________________________________________
1627 TH1D * AliHFPtSpectrum::ReweightHisto(TH1D *hToReweight, TH1D *hReference){
1628   //
1629   // Function to  reweight histograms for testing purposes:
1630   // This function takes the histo hToReweight and reweights
1631   //  it (its pt shape) with respect to hReference
1632   //
1633
1634   // check histograms consistency
1635   Bool_t areconsistent=kTRUE;
1636   areconsistent &= CheckHistosConsistency(hToReweight,hReference);
1637   if (!areconsistent) {
1638     AliInfo("the histograms to reweight are not consistent (bin width, bounds)");
1639     return NULL;
1640   }
1641
1642   // define a new empty histogram
1643   TH1D *hReweighted = (TH1D*)hToReweight->Clone("hReweighted");
1644   hReweighted->Reset();
1645   Double_t weight=1.0;
1646   Double_t yvalue=1.0;
1647   Double_t integralRef = hReference->Integral();
1648   Double_t integralH = hToReweight->Integral();
1649
1650   // now reweight the spectra
1651   //
1652   // the weight is the relative probability of the given pt bin in the reference histo
1653   //  divided by its relative probability (to normalize it) on the histo to re-weight
1654   for (Int_t i=0; i<=hToReweight->GetNbinsX(); i++) {
1655     weight = (hReference->GetBinContent(i)/integralRef) / (hToReweight->GetBinContent(i)/integralH) ;
1656     yvalue = hToReweight->GetBinContent(i);
1657     hReweighted->SetBinContent(i,yvalue*weight);
1658   }
1659
1660   return (TH1D*)hReweighted;
1661 }
1662
1663 //_________________________________________________________________________________________________________
1664 TH1D * AliHFPtSpectrum::ReweightRecHisto(TH1D *hRecToReweight, TH1D *hMCToReweight, TH1D *hMCReference){
1665   //
1666   // Function to  reweight histograms for testing purposes:
1667   // This function takes the histo hToReweight and reweights
1668   //  it (its pt shape) with respect to hReference /hMCToReweight
1669   //
1670
1671   // check histograms consistency
1672   Bool_t areconsistent=kTRUE;
1673   areconsistent &= CheckHistosConsistency(hMCToReweight,hMCReference);
1674   areconsistent &= CheckHistosConsistency(hRecToReweight,hMCReference);
1675   if (!areconsistent) {
1676     AliInfo("the histograms to reweight are not consistent (bin width, bounds)");
1677     return NULL;
1678   }
1679
1680   // define a new empty histogram
1681   TH1D *hReweighted = (TH1D*)hMCToReweight->Clone("hReweighted");
1682   hReweighted->Reset();
1683   TH1D *hRecReweighted = (TH1D*)hRecToReweight->Clone("hRecReweighted");
1684   hRecReweighted->Reset();
1685   Double_t weight=1.0;
1686   Double_t yvalue=1.0, yrecvalue=1.0;
1687   Double_t integralRef = hMCReference->Integral();
1688   Double_t integralH = hMCToReweight->Integral();
1689
1690   // now reweight the spectra
1691   //
1692   // the weight is the relative probability of the given pt bin
1693   //  that should be applied in the MC histo to get the reference histo shape
1694   //  Probabilities are properly normalized.
1695   for (Int_t i=0; i<=hMCToReweight->GetNbinsX(); i++) {
1696     weight = (hMCReference->GetBinContent(i)/integralRef) / (hMCToReweight->GetBinContent(i)/integralH) ;
1697     yvalue = hMCToReweight->GetBinContent(i);
1698     hReweighted->SetBinContent(i,yvalue*weight);
1699     yrecvalue = hRecToReweight->GetBinContent(i);
1700     hRecReweighted->SetBinContent(i,yrecvalue*weight);
1701   }
1702
1703   return (TH1D*)hRecReweighted;
1704 }
1705
1706
1707
1708 //_________________________________________________________________________________________________________
1709 Int_t AliHFPtSpectrum::FindTH2YBin(TH2D *histo, Float_t yvalue){
1710   //
1711   // Function to find the y-axis bin of a TH2 for a given y-value
1712   //
1713
1714   Int_t nbins = histo->GetNbinsY();
1715   Int_t ybin=0;
1716   for (int j=0; j<=nbins; j++) {
1717     Float_t value = histo->GetYaxis()->GetBinCenter(j);
1718     Float_t width = histo->GetYaxis()->GetBinWidth(j);
1719     //    if( TMath::Abs(yvalue-value)<= width/2. ) {
1720     if( TMath::Abs(yvalue-value)<= width ) {
1721       ybin =j;
1722       //      cout <<" value "<<value << ", yval "<< yvalue<<", bin width "<<width/2.<< " y ="<<ybin<<endl;
1723       break;
1724     }
1725   }
1726
1727   return ybin;
1728 }