compareTextLists.C

Go to the documentation of this file.

/* this macro takes the cmf and caf lists of events
 * and does a few things
 * -- it first checks that all events in caf are present in cmf and vice versa
 * -- it then produces a tree with any differences between 
 *    ereco, etrue, cvweights, pdgcodes and CCNC values
 */


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

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

struct EventInfo{

  EventInfo()
    : lepe(0.0)
    , hade(0.0)
    , ecalo(0.0)
    , ereco(0.)
    , etrue(0.)
    , pdg(0.)
    , ccnc(0.)
    , cvwgt(0.)
  {}

  EventInfo(double le,
            double he,
            double ec,
            double er,
            double et,
            double p,
            double cc,
            double ew)
    : lepe(le)
    , hade(he)
    , ecalo(ec)
    , ereco(er)
    , etrue(et)
    , pdg(p)
    , ccnc(cc)
    , cvwgt(ew)
  {}

  double lepe;
  double hade;
  double ecalo;
  double ereco;
  double etrue;
  double pdg;
  double ccnc;
  double cvwgt;

  // assume that ereco and cvwgt are the most likely
  // offenders so put them first
  bool operator<(EventInfo const& other) const {
    if(ereco < other.ereco   ) return true;
    else{
      if(ereco == other.ereco &&
         cvwgt < other.cvwgt) return true;
      else{
        if(ereco == other.ereco &&
           cvwgt == other.cvwgt &&
           etrue < other.etrue) return true;
        else{
          if(ereco == other.ereco &&
             cvwgt == other.cvwgt &&
             etrue == other.etrue &&
             pdg < other.pdg) return true;
          else{
            if(ereco == other.ereco &&
               cvwgt == other.cvwgt &&
               etrue == other.etrue &&
               pdg   == other.pdg   &&
               ccnc  < other.ccnc) return true;
          }
        }
      }
    }
    return false;
  }

  bool operator==(EventInfo const& other) const{
    if (ereco == other.ereco && 
        cvwgt == other.cvwgt &&
        etrue == other.etrue) return true;
    else return false;
  }
};

struct EventID{

  EventID()
      : run(0)
      , subrun(0)
      , event(0)
      , slice(0)
      , cycle(0)
      
  {}

  EventID(long r,
          long sr,
          long ev,
          long sl,
          long c)
      : run(r)
      , subrun(sr)
      , event(ev)
      , slice(sl)
      , cycle(c)
  {}

  long   run   ;
  long   subrun;
  long   event ;
  long   slice ;
  long   cycle ;

  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;
  }
};

// just a simple struct for holding the name and eventMap together
struct EventContainer{
  std::map<EventID, EventInfo> eventMap;
  std::string name;
};

void fillEventMapFromFile(EventContainer   & eventContainer,
                          std::string const& fileName)
{

  // read in the event lists list of events
  std::ifstream eventLists(fileName);

  std::cout << fileName
            << " is open "
            << eventLists.is_open()
            << std::endl;

  double run, subrun, event, slice, cycle;
  double lepe, hade, ecalo, ereco, etrue, pdg, ccnc, cvwgt;

  long ev;

  while( !eventLists.eof() && eventLists.is_open() ){
    eventLists >> run >> subrun >> event >> slice >> cycle 
               >> lepe >> hade >> ecalo >> ereco >> etrue >> pdg >> ccnc >> cvwgt;

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

    EventID evId((long)run,
                 (long)subrun,
                 (long)ev,
                 (long)slice,
                 (long)cycle);

    EventInfo evIn(lepe, hade, ecalo, ereco, etrue, pdg, ccnc, cvwgt);

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

    eventContainer.eventMap[evId] = evIn;

  } // end loop to load event lists events

  std::cout << "There are "
            << eventContainer.eventMap.size()
            << " events "
            << std::endl;
}

EventContainer findMissingFiles(EventContainer &baseSample, 
                                EventContainer &searchSample)
{

  EventContainer differentEventSample;
  differentEventSample.name = baseSample.name + "NotIn" + searchSample.name;

  std::map<EventID, EventInfo> baseMap   = baseSample.eventMap;
  std::map<EventID, EventInfo> searchMap = searchSample.eventMap;

  int nEvents=0;

  for (auto& itr : baseMap){
    
    // first, check if the event exists in the searchMap
    if (searchMap.find(itr.first) != searchMap.end()){
      // if it exists, check that the information is correct
      if (itr.second == searchMap[itr.first]){
      }
      else{
        // event exists, but information is wrong
        differentEventSample.eventMap[itr.first] = itr.second;
      }

    }
    else{
      // event doesn't exist
      std::cout 
        << "Event " 
        << itr.first.run << " " 
        << itr.first.subrun << " "
        << itr.first.event << " "
        << itr.first.slice << " "
        << itr.first.cycle 
        << " exists in "
        << baseSample.name
        << " but not in "
        << searchSample.name
        << ". Printing information."
        << "\n ereco " << itr.second.ereco
        << "\n etrue " << itr.second.etrue
        << "\n pdg   " << itr.second.pdg
        << "\n ccnc  " << itr.second.ccnc
        << "\n cvwgt " << itr.second.cvwgt
        << std::endl;
      ++nEvents;
    }
  }

  std::cout
    << "there are " 
    << nEvents 
    << " events in base sample not in search sample"
    << std::endl;

  return differentEventSample;
}

