demo_CPT.C

Go to the documentation of this file.

#include "CAFAna/Vars/Vars.h"
#include "CAFAna/Vars/FitVarsCPT.h"
#include "CAFAna/Experiment/MultiExperiment.h"
#include "3FlavorAna/Cuts/NumuCuts.h"
#include "CAFAna/Prediction/PredictionNoExtrap.h"
#include "CAFAna/Experiment/SingleSampleExperiment.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Fit/FrequentistSurface.h"
#include "3FlavorAna/Vars/Binnings.h"
#include "OscLib/OscCalcCPT.h"
#include "TCanvas.h"
#include "TGaxis.h"
#include "TMath.h"

using namespace ana;

void demo_CPT()
{

  const std::string fname =
    "$NOVA_PROD/mc/FA14-10-28/genie/fd/caf/010000/1000000/*_fhc_nonswap_*";
  const std::string fnameSwap =
    "$NOVA_PROD/mc/FA14-10-28/genie/fd/caf/010000/1000000/*_fhc_swap_*";
  const std::string fnameRHC =
    "$NOVA_PROD/mc/FA14-10-28/genie/fd/caf/010000/1000000/*_rhc_nonswap_*";
  const std::string fnameRHCSwap =
    "$NOVA_PROD/mc/FA14-10-28/genie/fd/caf/010000/1000000/*_rhc_swap_*";

  SpectrumLoader loader(fname);
  SpectrumLoader loaderSwap(fnameSwap);
  SpectrumLoader loaderRHC(fnameRHC);
  SpectrumLoader loaderRHCSwap(fnameRHCSwap);

  PredictionNoExtrap predCC(   loader, loaderSwap,
                               "Simple Energy (GeV)", kNumuEnergyBinning,
                               kCCE, kNumuNonQEFD );
  PredictionNoExtrap predQE(   loader, loaderSwap,
                               "QE Energy (GeV)",     kNumuEnergyBinning,
                               kQEE, kNumuQEFD );
  PredictionNoExtrap predCCRHC(loaderRHC, loaderRHCSwap,
                               "Simple Energy (GeV)", kNumuEnergyBinning,
                               kCCE, kNumuNonQEFD );
  PredictionNoExtrap predQERHC(loaderRHC, loaderRHCSwap,
                               "QE Energy (GeV)", kNumuEnergyBinning,
                               kQEE, kNumuQEFD );


  loader.Go();
  loaderSwap.Go();
  loaderRHC.Go();
  loaderRHCSwap.Go();


  const double pot = 18e20;
  const double potrhc = 18e20;

  //Set up the CPT-violating oscillation calculators
  osc::OscCalcCPT calcR;
  calcR.SetL(810);
  calcR.SetRho(2.75);
  calcR.SetDmsq21(7.59e-5);
  calcR.SetDmsq32(2.41e-3);
  calcR.SetTh12(.601);           //sin2(2t12) = .870
  calcR.SetTh13(.1567);          //sin2(2t13) = .095
  calcR.SetTh23(TMath::Pi()*.2141);  //sin2(2t23) = .95

  osc::OscCalcCPT calcB;
  calcB.SetL(810);
  calcB.SetRho(2.75);
  calcB.SetDmsq21(7.59e-5);
  calcB.SetDmsq32(2.41e-3, osc::ENuSign::kNu);
  calcB.SetDmsq32(3.00e-3, osc::ENuSign::kNuBar);
  calcB.SetTh12(.601);           //sin2(2t12) = .870
  calcB.SetTh13(.1567);          //sin2(2t13) = .095
  calcB.SetTh23(TMath::Pi()*.2141, osc::ENuSign::kNu);    //sin2(2t23) = .95
  calcB.SetTh23(TMath::Pi()*.1867, osc::ENuSign::kNuBar); //sin2(2t23) = .85

  //Make fake data
  Spectrum obsCCR = predCC.Predict(&calcR);
  Spectrum obsQER = predQE.Predict(&calcR);
  Spectrum obsCCRHCR = predCCRHC.Predict(&calcR);
  Spectrum obsQERHCR = predQERHC.Predict(&calcR);

  Spectrum fakeCCR = obsCCR.FakeData(pot);
  Spectrum fakeQER = obsQER.FakeData(pot);
  Spectrum fakeCCRHCR = obsCCRHCR.FakeData(potrhc);
  Spectrum fakeQERHCR = obsQERHCR.FakeData(potrhc);


  Spectrum obsCCB = predCC.Predict(&calcB);
  Spectrum obsQEB = predQE.Predict(&calcB);
  Spectrum obsCCRHCB = predCCRHC.Predict(&calcB);
  Spectrum obsQERHCB = predQERHC.Predict(&calcB);

  Spectrum fakeCCB = obsCCB.FakeData(pot);
  Spectrum fakeQEB = obsQEB.FakeData(pot);
  Spectrum fakeCCRHCB = obsCCRHCB.FakeData(potrhc);
  Spectrum fakeQERHCB = obsQERHCB.FakeData(potrhc);

  //Make two experiments
  SingleSampleExperiment exptCCR(&predCC, fakeCCR);
  SingleSampleExperiment exptQER(&predQE, fakeQER);
  SingleSampleExperiment exptCCRHCR(&predCCRHC, fakeCCRHCR);
  SingleSampleExperiment exptQERHCR(&predQERHC, fakeQERHCR);

  std::vector<const IExperiment*> exptsR;

  exptsR.push_back(&exptCCR);
  exptsR.push_back(&exptQER);
  exptsR.push_back(&exptCCRHCR);
  exptsR.push_back(&exptQERHCR);

  MultiExperiment exptR(exptsR);


  SingleSampleExperiment exptCCB(&predCC, fakeCCB);
  SingleSampleExperiment exptQEB(&predQE, fakeQEB);
  SingleSampleExperiment exptCCRHCB(&predCCRHC, fakeCCRHCB);
  SingleSampleExperiment exptQERHCB(&predQERHC, fakeQERHCB);

  std::vector<const IExperiment*> exptsB;

  exptsB.push_back(&exptCCB);
  exptsB.push_back(&exptQEB);
  exptsB.push_back(&exptCCRHCB);
  exptsB.push_back(&exptQERHCB);

  MultiExperiment exptB(exptsB);

  //Draw the surfaces
  TGaxis::SetMaxDigits(3);

  const int res = 15; //10 is fast but ugly, 100 is slow and pretty
                      //pick something in between
  const double xmin = -.15;
  const double xmax = .25;
  const double ymin = -1.5e-3;
  const double ymax = 1e-3;

  TCanvas* c1 = new TCanvas("c1", "cptPlot", 1024, 768);


  SigmaDelta sdtB(&calcB, kFitSinSq2Theta23nu, kFitSinSq2Theta23bar);
  SigmaDelta sdmB(&calcB, kFitDmSq32nuNorm, kFitDmSq32barNorm);
  FitDeltaCPT kFitDeltaCPTSinSq2Theta23B(&sdtB);
  FitSigmaCPT kFitSigmaCPTSinSq2Theta23B(&sdtB);
  FitDeltaCPT kFitDeltaCPTDmSq32NormB(&sdmB);
  FitSigmaCPT kFitSigmaCPTDmSq32NormB(&sdmB);

  FrequentistSurface surfB2(&exptB, &calcB,
                 &kFitDeltaCPTSinSq2Theta23B, res, xmin, xmax,
                 &kFitDeltaCPTDmSq32NormB,    res, ymin, ymax,
                 {&kFitSigmaCPTSinSq2Theta23B,&kFitSigmaCPTDmSq32NormB}
                );

  surfB2.DrawBestFit(kBlue);
  surfB2.DrawContour(Gaussian68Percent2D(surfB2), 1, kBlue);
  surfB2.DrawContour(Gaussian2Sigma2D(surfB2),    2, kBlue);
  surfB2.DrawContour(Gaussian3Sigma2D(surfB2),    3, kBlue);


  TLine *xaxis = new TLine(xmin,0,xmax,0);
  xaxis->SetLineColor(kGray);
  xaxis->SetLineStyle(7);
  xaxis->Draw();

  TLine *yaxis = new TLine(0,ymin,0,ymax);
  yaxis->SetLineColor(kGray);
  yaxis->SetLineStyle(7);
  yaxis->Draw();


  SigmaDelta sdtR(&calcR, kFitSinSq2Theta23nu, kFitSinSq2Theta23bar);
  SigmaDelta sdmR(&calcR, kFitDmSq32nuNorm, kFitDmSq32barNorm);
  FitDeltaCPT kFitDeltaCPTSinSq2Theta23R(&sdtR);
  FitSigmaCPT kFitSigmaCPTSinSq2Theta23R(&sdtR);
  FitDeltaCPT kFitDeltaCPTDmSq32NormR(&sdmR);
  FitSigmaCPT kFitSigmaCPTDmSq32NormR(&sdmR);

  FrequentistSurface surfR2(&exptR, &calcR,
                 &kFitDeltaCPTSinSq2Theta23R, res, xmin, xmax,
                 &kFitDeltaCPTDmSq32NormR,    res, ymin, ymax,
                 {&kFitSigmaCPTSinSq2Theta23R,&kFitSigmaCPTDmSq32NormR}
                );

  surfR2.DrawBestFit(kRed);
  surfR2.DrawContour(Gaussian68Percent2D(surfR2), 1, kRed);
  surfR2.DrawContour(Gaussian2Sigma2D(surfR2),    2, kRed);
  surfR2.DrawContour(Gaussian3Sigma2D(surfR2),    3, kRed);


  c1->SaveAs("CPT-sesitivity.pdf");

}