EvtTimeFNALBeam.cxx

Go to the documentation of this file.

#include "EvtTimeFNALBeam.h"
#include "EvtTimeShiftFactory.h"
TIMESHIFTREG3(evgb,EvtTimeFNALBeam,evgb::EvtTimeFNALBeam)

#include <iostream>
#include <iomanip>
#include <algorithm>

//GENIE includes
#ifdef GENIE_PRE_R3
  #include "GENIE/Utils/StringUtils.h"
#else
  #include "GENIE/Framework/Utils/StringUtils.h"
#endif
#include "messagefacility/MessageLogger/MessageLogger.h"
#include "cetlib_except/exception.h"

const double ksigma2fwhm = 2.0*std::sqrt(2.0*std::log(2.0));

namespace evgb {

  EvtTimeFNALBeam::EvtTimeFNALBeam(const std::string& config)
    // default to a NuMI config
    : EvtTimeShiftI(config)
    , fTimeBetweenBuckets(1e9/53.103e6)
    , fBucketTimeSigma(0.750)
    , fNBucketsPerBatch(84)        // NOvA-era 81+3, MINOS-era 81+5
    , fNFilledBucketsPerBatch(81)  // 81 for both eras
    , fDisallowedBatchMask(6,0)    // don't disallow any
    , fGlobalOffset(0)
  {
    std::vector<double> bi(6,1.0); // 6 equal batches
    SetBatchIntensities(bi);
    Config(config);
  }

  EvtTimeFNALBeam::~EvtTimeFNALBeam() { ; }

  void EvtTimeFNALBeam::Config(const std::string& config)
  {
    // parse config string
    if ( config == "" ) return;
    // GENIEHelper does a PrintConfig() when it gets this object ...

    // not the most sophisticated of parsing ... but FHICL would be overkill

    std::vector<std::string> strs = GetConfigTokens(config);
    // strings have all been tokenized and made lowercase

    size_t nstrs = strs.size();
    for (size_t i=0; i<nstrs; ++i) {
      if ( strs[i] == "numi" ) {
        fTimeBetweenBuckets     = 1e9/53.103e6;
        fBucketTimeSigma        = 0.750;
        fNBucketsPerBatch       = 84;  // NOvA-era 81+3, MINOS-era 81+5
        fNFilledBucketsPerBatch = 81;  // 81 for both eras
        fDisallowedBatchMask    = std::vector<int>(6,0); // don't disallow any
        fGlobalOffset           = 0;
        std::vector<double> bi(6,1.0); // 6 equal batches
        SetBatchIntensities(bi);
      } else
      if ( strs[i] == "booster" ) {
        fTimeBetweenBuckets     = 1e9/53.103e6;
        fBucketTimeSigma        = 2.0;
        fNBucketsPerBatch       = 84;  //
        fNFilledBucketsPerBatch = 81;  //
        fDisallowedBatchMask    = std::vector<int>(1,0); // don't disallow any
        fGlobalOffset           = 0;
        std::vector<double> bi(1,1.0); // 1 batch
        SetBatchIntensities(bi);
      } else
      if ( strs[i].find("intensity") != std::string::npos ) {
        // a list of batch intensities ... sets # of batches
        // eat values up until we see the end, or a word
        std::vector<double> bi;
        // count how many of next tokens are numbers (crude)
        for (size_t jj=i+1; jj<nstrs; ++jj) {
          // very crude check of being a number
          // can be  fooled by strange text
          size_t pos = strs[jj].find_first_not_of("0123456789.-+eE");
          if ( pos != std::string::npos ) break;
          // looks like a numeric value
          double val = atof(strs[jj].c_str());
          if ( val < 0 ) {
            mf::LogError("EvtTime")
              << "EvtTimeFNALBeam 'intensity' value [" << (jj-i-1)
              << "]=" << val << " '" << strs[jj] << "' "
              << "can't be less than zero, setting to zero";
            val = 0;
          }
          bi.push_back(val);
        }
        // ate up some strings ... move loop index
        i += bi.size();
        if ( bi.empty() ) {
          mf::LogError("EvtTime")
            << "EvtTimeFNALBeam error 'intensity' option didn't seem to have values";
        } else {
          SetBatchIntensities(bi);
        }
      } else
      if ( strs[i] == "bdisallowed" ) {
        mf::LogError("EvtTime")
          << "EvtTimeFNALBeam sorry 'bdisallowed' option not yet implemented";
      } else {
        // all the rest take one numeric value
        if ( i+1 >= nstrs ) {
          mf::LogError("EvtTime")
            << "EvtTimeFNALBeam sorry too few values for '" << strs[i] << "'";
          continue;
        }
        const char* arg = strs[i+1].c_str();
        if      ( strs[i] == "sigma"     ) fBucketTimeSigma        = atof(arg);
        else if ( strs[i] == "fwhm"      ) fBucketTimeSigma        = atof(arg)/ksigma2fwhm;
        else if ( strs[i] == "dtbucket"  ) fTimeBetweenBuckets     = atof(arg);
        else if ( strs[i] == "nperbatch" ) fNBucketsPerBatch       = atoi(arg);
        else if ( strs[i] == "nfilled"   ) fNFilledBucketsPerBatch = atoi(arg);
        else if ( strs[i] == "global"    ) fGlobalOffset           = atof(arg);
        else if ( strs[i] == "seed"      ) { ; } // handled in base
        else {
          mf::LogError("EvtTime")
            << "unknown EvtTimeFNALBeam config key '" << strs[i] << "'";
          continue;
        }
        ++i; // used up an argument
      }
    }
    // consistency check
    if ( fNFilledBucketsPerBatch > fNBucketsPerBatch ) {
      mf::LogError("EvtTime")
        << "EvtTimeFNALBeam nfilled " << fNFilledBucketsPerBatch
        << " of " << fNBucketsPerBatch << " buckets per batch,\n"
        << "set nfilled to match buckets per batch";
      fNFilledBucketsPerBatch = fNBucketsPerBatch;
    }

  }

