ReadoutSim_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
/// \brief   TODO
/// \author  Ryan Patterson - rbpatter@caltech.edu
////////////////////////////////////////////////////////////////////////

#include <algorithm>
#include "sys/stat.h"

// Framework includes
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Core/EDProducer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "art/Framework/Services/Optional/TFileService.h"
#include "art/Framework/Services/Optional/TFileDirectory.h"
#include "art/Framework/Services/Optional/RandomNumberGenerator.h"

#include "cetlib/search_path.h"
#include "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

#include "CLHEP/Random/RandFlat.h"
#include "CLHEP/Random/RandExponential.h"
#include "CLHEP/Random/RandPoissonQ.h"
#include "CLHEP/Random/RandGaussQ.h"

// NOvA includes
#include "CMap/service/CMapService.h"
#include "Geometry/Geometry.h"
#include "GeometryObjects/CellGeo.h"
#include "MCCheater/BackTracker.h"
#include "NovaDAQConventions/DAQConventions.h"
#include "RawData/RawTrigger.h"
#include "ReadoutSim/FPGAAlgorithms.h"
#include "ReadoutSim/NoiseMaker.h"
#include "ReadoutSim/ExcessNoiseMaker.h"
#include "Simulation/PhotonSignal.h"
#include "Simulation/Simulation.h"
#include "Utilities/func/RandHisto.h"
#include "Utilities/func/MathUtil.h"
#include "ReadoutSim/IPulseShaper.h"
#include "ReadoutSim/PulseShaper.h"
#include "ReadoutSim/LegacyPulseShaper.h"

// ROOT includes
#include "TFile.h"
#include "TH2.h"
#include "TH3.h"
#include "TMath.h"
#include "TSystem.h"
#include "TGraph.h"
#include "TString.h"

#include "TH1D.h"

#include <set>
#include <utility>
#include <iomanip>
#include <iostream>
#include <vector>

namespace rsim
{
  class ReadoutSim: public art::EDProducer
  {
  public:
    explicit ReadoutSim(const fhicl::ParameterSet& pset);
    virtual ~ReadoutSim();

    void produce(art::Event& evt);
    void beginRun(art::Run& run);

  protected:
    // internal methods
    double GetSmearedPE(int nPE);
    double PE2ADC(double smearedPE, double gain);

    //double PerfectASICCurve(double t);
    bool   CreateRawDigits(const std::vector<sim::PhotonSignal>& hitlist,
                           std::unique_ptr< std::vector<rawdata::RawDigit> > &rawdigitlist,
                           rawdata::RawTrigger rawtrig);

    // Draw a trace and save it as a TGraph using the TFileService 
    void   DrawASICOutput(double *ASIC_Output, int offset, int stride, 
                          int plane, int cell, std::string name="trace_");

    // configuration settings
    std::string fPhotonModuleLabel;    ///< photons made in the module with this label
    std::string fCopyNoTruthHitsLabel; ///< hits with no truth (i.e., overlay) from this label will be copied into the digit vector after simulation is done

    fhicl::ParameterSet fPSetNDOS, fPSetND, fPSetFD, fPSetTB;

    CommonParameters* fParams;

    double fScaleEmptyCellNoiseRate; ///< scale factor for empty cell noise
    bool   fGeneratePhysicsCellNoise;

    std::string fEmptyCellNoiseFile_FD; ///< File holding FD noise PDF
    std::string fEmptyCellNoiseFile_ND; ///< File holding ND noise PDF
    std::string fEmptyCellNoiseFile_NDOS; ///< File holding NDOS noise PDF
    std::string fEmptyCellNoiseFile_TB; ///< File holding TB noise PDF

    // internal utilities
    rsim::IFPGAAlgorithm* fFPGAAlgo;
    rsim::NoiseMaker* fNoiseMaker;    ///< NoiseMaker object for making correlated noise
    rsim::ExcessNoiseMaker* fExcessNoiseMaker;    ///< ExcessNoiseMaker object for modeling 
    rsim::IPulseShaper*      fPulseShaper;         ///< PulseShaper object for generating traces

    TH3D*     fEmptyCellNoisePDF;     ///< for fast generation of empty-cell noise hits
    double    fEmptyCellNoisePDFIntegral;

    std::vector< TH1* >      fGainPDF;
    std::vector< std::vector< float > > fGainMap;
    bool      fUseRandomGains;
    bool      fVaryNoiseByThreshold;
    bool      fUseLegacyShaper;
    int       fNanosliceVersion;
    int       fLookbackSamples;
    bool      fVaryBaseline;
    bool      fUseNewExcessNoise;
    double    fFEBFlashThreshold;
    std::vector< std::vector< double > > fBaselineMap;

    int       LoadEmptyCellNoisePDF(const char* fname);
    void      PopulateRawTrig(rawdata::RawTrigger &rawtrig, art::Event const &evt);
    TH1D     *fADC;
    TH1D     *fMaxFebADC;
    TH2D     *fPlaneCellPE;
    TH2D     *fPlaneCellADC;
    TH1D     *fNoiseADC;

    bool fSag;
    double fSagFrac;

    /// Vector of plane, cell pairs for which to write traces to the TFile 
    /// for debugging purposes. 
    std::vector<std::pair<int, int> > fTracePlaneCells;
  };

