timingPeak.C

Go to the documentation of this file.

#ifdef __CINT__
void timingPeak(std::string samdef, unsigned int timestamp, unsigned int nDays, std::string naming){
  std::cout << "Sorry, you must run in compiled mode" << std::endl;
}
#else

#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Analysis/Plots.h"
#include "CAFAna/Core/EventList.h"
#include "CAFAna/Core/ReweightableSpectrum.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Cuts/Cuts.h"
#include "CAFAna/Prediction/PredictionNoExtrap.h"
#include "CAFAna/Vars/Vars.h"
#include "CustomFunctions.h"
#include "VarsAndCuts.h"

using namespace ana;

#include "TCanvas.h"
#include "TF2.h"
#include "TFile.h"
#include "TGraph.h"
#include "TGraphAsymmErrors.h"
#include "TGraphErrors.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TLatex.h"
#include "TLegend.h"
#include "TMath.h"
#include "TPad.h"
#include "TPaveLabel.h"
#include "TProfile.h"
#include "TStyle.h"
#include "TTimeStamp.h"

void lastUpdated(unsigned int timestamp);
                                                                                                     
void timingPeak(std::string samdef, unsigned int timestamp, unsigned int nDays, std::string naming)
{
  const std::string fFiles = "defname: " + samdef + Form(" and online.subrunendtime >= '%u' and online.subrunendtime <= '%u'", timestamp - 86400*nDays, timestamp);
  SpectrumLoader loader(fFiles);

  // Hacky way of getting the livetime. Do as I say...
  TH1D* hLiveTime = new TH1D("hLiveTime",";Livetime (s)", 1, 0, 1);
  loader.AddSpillHistogram(hLiveTime, kSpillLivetime, kNoSpillCut, kSpillUnweighted);

  struct Plot
  {
    std::string label;
    std::string fname;
    Binning bins;
    Var var;
  };

  const int kNumPlots = 1;
  Plot plots[kNumPlots] =
    {
      {"Slice time [#mus]", "tus", Binning::Simple(45,25,475), kSliceTime}
    };

  Spectrum* hData[kNumPlots];
  Spectrum* hBkg[kNumPlots];

  for(int i = 0; i < kNumPlots; ++i){
    Plot p = plots[i];
    // Don't apply time cut for timing and runs plots, do for everything else
    hData[i] = new Spectrum(p.label, p.bins, loader, p.var, kCut);
  }

  loader.Go();

  // MakeEventListFile(fFiles, kCut && kInTimeTight, "eventlists/eventlist_"+naming+Form("_%u.txt",timestamp), true, true);
  // MakeEventListFile(fFiles, kCut && !kInTimeTight, Form("day/background/bkg_evtlist_%u.txt",timestamp), true);

  for(int i = 0; i < kNumPlots; ++i){

    hData[0]->OverrideLivetime(hLiveTime->Integral());

    std::cout << "Calculated livetime: " << hData[0]->Livetime() << std::endl;

    new TCanvas;
    TH1* hSig = hData[i]->ToTH1(hData[i]->Livetime(), kLivetime);
    hSig->SetName(("sig_"+plots[i].fname).c_str());
    // hSig->GetYaxis()->SetNdivisions(0);

    // TGraphAsymmErrors* grSig = GraphWithPoissonErrors(hSig, true, false);
    
    hSig->SetLineColor(kBlue);
    hSig->SetLineWidth(1);
    hSig->SetMarkerColor(kBlue);
    // hSig->GetYaxis()->SetRangeUser(0,hSig->GetMaximum()*1.12+4.);
    hSig->Draw("hist");
    // grSig->SetMarkerColor(kBlack);
    // grSig->SetLineColor(kBlack);
    // grSig->SetLineWidth(2);
    // grSig->SetName(("grSig_"+plots[i].fname).c_str());
    // grSig->Draw("P");
    
    int fSignal = (int) hSig->GetBinContent(30);
    double fBKG = (hSig->Integral() - fSignal)*1./(44);
    double prob = TMath::PoissonI(fSignal,fBKG);

    TLatex p1(0.15,0.80,Form("Evts. in time peak: %i",fSignal));
    TLatex p2(0.15,0.75,Form("BKG. estimate: %6.2g", fBKG));
    TLatex p3(0.15,0.70,Form("Probability: %6.2g", prob));

    p1.SetNDC();
    p1.SetTextSize(0.04);
    p1.Draw();
    p2.SetNDC();
    p2.SetTextSize(0.04);
    p2.Draw();
    p3.SetNDC();
    p3.SetTextSize(0.04);
    p3.Draw();

    TGraph g1,g2,g3;
    g1.SetPoint(0,315,0);
    g1.SetPoint(1,315,1e99);
    g2.SetPoint(0,325,0);
    g2.SetPoint(1,325,1e99);
    g3.SetPoint(0,0,fBKG);
    g3.SetPoint(1,500,fBKG);

    g1.SetLineColor(kRed);
    g2.SetLineColor(kRed);
    g3.SetLineColor(kPink);
    g1.SetLineStyle(7);
    g2.SetLineStyle(7);
    g3.SetLineStyle(2);
    g1.Draw("L");
    g2.Draw("L");
    g3.Draw("L");

    lastUpdated(timestamp);
    gPad->SaveAs(("/nova/ana/users/novadq/TimePeak/BNB/plots/" + plots[i].fname + "_" + naming + ".root").c_str());
    gPad->SaveAs(("/nova/ana/users/novadq/TimePeak/BNB/plots/" + plots[i].fname + "_" + naming + ".png").c_str());
    
  } // End loop over i

}

#endif