FCContour.C

Go to the documentation of this file.

#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Analysis/CountingExperimentStatistics.h"
#include "CAFAna/Analysis/Exposures.h"
#include "CAFAna/Fit/MinuitFitter.h"
#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/SystShifts.h"
#include "CAFAna/Experiment/CountingExperiment.h"
#include "CAFAna/Experiment/GaussianConstraint.h"
#include "CAFAna/Experiment/MultiExperiment.h"
#include "CAFAna/Prediction/IPrediction.h"
#include "CAFAna/Prediction/PredictionCombinePeriods.h"
#include "CAFAna/Prediction/PredictionInterp.h"
#include "NuXAna/Systs/NusSysts.h"
#include "CAFAna/Vars/FitVarsSterile.h"

#include "OscLib/OscCalcSterile.h"

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

using namespace ana;

double CalcSignificance(const std::vector<double>& num,
                        const std::vector<double>& den,
                        int nData,
                        double nPred);

void ResetAngles(osc::OscCalcSterile* calc); //should just make this globally somewhere

void FCContour(std::string outfile)
{
  // Get the data spectrum
  // Load in the NuMI data spectrum
  Spectrum* data = LoadFromFile<Spectrum>("$NOVA_ANA/steriles/Ana01/fout_nus_ana01_box_opening_restricted.root", "spec_nus_cale_numi").release();
  int nData = data->Integral(data->POT());
  std::cout << "Data spectrum: " << data->POT()
            << " and " << nData << " events" << std::endl;
  
  // Load in the out-of-NuMI-time weighted spectrum (cosmics)
  Spectrum* cosmic_scale = LoadFromFile<Spectrum>("$NOVA_ANA/steriles/Ana01/AnaResults.root", "cosmic__CalE").release();
  TH1D* hCosmic = cosmic_scale->ToTH1(cosmic_scale->Livetime(), kLivetime);
  Spectrum sCos(hCosmic, data->POT(), data->Livetime());
  double nCos = sCos.ToTH1(data->POT())->Integral(1, hCosmic->GetNbinsX());
  std::cout << "Total cosmic count: " << nCos << std::endl;
  
  // Load in the MC Prediction
  //  TFile* fileMC = new TFile("$NOVA_ANA/steriles/Ana01/nus_ana01_prediction.root", "READ");
  TFile* fileMC = new TFile("/nova/app/users/gsdavies/nus/ana.a/NuXAna/macros/nus_pred_osc_syst_incl_ana01.root");
  
  std::unique_ptr<PredictionInterp> pred123b =
    PredictionInterp::LoadFrom((TDirectory*)fileMC->Get("nus_pred_p1p2e3b"));
  
  std::unique_ptr<PredictionInterp> pred3c3d =
    PredictionInterp::LoadFrom((TDirectory*)fileMC->Get("nus_pred_e3c"));

  fileMC->Close();

  const double pot123b =
    kSecondAnaPeriod1POT
    + kSecondAnaPeriod2POT
    + kSecondAnaEpoch3bPOT;
  
  const double pot3c3d = kSecondAnaEpoch3cPOT + kSecondAnaEpoch3dPOT;
  
  // Combine period1/period2/epoch3b + epoch3c, weighted appropriately
  PredictionCombinePeriods pred({
    {&(*pred123b), pot123b},
    {&(*pred3c3d), pot3c3d}
  });
  
  // Initialize FC histogram
  std::string title = ";#theta_{34} (Deg.); #theta_{24} (Deg.)";
  int nbins34 = 200;
  double min34 = 0.;
  double max34 = 50.;
  int nbins24 = 180;
  double min24 = 0.;
  double max24 = 45.;
  TH2D* hFCSig = new TH2D("hFCSig", title.c_str(),
                          nbins34, min34, max34,
                          nbins24, min24, max24);
  TH2D* hFCChi = new TH2D("hFCChi", title.c_str(),
                          nbins34, min34, max34,
                          nbins24, min24, max24);

  osc::OscCalcSterile* calc4f = DefaultSterileCalc(4);
  calc4f->SetNFlavors(4);
  calc4f->SetDm(4, 0.5); // #Delta m^{2}_{41} = 0.5 eV^2
  calc4f->SetAngle(2, 4, 0.); // #theta_{24} = 0 degrees
  calc4f->SetAngle(3, 4, 0.); // #theta_{34} = 0 degrees
  
  std::vector<const IFitVar*> fitvars_num;
  fitvars_num.push_back(&kFitTheta23InDegreesSterile);
  fitvars_num.push_back(&kFitDelta24InPiUnitsSterile);

  std::vector<const IFitVar*> fitvars_den;
  fitvars_den.push_back(&kFitTheta23InDegreesSterile);
  fitvars_den.push_back(&kFitTheta24InDegreesSterile);
  fitvars_den.push_back(&kFitTheta34InDegreesSterile);
  fitvars_den.push_back(&kFitDelta24InPiUnitsSterile);

  std::vector<double> rank_den;

//  std::cout << "Rank denominator: " << std::endl;
  for(int i = 0; i <= nData; ++i) {

    ResetAngles(calc4f);

    // Create a fake data spectrum that only contains i events
    TH1* h = data->ToTH1(kSecondAnaPOT);
    h->Reset();
    h->SetBinContent(1, i);
    Spectrum fake(h, data->POT(), data->Livetime());

    GaussianConstraint th23Constraint(&kFitTheta23InDegreesSterile,
                                      45.8, 3.2);
    CountingExperiment nus_expt(&pred, fake, sCos);
    MultiExperiment    expt({&th23Constraint, &nus_expt});

    std::vector<const ISyst*> systs = getAllNusSysts();
    SystShifts shifts;
    for(const auto& syst : systs) {
      shifts.SetShift(syst, 0);
    }

    MinuitFitter fit(&expt, fitvars_den, systs);
    double chi2 = fit.Fit(calc4f, shifts, IFitter::kQuiet)->EvalMetricVal();

    double nPred = pred.PredictSyst(calc4f, shifts).Integral(kSecondAnaPOT);
    nPred += sCos.Integral(kSecondAnaPOT);

    rank_den.push_back(TMath::Poisson(i, nPred));
    std::cout << "Data value: " << i << std::endl;
//    std::cout << "Data value: " << i << ", " << nPred << ", "
//              << chi2 << ", " << rank_den[i] << std::endl;
  }
  std::cout << "Finished calculating R denominator values." << std::endl;

  for(int i34 = 1, n34 = hFCSig->GetNbinsX(); i34 <= n34; ++i34) {
    std::cout << "On loop over theta_34 = " << i34 << std::endl;
    
    for(int i24 = 1, n24 = hFCSig->GetNbinsY(); i24 <= n24; ++i24) {
      std::cout << "On loop over theta_24 = " << i24 << std::endl;

      ResetAngles(calc4f);
      calc4f->SetAngle(3, 4, hFCSig->GetXaxis()->GetBinCenter(i34)*M_PI/180.);
      calc4f->SetAngle(2, 4, hFCSig->GetYaxis()->GetBinCenter(i24)*M_PI/180.);
      std::cout << "t34: " << hFCSig->GetXaxis()->GetBinCenter(i34)*M_PI/180. << std::endl;
      std::cout << "t24: " << hFCSig->GetXaxis()->GetBinCenter(i24)*M_PI/180. << std::endl;

      std::vector<double> rank_num;
      
      GaussianConstraint th23Constraint(&kFitTheta23InDegreesSterile,
                                        45.8, 3.2);
      CountingExperiment nus_expt(&pred, *data, sCos);
      MultiExperiment    expt({&th23Constraint, &nus_expt});
      
      std::vector<const ISyst*> systs = getAllNusSysts();
      SystShifts shifts;
      for(const auto& syst : systs) {
        shifts.SetShift(syst, 0);
      }
      
      MinuitFitter fit(&expt, fitvars_num, systs);
      double chi2 = fit.Fit(calc4f, shifts, IFitter::kQuiet)->EvalMetricVal();
      std::cout << " chi2: " << chi2 << std::endl;      
      double nPred = pred.PredictSyst(calc4f, shifts).Integral(kSecondAnaPOT);
      std::cout << "npred: " << nPred << std::endl;
      nPred += sCos.Integral(kSecondAnaPOT);
      std::cout << "new nPred: " << nPred << std::endl;
      for(int i = 0; i <= nData; ++i) {
        chi2 += LogLikelihood(nPred, i);
        rank_num.push_back(TMath::Poisson(i, nPred));
//        std::cout << "Data value: " << i << ", " << nPred << ", "
//                  << chi2 << ", " << rank_num[i] << ", "
//                  << rank_num[i]/rank_den[i] << std::endl;
        chi2 -= LogLikelihood(nPred, i);
      }
      std::cout << "new chi2: " << chi2 << std::endl;
      double bin_significance = CalcSignificance(rank_num, rank_den, nData, nPred);
      std::cout << "bin_sig: " << bin_significance << std::endl;
      // Converting to “Gaussian” sigma
      const double sigma = PValueToSigma(1.0-bin_significance);

//      std::cout << "Bin significance: " << bin_significance << ", "
//                << "Sigma: " << sigma << ", "
//                << "Chi2: " << sigma*sigma << std::endl;
      hFCSig->SetBinContent(i34,i24, bin_significance);
      hFCChi->SetBinContent(i34,i24, sigma*sigma);
    } // end of loop over theta 24 bins
  } // end of loop over theta 34 bins

  TFile* rootF = new TFile(outfile.c_str(), "RECREATE");

  rootF->WriteTObject(hFCSig);
  rootF->WriteTObject(hFCChi);

  TCanvas* cSig = new TCanvas("cFCSig", "cFCSig", 800, 500);
  //CenterTitles(hFCSig);
  hFCSig->Draw("colz");
  gPad->RedrawAxis();
  //Preliminary();

  gPad->Print("FC2DSignificance_profd24.png");
  rootF->WriteTObject(cSig);

  TCanvas* cChi = new TCanvas("cFCChi", "cFCChi", 800, 500);
  //CenterTitles(hFCChi);
  hFCChi->Draw("colz");
  gPad->RedrawAxis();
  //Preliminary();

  gPad->Print("FC2DChiSq_profd24.png");
  rootF->WriteTObject(cChi);

  rootF->Close();

  std::cout << "Finished. Ctrl+Z should be safe if the macro didn't quit." << std::endl;
}

