ExposureSource.cxx

Go to the documentation of this file.

#include "CAFAna/Core/ExposureSource.h"

#include "CAFAna/Core/IExposureSink.h"

#include "StandardRecord/Proxy/SRProxy.h"

#include <iostream>
#include <map>

namespace ana
{
  // --------------------------------------------------------------------------
  ExposureSource::~ExposureSource()
  {
    if(!fFinalized && !fCuts.empty()){
      std::cout << "ExposureSource destroyed without Finalize() being called" << std::endl;
      abort();
    }
  }

  // --------------------------------------------------------------------------
  void ExposureSource::AddSink(const Cut& cut, IExposureSink* s)
  {
    assert(!fFinalized);

    // It's more important that we can iterate fast in HandleSpill() than that
    // this function is fast.

    for(auto& it: fCuts){
      if(it.first.ID() == cut.ID()){
        it.second.sinks.push_back(s);
        return;
      }
    }

    fCuts.push_back(std::make_pair(cut, CutDetails()));
    fCuts.back().second.sinks.push_back(s);
  }

  // --------------------------------------------------------------------------
  void ExposureSource::HandleSpill(const caf::SRSpillProxy* spill)
  {
    assert(!fFinalized);

    for(auto& it: fCuts){
      const Cut& cut = it.first;
      double pot = cut.POT(spill);
      // Not implemented -> fall back to usual counting
      if(pot < 0) pot = spill->spillpot;


      // FD MC is 1 neutrino-per-event, so the livetime count will be all
      // screwed up there. Just go with 0, ie "unknown".
      double livetime = 0;
      if(spill->det != caf::kFARDET || !spill->ismc) // TODO || fSource == kCosmic)
        livetime = cut.Livetime(spill);

      // Not implemented -> no good way to count
      if(livetime < 0){
        if(spill->det == caf::kFARDET && !spill->ismc){
          // People might well expect FD data to have a livetime. Warn.
          static bool once = true;
          if(once){
            once = false;
            std::cout << "WARNING: no way to compute livetime for FD data spectrum. If you want a livetime you need to be applying one of the cuts from TimingCuts.h or similar. You probably should be anyway to remove bad data near the spill ends." << std::endl;
          }
        }
        // Otherwise the user probably doesn't care. Just don't sum it into the
        // livetime
        livetime = 0;
      }

      it.second.pot += pot;
      it.second.livetime += livetime;
    }
  }

  //----------------------------------------------------------------------
  void ExposureSource::Finalize()
  {
    if(fFinalized){
      std::cout << "ExposureSource::Finalize() called twice!" << std::endl;
      abort();
    }

    for(auto& it: fCuts){
      for(IExposureSink* s: it.second.sinks){
        s->HandlePOT(it.second.pot);
        s->HandleLivetime(it.second.livetime);
      }
    }

    fFinalized = true;
  }

  //----------------------------------------------------------------------
  void POTPlusPercent(std::ostream& os, double pot, double denom)
  {
    os << pot << " POT";
    if(denom != 0) os << " (" << int(100*pot/denom+.5)<< "%)";
  }

  //----------------------------------------------------------------------
  void ExposureSource::ReportExposures(double denom) const
  {
    std::cout << "Total exposure:" << std::endl;

    // Normally using a floating-point number as a key is a bad idea. But this
    // seems to work, and if it goes wrong all it breaks is this printout.
    std::map<double, int> livetimes;
    std::map<double, int> pots;
    for(auto it: fCuts){
      livetimes[it.second.livetime] += it.second.sinks.size();
      pots[it.second.pot] += it.second.sinks.size();
    }

    if(pots.size() == 1){
      std::cout << "  ";
      POTPlusPercent(std::cout, pots.begin()->first, denom);
      std::cout << std::endl;
    }
    else{
      for(auto it: pots){
        std::cout << "  ";
        POTPlusPercent(std::cout, pots.begin()->first, denom);
        std::cout << " for " << it.second << " spectra" << std::endl;
      }
    }

    if(livetimes.size() == 1){
      if(livetimes.begin()->first == 0){
        std::cout << "  Unknown livetime" << std::endl;
      }
      else{
        std::cout << "  " << livetimes.begin()->first << " seconds" << std::endl;
      }
    }
    else{
      for(auto it: livetimes){
        std::cout << "  " << it.first << " seconds for " << it.second << " spectra" << std::endl;
      }
    }
  }
}