  //----------------------------------------------------------------------
  ReadoutSim::ReadoutSim(const fhicl::ParameterSet& pset) :
    fPhotonModuleLabel(pset.get<std::string>("PhotonModuleLabel")),
    fCopyNoTruthHitsLabel(pset.get<std::string>("CopyNoTruthHitsFromLabel")),
    fPSetNDOS(pset.get<fhicl::ParameterSet>("ndos")),
    fPSetND  (pset.get<fhicl::ParameterSet>("nd")),
    fPSetFD  (pset.get<fhicl::ParameterSet>("fd")),
    fPSetTB  (pset.get<fhicl::ParameterSet>("tb")),
    fParams(0),
    fFPGAAlgo(0),
    fNoiseMaker(0),
    fExcessNoiseMaker(0),
    fPulseShaper(0),
    fEmptyCellNoisePDF(0),
    fEmptyCellNoisePDFIntegral(0),
    fGainPDF(), fUseRandomGains(true),
    fVaryNoiseByThreshold(true), fVaryBaseline(false), fUseNewExcessNoise(false), 
    fFEBFlashThreshold(10000.0),
    fADC(0), fPlaneCellPE(0), fPlaneCellADC(0)
  {

    // get the random number seed, use a random default if not specified    
    // in the configuration file.  
    unsigned int seed = pset.get< unsigned int >("Seed", sim::GetRandomNumberSeed());
  
    createEngine( seed );

    produces< std::vector<rawdata::RawDigit>   >();
    produces< std::vector<rawdata::RawTrigger> >();

    std::vector<int> tracePlanes = pset.get<std::vector<int> >("TracePlanes");
    std::vector<int> traceCells  = pset.get<std::vector<int> >("TraceCells");
    if(traceCells.size() != tracePlanes.size())
      throw cet::exception("InvalidTracePlaneCell") <<  
      "FCL configurable TracePlanes and TraceCells must have " <<
      "same number of elements \n";

    for(size_t i = 0; i < traceCells.size(); ++i){
      int plane = tracePlanes[i];
      int cell  = traceCells[i];
      fTracePlaneCells.push_back(std::make_pair(plane, cell)); 
    }

  }

  //----------------------------------------------------------------------
  ReadoutSim::~ReadoutSim()
  {
    if (fParams)            delete fParams;
    if (fFPGAAlgo)          delete fFPGAAlgo;
    if (fNoiseMaker)        delete fNoiseMaker;
    if (fExcessNoiseMaker)  delete fExcessNoiseMaker;
    if (fEmptyCellNoisePDF) delete fEmptyCellNoisePDF;
    if (fPulseShaper)       delete fPulseShaper;
  }

  //----------------------------------------------------------------------

