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