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