NuHistInterpolator.cxx

Go to the documentation of this file.

#include <cassert>

#include "TF1.h"
#include "TGraph.h"

#include "MessageService/MsgService.h"
#include "NtupleUtils/NuHistInterpolator.h"

ClassImp(NuHistInterpolator)

CVSID("$Id: NuHistInterpolator.cxx,v 1.5 2009/04/08 20:16:46 gmieg Exp $");

//____________________________________________________________________72
void NuHistInterpolator::push_back(const Double_t shift,
                                   const TH1D& spectrum,
                                   const Double_t pot)
{
  fSpectra[shift] = spectrum;
  fPoT = pot;
}

//____________________________________________________________________72
void NuHistInterpolator::push_back(const Double_t shift,
                                   const TH1D& spectrum)
{
  fSpectra[shift] = spectrum;
}
/*
//____________________________________________________________________72
const TH1D NuHistInterpolator::ProvideSpectrum(const Double_t shift) const
{
  std::map<Double_t,TH1D>::const_iterator spectraIt = fSpectra.begin();
  while (spectraIt->first < shift && spectraIt != fSpectra.end()){
    ++spectraIt;
  }
  if (fSpectra.begin() == spectraIt){
    MSG("NuHistInterpolator",Msg::kWarning)
      << "Systematic shift " << shift << " below provided range" << endl;
    assert(false);
  }
  if (fSpectra.end() == spectraIt){
    MSG("NuHistInterpolator",Msg::kWarning)
      << "Systematic shift " << shift << " above provided range" << endl;
    assert(false);
  }

  TH1D spectrumAbove(spectraIt->second);
  const Double_t shiftAbove = spectraIt->first;

  ++spectraIt;
  if (fSpectra.end() == spectraIt){
    cout << "Dying" << endl;
    assert(false);
}
  TH1D spectrumTwoAbove(spectraIt->second);
  const Double_t shiftTwoAbove = spectraIt->first;

  --spectraIt;
  --spectraIt;
  TH1D spectrumBelow(spectraIt->second);
  const Double_t shiftBelow = spectraIt->first;

  if (fSpectra.begin() == spectraIt){
    cout << "Dying" << endl;
    assert(false);
  }

  --spectraIt;
  if (fSpectra.begin() == spectraIt){
    cout << "Dying" << endl;
    assert(false);
  }

  TH1D spectrumTwoBelow(spectraIt->second);
  const Double_t shiftTwoBelow = spectraIt->first;

  TH1D middleGradient(spectrumAbove);
  middleGradient.Add(&spectrumBelow,-1.0);
  middleGradient.Scale(1.0/(shiftAbove-shiftBelow));

  TH1D lowGradient(spectrumBelow);
  lowGradient.Add(&spectrumTwoBelow,-1.0);
  lowGradient.Scale(1.0/(shiftBelow-shiftTwoBelow));

  TH1D highGradient(spectrumTwoAbove);
  highGradient.Add(&spectrumAbove,-1.0);
  highGradient.Scale(1.0/(shiftTwoAbove-shiftAbove));

  if (shift<(shiftAbove+shiftBelow)/2.0){
    TH1D spectrumOut(lowGradient);
    spectrumOut.Scale(shift - shiftTwoBelow);
    spectrumOut.Add(&spectrumTwoBelow);
    for (Int_t bin=1; bin<=spectrumOut.GetNbinsX();++bin){
      spectrumOut.SetBinError(bin,1.0);
      if (spectrumOut.GetBinContent(bin)<0){
        spectrumOut.SetBinContent(bin,0);
      }
    }
    return spectrumOut;
  }
  else{
    TH1D spectrumOut(highGradient);
    spectrumOut.Scale(shift - shiftAbove);
    spectrumOut.Add(&spectrumAbove);
    for (Int_t bin=1; bin<=spectrumOut.GetNbinsX();++bin){
      if (spectrumOut.GetBinContent(bin)<0){
        spectrumOut.SetBinContent(bin,0);
      }
      spectrumOut.SetBinError(bin,1.0);
    }
    return spectrumOut;
  }
}
*/

