PWGHF/vertexingHF/AliHFPtSpectrum.cxx

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