compareEvents.C

Go to the documentation of this file.

#include <string>
#include <set>
#include <list>
#include <unordered_map>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <math.h>

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

struct EventID{

  EventID()
      : run(0)
      , subrun(0)
      , event(0)
      , slice(0)
      , cycle(0)
      , energy(0.)
      , cvn(0.)
      , wgt(0.)
      , select("")
  {}

  EventID(long r,
          long sr,
          long ev,
          long sl,
          long c,
          double en,
          double id,
          double weight,
          std::string sel)
      : run(r)
      , subrun(sr)
      , event(ev)
      , slice(sl)
      , cycle(c)
      , energy(en)
      , cvn(id)
      , wgt(weight)
      , select(sel)
  {}

  long   run   ;
  long   subrun;
  long   event ;
  long   slice ;
  long   cycle ;
  double energy;
  double cvn;
  double wgt;
  std::string select;

  bool operator<(EventID const& other) const {
    if(run < other.run   ) return true;
    else{
      if(run == other.run &&
         subrun < other.subrun) return true;
      else{
        if(run == other.run &&
           subrun == other.subrun &&
           event < other.event ) return true;
        else {
          if(run == other.run &&
             subrun == other.subrun &&
             event == other.event &&
             slice  < other.slice ) return true;
          else{
            if(run == other.run &&
               subrun == other.subrun &&
               event == other.event &&
               slice  == other.slice &&
               cycle  < other.cycle ) return true;
          }
        }
      }
    }
    return false;
  }

};

//-----------------------------------------------------------------------------
void fillEventIDsFromFile(std::string       const& inputName,
                          std::set<EventID>      & eventList,
                          std::string       const& selection)
{
  std::cout << inputName
            << std::endl;

  // read in the event lists list of events
  std::ifstream eventListFile(inputName);

  std::cout << inputName
            << " is open "
            << eventListFile.is_open()
            << std::endl;

  double energy, cvn, wgt;
  double run, subrun, event, slice, cycle;

  long ev;

  int highEventNumCnt = 0;

  std::string select;
  std::string selectToUse = "NuMu";
  if     (selection == "nue") selectToUse = "NuE";
  else if(selection == "nus") selectToUse = "NC";

  while( !eventListFile.eof() && eventListFile.is_open() ){

    if(inputName.find("caf") == std::string::npos){
      eventListFile >> select >> run >> subrun >> event >> slice >> cycle >> energy >> cvn >> wgt;
      if(select.find(selectToUse) == std::string::npos) continue;
    }
    else{
      eventListFile >> run >> subrun >> event >> slice >> cycle >> energy >> cvn >> wgt;
      select = "caf_" + selection;
    }

    if(event > 1.e6){
      ev = event / 10;
      ev *= 10.;
      ++highEventNumCnt;
    }
    else ev = event;

    EventID evId((long)run,
                 (long)subrun,
                 (long)ev,
                 (long)slice,
                 (long)cycle,
                 energy,
                 cvn,
                 wgt,
                 select);

    if(eventList.size() < 10)
      std::cout << select
                << " "
                << evId.run
                << " "
                << evId.subrun
                << " "
                << evId.event
                << " / "
                << event
                << " "
                << evId.slice
                << " "
                << evId.cycle
                << " "
                << eventList.size()
                << std::endl;

    eventList.insert(evId);

  } // end loop to load event lists events

  std::cout << "There are "
            << eventList.size()
            << " events from "
            << inputName
            << " and "
            << highEventNumCnt
            << " high event numbers"
            << std::endl;

}