  double EvtTimeFNALBeam::TimeOffset()
  {
    // calculate in small to large

    // pick a time within a bucket
    double offset = fRndmGen->Gaus(0.0,fBucketTimeSigma);

    // pick a bucket within a batch
    // assume all ~ buckets constant in batch until we have another model
    offset +=  fTimeBetweenBuckets *
               (double)fRndmGen->Integer(fNFilledBucketsPerBatch);

    // pick a bucket
    bool   disallowed = true;
    size_t ibatch = 0;
    size_t nbatch = fCummulativeBatchPDF.size();
    double r = 2;
    while ( disallowed ) {
      r = fRndmGen->Uniform();
      for (ibatch=0; ibatch<nbatch; ++ibatch) {
        if ( r <= fCummulativeBatchPDF[ibatch] ) break;
      }
      disallowed = ( fDisallowedBatchMask[ibatch] != 0 );
    }
    offset += fTimeBetweenBuckets*(double)fNBucketsPerBatch*(double)ibatch;

    // finally the global offset
    return offset + fGlobalOffset;
  }

  double EvtTimeFNALBeam::TimeOffset(std::vector<double> bi)
  {
    CalculateCPDF(bi);
    return TimeOffset();
  }

  void EvtTimeFNALBeam::PrintConfig(bool /* verbose */)
  {

    std::ostringstream msg;
    msg  << "EvtTimeFNALBeam config: \n"
         << "  TimeBetweenBuckets:     " << fTimeBetweenBuckets << " ns\n"
         << "  BucketTimeSigma:        " << fBucketTimeSigma << " ns "
         << "[FWHM " << fBucketTimeSigma*ksigma2fwhm << "]\n"
         << "  NBucketsPerBatch:       " << fNBucketsPerBatch << "\n"
         << "  NFilledBucketsPerBatch: " << fNFilledBucketsPerBatch << "\n"
         << "  Relative Fractions:    ";
    double prev=0;
    for (size_t i=0; i < fCummulativeBatchPDF.size(); ++i) {
      double frac = fCummulativeBatchPDF[i] - prev;
      bool skip = fDisallowedBatchMask[i];
      msg << " ";
      if (skip) msg << "{{";
      msg << frac;
      if (skip) msg << "}}";
      prev = fCummulativeBatchPDF[i];
    }
    msg << "\n"
        << "  GlobalOffset:           " << fGlobalOffset << " ns\n";

    mf::LogInfo("EvtTime") << msg.str();
  }


  void EvtTimeFNALBeam::SetBatchIntensities(std::vector<double> bi)
  {
    CalculateCPDF(bi);
  }

  void EvtTimeFNALBeam::SetDisallowedBatchMask(std::vector<int> disallow)
  {
    size_t ndis = disallow.size();
    size_t nbi  = fCummulativeBatchPDF.size();
    fDisallowedBatchMask = disallow;
    // expand it so it's mirrors # of batch intensities
    // but allow all that haven't been set
    if ( nbi > ndis ) fDisallowedBatchMask.resize(nbi,0);
  }

  void EvtTimeFNALBeam::CalculateCPDF(std::vector<double> bi)
  {
    fCummulativeBatchPDF.clear();
    double sum = 0;
    size_t nbi = bi.size();
    for (size_t i=0; i < nbi; ++i) {
      sum += bi[i];
      fCummulativeBatchPDF.push_back(sum);
    }
    // normalize to unit probability
    for (size_t i=0; i < nbi; ++i) fCummulativeBatchPDF[i] /= sum;
    // make sure the mask vector keeps up (but never make it smaller)
    // allowing all new batches
    if ( nbi > fDisallowedBatchMask.size() )
      fDisallowedBatchMask.resize(nbi,0);

    /*
    for (size_t j=0; j<nbi; ++j) {
      std::cout << " CPDF[" << j << "] " << fCummulativeBatchPDF[j]
                << "  " << ((fDisallowedBatchMask[j])?"dis":"") << "allowed"
                << std::endl;
    }
    */

  }


} // namespace evgb