DrawCCOscillations.C

Go to the documentation of this file.

/// DrawCCOscillations by J. Hewes <jhewes15@fnal.gov>

// CAFAna headers
#include "CAFAna/Core/Progress.h"
#include "CAFAna/Analysis/Calcs.h"
#include "NuXAna/macros/Nus20/Utilities.h"

// ROOT headers
#include "TCanvas.h"
#include "TGraph.h"
#include "TMarker.h"
#include "TLegend.h"

using namespace ana;

void DrawCCOscillations() {

  DontAddDirectory guard;
  gErrorIgnoreLevel = kError;
  std::vector<float> margin = { 0.15, 0.05 };

  auto calc = DefaultSterileCalc(4);
  SetNus20Params(calc);
  calc->SetL(1);

  // Flavours
  std::vector<int> initFlav = { 12, 14 };
  std::vector<int> finalFlav = { 12, 14, 16 };

  // Flavour names
  std::map<int, std::string> flavString = { {12, "#nu_{e}"}, {14, "#nu_{#mu}"},
    {16, "#nu_{#tau}"}, {-12, "#bar{#nu}_{e}"},
    {-14, "#bar{#nu}_{#mu}"}, {-16, "#bar{#nu}_{#tau}"} };


  Progress* p = new Progress("Drawing CC backgrounds");
  size_t counter = 0;

  Binning xBins = Binning::Simple(50, 0, 1);
  Binning yBins = Binning::LogUniform(50, 1e-2, 100);

  TH2D* hSurface = new TH2D(UniqueName().c_str(), ";sin^{2} #theta_{24};#Delta m_{41}^{2}",
    xBins.NBins(), &xBins.Edges()[0], yBins.NBins(), &yBins.Edges()[0]);

  std::vector<float> th34Vals = { 0., 20., 45. };
  std::vector<size_t> colours = { kRed-6, kGreen-5, kBlue-6 };

  // Loop over flavours
  size_t nPoints = 10000;
  for (int init : initFlav) {
    for (int final : finalFlav) {
      for (int sign = 0; sign <= 1; ++sign) {
        int i = sign? init : -1 * init;
        int f = sign? final : -1 * final;

        // Get the three-flavour oscillation probability
        SetNus20Params(calc);
        std::vector<float> x3f(nPoints);
        std::vector<float> y3f(nPoints);
        for (size_t it = 1; it <= nPoints; ++it) {
    float energy = (20./(float)nPoints) * (float)it;
    x3f[it] = energy;
    y3f[it] = calc->P(i, f, energy);
        }
        TGraph g3f(x3f.size(), &x3f[0], &y3f[0]);
        std::ostringstream title;
        title << "P(" << flavString[i] << " #rightarrow " << flavString[f] << ")";
        g3f.SetTitle(title.str().c_str());
        g3f.GetXaxis()->SetTitle("True E (GeV)");
        g3f.GetYaxis()->SetTitle("Osc. prob.");
        g3f.SetLineWidth(3);
        g3f.SetLineColor(kGray);

        // Loop over surface points
        for (int yVal = -2; yVal < 3; ++yVal) {
    double dm = pow(10, yVal);
    kFitDmSq41Sterile.SetValue(calc, dm);
    std::vector<std::string> images;
    for (size_t xBin = 1; xBin <= (size_t)xBins.NBins(); ++xBin) {
      double ssth24 = hSurface->GetXaxis()->GetBinCenter(xBin);
      kFitSinSqTheta24Sterile.SetValue(calc, ssth24);
      std::vector<float> x(nPoints);
      std::vector<std::vector<float>> y(th34Vals.size(), std::vector<float>(nPoints));
      std::vector<std::vector<float>> yDiff(th34Vals.size(), std::vector<float>(nPoints));
      for (size_t it = 1; it <= nPoints; ++it) {
        float energy = (20./(float)nPoints) * (float)it;
        x[it] = energy;
        for (size_t itTh34 = 0; itTh34 < th34Vals.size(); ++itTh34) {
          kFitTheta34InDegreesSterile.SetValue(calc, th34Vals[itTh34]);
          y[itTh34][it] = calc->P(i, f, energy);
          yDiff[itTh34][it] = y[itTh34][it] - y3f[it];
        }
      }
      TCanvas c("c", "c", 1600, 900);
      c.Divide(3);
      // Draw surface pad
      c.cd(1);
      gPad->SetPad(0.01, 0.01, 0.49, 0.99);
      gPad->SetLogy();
      hSurface->Draw("colz");
      TMarker m(ssth24, dm, 29);
      m.SetMarkerSize(3);
      m.Draw("same");
      // Draw oscillation plot
      c.cd(2);
      gPad->SetPad(0.51, 0.51, 0.99, 0.99);
      g3f.Draw("al");
      std::vector<TGraph*> g;
      TLegend l1(0.6, 0.15, 0.85, 0.4);
      for (size_t it = 0; it < y.size(); ++it) {
        g.push_back(new TGraph(x.size(), &x[0], &y[it][0]));
        g[it]->SetLineWidth(3);
        g[it]->SetLineColor(colours[it]);
        g[it]->Draw("l same");
        std::ostringstream name;
        name << "#theta_{34} = " << th34Vals[it] << "#circ";
        l1.AddEntry(g[it], name.str().c_str(), "l");
      }
      l1.Draw();
      // Draw ratio plot
      c.cd(3);
      gPad->SetPad(0.51, 0.01, 0.99, 0.49);
      std::vector<TGraph*> gDiff;
      TLegend l2(0.6, 0.15, 0.85, 0.4);
      for (size_t it = 0; it < yDiff.size(); ++it) {
        gDiff.push_back(new TGraph(x.size(), &x[0], &yDiff[it][0]));
        gDiff[it]->SetLineWidth(3);
        gDiff[it]->SetLineColor(colours[it]);
        gDiff[it]->SetTitle("");
        gDiff[it]->GetXaxis()->SetTitle("True E (GeV)");
        gDiff[it]->GetYaxis()->SetTitle("P(4flav) - P(3flav)");
        std::string drawOpts = it == 0 ? "al" : "l same";
        gDiff[it]->Draw(drawOpts.c_str());
        std::ostringstream name;
        name << "#theta_{34} = " << th34Vals[it] << "#circ";
        l2.AddEntry(gDiff[it], name.str().c_str(), "l");
      }
      l2.Draw();
      // Save plot
      std::ostringstream fileName;
      fileName << "plots/cc_oscprob/prob_flav" << i << "to" << f
        << "_dm" << dm << "_x" << xBin << ".png";
      c.cd(0);
      c.SaveAs(fileName.str().c_str());
      images.push_back(fileName.str());
      c.Clear();
      for (TGraph* tmp : g) delete tmp;
      for (TGraph* tmp : gDiff) delete tmp;
    } // for ssth24 value
    std::ostringstream gifName;
    gifName << "plots/cc_oscprob/prob_flav" << i << "to" << f
      << "_dm" << dm << ".gif";
    MakeGif(images, gifName.str(), 0.1);
    for (std::string pngFile : images) {
      std::ostringstream delCmd;
      delCmd << "rm " << pngFile;
      gSystem->Exec(delCmd.str().c_str());
    }
    ++counter;
    p->SetProgress(float(counter)/float(initFlav.size()*finalFlav.size()*10));
        } // for dm41 value
      } // for sign
    } // for final flavour
  } // for initial flavour

  p->Done();

} // Main macro function