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   Int_t nbins = fhRECpt->GetNbinsX();
 572   if (fFeedDownOption==1) {
 573     // Compute the feed-down correction via fc-method
 574     CalculateFeedDownCorrectionFc();
 575     // Correct the yield for feed-down correction via fc-method
 576     CalculateFeedDownCorrectedSpectrumFc();
 577   }
 578   else if (fFeedDownOption==2) {
 579     // Correct the yield for feed-down correction via Nb-method
 580     CalculateFeedDownCorrectedSpectrumNb(deltaY,branchingRatioBintoFinalDecay);
 581   }
 582   else if (fFeedDownOption==0) {
 583     // If there is no need for feed-down correction,
 584     //    the "corrected" yield is equal to the raw yield
 585     fhYieldCorr = (TH1D*)fhRECpt->Clone();
 586     fhYieldCorr->SetNameTitle("fhYieldCorr","un-corrected yield");
 587     fhYieldCorrMax = (TH1D*)fhRECpt->Clone();
 588     fhYieldCorrMin = (TH1D*)fhRECpt->Clone();
 589     fhYieldCorrMax->SetNameTitle("fhYieldCorrMax","un-corrected yield");
 590     fhYieldCorrMin->SetNameTitle("fhYieldCorrMin","un-corrected yield");
 591     fgYieldCorr = new TGraphAsymmErrors(nbins+1);
 592     fAsymUncertainties=kFALSE;
 593   }
 594   else {
 595     AliInfo(" Are you sure the feed-down correction option is right ?");
 596   }
 597
 598
 599   // Print out information
 600   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]);
 601   if (fPbPbElossHypothesis)  printf("\n\n     The considered Tab is  %4.2e +- %2.2e \n\n",fTab[0],fTab[1]);
 602
 603   //
 604   // Finally: Correct from yields to cross-section
 605   //
 606   Double_t binwidth = fhRECpt->GetBinWidth(1);
 607   Double_t *limits = new Double_t[nbins+1];
 608   Double_t *binwidths = new Double_t[nbins];
 609   Double_t xlow=0.;
 610   for (Int_t i=1; i<=nbins; i++) {
 611     binwidth = fhRECpt->GetBinWidth(i);
 612     xlow = fhRECpt->GetBinLowEdge(i);
 613     limits[i-1] = xlow;
 614     binwidths[i-1] = binwidth;
 615   }
 616   limits[nbins] = xlow + binwidth;
 617
 618
 619   // declare the output histograms
 620   fhSigmaCorr = new TH1D("fhSigmaCorr","corrected sigma",nbins,limits);
 621   fgSigmaCorr = new TGraphAsymmErrors(nbins+1);
 622
 623   fhSigmaCorrMax = new TH1D("fhSigmaCorrMax","max corrected sigma",nbins,limits);
 624   fhSigmaCorrMin = new TH1D("fhSigmaCorrMin","min corrected sigma",nbins,limits);
 625   fhSigmaCorrDataSyst = new TH1D("fhSigmaCorrDataSyst","data syst uncertainties on the corrected sigma",nbins,limits);
 626
 627   if (fPbPbElossHypothesis && fFeedDownOption==1) {
 628     fhSigmaCorrRcb = new TH2D("fhSigmaCorrRcb","corrected sigma vs Rcb Eloss hypothesis; p_{T} [GeV/c] ; Rcb Eloss hypothesis ; #sigma",nbins,limits,800,0.,4.);
 629     fnSigma = new TNtuple("fnSigma"," Sigma ntuple calculation","pt:Signal:Rcb:fc:Yield:Sigma:SigmaStatUnc:SigmaMax:SigmaMin");
 630   }
 631   if (fPbPbElossHypothesis && fFeedDownOption==2) {
 632     fhSigmaCorrRcb = new TH2D("fhSigmaCorrRcb","corrected sigma vs Rb Eloss hypothesis; p_{T} [GeV/c] ; Rb Eloss hypothesis ; #sigma",nbins,limits,800,0.,4.);
 633     fnSigma = new TNtuple("fnSigma"," Sigma ntuple calculation","pt:Signal:Rb:fc:Yield:Sigma:SigmaStatUnc:SigmaMax:SigmaMin");
 634   }
 635   // and the output TGraphAsymmErrors
 636   if (fAsymUncertainties){
 637     fgSigmaCorrExtreme = new TGraphAsymmErrors(nbins+1);
 638     fgSigmaCorrConservative = new TGraphAsymmErrors(nbins+1);
 639   }
 640   fhStatUncEffcSigma = new TH1D("fhStatUncEffcSigma","direct charm stat unc on the cross section",nbins,limits);
 641   fhStatUncEffbSigma = new TH1D("fhStatUncEffbSigma","secondary charm stat unc on the cross section",nbins,limits);
 642
 643
 644   // protect against null denominator
 645   if (deltaY==0. || fLuminosity[0]==0. || fTrigEfficiency[0]==0. || branchingRatioC==0.) {
 646     AliError(" Hey you ! Why luminosity or trigger-efficiency or the c-BR or delta_y are set to zero ?! ");
 647     delete [] limits;
 648     delete [] binwidths;
 649     return ;
 650   }
 651
 652   Double_t value=0, errValue=0, errvalueMax=0., errvalueMin=0.;
 653   Double_t errvalueExtremeMax=0., errvalueExtremeMin=0.;
 654   Double_t errvalueConservativeMax=0., errvalueConservativeMin=0.;
 655   Double_t errvalueStatUncEffc=0.;
 656   for(Int_t ibin=1; ibin<=nbins; ibin++){
 657
 658     // Variables initialization
 659     value=0.; errValue=0.; errvalueMax=0.; errvalueMin=0.;
 660     errvalueExtremeMax=0.; errvalueExtremeMin=0.;
 661     errvalueConservativeMax=0.; errvalueConservativeMin=0.;
 662     errvalueStatUncEffc=0.;
 663
 664     // Sigma calculation
 665     //   Sigma = ( 1. / (lumi * delta_y * BR_c * ParticleAntiPartFactor * eff_trig * eff_c ) ) * spectra (corrected for feed-down)
 666     value = (fhDirectEffpt->GetBinContent(ibin) && fhDirectEffpt->GetBinContent(ibin)!=0. && fhRECpt->GetBinContent(ibin)>0.) ?
 667       ( fhYieldCorr->GetBinContent(ibin) / ( deltaY * branchingRatioC * fParticleAntiParticle * fLuminosity[0] * fTrigEfficiency[0] * fhDirectEffpt->GetBinContent(ibin) ) )
 668       : 0. ;
 669
 670     // Sigma statistical uncertainty:
 671     //   delta_sigma = sigma * sqrt ( (delta_spectra/spectra)^2 )
 672     errValue = (value!=0.) ?  value * (fhYieldCorr->GetBinError(ibin)/fhYieldCorr->GetBinContent(ibin))  : 0. ;
 673
 674     //    cout<< " x "<< fhRECpt->GetBinCenter(ibin) << " sigma " << value << " +- "<< errValue << " (stat)"<<endl;
 675
 676     //
 677     // Sigma systematic uncertainties
 678     //
 679     if (fAsymUncertainties && value>0.) {
 680
 681       //  (syst but feed-down) delta_sigma = sigma * sqrt ( (delta_spectra_syst/spectra)^2 +
 682       //                                     (delta_lumi/lumi)^2 + (delta_eff_trig/eff_trig)^2 + (delta_eff/eff)^2  + (global_eff)^2 )
 683       errvalueMax = value * TMath::Sqrt( (fgYieldCorr->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin))*(fgYieldCorr->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin)) +
 684                                          (fLuminosity[1]/fLuminosity[0])*(fLuminosity[1]/fLuminosity[0]) +
 685                                          (fTrigEfficiency[1]/fTrigEfficiency[0])*(fTrigEfficiency[1]/fTrigEfficiency[0])  +
 686                                          (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin)) +
 687                                          fGlobalEfficiencyUncertainties[0]*fGlobalEfficiencyUncertainties[0] );
 688       errvalueMin = value * TMath::Sqrt( (fgYieldCorr->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin))*(fgYieldCorr->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin)) +
 689                                          (fLuminosity[1]/fLuminosity[0])*(fLuminosity[1]/fLuminosity[0]) +
 690                                          (fTrigEfficiency[1]/fTrigEfficiency[0])*(fTrigEfficiency[1]/fTrigEfficiency[0])  +
 691                                          (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))  +
 692                                          fGlobalEfficiencyUncertainties[0]*fGlobalEfficiencyUncertainties[0] );
 693
 694       // Uncertainties from feed-down
 695       //      (feed-down syst) delta_sigma = sigma * sqrt ( (delta_spectra_fd/spectra_fd)^2 )
 696       //   extreme case
 697       errvalueExtremeMax = value * (fgYieldCorrExtreme->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin));
 698       errvalueExtremeMin =  value * (fgYieldCorrExtreme->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin));
 699       //
 700       //   conservative case
 701       errvalueConservativeMax = value * (fgYieldCorrConservative->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin));
 702       errvalueConservativeMin =  value * (fgYieldCorrConservative->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin));
 703
 704
 705       // stat unc of the efficiencies, separately
 706       errvalueStatUncEffc = value * (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin)) ;
 707
 708     }
 709     else {
 710       // protect against null denominator
 711       errvalueMax = (value!=0.) ?
 712         value * TMath::Sqrt( (fLuminosity[1]/fLuminosity[0])*(fLuminosity[1]/fLuminosity[0]) +
 713                              (fTrigEfficiency[1]/fTrigEfficiency[0])*(fTrigEfficiency[1]/fTrigEfficiency[0])  +
 714                              (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))  +
 715                              fGlobalEfficiencyUncertainties[0]*fGlobalEfficiencyUncertainties[0] )
 716         : 0. ;
 717       errvalueMin = errvalueMax;
 718     }
 719
 720     //
 721     // Fill the histograms
 722     //
 723     fhSigmaCorr->SetBinContent(ibin,value);
 724     fhSigmaCorr->SetBinError(ibin,errValue);
 725
 726     Double_t x = fhYieldCorr->GetBinCenter(ibin);
 727     fgSigmaCorr->SetPoint(ibin,x,value); // i,x,y
 728     fgSigmaCorr->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueMin,errvalueMax); // i,xl,xh,yl,yh
 729     fhSigmaCorrMax->SetBinContent(ibin,value+errvalueMax);
 730     fhSigmaCorrMin->SetBinContent(ibin,value-errvalueMin);
 731
 732     //
 733     // Fill the histos and ntuple vs the Eloss hypothesis
 734     //
 735     if (fPbPbElossHypothesis) {
 736
 737       // Loop over the Eloss hypothesis
 738       if(!fnHypothesis) {
 739         AliError("Error reading the fc hypothesis no ntuple, please check !!");
 740         delete [] limits;
 741         delete [] binwidths;
 742         return ;
 743       }
 744       Int_t entriesHypo = fnHypothesis->GetEntries();
 745       Float_t pt=0., Rhy=0., fc=0., fcMin=0., fcMax=0.;
 746       fnHypothesis->SetBranchAddress("pt",&pt);
 747       if (fFeedDownOption==2) fnHypothesis->SetBranchAddress("Rb",&Rhy);
 748       else if (fFeedDownOption==1) fnHypothesis->SetBranchAddress("Rcb",&Rhy);
 749       fnHypothesis->SetBranchAddress("fc",&fc);
 750       fnHypothesis->SetBranchAddress("fcMax",&fcMax);
 751       fnHypothesis->SetBranchAddress("fcMin",&fcMin);
 752
 753       for (Int_t item=0; item<entriesHypo; item++) {
 754
 755         fnHypothesis->GetEntry(item);
 756         if ( TMath::Abs( pt - fhDirectEffpt->GetBinCenter(ibin) ) > 0.15 ) continue;
 757
 758         Double_t yieldRcbvalue = (fhRECpt->GetBinContent(ibin) ) ? fhRECpt->GetBinContent(ibin) * fc : 0. ;
 759         yieldRcbvalue /= fhRECpt->GetBinWidth(ibin) ;
 760         Double_t yieldRcbvalueMax = (fhRECpt->GetBinContent(ibin) ) ? fhRECpt->GetBinContent(ibin) * fcMax : 0. ;
 761         yieldRcbvalueMax /= fhRECpt->GetBinWidth(ibin) ;
 762         Double_t yieldRcbvalueMin = (fhRECpt->GetBinContent(ibin) ) ? fhRECpt->GetBinContent(ibin) * fcMin : 0. ;
 763         yieldRcbvalueMin /= fhRECpt->GetBinWidth(ibin) ;
 764
 765         // Sigma calculation
 766         //   Sigma = ( 1. / (lumi * delta_y * BR_c * ParticleAntiPartFactor * eff_trig * eff_c ) ) * spectra (corrected for feed-down)
 767         Double_t sigmaRcbvalue = (fhDirectEffpt->GetBinContent(ibin) && fhDirectEffpt->GetBinContent(ibin)>0.) ?
 768           ( yieldRcbvalue / ( deltaY * branchingRatioC * fParticleAntiParticle * fLuminosity[0] * fTrigEfficiency[0] * fhDirectEffpt->GetBinContent(ibin) ) )
 769           : 0. ;
 770         Double_t sigmaRcbvalueMax = (sigmaRcbvalue!=0.) ?
 771           ( yieldRcbvalueMax / ( deltaY * branchingRatioC * fParticleAntiParticle * fLuminosity[0] * fTrigEfficiency[0] * fhDirectEffpt->GetBinContent(ibin) ) )
 772           : 0. ;
 773         Double_t sigmaRcbvalueMin = (sigmaRcbvalue!=0.) ?
 774           ( yieldRcbvalueMin / ( deltaY * branchingRatioC * fParticleAntiParticle * fLuminosity[0] * fTrigEfficiency[0] * fhDirectEffpt->GetBinContent(ibin) ) )
 775           : 0. ;
 776         // Sigma statistical uncertainty:
 777         //   delta_sigma = sigma * sqrt ( (delta_spectra/spectra)^2 ) = sigma * ( delta_spectra / (spectra-corr * binwidth) )
 778         Double_t sigmaRcbvalueStatUnc = (sigmaRcbvalue!=0.) ?
 779           sigmaRcbvalue * ( fhRECpt->GetBinError(ibin) / ( yieldRcbvalue * fhRECpt->GetBinWidth(ibin) ) )  : 0. ;
 780
 781         fhSigmaCorrRcb->Fill( fhSigmaCorr->GetBinCenter(ibin) , Rhy, sigmaRcbvalue );
 782         //      if(ibin==3)
 783         //        cout << " pt "<< fhRECpt->GetBinCenter(ibin) <<" bin "<< ibin<<" rval="<<rval<<", rbin="<<rbin<<" fc-value="<< fhFcRcb->GetBinContent(ibin,rbin) <<", yield-fcRbvalue="<<yieldRcbvalue<<", sigma-fcRbvalue="<<sigmaRcbvalue<<endl;
 784         fnSigma->Fill(fhRECpt->GetBinCenter(ibin), fhRECpt->GetBinContent(ibin),
 785                       Rhy, fc, yieldRcbvalue, sigmaRcbvalue, sigmaRcbvalueStatUnc, sigmaRcbvalueMax, sigmaRcbvalueMin );
 786       }
 787     }
 788     //
 789     // Fill the TGraphAsymmErrors
 790     if (fAsymUncertainties) {
 791       fgSigmaCorrExtreme->SetPoint(ibin,x,value); // i,x,y
 792       fgSigmaCorrExtreme->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueExtremeMin,errvalueExtremeMax); // i,xl,xh,yl,yh
 793       fgSigmaCorrConservative->SetPoint(ibin,x,value); // i,x,y
 794       fgSigmaCorrConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueConservativeMin,errvalueConservativeMax); // i,xl,xh,yl,yh
 795
 796       fhStatUncEffcSigma->SetBinContent(ibin,0.);
 797       if(value>0.) fhStatUncEffcSigma->SetBinError(ibin,((errvalueStatUncEffc/value)*100.));
 798       fhStatUncEffbSigma->SetBinContent(ibin,0.); fhStatUncEffbSigma->SetBinError(ibin,0.);
 799       //      cout << " pt "<< fhRECpt->GetBinCenter(ibin) <<" bin "<< ibin<<" stat-unc-c-sigma "<< errvalueStatUncEffc/value << endl;
 800     }
 801
 802   }
 803   delete [] binwidths;
 804   delete [] limits;
 805
 806 }
 807
 808 //_________________________________________________________________________________________________________
 809 TH1D * AliHFPtSpectrum::EstimateEfficiencyRecoBin(TH1D *hSimu, TH1D *hReco, const char *name) {
 810   //
 811   // Function that computes the acceptance and efficiency correction
 812   //  based on the simulated and reconstructed spectra
 813   //  and using the reconstructed spectra bin width
 814   //
 815   //  eff = reco/sim ; err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
 816   //
 817
 818   if(!fhRECpt){
 819     AliInfo("Hey, the reconstructed histogram was not set yet !");
 820     return NULL;
 821   }
 822
 823   Int_t nbins = fhRECpt->GetNbinsX();
 824   Double_t *limits = new Double_t[nbins+1];
 825   Double_t xlow=0.,binwidth=0.;
 826   for (Int_t i=1; i<=nbins; i++) {
 827     binwidth = fhRECpt->GetBinWidth(i);
 828     xlow = fhRECpt->GetBinLowEdge(i);
 829     limits[i-1] = xlow;
 830   }
 831   limits[nbins] = xlow + binwidth;
 832
 833   TH1D * hEfficiency = new TH1D(name," acceptance #times efficiency",nbins,limits);
 834
 835   Double_t *sumSimu=new Double_t[nbins];
 836   Double_t *sumReco=new Double_t[nbins];
 837   for (Int_t ibin=0; ibin<nbins; ibin++){
 838     sumSimu[ibin]=0.;  sumReco[ibin]=0.;
 839   }
 840   for (Int_t ibin=0; ibin<=hSimu->GetNbinsX(); ibin++){
 841
 842     for (Int_t ibinrec=0; ibinrec<nbins; ibinrec++){
 843       if ( hSimu->GetBinCenter(ibin)>limits[ibinrec] &&
 844            hSimu->GetBinCenter(ibin)<limits[ibinrec+1] ) {
 845         sumSimu[ibinrec]+=hSimu->GetBinContent(ibin);
 846       }
 847       if ( hReco->GetBinCenter(ibin)>limits[ibinrec] &&
 848            hReco->GetBinCenter(ibin)<limits[ibinrec+1] ) {
 849         sumReco[ibinrec]+=hReco->GetBinContent(ibin);
 850       }
 851     }
 852
 853   }
 854
 855
 856   // the efficiency is computed as reco/sim (in each bin)
 857   //  its uncertainty is err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
 858   Double_t eff=0., erreff=0.;
 859   for (Int_t ibinrec=0; ibinrec<nbins; ibinrec++) {
 860     if (sumSimu[ibinrec]!= 0. && sumReco[ibinrec]!=0.) {
 861       eff = sumReco[ibinrec] / sumSimu[ibinrec] ;
 862       // protection in case eff > 1.0
 863       // test calculation (make the argument of the sqrt positive)
 864       erreff = TMath::Sqrt( eff * TMath::Abs(1.0 - eff) ) / TMath::Sqrt( sumSimu[ibinrec] );
 865     }
 866     else { eff=0.0; erreff=0.; }
 867     hEfficiency->SetBinContent(ibinrec+1,eff);
 868     hEfficiency->SetBinError(ibinrec+1,erreff);
 869   }
 870
 871   delete [] sumSimu;
 872   delete [] sumReco;
 873
 874   return (TH1D*)hEfficiency;
 875 }
 876
 877 //_________________________________________________________________________________________________________
 878 void AliHFPtSpectrum::EstimateAndSetDirectEfficiencyRecoBin(TH1D *hSimu, TH1D *hReco) {
 879   //
 880   // Function that computes the Direct  acceptance and efficiency correction
 881   //  based on the simulated and reconstructed spectra
 882   //  and using the reconstructed spectra bin width
 883   //
 884   //  eff = reco/sim ; err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
 885   //
 886
 887   if(!fhRECpt || !hSimu || !hReco){
 888     AliError("Hey, the reconstructed histogram was not set yet !");
 889     return;
 890   }
 891
 892   fhDirectEffpt = EstimateEfficiencyRecoBin(hSimu,hReco,"fhDirectEffpt");
 893   fhDirectEffpt->SetNameTitle("fhDirectEffpt"," direct acceptance #times efficiency");
 894
 895 }
 896
 897 //_________________________________________________________________________________________________________
 898 void AliHFPtSpectrum::EstimateAndSetFeedDownEfficiencyRecoBin(TH1D *hSimu, TH1D *hReco) {
 899   //
 900   // Function that computes the Feed-Down acceptance and efficiency correction
 901   //  based on the simulated and reconstructed spectra
 902   //  and using the reconstructed spectra bin width
 903   //
 904   //  eff = reco/sim ; err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
 905   //
 906
 907   if(!fhRECpt || !hSimu || !hReco){
 908     AliError("Hey, the reconstructed histogram was not set yet !");
 909     return;
 910   }
 911
 912   fhFeedDownEffpt = EstimateEfficiencyRecoBin(hSimu,hReco,"fhFeedDownEffpt");
 913   fhFeedDownEffpt->SetNameTitle("fhFeedDownEffpt"," feed-down acceptance #times efficiency");
 914
 915 }
 916
 917 //_________________________________________________________________________________________________________
 918 Bool_t AliHFPtSpectrum::Initialize(){
 919   //
 920   // Initialization of the variables (histograms)
 921   //
 922
 923   if (fFeedDownOption==0) {
 924     AliInfo("Getting ready for the corrections without feed-down consideration");
 925   } else if (fFeedDownOption==1) {
 926     AliInfo("Getting ready for the fc feed-down correction calculation");
 927   } else if (fFeedDownOption==2) {
 928     AliInfo("Getting ready for the Nb feed-down correction calculation");
 929   } else { AliError("The calculation option must be <=2"); return kFALSE; }
 930
 931   // Start checking the input histograms consistency
 932   Bool_t areconsistent=kTRUE;
 933
 934   // General checks
 935   if (!fhDirectEffpt || !fhRECpt) {
 936     AliError(" Reconstructed spectra and/or the Nc efficiency distributions are not defined");
 937     return kFALSE;
 938   }
 939   areconsistent &= CheckHistosConsistency(fhRECpt,fhDirectEffpt);
 940   if (!areconsistent) {
 941     AliInfo("Histograms required for Nb correction are not consistent (bin width, bounds)");
 942     return kFALSE;
 943   }
 944   if (fFeedDownOption==0) return kTRUE;
 945
 946   //
 947   // Common checks for options 1 (fc) & 2(Nb)
 948   if (!fhFeedDownMCpt || !fhFeedDownEffpt) {
 949     AliError(" Theoretical Nb and/or the Nb efficiency distributions are not defined");
 950     return kFALSE;
 951   }
 952   areconsistent &= CheckHistosConsistency(fhRECpt,fhFeedDownMCpt);
 953   areconsistent &= CheckHistosConsistency(fhFeedDownMCpt,fhFeedDownEffpt);
 954   if (fAsymUncertainties) {
 955     if (!fhFeedDownMCptMax || !fhFeedDownMCptMin) {
 956       AliError(" Max/Min theoretical Nb distributions are not defined");
 957       return kFALSE;
 958     }
 959     areconsistent &= CheckHistosConsistency(fhFeedDownMCpt,fhFeedDownMCptMax);
 960   }
 961   if (!areconsistent) {
 962     AliInfo("Histograms required for Nb correction are not consistent (bin width, bounds)");
 963     return kFALSE;
 964   }
 965   if (fFeedDownOption>1) return kTRUE;
 966
 967   //
 968   // Now checks for option 1 (fc correction)
 969   if (!fhDirectMCpt) {
 970     AliError("Theoretical Nc distributions is not defined");
 971     return kFALSE;
 972   }
 973   areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhFeedDownMCpt);
 974   areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhDirectEffpt);
 975   if (fAsymUncertainties) {
 976     if (!fhDirectMCptMax || !fhDirectMCptMin) {
 977       AliError(" Max/Min theoretical Nc distributions are not defined");
 978       return kFALSE;
 979     }
 980     areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhDirectMCptMax);
 981   }
 982   if (!areconsistent) {
 983     AliInfo("Histograms required for fc correction are not consistent (bin width, bounds)");
 984     return kFALSE;
 985   }
 986
 987   return kTRUE;
 988 }
 989
 990 //_________________________________________________________________________________________________________
 991 Bool_t AliHFPtSpectrum::CheckHistosConsistency(TH1D *h1, TH1D *h2){
 992   //
 993   // Check the histograms consistency (bins, limits)
 994   //
 995
 996   if (!h1 || !h2) {
 997     AliError("One or both histograms don't exist");
 998     return kFALSE;
 999   }
1000
1001   Double_t binwidth1 = h1->GetBinWidth(1);
1002   Double_t binwidth2 = h2->GetBinWidth(1);
1003   Double_t min1 = h1->GetBinCenter(1) - (binwidth1/2.) ;
1004 //   Double_t max1 = h1->GetBinCenter(nbins1) + (binwidth1/2.) ;
1005   Double_t min2 = h2->GetBinCenter(1) - (binwidth2/2.) ;
1006 //   Double_t max2 = h2->GetBinCenter(nbins2) + (binwidth2/2.) ;
1007
1008   if (binwidth1!=binwidth2) {
1009     AliInfo(" histograms with different bin width");
1010     return kFALSE;
1011   }
1012   if (min1!=min2) {
1013     AliInfo(" histograms with different minimum");
1014     return kFALSE;
1015   }
1016 //   if (max1!=max2) {
1017 //     AliInfo(" histograms with different maximum");
1018 //     return kFALSE;
1019 //   }
1020
1021   return kTRUE;
1022 }
1023
1024 //_________________________________________________________________________________________________________
1025 void AliHFPtSpectrum::CalculateFeedDownCorrectionFc(){
1026   //
1027   // Compute fc factor and its uncertainties bin by bin
1028   //   fc = 1 / ( 1 + (eff_b/eff_c)*(N_b/N_c) )
1029   //
1030   // uncertainties: (conservative) combine the upper/lower N_b & N_c predictions together
1031   //                (extreme) combine the upper N_b predictions with the lower N_c predictions & viceversa
1032   //                systematic uncertainty on the acceptance x efficiency b/c ratio are included
1033   //
1034   //  In addition, in HIC the feed-down correction varies with an energy loss hypothesis: Raa(c-->D) / Raa(b-->D) = Rcb
1035   //           fc (Rcb) = ( 1. / ( 1 + (eff_b/eff_c)*(N_b/N_c)* (1/Rcb) ) );
1036   //
1037   AliInfo("Calculating the feed-down correction factor (fc method)");
1038
1039   // define the variables
1040   Int_t nbins = fhRECpt->GetNbinsX();
1041   Double_t binwidth = fhRECpt->GetBinWidth(1);
1042   Double_t *limits = new Double_t[nbins+1];
1043   Double_t *binwidths = new Double_t[nbins];
1044   Double_t xlow=0.;
1045   for (Int_t i=1; i<=nbins; i++) {
1046     binwidth = fhRECpt->GetBinWidth(i);
1047     xlow = fhRECpt->GetBinLowEdge(i);
1048     limits[i-1] = xlow;
1049     binwidths[i-1] = binwidth;
1050   }
1051   limits[nbins] = xlow + binwidth;
1052
1053   Double_t correction=1.;
1054   Double_t theoryRatio=1.;
1055   Double_t effRatio=1.;
1056   Double_t correctionExtremeA=1., correctionExtremeB=1.;
1057   Double_t theoryRatioExtremeA=1., theoryRatioExtremeB=1.;
1058   Double_t correctionConservativeA=1., correctionConservativeB=1.;
1059   Double_t theoryRatioConservativeA=1., theoryRatioConservativeB=1.;
1060   //  Double_t correctionUnc=0.;
1061   Double_t correctionExtremeAUnc=0., correctionExtremeBUnc=0.;
1062   Double_t correctionConservativeAUnc=0., correctionConservativeBUnc=0.;
1063
1064   // declare the output histograms
1065   fhFc = new TH1D("fhFc","fc correction factor",nbins,limits);
1066   fhFcMax = new TH1D("fhFcMax","max fc correction factor",nbins,limits);
1067   fhFcMin = new TH1D("fhFcMin","min fc correction factor",nbins,limits);
1068   if(fPbPbElossHypothesis) {
1069     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.);
1070     fnHypothesis = new TNtuple("fnHypothesis"," Feed-down correction vs hypothesis (fc)","pt:Rcb:fc:fcMax:fcMin");
1071   }
1072   // two local control histograms
1073   TH1D *hTheoryRatio = new TH1D("hTheoryRatio","Theoretical B-->D over c-->D (feed-down/direct) ratio",nbins,limits);
1074   TH1D *hEffRatio = new TH1D("hEffRatio","Efficiency B-->D over c-->D (feed-down/direct) ratio",nbins,limits);
1075   // and the output TGraphAsymmErrors
1076   if (fAsymUncertainties) {
1077     fgFcExtreme = new TGraphAsymmErrors(nbins+1);
1078     fgFcExtreme->SetNameTitle("fgFcExtreme","fgFcExtreme");
1079     fgFcConservative = new TGraphAsymmErrors(nbins+1);
1080     fgFcConservative->SetNameTitle("fgFcConservative","fgFcConservative");
1081   }
1082
1083   fhStatUncEffcFD = new TH1D("fhStatUncEffcFD","direct charm stat unc on the feed-down correction",nbins,limits);
1084   fhStatUncEffbFD = new TH1D("fhStatUncEffbFD","secondary charm stat unc on the feed-down correction",nbins,limits);
1085   Double_t correctionConservativeAUncStatEffc=0., correctionConservativeBUncStatEffc=0.;
1086   Double_t correctionConservativeAUncStatEffb=0., correctionConservativeBUncStatEffb=0.;
1087
1088   //
1089   // Compute fc
1090   //
1091   for (Int_t ibin=1; ibin<=nbins; ibin++) {
1092
1093     // Variables initialization
1094     correction=1.; theoryRatio=1.; effRatio=1.;
1095     correctionExtremeA=1.; correctionExtremeB=1.;
1096     theoryRatioExtremeA=1.; theoryRatioExtremeB=1.;
1097     correctionConservativeA=1.; correctionConservativeB=1.;
1098     theoryRatioConservativeA=1.; theoryRatioConservativeB=1.;
1099     //    correctionUnc=0.;
1100     correctionExtremeAUnc=0.; correctionExtremeBUnc=0.;
1101     correctionConservativeAUnc=0.; correctionConservativeBUnc=0.;
1102     correctionConservativeAUncStatEffc=0.; correctionConservativeBUncStatEffc=0.;
1103     correctionConservativeAUncStatEffb=0.; correctionConservativeBUncStatEffb=0.;
1104
1105     //  theory_ratio = (N_b/N_c)
1106     theoryRatio = (fhDirectMCpt->GetBinContent(ibin)>0. && fhFeedDownMCpt->GetBinContent(ibin)>0.) ?
1107       fhFeedDownMCpt->GetBinContent(ibin) / fhDirectMCpt->GetBinContent(ibin) : 1.0 ;
1108
1109     //
1110     // Calculate the uncertainty [ considering only the theoretical uncertainties on Nb & Nc for now !!! ]
1111     //
1112     // extreme A = direct-max, feed-down-min
1113     theoryRatioExtremeA = (fhDirectMCptMax->GetBinContent(ibin)>0. && fhFeedDownMCptMin->GetBinContent(ibin)>0.) ?
1114       fhFeedDownMCptMin->GetBinContent(ibin) / fhDirectMCptMax->GetBinContent(ibin) : 1.0 ;
1115     // extreme B = direct-min, feed-down-max
1116     theoryRatioExtremeB = (fhDirectMCptMin->GetBinContent(ibin)>0. && fhDirectMCptMax->GetBinContent(ibin)>0.) ?
1117       fhFeedDownMCptMax->GetBinContent(ibin) / fhDirectMCptMin->GetBinContent(ibin) : 1.0 ;
1118     // conservative A = direct-max, feed-down-max
1119     theoryRatioConservativeA = (fhDirectMCptMax->GetBinContent(ibin)>0. && fhFeedDownMCptMin->GetBinContent(ibin)>0.) ?
1120       fhFeedDownMCptMax->GetBinContent(ibin) / fhDirectMCptMax->GetBinContent(ibin) : 1.0 ;
1121     // conservative B = direct-min, feed-down-min
1122     theoryRatioConservativeB = (fhDirectMCptMin->GetBinContent(ibin)>0. && fhDirectMCptMax->GetBinContent(ibin)>0.) ?
1123       fhFeedDownMCptMin->GetBinContent(ibin) / fhDirectMCptMin->GetBinContent(ibin) : 1.0 ;
1124
1125     //  eff_ratio = (eff_b/eff_c)
1126     effRatio = (fhDirectEffpt->GetBinContent(ibin) && fhDirectEffpt->GetBinContent(ibin)!=0.) ?
1127       fhFeedDownEffpt->GetBinContent(ibin) / fhDirectEffpt->GetBinContent(ibin) : 1.0 ;
1128
1129     //   fc = 1 / ( 1 + (eff_b/eff_c)*(N_b/N_c) )
1130     if( TMath::Abs(effRatio - 1.0)<0.0001 || TMath::Abs(theoryRatio - 1.0)<0.0001 ) {
1131       correction = 1.0;
1132       correctionExtremeA = 1.0;
1133       correctionExtremeB = 1.0;
1134       correctionConservativeA = 1.0;
1135       correctionConservativeB = 1.0;
1136     }
1137     else {
1138       correction = ( 1. / ( 1 + ( effRatio * theoryRatio ) ) );
1139       correctionExtremeA = ( 1. / ( 1 + ( effRatio * theoryRatioExtremeA ) ) );
1140       correctionExtremeB = ( 1. / ( 1 + ( effRatio * theoryRatioExtremeB ) ) );
1141       correctionConservativeA = ( 1. / ( 1 + ( effRatio * theoryRatioConservativeA ) ) );
1142       correctionConservativeB = ( 1. / ( 1 + ( effRatio * theoryRatioConservativeB ) ) );
1143     }
1144
1145
1146     // fc uncertainty from (eff_b/eff_c) = fc^2 * (N_b/N_c) * delta(eff_b/eff_c)
1147     //  delta(eff_b/eff_c) is a percentage = effRatio * sqrt( fGlobalEfficiencyUncertainties[1]^2 + unc_eff_c ^2 + unc_eff_b ^2 )
1148     Double_t relEffUnc =  TMath::Sqrt( fGlobalEfficiencyUncertainties[1]*fGlobalEfficiencyUncertainties[1] +
1149                                        (fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) +
1150                                        (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))
1151                                        );
1152
1153     //    correctionUnc = correction*correction * theoryRatio * effRatio * relEffUnc;
1154     correctionExtremeAUnc = correctionExtremeA*correctionExtremeA * theoryRatioExtremeA  * effRatio * relEffUnc;
1155     correctionExtremeBUnc = correctionExtremeB*correctionExtremeB * theoryRatioExtremeB  * effRatio * relEffUnc;
1156
1157     correctionConservativeAUnc = correctionConservativeA*correctionConservativeA * theoryRatioConservativeA  *effRatio * relEffUnc;
1158     //
1159     correctionConservativeAUncStatEffc = correctionConservativeA*correctionConservativeA * theoryRatioConservativeA  *effRatio *
1160       (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin));
1161     correctionConservativeAUncStatEffb = correctionConservativeA*correctionConservativeA * theoryRatioConservativeA  *effRatio *
1162       (fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin));
1163
1164     correctionConservativeBUnc = correctionConservativeB*correctionConservativeB * theoryRatioConservativeB  *effRatio * relEffUnc;
1165
1166     correctionConservativeBUncStatEffb = correctionConservativeB*correctionConservativeB * theoryRatioConservativeB  *effRatio *
1167       (fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin));
1168     correctionConservativeBUncStatEffc = correctionConservativeB*correctionConservativeB * theoryRatioConservativeB  *effRatio *
1169       (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin));
1170
1171
1172     // Fill in the histograms
1173     hTheoryRatio->SetBinContent(ibin,theoryRatio);
1174     hEffRatio->SetBinContent(ibin,effRatio);
1175     fhFc->SetBinContent(ibin,correction);
1176     //
1177     // Estimate how the result varies vs charm/beauty Eloss hypothesis
1178     //
1179     if ( correction>1.0e-16 && fPbPbElossHypothesis){
1180       // Loop over the Eloss hypothesis
1181       //      Int_t rbin=0;
1182       for (Float_t rval=0.0025; rval<4.0; rval+=0.005){
1183         // Central fc with Eloss hypothesis
1184         Double_t correctionRcb = ( 1. / ( 1 + ( effRatio * theoryRatio * (1/rval) ) ) );
1185         fhFcRcb->Fill( fhFc->GetBinCenter(ibin) , rval, correctionRcb );
1186         //      if(ibin==3){
1187         //        cout << " pt "<< fhFc->GetBinCenter(ibin) <<" bin "<< ibin<<" rval="<<rval<<", rbin="<<rbin<<", fc-Rcb-value="<<correctionRcb<<endl;
1188         //        rbin++;
1189         //      }
1190         // Upper / lower fc with up / low FONLL bands and Eloss hypothesis
1191         Double_t correctionConservativeARcb = ( 1. / ( 1 + ( effRatio * theoryRatioConservativeA * (1/rval) ) ) );
1192         Double_t correctionConservativeBRcb = ( 1. / ( 1 + ( effRatio * theoryRatioConservativeB * (1/rval) ) ) );
1193         Double_t correctionConservativeARcbUnc = correctionConservativeARcb*correctionConservativeARcb * theoryRatioConservativeA * (1/rval) *effRatio * relEffUnc;
1194         Double_t correctionConservativeBRcbUnc = correctionConservativeBRcb*correctionConservativeBRcb * theoryRatioConservativeB * (1/rval) *effRatio * relEffUnc;
1195         Double_t consvalRcb[4] = { correctionConservativeARcb - correctionConservativeARcbUnc, correctionConservativeARcb + correctionConservativeARcbUnc,
1196                                    correctionConservativeBRcb - correctionConservativeBRcbUnc, correctionConservativeBRcb + correctionConservativeBRcbUnc};
1197         Double_t uncConservativeRcbMin = correctionRcb - TMath::MinElement(4,consvalRcb);
1198         Double_t uncConservativeRcbMax =  TMath::MaxElement(4,consvalRcb) - correctionRcb;
1199 //      if(ibin==3)
1200 //        cout << " pt="<<fhDirectEffpt->GetBinCenter(ibin)<<", hypo="<<rval<<", fc="<<correctionRcb<<" +"<<uncConservativeRcbMax<<" -"<<uncConservativeRcbMin<<endl;
1201         fnHypothesis->Fill( fhDirectEffpt->GetBinCenter(ibin), rval, correctionRcb, correctionRcb+uncConservativeRcbMax, correctionRcb-uncConservativeRcbMin);
1202       }
1203     }
1204     //
1205     // Fill the rest of (asymmetric) histograms
1206     //
1207     if (fAsymUncertainties) {
1208       Double_t x = fhDirectMCpt->GetBinCenter(ibin);
1209       Double_t val[4] = { correctionExtremeA + correctionExtremeAUnc, correctionExtremeA - correctionExtremeAUnc,
1210                           correctionExtremeB + correctionExtremeBUnc, correctionExtremeB - correctionExtremeBUnc };
1211       Double_t uncExtremeMin = correction - TMath::MinElement(4,val);
1212       Double_t uncExtremeMax = TMath::MaxElement(4,val) - correction;
1213       fgFcExtreme->SetPoint(ibin,x,correction); // i,x,y
1214       fgFcExtreme->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),uncExtremeMin,uncExtremeMax); // i,xl,xh,yl,yh
1215       fhFcMax->SetBinContent(ibin,correction+uncExtremeMax);
1216       fhFcMin->SetBinContent(ibin,correction-uncExtremeMin);
1217       Double_t consval[4] = { correctionConservativeA - correctionConservativeAUnc, correctionConservativeA + correctionConservativeAUnc,
1218                               correctionConservativeB - correctionConservativeBUnc, correctionConservativeB + correctionConservativeBUnc};
1219       Double_t uncConservativeMin = correction - TMath::MinElement(4,consval);
1220       Double_t uncConservativeMax =  TMath::MaxElement(4,consval) - correction;
1221       fgFcConservative->SetPoint(ibin,x,correction); // i,x,y
1222       fgFcConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),uncConservativeMin,uncConservativeMax); // i,xl,xh,yl,yh
1223       if( !(correction>0.) ){
1224         fgFcExtreme->SetPoint(ibin,x,0.); // i,x,y
1225         fgFcExtreme->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),0.,0.); // i,xl,xh,yl,yh
1226         fgFcConservative->SetPoint(ibin,x,0.); // i,x,y
1227         fgFcConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),0.,0.); // i,xl,xh,yl,yh
1228       }
1229
1230       Double_t valStatEffc[2] = { correctionConservativeAUncStatEffc/correctionConservativeA,
1231                                   correctionConservativeBUncStatEffc/correctionConservativeB };
1232       Double_t valStatEffb[2] = { correctionConservativeAUncStatEffb/correctionConservativeA,
1233                                   correctionConservativeBUncStatEffb/correctionConservativeB };
1234       Double_t uncConservativeStatEffc = TMath::MaxElement(2,valStatEffc);
1235       Double_t uncConservativeStatEffb = TMath::MaxElement(2,valStatEffb);
1236       fhStatUncEffcFD->SetBinContent(ibin,0.); fhStatUncEffcFD->SetBinError(ibin,uncConservativeStatEffc*100.);
1237       fhStatUncEffbFD->SetBinContent(ibin,0.); fhStatUncEffbFD->SetBinError(ibin,uncConservativeStatEffb*100.);
1238       //      cout << " pt "<< fhStatUncEffcFD->GetBinCenter(ibin) <<" bin "<< ibin<<" fc-stat-c ="<<uncConservativeStatEffc<<" fc-stat-b ="<<uncConservativeStatEffb<<endl;
1239     }
1240
1241   }
1242   delete [] binwidths;
1243   delete [] limits;
1244
1245 }
1246
1247 //_________________________________________________________________________________________________________
1248 void AliHFPtSpectrum::CalculateFeedDownCorrectedSpectrumFc(){
1249   //
1250   // Compute the feed-down corrected spectrum if feed-down correction is done via fc factor (bin by bin)
1251   //    physics = reco * fc / bin-width
1252   //
1253   //    uncertainty:             (stat) delta_physics = physics * sqrt ( (delta_reco/reco)^2 )
1254   //               (syst but feed-down) delta_physics = physics * sqrt ( (delta_reco_syst/reco)^2 )
1255   //                   (feed-down syst) delta_physics = physics * sqrt ( (delta_fc/fc)^2 )
1256   //
1257   //    ( Calculation done bin by bin )
1258   //
1259   //  In addition, in HIC the feed-down correction varies with an energy loss hypothesis: Raa(c-->D) / Raa(b-->D) = Rcb
1260
1261   AliInfo(" Calculating the feed-down corrected spectrum (fc method)");
1262
1263   if (!fhFc || !fhRECpt) {
1264     AliError(" Reconstructed or fc distributions are not defined");
1265     return;
1266   }
1267
1268   Int_t nbins = fhRECpt->GetNbinsX();
1269   Double_t value = 0., errvalue = 0., errvalueMax= 0., errvalueMin= 0.;
1270   Double_t valueExtremeMax= 0., valueExtremeMin= 0.;
1271   Double_t valueConservativeMax= 0., valueConservativeMin= 0.;
1272   Double_t binwidth = fhRECpt->GetBinWidth(1);
1273   Double_t *limits = new Double_t[nbins+1];
1274   Double_t *binwidths = new Double_t[nbins];
1275   Double_t xlow=0.;
1276   for (Int_t i=1; i<=nbins; i++) {
1277     binwidth = fhRECpt->GetBinWidth(i);
1278     xlow = fhRECpt->GetBinLowEdge(i);
1279     limits[i-1] = xlow;
1280     binwidths[i-1] = binwidth;
1281   }
1282   limits[nbins] = xlow + binwidth;
1283
1284   // declare the output histograms
1285   fhYieldCorr = new TH1D("fhYieldCorr","corrected yield (by fc)",nbins,limits);
1286   fhYieldCorrMax = new TH1D("fhYieldCorrMax","max corrected yield (by fc)",nbins,limits);
1287   fhYieldCorrMin = new TH1D("fhYieldCorrMin","min corrected yield (by fc)",nbins,limits);
1288   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.);
1289   // and the output TGraphAsymmErrors
1290   fgYieldCorr = new TGraphAsymmErrors(nbins+1);
1291   if (fAsymUncertainties){
1292     fgYieldCorrExtreme = new TGraphAsymmErrors(nbins+1);
1293     fgYieldCorrConservative = new TGraphAsymmErrors(nbins+1);
1294   }
1295
1296   //
1297   // Do the calculation
1298   //
1299   for (Int_t ibin=1; ibin<=nbins; ibin++) {
1300
1301     // Variables initialization
1302     value = 0.; errvalue = 0.; errvalueMax= 0.; errvalueMin= 0.;
1303     valueExtremeMax= 0.; valueExtremeMin= 0.;
1304     valueConservativeMax= 0.; valueConservativeMin= 0.;
1305
1306
1307     // calculate the value
1308     //    physics = reco * fc / bin-width
1309     value = (fhRECpt->GetBinContent(ibin) && fhFc->GetBinContent(ibin)) ?
1310       fhRECpt->GetBinContent(ibin) * fhFc->GetBinContent(ibin) : 0. ;
1311     value /= fhRECpt->GetBinWidth(ibin) ;
1312
1313     // Statistical uncertainty
1314     //    (stat) delta_physics = physics * sqrt ( (delta_reco/reco)^2 )
1315     errvalue = (value!=0. && fhRECpt->GetBinContent(ibin) && fhRECpt->GetBinContent(ibin)!=0.) ?
1316       value * (fhRECpt->GetBinError(ibin)/fhRECpt->GetBinContent(ibin)) : 0. ;
1317
1318     // Calculate the systematic uncertainties
1319     //    (syst but feed-down) delta_physics = physics * sqrt ( (delta_reco_syst/reco)^2 )
1320     //        (feed-down syst) delta_physics = physics * sqrt ( (delta_fc/fc)^2 )
1321     //
1322     // Protect against null denominator. If so, define uncertainty as null
1323     if (fhRECpt->GetBinContent(ibin) && fhRECpt->GetBinContent(ibin)!=0.) {
1324
1325       if (fAsymUncertainties) {
1326
1327         // Systematics but feed-down
1328         if (fgRECSystematics) {
1329           errvalueMax = value * ( fgRECSystematics->GetErrorYhigh(ibin) / fhRECpt->GetBinContent(ibin) );
1330           errvalueMin = value * ( fgRECSystematics->GetErrorYlow(ibin) / fhRECpt->GetBinContent(ibin) );
1331         }
1332         else { errvalueMax = 0.; errvalueMin = 0.; }
1333
1334         // Extreme feed-down systematics
1335         valueExtremeMax = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) + fgFcExtreme->GetErrorYhigh(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
1336         valueExtremeMin = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) - fgFcExtreme->GetErrorYlow(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
1337
1338         // Conservative feed-down systematics
1339         valueConservativeMax = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) + fgFcConservative->GetErrorYhigh(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
1340         valueConservativeMin = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) - fgFcConservative->GetErrorYlow(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
1341
1342       }
1343
1344     }
1345     else { errvalueMax = 0.; errvalueMin = 0.; }
1346
1347     //
1348     // Fill in the histograms
1349     //
1350     fhYieldCorr->SetBinContent(ibin,value);
1351     fhYieldCorr->SetBinError(ibin,errvalue);
1352     //
1353     // Fill the histos and ntuple vs the Eloss hypothesis
1354     //
1355     if (fPbPbElossHypothesis) {
1356       // Loop over the Eloss hypothesis
1357       for (Float_t rval=0.0025; rval<4.0; rval+=0.005){
1358         Int_t rbin = FindTH2YBin(fhYieldCorrRcb,rval);
1359         Double_t fcRcbvalue = fhFcRcb->GetBinContent(ibin,rbin);
1360         //    physics = reco * fcRcb / bin-width
1361         Double_t Rcbvalue = (fhRECpt->GetBinContent(ibin) && fcRcbvalue) ?
1362           fhRECpt->GetBinContent(ibin) * fcRcbvalue : 0. ;
1363         Rcbvalue /= fhRECpt->GetBinWidth(ibin) ;
1364         fhYieldCorrRcb->Fill( fhYieldCorr->GetBinCenter(ibin) , rval, Rcbvalue );
1365         //        cout << " pt "<< fhRECpt->GetBinCenter(ibin) <<" bin "<< ibin<<" rval="<<rval<<", rbin="<<rbin<<" fc-fcRbvalue="<<fcRcbvalue<<", yield="<<Rcbvalue<<endl;
1366       }
1367     }
1368     if (fAsymUncertainties) {
1369       Double_t center = fhYieldCorr->GetBinCenter(ibin);
1370       fgYieldCorr->SetPoint(ibin,center,value); // i,x,y
1371       fgYieldCorr->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueMin,errvalueMax); // i,xl,xh,yl,yh
1372       fhYieldCorrMax->SetBinContent(ibin,value+errvalueMax);
1373       fhYieldCorrMin->SetBinContent(ibin,value-errvalueMin);
1374       fgYieldCorrExtreme->SetPoint(ibin,center,value);
1375       fgYieldCorrExtreme->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),value-valueExtremeMin,valueExtremeMax-value);
1376       fgYieldCorrConservative->SetPoint(ibin,center,value);
1377       fgYieldCorrConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),value-valueConservativeMin,valueConservativeMax-value);
1378     }
1379
1380   }
1381   delete [] binwidths;
1382   delete [] limits;
1383
1384 }
1385
1386 //_________________________________________________________________________________________________________
1387 void AliHFPtSpectrum::CalculateFeedDownCorrectedSpectrumNb(Double_t deltaY, Double_t branchingRatioBintoFinalDecay) {
1388   //
1389   // Compute the feed-down corrected spectrum if feed-down correction is done via Nb (bin by bin)
1390   //    physics =  [ reco  - (lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th) ] / bin-width
1391   //
1392   //    uncertainty:   (stat)  delta_physics = sqrt ( (delta_reco)^2 )  / bin-width
1393   //     (syst but feed-down)  delta_physics = sqrt ( (delta_reco_syst)^2 )  / bin-width
1394   //         (feed-down syst)  delta_physics = sqrt ( (k*delta_lumi/lumi)^2 + (k*delta_eff_trig/eff_trig)^2
1395   //                                                   + (k*delta_Nb/Nb)^2 + (k*delta_eff/eff)^2  + (k*global_eff_ratio)^2 ) / bin-width
1396   //                    where k = lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th
1397   //
1398   //  In addition, in HIC the feed-down correction varies with an energy loss hypothesis: Raa(b-->D) = Rb
1399   //    physics =  [ reco  - ( Tab * Nevt * delta_y * BR_b * eff_trig * eff_b * Nb_th * Rb ) ] / bin-width
1400   //
1401   AliInfo("Calculating the feed-down correction factor and spectrum (Nb method)");
1402
1403   Int_t nbins = fhRECpt->GetNbinsX();
1404   Double_t binwidth = fhRECpt->GetBinWidth(1);
1405   Double_t value = 0., errvalue = 0., errvalueMax = 0., errvalueMin = 0., kfactor = 0.;
1406   Double_t errvalueExtremeMax = 0., errvalueExtremeMin = 0.;
1407   Double_t *limits = new Double_t[nbins+1];
1408   Double_t *binwidths = new Double_t[nbins];
1409   Double_t xlow=0.;
1410   for (Int_t i=1; i<=nbins; i++) {
1411     binwidth = fhRECpt->GetBinWidth(i);
1412     xlow = fhRECpt->GetBinLowEdge(i);
1413     limits[i-1] = xlow;
1414     binwidths[i-1] = binwidth;
1415   }
1416   limits[nbins] = xlow + binwidth;
1417
1418   // declare the output histograms
1419   fhYieldCorr = new TH1D("fhYieldCorr","corrected yield (by Nb)",nbins,limits);
1420   fhYieldCorrMax = new TH1D("fhYieldCorrMax","max corrected yield (by Nb)",nbins,limits);
1421   fhYieldCorrMin = new TH1D("fhYieldCorrMin","min corrected yield (by Nb)",nbins,limits);
1422   if(fPbPbElossHypothesis) {
1423     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.);
1424     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.);
1425     fnHypothesis = new TNtuple("fnHypothesis"," Feed-down correction vs hypothesis (Nb)","pt:Rb:fc:fcMax:fcMin");
1426   }
1427   // and the output TGraphAsymmErrors
1428   fgYieldCorr = new TGraphAsymmErrors(nbins+1);
1429   if (fAsymUncertainties){
1430     fgYieldCorrExtreme = new TGraphAsymmErrors(nbins+1);
1431     fgYieldCorrConservative = new TGraphAsymmErrors(nbins+1);
1432     // Define fc-conservative
1433     fgFcConservative = new TGraphAsymmErrors(nbins+1);
1434     AliInfo(" Beware the conservative & extreme uncertainties are equal by definition !");
1435   }
1436
1437   // variables to define fc-conservative
1438   double correction=0, correctionMax=0., correctionMin=0.;
1439
1440   fhStatUncEffcFD = new TH1D("fhStatUncEffcFD","direct charm stat unc on the feed-down correction",nbins,limits);
1441   fhStatUncEffbFD = new TH1D("fhStatUncEffbFD","secondary charm stat unc on the feed-down correction",nbins,limits);
1442   Double_t correctionUncStatEffc=0.;
1443   Double_t correctionUncStatEffb=0.;
1444
1445
1446   //
1447   // Do the calculation
1448   //
1449   for (Int_t ibin=1; ibin<=nbins; ibin++) {
1450
1451     // Calculate the value
1452     //    physics =  [ reco  - (lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th) ] / bin-width
1453     //  In HIC :   physics =  [ reco  - ( Tab * Nevt * delta_y * BR_b * eff_trig * eff_b * Nb_th * Rb ) ] / bin-width
1454     //
1455     //
1456     Double_t frac = 1.0, errfrac =0.;
1457
1458     // Variables initialization
1459     value = 0.; errvalue = 0.; errvalueMax = 0.; errvalueMin = 0.; kfactor = 0.;
1460     errvalueExtremeMax = 0.; errvalueExtremeMin = 0.;
1461     correction=0; correctionMax=0.; correctionMin=0.;
1462     correctionUncStatEffc=0.; correctionUncStatEffb=0.;
1463
1464     if(fPbPbElossHypothesis) {
1465       frac = fTab[0]*fNevts;
1466       if(fIsEventPlane) frac = frac/2.0;
1467       errfrac = frac * TMath::Sqrt( (fTab[1]/fTab[0])*(fTab[1]/fTab[0]) + (1/fNevts) );
1468     } else {
1469       frac = fLuminosity[0];
1470       errfrac = fLuminosity[1];
1471     }
1472
1473     value = ( fhRECpt->GetBinContent(ibin)>0. && fhRECpt->GetBinContent(ibin)!=0. &&
1474               fhFeedDownMCpt->GetBinContent(ibin)>0. && fhFeedDownEffpt->GetBinContent(ibin)>0. ) ?
1475       fhRECpt->GetBinContent(ibin) - frac*(deltaY*branchingRatioBintoFinalDecay*fParticleAntiParticle*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) * fhRECpt->GetBinWidth(ibin) )
1476       : 0. ;
1477     value /= fhRECpt->GetBinWidth(ibin);
1478     if (value<0.) value =0.;
1479
1480     //  Statistical uncertainty:   delta_physics = sqrt ( (delta_reco)^2 )  / bin-width
1481     errvalue = (value!=0. && fhRECpt->GetBinError(ibin) && fhRECpt->GetBinError(ibin)!=0.)  ?
1482       fhRECpt->GetBinError(ibin) : 0.;
1483     errvalue /= fhRECpt->GetBinWidth(ibin);
1484
1485     // Correction (fc) : Estimate of the relative amount feed-down subtracted
1486     // correction =  [ 1  - (lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th)/reco ]
1487     // in HIC: correction =  [ 1  - ( Tab * Nevt * delta_y * BR_b * eff_trig * eff_b * Nb_th)/reco ]
1488     correction = (value>0.) ?
1489       1 - (frac*deltaY*branchingRatioBintoFinalDecay*fParticleAntiParticle*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) * fhRECpt->GetBinWidth(ibin) ) / fhRECpt->GetBinContent(ibin)  : 0. ;
1490     if (correction<0.) correction = 0.;
1491
1492     //    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;
1493
1494     // Systematic uncertainties
1495     //     (syst but feed-down)  delta_physics = sqrt ( (delta_reco_syst)^2 )  / bin-width
1496     //         (feed-down syst)  delta_physics = sqrt ( (k*delta_lumi/lumi)^2 + (k*delta_eff_trig/eff_trig)^2
1497     //                                                   + (k*delta_Nb/Nb)^2 + (k*delta_eff/eff)^2 + (k*global_eff_ratio)^2 ) / bin-width
1498     //                    where k = lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th * bin-width
1499     kfactor = frac*deltaY*branchingRatioBintoFinalDecay*fParticleAntiParticle*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) * fhRECpt->GetBinWidth(ibin) ;
1500     //
1501     if (fAsymUncertainties && value>0.) {
1502       Double_t nb =  fhFeedDownMCpt->GetBinContent(ibin);
1503       Double_t nbDmax = fhFeedDownMCptMax->GetBinContent(ibin) - fhFeedDownMCpt->GetBinContent(ibin);
1504       Double_t nbDmin = fhFeedDownMCpt->GetBinContent(ibin) - fhFeedDownMCptMin->GetBinContent(ibin);
1505
1506       // Systematics but feed-down
1507       if (fgRECSystematics){
1508         errvalueMax = fgRECSystematics->GetErrorYhigh(ibin) / fhRECpt->GetBinWidth(ibin) ;
1509         errvalueMin = fgRECSystematics->GetErrorYlow(ibin) / fhRECpt->GetBinWidth(ibin);
1510       }
1511       else { errvalueMax = 0.; errvalueMin = 0.; }
1512
1513       // Feed-down systematics
1514       // min value with the maximum Nb
1515       Double_t errCom =  ( (kfactor*errfrac/frac)*(kfactor*errfrac/frac) ) +
1516         ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) +
1517         ( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) ) +
1518         ( (kfactor*fGlobalEfficiencyUncertainties[1])*(kfactor*fGlobalEfficiencyUncertainties[1]) ) ;
1519       errvalueExtremeMin = TMath::Sqrt( errCom + ( (kfactor*nbDmax/nb)*(kfactor*nbDmax/nb) ) ) / fhRECpt->GetBinWidth(ibin);
1520       // max value with the minimum Nb
1521       errvalueExtremeMax =  TMath::Sqrt( errCom + ( (kfactor*nbDmin/nb)*(kfactor*nbDmin/nb) ) ) / fhRECpt->GetBinWidth(ibin);
1522
1523       // Correction systematics (fc)
1524       // min value with the maximum Nb
1525       correctionMin = TMath::Sqrt( errCom + ( (kfactor*nbDmax/nb)*(kfactor*nbDmax/nb) ) ) / fhRECpt->GetBinContent(ibin) ;
1526       // max value with the minimum Nb
1527       correctionMax =  TMath::Sqrt( errCom + ( (kfactor*nbDmin/nb)*(kfactor*nbDmin/nb) ) ) / fhRECpt->GetBinContent(ibin) ;
1528       //
1529       correctionUncStatEffb = TMath::Sqrt(  ( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))       )
1530                                             ) / fhRECpt->GetBinContent(ibin) ;
1531       correctionUncStatEffc = 0.;
1532     }
1533     else{ // Don't consider Nb uncertainty in this case [ to be tested!!! ]
1534       errvalueExtremeMax =  TMath::Sqrt( ( (kfactor*errfrac/frac)*(kfactor*errfrac/frac) ) +
1535                                          ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) )  +
1536                                          ( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))  )  +
1537                                          ( (kfactor*fGlobalEfficiencyUncertainties[1])*(kfactor*fGlobalEfficiencyUncertainties[1]) )
1538                                          ) / fhRECpt->GetBinWidth(ibin);
1539       errvalueExtremeMin =  errvalueExtremeMax ;
1540     }
1541
1542
1543     // fill in histograms
1544     fhYieldCorr->SetBinContent(ibin,value);
1545     fhYieldCorr->SetBinError(ibin,errvalue);
1546     //
1547     // Estimate how the result varies vs charm/beauty Eloss hypothesis
1548     //
1549     if ( correction>1.0e-16 && fPbPbElossHypothesis){
1550       // Loop over the Eloss hypothesis
1551       //      Int_t rbin=0;
1552       for (Float_t rval=0.0025; rval<4.0; rval+=0.005){
1553         // correction =  [ 1  - (Tab *Nevt * delta_y * BR_b * eff_trig * eff_b * Nb_th *binwidth )* (rval) /reco ]
1554         Double_t fcRcbvalue = 1 - (frac*deltaY*branchingRatioBintoFinalDecay*fParticleAntiParticle*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin)*fhRECpt->GetBinWidth(ibin) * rval ) / fhRECpt->GetBinContent(ibin) ;
1555         //      cout << " rval="<<rval <<" , fc="<<fcRcbvalue<<"    ";
1556         if(fcRcbvalue<1.e-16) fcRcbvalue=0.;
1557         fhFcRcb->Fill( fhRECpt->GetBinCenter(ibin) , rval, fcRcbvalue );
1558         //    physics = reco * fcRcb / bin-width
1559         Double_t Rcbvalue = (fhRECpt->GetBinContent(ibin) && fcRcbvalue) ?
1560           fhRECpt->GetBinContent(ibin) * fcRcbvalue : 0. ;
1561         Rcbvalue /= fhRECpt->GetBinWidth(ibin) ;
1562         fhYieldCorrRcb->Fill( fhYieldCorr->GetBinCenter(ibin) , rval, Rcbvalue );
1563         //      if(ibin==3){
1564         //        cout << " pt "<< fhFcRcb->GetBinCenter(ibin) <<" bin "<< ibin<<" rval="<<rval<<", rbin="<<rbin<<", fc-Rb-value="<< fcRcbvalue << ", yield-Rb-value="<< Rcbvalue <<endl;
1565         //      cout << " pt "<< fhFcRcb->GetBinCenter(ibin) <<" bin "<< ibin<<" rval="<<rval<<", fc-Rb-value="<< fcRcbvalue << ", yield-Rb-value="<< Rcbvalue <<endl;
1566         //        rbin++;
1567         //      }
1568         Double_t correctionMaxRcb = correctionMax*rval;
1569         Double_t correctionMinRcb = correctionMin*rval;
1570         fnHypothesis->Fill( fhYieldCorr->GetBinCenter(ibin), rval, fcRcbvalue, fcRcbvalue + correctionMaxRcb, fcRcbvalue - correctionMinRcb);
1571 //      if(ibin==3){
1572 //        cout << " pt="<< fhFcRcb->GetBinCenter(ibin) <<", rval="<<rval<<", fc="<<fcRcbvalue<<" +"<<correctionMaxRcb<<" -"<<correctionMinRcb<<endl;}
1573       }
1574     }
1575     //
1576     // Fill the rest of (asymmetric) histograms
1577     //
1578     if (fAsymUncertainties) {
1579       Double_t x = fhYieldCorr->GetBinCenter(ibin);
1580       fgYieldCorr->SetPoint(ibin,x,value); // i,x,y
1581       fgYieldCorr->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueMin,errvalueMax); // i,xl,xh,yl,yh
1582       fhYieldCorrMax->SetBinContent(ibin,value+errvalueMax);
1583       fhYieldCorrMin->SetBinContent(ibin,value-errvalueMin);
1584       fgYieldCorrExtreme->SetPoint(ibin,x,value); // i,x,y
1585       fgYieldCorrExtreme->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueExtremeMin,errvalueExtremeMax); // i,xl,xh,yl,yh
1586       fgYieldCorrConservative->SetPoint(ibin,x,value); // i,x,y
1587       fgYieldCorrConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueExtremeMin,errvalueExtremeMax); // i,xl,xh,yl,yh
1588       //      cout << " bin " << ibin << ", correction " << correction << ", min correction unc " << correctionMin << ", max correction unc " << correctionMax << endl;
1589       if(correction>0.){
1590         fgFcConservative->SetPoint(ibin,x,correction);
1591         fgFcConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),correctionMin,correctionMax);
1592
1593         fhStatUncEffbFD->SetBinContent(ibin,0.); fhStatUncEffbFD->SetBinError(ibin,correctionUncStatEffb/correction*100.);
1594         fhStatUncEffcFD->SetBinContent(ibin,0.); fhStatUncEffcFD->SetBinError(ibin,correctionUncStatEffc/correction*100.);
1595         //      cout << " pt "<< fhStatUncEffcFD->GetBinCenter(ibin) <<" bin "<< ibin<<" fc-stat-c ="<< correctionUncStatEffc/correction <<" fc-stat-b ="<< correctionUncStatEffb/correction <<endl;
1596       }
1597       else{
1598         fgFcConservative->SetPoint(ibin,x,0.);
1599         fgFcConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),0.,0.);
1600       }
1601     }
1602
1603   }
1604   delete [] binwidths;
1605   delete [] limits;
1606
1607 }
1608
1609
1610 //_________________________________________________________________________________________________________
1611 void AliHFPtSpectrum::ComputeSystUncertainties(AliHFSystErr *systematics, Bool_t combineFeedDown) {
1612   //
1613   // Function that re-calculates the global systematic uncertainties
1614   //   by calling the class AliHFSystErr and combining those
1615   //   (in quadrature) with the feed-down subtraction uncertainties
1616   //
1617
1618   // Estimate the feed-down uncertainty in percentage
1619   Int_t nentries = 0;
1620   TGraphAsymmErrors *grErrFeeddown = 0;
1621   Double_t x=0., y=0., errx=0., erryl=0., erryh=0;
1622   if(fFeedDownOption!=0) {
1623     nentries = fgSigmaCorrConservative->GetN();
1624     grErrFeeddown = new TGraphAsymmErrors(nentries);
1625     for(Int_t i=0; i<nentries; i++) {
1626       x=0.; y=0.; errx=0.; erryl=0.; erryh=0.;
1627       fgSigmaCorrConservative->GetPoint(i,x,y);
1628       if(y>0.){
1629         errx = fgSigmaCorrConservative->GetErrorXlow(i) ;
1630         erryl = fgSigmaCorrConservative->GetErrorYlow(i) / y ;
1631         erryh = fgSigmaCorrConservative->GetErrorYhigh(i) / y ;
1632       }
1633       //    cout << " x "<< x << " +- "<<errx<<" , y "<<y<<" + "<<erryh<<" - "<<erryl<<endl;
1634       grErrFeeddown->SetPoint(i,x,0.);
1635       grErrFeeddown->SetPointError(i,errx,errx,erryl,erryh); //i, xl, xh, yl, yh
1636     }
1637   }
1638
1639   // Draw all the systematics independently
1640   systematics->DrawErrors(grErrFeeddown);
1641
1642   // Set the sigma systematic uncertainties
1643   // possibly combine with the feed-down uncertainties
1644   Double_t errylcomb=0., erryhcomb=0;
1645   for(Int_t i=1; i<nentries; i++) {
1646     fgSigmaCorr->GetPoint(i,x,y);
1647     errx = grErrFeeddown->GetErrorXlow(i) ;
1648     erryl = grErrFeeddown->GetErrorYlow(i);
1649     erryh = grErrFeeddown->GetErrorYhigh(i);
1650     if (combineFeedDown) {
1651       errylcomb = systematics->GetTotalSystErr(x,erryl) * y ;
1652       erryhcomb = systematics->GetTotalSystErr(x,erryh) * y ;
1653     } else {
1654       errylcomb = systematics->GetTotalSystErr(x) * y ;
1655       erryhcomb = systematics->GetTotalSystErr(x) * y ;
1656     }
1657     fgSigmaCorr->SetPointError(i,errx,errx,errylcomb,erryhcomb);
1658     //
1659     fhSigmaCorrDataSyst->SetBinContent(i,y);
1660     erryl = systematics->GetTotalSystErr(x) * y ;
1661     fhSigmaCorrDataSyst->SetBinError(i,erryl);
1662   }
1663
1664 }
1665
1666
1667 //_________________________________________________________________________________________________________
1668 void AliHFPtSpectrum::DrawSpectrum(TGraphAsymmErrors *gPrediction) {
1669   //
1670   // Example method to draw the corrected spectrum & the theoretical prediction
1671   //
1672
1673   TCanvas *csigma = new TCanvas("csigma","Draw the corrected cross-section & the prediction");
1674   csigma->SetFillColor(0);
1675   gPrediction->GetXaxis()->SetTitleSize(0.05);
1676   gPrediction->GetXaxis()->SetTitleOffset(0.95);
1677   gPrediction->GetYaxis()->SetTitleSize(0.05);
1678   gPrediction->GetYaxis()->SetTitleOffset(0.95);
1679   gPrediction->GetXaxis()->SetTitle("p_{T}  [GeV]");
1680   gPrediction->GetYaxis()->SetTitle("BR #times #frac{d#sigma}{dp_{T}} |_{|y|<0.5}   [pb/GeV]");
1681   gPrediction->SetLineColor(kGreen+2);
1682   gPrediction->SetLineWidth(3);
1683   gPrediction->SetFillColor(kGreen+1);
1684   gPrediction->Draw("3CA");
1685   fgSigmaCorr->SetLineColor(kRed);
1686   fgSigmaCorr->SetLineWidth(1);
1687   fgSigmaCorr->SetFillColor(kRed);
1688   fgSigmaCorr->SetFillStyle(0);
1689   fgSigmaCorr->Draw("2");
1690   fhSigmaCorr->SetMarkerColor(kRed);
1691   fhSigmaCorr->Draw("esame");
1692   csigma->SetLogy();
1693   TLegend * leg = new TLegend(0.7,0.75,0.87,0.5);
1694   leg->SetBorderSize(0);
1695   leg->SetLineColor(0);
1696   leg->SetFillColor(0);
1697   leg->SetTextFont(42);
1698   leg->AddEntry(gPrediction,"FONLL ","fl");
1699   leg->AddEntry(fhSigmaCorr,"data stat. unc.","pl");
1700   leg->AddEntry(fgSigmaCorr,"data syst. unc.","f");
1701   leg->Draw();
1702   csigma->Draw();
1703
1704 }
1705
1706 //_________________________________________________________________________________________________________
1707 TH1D * AliHFPtSpectrum::ReweightHisto(TH1D *hToReweight, TH1D *hReference){
1708   //
1709   // Function to  reweight histograms for testing purposes:
1710   // This function takes the histo hToReweight and reweights
1711   //  it (its pt shape) with respect to hReference
1712   //
1713
1714   // check histograms consistency
1715   Bool_t areconsistent=kTRUE;
1716   areconsistent &= CheckHistosConsistency(hToReweight,hReference);
1717   if (!areconsistent) {
1718     AliInfo("the histograms to reweight are not consistent (bin width, bounds)");
1719     return NULL;
1720   }
1721
1722   // define a new empty histogram
1723   TH1D *hReweighted = (TH1D*)hToReweight->Clone("hReweighted");
1724   hReweighted->Reset();
1725   Double_t weight=1.0;
1726   Double_t yvalue=1.0;
1727   Double_t integralRef = hReference->Integral();
1728   Double_t integralH = hToReweight->Integral();
1729
1730   // now reweight the spectra
1731   //
1732   // the weight is the relative probability of the given pt bin in the reference histo
1733   //  divided by its relative probability (to normalize it) on the histo to re-weight
1734   for (Int_t i=0; i<=hToReweight->GetNbinsX(); i++) {
1735     weight = (hReference->GetBinContent(i)/integralRef) / (hToReweight->GetBinContent(i)/integralH) ;
1736     yvalue = hToReweight->GetBinContent(i);
1737     hReweighted->SetBinContent(i,yvalue*weight);
1738   }
1739
1740   return (TH1D*)hReweighted;
1741 }
1742
1743 //_________________________________________________________________________________________________________
1744 TH1D * AliHFPtSpectrum::ReweightRecHisto(TH1D *hRecToReweight, TH1D *hMCToReweight, TH1D *hMCReference){
1745   //
1746   // Function to  reweight histograms for testing purposes:
1747   // This function takes the histo hToReweight and reweights
1748   //  it (its pt shape) with respect to hReference /hMCToReweight
1749   //
1750
1751   // check histograms consistency
1752   Bool_t areconsistent=kTRUE;
1753   areconsistent &= CheckHistosConsistency(hMCToReweight,hMCReference);
1754   areconsistent &= CheckHistosConsistency(hRecToReweight,hMCReference);
1755   if (!areconsistent) {
1756     AliInfo("the histograms to reweight are not consistent (bin width, bounds)");
1757     return NULL;
1758   }
1759
1760   // define a new empty histogram
1761   TH1D *hReweighted = (TH1D*)hMCToReweight->Clone("hReweighted");
1762   hReweighted->Reset();
1763   TH1D *hRecReweighted = (TH1D*)hRecToReweight->Clone("hRecReweighted");
1764   hRecReweighted->Reset();
1765   Double_t weight=1.0;
1766   Double_t yvalue=1.0, yrecvalue=1.0;
1767   Double_t integralRef = hMCReference->Integral();
1768   Double_t integralH = hMCToReweight->Integral();
1769
1770   // now reweight the spectra
1771   //
1772   // the weight is the relative probability of the given pt bin
1773   //  that should be applied in the MC histo to get the reference histo shape
1774   //  Probabilities are properly normalized.
1775   for (Int_t i=0; i<=hMCToReweight->GetNbinsX(); i++) {
1776     weight = (hMCReference->GetBinContent(i)/integralRef) / (hMCToReweight->GetBinContent(i)/integralH) ;
1777     yvalue = hMCToReweight->GetBinContent(i);
1778     hReweighted->SetBinContent(i,yvalue*weight);
1779     yrecvalue = hRecToReweight->GetBinContent(i);
1780     hRecReweighted->SetBinContent(i,yrecvalue*weight);
1781   }
1782
1783   return (TH1D*)hRecReweighted;
1784 }
1785
1786
1787
1788 //_________________________________________________________________________________________________________
1789 Int_t AliHFPtSpectrum::FindTH2YBin(TH2D *histo, Float_t yvalue){
1790   //
1791   // Function to find the y-axis bin of a TH2 for a given y-value
1792   //
1793
1794   Int_t nbins = histo->GetNbinsY();
1795   Int_t ybin=0;
1796   for (int j=0; j<=nbins; j++) {
1797     Float_t value = histo->GetYaxis()->GetBinCenter(j);
1798     Float_t width = histo->GetYaxis()->GetBinWidth(j);
1799     //    if( TMath::Abs(yvalue-value)<= width/2. ) {
1800     if( TMath::Abs(yvalue-value)<= width ) {
1801       ybin =j;
1802       //      cout <<" value "<<value << ", yval "<< yvalue<<", bin width "<<width/2.<< " y ="<<ybin<<endl;
1803       break;
1804     }
1805   }
1806
1807   return ybin;
1808 }
1809
1810 //_________________________________________________________________________________________________________
1811 void AliHFPtSpectrum::ResetStatUncEff(){
1812
1813   Int_t entries = fhDirectEffpt->GetNbinsX();
1814   for(Int_t i=0; i<=entries; i++){
1815     fhDirectEffpt->SetBinError(i,0.);
1816     fhFeedDownEffpt->SetBinError(i,0.);
1817   }
1818
1819 }