void saveToFile(EventContainer cmf, EventContainer caf){

  // setup output histograms, canvas
  TFile* eventDifferences = new TFile("eventDifferences.root", "recreate");
  TTree* t = new TTree("ed", "Event Differences");
  int    run;
  int    subrun;
  int    event;
  int    slice;
  int    cycle;
  double cmf_lepe;
  double cmf_hade;
  double cmf_ecalo;
  double cmf_ereco;
  double cmf_etrue;
  double cmf_pdg;
  double cmf_ccnc;
  double cmf_cvwgt;
  double caf_lepe;
  double caf_hade;
  double caf_ecalo;
  double caf_ereco;
  double caf_etrue;
  double caf_pdg;
  double caf_ccnc;
  double caf_cvwgt;
  t->Branch("run"       , &run);
  t->Branch("subrun"    , &subrun);
  t->Branch("event"     , &event);
  t->Branch("slice"     , &slice);
  t->Branch("cycle"     , &cycle);
  t->Branch("cmf_lepe"  , &cmf_lepe);
  t->Branch("cmf_hade"  , &cmf_hade);
  t->Branch("cmf_ecalo" , &cmf_ecalo);
  t->Branch("cmf_ereco" , &cmf_ereco);
  t->Branch("cmf_etrue" , &cmf_etrue);
  t->Branch("cmf_pdg"   , &cmf_pdg);
  t->Branch("cmf_ccnc"  , &cmf_ccnc);
  t->Branch("cmf_cvwgt" , &cmf_cvwgt);
  t->Branch("caf_lepe"  , &caf_lepe);
  t->Branch("caf_hade"  , &caf_hade);
  t->Branch("caf_ecalo" , &caf_ecalo);
  t->Branch("caf_ereco" , &caf_ereco);
  t->Branch("caf_etrue" , &caf_etrue);
  t->Branch("caf_pdg"   , &caf_pdg);
  t->Branch("caf_ccnc"  , &caf_ccnc);
  t->Branch("caf_cvwgt" , &caf_cvwgt);

  for (auto& itr : cmf.eventMap){
    run    = itr.first.run;
    subrun = itr.first.subrun;
    event  = itr.first.event;
    cycle  = itr.first.cycle;
    slice  = itr.first.slice;

    cmf_lepe  = itr.second.lepe;
    cmf_hade  = itr.second.hade;
    cmf_ecalo = itr.second.ecalo;
    cmf_ereco = itr.second.ereco;
    cmf_etrue = itr.second.etrue;
    cmf_pdg   = itr.second.pdg;
    cmf_ccnc  = itr.second.ccnc;
    cmf_cvwgt = itr.second.cvwgt;

    caf_lepe  = caf.eventMap[itr.first].lepe;
    caf_hade  = caf.eventMap[itr.first].hade;
    caf_ecalo = caf.eventMap[itr.first].ecalo;
    caf_ereco = caf.eventMap[itr.first].ereco;
    caf_etrue = caf.eventMap[itr.first].etrue;
    caf_pdg   = caf.eventMap[itr.first].pdg;
    caf_ccnc  = caf.eventMap[itr.first].ccnc;
    caf_cvwgt = caf.eventMap[itr.first].cvwgt;

    t->Fill();
  }
  t->Write();
}

void compareTextLists(std::string cmfFile, std::string cafFile)
{

  EventContainer cmfEvents;
  EventContainer cafEvents;

  cmfEvents.name = "cmf";
  cafEvents.name = "caf";

  std::cout 
    << "Configured with "
    << "\n -- cmf file: " << cmfFile
    << "\n -- caf file: " << cafFile
    << std::endl;

  fillEventMapFromFile(cmfEvents, cmfFile);
  fillEventMapFromFile(cafEvents, cafFile);

  EventContainer CMFEventDiff = findMissingFiles(cmfEvents, cafEvents);
  EventContainer CAFEventDiff = findMissingFiles(cafEvents, cmfEvents);

  saveToFile(CMFEventDiff, CAFEventDiff);

}