do_extrap.h

Go to the documentation of this file.

#ifndef __CINT__

#include "MCReweight/MCReweight.h"
#include "MCReweight/NeugenWeightCalculator.h"
#include "NCUtils/Extrapolation/NCExtrapolationModule.h"

#include "NCUtils/NCOscProb.h"
#include "NCUtils/NCType.h"

#include "TString.h"

#include <iostream>

enum EFitType {
  kRealData,
  kFakeData,
  kFakeDataSplitMC,
  kMockData
};

// The Module and its Registry as globals so they don't need to be
// passed to all the functions
NCExtrapolationModule extrapMod;
Registry reg;

// =====================================================================
//         Utility functions
// =====================================================================

TString add_to_output_name(TString s, bool appendSep=true)
{
  static TString outputName;
  outputName.Append(s);
  if (appendSep) outputName.Append("_");

  return outputName;
}

// =====================================================================
//        Basics
// =====================================================================

void set_msg_levels()
{
  MsgService* m=MsgService::Instance();
  m->GetStream("NCExtrapolationModule")->SetLogLevel(Msg::kInfo);
  m->GetStream("NCAnalysisCuts")->SetLogLevel(Msg::kInfo);
  m->GetStream("NCUtils")->SetLogLevel(Msg::kInfo);
  //m->GetStream("NCExtrapolation"->SetLogLevel(Msg::kInfo);
  //m->GetStream("MinosMinimization"->SetLogLevel(Msg::kDebug);
  m->GetStream("Registry")->SetLogLevel(Msg::kFatal);
}

void set_beam_types(TString beamTypes)
{
  reg.Set("BeamType", beamTypes);
}

void set_weight_config()
{
  //default is PiMinus_CedarDaikon
  reg.Set("MEGAWeightConfig", "PiMinus_CedarDaikon");
}


void set_near_run_limits(int dataLimit, int mcLimit)
{
  reg.Set("RunLimitNearData", dataLimit);
  reg.Set("RunLimitNearMC",   mcLimit);
}

void set_far_run_limits(int dataLimit, int mcLimit)
{
  reg.Set("RunLimitFarData", dataLimit);
  reg.Set("RunLimitFarMC",   mcLimit);
}

void set_run_to_use(int run)
{
  //set which runs to use - see NCRunUtil.h for the numbering
  reg.Set("RunToUse", run);
}

// =====================================================================
//        General settings
// =====================================================================
void set_data_mc_type(EFitType fitType,
                      int mockDataRun=-1,
                      /* -1 mock data subrun means "use all subruns",
                         so use -2 as default */
                      int mockDataSubRun=-2)
{

  if (fitType!=kMockData &&
      (mockDataRun!=-1 || mockDataSubRun!=-2))
    assert(0 && "Mock data parameters set when fit type is not mock data");

  if (fitType==kMockData &&
      (mockDataRun==-1 || mockDataSubRun==-2))
    assert(0 && "Mock data parameters not set when fit type is mock data");

  switch (fitType) {
  case kRealData:
    reg.Set("UseMCAsData", false);
    reg.Set("UseMockData", false);
    add_to_output_name("data");
    break;
  case kFakeData:
    reg.Set("UseMCAsData", true);
    reg.Set("UseMockData", false);
    reg.Set("SplitMC",     false);
    add_to_output_name("mconly");
    break;
  case kFakeDataSplitMC:
    reg.Set("UseMCAsData", true);
    reg.Set("UseMockData", false);
    reg.Set("SplitMC",     true);
    add_to_output_name("splitmc");
    break;
  case kMockData:
    reg.Set("UseMCAsData", false);
    reg.Set("UseMockData", true);
    reg.Set("MockDataSet", Form("F%d", mockDataRun));
    reg.Set("MockDataSubRun", mockDataSubRun);
    add_to_output_name("mock");
    add_to_output_name(Form("run%d", mockDataRun));
    add_to_output_name(Form("%d", mockDataSubRun));
    break;
  default:
    assert(0 && "Invalid fit type");
  }
}