  void ReadoutSim::beginRun(art::Run& run)
  {
    // Figure out which parameter set we're using
    fhicl::ParameterSet pset;
    art::ServiceHandle<geo::Geometry> geom;
    switch(geom->DetId()){
    case novadaq::cnv::kNDOS:
      pset = fPSetNDOS;
      break;
    case novadaq::cnv::kNEARDET:
      pset = fPSetND;
      break;
    case novadaq::cnv::kFARDET:
      pset = fPSetFD;
      break;
    case novadaq::cnv::kTESTBEAM:
      pset = fPSetTB;
      break;
    default:
      assert(0 && "Unknown detector");
    }

    fNanosliceVersion     = pset.get<int>("FPGA_DCS_NanosliceVersion");
    fUseRandomGains       = pset.get<bool>("UseRandomGains");
    fVaryNoiseByThreshold = pset.get<bool>("VaryNoiseByThreshold");
    fVaryBaseline         = pset.get<bool>("VaryBaseline");
    fUseNewExcessNoise    = pset.get<bool>("UseNewExcessNoise");
    fFEBFlashThreshold    = pset.get<double>("FEBFlashThreshold");
    fLookbackSamples      = pset.get<int>("FPGA_DCS_LookbackSamples");
    fUseLegacyShaper      = pset.get<bool>("UseLegacyShaper");

    // Noise histograms
    fScaleEmptyCellNoiseRate = pset.get<double>("ScaleEmptyCellNoiseRate");
    fGeneratePhysicsCellNoise = pset.get<bool>("GeneratePhysicsCellNoise");

    // Create helpers
    if(!fParams)           fParams = new CommonParameters(pset);
    if(!fFPGAAlgo)         fFPGAAlgo = CreateFPGAAlgorithm(pset);
    if(!fNoiseMaker)       fNoiseMaker = new rsim::NoiseMaker(pset);
    if(!fExcessNoiseMaker) fExcessNoiseMaker = new rsim::ExcessNoiseMaker(pset);
    if(!fPulseShaper) 
      {
        if (fUseLegacyShaper) fPulseShaper = new rsim::LegacyPulseShaper(pset);
        else fPulseShaper = new rsim::PulseShaper(pset);
      }
    
    int runNumber = run.run();
    if(runNumber >= 1000000 || runNumber < 10000) {
      int forceGain = pset.get<int>("ForceGain");
      if (forceGain == 0) {
        throw cet::exception("NoGainSetting")
          << "Ideal MC run number detected, but ForceGain == 0. Set ForceGain to the value of gain you want to simulate.\n"
          << __FILE__ << ":" << __LINE__ << "\n";
      }
    }

    // constructor decides if initialized value is a path or an environment variable
    cet::search_path sp("FW_SEARCH_PATH");
    
    std::string noiseFile;
    sp.find_file(fParams->fEmptyCellNoiseFile, noiseFile);
    
    std::string gainFile;
    sp.find_file(fParams->fGainFile, gainFile);
    
    struct stat sb;
    
    if(noiseFile.empty() || stat(noiseFile.c_str(), &sb) != 0){
      // failed to resolve the file name
      throw cet::exception("NoCellNoiseFile")
        << "noise file " << noiseFile << " not found!\n"
        << __FILE__ << ":" << __LINE__ << "\n";
    }
    
    if(gainFile.empty() || stat(gainFile.c_str(), &sb) != 0){
      // failed to resolve the file name
      throw cet::exception("NoGainFile")
        << "gain file " << gainFile << " not found!\n"
        << __FILE__ << ":" << __LINE__ << "\n";
    }
    
    if(fEmptyCellNoisePDF == 0){
      if( LoadEmptyCellNoisePDF(noiseFile.c_str()) ){
        mf::LogError("ReadoutSimBeginJob") << "failed to load cell noise pdf";
      }
    }
    
    if (fGainPDF.size() == 0){
      if (gSystem->AccessPathName(gainFile.c_str())){
        throw cet::exception("FileNotOpen")
          << "Unable to open " << gainFile << "\n"
          << __FILE__ << ":" << __LINE__ << "\n";
      }
      
      TFile f(gainFile.c_str());
      for (int ipix = 0; ipix < 32; ++ipix)
        {
          TString gainName = "hGains_PIX";
          gainName += ipix;
          TH1D *htemp = (TH1D*)f.Get(gainName);
          if (!htemp){
            throw cet::exception("FatalRootError") << "Unable to get gain histogram from " << gainFile << "\n"
                                                   << __FILE__ << ":" << __LINE__ << "\n";
          }
          gainName = "GainPDF_PIX";
          gainName += ipix;
          fGainPDF.push_back( (TH1D*)htemp->Clone(gainName) );
          fGainPDF.back()->SetDirectory(0);
        }
      f.Close();
      
      //first time through, fill with nominal gains
      for (unsigned int plane = 0; plane < geom->NPlanes(); plane++){
        std::vector< float > plane_gains;
        for (unsigned int cell = 0; cell < geom->Plane(plane)->Ncells(); cell++){
          plane_gains.push_back(fParams->fGain);
        }
        fGainMap.push_back(plane_gains);
      }
    }
    
    if (fVaryBaseline) {
      //Random baseline                                                                                
      fBaselineMap.clear();
      double baselineMean  = pset.get<double>("BaselineMean");
      double baselineSigma = pset.get<double>("BaselineSigma");
      art::ServiceHandle<art::RandomNumberGenerator> rng;
      CLHEP::HepRandomEngine &engine = rng->getEngine();
      CLHEP::RandGaussQ       gauss(engine);
      for (unsigned int plane = 0; plane < geom->NPlanes(); plane++){
        std::vector< double > plane_baselines;
        for (unsigned int cell = 0; cell < geom->Plane(plane)->Ncells(); cell++){
          double baseline = gauss.fire(baselineMean, baselineSigma);
          plane_baselines.push_back(baseline);
        }
        fBaselineMap.push_back(plane_baselines);
      }
    }

    fFPGAAlgo->FetchThresholds();

    fSag = pset.get<bool>("Sag");
    fSagFrac = pset.get<double>("SagFrac");
    
    // Histograms
    if(fPlaneCellPE) return; // Don't make histograms twice

    art::ServiceHandle<art::TFileService> tfs;

    fADC         = tfs->make<TH1D>("ADC", ";ADC", 4095, 0, 4095);
    fMaxFebADC   = tfs->make<TH1D>("MaxFebADC", ";MaxFebADC", 2000, 0, 20000);
    fNoiseADC    = tfs->make<TH1D>("NoiseADC",";ADC3-ADC0", 150, 0, 300);
    fPlaneCellPE = tfs->make<TH2D>("planeCellPE", ";Plane;Cell",
                                   geom->NPlanes(), 0, geom->NPlanes(),
                                   geom->Plane(0)->Ncells(), 0, geom->Plane(0)->Ncells());
    fPlaneCellADC = tfs->make<TH2D>("planeCellADC", ";Plane;Cell",
                                    geom->NPlanes(), 0., geom->NPlanes(),
                                    geom->Plane(0)->Ncells(), 0, geom->Plane(0)->Ncells());
    
  }

