MakeFakeDataLLTest.C

Go to the documentation of this file.

/// MakeFakeData.C by J. Hewes <jhewes15@fnal.gov>
/// Modified by   David Duenas <dueasdf@mail.uc.edu>
//#include "CAFAna/Prediction/CovarianceMatrix.h"
#include "CAFAna/Analysis/Exposures.h"
#include "CAFAna/Analysis/CovMxManager.h"
#include "NuXAna/macros/Nus20/Utilities.h"
#include "Utilities/rootlogon.C"
#include "CAFAna/Core/Progress.h"
#include "CAFAna/Analysis/Style.h"

#include "TH1.h"
#include "TCanvas.h"
#include "TMath.h"
#include "TLegend.h"
#include "TError.h"

using namespace ana;

void FakeSignal(osc::OscCalcSterile* calc, size_t n) {

  SetNus20Params(calc);

  if (n == 1) {
    kFitSinSqTheta24Sterile.SetValue(calc, 0.1);
    kFitSinSqTheta34Sterile.SetValue(calc, 0.5);
    kFitDmSq41Sterile.SetValue(calc, 0.1);
  } else if (n == 2) {
    kFitSinSqTheta24Sterile.SetValue(calc, 0.01);
    kFitSinSqTheta34Sterile.SetValue(calc, 0.5);
    kFitDmSq41Sterile.SetValue(calc, 1.0);
  } else if (n == 3) {
    kFitSinSqTheta24Sterile.SetValue(calc, 0.8);
    kFitSinSqTheta34Sterile.SetValue(calc, 0.5);
    kFitDmSq41Sterile.SetValue(calc, 10.0);
  } else {
    throw std::runtime_error("Fake signal number not valid!");
  }
} // function FakeSignal

void MiniBooNEParameters(osc::OscCalcSterile* calc)
{ 
  //MiniBooNe Best Fit values
  Double_t deltm2_41=0.041;
  Double_t theta_14=TMath::ASin(TMath::Sqrt(0.5000));
  Double_t theta_24=TMath::ASin(TMath::Sqrt(0.9580));
  Double_t theta_34=TMath::ASin(TMath::Sqrt(0.9580));
  calc->SetDm(4, deltm2_41);
  calc->SetAngle(1, 4, theta_14);
  calc->SetAngle(2, 4, theta_24);
  calc->SetAngle(3, 4, theta_34);
 
}

