plot_kinematics_cafana.C

Go to the documentation of this file.

#ifdef __CINT__
void plot_kinematics_cafana(){}
#else

#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Vars/Vars.h"
using namespace ana;

#include "StandardRecord/StandardRecord.h"

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

void Compare(const Spectrum* gen, const Spectrum* gib, const std::string& elem)
{
  TH1* hgen = gen->ToTH1(1e20, kBlue);
  TH1* hgib = gib->ToTH1(1e20, kRed);

  if(hgen->GetBinContent(hgen->GetMaximumBin()) > hgib->GetBinContent(hgib->GetMaximumBin())){
    hgen->SetTitle(("#nu_{#mu} CC on "+elem).c_str());
    hgen->GetYaxis()->SetTitle("Events / 10^{20} POT");
    hgen->Draw("hist");
    hgib->Draw("hist same");
  }
  else{
    hgib->SetTitle(("#nu_{#mu} CC on "+elem).c_str());
    hgib->GetYaxis()->SetTitle("Events / 10^{20} POT");
    hgib->Draw("hist");
    hgen->Draw("hist same");
  }
}

void plot_kinematics_cafana()
{
  SpectrumLoader lgenie("/nfs/raid5/bckhouse/gibuu/reprocess/*.orig.caf.root");
  SpectrumLoader lgibuu("/nfs/raid5/bckhouse/gibuu/reprocess/*.gibuu.caf.root");

  const Cut kPreselAny({"mc.nnu", "mc.nu.iscc", "mc.nu.pdg",
                        "mc.nu.vtx.x", "mc.nu.vtx.y", "mc.nu.vtx.z"},
    [](const caf::StandardRecord* sr)
    {
      if(sr->mc.nnu == 0) return false;
      if(!sr->mc.nu[0].iscc || !sr->mc.nu[0].pdg == 14) return false;
      if(sr->mc.nu[0].vtx.x < -200 || sr->mc.nu[0].vtx.x > +200 ||
         sr->mc.nu[0].vtx.y < -200 || sr->mc.nu[0].vtx.y > +200 ||
         sr->mc.nu[0].vtx.z < -100 || sr->mc.nu[0].vtx.z > 1600) return false;
      return true;
    });

  const Var kGiBUUWgt({"mc.nnu", "mc.nu.mode", "mc.nu.genweight"},
                      [](const caf::StandardRecord* sr)
                      {
                        if(sr->mc.nnu == 0) return 0.f;
                        if(sr->mc.nu[0].mode >= 0) return 1.f; // GENIE
                        return sr->mc.nu[0].genweight;
                      });

  const Var kLogWgt({"mc.nnu", "mc.nu.mode", "mc.nu.genweight"},
                    [](const caf::StandardRecord* sr)
                    {
                      if(sr->mc.nnu == 0) return -100.f;
                      if(sr->mc.nu[0].mode >= 0) return -100.f; // GENIE
                      return log10(sr->mc.nu[0].genweight);
                    });

  const Var kTgtZ({"mc.nnu", "mc.nu.tgtZ"},
                  [](const caf::StandardRecord* sr)
                  {
                    if(sr->mc.nnu == 0) return -1;
                    return sr->mc.nu[0].tgtZ;
                  });

  const Var kTrueX({"mc.nnu", "mc.nu.x"},
                  [](const caf::StandardRecord* sr)
                  {
                    if(sr->mc.nnu == 0) return -1.f;
                    return sr->mc.nu[0].x;
                  });

  const Var kTrueY({"mc.nnu", "mc.nu.y"},
                  [](const caf::StandardRecord* sr)
                  {
                    if(sr->mc.nnu == 0) return -1.f;
                    return sr->mc.nu[0].y;
                  });

  const Var kTrueQ2({"mc.nnu", "mc.nu.q2"},
                    [](const caf::StandardRecord* sr)
                    {
                      if(sr->mc.nnu == 0) return -1.f;
                      return sr->mc.nu[0].q2;
                    });

  const Var kTrueW2({"mc.nnu", "mc.nu.W2"},
                    [](const caf::StandardRecord* sr)
                    {
                      if(sr->mc.nnu == 0) return -1.f;
                      return sr->mc.nu[0].W2;
                    });

  struct Spects
  {
    Spectrum *E, *x, *y, *q2, *W2;
  };

  std::map<int, Spects> ssgen, ssgib;

  Spectrum sZgen("Target Z", Binning::Simple(30,0,30), lgenie, kTgtZ, kPreselAny);
  Spectrum sZgib("Target Z", Binning::Simple(30,0,30), lgibuu, kTgtZ, kPreselAny, kNoShift, kGiBUUWgt);

  Spectrum sWtgib("log_{10}(weight)", Binning::Simple(60,-2,+4), lgibuu, kLogWgt, kPreselAny, kNoShift, kGiBUUWgt);

  std::map<int, std::string> names = {{1, "H"},
                                      {6, "C"},
                                      {8, "O"},
                                      {17, "Cl"},
                                      {22, "Ti"},
                                      {26, "Fe"}};
  for(int Z: {1, 6, 8, 17, 22, 26}){
    const Cut kPresel = kPreselAny && (kTgtZ == Z);

    ssgen[Z].E = new Spectrum("True neutrino energy", Binning::Simple(40,0,10), lgenie, kTrueE, kPresel);
    ssgib[Z].E = new Spectrum("True neutrino energy", Binning::Simple(40,0,10), lgibuu, kTrueE, kPresel, kNoShift, kGiBUUWgt);

    ssgen[Z].x = new Spectrum("x", Binning::Simple(40,0,2), lgenie, kTrueX, kPresel);
    ssgib[Z].x = new Spectrum("x", Binning::Simple(40,0,2), lgibuu, kTrueX, kPresel, kNoShift, kGiBUUWgt);

    ssgen[Z].y = new Spectrum("y", Binning::Simple(50,0,1), lgenie, kTrueY, kPresel);
    ssgib[Z].y = new Spectrum("y", Binning::Simple(50,0,1), lgibuu, kTrueY, kPresel, kNoShift, kGiBUUWgt);

    ssgen[Z].q2 = new Spectrum("Q^{2}", Binning::Simple(40,0,2), lgenie, kTrueQ2, kPresel);
    ssgib[Z].q2 = new Spectrum("Q^{2}", Binning::Simple(40,0,2), lgibuu, kTrueQ2, kPresel, kNoShift, kGiBUUWgt);

    ssgen[Z].W2 = new Spectrum("W^{2}", Binning::Simple(40,0,2), lgenie, kTrueW2, kPresel);
    ssgib[Z].W2 = new Spectrum("W^{2}", Binning::Simple(40,0,2), lgibuu, kTrueW2, kPresel, kNoShift, kGiBUUWgt);
  } // end for Z

  lgenie.Go();
  lgibuu.Go();

  sZgen.ToTH1(1e20, kBlue)->Draw("hist");
  sZgib.ToTH1(1e20, kRed)->Draw("hist same");

  TLegend* leg = new TLegend(.725, .725, .875, .875);
  leg->SetFillStyle(0);
  leg->AddEntry(sZgen.ToTH1(1e20, kBlue), "GENIE", "l");
  leg->AddEntry(sZgib.ToTH1(1e20, kRed),  "GiBUU", "l");
  leg->Draw();

  gPad->Print("tgtZ.pdf");

  new TCanvas;

  sWtgib.ToTH1(1e20)->Draw("hist");
  gPad->Print("genweight.pdf");

  for(int Z: {1, 6, 8, 17, 22, 26}){
    new TCanvas;
    Compare(ssgen[Z].E, ssgib[Z].E, names[Z]);
    leg->Draw();
    gPad->Print(("Enu_"+names[Z]+".pdf").c_str());

    new TCanvas;
    Compare(ssgen[Z].x, ssgib[Z].x, names[Z]);
    leg->Draw();
    gPad->Print(("x_"+names[Z]+".pdf").c_str());

    new TCanvas;
    Compare(ssgen[Z].y, ssgib[Z].y, names[Z]);
    leg->Draw();
    gPad->Print(("y_"+names[Z]+".pdf").c_str());

    new TCanvas;
    Compare(ssgen[Z].q2, ssgib[Z].q2, names[Z]);
    leg->Draw();
    gPad->Print(("q2_"+names[Z]+".pdf").c_str());

    new TCanvas;

    Compare(ssgen[Z].W2, ssgib[Z].W2, names[Z]);
    leg->Draw();
    gPad->Print(("W2_"+names[Z]+".pdf").c_str());
  }
}

#endif