  //----------------------------------------------------------------------
  void ReadoutSim::produce(art::Event& evt)
  {
    std::unique_ptr< std::vector<rawdata::RawDigit> > rawdigitlist( new std::vector<rawdata::RawDigit>);
    std::unique_ptr< std::vector<rawdata::RawTrigger> > rawtriglist( new std::vector<rawdata::RawTrigger>);

    art::ServiceHandle<geo::Geometry> geom;

    //needs to be here for CMap to be initialized to get pixel-by-pixel gains
    //random gains will change event-by-event, but that shouldn't be an issue
    if (fUseRandomGains)
      {
        // get the RandomNumberGenerator service
        art::ServiceHandle<art::RandomNumberGenerator> rng;
        art::ServiceHandle<cmap::CMap> cmap;
        
        CLHEP::HepRandomEngine &engine = rng->getEngine();
        RandHisto randHisto(engine);
        //double gainScale(0.0), dummy1(0.0), dummy2(0.0);
        for (unsigned int plane = 0; plane < geom->NPlanes(); plane++) {
          for (unsigned int cell = 0; cell < geom->Plane(plane)->Ncells(); cell++) {
            const daqchannelmap::lchan lc = cmap->Map()->encodeLChan(geom->DetId(), plane, cell);
            const daqchannelmap::dchan dc = cmap->Map()->encodeDChan(lc);
            const int pix = cmap->Map()->getPixel(dc);
            
            //randHisto.GetRandom( fGainPDF[pix], gainScale, dummy1, dummy2 );
            //for now, just use the mean of the distibution for each pixel position
            double gainScale = fGainPDF[pix]->GetMean(1);
            fGainMap[plane][cell] = gainScale*fParams->fGain;
          }
        }
      }
    
    // Load hit list (PhotonSignal objects)
    art::Handle< std::vector<sim::PhotonSignal> > hitlistHandle;
    evt.getByLabel(fPhotonModuleLabel, hitlistHandle);

    std::vector<sim::PhotonSignal> hitlist;
    for(unsigned int i = 0; i < hitlistHandle->size(); ++i){
      if((*hitlistHandle)[i].NPhoton() > 0)
        hitlist.push_back((*hitlistHandle)[i]);
    }

    // Need to ensure that hits on the same channel are adjacent in the vector.
    // The efficiency optimizations in CreateRawDigits assume this to be true
    // so it's required for correctness.
    sim::SortByPlaneAndCell(hitlist);
  
    // Make a "trigger" happen first.  (We'll need the trigger
    // time when setting hits' TDCs.)
    // Write the raw trigger into the event.
    rawdata::RawTrigger rawtrig;
    PopulateRawTrig(rawtrig, evt);
    rawtriglist->push_back(rawtrig);
    evt.put(std::move(rawtriglist));

    // Create raw digit list from hit list (i.e., do all the work)
    int failed = CreateRawDigits(hitlist, rawdigitlist, rawtrig);
    if (failed){
      mf::LogWarning("FailedDigitCreation") << "failed to create RawDigits, \n"
                                            << "putting empty RawDigit list into event";
      evt.put(std::move(rawdigitlist));
      return;
    }

    art::ServiceHandle<cmap::CMap> cmap;
    for(unsigned int i = 0; i < rawdigitlist->size(); ++i){
      TString debugOut = cmap->GetPlane(&rawdigitlist->at(i));
      debugOut += " ";
      debugOut += cmap->GetCell(&rawdigitlist->at(i));
      debugOut += "  ";
      for (int iADC = 0; iADC < (int)rawdigitlist->at(i).NADC(); ++ iADC) {
        debugOut += rawdigitlist->at(i).ADC(iADC);
        debugOut += " ";
      }
      LOG_DEBUG("ReadoutSim") << debugOut;

      if (fNanosliceVersion == 3) {
      fADC->Fill(rawdigitlist->at(i).ADC(3) - rawdigitlist->at(i).ADC(0));
      fPlaneCellADC->Fill(cmap->GetPlane(&rawdigitlist->at(i)), cmap->GetCell(&rawdigitlist->at(i)), rawdigitlist->at(i).ADC(3) - rawdigitlist->at(i).ADC(0)  ); 
      }
    }

    // Ensure all the digits are flagged as MC
    for(unsigned int i = 0; i < rawdigitlist->size(); ++i) 
      rawdigitlist->at(i).SetMC();

    // If we want to copy unsimulated hits from a preceding DAQ run, do that now before storing them.
    // (Primary use case: cosmic overlay hits in light level variations where we're resimulating the
    //  photons in the MC hits but don't want to drop the data cosmic hits.)
    if (!fCopyNoTruthHitsLabel.empty())
    {
      art::Handle< std::vector<rawdata::RawDigit>> digitcol;
      evt.getByLabel(fCopyNoTruthHitsLabel, digitcol);

      if (digitcol->empty())
        mf::LogWarning("ReadoutSim") << "Found no digits in fCopyNoTruthHitsLabel label '" << fCopyNoTruthHitsLabel << "' to copy!";

      for (const auto & srcDigit : *digitcol)
      {
        // only want digits NOT from simulation
        if (!srcDigit.IsMC())
          rawdigitlist->push_back(srcDigit);
      }
    }

    evt.put(std::move(rawdigitlist));

    try{
      art::ServiceHandle<cheat::BackTracker> bt;
      // This is the first point in MC production that we have all the products
      // that BackTracker wants to build maps of. Ensure they're built
      // correctly for all downstream BackTracker users.
      bt->Rebuild(evt);
    }
    catch(...){
      // There's no BackTracker service in this job. That's fine, it just
      // means we don't have to do anything here.
    }

    return;
  }


