test_nue2017Prediction.C

Go to the documentation of this file.

#include "OscLib/OscCalcDumb.h"
#include "CAFAna/Analysis/Calcs.h"

#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/LoadFromFile.h"

#include "CAFAna/Vars/FitVars.h"

#include "3FlavorAna/Prediction/PredictionExtendToPeripheral.h"
#include "3FlavorAna/Prediction/PredictionAddRock.h"
#include "CAFAna/Prediction/PredictionNoExtrap.h"

#include "Utilities/func/MathUtil.h"

#include "TCanvas.h"
#include "TH1.h"
#include "TFile.h"

#include <iostream>
#include <cmath>

using namespace ana;

void test_nue2017Prediction()
{
  double pot = 9.48e20; // From Nitish

  std::string dirName = "/nova/ana/users/dpershey/Extraps2017/";
  std::string fName1 = "out_nue2017_extrap.root";
  std::string fName2 = "out_fdrock.root";

  // First step, take the extrapolated core sample and extend the dcmp result
  // to the peripheral sample
  // The extrapolated result
  PredictionExtrap *core_merge =
    LoadFromFile<PredictionExtrap>(dirName+fName1,"predProp").release();
  // The FD MC
  PredictionNoExtrap *fdmc_merge =
    LoadFromFile<PredictionNoExtrap>(dirName+fName1,"fdmcNom").release();
  // Those together give the fiducial prediction
  PredictionExtendToPeripheral pred_merge(core_merge,fdmc_merge,true);


  PredictionExtrap *core =
    LoadFromFile<PredictionExtrap>(dirName+fName1,"predProp").release();
  // The FD MC
  PredictionNoExtrap *fdmc =
    LoadFromFile<PredictionNoExtrap>(dirName+fName1,"fdmcNom").release();
  PredictionExtendToPeripheral pred(core,fdmc,false);

  // Now, make a plot
  // Close to one of the SA best fit points
  double dCP = 1.5;
  double ssth23 = 0.402;
  double dmsq32 = 0.0027;

  osc::IOscCalcAdjustable* calc = DefaultOscCalc();

  kFitDeltaInPiUnits.SetValue(calc,dCP);
  kFitDmSq32.SetValue(calc,dmsq32);
  kFitSinSqTheta23.SetValue(calc,ssth23);

  Spectrum tot = pred.Predict(calc);
  Spectrum sig = pred.PredictComponent(calc, Flavors::kNuMuToNuE,
                                       Current::kCC, Sign::kNu);
  Spectrum bb  = pred.Predict(calc);
  bb -= sig;

  TH1 *htot = tot.ToTH1(pot);
  TH1 *hsig = sig.ToTH1(pot);
  TH1 *hbb  = bb. ToTH1(pot);

  htot->SetLineColor(kRed);
  hsig->SetLineColor(kPink+9);
  hbb ->SetLineColor(kAzure);



  TCanvas *can = new TCanvas();
  htot->Draw("hist");
  hsig->Draw("hist,same");
  hbb ->Draw("hist,same");


  Spectrum tot_merge = pred_merge.Predict(calc);
  Spectrum sig_merge = pred_merge.PredictComponent(calc, Flavors::kNuMuToNuE,
                                                   Current::kCC, Sign::kNu);
  Spectrum bb_merge  = pred_merge.Predict(calc);
  bb_merge -= sig_merge;

  TH1 *htot_merge = tot_merge.ToTH1(pot);
  TH1 *hsig_merge = sig_merge.ToTH1(pot);
  TH1 *hbb_merge  = bb_merge. ToTH1(pot);

  htot_merge->SetLineColor(kRed);
  hsig_merge->SetLineColor(kPink+9);
  hbb_merge ->SetLineColor(kAzure);

  TCanvas *can_merge = new TCanvas();
  htot_merge->Draw("hist");
  hsig_merge->Draw("hist,same");
  hbb_merge ->Draw("hist,same");

  std::cout << "Compare the event counts in the peripheral: " << std::endl;
  std::cout << "Total:      Merged - " << htot_merge->Integral(28,28) << " and Not Merged - " << htot->Integral(28,36) << std::endl;
  std::cout << "Signal:     Merged - " << hsig_merge->Integral(28,28) << " and Not Merged - " << hsig->Integral(28,36) << std::endl;
  std::cout << "Beam Bkgd:  Merged - " << hbb_merge->Integral(28,28) << " and Not Merged - " << hbb->Integral(28,36) << std::endl;


}