void set_data_mc_path(TString uDSTDir)
{
  reg.Set("DataMCPath", uDSTDir);
}

void set_mc_exposure_for_beam(BeamType::BeamType_t beamType, NC::RunUtil::ERunType run, double pot)
{
  extrapMod.SetMCExposureForBeam(NCBeam::Info(beamType, run), pot);
}

void set_extraction(TString name)
{
  // Add a dot at the end if we missed it
  if (!name.EndsWith(".")) name+=".";

  // pick the extraction to use - use only 1 at a time
  // in this module.
  //Loop through extractions and use the one we want
  for (Int_t i = 0; i < NCType::kNumExtractions; ++i){
    if(name == NCType::kExtractionNames[i]){
      reg.Set("ExtractionType", i);
      return;
    }
  }

  assert(0 && "Specified extraction not found");
}

void set_use_custom_pid_cut_and_value(bool useCustom,
                                      float value=-9999)
{
  reg.Set("UsePIDCustomCut", useCustom);
  if (useCustom)
    reg.Set("PIDCustomCut", value);
}

void set_prediction(bool makePred, bool writeFDData=false)
{
  reg.Set("Prediction", makePred);
  if(writeFDData) reg.Set("PredictionWriteData", writeFDData);
}


// =====================================================================
//          Systematics-related stuff
// =====================================================================

void set_syst_study_params(bool doSystStudy, float shiftInSigma)
{
  //Whether or not to do syst study.
  //If true, Adjust conditions are ignored in favor of SystShift(sigma).
  reg.Set("DoSystStudy", doSystStudy);
  reg.Set("SystShift", shiftInSigma);
}

void turn_off_systematic_params()
{
  // This is mostly copy-pasted from NCExtrapolationModule, so we don't
  // actually *have* to do it. Anyway, it doesn't hurt to make it
  // explicit in the macro

  TString adjust   = "Adjust";
  TString mcAsData = "ChangeMCAsData";
  TString fit      = "Fit";
  for(int i = 0; i < NCType::kNumSystematicParameters; ++i){
    reg.Set(NCType::kParams[i].name,         false);
    reg.Set(mcAsData + NCType::kParams[i].name, false);
    reg.Set(fit + NCType::kParams[i].name,      false);
  }
}

// No adjustment value given, so will use the value of "SystShift"
void adjust_systematic(NCType::EFitParam syst)
{
  reg.Set("ChangeMCAsData"+NCType::kParams[syst].name, true);
}

void adjust_systematic(NCType::EFitParam syst, float adjVal)
{
  reg.Set("ChangeMCAsData"+NCType::kParams[syst].name, true);
  reg.Set("Adjust"+NCType::kParams[syst].name, adjVal);
}

void fit_systematic(NCType::EFitParam syst)
{
  reg.Set("Fit"+NCType::kParams[syst].name, true);
}

// =====================================================================
//        Oscillation and fitter things
// =====================================================================

TString get_model_short_name(NCType::EOscModel m)
{
  switch(m) {
  case NCType::kThreeFlavor: return "3flavor";
  case NCType::kNoOscillations: return "noOsc";
  case NCType::kFourFlavorGeneral: return "4flavorGen";
  case NCType::kFourFlavorDelta43IsBig: return "delta43IsBig";
  case NCType::kFourFlavorDelta41Is0: return "delta41Is0";
  case NCType::kFourFlavorDelta43Is0: return "delta43Is0";
  case NCType::kSterileFraction: return "sFrac";
  case NCType::kDecay: return "decay";
  case NCType::kDecoherence: return "decoherence";
  default:
    assert(0 && "Unknown oscmodel");
    return "Unknown";
  }
}