  //---------------------------------------------
  //-------------- CreateRawDigits --------------
  //---------------------------------------------
  bool ReadoutSim::CreateRawDigits
  (const std::vector<sim::PhotonSignal>& hitlist,
   std::unique_ptr< std::vector<rawdata::RawDigit> > &rawdigitlist,
   rawdata::RawTrigger rawtrig)
  {
    // This will return (via the argument) a vector of RawDigits
    // that holds the simulated readout for the hits in hitlist.
    // There will be one RawDigit object per active channel in the
    // returned vector.
    //
    // The steps:
    //   (0) Init (e.g., build template ASIC output curve)
    //   (1) For each cell, lay down scaled/jittered ASIC output curves
    //       on an array of clockticks (all clockticks for simplicity)
    //   (2) Pick out the ones the MUXer will actually pass through
    //       (depends on detector and channel)
    //   (3) Pass these digitizations to the appropriate digital signal
    //       processor.  For readibility, also pass a time array even
    //       though it can always be back-computed from the array indices.
    //
    // The DSP steps are discussed in the DSP code.
    //

    // get the RandomNumberGenerator service
    art::ServiceHandle<art::RandomNumberGenerator> rng;
    CLHEP::HepRandomEngine &engine = rng->getEngine();
    CLHEP::RandPoissonQ     poisson(engine);
    CLHEP::RandFlat         flat(engine);
    CLHEP::RandGaussQ       gauss(engine);
    CLHEP::RandExponential  expo(engine);
    
    RandHisto randHisto(engine);

    art::ServiceHandle<geo::Geometry> geom;
    art::ServiceHandle<cmap::CMap>    cmap;
    
    // Things we'll need:
    //const int nHits = hitlist.size();

    // Not all of the clockticks result in a digitization, thanks to the MUXing.
    // Which ones get seen by the FPGA depends on the channel and the detector.
    // Thus, we must pass to the FPGA algorithm two things: (1) the first clocktick
    // to consider, and (2) the number of clockticks that elapse between relevant
    // digitizations.  (The algorithm may also have other parameters that can be
    // passed in to create algo variants without duplicating code.)

    // Loop over hits.
    // In order to handle all the hits on each channel together, simply
    // loop twice.  The inner loop is to look for more hits on the
    // same channel. Because the list had been sorted, this code
    // runs with O(N) complexity.

    // If APD sag is enabled, create a sagged baseline trace for every FEB that
    // has physics in it. Indexing here is first by plane and then by "FEB
    // number", which for these purposes is just cell/32.
    //key is (plane, cell)
    std::map< std::pair<int, int>, std::list< std::pair<double, double> > > ADCPulseMap;
    //key is (plane, feb)
    std::map< std::pair<int, int>, std::list< std::pair<double, double> > > SaggedPulseMap;
    
    std::vector<bool> channelHasHits[geom->NPlanes()];
    for (unsigned int plane=0;plane<geom->NPlanes();plane++) {
      for (unsigned int cell=0;cell<geom->Plane(plane)->Ncells();cell++) {
        channelHasHits[plane].push_back(false);
      }
    }
    
    for(const sim::PhotonSignal& phot: hitlist){
      const int plane = phot.Plane();
      const int cell = phot.Cell();
      
      std::pair<int, int> hitCoord(plane, cell);
      std::pair<int, int> febCoord(plane, cell/32);
      
      int nPE = phot.NPhoton();
      double smearedPE = GetSmearedPE(nPE);
      double ADC = PE2ADC(smearedPE, fGainMap[plane][cell]);
      double time = phot.TimeMean();

      channelHasHits[plane][cell] = true;

      if (fSag) {
        //sagged pulses originating from the same channel as the unsagged version will cancel to give the original ADC after clustering
        SaggedPulseMap[febCoord].push_back( std::make_pair(time, -fSagFrac*ADC) );
        ADCPulseMap[hitCoord].push_back( std::make_pair(time, (1.0 + fSagFrac)*ADC) );
      }
      else ADCPulseMap[hitCoord].push_back( std::make_pair(time, ADC) );
    }
    
    for (auto&& it: ADCPulseMap) it.second.sort();    
    for (auto&& it: SaggedPulseMap) it.second.sort();
    
    //cluster the sagged hits within 15 ns to determine if enough charge was collected to trigger an FEB flash
    std::vector<bool> febHasFlashHit[geom->NPlanes()];
    for (unsigned int plane=0;plane<geom->NPlanes();plane++) {
      for (unsigned int feb=0;feb<geom->Plane(plane)->Ncells()/32;feb++) {
        std::list< std::pair<double,double> > SaggedPulses = SaggedPulseMap[ std::make_pair(plane, feb) ];
        double meanTime(0.0);
        double ADC(0.0);
        double maxADC(0.0);
        for (auto it = SaggedPulses.begin(); it != SaggedPulses.end(); ++it) {
          double currTime = it->first;
          double currADC  = fabs(it->second);
          if (ADC == 0) {
            meanTime = currTime;
            ADC      = currADC;
          }
          else if ( std::fabs(meanTime - currTime) < 15.0 ) {
            meanTime = meanTime*ADC + currTime*currADC;
            ADC += currADC;
            meanTime /= ADC;
          }
          else {
            if (ADC > maxADC) maxADC = ADC;
            meanTime = currTime;
            ADC      = currADC;
          }
        }
        if (ADC > maxADC) maxADC = ADC;
        maxADC /= fSagFrac;
        fMaxFebADC->Fill(maxADC);
        if (maxADC > fFEBFlashThreshold) febHasFlashHit[plane].push_back(true);
        else febHasFlashHit[plane].push_back(false);
      }
    }
    
    // For keeping track of which hits were treated already (by
    // having had a previous hit on the same channel)...
    //bool alreadyDone[nHits];
    //for (int iH=0; iH<nHits; iH++) alreadyDone[iH] = false;

    // Keeping it simple here.  Channels without physics activity need
    // noise manually applied.  That's most of them.  Here's how I'll
    // keep track.

    // Each FEB, in each plane, has a readout phase set at random. Precalculate
    // them all here.
    //
    // One theoretical problem with this method is that each new spill gets a
    // new phasing pattern, but I'm hard pressed to think why you'd be doing a
    // study that could notice this.
    std::vector<std::vector<int>> phases;
    phases.resize(geom->NPlanes());
    for(unsigned int i = 0; i < geom->NPlanes(); ++i){
      for(unsigned int j = 0; j < geom->Plane(i)->Ncells()/32; ++j){
        phases[i].push_back(flat.fireInt(fParams->fNumClockticksPerDigitization));
      }
    }

    //will need this for making noise on empty cells
    if (fEmptyCellNoisePDF == 0) return 1; // return if error

    //loop over all planes and cells
    for (unsigned int plane=0;plane<geom->NPlanes();plane++) {
      for (unsigned int cell=0;cell<geom->Plane(plane)->Ncells();cell++) {    
        const daqchannelmap::lchan lchan = cmap->Map()->encodeLChan(geom->DetId(), plane, cell);
        const daqchannelmap::dchan dchan = cmap->Map()->encodeDChan(lchan);
        const int pix = cmap->Map()->getPixel(dchan);
        // Within the FEB the phase of the digitization proceeds round-robin
        // around the pixels: 0,8,16,24 together then 1,9,17,25, etc...
        const int phase = phases[plane][cell/32];
        int firstDigitizationClocktick = ((pix+phase)%fParams->fNumClockticksPerDigitization);

        //Set plane and cell so the proper threshold is used
        fFPGAAlgo->SetPlaneAndCell(plane,cell);
        //Get threshold for noise in this channel
        double thresh = fFPGAAlgo->ThresholdADC();
        //Get the baseline
        double baseline = fParams->fASICBaselineInADCCounts;
        if (fVaryBaseline) baseline = fBaselineMap[plane][cell];
        //Get the noise scaling
        double noiseFactor = 1.0;
        if (fVaryNoiseByThreshold) noiseFactor = fFPGAAlgo->ThresholdADC()/fNoiseMaker->GetNominalThreshold();

        //normally only simulate a trace if a channel has a real hit on it
        //if any channel on the feb has a hit larger than the flasher threshold, simulate a trace for all channels on that feb
        if ( channelHasHits[plane][cell] || febHasFlashHit[plane][cell/32] ) {
          std::list< std::pair<double,double> > ADCPulses = ADCPulseMap[ std::make_pair(plane, cell) ];
          std::list< std::pair<double,double> > SaggedPulses = SaggedPulseMap[ std::make_pair(plane, cell/32) ];

          fPulseShaper->SetPhase(firstDigitizationClocktick);
          fPulseShaper->SetBaseline(baseline);
          //construct the trace
          // This holds the real-valued "voltages" coming out of the ASIC
          // (post-shaping).  These will be digitized (i.e., integer-ized)
          // after all the hits on the channel have contributed.
          std::vector<double> ASIC_Output = fPulseShaper->CreateTrace( ADCPulses, SaggedPulses );
          // Both "current" and "voltage" noise are being handled by the NoiseMaker
          // object.  See doc-4401 for background.  This eliminates the need for
          // explicit sprinkling of dark noise.
          //
          // An annoyance: the noise formalism expects that I will be needing noise
          // only every digitization cycle, not every clocktick.  The correlations
          // will be wrong if I draw a noise value every clocktick.  So, I'll only
          // put noise on the to-be-used clockticks, which is kind of weird and
          // will make the DumpASICOutput routine less useful.  I can maybe clean
          // this up later.
          //
          // In the same loop, pass the ASIC output through the ADC.  This amounts to:
          //   (a) lifting the output by a baseline
          //   (b) flooring the real-valued numbers to make integers
          //   (c) truncating the range to [0, (1<<12)-1]
          //
          for (int iClocktick=firstDigitizationClocktick; 
               iClocktick<fParams->fNumClockticksInSpill; 
               iClocktick += fParams->fNumClockticksPerDigitization) {
            if (fGeneratePhysicsCellNoise) ASIC_Output[iClocktick] += noiseFactor*fNoiseMaker->GetNoise(&gauss)/fParams->fADCMaxPE * ( (1<<12) - 1 );
            ASIC_Output[iClocktick] = int(ASIC_Output[iClocktick]);
            if      ( ASIC_Output[iClocktick] < 0         ) ASIC_Output[iClocktick] = 0        ;
            else if ( ASIC_Output[iClocktick] > (1<<12)-1 ) ASIC_Output[iClocktick] = (1<<12)-1;
          } //end iClocktick

          // Draw traces if requested in fcl                                                                
          if(!fTracePlaneCells.empty()) {
            for(auto tracePlaneCell:fTracePlaneCells) {
              if(plane != (unsigned int)tracePlaneCell.first) continue;
              // Draw everything in the same FEB, important in the sag case                                 
              if(cell/32 != (unsigned int)tracePlaneCell.second/32) continue;
              
              DrawASICOutput(&ASIC_Output[0], firstDigitizationClocktick,
                             fParams->fNumClockticksPerDigitization,
                             plane, cell, "trace_");
              
            } // for tracePlaneCell                                                          
          }

          // Process the now-complete signal through the FPGA.
          //
          // The list of FPGA algorithms can grow.  Add new ones by following the
          // examples already present.
          //
          // Different detectors may have different signal processing
          // algorithms. At the moment, though, all the detectors are controlled by
          // a single string parameter
          
          // Call the requested algorithm to turn the ASIC output into a list of
          // RawDigits to be added to the outgoing list
          std::vector<rawdata::RawDigit> local_digitlist =
            fFPGAAlgo->ASICToDigits(&ASIC_Output[0], firstDigitizationClocktick);
          LOG_DEBUG("ReadoutSim") << "    == size of digilist == "
                                  << local_digitlist.size();

          // Store any rawdigits (if there are actually triggered hits on the
          // channel)
          for (unsigned int iD=0; iD<local_digitlist.size(); iD++) {
            rawdata::RawDigit digit = local_digitlist[iD];

            if ( ! digit.NADC() ) {
              throw cet::exception("NoADC")
                << "ReadoutSim::CreateRawDigits() -- FPGA algorithm returned "
                << "a digit with no ADC values.  This should never happen. Fatal!\n"
                << __FILE__ << ":" << __LINE__ << "\n";
              return 1;
            }
            // Finalize the digit
            digit.SetChannel(lchan);
            digit.SetDaqChannel(dchan);

            // Store the digit!
            LOG_DEBUG("ReadoutSim") << " -- digit -- " << plane << " " << cell
                                    << " " << digit.TDC() << " " << digit.ADC(0);
            
            rawdigitlist->push_back(digit);
          } // end iD
        
        } //end for channels with hits
        else {
          if( fScaleEmptyCellNoiseRate == 0) continue;
          double rate = fScaleEmptyCellNoiseRate*fEmptyCellNoisePDFIntegral; // prob. per cell per ns
          double mean = rate*fParams->fClocktick*fParams->fNumClockticksPerDigitization; //prob. per cell per digitization
          double mean_arrival = 0;
          if (mean != 0)
            mean_arrival = 1.0/mean;
          
          //start at the offset first digitization clocktick
          for (Long64_t iTick = firstDigitizationClocktick;
               iTick < fParams->fNumClockticksInSpill;
               iTick += fParams->fNumClockticksPerDigitization) {          
            //std::floor( expo.fire(mean_arrival) ) is the number of digitization periods until a noise hit
            //It's rounded down since any hit within a period should be assigned to the first clocktick of the period
            iTick += std::floor( expo.fire(mean_arrival) )*fParams->fNumClockticksPerDigitization;
            if (iTick < fParams->fNumClockticksInSpill) {
              // random adc values
              double adc0, adc1, adc2;
              randHisto.GetRandom( fEmptyCellNoisePDF, adc0, adc1, adc2 );
              // make the noise RawDigit
              rawdata::RawDigit digit;
              digit.SetVersion(fNanosliceVersion);
            
              //need to store peak time, not baseline time
              digit.SetTDC((iTick+fLookbackSamples*fParams->fNumClockticksPerDigitization)<<2); // shift by two bits to move 16 MHz-->64 MHz
              //digit.SetTDC(iTick<<2); // shift by two bits to move 16 MHz-->64 MHz
              // GetRandom will have given us a value between binLowEdge and
              // binHighEdge. All that bin's contents represent binLowEdge, so
              // rounding down is correct.
              if (fNanosliceVersion == 0) {
                if (fParams->fUseNewEmptyCellNoise) digit.SetADC(0,uint16_t(adc2));
                else                                digit.SetADC(0,uint16_t(adc0));
              }
              else {
                digit.SetADC(0,uint16_t(baseline));
                digit.SetADC(1,uint16_t(adc0 + baseline));
                digit.SetADC(2,uint16_t(adc1 + baseline));
                digit.SetADC(3,uint16_t(adc2 + baseline));
              }
              digit.SetChannel(lchan);
              digit.SetDaqChannel(dchan);
              
              // Store the noise digit!
              // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
              if (digit.ADC() >= thresh) {
                fNoiseADC->Fill(digit.ADC());
                rawdigitlist->push_back(digit);
              }
            }
          }
        } //end for channels without hits
      } //end for cell
    } //end for plane
    return 0;
  }
  
