PWG3/vertexingHF/AliHFPtSpectrum.cxx

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