getTimePeakPlots.C

Go to the documentation of this file.

#pragma once

#include "CAFAna/Vars/GenieWeights.h"
#include "CAFAna/Vars/PPFXWeights.h"


#include "Includes.h"
#include "VarsAndCuts.h"
#include "CustomFunctions.h"
#include "GetSpectra.h"

using namespace ana;

void getTimePeakPlots(std::string samdef, unsigned int timestamp, unsigned int nDays, std::string naming)
{

  //this should make us more economic and efficient
  //only data spectra will be made every time
  bool makePred = 0;  // false because the predictions already exist.
  if(makePred){
    getPredictions(timestamp, nDays, naming);
  }
  bool makeData = 0;
  if(makeData){
    getData(timestamp, nDays, naming, samdef);
  }

  std::cout << "Getting files ..." << std::endl;
  std::string mcFile = "predictions_" + naming + ".root";
  std::string dtFile = "dataspectra_" + naming + ".root";



  TFile* mcPred      = new TFile(mcFile.c_str());
  TFile* dtSpec      = new TFile(dtFile.c_str());
  std::cout << "Setting osc pars" << std::endl;
  osc::OscCalcPMNSOpt calcMin;
  osc::OscCalcPMNSOpt calcMax;
  osc::OscCalcPMNSOpt calcDef;
  ana::ResetOscCalcToDefault(&calcMin);
  ana::ResetOscCalcToDefault(&calcMax);
  ana::ResetOscCalcToDefault(&calcDef);

  std::cout << "PASSED" << std::endl;
  std::cout << " "      << std::endl;

  calcMin.SetTh13(asin(sqrt(.084-0.008))/2); // (0.084) 0.076 to 0.092
  calcMin.SetTh12(asin(sqrt(.846+0.021))/2); // (0.846) 0.825 to 0.867
  calcMin.SetTh23(asin(sqrt(0.437))); // deltam2 > 0 -> sin2theta23 = 0.437 (PDG2016) 
  calcMin.SetDmsq21(7.71e-5); // (7.53) 7.35 to 7.71
  calcMin.SetDmsq32(2.46e-3); // (2.40) 2.34 to 2.46

  calcMax.SetTh13(asin(sqrt(.084+0.008))/2); // (0.084) 0.076 to 0.092
  calcMax.SetTh12(asin(sqrt(.846+0.021))/2); // (0.846) 0.825 to 0.867
  calcMax.SetTh23(asin(sqrt(0.569))); // deltam2 < 0 -> sin2theta23 = 0.569 (PDG2016)
  calcMax.SetDmsq21(7.35e-5); // (7.53) 7.35 to 7.71
  calcMax.SetDmsq32(2.34e-3); // (2.40) 2.34 to 2.46 
  

  std::cout << "Defining plots ..." << std::endl;
  Plot plots[kNumPlots] =
    {
      //      {"Slice time [#mus]", "tus", Binning::Simple(35,13,433), ana::kSliceTime},
      {"Slice time [#mus]", "tus", Binning::Simple(45,1,541), ana::kSliceTime}, 
      {"Run", "run", Binning::Simple(15000, 15000, 30000), SIMPLEVAR(hdr.run)},
      {"Time", "time", Binning::Simple(nDays<3 ? 24*nDays : nDays,timestamp - 86400*nDays,timestamp), SIMPLEVAR(hdr.unixtime)},
      {"Max y [m]", "maxy", kXYBins, kBoxMaxY},
      {"Slice time [#mus]", "tighttus", Binning::Simple(360,205,241), ana::kSliceTime},
      {"No. Kalman tracks", "nkal", Binning::Simple(10, 0, 10), SIMPLEVAR(trk.kalman.ntracks)},
      {"No. Hits", "nhit", Binning::Simple(10, 0, 500), kNHit},
      {"Kalman angle", "angkal", Binning::Simple(5, 0, 1), SIMPLEVAR(sel.cosrej.anglekal)},
      {"Z direction", "dirz", Binning::Simple(5, 0.5, 1), ana::kDirZ},
      {"BDT", "bdt", Binning::Simple(20, 0, 1), ana::kNumuContPID},
      {"p_{T}/p", "ptp", Binning::Simple(20, 0, 1), ana::kPtP}
    };
  


  std::cout << "Defining spectrums ..." << std::endl;
  Spectrum* hData[kNumPlots];
  Spectrum* hBkg[kNumPlots];
  PredictionNoExtrap* hPred[kNumPlots];


  for(int i = 0; i < kNumPlots; ++i){
    std::cout << "" << i << std::endl;
    Plot p = plots[i];
    std::string dataName = Form("data_%s", (p.fname).c_str());
    std::string bkgName =  Form("bkg_%s", (p.fname).c_str());
    std::cout << "Data name :" << dataName << " and the Bkg Name :" << bkgName << std::endl;
    hData[i] = LoadFrom<Spectrum>(dtSpec, dataName).release(); 
    hBkg[i] = LoadFrom<Spectrum>(dtSpec, bkgName).release();
    hPred[i] = LoadFrom<PredictionNoExtrap>(mcPred, (p.fname).c_str() ).release();
  }
  
  
  std::cout << "Getting miscellanius spectra ..." << std::endl;
  Spectrum* potTime   = LoadFrom<Spectrum>(dtSpec, "potTime").release(); 
  Spectrum* massDenom = LoadFrom<Spectrum>(dtSpec, "massDenom").release();
  Spectrum* massNum   = LoadFrom<Spectrum>(dtSpec, "massNum").release();
  Spectrum* nonSwapDenom = LoadFrom<Spectrum>(mcPred, "nonSwapDenom").release();
  Spectrum* swapDenom = LoadFrom<Spectrum>(mcPred, "swapDenom").release();

  const double massFactor = massNum->ToTH1(1)->Integral(0, -1) / massDenom->ToTH1(1)->Integral(0, -1);
  double signal = 0, mcsig = 0, mcsig_min = 0, mcsig_max = 0, bkg = 0;
  double mc_numu = 0, mc_nue = 0, mc_nc = 0;

  std::cout << "Effective mass: " << massFactor*14 << " diblocks" << std::endl;

  for(int i = 0; i < kNumPlots; ++i){
    new TCanvas;
    TH1* hSig = hData[i]->ToTH1(hData[i]->POT());
    hSig->SetName(("sig_"+plots[i].fname).c_str());
  
    TGraphAsymmErrors* grSig = GraphWithPoissonErrors(hSig, true, false);
   
    if(i==1 || i==2) // no errors in run or time plots
      hSig->Draw("hist");
    else
      {
        hSig->SetLineColor(0);
        hSig->SetLineWidth(1);
        hSig->SetMarkerColor(0);
        hSig->GetXaxis()->SetTitle((plots[i].label).c_str());
        hSig->GetYaxis()->SetRangeUser(0,hSig->GetMaximum()*1.12+4.);
        hSig->Draw("hist");
        grSig->SetMarkerStyle(20);
        grSig->SetMarkerColor(kBlack);
        grSig->SetLineColor(kBlack);
        grSig->SetLineWidth(2);
        grSig->SetName(("grSig_"+plots[i].fname).c_str());
        grSig->GetXaxis()->SetTitle((plots[i].label).c_str());
        grSig->Draw("AP");
      }

    TH1* hMC = hPred[i]->Predict(&calcDef).ToTH1(hData[i]->POT()*massFactor);
    TH1* hMC_min = hPred[i]->Predict(&calcMin).ToTH1(hData[i]->POT()*massFactor);
    hMC_min->SetName(("MC_min_"+plots[i].fname).c_str());
    hMC_min->SetLineColor(kRed-9);
    hMC_min->SetLineWidth(1);
    TH1* hMC_max = hPred[i]->Predict(&calcMax).ToTH1(hData[i]->POT()*massFactor);
    hMC_max->SetName(("MC_max_"+plots[i].fname).c_str());
    hMC_max->SetLineColor(kRed-9);
    hMC_max->SetLineWidth(1);

    TH1* hMC_numu = hPred[i]->PredictComponent(&calcDef, Flavors::kAllNuMu, Current::kCC, Sign::kBoth).ToTH1(hData[i]->POT()*massFactor);
    hMC_numu->SetLineColor(kBlue);
    TH1* hMC_nue  = hPred[i]->PredictComponent(&calcDef, Flavors::kAllNuE, Current::kCC, Sign::kBoth).ToTH1(hData[i]->POT()*massFactor);
    hMC_nue->SetLineColor(kGreen+2);
    TH1* hMC_NC   = hPred[i]->PredictComponent(&calcDef, Flavors::kAll, Current::kNC, Sign::kBoth).ToTH1(hData[i]->POT()*massFactor);
    hMC_NC->SetLineColor(kGray);

    TH1* hBKG = hBkg[i]->ToTH1(hData[i]->POT()*timeWindowFactor);
    hBKG->SetName(("BKG_"+plots[i].fname).c_str());
    hBKG->SetLineColor(kBlue);

    TH1* hExp = hPred[i]->Predict(&calcDef).ToTH1(hData[i]->POT()*massFactor);
    TH1* hExp_min = hPred[i]->Predict(&calcMin).ToTH1(hData[i]->POT()*massFactor);
    TH1* hExp_max = hPred[i]->Predict(&calcMax).ToTH1(hData[i]->POT()*massFactor);

    if(i == 0){
      // But for which we should print s/b stats
      signal = hSig->GetBinContent(19);
      mcsig = hPred[i]->Predict(&calcDef).ToTH1(hData[i]->POT()*massFactor)->Integral(0, -1);
      mc_numu = hMC_numu->Integral(0, -1);
      mc_nue = hMC_nue->Integral(0, -1);
      mc_nc = hMC_NC->Integral(0, -1);
      mcsig_min = hPred[i]->Predict(&calcMin).ToTH1(hData[i]->POT()*massFactor)->Integral(0, -1);
      mcsig_max = hPred[i]->Predict(&calcMax).ToTH1(hData[i]->POT()*massFactor)->Integral(0, -1);
      bkg = (hSig->Integral(1, hSig->GetNbinsX())-signal)/(hSig->GetNbinsX()-1);
      std::cout << "Signal: " << signal << std::endl;
      std::cout << "Bkg estimate: " << bkg << std::endl;
      std::cout << "MC signal min expectation: " << mcsig_min << std::endl;
      std::cout << "MC signal expectation: " << mcsig << std::endl;
      std::cout << "MC signal max expectation: " << mcsig_max << std::endl;
      std::cout << "MC before cuts " << (*nonSwapDenom + *swapDenom).ToTH1(hData[i]->POT())->Integral(0, -1)*massFactor << std::endl;
      std::cout << "MC expected numu : " << mc_numu << std::endl;
      std::cout << "MC expected nue  : " << mc_nue << std::endl;
      std::cout << "MC expected NC   : " << mc_nc << std::endl;

      hExp->SetBinContent(19,mcsig+bkg);
      hExp_min->SetBinContent(19,mcsig_min+bkg);
      hExp_max->SetBinContent(19,mcsig_max+bkg);

      TGraph* g = new TGraph;
      g->SetPoint(0, -1000e3, bkg);
      g->SetPoint(1, +1000e3, bkg);
      g->SetLineColor(kBlue);
      g->SetLineWidth(2);
      g->SetLineStyle(7);
      g->Draw("L");
    }
    //for time so we can reuse the predictions
    //first clone BKG so the bins coincide
    if(i==2){
      TH1* hHelp = (TH1*)hBKG->Clone();
      TH1* hHelpMin = (TH1*)hBKG->Clone();
      TH1* hHelpMax = (TH1*)hBKG->Clone();

      for(int binId=0; binId<hExp->GetNbinsX(); binId++){
        hHelp->SetBinContent(binId+1, hExp->GetBinContent(binId+1));
        hHelpMin->SetBinContent(binId+1, hExp_min->GetBinContent(binId+1));
        hHelpMax->SetBinContent(binId+1, hExp_max->GetBinContent(binId+1));
      }

      hExp->Clear();
      hExp = (TH1*)hHelp->Clone();
      hExp->Add(hBKG);

      hExp_min->Clear();
      hExp_min = (TH1*)hHelp->Clone();
      hExp_min->Add(hBKG);

      hExp_max->Clear();
      hExp_max = (TH1*)hHelp->Clone();
      hExp_max->Add(hBKG);
    }
    else
    {
      hExp->Add(hBKG);
      hExp_min->Add(hBKG);
      hExp_max->Add(hBKG);
    }

    hExp_max->SetName(("Exp_max_"+plots[i].fname).c_str());
    hExp_max->SetLineColor(kRed-9);
    hExp_min->SetName(("Exp_min_"+plots[i].fname).c_str());
    hExp_max->SetLineColor(kRed-9);

    TGraph* grShade = ShadeBetweenHistograms(hExp_min, hExp_max);
    grShade->SetName(("shade_"+plots[i].fname).c_str());
    grShade->SetFillColor(kRed-9);
    grShade->SetLineColor(kRed-9);
    grShade->Draw("lf");
  
    double ymax = 1.3* std::max( hSig->GetMaximum(), hExp->GetMaximum());
    hExp->GetYaxis()->SetRangeUser(0,ymax);
    gPad->Update();
    hExp->GetXaxis()->SetTitle((plots[i].label).c_str());
    hExp->SetName(("Exp_"+plots[i].fname).c_str());
    hExp->SetLineColor(kRed);
    hExp->SetLineStyle(1);
    hExp->Draw("hist same");
    hMC_numu->Draw("hist same");
    hMC_nue->Draw("hist same");
    hMC_NC->Draw("hist same");
    hExp->Draw("hist same");
    if(i!=1 && i!=2){
      hSig->Draw("hist same");
      grSig->Draw("P same");
    }

    TLegend* ll = new TLegend(0.15,0.60,0.45,0.89);
    ll->AddEntry(grSig,"Data","l");
    ll->AddEntry(grShade,"Osc. Uncty.","lf");
    ll->AddEntry(hExp,"Expectation","l");
    ll->AddEntry(hMC_numu,"#nu_{#mu}","l");
    ll->AddEntry(hMC_nue,"#nu_{e}","l");
    ll->AddEntry(hMC_NC,"NC","l");
    ll->SetFillColor(0);
    ll->Draw();


    hSig->SaveAs(("/nova/ana/users/novadq/TimePeak/Plots/" + naming + "/Tools/sig_"+plots[i].fname+".root").c_str());
    hMC->SaveAs(("/nova/ana/users/novadq/TimePeak/Plots/" + naming + "/Tools/MC_"+plots[i].fname+".root").c_str());
    hBKG->SaveAs(("/nova/ana/users/novadq/TimePeak/Plots/" + naming + "/Tools/BKG_"+plots[i].fname+".root").c_str());

    lastUpdated(timestamp);
    gPad->SaveAs(("/nova/ana/users/novadq/TimePeak/Plots/" + naming + "/Plots/"+plots[i].fname+".png").c_str());

    if(i == 1){
      std::cout << "Bad runs:" << std::endl;
      TH1* h = hData[i]->ToTH1(hData[i]->POT());
      for(int j = 1; j <= h->GetNbinsX(); ++j){
        if(h->GetBinContent(j) > 5)
          std::cout << h->GetXaxis()->GetBinLowEdge(j) << std::endl;
      }
    }    
  } // End loop over i


  TH2* hPOT = potTime->ToTH2(mcsig);
  hPOT->SetName("hPOT");
  hPOT->SaveAs(("/nova/ana/users/novadq/TimePeak/Plots/"+naming+"/Tools/hPOT.root").c_str());
  if(naming == "All") 
    getSigmaPlotsAll(mcsig, timestamp, nDays, naming);
  else 
    getSigmaPlots(mcsig, timestamp, nDays, naming);

}