  //---------------------------------------------
  //-------------- PE2ADC --------------
  //---------------------------------------------
  double ReadoutSim::PE2ADC(double smearedPE, double gain)
  {
    //Converts photoelectrons to ADC

    return smearedPE / fParams->fADCMaxPE * ( (1<<12) - 1 ) * (gain/fParams->fGain);
  }
  
  //------------------------------------------------
  //-------------- GetSmearedPE------ --------------
  //------------------------------------------------
  double ReadoutSim::GetSmearedPE(int nPE) 
  {
    // Smear PEs to account for excess noise such that the mean PE remains the same
    // and the spread is sqrt(fAPDExcessNoiseFactor*nPE)
    
    // For smaller numbers of PE, use the ExcessNoiseMaker
    // For larger numbers, approximate with the log-normal distribution

    // get the RandomNumberGenerator service
    art::ServiceHandle<art::RandomNumberGenerator> rng;
    CLHEP::HepRandomEngine &engine = rng->getEngine();
    CLHEP::RandFlat         flat(engine);
    CLHEP::RandGaussQ       gauss(engine);

    if (fUseNewExcessNoise && nPE < 250) return fExcessNoiseMaker->DrawSmearedPE(nPE, &flat);
    else {
      if ( nPE == 0 ) {
        throw cet::exception("NPEEq0")
          << "ReadoutSim::GetASICScaleSpread() -- nPE==0.  This should never happen.\n"
          << __FILE__ << ":" << __LINE__ << "\n";
        return 1;
      }
      
      // This is the variance we need, to be combined with the preexisting poisson variance
      const double var = (fParams->fAPDExcessNoiseFactor-1)/nPE;
      // These paramters define a log-normal distribution with mean 1 and
      // variance given by 'var'
      const double mu = -log(1+var)/2;
      const double sigma = sqrt(log(1+var));
      
      // Generate a log-normal number with parameters mu and sigma
      return nPE*exp(mu+sigma*gauss.fire(0, 1));
    }
  }

