PWGHF/vertexingHF/AliHFPtSpectrum.cxx

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