FluxWeightSystematics.cxx

#include "FluxWeightSystematics.hxx"
#include "ND280AnalysisUtils.hxx"
#include "ToyBoxND280.hxx"
#include "Parameters.hxx"

bool debug=false;

//********************************************************************
FluxWeightSystematics::FluxWeightSystematics(const std::string& name):BinnedParams(name,k2D_SYMMETRIC){
//********************************************************************
   SetNParameters(GetNBins());

   _flux = new FluxWeighting(ND::params().GetParameterS("psycheSystematics.FluxWeighting.Folder"),
                             ND::params().GetParameterS("psycheSystematics.FluxWeighting.Version"),
                             ND::params().GetParameterS("psycheSystematics.FluxWeighting.Tuning"));
   
}

//********************************************************************
Weight_h FluxWeightSystematics::ComputeWeight(const ToyExperiment& toy, const AnaEventC& eventBB, const ToyBoxB& boxB){
//********************************************************************

  const AnaEventB& event = *static_cast<const AnaEventB*>(&eventBB); 

  // Cast the ToyBox to the appropriate type
  const ToyBoxND280& box = *static_cast<const ToyBoxND280*>(&boxB); 

  Weight_h eventWeight=1;

  if (event.GetIsSandMC()) return eventWeight; // No weight for sand MC
  if (!this->IsCorrectRunPeriod(anaUtils::GetRunPeriod(event.EventInfo.Run))) return eventWeight; //should be in the correct run
  if (!box.TrueVertex) return eventWeight;         // True vertex associated to the recon vertex should exist

  // Get the true neutrino energy
  Float_t enu = box.TrueVertex->NuEnergy;
  int type = box.TrueVertex->NuPDG;

  // Get the flux error values for this energy
  Float_t mean_dummy;
  Float_t sigma;
  Int_t index;
  if(!GetBinValues(enu, type, mean_dummy, sigma, index)) return eventWeight;
 
  // Actual mean from FluxWeighting class
  Float_t mean=  _flux->GetWeight(box.TrueVertex,anaUtils::GetRunPeriod(event.EventInfo.Run));

  // Compute the weight
  eventWeight.Systematic = 1 + (mean-1) + sigma * toy.GetToyVariations(_index)->Variations[index];
  eventWeight.Correction = 1 + (mean-1);
  if(debug)
  {
    std::cout << "Run: " << event.EventInfo.Run << " Enu: " << enu << " type: " << type << " mean " << mean << " sigma " << sigma 
               << " index " << index << " variation " << toy.GetToyVariations(_index)->Variations[index] << " weight " << eventWeight << std::endl;
  }

  return eventWeight;

}

//********************************************************************
FluxWeightSystematicsNeutrino::FluxWeightSystematicsNeutrino():FluxWeightSystematics("FluxWeightNeutrino"){
//********************************************************************

}

//********************************************************************
bool FluxWeightSystematicsNeutrino::IsCorrectRunPeriod(int runPeriod){
//********************************************************************
   return (runPeriod<8);

}


//********************************************************************
FluxWeightSystematicsAntiNeutrino::FluxWeightSystematicsAntiNeutrino():FluxWeightSystematics("FluxWeightAntiNeutrino"){
//********************************************************************

}

//********************************************************************
bool FluxWeightSystematicsAntiNeutrino::IsCorrectRunPeriod(int runPeriod){
//********************************************************************

   return (runPeriod>=8);

}