PWG3/vertexingHF/AliHFPtSpectrum.cxx

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