  //---------------------------------------------
  //----------- LoadEmptyCellNoisePDF -----------
  //---------------------------------------------
  int ReadoutSim::LoadEmptyCellNoisePDF(const char* fname) {

    // Load PDF for generating noise on physics-less cell.
    // I'd love for these files to stay as root files.  Please
    // don't change that unless absolutely necessary.

    if (gSystem->AccessPathName(fname)) {
      cet::exception("FileNotOpen")  << "ReadoutSim: Unable to open " << fname
                                     << "   Cannot get noise PDFs!";
      return 1;
    }

    TFile f(fname);

    TH3D *htemp = (TH3D*)f.Get("noise/h3_adc0adc1adc2");
    if (!htemp) {
      std::cerr << "ReadoutSim: Histogram h3_adc0adc1adc2 not found in " << fname << std::endl;
      std::cerr << "   Cannot get fine timing!  " << std::endl;
      return 1;
    }
    fEmptyCellNoisePDF = (TH3D*)htemp->Clone();
    fEmptyCellNoisePDF->SetDirectory(0);

    f.Close();

    // Apply threshold to noise histogram so that it is consistent with the
    // parameters we're currently simulating.
    //const double thresh = fFPGAAlgo->ThresholdADC();
    //for(int i = 0; i < fEmptyCellNoisePDF->GetNbinsX()+2; ++i){
    //for(int j = 0; j < fEmptyCellNoisePDF->GetNbinsY()+2; ++j){
    //for(int k = 0; k < fEmptyCellNoisePDF->GetNbinsZ()+2; ++k){
    //  if(fEmptyCellNoisePDF->GetXaxis()->GetBinCenter(i) < thresh)
    //    fEmptyCellNoisePDF->SetBinContent(i, j, k, 0);
    //}
    //}
    //}

    fEmptyCellNoisePDFIntegral = fEmptyCellNoisePDF->Integral();

    return 0;
  }