//-----------------------------------------------------------------------------
int findMissingEvents(bool firstCAF,
                      std::string const& missingListName,
                      std::set<EventID> const& firstList,
                      std::set<EventID> const& secondList,
                      TH1D* energyDiff,
                      TH1D* cvnDiff,
                      TH1D* wgtDiff,
                      std::unordered_map<std::string, TH1D*> & weightedSpect,
                      std::unordered_map<std::string, TH1D*> & unweightedSpect)
{
  int cntMissing = 0;

  std::ofstream missingList(missingListName);

  for(auto const& itr : firstList){
    if(secondList.find(itr) == secondList.end()){
      ++cntMissing;

      missingList << itr.run
                  << " "
                  << itr.subrun
                  << " "
                  << itr.event
                  << " "
                  << itr.slice
                  << " "
                  << itr.cycle
                  << std::endl;
    }
    else{
      if(!firstCAF){
        energyDiff                 ->Fill(itr.energy - secondList.find(itr)->energy);
        cvnDiff                    ->Fill(itr.cvn - secondList.find(itr)->cvn);
        wgtDiff                    ->Fill(itr.wgt - secondList.find(itr)->wgt);
        weightedSpect[itr.select]  ->Fill(itr.energy, itr.wgt * 1.e-3);
        unweightedSpect[itr.select]->Fill(itr.energy, 1.e-3);
      }
    } // end we found the event in the second list
  } // end loop to look for missing caf evnets

  missingList.close();

  return cntMissing;
}

//-----------------------------------------------------------------------------
void drawHistograms(std::string const& fileName,
                    TH1D* energyDiff,
                    TH1D* cvnDiff,
                    TH1D* wgtDiff,
                    TH1D* missing,
                    std::unordered_map<std::string, TH1D*> weightedSpect,
                    std::unordered_map<std::string, TH1D*> unweightedSpect,
                    std::unordered_map<std::string, TH1D*> spectRatio)
{
  // set how many pads and which pad is which
  // plot now
  int xPads       = 4;
  int yPads       = 3;
  int countPad    = 1;
  int deltaEPad   = 2;
  int deltaWgtPad = 3;
  int deltaCVNPad = 4;
  int spectPadBeg = 5;
  int ratioPadBeg = 9;

  if(fileName.find("nus") != std::string::npos){
    xPads       = 3;
    yPads       = 2;
    spectPadBeg = 4;
    ratioPadBeg = 5;
  }

  TCanvas *canv = new TCanvas("compareCanv",
                              "Compare CMF vs CAF",
                              300 * xPads,
                              300 * yPads);

  // the number of divisions is equal to
  // the number of selections in the x and 4 in y
  canv->Divide(xPads, yPads);

  canv->cd(countPad)->SetLogy();
  missing->Draw();

  canv->cd(deltaEPad)->SetLogy();
  energyDiff->GetXaxis()->SetRangeUser(-0.5, 0.5);
  energyDiff->Draw();

  canv->cd(deltaWgtPad)->SetLogy();
  wgtDiff->GetXaxis()->SetRangeUser(-0.5, 0.5);
  wgtDiff->Draw();

  if(fileName.find("nue") != std::string::npos){
    canv->cd(deltaCVNPad)->SetLogy();
    cvnDiff->GetXaxis()->SetRangeUser(-0.5, 0.5);
    cvnDiff->Draw();
  }

  size_t i = 0;
  for(auto const& itr : weightedSpect){
    canv->cd(i + spectPadBeg);

    itr.second->SetMaximum(std::max(itr.second->GetMaximum(), unweightedSpect[itr.first]->GetMaximum()) * 1.75);
    itr.second->Draw();
    unweightedSpect[itr.first]->Draw("same");
    TLegend *leg = new TLegend(0.5, 0.5, 0.8, 0.8);
    leg->AddEntry(itr.second, itr.first.c_str(), "");
    leg->AddEntry(itr.second, "Weighted Spectrum", "l");
    leg->AddEntry(unweightedSpect[itr.first], "Unweighted Spectrum", "l");
    leg->SetFillStyle(0);
    leg->SetBorderSize(0);
    leg->Draw();

    canv->cd(i + ratioPadBeg)->SetGridy();
    spectRatio[itr.first]->GetYaxis()->SetRangeUser(0.5, 1.5);
    spectRatio[itr.first]->Draw();

    ++i;
  } // end loop to draw spectra

  canv->Update();
  canv->Write();
  canv->Print(fileName.c_str(), "pdf");

  energyDiff->Write();
  cvnDiff   ->Write();
  missing   ->Write();

  for(auto const& itr : weightedSpect){
    itr.second                ->Write();
    unweightedSpect[itr.first]->Write();
    spectRatio     [itr.first]->Write();
  }

}