//____________________________________________________________________72
const TH1D NuHistInterpolator::ProvideSpectrum(const Double_t shift) const
{
  //This is the simple interpolator
  std::map<Double_t,TH1D>::const_iterator spectraIt = fSpectra.begin();
  while (spectraIt->first < shift && spectraIt != fSpectra.end()){
    ++spectraIt;
  }
  if (fSpectra.begin() == spectraIt){
    MSG("NuHistInterpolator",Msg::kWarning)
      << "Systematic shift " << shift << " below provided range" << endl;
    assert(false);
  }
  if (fSpectra.end() == spectraIt){
    MSG("NuHistInterpolator",Msg::kWarning)
      << "Systematic shift " << shift << " above provided range" << endl;
    assert(false);
  }

  TH1D spectrumAbove(spectraIt->second);
  const Double_t shiftAbove = spectraIt->first;
  --spectraIt;
  TH1D spectrumBelow(spectraIt->second);
  const Double_t shiftBelow = spectraIt->first;
//   cout << "Shift above = " << shiftAbove
//        << "; shift below = " << shiftBelow << endl;

  if (shiftAbove == shiftBelow){
    MSG("NuHistInterpolator",Msg::kWarning)
      << "Two shifted spectra have the identical systematic shifts" << endl;
    return spectrumAbove;
  }

  spectrumAbove.Scale(shift-shiftBelow);
  spectrumAbove.Add(&spectrumBelow,shiftAbove-shift);
  spectrumAbove.Scale(1.0/(shiftAbove-shiftBelow));
//    static Int_t count=0;
//     if (!(count%6)){
//    cout << "Being asked for spectrum with shift " << shift << endl;
//   cout << "Bin 10 " << spectrumAbove.GetBinContent(10) << endl;
//     }
//     ++count;
  return spectrumAbove;
}

