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     delete [] limits;
 565     delete [] binwidths;
 566     return ;
 567   }
 568
 569   Double_t value=0, errValue=0, errvalueMax=0., errvalueMin=0.;
 570   Double_t errvalueExtremeMax=0., errvalueExtremeMin=0.;
 571   Double_t errvalueConservativeMax=0., errvalueConservativeMin=0.;
 572   for(Int_t ibin=1; ibin<=nbins; ibin++){
 573
 574     // Sigma calculation
 575     //   Sigma = ( 1. / (lumi * delta_y * BR_c * eff_trig * eff_c ) ) * spectra (corrected for feed-down)
 576     value = (fhDirectEffpt->GetBinContent(ibin) && fhDirectEffpt->GetBinContent(ibin)!=0.) ?
 577       ( fhYieldCorr->GetBinContent(ibin) / ( deltaY * branchingRatioC * fLuminosity[0] * fTrigEfficiency[0] * fhDirectEffpt->GetBinContent(ibin) ) )
 578       : 0. ;
 579
 580     // Sigma statistical uncertainty:
 581     //   delta_sigma = sigma * sqrt ( (delta_spectra/spectra)^2 )
 582     errValue = (value!=0.) ?  value * (fhYieldCorr->GetBinError(ibin)/fhYieldCorr->GetBinContent(ibin))  : 0. ;
 583
 584
 585     //
 586     // Sigma systematic uncertainties
 587     //
 588     if (fAsymUncertainties) {
 589
 590       //  (syst but feed-down) delta_sigma = sigma * sqrt ( (delta_spectra_syst/spectra)^2 +
 591       //                                     (delta_lumi/lumi)^2 + (delta_eff_trig/eff_trig)^2 + (delta_eff/eff)^2  + (global_eff)^2 )
 592       errvalueMax = value * TMath::Sqrt( (fgYieldCorr->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin))*(fgYieldCorr->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin)) +
 593                                          (fLuminosity[1]/fLuminosity[0])*(fLuminosity[1]/fLuminosity[0]) +
 594                                          (fTrigEfficiency[1]/fTrigEfficiency[0])*(fTrigEfficiency[1]/fTrigEfficiency[0])  +
 595                                          (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin)) +
 596                                          fGlobalEfficiencyUncertainties[0]*fGlobalEfficiencyUncertainties[0] );
 597       errvalueMin = value * TMath::Sqrt( (fgYieldCorr->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin))*(fgYieldCorr->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin)) +
 598                                          (fLuminosity[1]/fLuminosity[0])*(fLuminosity[1]/fLuminosity[0]) +
 599                                          (fTrigEfficiency[1]/fTrigEfficiency[0])*(fTrigEfficiency[1]/fTrigEfficiency[0])  +
 600                                          (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))  +
 601                                          fGlobalEfficiencyUncertainties[0]*fGlobalEfficiencyUncertainties[0] );
 602
 603       // Uncertainties from feed-down
 604       //      (feed-down syst) delta_sigma = sigma * sqrt ( (delta_spectra_fd/spectra_fd)^2 )
 605       //   extreme case
 606       errvalueExtremeMax = value * (fgYieldCorrExtreme->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin));
 607       errvalueExtremeMin =  value * (fgYieldCorrExtreme->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin));
 608       //
 609       //   conservative case
 610       errvalueConservativeMax = value * (fgYieldCorrConservative->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin));
 611       errvalueConservativeMin =  value * (fgYieldCorrConservative->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin));
 612
 613     }
 614     else {
 615       // protect against null denominator
 616       errvalueMax = (value!=0.) ?
 617         value * TMath::Sqrt( (fLuminosity[1]/fLuminosity[0])*(fLuminosity[1]/fLuminosity[0]) +
 618                              (fTrigEfficiency[1]/fTrigEfficiency[0])*(fTrigEfficiency[1]/fTrigEfficiency[0])  +
 619                              (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))  +
 620                              fGlobalEfficiencyUncertainties[0]*fGlobalEfficiencyUncertainties[0] )
 621         : 0. ;
 622       errvalueMin = errvalueMax;
 623     }
 624
 625     // Fill the histograms
 626     fhSigmaCorr->SetBinContent(ibin,value);
 627     fhSigmaCorr->SetBinError(ibin,errValue);
 628     // Fill the TGraphAsymmErrors
 629     if (fAsymUncertainties) {
 630       Double_t x = fhYieldCorr->GetBinCenter(ibin);
 631       fgSigmaCorr->SetPoint(ibin,x,value); // i,x,y
 632       fgSigmaCorr->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueMin,errvalueMax); // i,xl,xh,yl,yh
 633       fhSigmaCorrMax->SetBinContent(ibin,value+errvalueMax);
 634       fhSigmaCorrMin->SetBinContent(ibin,value-errvalueMin);
 635       fgSigmaCorrExtreme->SetPoint(ibin,x,value); // i,x,y
 636       fgSigmaCorrExtreme->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueExtremeMin,errvalueExtremeMax); // i,xl,xh,yl,yh
 637       fgSigmaCorrConservative->SetPoint(ibin,x,value); // i,x,y
 638       fgSigmaCorrConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueConservativeMin,errvalueConservativeMax); // i,xl,xh,yl,yh
 639
 640     }
 641
 642   }
 643   delete [] binwidths;
 644   delete [] limits;
 645
 646 }
 647
 648 //_________________________________________________________________________________________________________
 649 TH1D * AliHFPtSpectrum::EstimateEfficiencyRecoBin(TH1D *hSimu, TH1D *hReco, const char *name) {
 650   //
 651   // Function that computes the acceptance and efficiency correction
 652   //  based on the simulated and reconstructed spectra
 653   //  and using the reconstructed spectra bin width
 654   //
 655   //  eff = reco/sim ; err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
 656   //
 657
 658   if(!fhRECpt){
 659     AliInfo("Hey, the reconstructed histogram was not set yet !");
 660     return NULL;
 661   }
 662
 663   Int_t nbins = fhRECpt->GetNbinsX();
 664   Double_t *limits = new Double_t[nbins+1];
 665   Double_t xlow=0.,binwidth=0.;
 666   for (Int_t i=1; i<=nbins; i++) {
 667     binwidth = fhRECpt->GetBinWidth(i);
 668     xlow = fhRECpt->GetBinLowEdge(i);
 669     limits[i-1] = xlow;
 670   }
 671   limits[nbins] = xlow + binwidth;
 672
 673   TH1D * hEfficiency = new TH1D(name," acceptance #times efficiency",nbins,limits);
 674
 675   Double_t sumSimu[nbins], sumReco[nbins];
 676   for (Int_t ibin=0; ibin<nbins; ibin++){
 677     sumSimu[ibin]=0.;  sumReco[ibin]=0.;
 678   }
 679   for (Int_t ibin=0; ibin<=hSimu->GetNbinsX(); ibin++){
 680
 681     for (Int_t ibinrec=0; ibinrec<nbins; ibinrec++){
 682       if ( hSimu->GetBinCenter(ibin)>limits[ibinrec] &&
 683            hSimu->GetBinCenter(ibin)<limits[ibinrec+1] ) {
 684         sumSimu[ibinrec]+=hSimu->GetBinContent(ibin);
 685       }
 686       if ( hReco->GetBinCenter(ibin)>limits[ibinrec] &&
 687            hReco->GetBinCenter(ibin)<limits[ibinrec+1] ) {
 688         sumReco[ibinrec]+=hReco->GetBinContent(ibin);
 689       }
 690     }
 691
 692   }
 693
 694
 695   // the efficiency is computed as reco/sim (in each bin)
 696   //  its uncertainty is err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
 697   Double_t eff=0., erreff=0.;
 698   for (Int_t ibinrec=0; ibinrec<=nbins; ibinrec++) {
 699     if (sumSimu[ibinrec]!= 0. && sumReco[ibinrec]!=0.) {
 700       eff = sumReco[ibinrec] / sumSimu[ibinrec] ;
 701       // protection in case eff > 1.0
 702       // test calculation (make the argument of the sqrt positive)
 703       erreff = TMath::Sqrt( eff * TMath::Abs(1.0 - eff) ) / TMath::Sqrt( sumSimu[ibinrec] );
 704     }
 705     else { eff=0.0; erreff=0.; }
 706     hEfficiency->SetBinContent(ibinrec+1,eff);
 707     hEfficiency->SetBinError(ibinrec+1,erreff);
 708   }
 709
 710   return (TH1D*)hEfficiency;
 711 }
 712
 713 //_________________________________________________________________________________________________________
 714 void AliHFPtSpectrum::EstimateAndSetDirectEfficiencyRecoBin(TH1D *hSimu, TH1D *hReco) {
 715   //
 716   // Function that computes the Direct  acceptance and efficiency correction
 717   //  based on the simulated and reconstructed spectra
 718   //  and using the reconstructed spectra bin width
 719   //
 720   //  eff = reco/sim ; err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
 721   //
 722
 723   if(!fhRECpt || !hSimu || !hReco){
 724     AliError("Hey, the reconstructed histogram was not set yet !");
 725     return;
 726   }
 727
 728   fhDirectEffpt = EstimateEfficiencyRecoBin(hSimu,hReco,"fhDirectEffpt");
 729   fhDirectEffpt->SetNameTitle("fhDirectEffpt"," direct acceptance #times efficiency");
 730
 731 }
 732
 733 //_________________________________________________________________________________________________________
 734 void AliHFPtSpectrum::EstimateAndSetFeedDownEfficiencyRecoBin(TH1D *hSimu, TH1D *hReco) {
 735   //
 736   // Function that computes the Feed-Down acceptance and efficiency correction
 737   //  based on the simulated and reconstructed spectra
 738   //  and using the reconstructed spectra bin width
 739   //
 740   //  eff = reco/sim ; err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
 741   //
 742
 743   if(!fhRECpt || !hSimu || !hReco){
 744     AliError("Hey, the reconstructed histogram was not set yet !");
 745     return;
 746   }
 747
 748   fhFeedDownEffpt = EstimateEfficiencyRecoBin(hSimu,hReco,"fhFeedDownEffpt");
 749   fhFeedDownEffpt->SetNameTitle("fhFeedDownEffpt"," feed-down acceptance #times efficiency");
 750
 751 }
 752
 753 //_________________________________________________________________________________________________________
 754 Bool_t AliHFPtSpectrum::Initialize(){
 755   //
 756   // Initialization of the variables (histograms)
 757   //
 758
 759   if (fFeedDownOption==0) {
 760     AliInfo("Getting ready for the corrections without feed-down consideration");
 761   } else if (fFeedDownOption==1) {
 762     AliInfo("Getting ready for the fc feed-down correction calculation");
 763   } else if (fFeedDownOption==2) {
 764     AliInfo("Getting ready for the Nb feed-down correction calculation");
 765   } else { AliError("The calculation option must be <=2"); return kFALSE; }
 766
 767   // Start checking the input histograms consistency
 768   Bool_t areconsistent=kTRUE;
 769
 770   // General checks
 771   if (!fhDirectEffpt || !fhRECpt) {
 772     AliError(" Reconstructed spectra and/or the Nc efficiency distributions are not defined");
 773     return kFALSE;
 774   }
 775   areconsistent &= CheckHistosConsistency(fhRECpt,fhDirectEffpt);
 776   if (!areconsistent) {
 777     AliInfo("Histograms required for Nb correction are not consistent (bin width, bounds)");
 778     return kFALSE;
 779   }
 780   if (fFeedDownOption==0) return kTRUE;
 781
 782   //
 783   // Common checks for options 1 (fc) & 2(Nb)
 784   if (!fhFeedDownMCpt || !fhFeedDownEffpt) {
 785     AliError(" Theoretical Nb and/or the Nb efficiency distributions are not defined");
 786     return kFALSE;
 787   }
 788   areconsistent &= CheckHistosConsistency(fhRECpt,fhFeedDownMCpt);
 789   areconsistent &= CheckHistosConsistency(fhFeedDownMCpt,fhFeedDownEffpt);
 790   if (fAsymUncertainties) {
 791     if (!fhFeedDownMCptMax || !fhFeedDownMCptMin) {
 792       AliError(" Max/Min theoretical Nb distributions are not defined");
 793       return kFALSE;
 794     }
 795     areconsistent &= CheckHistosConsistency(fhFeedDownMCpt,fhFeedDownMCptMax);
 796   }
 797   if (!areconsistent) {
 798     AliInfo("Histograms required for Nb correction are not consistent (bin width, bounds)");
 799     return kFALSE;
 800   }
 801   if (fFeedDownOption>1) return kTRUE;
 802
 803   //
 804   // Now checks for option 1 (fc correction)
 805   if (!fhDirectMCpt) {
 806     AliError("Theoretical Nc distributions is not defined");
 807     return kFALSE;
 808   }
 809   areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhFeedDownMCpt);
 810   areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhDirectEffpt);
 811   if (fAsymUncertainties) {
 812     if (!fhDirectMCptMax || !fhDirectMCptMin) {
 813       AliError(" Max/Min theoretical Nc distributions are not defined");
 814       return kFALSE;
 815     }
 816     areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhDirectMCptMax);
 817   }
 818   if (!areconsistent) {
 819     AliInfo("Histograms required for fc correction are not consistent (bin width, bounds)");
 820     return kFALSE;
 821   }
 822
 823   return kTRUE;
 824 }
 825
 826 //_________________________________________________________________________________________________________
 827 Bool_t AliHFPtSpectrum::CheckHistosConsistency(TH1D *h1, TH1D *h2){
 828   //
 829   // Check the histograms consistency (bins, limits)
 830   //
 831
 832   if (!h1 || !h2) {
 833     AliError("One or both histograms don't exist");
 834     return kFALSE;
 835   }
 836
 837   Double_t binwidth1 = h1->GetBinWidth(1);
 838   Double_t binwidth2 = h2->GetBinWidth(1);
 839   Double_t min1 = h1->GetBinCenter(1) - (binwidth1/2.) ;
 840 //   Double_t max1 = h1->GetBinCenter(nbins1) + (binwidth1/2.) ;
 841   Double_t min2 = h2->GetBinCenter(1) - (binwidth2/2.) ;
 842 //   Double_t max2 = h2->GetBinCenter(nbins2) + (binwidth2/2.) ;
 843
 844   if (binwidth1!=binwidth2) {
 845     AliInfo(" histograms with different bin width");
 846     return kFALSE;
 847   }
 848   if (min1!=min2) {
 849     AliInfo(" histograms with different minimum");
 850     return kFALSE;
 851   }
 852 //   if (max1!=max2) {
 853 //     AliInfo(" histograms with different maximum");
 854 //     return kFALSE;
 855 //   }
 856
 857   return kTRUE;
 858 }
 859
 860 //_________________________________________________________________________________________________________
 861 void AliHFPtSpectrum::CalculateFeedDownCorrectionFc(){
 862   //
 863   // Compute fc factor and its uncertainties bin by bin
 864   //   fc = 1 / ( 1 + (eff_b/eff_c)*(N_b/N_c) )
 865   //
 866   // uncertainties: (conservative) combine the upper/lower N_b & N_c predictions together
 867   //                (extreme) combine the upper N_b predictions with the lower N_c predictions & viceversa
 868   //                systematic uncertainty on the acceptance x efficiency b/c ratio are included
 869   //
 870
 871   // define the variables
 872   Int_t nbins = fhRECpt->GetNbinsX();
 873   Double_t binwidth = fhRECpt->GetBinWidth(1);
 874   Double_t *limits = new Double_t[nbins+1];
 875   Double_t *binwidths = new Double_t[nbins];
 876   Double_t xlow=0.;
 877   for (Int_t i=1; i<=nbins; i++) {
 878     binwidth = fhRECpt->GetBinWidth(i);
 879     xlow = fhRECpt->GetBinLowEdge(i);
 880     limits[i-1] = xlow;
 881     binwidths[i-1] = binwidth;
 882   }
 883   limits[nbins] = xlow + binwidth;
 884
 885   Double_t correction=1.;
 886   Double_t theoryRatio=1.;
 887   Double_t effRatio=1.;
 888   Double_t correctionExtremeA=1., correctionExtremeB=1.;
 889   Double_t theoryRatioExtremeA=1., theoryRatioExtremeB=1.;
 890   Double_t correctionConservativeA=1., correctionConservativeB=1.;
 891   Double_t theoryRatioConservativeA=1., theoryRatioConservativeB=1.;
 892   Double_t correctionUnc=0.;
 893   Double_t correctionExtremeAUnc=0., correctionExtremeBUnc=0.;
 894   Double_t correctionConservativeAUnc=0., correctionConservativeBUnc=0.;
 895
 896   // declare the output histograms
 897   fhFc = new TH1D("fhFc","fc correction factor",nbins,limits);
 898   fhFcMax = new TH1D("fhFcMax","max fc correction factor",nbins,limits);
 899   fhFcMin = new TH1D("fhFcMin","min fc correction factor",nbins,limits);
 900   // two local control histograms
 901   TH1D *hTheoryRatio = new TH1D("hTheoryRatio","Theoretical B-->D over c-->D (feed-down/direct) ratio",nbins,limits);
 902   TH1D *hEffRatio = new TH1D("hEffRatio","Efficiency B-->D over c-->D (feed-down/direct) ratio",nbins,limits);
 903   // and the output TGraphAsymmErrors
 904   if (fAsymUncertainties) {
 905     fgFcExtreme = new TGraphAsymmErrors(nbins+1);
 906     fgFcExtreme->SetNameTitle("fgFcExtreme","fgFcExtreme");
 907     fgFcConservative = new TGraphAsymmErrors(nbins+1);
 908     fgFcConservative->SetNameTitle("fgFcConservative","fgFcConservative");
 909   }
 910
 911
 912   //
 913   // Compute fc
 914   //
 915   for (Int_t ibin=1; ibin<=nbins; ibin++) {
 916
 917     //  theory_ratio = (N_b/N_c)
 918     theoryRatio = (fhDirectMCpt->GetBinContent(ibin)>0. && fhFeedDownMCpt->GetBinContent(ibin)>0.) ?
 919       fhFeedDownMCpt->GetBinContent(ibin) / fhDirectMCpt->GetBinContent(ibin) : 1.0 ;
 920
 921     //
 922     // Calculate the uncertainty [ considering only the theoretical uncertainties on Nb & Nc for now !!! ]
 923     //
 924     // extreme A = direct-max, feed-down-min
 925     theoryRatioExtremeA = (fhDirectMCptMax->GetBinContent(ibin)>0. && fhFeedDownMCptMin->GetBinContent(ibin)>0.) ?
 926       fhFeedDownMCptMin->GetBinContent(ibin) / fhDirectMCptMax->GetBinContent(ibin) : 1.0 ;
 927     // extreme B = direct-min, feed-down-max
 928     theoryRatioExtremeB = (fhDirectMCptMin->GetBinContent(ibin)>0. && fhDirectMCptMax->GetBinContent(ibin)>0.) ?
 929       fhFeedDownMCptMax->GetBinContent(ibin) / fhDirectMCptMin->GetBinContent(ibin) : 1.0 ;
 930     // conservative A = direct-max, feed-down-max
 931     theoryRatioConservativeA = (fhDirectMCptMax->GetBinContent(ibin)>0. && fhFeedDownMCptMin->GetBinContent(ibin)>0.) ?
 932       fhFeedDownMCptMax->GetBinContent(ibin) / fhDirectMCptMax->GetBinContent(ibin) : 1.0 ;
 933     // conservative B = direct-min, feed-down-min
 934     theoryRatioConservativeB = (fhDirectMCptMin->GetBinContent(ibin)>0. && fhDirectMCptMax->GetBinContent(ibin)>0.) ?
 935       fhFeedDownMCptMin->GetBinContent(ibin) / fhDirectMCptMin->GetBinContent(ibin) : 1.0 ;
 936
 937     //  eff_ratio = (eff_b/eff_c)
 938     effRatio = (fhDirectEffpt->GetBinContent(ibin) && fhDirectEffpt->GetBinContent(ibin)!=0.) ?
 939       fhFeedDownEffpt->GetBinContent(ibin) / fhDirectEffpt->GetBinContent(ibin) : 1.0 ;
 940
 941     //   fc = 1 / ( 1 + (eff_b/eff_c)*(N_b/N_c) )
 942     if( TMath::Abs(effRatio - 1.0)<0.01 || TMath::Abs(theoryRatio - 1.0)<0.01 ) {
 943       correction = 1.0;
 944       correctionExtremeA = 1.0;
 945       correctionExtremeB = 1.0;
 946       correctionConservativeA = 1.0;
 947       correctionConservativeB = 1.0;
 948     }
 949     else {
 950       correction = ( 1. / ( 1 + ( effRatio * theoryRatio ) ) );
 951       correctionExtremeA = ( 1. / ( 1 + ( effRatio * theoryRatioExtremeA ) ) );
 952       correctionExtremeB = ( 1. / ( 1 + ( effRatio * theoryRatioExtremeB ) ) );
 953       correctionConservativeA = ( 1. / ( 1 + ( effRatio * theoryRatioConservativeA ) ) );
 954       correctionConservativeB = ( 1. / ( 1 + ( effRatio * theoryRatioConservativeB ) ) );
 955     }
 956
 957
 958     // fc uncertainty from (eff_b/eff_c) = fc^2 * (N_b/N_c) * delta(eff_b/eff_c)
 959     //  delta(eff_b/eff_c) is a percentage = fGlobalEfficiencyUncertainties[1]*effRatio
 960     //    correctionUnc = correction*correction * theoryRatio * fGlobalEfficiencyUncertainties[1]*effRatio;
 961     //     correctionExtremeAUnc = correctionExtremeA*correctionExtremeA * theoryRatioExtremeA  * fGlobalEfficiencyUncertainties[1]*effRatio;
 962     //     correctionExtremeBUnc = correctionExtremeB*correctionExtremeB * theoryRatioExtremeB  * fGlobalEfficiencyUncertainties[1]*effRatio;
 963     //     correctionConservativeAUnc = correctionConservativeA*correctionConservativeA * theoryRatioConservativeA  * fGlobalEfficiencyUncertainties[1]*effRatio;
 964     //     correctionConservativeBUnc = correctionConservativeB*correctionConservativeB * theoryRatioConservativeB  * fGlobalEfficiencyUncertainties[1]*effRatio;
 965     correctionUnc = correction*correction * theoryRatio * effRatio *
 966       TMath::Sqrt( fGlobalEfficiencyUncertainties[1]*fGlobalEfficiencyUncertainties[1] +
 967                    (fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) +
 968                    (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))
 969                    );
 970     correctionExtremeAUnc = correctionExtremeA*correctionExtremeA * theoryRatioExtremeA  * effRatio *
 971       TMath::Sqrt( fGlobalEfficiencyUncertainties[1]*fGlobalEfficiencyUncertainties[1] +
 972                    (fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) +
 973                    (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))
 974                    );
 975     correctionExtremeBUnc = correctionExtremeB*correctionExtremeB * theoryRatioExtremeB  * effRatio *
 976       TMath::Sqrt( fGlobalEfficiencyUncertainties[1]*fGlobalEfficiencyUncertainties[1] +
 977                    (fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) +
 978                    (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))
 979                    );
 980     correctionConservativeAUnc = correctionConservativeA*correctionConservativeA * theoryRatioConservativeA  *effRatio *
 981       TMath::Sqrt( fGlobalEfficiencyUncertainties[1]*fGlobalEfficiencyUncertainties[1] +
 982                    (fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) +
 983                    (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))
 984                    );
 985     correctionConservativeBUnc = correctionConservativeB*correctionConservativeB * theoryRatioConservativeB  *effRatio *
 986       TMath::Sqrt( fGlobalEfficiencyUncertainties[1]*fGlobalEfficiencyUncertainties[1] +
 987                    (fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) +
 988                    (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))
 989                    );
 990
 991
 992     // Fill in the histograms
 993     hTheoryRatio->SetBinContent(ibin,theoryRatio);
 994     hEffRatio->SetBinContent(ibin,effRatio);
 995     fhFc->SetBinContent(ibin,correction);
 996     if (fAsymUncertainties) {
 997       Double_t x = fhDirectMCpt->GetBinCenter(ibin);
 998       Double_t val[4] = { correctionExtremeA + correctionExtremeAUnc, correctionExtremeA - correctionExtremeAUnc,
 999                           correctionExtremeB + correctionExtremeBUnc, correctionExtremeB - correctionExtremeBUnc };
1000       Double_t uncExtremeMin = correction - TMath::MinElement(4,val);
1001       Double_t uncExtremeMax = TMath::MaxElement(4,val) - correction;
1002       fgFcExtreme->SetPoint(ibin,x,correction); // i,x,y
1003       fgFcExtreme->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),uncExtremeMin,uncExtremeMax); // i,xl,xh,yl,yh
1004       fhFcMax->SetBinContent(ibin,correction+uncExtremeMax);
1005       fhFcMin->SetBinContent(ibin,correction-uncExtremeMin);
1006       Double_t consval[4] = { correctionConservativeA - correctionConservativeAUnc, correctionConservativeA + correctionConservativeAUnc,
1007                               correctionConservativeB - correctionConservativeBUnc, correctionConservativeB + correctionConservativeBUnc};
1008       Double_t uncConservativeMin = correction - TMath::MinElement(4,consval);
1009       Double_t uncConservativeMax =  TMath::MaxElement(4,consval) - correction;
1010       fgFcConservative->SetPoint(ibin,x,correction); // i,x,y
1011       fgFcConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),uncConservativeMin,uncConservativeMax); // i,xl,xh,yl,yh
1012       if( !(correction>0.) ){
1013         fgFcExtreme->SetPoint(ibin,x,0.); // i,x,y
1014         fgFcExtreme->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),0.,0.); // i,xl,xh,yl,yh
1015         fgFcConservative->SetPoint(ibin,x,0.); // i,x,y
1016         fgFcConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),0.,0.); // i,xl,xh,yl,yh
1017       }
1018     }
1019
1020   }
1021   delete [] binwidths;
1022   delete [] limits;
1023
1024 }
1025
1026 //_________________________________________________________________________________________________________
1027 void AliHFPtSpectrum::CalculateFeedDownCorrectedSpectrumFc(){
1028   //
1029   // Compute the feed-down corrected spectrum if feed-down correction is done via fc factor (bin by bin)
1030   //    physics = reco * fc / bin-width
1031   //
1032   //    uncertainty:             (stat) delta_physics = physics * sqrt ( (delta_reco/reco)^2 )
1033   //               (syst but feed-down) delta_physics = physics * sqrt ( (delta_reco_syst/reco)^2 )
1034   //                   (feed-down syst) delta_physics = physics * sqrt ( (delta_fc/fc)^2 )
1035   //
1036   //    ( Calculation done bin by bin )
1037
1038   if (!fhFc || !fhRECpt) {
1039     AliError(" Reconstructed or fc distributions are not defined");
1040     return;
1041   }
1042
1043   Int_t nbins = fhRECpt->GetNbinsX();
1044   Double_t value = 0., errvalue = 0., errvalueMax= 0., errvalueMin= 0.;
1045   Double_t valueExtremeMax= 0., valueExtremeMin= 0.;
1046   Double_t valueConservativeMax= 0., valueConservativeMin= 0.;
1047   Double_t binwidth = fhRECpt->GetBinWidth(1);
1048   Double_t *limits = new Double_t[nbins+1];
1049   Double_t *binwidths = new Double_t[nbins];
1050   Double_t xlow=0.;
1051   for (Int_t i=1; i<=nbins; i++) {
1052     binwidth = fhRECpt->GetBinWidth(i);
1053     xlow = fhRECpt->GetBinLowEdge(i);
1054     limits[i-1] = xlow;
1055     binwidths[i-1] = binwidth;
1056   }
1057   limits[nbins] = xlow + binwidth;
1058
1059   // declare the output histograms
1060   fhYieldCorr = new TH1D("fhYieldCorr","corrected yield (by fc)",nbins,limits);
1061   fhYieldCorrMax = new TH1D("fhYieldCorrMax","max corrected yield (by fc)",nbins,limits);
1062   fhYieldCorrMin = new TH1D("fhYieldCorrMin","min corrected yield (by fc)",nbins,limits);
1063   // and the output TGraphAsymmErrors
1064   if (fAsymUncertainties){
1065     fgYieldCorr = new TGraphAsymmErrors(nbins+1);
1066     fgYieldCorrExtreme = new TGraphAsymmErrors(nbins+1);
1067     fgYieldCorrConservative = new TGraphAsymmErrors(nbins+1);
1068   }
1069
1070   //
1071   // Do the calculation
1072   //
1073   for (Int_t ibin=1; ibin<=nbins; ibin++) {
1074
1075     // calculate the value
1076     //    physics = reco * fc / bin-width
1077     value = (fhRECpt->GetBinContent(ibin) && fhFc->GetBinContent(ibin)) ?
1078       fhRECpt->GetBinContent(ibin) * fhFc->GetBinContent(ibin) : 0. ;
1079     value /= fhRECpt->GetBinWidth(ibin) ;
1080
1081     // Statistical uncertainty
1082     //    (stat) delta_physics = physics * sqrt ( (delta_reco/reco)^2 )
1083     errvalue = (value!=0. && fhRECpt->GetBinContent(ibin) && fhRECpt->GetBinContent(ibin)!=0.) ?
1084       value * (fhRECpt->GetBinError(ibin)/fhRECpt->GetBinContent(ibin)) : 0. ;
1085
1086     // Calculate the systematic uncertainties
1087     //    (syst but feed-down) delta_physics = physics * sqrt ( (delta_reco_syst/reco)^2 )
1088     //        (feed-down syst) delta_physics = physics * sqrt ( (delta_fc/fc)^2 )
1089     //
1090     // Protect against null denominator. If so, define uncertainty as null
1091     if (fhRECpt->GetBinContent(ibin) && fhRECpt->GetBinContent(ibin)!=0.) {
1092
1093       if (fAsymUncertainties) {
1094
1095         // Systematics but feed-down
1096         if (fgRECSystematics) {
1097           errvalueMax = value * ( fgRECSystematics->GetErrorYhigh(ibin) / fhRECpt->GetBinContent(ibin) );
1098           errvalueMin = value * ( fgRECSystematics->GetErrorYlow(ibin) / fhRECpt->GetBinContent(ibin) );
1099         }
1100         else { errvalueMax = 0.; errvalueMin = 0.; }
1101
1102         // Extreme feed-down systematics
1103         valueExtremeMax = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) + fgFcExtreme->GetErrorYhigh(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
1104         valueExtremeMin = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) - fgFcExtreme->GetErrorYlow(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
1105
1106         // Conservative feed-down systematics
1107         valueConservativeMax = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) + fgFcConservative->GetErrorYhigh(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
1108         valueConservativeMin = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) - fgFcConservative->GetErrorYlow(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
1109
1110       }
1111
1112     }
1113     else { errvalueMax = 0.; errvalueMin = 0.; }
1114
1115     // fill in the histograms
1116     fhYieldCorr->SetBinContent(ibin,value);
1117     fhYieldCorr->SetBinError(ibin,errvalue);
1118     if (fAsymUncertainties) {
1119       Double_t center = fhYieldCorr->GetBinCenter(ibin);
1120       fgYieldCorr->SetPoint(ibin,center,value); // i,x,y
1121       fgYieldCorr->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueMin,errvalueMax); // i,xl,xh,yl,yh
1122       fhYieldCorrMax->SetBinContent(ibin,value+errvalueMax);
1123       fhYieldCorrMin->SetBinContent(ibin,value-errvalueMin);
1124       fgYieldCorrExtreme->SetPoint(ibin,center,value);
1125       fgYieldCorrExtreme->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),value-valueExtremeMin,valueExtremeMax-value);
1126       fgYieldCorrConservative->SetPoint(ibin,center,value);
1127       fgYieldCorrConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),value-valueConservativeMin,valueConservativeMax-value);
1128     }
1129
1130   }
1131   delete [] binwidths;
1132   delete [] limits;
1133
1134 }
1135
1136 //_________________________________________________________________________________________________________
1137 void AliHFPtSpectrum::CalculateFeedDownCorrectedSpectrumNb(Double_t deltaY, Double_t branchingRatioBintoFinalDecay) {
1138   //
1139   // Compute the feed-down corrected spectrum if feed-down correction is done via Nb (bin by bin)
1140   //    physics =  [ reco  - (lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th) ] / bin-width
1141   //
1142   //    uncertainty:   (stat)  delta_physics = sqrt ( (delta_reco)^2 )  / bin-width
1143   //     (syst but feed-down)  delta_physics = sqrt ( (delta_reco_syst)^2 )  / bin-width
1144   //         (feed-down syst)  delta_physics = sqrt ( (k*delta_lumi/lumi)^2 + (k*delta_eff_trig/eff_trig)^2
1145   //                                                   + (k*delta_Nb/Nb)^2 + (k*delta_eff/eff)^2  + (k*global_eff_ratio)^2 ) / bin-width
1146   //                    where k = lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th
1147   //
1148
1149   Int_t nbins = fhRECpt->GetNbinsX();
1150   Double_t binwidth = fhRECpt->GetBinWidth(1);
1151   Double_t value = 0., errvalue = 0., errvalueMax = 0., errvalueMin = 0., kfactor = 0.;
1152   Double_t errvalueExtremeMax = 0., errvalueExtremeMin = 0.;
1153   Double_t *limits = new Double_t[nbins+1];
1154   Double_t *binwidths = new Double_t[nbins];
1155   Double_t xlow=0.;
1156   for (Int_t i=1; i<=nbins; i++) {
1157     binwidth = fhRECpt->GetBinWidth(i);
1158     xlow = fhRECpt->GetBinLowEdge(i);
1159     limits[i-1] = xlow;
1160     binwidths[i-1] = binwidth;
1161   }
1162   limits[nbins] = xlow + binwidth;
1163
1164   // declare the output histograms
1165   fhYieldCorr = new TH1D("fhYieldCorr","corrected yield (by Nb)",nbins,limits);
1166   fhYieldCorrMax = new TH1D("fhYieldCorrMax","max corrected yield (by Nb)",nbins,limits);
1167   fhYieldCorrMin = new TH1D("fhYieldCorrMin","min corrected yield (by Nb)",nbins,limits);
1168   // and the output TGraphAsymmErrors
1169   if (fAsymUncertainties){
1170     fgYieldCorr = new TGraphAsymmErrors(nbins+1);
1171     fgYieldCorrExtreme = new TGraphAsymmErrors(nbins+1);
1172     fgYieldCorrConservative = new TGraphAsymmErrors(nbins+1);
1173     // Define fc-conservative
1174     fgFcConservative = new TGraphAsymmErrors(nbins+1);
1175     AliInfo(" Beware the conservative & extreme uncertainties are equal by definition !");
1176   }
1177
1178   // variables to define fc-conservative
1179   double correction=0, correctionMax=0., correctionMin=0.;
1180
1181   //
1182   // Do the calculation
1183   //
1184   for (Int_t ibin=1; ibin<=nbins; ibin++) {
1185
1186     // Calculate the value
1187     //    physics =  [ reco  - (lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th) ] / bin-width
1188     value = ( fhRECpt->GetBinContent(ibin) && fhRECpt->GetBinContent(ibin)!=0. &&
1189               fhFeedDownMCpt->GetBinContent(ibin)>0. && fhFeedDownEffpt->GetBinContent(ibin)>0. ) ?
1190       fhRECpt->GetBinContent(ibin) - (deltaY*branchingRatioBintoFinalDecay*fLuminosity[0]*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) )
1191       : 0. ;
1192     value /= fhRECpt->GetBinWidth(ibin);
1193
1194     //  Statistical uncertainty:   delta_physics = sqrt ( (delta_reco)^2 )  / bin-width
1195     errvalue = (value!=0. && fhRECpt->GetBinError(ibin) && fhRECpt->GetBinError(ibin)!=0.)  ?
1196       fhRECpt->GetBinError(ibin) : 0.;
1197     errvalue /= fhRECpt->GetBinWidth(ibin);
1198
1199     // Correction (fc) : Estimate of the relative amount feed-down subtracted
1200     // correction =  [ 1  - (lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th)/reco ]
1201     correction = 1 - (deltaY*branchingRatioBintoFinalDecay*fLuminosity[0]*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) ) / fhRECpt->GetBinContent(ibin) ;
1202
1203     // Systematic uncertainties
1204     //     (syst but feed-down)  delta_physics = sqrt ( (delta_reco_syst)^2 )  / bin-width
1205     //         (feed-down syst)  delta_physics = sqrt ( (k*delta_lumi/lumi)^2 + (k*delta_eff_trig/eff_trig)^2
1206     //                                                   + (k*delta_Nb/Nb)^2 + (k*delta_eff/eff)^2 + (k*global_eff_ratio)^2 ) / bin-width
1207     //                    where k = lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th
1208     kfactor = deltaY*branchingRatioBintoFinalDecay*fLuminosity[0]*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) ;
1209
1210     //
1211     if (fAsymUncertainties) {
1212       Double_t nb =  fhFeedDownMCpt->GetBinContent(ibin);
1213       Double_t nbDmax = fhFeedDownMCptMax->GetBinContent(ibin) - fhFeedDownMCpt->GetBinContent(ibin);
1214       Double_t nbDmin = fhFeedDownMCpt->GetBinContent(ibin) - fhFeedDownMCptMin->GetBinContent(ibin);
1215
1216       // Systematics but feed-down
1217       if (fgRECSystematics){
1218         errvalueMax = fgRECSystematics->GetErrorYhigh(ibin) / fhRECpt->GetBinWidth(ibin) ;
1219         errvalueMin = fgRECSystematics->GetErrorYlow(ibin) / fhRECpt->GetBinWidth(ibin);
1220       }
1221       else { errvalueMax = 0.; errvalueMin = 0.; }
1222
1223       // Feed-down systematics
1224       // min value with the maximum Nb
1225       errvalueExtremeMin = TMath::Sqrt( ( (kfactor*fLuminosity[1]/fLuminosity[0])*(kfactor*fLuminosity[1]/fLuminosity[0]) ) +
1226                                         ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) +
1227                                         ( (kfactor*nbDmax/nb)*(kfactor*nbDmax/nb) )  +
1228                                         ( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) ) +
1229                                         ( (kfactor*fGlobalEfficiencyUncertainties[1])*(kfactor*fGlobalEfficiencyUncertainties[1]) )
1230                                         ) / fhRECpt->GetBinWidth(ibin);
1231       // max value with the minimum Nb
1232       errvalueExtremeMax =  TMath::Sqrt( ( (kfactor*fLuminosity[1]/fLuminosity[0])*(kfactor*fLuminosity[1]/fLuminosity[0]) ) +
1233                                          ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) +
1234                                          ( (kfactor*nbDmin/nb)*(kfactor*nbDmin/nb) )  +
1235                                          ( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))  ) +
1236                                          ( (kfactor*fGlobalEfficiencyUncertainties[1])*(kfactor*fGlobalEfficiencyUncertainties[1]) )
1237                                          ) / fhRECpt->GetBinWidth(ibin);
1238
1239       // Correction systematics (fc)
1240       // min value with the maximum Nb
1241       correctionMin = TMath::Sqrt( ( (kfactor*fLuminosity[1]/fLuminosity[0])*(kfactor*fLuminosity[1]/fLuminosity[0]) ) +
1242                                    ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) +
1243                                    ( (kfactor*nbDmax/nb)*(kfactor*nbDmax/nb) )  +
1244                                    ( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) ) +
1245                                    ( (kfactor*fGlobalEfficiencyUncertainties[1])*(kfactor*fGlobalEfficiencyUncertainties[1]) )
1246                                    ) / fhRECpt->GetBinContent(ibin) ;
1247       // max value with the minimum Nb
1248       correctionMax =  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*nbDmin/nb)*(kfactor*nbDmin/nb) )  +
1251                                     ( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))       ) +
1252                                     ( (kfactor*fGlobalEfficiencyUncertainties[1])*(kfactor*fGlobalEfficiencyUncertainties[1]) )
1253                                     ) / fhRECpt->GetBinContent(ibin) ;
1254     }
1255     else{ // Don't consider Nb uncertainty in this case [ to be tested!!! ]
1256       errvalueExtremeMax =  TMath::Sqrt( ( (kfactor*fLuminosity[1]/fLuminosity[0])*(kfactor*fLuminosity[1]/fLuminosity[0]) ) +
1257                                          ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) )  +
1258                                          ( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))  )  +
1259                                          ( (kfactor*fGlobalEfficiencyUncertainties[1])*(kfactor*fGlobalEfficiencyUncertainties[1]) )
1260                                          ) / fhRECpt->GetBinWidth(ibin);
1261       errvalueExtremeMin =  errvalueExtremeMax ;
1262     }
1263
1264
1265     // fill in histograms
1266     fhYieldCorr->SetBinContent(ibin,value);
1267     fhYieldCorr->SetBinError(ibin,errvalue);
1268     if (fAsymUncertainties) {
1269       Double_t x = fhYieldCorr->GetBinCenter(ibin);
1270       fgYieldCorr->SetPoint(ibin,x,value); // i,x,y
1271       fgYieldCorr->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueMin,errvalueMax); // i,xl,xh,yl,yh
1272       fhYieldCorrMax->SetBinContent(ibin,value+errvalueMax);
1273       fhYieldCorrMin->SetBinContent(ibin,value-errvalueMin);
1274       fgYieldCorrExtreme->SetPoint(ibin,x,value); // i,x,y
1275       fgYieldCorrExtreme->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueExtremeMin,errvalueExtremeMax); // i,xl,xh,yl,yh
1276       fgYieldCorrConservative->SetPoint(ibin,x,value); // i,x,y
1277       fgYieldCorrConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueExtremeMin,errvalueExtremeMax); // i,xl,xh,yl,yh
1278       //      cout << " bin " << ibin << ", correction " << correction << ", min correction unc " << correctionMin << ", max correction unc " << correctionMax << endl;
1279       if(correction>0.){
1280         fgFcConservative->SetPoint(ibin,x,correction);
1281         fgFcConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),correctionMin,correctionMax);
1282       }
1283       else{
1284         fgFcConservative->SetPoint(ibin,x,0.);
1285         fgFcConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),0.,0.);
1286       }
1287     }
1288
1289   }
1290   delete [] binwidths;
1291   delete [] limits;
1292
1293 }
1294
1295
1296 //_________________________________________________________________________________________________________
1297 void AliHFPtSpectrum::ComputeSystUncertainties(Int_t decay, Bool_t combineFeedDown) {
1298   //
1299   // Function that re-calculates the global systematic uncertainties
1300   //   by calling the class AliHFSystErr and combining those
1301   //   (in quadrature) with the feed-down subtraction uncertainties
1302   //
1303
1304   // Call the systematics uncertainty class for a given decay
1305   AliHFSystErr systematics(decay);
1306
1307   // Estimate the feed-down uncertainty in percentage
1308   Int_t nentries = fgSigmaCorrConservative->GetN();
1309   TGraphAsymmErrors *grErrFeeddown = new TGraphAsymmErrors(nentries);
1310   Double_t x=0., y=0., errx=0., erryl=0., erryh=0;
1311   for(Int_t i=0; i<nentries; i++) {
1312     fgSigmaCorrConservative->GetPoint(i,x,y);
1313     errx = fgSigmaCorrConservative->GetErrorXlow(i) ;
1314     erryl = fgSigmaCorrConservative->GetErrorYlow(i) / y ;
1315     erryh = fgSigmaCorrConservative->GetErrorYhigh(i) / y ;
1316     grErrFeeddown->SetPoint(i,x,0.);
1317     grErrFeeddown->SetPointError(i,errx,errx,erryl,erryh); //i, xl, xh, yl, yh
1318   }
1319
1320   // Draw all the systematics independently
1321   systematics.DrawErrors(grErrFeeddown);
1322
1323   // Set the sigma systematic uncertainties
1324   // possibly combine with the feed-down uncertainties
1325   Double_t errylcomb=0., erryhcomb=0;
1326   for(Int_t i=1; i<nentries; i++) {
1327     fgSigmaCorr->GetPoint(i,x,y);
1328     errx = grErrFeeddown->GetErrorXlow(i) ;
1329     erryl = grErrFeeddown->GetErrorYlow(i);
1330     erryh = grErrFeeddown->GetErrorYhigh(i);
1331     if (combineFeedDown) {
1332       errylcomb = systematics.GetTotalSystErr(x,erryl) * y ;
1333       erryhcomb = systematics.GetTotalSystErr(x,erryh) * y ;
1334     } else {
1335       errylcomb = systematics.GetTotalSystErr(x) * y ;
1336       erryhcomb = systematics.GetTotalSystErr(x) * y ;
1337     }
1338     fgSigmaCorr->SetPointError(i,errx,errx,errylcomb,erryhcomb);
1339   }
1340
1341 }
1342
1343
1344 //_________________________________________________________________________________________________________
1345 void AliHFPtSpectrum::DrawSpectrum(TGraphAsymmErrors *gPrediction) {
1346   //
1347   // Example method to draw the corrected spectrum & the theoretical prediction
1348   //
1349
1350   TCanvas *csigma = new TCanvas("csigma","Draw the corrected cross-section & the prediction");
1351   csigma->SetFillColor(0);
1352   gPrediction->GetXaxis()->SetTitleSize(0.05);
1353   gPrediction->GetXaxis()->SetTitleOffset(0.95);
1354   gPrediction->GetYaxis()->SetTitleSize(0.05);
1355   gPrediction->GetYaxis()->SetTitleOffset(0.95);
1356   gPrediction->GetXaxis()->SetTitle("p_{T}  [GeV]");
1357   gPrediction->GetYaxis()->SetTitle("BR #times #frac{d#sigma}{dp_{T}} |_{|y|<0.5}   [pb/GeV]");
1358   gPrediction->SetLineColor(kGreen+2);
1359   gPrediction->SetLineWidth(3);
1360   gPrediction->SetFillColor(kGreen+1);
1361   gPrediction->Draw("3CA");
1362   fgSigmaCorr->SetLineColor(kRed);
1363   fgSigmaCorr->SetLineWidth(1);
1364   fgSigmaCorr->SetFillColor(kRed);
1365   fgSigmaCorr->SetFillStyle(0);
1366   fgSigmaCorr->Draw("2");
1367   fhSigmaCorr->SetMarkerColor(kRed);
1368   fhSigmaCorr->Draw("esame");
1369   csigma->SetLogy();
1370   TLegend * leg = new TLegend(0.7,0.75,0.87,0.5);
1371   leg->SetBorderSize(0);
1372   leg->SetLineColor(0);
1373   leg->SetFillColor(0);
1374   leg->SetTextFont(42);
1375   leg->AddEntry(gPrediction,"FONLL ","fl");
1376   leg->AddEntry(fhSigmaCorr,"data stat. unc.","pl");
1377   leg->AddEntry(fgSigmaCorr,"data syst. unc.","f");
1378   leg->Draw();
1379   csigma->Draw();
1380
1381 }
1382
1383 //_________________________________________________________________________________________________________
1384 TH1D * AliHFPtSpectrum::ReweightHisto(TH1D *hToReweight, TH1D *hReference){
1385   //
1386   // Function to  reweight histograms for testing purposes:
1387   // This function takes the histo hToReweight and reweights
1388   //  it (its pt shape) with respect to hReference
1389   //
1390
1391   // check histograms consistency
1392   Bool_t areconsistent=kTRUE;
1393   areconsistent &= CheckHistosConsistency(hToReweight,hReference);
1394   if (!areconsistent) {
1395     AliInfo("the histograms to reweight are not consistent (bin width, bounds)");
1396     return NULL;
1397   }
1398
1399   // define a new empty histogram
1400   TH1D *hReweighted = (TH1D*)hToReweight->Clone("hReweighted");
1401   hReweighted->Reset();
1402   Double_t weight=1.0;
1403   Double_t yvalue=1.0;
1404   Double_t integralRef = hReference->Integral();
1405   Double_t integralH = hToReweight->Integral();
1406
1407   // now reweight the spectra
1408   //
1409   // the weight is the relative probability of the given pt bin in the reference histo
1410   //  divided by its relative probability (to normalize it) on the histo to re-weight
1411   for (Int_t i=0; i<=hToReweight->GetNbinsX(); i++) {
1412     weight = (hReference->GetBinContent(i)/integralRef) / (hToReweight->GetBinContent(i)/integralH) ;
1413     yvalue = hToReweight->GetBinContent(i);
1414     hReweighted->SetBinContent(i,yvalue*weight);
1415   }
1416
1417   return (TH1D*)hReweighted;
1418 }
1419
1420 //_________________________________________________________________________________________________________
1421 TH1D * AliHFPtSpectrum::ReweightRecHisto(TH1D *hRecToReweight, TH1D *hMCToReweight, TH1D *hMCReference){
1422   //
1423   // Function to  reweight histograms for testing purposes:
1424   // This function takes the histo hToReweight and reweights
1425   //  it (its pt shape) with respect to hReference /hMCToReweight
1426   //
1427
1428   // check histograms consistency
1429   Bool_t areconsistent=kTRUE;
1430   areconsistent &= CheckHistosConsistency(hMCToReweight,hMCReference);
1431   areconsistent &= CheckHistosConsistency(hRecToReweight,hMCReference);
1432   if (!areconsistent) {
1433     AliInfo("the histograms to reweight are not consistent (bin width, bounds)");
1434     return NULL;
1435   }
1436
1437   // define a new empty histogram
1438   TH1D *hReweighted = (TH1D*)hMCToReweight->Clone("hReweighted");
1439   hReweighted->Reset();
1440   TH1D *hRecReweighted = (TH1D*)hRecToReweight->Clone("hRecReweighted");
1441   hRecReweighted->Reset();
1442   Double_t weight=1.0;
1443   Double_t yvalue=1.0, yrecvalue=1.0;
1444   Double_t integralRef = hMCReference->Integral();
1445   Double_t integralH = hMCToReweight->Integral();
1446
1447   // now reweight the spectra
1448   //
1449   // the weight is the relative probability of the given pt bin
1450   //  that should be applied in the MC histo to get the reference histo shape
1451   //  Probabilities are properly normalized.
1452   for (Int_t i=0; i<=hMCToReweight->GetNbinsX(); i++) {
1453     weight = (hMCReference->GetBinContent(i)/integralRef) / (hMCToReweight->GetBinContent(i)/integralH) ;
1454     yvalue = hMCToReweight->GetBinContent(i);
1455     hReweighted->SetBinContent(i,yvalue*weight);
1456     yrecvalue = hRecToReweight->GetBinContent(i);
1457     hRecReweighted->SetBinContent(i,yrecvalue*weight);
1458   }
1459
1460   return (TH1D*)hRecReweighted;
1461 }
1462