//-----------------------------------------------------------------------------
void compareEvents(std::string det="DETECTOR",
                   std::string sel="SELECT",
                   std::string sample="FILETYPE",
                   std::string horncur="HORNCUR",
                   std::string tag="TAG",
                   std::string syst="SYST")
{
  if( !syst.empty() ) syst.insert(0, "_");

  std::set<EventID> eventListEvents;
  std::set<EventID> cafEvents;

  TFile *file = new TFile(("compare_hist_" + det + "_" + sel + "_" + horncur + "_" + sample + syst + "_" + tag + ".root").c_str(), "RECREATE");

  TH1D* energyDiff = new TH1D("energyDiff", ";#Delta E (GeV);Events", 2000, -1.e-4, 1.e-4);
  energyDiff->GetXaxis()->CenterTitle();
  energyDiff->GetYaxis()->CenterTitle();

  TH1D* cvnDiff = new TH1D("cvnDiff", ";#Delta CVN;Events", 2000, -1.e-4, 1.e-4);
  cvnDiff->GetXaxis()->CenterTitle();
  cvnDiff->GetYaxis()->CenterTitle();

  TH1D* wgtDiff = new TH1D("wgtDiff", ";#Delta (PPFX #times XSecCV);Events", 2000, -1.e-4, 1.e-4);
  wgtDiff->GetXaxis()->CenterTitle();
  wgtDiff->GetYaxis()->CenterTitle();

  // Histogram to visualize the number of missing events and totals
  TH1D* missing = new TH1D("missing", ";Missing in CAF or CMF;Events", 4, 0., 4.);
  missing->GetXaxis()->SetBinLabel(1, "Total CMF");
  missing->GetXaxis()->SetBinLabel(2, "CMF not in CAF");
  missing->GetXaxis()->SetBinLabel(3, "Total CAF");
  missing->GetXaxis()->SetBinLabel(4, "CAF not in CMF");
  missing->GetXaxis()->CenterTitle();
  missing->GetYaxis()->CenterTitle();

  std::list<std::string> sels;
  std::vector<double> bins;
  std::vector<double> nueBins ({1., 1.5, 2., 2.5, 3., 3.5, 4., 4.5, 5.});
  std::vector<double> numuBins({0., 0.75, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2., 2.25, 2.5, 2.75, 3., 3.5, 4., 5.});
  std::vector<double> nusBins;
  for(double b = 0.75; b < 20.25; b += 0.25) nusBins.push_back(b);

  if(sel == "numu"){
    bins.swap(numuBins);
    sels.emplace_back("NuMuSelQ1");
    sels.emplace_back("NuMuSelQ2");
    sels.emplace_back("NuMuSelQ3");
    sels.emplace_back("NuMuSelQ4");
  }
  else if(sel == "nue"){
    bins.swap(nueBins);
    sels.emplace_back("NuESel_LowPID");
    sels.emplace_back("NuESel_HighPID");
    sels.emplace_back("NuESel_Peripheral");
  }
  else if(sel == "nus"){
    bins.swap(nusBins);
    sels.emplace_back("NCSel");
  }

  std::unordered_map<std::string, TH1D*> weightedSpect;
  std::unordered_map<std::string, TH1D*> unweightedSpect;
  std::unordered_map<std::string, TH1D*> spectRatio;

  for(auto const& itr : sels){
    //std::cout << itr << std::endl;

    weightedSpect[itr] = new TH1D(("weighted" + itr).c_str(),
                                  ";E (GeV);Events / GeV (#times 1000)",
                                  bins.size() - 1,
                                  &bins[0]);
    weightedSpect[itr]->GetXaxis()->CenterTitle();
    weightedSpect[itr]->GetYaxis()->CenterTitle();
    weightedSpect[itr]->GetXaxis()->SetDecimals();
    weightedSpect[itr]->GetYaxis()->SetDecimals();
    weightedSpect[itr]->Sumw2();

    unweightedSpect[itr] = new TH1D(("unweighted" + itr).c_str(),
                                    ";E (GeV);Events / GeV (#times 1000)",
                                    bins.size() - 1,
                                    &bins[0]);
    unweightedSpect[itr]->GetXaxis()->CenterTitle();
    unweightedSpect[itr]->GetYaxis()->CenterTitle();
    unweightedSpect[itr]->GetXaxis()->SetDecimals();
    unweightedSpect[itr]->GetYaxis()->SetDecimals();
    unweightedSpect[itr]->Sumw2();
    unweightedSpect[itr]->SetLineColor(2);
    unweightedSpect[itr]->SetMarkerColor(2);

    spectRatio[itr] = new TH1D(("spectRatio" + itr).c_str(),
                               ";E (GeV);Weighted / Unweighted",
                               bins.size() - 1,
                               &bins[0]);
    spectRatio[itr]->GetXaxis()->CenterTitle();
    spectRatio[itr]->GetYaxis()->CenterTitle();
    spectRatio[itr]->GetXaxis()->SetDecimals();
    spectRatio[itr]->GetYaxis()->SetDecimals();
    spectRatio[itr]->SetMarkerStyle(22);
    spectRatio[itr]->SetMarkerSize(0.8);
  } // end loop to make spectra for each selection

  std::string inputName("selectedEvents_" + det + "_" + sel + "_" + horncur + "_" + sample + syst + "_" + tag + ".txt");

  fillEventIDsFromFile(inputName,
                       eventListEvents,
                       sel);

  inputName = "cafEvents_" + det + "_" + sel + "_" + horncur + "_" + sample + syst + "_" + tag + ".txt";

  size_t pos = inputName.find("-");
  while( pos != std::string::npos){
    inputName.replace(pos, 1, "_");
    pos = inputName.find("-");
  }

  fillEventIDsFromFile(inputName, cafEvents, sel);

  // loop over the event list and look for missing caf events
  int cntMissingCAF = findMissingEvents(false,
                                        "missing_caf_list_" + det + "_" + sel + "_" + horncur + "_" + sample + syst + "_" + tag + ".txt",
                                        eventListEvents,
                                        cafEvents,
                                        energyDiff,
                                        cvnDiff,
                                        wgtDiff,
                                        weightedSpect,
                                        unweightedSpect);

  int cntMissingEventList = findMissingEvents(true,
                                              "missing_event_list_" + det + "_" + sel + "_" + horncur + "_" + sample + syst + "_" + tag + ".txt",
                                              cafEvents,
                                              eventListEvents,
                                              energyDiff,
                                              cvnDiff,
                                              wgtDiff,
                                              weightedSpect,
                                              unweightedSpect);

  std::cout << "There are "
            << cntMissingCAF
            << " event list events not in the caf events and "
            << cntMissingEventList
            << " caf events not in the event list events"
            << std::endl;

  // make the ratios of the weighted and unweighted histograms
  double norm;
  for(auto const& itr : weightedSpect){
    spectRatio[itr.first]->Divide(weightedSpect[itr.first], unweightedSpect[itr.first]);

    // normalize the bin contents to be per 0.1 GeV
    for(int i = 1; i < itr.second->GetNbinsX() + 1; ++i){
      norm = 1. / itr.second->GetBinWidth(i);
      itr.second->SetBinContent(i, itr.second->GetBinContent(i) * norm);
      itr.second->SetBinError(i, itr.second->GetBinError(i) * norm);

      unweightedSpect[itr.first]->SetBinContent(i, unweightedSpect[itr.first]->GetBinContent(i) * norm);
      unweightedSpect[itr.first]->SetBinError(i, unweightedSpect[itr.first]->GetBinError(i) * norm);
    }
  }

  missing->Fill(0.5, eventListEvents.size());
  missing->Fill(1.5, cntMissingCAF);
  missing->Fill(2.5, cafEvents.size());
  missing->Fill(3.5, cntMissingEventList);

  drawHistograms("compare_" + det + "_" + sel + "_" + horncur + "_" + sample + "_" + tag + ".pdf",
                 energyDiff,
                 cvnDiff,
                 wgtDiff,
                 missing,
                 weightedSpect,
                 unweightedSpect,
                 spectRatio);

  file->Write();
  file->Close();
}