NuMMRun.cxx

Go to the documentation of this file.

#include <cassert>
#include <cmath>

#include "Rtypes.h"
#include "TH1D.h"

#include "Conventions/Detector.h"
#include "Conventions/SimFlag.h"
#include "MessageService/MsgService.h"
#include "NtupleUtils/NuCuts.h"
#include "NtupleUtils/NuHistInterpolator.h"
#include "NtupleUtils/NuMMHelperPRL.h"
#include "NtupleUtils/NuMMRun.h"

ClassImp(NuMMRun)

CVSID("$Id: NuMMRun.cxx,v 1.16 2010/01/17 18:15:45 bckhouse Exp $");

//____________________________________________________________________72
NuMMRun::NuMMRun()
  : fFakeDataStudy(false),
    fQuietMode(false),
    fDisappearanceModel(0),
    fConstrainEnergy(false),
    fMaximumEnergy(0.0)
//    fMaximumBin(0)
{
  NuCuts nuCuts;
  fNDFidMass = nuCuts.FiducialMass(Detector::kNear,
                                   SimFlag::kData,
                                   NuCuts::kJJE1);
  fFDFidMass = nuCuts.FiducialMass(Detector::kFar,
                                   SimFlag::kData,
                                   NuCuts::kJJE1);
}

//____________________________________________________________________72
NuMMRun::NuMMRun(NuMMHelperPRL* helper,
                 NuHistInterpolator* ndData,
                 NuHistInterpolator* fdData)
  : fFakeDataStudy(false),
    fQuietMode(false),
    fDisappearanceModel(0)
{
  NuCuts nuCuts;
  fNDFidMass = nuCuts.FiducialMass(Detector::kNear,
                                   SimFlag::kData,
                                   NuCuts::kJJE1);
  fFDFidMass = nuCuts.FiducialMass(Detector::kFar,
                                   SimFlag::kData,
                                   NuCuts::kJJE1);
  fHelper = helper;
  fNDData = ndData;
  fFDData = fdData;
}

//____________________________________________________________________72
void NuMMRun::ConfigForDisappearanceModel(int model)
{
  assert(model == 0 || model == 1 || model == 2);
  fDisappearanceModel = model;
}

//____________________________________________________________________72
const Double_t NuMMRun::StatsLikelihood(const TH1D* fdPred,
                                        const TH1D* fdData) const
{
/*
  Bool_t    fConstrainEnergy;
  Double_t  fMaximumEnergy;
  Int_t     fMaximumBin;*/
  Int_t MaximumBin = 0;

  // Are we constraining?
  if (fConstrainEnergy) {
    // Calculate the upper bin (not sure why using findfixbin... is this old?)
    MaximumBin = fdPred->GetXaxis()->FindFixBin(fMaximumEnergy);
    MAXMSG("NuMMRun", 5, Msg::kDebug) << "    Calculated highest bin as "
      << MaximumBin << endl;
  } else {
    // In this case, we are not constraining - use the full range!
    MAXMSG("NuMMRun", 5, Msg::kDebug) << "    Not constraining bin comparison" << endl;
    MaximumBin = fdPred->GetNbinsX();
  }

  Double_t like = 0;

  //Aim to minimise -2lnL. That's what I'm returning.
  for (Int_t i=1; i<=MaximumBin; ++i){
    //Bizarre limits because root histograms are silly
    Double_t mnu = fdPred->GetBinContent(i);
    Double_t dnu = fdData->GetBinContent(i);

    // Replace the previous Log-likelihood with one that is a bit
    // more robust for fringe cases.
    like += NuUtilities::LogLikelihood(dnu,mnu);
  }
  return like;
}

---