void MakeFakeData(bool rhc = false, bool stats = true)
{
  try { // i know wrapping the entire macro in a try/catch is bad but here we are

    DontAddDirectory guard;
    rootlogon(); 
    
    //Number of Universes
    size_t nUniverses = 1000;

    osc::OscCalcSterile* calc = DefaultSterileCalc(4);

    // TODO: confirm arguments to SetPrediction
    covmx::Sample nd = kNusFHCNearDet;
    nd.SetPOT(kAna2018SensitivityFHCNDPOT);
    nd.SetLivetime(-1);
    nd.SetBinning(kNCNDBinning);
    SetPrediction(nd, false); // Get ND prediction without systs

    covmx::Sample fd = kNusFHCFarDet;
    fd.SetPOT(kAna2018FHCPOT);
    fd.SetLivetime(kAna2018FHCLivetime);
    fd.SetBinning(kNCFDBinning);
    SetPrediction(fd, false); // Get FD prediction without systs

    TFile* f = TFile::Open("fhc_fake_data_lltest.root", "recreate");

    Spectrum* cosmic = LoadCosmic(kNusFHCFarDet, kUseCVMFS).release();
    TH1D* hCosmic = cosmic->ToTH1(fd.GetLivetime(), kLivetime);

    // Predict three-flavour spectra here
    SetNus20Params(calc);

    TH1D* hND3Flav = nd.GetPrediction()->Predict(calc).ToTH1(nd.GetPOT());
    hND3Flav->GetYaxis()->SetMaxDigits(4);
    hND3Flav->SetLineWidth(2);
    hND3Flav->SetLineColor(kAzure-2);
    hND3Flav->GetXaxis()->CenterTitle();
    hND3Flav->GetYaxis()->CenterTitle();

    TH1D* hFD3Flav = fd.GetPrediction()->Predict(calc).ToTH1(fd.GetPOT());
    hFD3Flav->GetYaxis()->SetMaxDigits(4);
    hFD3Flav->SetLineWidth(2);
    hFD3Flav->SetLineColor(kGreen+2);
    hFD3Flav->GetXaxis()->CenterTitle();
    hFD3Flav->GetYaxis()->CenterTitle();
    hFD3Flav->Add(hCosmic);

    // Get list of active systematics
    std::vector<std::string> systs;
    for (const ISyst* syst : getAllXsecSysts_2018_RPAFix()) {
      systs.push_back(syst->ShortName());
    }
    for (const ISyst* syst : LoadSysts(nd, kUseCVMFS)) {
      systs.push_back(dynamic_cast<const NuISyst*>(syst)->BaseName());
    }

    // Load matrix from manager
    std::string mxPath = kNus20Path + "/covmx/covmxmanager.root";
    TFile* mxFile = TFile::Open(pnfs2xrootd(mxPath).c_str(), "read");
    CovMxManager* mgr = CovMxManager::LoadFrom(mxFile).release();
    covmx::CovarianceMatrix* mx = mgr->GetCovarianceMatrix({nd,fd}, systs);
    mx->Predict({nd,fd}, calc);

    // Get binning information
    unsigned int nBins = mx->GetBinning().NBins();

    // Loop over every set of fake oscillation parameters
    for (size_t u = 1; u < 4; ++u) {

      TDirectory* uDir = f->mkdir(Form("fs%zu", u));
   
      FakeSignal(calc, u);
      mx->Predict({nd,fd}, calc);

      TH1D* hND4Flav = nd.GetPrediction()->Predict(calc).ToTH1(nd.GetPOT());
      hND4Flav->GetYaxis()->SetMaxDigits(4);
      hND4Flav->SetLineWidth(2);
      hND4Flav->SetLineColor(kRed-3);
      hND4Flav->GetXaxis()->CenterTitle();
      hND4Flav->GetYaxis()->CenterTitle();

      TH1D* hFD4Flav = fd.GetPrediction()->Predict(calc).ToTH1(fd.GetPOT());
      hFD4Flav->GetYaxis()->SetMaxDigits(4);
      hFD4Flav->SetLineWidth(2);
      hFD4Flav->SetLineColor(kRed-3);
      hFD4Flav->GetXaxis()->CenterTitle();
      hFD4Flav->GetYaxis()->CenterTitle();
      hFD4Flav->Add(hCosmic);

      // Combine nominal 4-flavour spectra into a single histogram

      Progress* p = new Progress("Generating fake universes");
      for (size_t i = 1; i <= nUniverses; ++i) {

        TCanvas* c3vs4 = new TCanvas(Form("c%zu", i), "", 1600, 1200);  
        c3vs4->Divide(3, 2);
          
        //ND Plots (NO Fluctuations)                                      
        c3vs4->cd(1);
        TLegend l1(0.4,0.5,0.8,0.79);
        l1.AddEntry(hND3Flav,"#splitline{ND Fake Data}{#splitline{3 Flavor }{NO Fluctuations}}","l");
        hND3Flav->Draw("hist");
        l1.Draw();

        //FD Plots (No Fluctuations)
        c3vs4->cd(4);
        TLegend l2(0.4,0.13,0.8,0.42);
        l2.AddEntry(hFD3Flav,"#splitline{FD Fake Data}{#splitline{3 Flavor }{NO Fluctuations}}","l");
        hFD3Flav->Draw("hist");
        l2.Draw(); 
         
        TDirectory* dir = uDir->mkdir(Form("universe_%zu", i));
        std::vector<Spectrum> specs = mx->GetUniverse(calc, {nd,fd});
               
        TH1D* hsND = specs[0].ToTH1(nd.GetPOT());
        TH1D* hsFD = specs[1].ToTH1(fd.GetPOT());
        hsFD->Add(hCosmic);

        //ND Fluctuation Plots
        c3vs4->cd(2);
        TLegend l3(0.4,0.5,0.8,0.79);
        l3.AddEntry(hsND,"#splitline{ND Fake Data}{#splitline{4 Flavor }{Fluctuations}}","l");
        hsND->GetYaxis()->SetMaxDigits(4);
        hsND->SetLineWidth(2);
        hsND->SetLineColor(kAzure-2);
        hsND->Draw("hist");
        hsND->GetXaxis()->CenterTitle();
        hsND->GetYaxis()->CenterTitle();
        if (hND4Flav->GetMaximum() > hsND->GetMaximum())
          hsND->SetMaximum(hND4Flav->GetMaximum());
        l3.AddEntry(hND4Flav, "#splitline{4 Flavor}{No fluctuations}", "l");
        hND4Flav->Draw("hist same");
        l3.Draw();

        //FD Fluctuation Plots     
        c3vs4->cd(5);
        TLegend l4(0.4,0.13,0.8,0.42);
        l4.AddEntry(hsFD,"#splitline{FD Fake Data}{#splitline{4 Flavor }{Fluctuations}}","l");
        hsFD->SetLineWidth(2);
        hsFD->SetLineColor(kGreen+2);
        hsFD->Draw("hist");
        hsFD->GetXaxis()->CenterTitle();
        hsFD->GetYaxis()->CenterTitle();
        if (hFD4Flav->GetMaximum() > hsFD->GetMaximum())
          hsFD->SetMaximum(hFD4Flav->GetMaximum());
        l4.AddEntry(hFD4Flav, "#splitline{4 Flavor}{No fluctuations}", "l");
        hFD4Flav->Draw("hist same");
        l4.Draw();
                
        //ratio Plots
        c3vs4->cd(3);
        TH1D* hratioND = (TH1D*)hND3Flav->Clone("hratioND");
        hratioND->Divide(hsND);
        hratioND->GetYaxis()->SetRangeUser(0,2);
        hratioND->GetYaxis()->SetTitle("ND 3 Flavor/4 Flavor");
        hratioND->SetLineColor(kBlack);
        hratioND->Draw("hist");
              
        c3vs4->cd(6);
        TH1D* hratioFD = (TH1D*)hFD3Flav->Clone("hratioFD");
        hratioFD->Divide(hsFD);
        hratioFD->GetYaxis()->SetRangeUser(0,2.5);
        hratioFD->GetYaxis()->SetTitle("FD 3 Flavor/4 Flavor");
        hratioFD->SetLineColor(kBlack);
        hratioFD->Draw("hist");
        
        // saving Fake data with the SetNewParameters() values
        specs[0].SaveTo(dir->mkdir(kNusFHCNearDet.GetTag().c_str()));
        specs[1].SaveTo(dir->mkdir(kNusFHCFarDet.GetTag().c_str()));

        //Directory to store FakeData plots
        gErrorIgnoreLevel = kError;
        gSystem->Exec(Form("mkdir -p plots/FakeDataPlots/fs%zu", u));
        c3vs4->SaveAs(Form("plots/FakeDataPlots/fs%zu/universe_%zu.png", u, i));
        gErrorIgnoreLevel = kWarning;

        p->SetProgress(double(i)/double(nUniverses));

        delete c3vs4;
        delete hsND;
        delete hsFD;
        delete hratioND;
        delete hratioFD;
      }
      p->Done();

    }

    delete f;

  } catch(const std::exception& e) {
    std::cout << std::endl << "Exception occurred! " << e.what() << std::endl;
  }

}