/*
//____________________________________________________________________72
const TH1D NuHistInterpolator::ProvideSpectrum(const Double_t shift) const
{
  std::map<Double_t,TH1D>::const_iterator spectraIt = fSpectra.begin();
  while (spectraIt->first < shift && spectraIt != fSpectra.end()){
    ++spectraIt;
  }
  if (fSpectra.begin() == spectraIt){
    MSG("NuHistInterpolator",Msg::kWarning)
      << "Systematic shift " << shift << " below provided range" << endl;
    assert(false);
  }
  if (fSpectra.end() == spectraIt){
    MSG("NuHistInterpolator",Msg::kWarning)
      << "Systematic shift " << shift << " above provided range" << endl;
    assert(false);
  }
  --spectraIt;
  if (fSpectra.begin() == spectraIt){
    MSG("NuHistInterpolator",Msg::kWarning)
      << "Systematic shift " << shift << " below provided range" << endl;
    assert(false);
  }
  --spectraIt;
  TH1D outSpectrum(spectraIt->second);
//   outSpectrum.Clear();
//   outSpectrum.Reset();
  for (Int_t bin = 1; bin<=spectraIt->second.GetNbinsX(); ++bin){
    TGraph gInterpolate(4);
    gInterpolate.SetPoint(0,spectraIt->first,spectraIt->second.GetBinContent(bin));
    ++spectraIt;
    gInterpolate.SetPoint(1,spectraIt->first,spectraIt->second.GetBinContent(bin));
    ++spectraIt;
    gInterpolate.SetPoint(2,spectraIt->first,spectraIt->second.GetBinContent(bin));
    ++spectraIt;
    gInterpolate.SetPoint(3,spectraIt->first,spectraIt->second.GetBinContent(bin));
    gInterpolate.Fit("pol2","Q");
     Double_t par0 = gInterpolate.GetFunction("pol2")->GetParameter(0);
     Double_t par1 = gInterpolate.GetFunction("pol2")->GetParameter(1);
     Double_t par2 = gInterpolate.GetFunction("pol2")->GetParameter(2);
     Double_t binContent = par0 + shift*par1 + shift*shift*par2;
     if (binContent<0.0){binContent=0.0;}
     outSpectrum.SetBinContent(bin,binContent);
    --spectraIt;
    --spectraIt;
    --spectraIt;
  }
  return outSpectrum;
}
*/
/*
//____________________________________________________________________72
const TH1D NuHistInterpolator::ProvideSpectrum(const Double_t shift) const
{
//This is the good interpolator
  std::map<Double_t,TH1D>::const_iterator spectraIt = fSpectra.begin();
  TH1D outSpectrum(spectraIt->second);
  for (Int_t bin = 1; bin<=outSpectrum.GetNbinsX(); ++bin){
    TGraph gInterpolate(fSpectra.size());
    Int_t count=0;
    for (spectraIt = fSpectra.begin();
         spectraIt != fSpectra.end();
         ++spectraIt){
      gInterpolate.SetPoint(count,
                            spectraIt->first,
                            spectraIt->second.GetBinContent(bin));
      ++count;
//       cout << "Limit: " << outSpectrum.GetNbinsX() << endl;
//        cout << "Bin: " << bin << "Count: " << count << endl;
    }
    gInterpolate.Fit("pol2","Q");
    Double_t par0 = gInterpolate.GetFunction("pol2")->GetParameter(0);
    Double_t par1 = gInterpolate.GetFunction("pol2")->GetParameter(1);
    Double_t par2 = gInterpolate.GetFunction("pol2")->GetParameter(2);
    Double_t binContent = par0 + shift*par1 + shift*shift*par2;

//     gInterpolate.Fit("pol6","Q");
//     Double_t par0 = gInterpolate.GetFunction("pol6")->GetParameter(0);
//     Double_t par1 = gInterpolate.GetFunction("pol6")->GetParameter(1);
//     Double_t par2 = gInterpolate.GetFunction("pol6")->GetParameter(2);
//     Double_t par3 = gInterpolate.GetFunction("pol6")->GetParameter(3);
//     Double_t par4 = gInterpolate.GetFunction("pol6")->GetParameter(4);
//     Double_t par5 = gInterpolate.GetFunction("pol6")->GetParameter(5);
//     Double_t par6 = gInterpolate.GetFunction("pol6")->GetParameter(6);
//     Double_t binContent = par0 + shift*par1 + shift*shift*par2 + shift*shift*shift*par3 + shift*shift*shift*shift*par4
//       + shift*shift*shift*shift*shift*par5 + shift*shift*shift*shift*shift*shift*par6;

    if (binContent<0.0){binContent=0.0;}
    outSpectrum.SetBinContent(bin,binContent);
  }
  return outSpectrum;
}
*/
//____________________________________________________________________72
const TH1D NuHistInterpolator::ApplyShiftTo(const TH1D& spectrumToShift,
                                            const Double_t shift) const
{
  TH1D outSpectrum(spectrumToShift);
  string name = spectrumToShift.GetName();
  string title = spectrumToShift.GetTitle();
//   outSpectrum.Reset();
//   outSpectrum.Clear();

  std::map<Double_t,TH1D>::const_iterator spectraIt = fSpectra.begin();
  for (Int_t bin = 1; bin<=outSpectrum.GetNbinsX(); ++bin){
    TGraph gInterpolate(fSpectra.size());
    Int_t count=0;
    for (spectraIt = fSpectra.begin();
         spectraIt != fSpectra.end();
         ++spectraIt){
      gInterpolate.SetPoint(count,
                            spectraIt->first,
                            spectraIt->second.GetBinContent(bin));
      ++count;
//       cout << "Limit: " << outSpectrum.GetNbinsX() << endl;
//        cout << "Bin: " << bin << "Count: " << count << endl;
    }
    gInterpolate.Fit("pol2","Q");
    Double_t par0 = gInterpolate.GetFunction("pol2")->GetParameter(0);
    Double_t par1 = gInterpolate.GetFunction("pol2")->GetParameter(1);
    Double_t par2 = gInterpolate.GetFunction("pol2")->GetParameter(2);

    Double_t scaleFactor = par0 + shift*par1 + shift*shift*par2;
    scaleFactor /= par0;
    Double_t binContent = spectrumToShift.GetBinContent(bin);
//     if (10==bin){cout << "Scale factor: " << scaleFactor << endl;}
    binContent *= scaleFactor;
    if (binContent<0.0){binContent=0.0;}
    outSpectrum.SetBinContent(bin,binContent);
  }
//   outSpectrum.SetName(name.c_str());
//   outSpectrum.SetTitle(title.c_str());
  return outSpectrum;
}

---