  void ReadoutSim::PopulateRawTrig(rawdata::RawTrigger &rawtrig, art::Event const &evt) {
    // We may need to flesh out more of the trigger block in the
    // future (although not all of the items make sense in MC).

    // Say this is MC:
    rawtrig.fTriggerMask_MCBit = 1;

    // Units of this in data are 500ns
    rawtrig.fTriggerRange_TriggerLength = fParams->fClocktick*fParams->fNumClockticksInSpill/500.0;

    // Set the trigger time such that it increments by 0.1 seconds
    // every event (to mimic calibration triggers, though not beam
    // triggers).  Also, add one "day" per subrun and one "year" per
    // run, just to keep the times in different files different.
    // I'm not making these FHICL parameters because we will want
    // to do this properly (e.g., have each MC file representing a
    // specific timestamp drawn to match real data).
    //
    // The TDC ticks run at 64 MHz

    // using unsigned long long rather than uint64_t to ensure consistent behavior on OSX and LINUX
    static const unsigned long long ticksPerYear  = (unsigned long long)(64000000) * 31556926;
    static const unsigned long long ticksPerDay   = (unsigned long long)(64000000) * 86400;
    static const unsigned long long ticksPerEvent = (unsigned long long)(64000000) / 10; // 10 Hz trigger

    rawtrig.fTriggerTimingMarker_TimeStart =
      evt.run() * ticksPerYear +
      evt.subRun() * ticksPerDay +
      evt.event() * ticksPerEvent;

  }

  //---------------------------------------------
  //----------- DrawASICOutput        -----------
  //---------------------------------------------
  void   ReadoutSim::DrawASICOutput(double *ASIC_Output, int offset, int stride, 
                        int plane, int cell, std::string name){

    
    // Make this function sort of noisy so that it doesn't get accidentally 
    // get called in production or when folks don't expect it is being called. 
    std::cout << "Saving a trace for plane: " << plane 
              << " cell: " << cell << std::endl;
    
    std::vector<double> traceTimes;
    std::vector<double> trace;
    
    for(int iClocktick = offset; 
        iClocktick < fParams->fNumClockticksInSpill; 
        iClocktick += stride)
      {
        
        traceTimes.push_back(iClocktick * fParams->fClocktick);
        trace.push_back(ASIC_Output[iClocktick]);
      }
    
    
    art::ServiceHandle<art::TFileService> tfs;
    TGraph* traceGraph = tfs->make<TGraph>(traceTimes.size(), 
                                           &traceTimes[0], &trace[0]);     
    std::stringstream traceName;
    traceName<<name<<"_plane" <<plane<<"_feb"<<cell/32<< "_cell"<<cell; 
    traceGraph->SetName(traceName.str().c_str());
    traceGraph->SetTitle(";Time [ns]; ADC");
    traceGraph->Write();
    delete traceGraph;
  }
  
  
  
  DEFINE_ART_MODULE(ReadoutSim)

}//namespace