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