//------------------------------------------------------------------------------
double CalcSignificance(const std::vector<double>& num,
                        const std::vector<double>& den,
                        int nData,
                        double nPred)
{
  double prob_tot = 0.;

  std::map<double, std::vector<int> > rank;

  assert(num.size() == den.size());
  for(int i = 0, n = num.size(); i < n; ++i) {
    double r = num[i]/den[i];

    if(rank.find(r) == rank.end()) {
      rank[r] = std::vector<int>();
    }

    rank[r].push_back(i);
  }

  double rData = num[nData]/den[nData];
  //std::cout << "rData: " << rData << std::endl;

  std::map<double, std::vector<int> >::reverse_iterator top_rank = rank.rbegin();
  double rCurr = top_rank->first;

  while(rCurr >= rData && top_rank != rank.rend()) {
//    std::cout << "rCurr: " << rCurr << std::endl;
    for(const auto& nObs : top_rank->second) {
      prob_tot += TMath::Poisson(nObs, nPred);
    }

    rank.erase(rCurr);
    if(rank.size() == 0) { break; }
    top_rank = rank.rbegin();
    rCurr = top_rank->first;
  }

//  std::cout << "Rank Size after: " << rank.size() << std::endl;

  return prob_tot;
}

//------------------------------------------------------------------------------
void ResetAngles(osc::OscCalcSterile* calc)
{
  ResetSterileCalcToDefault(calc);
  calc->SetDm(4, 0.5);      // #deltam41 = 0.5 eV^2
  calc->SetAngle(2, 3, M_PI/4);
  calc->SetAngle(3, 4, 0.);
  calc->SetAngle(2, 4, 0.);
}