TString get_param_string(NCType::EFitParam p)
{
  switch(p) {
  case NCType::kUE3Sqr:  return "UE3Sqr";
  case NCType::kTheta13: return "Theta13";
  case NCType::kDelta1:  return "Delta1";
  case NCType::kDelta2:  return "Delta2";
  case NCType::kTheta14: return "Theta14";
  case NCType::kTheta24: return "Theta24";
  case NCType::kTheta23: return "Theta23";
  case NCType::kFs:      return "Fs";
  default:
    assert(0 && "No String for parameter");
    return "Unknown";
  }

}

void fit_osc_param(NCType::EFitParam p)
{
  const TString paramString=get_param_string(p);
  reg.Set("Fit"+paramString, true);
}

void fix_osc_param(NCType::EFitParam p, double value)
{
  const TString paramString=get_param_string(p);
  reg.Set("Fit"+paramString, false);
  reg.Set(paramString+"Val", value);
}

void set_osc_pars(NC::OscProb::OscPars* oscPars)
{
  assert(sizeof(oscPars) == sizeof(int));
  reg.Set("OscPars", int(oscPars));
  add_to_output_name(get_model_short_name(oscPars->OscillationModel()));
}

void set_model_to_fit(NCType::EOscModel m)
{
  reg.Set("OscillationModel", int(m));
  add_to_output_name(get_model_short_name(m)+"Fit");
}

void use_fitter_confidence(int confidenceLevel, bool wanted=true)
{
  assert(confidenceLevel==68 || confidenceLevel==90 || confidenceLevel==99);
  reg.Set(Form("LazyContoursWant%d", confidenceLevel), wanted);
}

void set_fitter_params(float precScale=1.)
{
  use_fitter_confidence(68);
  use_fitter_confidence(90);
  use_fitter_confidence(99, false);

  reg.Set("PrecScale", precScale);
}

void add_contour(NCType::EFitParam x, NCType::EFitParam y)
{
  static TString contoursList("");
  contoursList+=Form("(%d,%d) ", x, y);

  reg.Set("ContoursList", contoursList);
}

void add_projection(NCType::EFitParam proj)
{
  static TString projsList;
  projsList += TString::Format("%d ", proj);
  reg.Set("ProjectionsList", projsList);
}

void do_syst_interpolation(bool interp)
{
  reg.Set("DoSystematicInterpolation", interp);
}

void make_shifted_beams(bool make)
{
  reg.Set("MakeShiftedBeams", make);
  reg.Set("ShiftedBeamsSimpleOnly", make);
}

void force_best_fit_to_truth(bool force)
{
  reg.Set("ForceBestFitToTruth", force);
}

// =====================================================================
//       Extrapolation settings
// =====================================================================

void farnear_smooth_width(double smoothingWidth)
{
  reg.Set("FarNearSmoothingWidth", smoothingWidth);
}

void pid_emulate_farnear(bool emulate)
{
  reg.Set("PIDFitEmulateFarNear", emulate);
}

void pid_true_bin_factor(int factor)
{
  reg.Set("PIDFitTrueBinFactor", factor);
}

void pid_binning(int nbins, double minpid, double maxpid)
{
  // Number of bins in PID to use
  reg.Set("PIDFitPIDNbins",          nbins);
  // Minimum, maximum value of PID
  reg.Set("PIDFitPIDmin",          minpid);
  reg.Set("PIDFitPIDmax",          maxpid);
}

// =====================================================================
//       Miscellaneous
// =====================================================================
void add_neugen_weight_calculator()
{
  MCReweight &mcReweight = MCReweight::Instance();
  NeugenWeightCalculator *fNeugenWeightCal = new NeugenWeightCalculator;
  mcReweight.AddWeightCalculator(fNeugenWeightCal);
}

void run()
{

  // Add the .root to the end of the output filename and set it in the module
  reg.Set("FileName", add_to_output_name(".root", false));

  add_neugen_weight_calculator();

  Registry def = extrapMod.DefaultConfig();
  def.UnLockValues();

  def.Merge(reg);
  extrapMod.Config(def);

  extrapMod.Run();
}

#endif

---