PWG3/vertexingHF/AliHFPtSpectrum.cxx

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