MdOptoElectronicSimulator.cc

Go to the documentation of this file.

#include "MdOptoElectronicSimulator.h"

#include <evt/Event.h>
//
#include <sevt/SEvent.h>
#include <sevt/Station.h>
#include <sevt/StationTriggerData.h>
#include <sevt/StationSimData.h>
#include <sdet/SDetector.h>
//
#include <mdet/MDetector.h>
#include <mdet/Counter.h>
#include <mdet/Module.h>
#include <mdet/Scintillator.h>
#include <mdet/Fiber.h>
#include <mdet/PMT.h>
#include <mdet/FrontEnd.h>
#include <mdet/Pixel.h>
#include <mdet/SiPMArray.h>
#include <mdet/FrontEndSiPM.h>
#include <mdet/SiPM.h>
#include <mdet/ChannelSiPM.h>
//
#include <mevt/MEvent.h>
#include <mevt/Counter.h>
#include <mevt/Module.h>
#include <mevt/Scintillator.h>
#include <mevt/ScintillatorSimData.h>
//
#include <utl/Particle.h>
#include <utl/PhysicalConstants.h>
#include <utl/Trace.h>
#include <utl/RandomEngine.h>
#include <utl/FFTDataContainer.h>
#include <utl/AugerUnits.h>
#include <utl/MathConstants.h>
#include <utl/ConsecutiveEnumFactory.h>
#include <utl/Branch.h>
#include <utl/ConfigUtils.h>
#include <utl/TimeStamp.h>
#include <utl/TimeInterval.h>
#include <utl/TabularStream.h>
#include <utl/Segment.h>
#include <utl/Plane.h>
#include <utl/Point.h>
#include <utl/GeometryUtilities.h>
//
#include <fwk/CentralConfig.h>
#include <fwk/RandomEngineRegistry.h>
//
#include <CLHEP/Random/RandFlat.h>
//
#include <cfloat>
#include <algorithm>
#include <vector>
#include <iostream>
// "Standardaly" this is _the_ place for endl and so (see http://www.research.att.com/~bs/3rd_issues.html):
#include <iomanip>
#include <sstream>
#include <string>
//
#include <boost/algorithm/minmax_element.hpp>
#include <boost/tokenizer.hpp>
//
#include <TF1.h>
#include <TCanvas.h>
#include <TGraph.h>
#include <TGaxis.h>
#include <TVector.h>
#include <TAxis.h>
#include <TLegend.h>
#include <TH1.h>
#include <TH2.h>
#include <TMath.h>
#include <TROOT.h>
#include <TStyle.h>


namespace {

class Tag {
        public:
          template<class Component>
          Tag(const std::string& prefix, const Component& c) {
                const mdet::Module&  mod = GetModule(c);
                const mdet::Counter& ctr = mod.GetCounter();
                fStream << prefix      << '_'
                                << ctr.GetId() << '_'
                                << mod.GetId() << '_';
          if( prefix != "outcomeIntegrator" )
                                fStream << c.GetId()  << '_'; // at least we're always putting the run number.
          }
          template<class S>
          Tag& operator()(const S& s) {
                fStream << s;
                return *this;
          }
          const std::stringstream& sstr() const {
                return fStream;
          }
        private:
          const mdet::Module& GetModule(const mdet::Pixel& p) const {
                return GetModule(p.GetPMT());
          }
          const mdet::Module& GetModule(const mdet::Channel& c) const {
                return GetModule(c.GetFrontEnd());
          }
          const mdet::Module& GetModule(const mdet::SiPM& s) const {
                return GetModule(s.GetSiPMArray());
          }
          const mdet::Module& GetModule(const mdet::ChannelSiPM& c) const {
                return GetModule(c.GetFrontEnd());
          }

          template<class Component>
          const mdet::Module& GetModule(const Component& c) const {
                return c.GetModule();
          }
          std::stringstream fStream;
};

template<class T1, class T2>
struct CanCopy {
  static void Constraints(T1 a, T2 b) {
    T2 c = a;
    b = a;
  }
  CanCopy() {
    void(*p)(T1,T2) = Constraints;
    ((void)p); // unused on purpose, avoid warnings.
  }
};

template<class StringContainer>
void PrintPlots(const TCanvas& c, const StringContainer& sc, const Tag& t)
{
  // Constrain the container to contain strings.
  // We do this after Stroustrup's FAQ answer
  // "Why can't I define constraints for my template parameters?"
  // http://www.research.att.com/~bs/bs_faq2.html#constraints
  // Using Boost's library may be too much.
  typedef typename StringContainer::value_type S;
  CanCopy<S, std::string>(); // trigger the check.
  const std::string& n = t.sstr().str();
  for (typename StringContainer::const_iterator i = sc.begin(), e = sc.end(); i != e; ++i)
    c.Print((n + *i).c_str());
}

template<typename T>
struct OwnerVector : public std::vector<T*> {
  void operator()(const T * const t) const {
    delete t;
  }
  ~OwnerVector() {
    typedef typename OwnerVector<T>::const_iterator It;
    typedef const OwnerVector<T>& ConstRef;
    // for_each works with copies by default, but we don't want
    // to mess around making copies of this vector (specially when
    // destroying the copies would lead to a double delete!)
    std::for_each<It,ConstRef>(this->begin(), this->end(), *this);
  }
};


/*
 * Several wrapping helper functors.
 */

template<class Container>
struct TotalPulseFunctor {
  TotalPulseFunctor(const Container& pulses, double ux, double uy) :
    fPulses(pulses), fUnitX(ux), fUnitY(uy), fAbs(false) { }
  double operator()(double t) const {
    double totalSignal = 0;
    // We accumulate the signal of all the pulses...
    for (typename Container::const_iterator p = fPulses.begin(), e = fPulses.end(); p != e; ++p)
      totalSignal += p->TimeSpaceEvaluate(t);
    if (fAbs)
      totalSignal = std::abs(totalSignal);
    return totalSignal;
  }
  double operator()(double *x, double * /*p*/) const {
    // For plottin', with units.
    return (*this)(x[0] * fUnitX) / fUnitY;
  }
  const Container& fPulses;//Love the abstraction.
  double fUnitX;
  double fUnitY;
  bool fAbs;
};

template<class Functor>
struct ProxyFunctor {
  ProxyFunctor(const Functor& func, double ux, double uy) :
    fFunctor(func), fUnitX(ux), fUnitY(uy) {}
  double operator()(double *x, double * /*p*/) const {
    return fFunctor(x[0] * fUnitX) / fUnitY;
  }
  const Functor& fFunctor;
  double fUnitX;
  double fUnitY;
};

struct PulseFunctor {
  PulseFunctor(const mdet::Pixel::SPE& spe, double ux, double uy) :
    fPulse(&spe), fUnitX(ux), fUnitY(uy) {}
  //overload for sipm
  PulseFunctor(const mdet::SiPM::PE& pe, double ux, double uy) :
    fPhotoEquivalent(&pe), fUnitX(ux), fUnitY(uy) {}
  double operator()(double *x, double * /*p*/) const {
    return fPulse->TimeSpaceEvaluate(x[0] * fUnitX) / fUnitY;
  }
  const mdet::Pixel::SPE * fPulse;
  const mdet::SiPM::PE * fPhotoEquivalent;
  double fUnitX;
  double fUnitY;
};

struct ThresholdFunctor {
  ThresholdFunctor(double t, double uy) :
    fThreshold(t), fUnitY(uy) {}
  double operator()(double * x = 0, double *p = 0) const {
    // Since the variables are not used, a warning may be issued,
    // but in the other hand they serve as placeholders for the default value so cannot
    // be commented out: cast to void to "use them". Note the cast in functional style
    // cannot be used: void(x) is taken as a new declaration.
    ((void)x);
    ((void)p);
    //
    // No need for unit in x.
    return fThreshold / fUnitY;
  }
  double fThreshold;
  double fUnitY;
};

struct TGraphFunctor {
  TGraphFunctor(const TGraph& g) : fGraph(g) {}
  double operator()(double *x, double * /*p*/) const {
    return fGraph.Eval(x[0]);
  }
  const TGraph& fGraph;
};

struct DerivativeFunctor {
  DerivativeFunctor(const TF1& f, double s = 1.0) : fF1(f), fScale(s) { }
  double operator()(double *x, double * /*p*/) const {
    return fScale * fF1.Derivative(x[0]);
  }
  const TF1& fF1;
  double fScale;
};

template<class T>
struct TraceFunctor {
  TraceFunctor(const T& t, double m, double ux, double uy) :
    fTrace(t), fMinTime(m), fUnitX(ux), fUnitY(uy), fBinning(t.GetBinning()), fAbs(false) {}
  double operator()(double *x, double * /*p*/) const {
    return ((*this)(x[0] * fUnitX)) / fUnitY;
  }
  double operator()(double t) const {
    // Do range checking since when running the operator gets called with
    // out of range values, and so we may access non-owned memory.
    double timeOffset = t - fMinTime;
    if (timeOffset < 0)
      return 0;
    double maxOffset = fTrace.GetSize() * fBinning;
    if (timeOffset > maxOffset)
      return 0;
    // We're in range, do a linear interpolation:
    typedef typename T::SizeType Size;
    /*
     * From the C++ standard ISO/IEC 14882:
     * 4.9) Floating-integral conversions
     * "An rvalue of a floating point type can be converted to an
     * rvalue of an integer type. The conversion truncates;
     * that is, the fractional part is discarded."
     */
    Size leftBin       = static_cast<Size>(timeOffset / fBinning);
    Size rightBin      = leftBin + 1;
    double leftOffset  = fBinning * leftBin;
    double rightOffset = fBinning * rightBin;
    double leftValue   = fTrace[ leftBin ];
    double rightValue  = fTrace[ rightBin ];
    double res = leftValue + (rightValue-leftValue) * (timeOffset - leftOffset) / (rightOffset - leftOffset);
    if (fAbs)
      res = std::abs(res);
    return res;
  }
  const T& fTrace;
  double fMinTime;
  double fUnitX;
  double fUnitY;
  double fBinning;
  bool fAbs;
};

}

using namespace fwk;

namespace MdOptoElectronicSimulatorAG {

const char* const
MdOptoElectronicSimulator::kSimulationTypeTags[] = { "FromMEvent", "FromMEventSimulatedScint"};

const char* const
MdOptoElectronicSimulator::kIntegratorSimulationTypeTags[] = { "StepByStep", "Simplified"};

MdOptoElectronicSimulator::MdOptoElectronicSimulator() :
  fRunNumber(1)
{
}

MdOptoElectronicSimulator::~MdOptoElectronicSimulator()
{
}

VModule::ResultFlag
MdOptoElectronicSimulator::Init()
{
  // Configuration reading.
  utl::Branch config = fwk::CentralConfig::GetInstance()->GetTopBranch("MdOptoElectronicSimulator");
  fLog.Configure(config);
  // Override defaults:
  fUnits.SetFrequencyDefault(utl::megahertz, "MHz");
  fUnits.SetElectricPotentialDefault(utl::milli*utl::volt, "mV");
  fUnits.SetElectricCurrentDefault(utl::milli*utl::ampere, "mA");
  fUnits.SetLengthDefault(utl::cm, "cm");
  fUnits.SetTimeDefault(utl::ns, "ns");
  fUnits.Configure(config);
  // Load the different file extensions
  std::string exts; // Load a raw string from configuration, and the fill the actual field.
  // XXX Check because VManager has, in fact, a GetData(list<string>) or vector<string>, which
  // ends calling >> upon a istringstream. And the use of a tokenizer is toooo-java like.
  LoadConfig(config, "plotFileExtensions", exts, ".eps .root");
  typedef boost::char_separator<std::string::value_type>  Sep;
  typedef boost::tokenizer<Sep> Tokenizer;
  Sep sep(" "); // separate with spaces
  Tokenizer tokens(exts, sep);
  fPlotFileExtensions.insert(fPlotFileExtensions.begin(), tokens.begin(), tokens.end());
  // Load the allowed particle types.
  // Load the defaults in a local variable...
  std::vector<ParticleType> particleTypesDefault;
  /*
   * particleTypesDefault.push_back(utl::Particle::eElectron);
   * particleTypesDefault.push_back(utl::Particle::eMuon);
   */
  // particleTypesDefault.insert(utl::Particle::eTau);
  // XXX Use a local variable of a type present in VManager,
  // and then pass the loaded data into the actual field.
  std::vector<ParticleType> particleTypes;
  LoadConfig(config, "allowedParticleTypes", particleTypes, particleTypesDefault);
  // Fill the actual field...
  fAllowedParticleTypes.insert(particleTypes.begin(), particleTypes.end());
  if (!fAllowedParticleTypes.size())
    fLog("No particle types indicated: allowing'em all", 1);
  // Load the kind of simulation to be performed.
  std::string tag;
  LoadConfig(config, "forcedSDTrigger", fForcedSDTrigger, false);
  LoadConfig(config, "nRepetitions", fNRepetitions, 10);
  LoadConfig(config, "nBinsHistograms", fNBinsHistograms, 100);
  LoadConfig(config, "minSPE", fMinSPE, 1);
  LoadConfig(config, "maxSPE", fMaxSPE, 5);
  LoadConfig(config, "stepSPE", fStepSPE, 1);
  LoadConfig(config, "nDiscretization", fNDiscretization, 2048);
  LoadConfig(config, "numPlotPoints", fNumPlotPoints, 500);
  LoadConfig(config, "ignoreCrossTalk", fIgnoreCrossTalk, false);
  LoadConfig(config, "pulseFilename", fPulseFilename, "");
  LoadConfig(config, "nPulseSamples", fNPulseSamples, 500);
  LoadConfig(config, "pulseSampleWindow", fPulseSampleWindow, 200 * fUnits.GetTimeUnit());
  LoadConfig(config, "generateSPEPulseOutput", fGenerateSPEPulseOutput, false);
  LoadConfig(config, "generatePreFETotalPulseOutput", fGeneratePreFETotalPulseOutput, false);
  LoadConfig(config, "generatePostFETotalPulseOutput", fGeneratePostFETotalPulseOutput, false);
  LoadConfig(config, "integratorSimType", tag, kIntegratorSimulationTypeTags[1]);
  fIntegratorSimType = utl::ConsecutiveEnumFactory<IntegratorSimulationType, eSimplified, kIntegratorSimulationTypeTags>::Create(tag);
  LoadConfig(config, "includeBaseLineFluctuationIntegrator", fIncludeBaseLineFluctuationIntegrator, true);
  LoadConfig(config, "injectNoiseBinary", fInjectNoiseBinary, false);
  // Check if a filename was loaded.
  if (! fPulseFilename.empty())
    fPulseFile.open(fPulseFilename.c_str());
  // Plot variables:
  LoadConfig(config, "togglePlotchannelPulses",     fTogglePlotchannelPulses, false);
  LoadConfig(config, "toggleOutcome",               fToggleOutcome, false);
  LoadConfig(config, "toggleInputOutput",           fToggleInputOutput, false);
  LoadConfig(config, "toggleTransferAmpResponse",   fToggleTransferAmpResponse, false);
  LoadConfig(config, "toggleTransferPhaseResponse", fToggleTransferPhaseResponse, false);
  LoadConfig(config, "toggleDelay",                 fToggleDelay, false);
  LoadConfig(config, "toggleFftPriorAmp",           fToggleFftPriorAmp, false);
  LoadConfig(config, "togglePlotFftPostAmp",        fTogglePlotFftPostAmp, false);
  // Now configure the defaults according the simulation type:

  fPlotChannelPulses         = false;
  fPlotOutcome               = false;
  fPlotInputOutput           = false;
  fPlotTransferAmpResponse   = false;
  fPlotTransferPhaseResponse = false;
  fPlotDelay                 = false;
  fPlotFftPriorAmp           = false;
  fPlotFftPostAmp            = false;
  // Now apply the toggle:
  if (fTogglePlotchannelPulses)
    fPlotChannelPulses         = !fPlotChannelPulses;
  if (fToggleOutcome)
    fPlotOutcome               = !fPlotOutcome;
  if (fToggleInputOutput)
    fPlotInputOutput           = !fPlotInputOutput;
  if (fToggleTransferAmpResponse)
    fPlotTransferAmpResponse   = !fPlotTransferAmpResponse;
  if (fToggleTransferPhaseResponse)
    fPlotTransferPhaseResponse = !fPlotTransferPhaseResponse;
  if (fToggleDelay)
    fPlotDelay                 = !fPlotDelay;
  if (fToggleFftPriorAmp)
    fPlotFftPriorAmp           = !fPlotFftPriorAmp;
  if (fTogglePlotFftPostAmp)
    fPlotFftPostAmp            = !fPlotFftPostAmp;
  Style_t font=102;
  Float_t fsize=0.03;
  Option_t *axis="XYZ";
  Float_t margin = 0.15;
  fStyle = new TStyle;
  //Canvas style
  fStyle->SetCanvasBorderMode(0);
  fStyle->SetCanvasColor(0);
  //Axis style
  fStyle->SetTickLength( 0.01, "XYZ" );
  fStyle->SetLabelFont( font , axis);
  fStyle->SetLabelSize( fsize, axis);
  fStyle->SetLabelOffset( 0.005, axis);
  fStyle->SetTitleFont( font , axis);
  fStyle->SetTitleSize( fsize, axis);
  fStyle->SetTitleOffset( 1.0, axis);
  fStyle->SetTitleXOffset( 1.5 );
  fStyle->SetTitleYOffset( 2.5 );
  fStyle->SetTitleAlign( 22 );
  fStyle->SetTitleBorderSize( 0 );
  fStyle->SetTitleColor( 1, axis );
  fStyle->SetTitleFillColor( 0 );
  fStyle->SetTitleH( 0 );
  fStyle->SetTitleW( 0 );
  fStyle->SetTitleX( 0.5 );
  fStyle->SetTitleY( 1-margin/2 );
  fStyle->SetTitleTextColor( kBlue );
  //fStyle->SetTitlePS(const char* pstitle);
  //fStyle->SetTitleStyle(Style_t style = 1001);
  //fStyle->SetTitleXSize(Float_t size = 0.02);
  //fStyle->SetTitleYSize(Float_t size = 0.02);
  fStyle->SetPalette(1,0);
  fStyle->SetStatColor(0);
  //Pad style
  fStyle->SetPadBottomMargin(margin);
  fStyle->SetPadLeftMargin  (margin);
  fStyle->SetPadRightMargin (margin);
  fStyle->SetPadTopMargin   (margin);
  fStyle->SetPadBorderMode( 0 );
  fStyle->SetPadBorderSize( 0 );
  fStyle->SetPadColor( 0 );
  fStyle->SetPadGridX(false);
  fStyle->SetPadGridY(false);
  fStyle->SetPadTickX(true);
  fStyle->SetPadTickY(false);
  //Frame style
  fStyle->SetFrameBorderMode( 0 );
  fStyle->SetFrameBorderSize( 0 );
  fStyle->SetFrameFillColor( 0 );
  fStyle->SetFrameFillStyle( 1 );
  fStyle->SetFrameLineColor( 0 );
  fStyle->SetFrameLineStyle( 0 );
  fStyle->SetFrameLineWidth( 0 );
  return eSuccess;
}

VModule::ResultFlag
MdOptoElectronicSimulator::Run(evt::Event& event)
{
  VModule::ResultFlag r = eFailure; // default.

  fStyle->cd();
  r = OptoElectronics(event);
  ++fRunNumber;
  return r;
}

/*********************************************************************************
 *                                                                               *
 * This method generate the PMT/SiPM pulses for the arriving photons             *
 * simulated in G4StationSimulator (loaded in ScintillatorSimData::PhotonTimes)  *
 * and runs the electronics simulation.                                          *
 * Replaces the RunFromMEvent method which first creates the photons with and    *
 * the optic fiber curve parametrization in the G4109 code.                      *
 *                                                                               *
 ********************************************************************************/
VModule::ResultFlag
MdOptoElectronicSimulator::OptoElectronics(evt::Event& event)
{
  if (event.HasMEvent()) {
    mevt::MEvent& me = event.GetMEvent();
    double deltaLatch = 0;
    const utl::TimeStamp& mEventTime = event.GetSEvent().GetHeader().GetTime();
    fLog("The event time for SiPM/PMT simulation ", 2)(mEventTime);
    const mdet::MDetector& detector = det::Detector::GetInstance().GetMDetector();
    for (mevt::MEvent::CounterIterator ic = me.CountersBegin(), ec = me.CountersEnd(); ic != ec; ++ic) {
      const mdet::Counter& counter = detector.GetCounter(ic->GetId());
      int triggerFound = GetTriggerTimeFromSD(event, counter, ic, deltaLatch);

      //Check trigger condition in WCD
      if (triggerFound) {
      //std::cout << "SD trigger found for counter " << ic->GetId() << ". MD: StartTrace - EventTime = " << deltaLatch/utl::ns << " ns " << std::endl;
        ic->SetT1(deltaLatch);
      } else if( fForcedSDTrigger) {
        //std::cout << "XXX NO SD trigger found for counter " << ic->GetId() << " but used anyway" << std::endl;
        //ic->SetT1(deltaLatch);
      } else {
            //std::cout << "XXX NO SD trigger found for counter " << ic->GetId() << std::endl;
            ic->SetRejected("SD associated un-triggered");
            continue;//do not go into modules of a non-triggered WCD
      }

      if (triggerFound || fForcedSDTrigger) {
        for (mdet::Counter::ModuleConstIterator mmDet = counter.ModulesBegin(); mmDet != counter.ModulesEnd(); ++mmDet) {
          if (!ic->HasModule(mmDet->GetId()))
            ic->MakeModule(mmDet->GetId());
        }
      }

      for (mevt::Counter::ModuleIterator im = ic->ModulesBegin(), em = ic->ModulesEnd(); im != em; ++im) {
        //Make all modules in the MEvent structure
        im->SetCandidate();
        const mdet::Module& module = counter.GetModule(im->GetId());
        SignalsMap signalsMap;
        for (mevt::Module::ScintillatorIterator is = im->ScintillatorsBegin(), es = im->ScintillatorsEnd();
              is != es; ++is) {
          if (is->HasSimData()) {
            const mdet::Scintillator& scint = module.GetScintillator(is->GetId());
            mevt::ScintillatorSimData& simData = is->GetSimData();

            for (mevt::ScintillatorSimData::ConstPhotonTimeIterator photon = simData.PhotonTimesBegin(), lastPhoton = simData.PhotonTimesEnd();
                photon != lastPhoton; ++photon) {

               const double delay = simData.GetPhotonTime(*photon);
               if ( !module.IsSiPM() ) { 
                const mdet::Pixel& pixel = module.GetPixelFor(scint);

                mdet::Pixel::SPE spe = pixel.MakeSPEAt(delay);

                int finalPixId = (pixel.GetPMT().ComputePulseDestination(pixel)).GetId();
                if (fIgnoreCrossTalk) {
                 finalPixId = pixel.GetId();
                }

                signalsMap[ finalPixId ].fPulses.push_back(spe);
                simData.AddSPEPulse(spe.GetMu(), spe.GetSigma(), spe.GetAmplitude(), finalPixId);

               } else {                  
                const mdet::SiPM& sipm = module.GetSiPMFor(scint);

                mdet::SiPM::PE pe = sipm.MakePEAt(delay);

                signalsMap[ sipm.GetId() ].fSiPMPulses.push_back(pe);
                simData.AddPEPulse(pe.GetStartTime(), pe.GetAmplitude1(), pe.GetAmplitude2(), pe.GetAmplitude3(),
                                pe.GetRiseTime(), pe.GetFallTime1(), pe.GetFallTime2(), pe.GetFallTime3(), sipm.GetId());
               }
            }
          }
        }

        VModule::ResultFlag f = SimulateElectronics(*im, module, signalsMap, mEventTime);
        if (f != eSuccess)
          return f;
      }
    }
  } else {
    INFO("Event without MEvent: nothing to be done.");
  }
  return eSuccess;
}


void
MdOptoElectronicSimulator::ProcessPulses(const mdet::Channel& channel,
                                         const SignalInformation& si,
                                         utl::TraceB& trace,
                                         double& span,
                                         const utl::TimeInterval& traceStart/*,
                                         const utl::TimeStamp& eventTime*/)

{
  //fLog("Processing pulses for channel", 1)(channel.GetId());
  const mdet::FrontEnd& frontEnd = channel.GetFrontEnd();
  // TODO No Jitter by now: so a constant sample rate.
  const double sRate = frontEnd.GetMeanSampleRatePeriod();
  // Clear the trace:
  trace.Clear();
  //
  // Set the appropiate binning equal to the (mean) sample rate.
  // Even in the case of having jitter, this would be the value because
  // the digital trace has no information of time, we just supose that that's
  // the period.
  trace.SetBinning(sRate);
  fLog("Cleaning trace and setting a binning of", 2)(sRate / fUnits.GetTimeUnit())(fUnits.GetTimeName());
  const unsigned int nBinsTrace = frontEnd.GetPreT1BufferLength() + frontEnd.GetPostT1BufferLength();
  fLog("Trace size", 2)(nBinsTrace);
  
  if (si.fPulses.empty()) {
    fLog("No pulses and other sim data found for channel", 1)(channel.GetId());
    // No pulses...so not a single 1 in the boolean trace...
    for (unsigned int i = 0; i < nBinsTrace; ++i)
      trace.PushBack(false);
    span = 0; // No pulse, no span.
    return; // done.
  }
  // First we determine the maximum extent of signal.
  // Init with the first so as to have something to compare with...
  // Now we are considering two cases times span before and after FrontEnd (see later also).
  double minTimePreFE = std::numeric_limits<double>::max();
  double maxTimePreFE = -std::numeric_limits<double>::max();
  std::unique_ptr<utl::TabularStream> spanTable;
  const unsigned int pulseSpanLogLevel = 3;
  if (fLog.GetLevel() >= pulseSpanLogLevel) {
    Init(spanTable, 3);
    *spanTable << utl::hline
               << "# spe" << utl::endc // 1
               << "min (" << fUnits.GetTimeName() << ")" << utl::endc // 2
               << "max (" << fUnits.GetTimeName() << ")" << utl::endr;// 3
  }
  unsigned int nSPE = 0; // only used when logging to the table...
  for (PulseContainer::const_iterator p = si.fPulses.begin(), e = si.fPulses.end(); p != e; ++p) {
    minTimePreFE = std::min(minTimePreFE, p->LowerLimit());
    maxTimePreFE = std::max(maxTimePreFE, p->UpperLimit());
    if (fLog.GetLevel() >= (pulseSpanLogLevel+1)) { // 4 each one, 1 more needed
      *spanTable << utl::hline
                 << ++nSPE << utl::endc
                 << fLog << p->LowerLimit() / fUnits.GetTimeUnit() << utl::endc
                 << fLog << p->UpperLimit() / fUnits.GetTimeUnit() << utl::endr;
    }
  }
  
  // The span covers both the range of pulses as well that of the analogic trace.
  const double pulseTimeSpan = maxTimePreFE - minTimePreFE;
  if (fLog.GetLevel() >= pulseSpanLogLevel) {
    *spanTable << utl::hline
               << "all" << utl::endc
               << fLog << minTimePreFE / fUnits.GetTimeUnit()  << utl::endc
               << fLog << maxTimePreFE / fUnits.GetTimeUnit()  << utl::endr
               << utl::hline;
    fLog("SPE pulses and other signals span:", pulseSpanLogLevel);
    fLog(spanTable->Str(), pulseSpanLogLevel);
  }
  // The pulses start, so take then into account...
  TotalPulseFunctor<PulseContainer> totalPulsePreFrontEnd(si.fPulses, fUnits.GetTimeUnit(),
                                                          fUnits.GetElectricPotentialUnit());
  fLog("Applying FFT...", 3);
  utl::FFTDataContainer<utl::Trace, TimeTrace::ValueType, FrequencyTrace::ValueType> fft;
  // Nyquist zone: first.
  // Discretization of total pulse (so to perform FFT).
  double binning = pulseTimeSpan / fNDiscretization;
  TimeTrace& timeTrace = fft.GetTimeSeries();
  fLog("Time discretization binning ", 4)(binning / fUnits.GetTimeUnit())(fUnits.GetTimeName());
  timeTrace.SetBinning(binning);
  fLog("Discretizing total pulses. Number of samples:", 3, false);
  // Fill the time series...
  for (unsigned int i = 0; i < fNDiscretization || fLog(i, 3)('.'); ++i) {
    double time  = minTimePreFE + i * binning;
    double value = totalPulsePreFrontEnd(time);

    fft.GetTimeSeries().PushBack(value);
  }
  // Now apply the circuit transfer function:
  // calling GetFrequencySpectrum triggers the update...
  fLog("Applying transfer function in frequency spectrum.", 3);
  FrequencyTrace& freqTrace            = fft.GetFrequencySpectrum();
  const FrequencyTrace::SizeType nFreq = fft.GetConstFrequencySpectrum().GetSize();
  const double freqBinning             = fft.GetConstFrequencySpectrum().GetBinning();
  TVectorD freqArray(nFreq);
  TVectorD ampArray(nFreq);
  TVectorD phaArray(nFreq);
  TVectorD preFFT(nFreq);
  TVectorD postFFT(nFreq);
  fLog("Generating frequency response.", 3);
  const unsigned int tabFFTLogLevel = 4;
  std::unique_ptr<utl::TabularStream> tabFFT;
  if (fLog.GetLevel() >= tabFFTLogLevel) {
    Init(tabFFT, 4);
    *tabFFT <<  utl::hline
           << "Size"        << utl::endc //1
           << "Binning ("   << fUnits.GetFrequencyName() << ")" << utl::endc //2
           << "1st freq ("  << fUnits.GetFrequencyName() << ")" << utl::endc //3
           << "Last freq (" << fUnits.GetFrequencyName() << ")" << utl::endr //4
           << nFreq         << utl::endc // 1
           << fLog << (freqBinning                      / fUnits.GetFrequencyUnit()) << utl::endc // 2
           << fLog << (fft.GetFrequencyOfBin(0)         / fUnits.GetFrequencyUnit()) << utl::endc // 3
           << fLog << (fft.GetFrequencyOfBin(nFreq - 1) / fUnits.GetFrequencyUnit()) << utl::endr // 4
           << utl::hline;
    fLog(tabFFT->Str(), tabFFTLogLevel);
  }
  for (FrequencyTrace::SizeType freqBin = 0; freqBin < nFreq; ++freqBin) {
    double freq = fft.GetFrequencyOfBin(freqBin);
    double preAmp = std::abs(freqTrace[ freqBin ]) / fUnits.GetElectricCurrentUnit();
    // This has resistance units..
    std::complex<double> c = channel.ComputeTransfer(freq);
    // Apply the transfer: now we have voltage...
    freqTrace[ freqBin ] *= c;
    // For plottin....
    double amp = std::abs(c) / fUnits.GetElectricResistanceUnit();
    double pha = std::arg(c) / fUnits.GetAngleUnit();
    freqArray(freqBin) = freq / fUnits.GetFrequencyUnit();
    ampArray(freqBin)  = amp;
    phaArray(freqBin)  = pha;
    preFFT(freqBin)    = preAmp;
    postFFT(freqBin)   = std::abs(freqTrace[ freqBin ]) / fUnits.GetElectricPotentialUnit();
  }
  fLog("Obtaining the amplified pulses in time-space...", 3);
  // Refresh back the time series: now we have the pulse widenned.
  const TimeTrace& totalPulsePostFrontEndTrace = fft.GetConstTimeSeries();
  /*
   * Is it correct to keep minTime and maxTime calculated prior front end?
   * 1st, we're not doing that anymore:
   *
   * a) Actually it isn't strictly so well; because the transfer function
   * has a non negligible phase effect: so it causes a shift in the output pulse.
   * Now, if the boundaries taken for the original pulse are quite far, it seems that
   * there's no actual problem.
   * If the limits, fall short with respect to the imposed shift, then the output pulse
   * may slip back entering from the other side of the axis as a periodicity effect due
   * to the transformation.
   *
   * b) In fact, in the real-life device there's an overall ~ 13 ns delay between the input and
   * output which is not being taken into account by now: ie we have only partially
   * described the real transfer function of the device.
   *
   * c) With the shift we're are taking into account b).
   */
  // Total analog pulse after the front-end, it also suffers a time shift
  // that rigidly shifts the total pulse:
  double shift = channel.ComputeSignalShift();
  fLog("Signal shift time", 3)(shift / fUnits.GetTimeUnit())(fUnits.GetTimeName())('.');
  // so the limits get shifted: compute new values, so to keep the original.
  const double minTimePostFE = minTimePreFE + shift;
  const double maxTimePostFE = maxTimePreFE + shift;
  // Note that here minTime is taken to be the reference start time, not
  // the "left" time window limit.
  TraceFunctor<TimeTrace> totalPulsePostFrontEnd(totalPulsePostFrontEndTrace,
                                                 minTimePostFE,
                                                 fUnits.GetTimeUnit(),
                                                 fUnits.GetElectricPotentialUnit());
  // Discrimination
  fLog("Creating discriminator object.", 3);
  const mdet::Channel::Discriminator discriminator = channel.MakeDiscriminator(
                                                      totalPulsePostFrontEnd,
                                                      minTimePostFE,
                                                      maxTimePostFE);
  // Loggin' of the times when the input reaches the threshold:
  // ie the roots of the equation threhols == input.
  const unsigned int rootLogLevel = 4;
  if (fLog.GetLevel() >= rootLogLevel) {
    typedef mdet::Channel::Discriminator::AtThresholdConstIterator ItRoot;
    ItRoot b = discriminator.AtThresholdBegin();
    ItRoot e = discriminator.AtThresholdEnd();
    if (b == e)
      fLog("Incoming pulse never at threshold level.", rootLogLevel);
    fLog("Incoming pulse at threshold level at ", rootLogLevel, false);
    for (ItRoot i = b; i != e; ++i) {
      fLog(*i / fUnits.GetTimeUnit(), rootLogLevel, false)(fUnits.GetTimeName())(' ');
    }
    fLog(".", rootLogLevel);
  }
  // Sampling
  fLog("Creating sampler object.", 3);
  mdet::FrontEnd::Sampler sampler = frontEnd.MakeSampler();
  /*
   * Note that we restrict the sampling to the maximun trace size according to buffer lenghts:
   * - First there's a all-false sampling trace segment, until reaching the zone of pulses,
   * - then the region with pulses,
   * - and lastly the trace is filled up with falses, until the max size is reached.
   */
  std::vector<double> sampleTimes;
  fLog("Sampling pulses.", 3);

  const double startTime = traceStart.GetInterval() - sRate* frontEnd.GetPreT1BufferLength();
  const double stopTime  = traceStart.GetInterval() + sRate* frontEnd.GetPostT1BufferLength();

  /*
   * XXX We could avoid everything if we can detect that there's no signal
   * inside the window frame defined by the trace.
   */
  // This is correct: we compare now with the updated minTime, which is the one
  // taken for the output signal.
  if (startTime < minTimePostFE) {
    // See inner body of the following loop (then condition).
    fLog("Samples false in the range", 4)
        (startTime / fUnits.GetTimeUnit())(fUnits.GetTimeName())
        ("to")
        (minTimePostFE / fUnits.GetTimeUnit())(fUnits.GetTimeName())
        ("as no signal in it.");
  }
  // The trace is empty: size is zero, so we start inserting values 'till t
  // he size gets _equal_ to the desired number of bins in the trace.
  // So when it's equal: don't eventer the for, cause one more will be added.
  int bin = 0;
  for (double sampleTime = startTime; trace.GetSize() < nBinsTrace || fLog('\n', 4);
    sampleTime += sRate, bin++) {
    //
    bool s;
    if (sampleTime < minTimePostFE || maxTimePostFE < sampleTime) {
      // If we are outside the time range where we have signal, then
      // take false as the certain sampled value.
      // See prev & next if condition and messages.
      s = false;
    } else {
      double o = discriminator(sampleTime);
      s = sampler(o);
      // Only log the relevant part:
      fLog(s, 4, false)('@')(sampleTime / fUnits.GetTimeUnit())(fUnits.GetTimeName())("(")(bin)(")")('|');
    }
    trace.PushBack(s);
    sampleTimes.push_back(sampleTime);
  }
  if (maxTimePostFE < stopTime) {
    // See inner body of the previous loop (then condition).
    fLog("Samples false in the range ", 4)
        (maxTimePostFE / fUnits.GetTimeUnit())(fUnits.GetTimeName())
        ("to")
        (stopTime / fUnits.GetTimeUnit())(fUnits.GetTimeName())
        ("as no signal in it.");
  }
  span = discriminator.GetOverThresholdTimeSpan();
  // PLOTS
  // If there are no extensions, then no plot will be performed at last:
  // so simply skip all this stuff.
  // On the other hand, every plot has an own flag to decide if it's
  // printed or not. Note that the flags are only used to do or not
  // the final print: ie all the logic is performed anyway, and just
  // the print may be avoided. This comes because in some cases there
  // dependencies between the variables used in different cases.
  //
  // But, if we are sampling the pulses, we need these stuff anyway.
  //
  if (fPlotFileExtensions.empty()     &&
      !fGenerateSPEPulseOutput        &&
      !fGeneratePreFETotalPulseOutput &&
      !fGeneratePostFETotalPulseOutput)
    return;
  // XXX Backing up the previous canvas may be needed (global variable)?
  fLog("Generating plots.", 3);
  // The canvas...
  std::stringstream canvasName;
  canvasName << "c"
             << channel.GetFrontEnd().GetModule().GetCounter().GetId() << "_"
             << channel.GetFrontEnd().GetModule().GetId() << "_"
             << channel.GetId();
  TCanvas canvas(canvasName.str().c_str() ,"Total SPE pulses", 200, 10, 700, 500);
  canvas.SetTickx();
  canvas.SetTicky();
  canvas.SetBorderMode(0);
  canvas.SetFillColor(0);
  canvas.SetFrameBorderMode(0);
  //canvas.SetGridx();
  //canvas.SetGridy();

  // The function for the total pulse prior front-end:
  std::stringstream totalPulseNamePre;
  totalPulseNamePre << "totalPulsePre_"
                    << channel.GetFrontEnd().GetModule().GetCounter().GetId() << "_"
                    << channel.GetFrontEnd().GetModule().GetId() << "_"
                    << channel.GetId();
  TF1 totalPulseFunPre(totalPulseNamePre.str().c_str(),
                       & totalPulsePreFrontEnd,
                       minTimePreFE / fUnits.GetTimeUnit(),
                       maxTimePreFE / fUnits.GetTimeUnit(),
                       1, "");
  totalPulseFunPre.SetLineColor(kRed);
  totalPulseFunPre.SetLineWidth(2);
  totalPulseFunPre.SetLineStyle(1); //  lines
  totalPulseFunPre.SetNpx(fNumPlotPoints);
  totalPulseFunPre.SetFillColor(canvas.GetFillColor()); // avoid a bounding box in the le
  totalPulseFunPre.Draw("X+");
  //
  // The individual pulses:
  unsigned int nPulse = 0;
  // Containers that own the pointers.
  OwnerVector<PulseFunctor> pFunctors;
  OwnerVector<TF1> pFunctions;
  pFunctors.reserve(si.fPulses.size());
  pFunctions.reserve(si.fPulses.size());
  for (PulseContainer::const_iterator p = si.fPulses.begin(), pe = si.fPulses.end(); p != pe; ++p) {
    PulseFunctor* pulse = new PulseFunctor(*p, fUnits.GetTimeUnit(), fUnits.GetElectricCurrentUnit());
    std::stringstream pulseName;
    pulseName << "pulse_" << nPulse++ << "_"
              << channel.GetFrontEnd().GetModule().GetCounter().GetId() << "_"
              << channel.GetFrontEnd().GetModule().GetId() << "_"
              << channel.GetId();
    TF1* pulseFun = new TF1(pulseName.str().c_str(),
                            pulse, minTimePreFE / fUnits.GetTimeUnit(),
                            maxTimePreFE / fUnits.GetTimeUnit(),
                            1, "");
    pulseFun->SetLineColor(kBlack);
    pulseFun->SetLineWidth(1);
    pulseFun->SetLineStyle(2); // lines
    pulseFun->SetNpx(fNumPlotPoints);
    pulseFun->Draw("same");
    pFunctors.push_back(pulse);
    pFunctions.push_back(pulseFun);
    if (fGenerateSPEPulseOutput) {
      std::stringstream s;
      s << "spe " << nPulse << " " << channel.GetIdMessage();
      Dump(*pulseFun, s.str());
    }
  }
  totalPulseFunPre.SetTitle("SPEs");
  totalPulseFunPre.GetXaxis()->SetTitle(("Time (" + fUnits.GetTimeName() + ")").c_str());
  totalPulseFunPre.GetYaxis()->SetTitle(("Current (" + fUnits.GetElectricCurrentName() + ")").c_str());
  canvas.Update();
  if (fGeneratePreFETotalPulseOutput) {
    std::stringstream s;
    s << "total_pre " << channel.GetIdMessage();
    Dump(totalPulseFunPre, s.str());
  }
  if (fPlotChannelPulses)
    PrintPlots(canvas, fPlotFileExtensions, Tag("channelPulses", channel)(fRunNumber));
  // ...and after the front-end pulse in voltage.
  fLog("Generating after-front-end total pulse.", 3);
  canvas.Clear();
  //
  const float rightMargin = 0.2;
  canvas.cd();
  TPad pad1("pad1", "", 0, 0, 1, 1);
  pad1.SetFrameBorderMode(0);
  pad1.SetFillColor(0);
  pad1.SetRightMargin(rightMargin);
  //pad1.SetGridx();
  //pad1.SetGridy();
  pad1.Draw();
  pad1.cd();
  /*
   **************************************************
   * XXX Plot the other signals prior the front-end.
   **************************************************
   */
  //
  // Total pulse after the front-end.
  fLog("Total pulse.", 3);
  std::stringstream totalPulseNamePost;
  totalPulseNamePost << "totalPulsePost_"
                    << channel.GetFrontEnd().GetModule().GetCounter().GetId() << "_"
                    << channel.GetFrontEnd().GetModule().GetId() << "_"
                    << channel.GetId();
  TF1 totalPulseFunPost(totalPulseNamePost.str().c_str(),
                        & totalPulsePostFrontEnd,
                        minTimePostFE / fUnits.GetTimeUnit(),
                        maxTimePostFE / fUnits.GetTimeUnit(),
                        1, "");
  totalPulseFunPost.SetLineColor(kBlack);
  totalPulseFunPost.SetLineWidth(2);
  totalPulseFunPost.SetLineStyle(1); // long dashed lines
  totalPulseFunPost.SetFillColor(canvas.GetFillColor()); // avoid a bounding box in the le
  totalPulseFunPost.SetNpx(fNumPlotPoints);
  totalPulseFunPost.GetYaxis()->SetTitleOffset(1.2);
  totalPulseFunPost.Draw();
  //
  // Threshold (on top of the previous):
  fLog("Discrimination level [", 3)
    (channel.GetThreshold() / fUnits.GetElectricPotentialUnit())(fUnits.GetElectricPotentialName())("].");
  ThresholdFunctor thresFunctor(channel.GetThreshold(), fUnits.GetElectricPotentialUnit());
  TF1 thresFun("Discrimination level",
               & thresFunctor,
               minTimePostFE / fUnits.GetTimeUnit(),
               maxTimePostFE / fUnits.GetTimeUnit(),
               0, "");
  thresFun.SetLineColor(kBlack);
  thresFun.SetLineStyle(2);
  thresFun.SetLineWidth(1); //
  thresFun.SetNpx(fNumPlotPoints);
  thresFun.SetFillColor(canvas.GetFillColor()); // avoid a bounding box in the le
  thresFun.Draw("same");
  //
  // Discriminator output
  fLog("Discriminator output.", 3);
  canvas.cd();
  TPad pad2("pad2", "", 0, 0, 1, 1);
  pad2.SetRightMargin(rightMargin);
  pad2.SetFillStyle(0);
  pad2.SetFillColor(0);
  pad2.SetFrameFillStyle(4000);
  pad2.SetFrameBorderMode(0);
  pad2.SetFrameFillColor(0);
  pad2.Draw();
  pad2.cd();
  ProxyFunctor<mdet::Channel::Discriminator> discriminatorFunctor(
    discriminator, fUnits.GetTimeUnit(), fUnits.GetElectricPotentialUnit());
  TF1 discrimFun("Discriminator output",
                 & discriminatorFunctor,
                 minTimePostFE / fUnits.GetTimeUnit(),
                 maxTimePostFE / fUnits.GetTimeUnit(),
                 1, "");
  discrimFun.SetLineColor(kBlue);
  discrimFun.SetNpx(fNumPlotPoints);
  discrimFun.SetLineStyle(2); //
  discrimFun.SetLineWidth(1); //
  discrimFun.SetFillColor(canvas.GetFillColor()); // avoid a bounding box in the legend.
   // The output is fundamentally two valued (it's a discrimination): few divisions.
  discrimFun.GetYaxis()->SetNdivisions(4);
  // Increase the distance.
  discrimFun.GetYaxis()->SetTitleOffset(1.4);
  /*
  // XXX
  // Attempt: the grid is shown according the main ticks in the left axis, it would
  // be nice to have the grid matching the main ticks of both axis.
  //
  // Force the axis range:
  // discrimFun.GetHistogram()->SetMinimum(discrimFun.GetMinimum()); // to the minimum..
  // Use static_cast instead of C-style cast for primitive type unsigned int since it has a compound
  // name, so "unsigned int(a)" needs to be written as "(unsigned int)(a)": I prefer static_cast.
  unsigned int nIntDigits =
    static_cast<unsigned int>(TMath::Log10(totalPulseFunPost.GetHistogram()->GetMaximum()));
  double spacing          = TMath::Power(10, nIntDigits);
  double grossRange       = spacing *
    static_cast<unsigned int>(totalPulseFunPost.GetHistogram()->GetMaximum() / spacing);
  unsigned int axisMult   = static_cast<unsigned int>(discrimFun.GetMaximum() / grossRange) + 1;
  std::cout << "FACT=" << nIntDigits << std::endl;
  std::cout << "FACT=" << spacing << std::endl;
  std::cout << "FACT=" << grossRange << std::endl;
  std::cout << "FACT=" << axisMult << std::endl;
  discrimFun.GetHistogram()->SetMaximum(axisMult * spacing);
  */
  // y on the right.
  discrimFun.Draw("Y+");
  // Hide x labels (once drawn)
  discrimFun.GetXaxis()->SetLabelSize(0);
  // Alternative A:
  // The samples equal to true are drawn at half the total function maximum.
  // const double trueLevel = totalPulseFun.GetMaximum(minTime, maxTime) / 2;
  // Alternative B(rian):
  // Put the trues over the threshold evaluating thresFunctor
  // Alternative C:
  // Having now an analog discriminator out, makes more to put the values
  // in the maximum of it.
  const double trueLevel = discrimFun.GetMaximum();
  // The samples.
  TVectorD samplesArray(nBinsTrace);
  // fLog("Samples: number of samples", 3)(nBinsTrace)(", number of times")(sampleTimes.size());
  // Draw the true samples with the same units as the drawn funcs.
  for (unsigned int i = 0; i < nBinsTrace; ++i)
    samplesArray[i] = trace[i] ? trueLevel : 0;
  // Samples
  /*
   * Regarding the to-C-array conversion implemented with the construct
   * &v[0] you have to see:
   * 6.2.3 - "Using Vectors as Ordinary Arrays" from
   * "The C++ Standard Library, The: A Tutorial and Reference"
   * by Nicolai M. Josuttis
   * Publisher: Addison Wesley
   *
   * and
   *
   * "C++ Standard Library Defect Report List"
   * (http://std.dkuug.dk/jtc1/sc22/wg21/docs/lwg-defects.html)
   * 69. Must elements of a vector be contiguous?
   * Section: 23.2.4 [lib.vector]  Status: TC  Submitter: Andrew Koenig  Date: 29 Jul 1998
   *
   */
  TVectorD sampleTimes_vec;
  sampleTimes_vec.Use(sampleTimes.size(),&(sampleTimes[0]));
  TGraph samples(sampleTimes_vec, samplesArray);
  samples.SetMarkerColor(kRed);
  samples.SetMarkerStyle(kFullCircle);
  samples.SetMarkerSize(1);
  samples.SetFillStyle(0); // hollow
  samples.SetFillColor(canvas.GetFillColor()); // avoid a bounding box in the legend.
  samples.SetLineColor(canvas.GetFillColor()); // avoid the line in the legend.
  samples.Draw("P");
  // Axis labels.
  totalPulseFunPost.GetXaxis()->SetTitle(("Time (" + fUnits.GetTimeName() + ")").c_str());
  totalPulseFunPost.GetYaxis()->SetTitle(
    ("Front-end output voltage (" + fUnits.GetElectricPotentialName() + ")").c_str());
  discrimFun.GetYaxis()->SetTitle(
    ("Discriminator output voltage (" + fUnits.GetElectricPotentialName() + ")").c_str());
  // Legend...
  canvas.cd(); // to make the legend's coordinates relative to the main pad.
  //TLegend leg(0.02 + 1.0 - rightMargin, 0.02, 0.18 + 1.0 - rightMargin, 0.30);
  //TLegend leg(0.02 + 1.0 - rightMargin, 0.02, 0.18 + 1.0 - rightMargin, 0.30);
  TLegend leg;
  leg.SetNColumns( 3 );
  leg.SetEntrySeparation( 0.8 );
  //leg.SetX1NDC( fStyle->GetPadLeftMargin() );
  //leg.SetX2NDC( 1 - fStyle->GetPadRightMargin() );
  leg.SetX1NDC( 0.1 );
  leg.SetX2NDC( 0.9 );
  leg.SetY1NDC( 1 - fStyle->GetPadTopMargin() );
  leg.SetY2NDC( 1 );
  leg.SetBorderSize( 0 );
  leg.SetTextFont( fStyle->GetTitleFont() );
  leg.SetTextSize( fStyle->GetTitleSize()+0.01 );
  //leg.AddEntry(&thresFun, thresFun.GetName());
  leg.AddEntry(&totalPulseFunPost, "PMT","L");
  leg.AddEntry(&discrimFun, "Discrim.", "L");
  leg.AddEntry(&samples, "Binary", "P");
  leg.SetFillColor(canvas.GetFillColor());
  leg.Draw();
  canvas.Update();
  if (fGeneratePostFETotalPulseOutput) {
    std::stringstream s;
    s << "total_post " << channel.GetIdMessage();
    Dump(totalPulseFunPost, s.str());
  }
  if (fPlotOutcome)
    PrintPlots(canvas, fPlotFileExtensions, Tag("outcome", channel)(fRunNumber));
  //
  // Input / output comparison
  totalPulsePreFrontEnd.fAbs  = true;
  totalPulsePostFrontEnd.fAbs = true;
  // Due to shift, and since we are plotting them overlapped here,
  // some commons begin/end values have to be calculated (if not two non-overlapping
  // x axis may be shown).
  const double minTimeMin = std::min(minTimePreFE, minTimePostFE) / fUnits.GetTimeUnit();
  const double maxTimeMax = std::max(maxTimePreFE, maxTimePostFE) / fUnits.GetTimeUnit();
  fLog("Generating input/output comparison plot.", 3);
  canvas.Clear();
  TPad pad3("pad3", "", 0, 0, 1, 1);
  pad3.SetFillColor(0);
  pad3.SetTickx();
  //pad3.SetGridx();
  //pad3.SetGridy();
  pad3.Draw();
  pad3.cd();
  //
  // New title:
  totalPulseFunPre.SetTitle("I/O");
  totalPulseFunPre.SetRange(minTimeMin, maxTimeMax);
  totalPulseFunPre.GetXaxis()->SetTitle(("Time (" + fUnits.GetTimeName() + ")").c_str());
  totalPulseFunPre.GetYaxis()->SetTitle(
    ("Current amplitude (" + fUnits.GetElectricCurrentName() + ")").c_str());
  totalPulseFunPre.Draw();
  TPad pad4("pad4", "", 0, 0, 1, 1);
  pad4.SetFillStyle(0);
  pad4.SetFillColor(0);
  pad4.SetFrameFillStyle(4000);
  pad4.SetFrameBorderMode(0);
  pad4.SetFrameFillColor(0);
  pad4.Draw();
  pad4.cd();
  totalPulseFunPost.SetRange(minTimeMin, maxTimeMax);
  totalPulseFunPost.GetYaxis()->SetTitle(
    ("Front-end output voltage amplitude (" + fUnits.GetElectricPotentialName() + ")").c_str());
  totalPulseFunPost.Draw("Y+");//X+ to put axis on top also.
  canvas.Update();
  if (fPlotInputOutput)
    PrintPlots(canvas, fPlotFileExtensions, Tag("input_output", channel)(fRunNumber));
  //
  // Transfer function (amplitude):
  fLog("Generating frequency response plots.", 3);
  canvas.Clear();
  canvas.SetLogx();
  canvas.SetLogy();
  TGraph ampFreqResponse(freqArray, ampArray);
  ampFreqResponse.SetMarkerStyle(21);
  ampFreqResponse.Draw("AP");
  ampFreqResponse.SetTitle("Front-end transfer function amplitude response");
  ampFreqResponse.GetXaxis()->SetTitle(("Frequency (" + fUnits.GetFrequencyName() + ")").c_str());
  ampFreqResponse.GetYaxis()->SetTitle(
    ("Amplitude (" + fUnits.GetElectricResistanceName() + ")").c_str());
  canvas.Update();
  if (fPlotTransferAmpResponse)
    PrintPlots(canvas, fPlotFileExtensions, Tag("transferAmpResponse", channel)(fRunNumber));
  //
  // Transfer function (phase):
  canvas.Clear();
  canvas.SetLogx();
  canvas.SetLogy(kFALSE);
  TGraph phaFreqResponse(freqArray, phaArray);
  phaFreqResponse.SetMarkerStyle(21);
  phaFreqResponse.Draw("AP");
  phaFreqResponse.SetTitle("Front-end transfer function phase response");
  phaFreqResponse.GetXaxis()->SetTitle(("Frequency (" + fUnits.GetFrequencyName()  + ")").c_str());
  phaFreqResponse.GetYaxis()->SetTitle(("Phase (" + fUnits.GetAngleName()  + ")").c_str());
  canvas.Update();
  if (fPlotTransferPhaseResponse)
    PrintPlots(canvas, fPlotFileExtensions, Tag("transferPhaseResponse", channel)(fRunNumber));
  //
  // Delay.
  if (phaFreqResponse.GetN() > 1) {
    canvas.Clear();
    TGraphFunctor phaFreRespFunctor(phaFreqResponse);
    const double bFreq = phaFreqResponse.GetX()[0];
    const double eFreq = phaFreqResponse.GetX()[phaFreqResponse.GetN() - 1];
    TF1 phaFreqRespTF1("Phase", &phaFreRespFunctor, bFreq, eFreq, 1, "");
    phaFreqRespTF1.SetNpx(fNumPlotPoints);
    phaFreqRespTF1.SetLineStyle(2); // long dash-dot
    // The group delay is Gd = - d Phase  / d w
    double factor = -1.0;  // So put the sign.
    factor /= utl::kTwoPi; // the former data is in terms of f (w = 2 pi f).
    // Now, since we are computing a derivative against f, the units will be time; but
    // not necesarily the common time unit used in the module: it will be the inverse
    // of the chosen frequency unit. So apply the conversion factor:
    // 1) now we have back in time units the values (since the data was divided by the
    // freq unit, as this comes from the TGraph used to plot). We have to divide because
    // the factor entered as
    //   d Phase / d ( f / fu)
    // that's the derivative against f scaled by the unit fu, so fu entered multiplying,
    // so we divide to get the value in Auger-units.
    factor /= fUnits.GetFrequencyUnit() ;
    // 2) now we divide by the desired unit to express the final value in the desired units
    // of time.
    factor /= fUnits.GetTimeUnit();
    DerivativeFunctor derivative(phaFreqRespTF1, factor);
    TF1 delayRespTF1("Delay", &derivative, bFreq, eFreq, 1, "");
    delayRespTF1.SetNpx(fNumPlotPoints);
    delayRespTF1.SetLineStyle(1); // short dash
    TPad pad5("pad5", "", 0, 0, 1, 1);
    pad5.SetFillColor(0);
    pad5.SetFrameBorderMode(0);
    //pad5.SetGridx();
    //pad5.SetGridy();
    pad5.SetLogx(kTRUE);
    pad5.SetLogy(kFALSE);
    pad5.Draw();
    pad5.cd();
    phaFreqRespTF1.Draw();
    TPad pad6("pad6", "", 0, 0, 1, 1);
    pad6.SetFillStyle(0);
    pad6.SetFillColor(0);
    pad6.SetFrameFillStyle(4000);
    pad6.SetFrameBorderMode(0);
    pad6.SetFrameFillColor(0);
    pad6.SetLogx(kTRUE);
    pad6.SetLogy(kFALSE);
    pad6.Draw();
    pad6.cd();
    delayRespTF1.Draw("Y+");
    phaFreqRespTF1.SetTitle("Phase and group delay");
    phaFreqRespTF1.GetXaxis()->SetTitle(("Frequency (" + fUnits.GetFrequencyName() + ")").c_str());
    phaFreqRespTF1.GetYaxis()->SetTitle(("Phase (" + fUnits.GetAngleName() + ")").c_str());
    delayRespTF1.GetYaxis()->SetTitle(("Time (" + fUnits.GetTimeName() + ")").c_str());
    TLegend leg(0.65, 0.70, 0.85, 0.9);
    leg.AddEntry(&phaFreqRespTF1, "Phase response");
    leg.AddEntry(&delayRespTF1, "Group delay");
    leg.Draw();
    canvas.Update();
    if (fPlotDelay)
      PrintPlots(canvas, fPlotFileExtensions, Tag("delay", channel)(fRunNumber));
  }
  //
  // FFT amplitude plot prior front-end.
  canvas.Clear();
  canvas.SetLogx();
  canvas.SetLogy();
  TGraph fftPlotPre(freqArray, preFFT);
  fftPlotPre.SetMarkerStyle(21);
  fftPlotPre.SetMarkerColor(kBlue);
  fftPlotPre.Draw("AP");
  fftPlotPre.SetTitle("Fourier amplitude spectrum prior front-end");
  fftPlotPre.GetXaxis()->SetTitle(("Frequency (" + fUnits.GetFrequencyName() + ")").c_str());
  fftPlotPre.GetYaxis()->SetTitle(("Current (" + fUnits.GetElectricCurrentName() + ")").c_str());
  canvas.Update();
  if (fPlotFftPriorAmp)
    PrintPlots(canvas, fPlotFileExtensions, Tag("fftPriorAmp", channel)(fRunNumber));
  //
  // FFT amplitude plot after front-end.
  canvas.Clear();
  canvas.SetLogx();
  canvas.SetLogy();
  TGraph fftPlotPost(freqArray, postFFT);
  fftPlotPost.SetMarkerStyle(21);
  fftPlotPost.SetMarkerColor(kRed);
  fftPlotPost.Draw("AP");
  fftPlotPost.SetTitle("Fourier amplitude spectrum after front-end.");
  fftPlotPost.GetXaxis()->SetTitle(("Frequency (" + fUnits.GetFrequencyName() + ")").c_str());
  fftPlotPost.GetYaxis()->SetTitle(("Voltage (" + fUnits.GetElectricPotentialName() + ")").c_str());
  canvas.Update();
  if (fPlotFftPostAmp)
    PrintPlots(canvas, fPlotFileExtensions, Tag("fftPostAmp", channel)(fRunNumber));
}

/*
 * Over-load for the SiPM simulations.
 *
 * We tried to make this code as similar as possible as it was for the PMTs.
 * However, the PMT method is quite a mess (> 700 lines!), and the code is way over-commented.
 *
 * Examples:
 *
 *                              "return; // done." (line 843)
 *                              ---------------------------------
 *                                // Clear the trace: (line 825)
 *                                trace.Clear();
 *
 * We are breaking the SiPM code in new methods that can (and should!) be used by the PMTs.
 * Still, we are not touching the PMT code. At this point, either the PMT code is removed or someone
 * should re-write it.
 *
 * */
void
MdOptoElectronicSimulator::ProcessPulses(const mdet::ChannelSiPM& channel,
                                         const SignalInformation& si,
                                         utl::TraceB& trace,
                                         utl::TraceD& traceAnalogical,
                                         double& span,
                                         const utl::TimeInterval& traceStart,
                                         double& minTimePreFE,
                                         double& maxTimePreFE
                                         )
{

  const mdet::FrontEndSiPM& frontEnd = channel.GetFrontEnd();

  const double sRate = frontEnd.GetMeanSampleRatePeriod();
  trace.Clear();
  traceAnalogical.Clear();
  trace.SetBinning(sRate);

  fLog("Cleaning trace and setting a binning of", 2)(sRate / fUnits.GetTimeUnit())(fUnits.GetTimeName());
  const unsigned int nBinsTrace = frontEnd.GetPreT1BufferLength() + frontEnd.GetPostT1BufferLength();
  fLog("Trace size", 2)(nBinsTrace);


  if (si.fSiPMPulses.empty() ) {
    fLog("No pulses and other sim data found for channel", 1)(channel.GetId());
    for (unsigned int i = 0; i < nBinsTrace; ++i)
      trace.PushBack(false);
    span = 0;
    return;
  }

  /************************
   *    Get trace timing  *
   ************************/
  const double pulseTimeSpan = maxTimePreFE-minTimePreFE;

  /*****************************
   *  Apply CITIROC Transfer   *
   *****************************/
  TimeTrace totalPulsePostFrontEndTrace;
  TimeTrace totalPulsePostDiscriminatorTrace;

  ApplyCITIROCTransfer(channel, si, pulseTimeSpan, minTimePreFE ,
                  totalPulsePostFrontEndTrace, totalPulsePostDiscriminatorTrace, traceAnalogical);

  /*****************************
   *  FPGA Sampling and plots  *
   *****************************/
  SampleTrace(minTimePreFE, maxTimePreFE, pulseTimeSpan / fNDiscretization, frontEnd, traceStart, totalPulsePostDiscriminatorTrace, trace);

  const double startTime = traceStart.GetInterval() - sRate*channel.GetFrontEnd().GetPreT1BufferLength();
  PlotChannel( startTime, minTimePreFE, pulseTimeSpan / fNDiscretization, channel, traceAnalogical,
                  totalPulsePostFrontEndTrace, totalPulsePostDiscriminatorTrace, trace);

}

/********************************************
 *                                          *
 *      Trace Timing                        *
 *                                          *
 ********************************************/

double
MdOptoElectronicSimulator::GetPulseTimeSpan(const SignalInformation& si, /*const utl::TimeStamp& eventTime,*/ double& minTimePreFE, double& maxTimePreFE
){

       

        minTimePreFE = std::numeric_limits<double>::max();
        maxTimePreFE = -std::numeric_limits<double>::max();


        std::unique_ptr<utl::TabularStream> spanTable;
        const unsigned int pulseSpanLogLevel = 3;
        if (fLog.GetLevel() >= pulseSpanLogLevel) {
          Init(spanTable, 3);
          *spanTable << utl::hline
                                 << "# pe" << utl::endc // 1
                                 << "min (" << fUnits.GetTimeName() << ")" << utl::endc // 2
                                 << "max (" << fUnits.GetTimeName() << ")" << utl::endr;// 3
        }
        unsigned int nPE = 0;
        for (SiPMPulseContainer::const_iterator p = si.fSiPMPulses.begin(), e = si.fSiPMPulses.end(); p != e; ++p) {
          minTimePreFE = std::min(minTimePreFE, p->LowerLimit());
          maxTimePreFE = std::max(maxTimePreFE, p->UpperLimit());
          if (fLog.GetLevel() >= (pulseSpanLogLevel+1)) {
                *spanTable << utl::hline
                                   << ++nPE << utl::endc
                                   << fLog << p->LowerLimit() / fUnits.GetTimeUnit() << utl::endc
                                   << fLog << p->UpperLimit() / fUnits.GetTimeUnit() << utl::endr;
          }
        }


        // The span covers both the range of pulses as well that of the analogic trace.
        const double pulseTimeSpan = maxTimePreFE - minTimePreFE;
        if (fLog.GetLevel() >= pulseSpanLogLevel) {
          *spanTable << utl::hline
                                 << "all" << utl::endc
                                 << fLog << minTimePreFE / fUnits.GetTimeUnit() << utl::endc
                                 << fLog << maxTimePreFE / fUnits.GetTimeUnit() << utl::endr
                                 << utl::hline;
          fLog("PE pulses and other signals span:", pulseSpanLogLevel);
          fLog(spanTable->Str(), pulseSpanLogLevel);
        }
        return pulseTimeSpan;
}

/****************************************
 *                                      *
 * CITIROC Transfer                     *
 *                                      *
 ****************************************/

void
MdOptoElectronicSimulator::ApplyCITIROCTransfer(const mdet::ChannelSiPM& channel, const SignalInformation& si, const double pulseTimeSpan, double minTimePreFE,
                 TimeTrace& totalPulsePostFrontEndTrace,
                TimeTrace& totalPulsePostDiscriminatorTrace, utl::TraceD& traceAnalogical){

        TotalPulseFunctor<SiPMPulseContainer> totalPulsePreFrontEnd(si.fSiPMPulses,
                                                                  fUnits.GetTimeUnit(),
                                                                  fUnits.GetElectricPotentialUnit());

        

        fLog("Applying FFT...", 3);

        /************************************************
         *First transform the trace into Fourier space  *
         ************************************************/

        utl::FFTDataContainer<utl::Trace, TimeTrace::ValueType, FrequencyTrace::ValueType> fft;

        double binning = pulseTimeSpan / fNDiscretization;
        TimeTrace& timeTrace = fft.GetTimeSeries();
        fLog("Time discretization binning ", 4)(binning / fUnits.GetTimeUnit())(fUnits.GetTimeName());
        timeTrace.SetBinning(binning);
        fLog("Discretizing total pulses. Number of samples:", 3, false);

        for (unsigned int i = 0; i < fNDiscretization || fLog(i, 3)('.'); ++i) {
          double time  = minTimePreFE + i * binning;
          double value = totalPulsePreFrontEnd(time);

          fft.GetTimeSeries().PushBack(value);
          traceAnalogical.PushBack(value); //Already prepare the Input trace for the integrator
        }

        /************************************************
         *Apply pre-amplifier and fast shaper transfer  *
         ************************************************/

        fLog("Applying transfer function in frequency spectrum.", 3);
        FrequencyTrace& freqTrace            = fft.GetFrequencySpectrum();
        const FrequencyTrace::SizeType nFreq = fft.GetConstFrequencySpectrum().GetSize();
        const double freqBinning             = fft.GetConstFrequencySpectrum().GetBinning();

        fLog("Generating frequency response.", 3);
        const unsigned int tabFFTLogLevel = 4;
        std::unique_ptr<utl::TabularStream> tabFFT;
        if (fLog.GetLevel() >= tabFFTLogLevel) {
          Init(tabFFT, 4);
          *tabFFT <<  utl::hline
                         << "Size"        << utl::endc //1
                         << "Binning ("   << fUnits.GetFrequencyName() << ")" << utl::endc //2
                         << "1st freq ("  << fUnits.GetFrequencyName() << ")" << utl::endc //3
                         << "Last freq (" << fUnits.GetFrequencyName() << ")" << utl::endr //4
                         << nFreq         << utl::endc // 1
                         << fLog << (freqBinning                      / fUnits.GetFrequencyUnit()) << utl::endc // 2
                         << fLog << (fft.GetFrequencyOfBin(0)         / fUnits.GetFrequencyUnit()) << utl::endc // 3
                         << fLog << (fft.GetFrequencyOfBin(nFreq - 1) / fUnits.GetFrequencyUnit()) << utl::endr // 4
                         << utl::hline;
          fLog(tabFFT->Str(), tabFFTLogLevel);
        }
        for (FrequencyTrace::SizeType freqBin = 0; freqBin < nFreq; ++freqBin) {
          double freq = fft.GetFrequencyOfBin(freqBin) / fUnits.GetFrequencyUnit();//From GHz(ns) to MHz
          std::complex<double> c = channel.ComputeTransfer(freq);
          freqTrace[ freqBin ] *= c;

        }
        fLog("Obtaining the amplified pulses in time-space...", 3);
        totalPulsePostFrontEndTrace = fft.GetTimeSeries();

        /**************************************************
         *      Simulate discriminator response   *
         **************************************************/

        const TimeTrace::SizeType nTime = totalPulsePostFrontEndTrace.GetSize();
        for (TimeTrace::SizeType timeBin = 0; timeBin < nTime; ++timeBin) {
                  double discriminatorOuput = channel.ComputeDiscriminator(totalPulsePostFrontEndTrace[timeBin], binning);
                  totalPulsePostDiscriminatorTrace.PushBack(discriminatorOuput);
        }
}

/********************************
 *                              *
 * FPGA Sampling                *
 *                              *
 ********************************/

void
MdOptoElectronicSimulator::SampleTrace(double minTimePostFE, double maxTimePostFE, double binning, const mdet::FrontEndSiPM& frontEnd,
                const utl::TimeInterval& traceStart, TimeTrace& totalPulsePostDiscriminatorTrace, utl::TraceB& trace){

          fLog("Creating sampler object.", 3);
          mdet::FrontEndSiPM::Sampler sampler = frontEnd.MakeSampler();
          const double sRate = frontEnd.GetMeanSampleRatePeriod();
          const unsigned int nBinsTrace = frontEnd.GetPreT1BufferLength() + frontEnd.GetPostT1BufferLength();

          std::vector<double> sampleTimes;
          // sampleTimes.reserve(nBinsTrace);
          fLog("Sampling pulses.", 3);

          const double startTime = traceStart.GetInterval() - sRate* frontEnd.GetPreT1BufferLength();
          const double stopTime  = traceStart.GetInterval() + sRate* frontEnd.GetPostT1BufferLength();

          //const double startTime  = traceStart.GetInterval();
          //const double stopTime   = startTime +
          //  sRate * (frontEnd.GetPreT1BufferLength() + frontEnd.GetPostT1BufferLength());

          if (startTime < minTimePostFE) {
            fLog("Samples false in the range (startTime < minTimePostFE)", 4)
                (startTime / fUnits.GetTimeUnit())(fUnits.GetTimeName())
                ("to")
                (minTimePostFE / fUnits.GetTimeUnit())(fUnits.GetTimeName())
                ("as no signal in it.");
          }
          //fLog("Samples:", 4);

          for (double sampleTime = startTime; trace.GetSize() < nBinsTrace || fLog('\n', 4);
            sampleTime += sRate) {

            bool s;
            if (sampleTime < minTimePostFE || maxTimePostFE < sampleTime) {
              s = false;
            } else {
              int bin = (int)((sampleTime-minTimePostFE)/binning);
              double o = totalPulsePostDiscriminatorTrace[bin];
              s = sampler(o);
              // Only log the relevant part:
              fLog(s, 4, false)('@')(sampleTime / fUnits.GetTimeUnit())(fUnits.GetTimeName())('|');
            }
            trace.PushBack(s);
            sampleTimes.push_back(sampleTime);
          }

          if (maxTimePostFE < stopTime) {
            // See inner body of the previous loop (then condition).
            fLog("Samples false in the range (maxTimePostFE < stopTime) ", 4)
                (maxTimePostFE / fUnits.GetTimeUnit())(fUnits.GetTimeName())
                ("to")
                (stopTime / fUnits.GetTimeUnit())(fUnits.GetTimeName())
                ("as no signal in it.");
          }

}

void
MdOptoElectronicSimulator::InjectDigitalNoise(const mdet::Module& module, mevt::Module& evtModule){

        const mdet::FrontEndSiPM& frontEnd = module.GetFrontEndSiPM();
        short bin = frontEnd.GetInjectDigitalNoiseBin();
        if(bin>=0){

                unsigned short width = frontEnd.GetInjectDigitalNoiseWidth();
                unsigned short channelId = frontEnd.GetInjectDigitalNoiseChannel();
                if(evtModule.HasChannel(channelId)){
                        mevt::Channel& channel = evtModule.GetChannel(channelId);
                        if (!channel.HasTrace())
                                channel.MakeTrace();
                        utl::TraceB& trace = channel.GetTrace();
                        unsigned short limit = std::min(width+bin, (int)trace.GetSize());

                        for(short i = bin; i<limit; i++){
                                trace[i] = true;
                        }
                }
        }

}

void
MdOptoElectronicSimulator::PlotChannel(const double traceStartTime, const double minTimePreFE, const double binning, const mdet::ChannelSiPM& channel,
                utl::TraceD& traceAnalogical, TimeTrace& totalPulsePostFrontEndTrace, TimeTrace& totalPulsePostDiscriminatorTrace, utl::TraceB& trace){

  if (fPlotFileExtensions.empty()     &&
      !fGenerateSPEPulseOutput        &&
      !fGeneratePreFETotalPulseOutput &&
      !fGeneratePostFETotalPulseOutput)
    return;

  TVectorD analogicalTime(traceAnalogical.GetSize());
  TVectorD traceAnalogicalVec(traceAnalogical.GetSize());
  TVectorD totalPulsePostFrontEndVec(traceAnalogical.GetSize());
  TVectorD totalPulsePostDiscriminatorVec(traceAnalogical.GetSize());

  TVectorD digitalTime(trace.GetSize());
  TVectorD traceFPGA(trace.GetSize());

  const double sRate = channel.GetFrontEnd().GetMeanSampleRatePeriod();
  const double discriminatorHiLevel = 3.0;//channel.GetDiscriminatorHiLevel();
  const double scaleFactor = 50.0;//To visualize the SiPM pulse;

  for(unsigned int i=0; i<traceAnalogical.GetSize(); i++){
          analogicalTime[i] = (i*binning + minTimePreFE);
          traceAnalogicalVec[i] = -1*traceAnalogical[i]*scaleFactor / fUnits.GetElectricPotentialUnit(); //Invert only for visual purpose.
          totalPulsePostFrontEndVec[i] = totalPulsePostFrontEndTrace[i] / fUnits.GetElectricPotentialUnit();
          totalPulsePostDiscriminatorVec[i] = totalPulsePostDiscriminatorTrace[i] / scaleFactor / fUnits.GetElectricPotentialUnit();
  }


  for(unsigned int i=0; i<trace.GetSize(); i++){
          digitalTime[i] = (i*sRate) + traceStartTime;
          traceFPGA[i] = ((double)trace[i]*discriminatorHiLevel) / scaleFactor / fUnits.GetElectricPotentialUnit();//rescale
  }

  fLog("Generating plots.", 3);
  // The canvas...
  std::stringstream canvasName;
  canvasName << "c"
             << channel.GetFrontEnd().GetModule().GetCounter().GetId() << "_"
             << channel.GetFrontEnd().GetModule().GetId() << "_"
             << channel.GetId();
  TCanvas canvas(canvasName.str().c_str() ,"Channel Output", 200, 10, 700, 500);
  canvas.SetTickx();
  canvas.SetTicky();
  canvas.SetBorderMode(0);
  canvas.SetFillColor(0);
  canvas.SetFrameBorderMode(0);
  //canvas.SetGridx();
  //canvas.SetGridy();

  TGraph totalPulsePrePlot(analogicalTime, traceAnalogicalVec);
  totalPulsePrePlot.SetLineWidth(3);
  totalPulsePrePlot.SetLineColor(kBlue);

  TGraph pulsePostFrontEndPlot(analogicalTime, totalPulsePostFrontEndVec);
  pulsePostFrontEndPlot.SetLineWidth(3);
  pulsePostFrontEndPlot.SetLineColor(kOrange);

  TGraph pulsePostDiscriminatorPlot(analogicalTime, totalPulsePostDiscriminatorVec);
  pulsePostDiscriminatorPlot.SetLineWidth(3);
  pulsePostDiscriminatorPlot.SetLineColor(kBlack);

  TGraph pulsePostFPGAPlot(digitalTime, traceFPGA);
  pulsePostFPGAPlot.SetMarkerColor(kRed);
  pulsePostFPGAPlot.SetMarkerStyle(kFullCircle);
  pulsePostFPGAPlot.SetMarkerSize(1);
  pulsePostFPGAPlot.SetTitle("");

  pulsePostFPGAPlot.GetXaxis()->SetTitle(("Time (" + fUnits.GetTimeName() + ")").c_str());
  pulsePostFPGAPlot.GetYaxis()->SetTitle(("Voltage (" + fUnits.GetElectricPotentialName() + ")").c_str());
  canvas.Update();
  canvas.cd();

  pulsePostFPGAPlot.Draw("AP");
  pulsePostFrontEndPlot.Draw("SAME");
  totalPulsePrePlot.Draw("SAME");
  pulsePostDiscriminatorPlot.Draw("SAME");

  pulsePostFPGAPlot.GetXaxis()->SetRangeUser(minTimePreFE-50,minTimePreFE+150);
  pulsePostFPGAPlot.GetYaxis()->SetRangeUser(0,discriminatorHiLevel / scaleFactor / fUnits.GetElectricPotentialUnit());

  TLegend leg;
  leg.SetNColumns( 2 );
  leg.SetEntrySeparation( 0.8 );
  //leg.SetX1NDC( fStyle->GetPadLeftMargin() );
  //leg.SetX2NDC( 1 - fStyle->GetPadRightMargin() );
  leg.SetX1NDC( 0.1 );
  leg.SetX2NDC( 0.9 );
  leg.SetY1NDC( 1 - fStyle->GetPadTopMargin() );
  leg.SetY2NDC( 1 );
  leg.SetBorderSize( 0 );
  leg.SetTextFont( fStyle->GetTitleFont() );
  leg.SetTextSize( fStyle->GetTitleSize()+0.01 );
  //leg.AddEntry(&thresFun, thresFun.GetName());
  leg.AddEntry(&totalPulsePrePlot, "SiPM (#times50)","L");
  leg.AddEntry(&pulsePostFrontEndPlot, "Shaper", "L");
  leg.AddEntry(&pulsePostDiscriminatorPlot, "Discrim. (/50)", "L");
  leg.AddEntry(&pulsePostFPGAPlot, "Binary ", "P");
  leg.SetFillColor(canvas.GetFillColor());
  leg.Draw();

  canvas.Update();

  if (fPlotOutcome)
    PrintPlots(canvas, fPlotFileExtensions, Tag("outcome", channel)(fRunNumber));

}

/*
 *
 *Simulate the Integrator electronics. 1 run per module
 *
 * */
void
MdOptoElectronicSimulator::ProcessPulsesIntegrator(const mdet::Module& module,
                                                   std::vector<utl::TraceD> analogicalTraces,
                                                   utl::TraceUSI& traceIntegratorA,
                                                   utl::TraceUSI& traceIntegratorB,
                                                   const utl::TimeInterval& traceStart,
                                                   double& minTimePreFE,
                                                   double& maxTimePreFE)
{

  const mdet::BackEndSiPM& backEnd   = module.GetBackEndSiPM();
  const mdet::FrontEndSiPM& frontEnd = module.GetFrontEndSiPM();

  const double sRateCounter    = frontEnd.GetMeanSampleRatePeriod();
  const double sRateIntegrator = frontEnd.GetSampleTimeADC();
  traceIntegratorA.Clear();
  traceIntegratorB.Clear();
  traceIntegratorA.SetBinning(sRateIntegrator);
  traceIntegratorB.SetBinning(sRateIntegrator);
  fLog("Cleaning Integrator Traces and setting a binning of", 2)(sRateIntegrator / fUnits.GetTimeUnit())(fUnits.GetTimeName());
  const unsigned int nBinsTrace = (frontEnd.GetPreT1BufferLength() + frontEnd.GetPostT1BufferLength())*sRateCounter/sRateIntegrator;
  fLog("Trace size", 2)(nBinsTrace);
  if (analogicalTraces.empty()) {
    fLog("No pulses data for this module", 2)(module.GetId());
    for (unsigned int i = 0; i < nBinsTrace; ++i){
        traceIntegratorA.PushBack(0);
        traceIntegratorB.PushBack(0);
    }
    return;
  }
  for(unsigned int i = 0; i<analogicalTraces.size(); i++){
        if(analogicalTraces[i].GetSize() != fNDiscretization){
                fLog("Channel traces are corrupted", 1)(module.GetId());
                return;
        }
  }

  /**********************************************
   *   Apply Transfer from first adder to ADC   *
   **********************************************/
  TimeTrace traceAfterADCLowGain;
  TimeTrace traceAfterADCHighGain;
  TVectorD totalAnalogicalInput(fNDiscretization);

  if(fIntegratorSimType == eSimplified){
          ApplyBackEndTransfer(backEnd, maxTimePreFE, minTimePreFE, traceAfterADCLowGain, traceAfterADCHighGain, totalAnalogicalInput, analogicalTraces);
  }else{
          ApplyBackEndTransferWStepSaturation(backEnd, maxTimePreFE, minTimePreFE, traceAfterADCLowGain, traceAfterADCHighGain, totalAnalogicalInput, analogicalTraces);

  }
   /*****************************
   *  ADC Sampling and plots  *
   *****************************/
  SampleTraceADC(minTimePreFE, maxTimePreFE, frontEnd, traceStart, traceAfterADCLowGain, traceAfterADCHighGain, traceIntegratorA, traceIntegratorB);

  //Pre and Post T1Buffers are given in terms of binary samples
  const double sRateFactor = frontEnd.GetMeanSampleRatePeriod()/frontEnd.GetSampleTimeADC();

  const double startTime = traceStart.GetInterval() - frontEnd.GetSampleTimeADC() * frontEnd.GetPreT1BufferLength() * sRateFactor;

  PlotIntegrator(startTime, minTimePreFE, (maxTimePreFE-minTimePreFE) / fNDiscretization,
                  frontEnd, totalAnalogicalInput, traceAfterADCLowGain, traceAfterADCHighGain, traceIntegratorA, traceIntegratorB, frontEnd.GetDelayBinaryADC());
}


/************************************************************************
 *                                                                      *
 * Integrator BackEnd Transfer Apply intermediate steps saturation      *
 *                                                                      *
 ***********************************************************************/

void
MdOptoElectronicSimulator::ApplyBackEndTransferWStepSaturation(const mdet::BackEndSiPM& backEnd, const double maxTimePreFE, const double minTimePreFE,
                TimeTrace& traceAfterADCLowGain, TimeTrace& traceAfterADCHighGain, TVectorD& totalAnalogicalInput, const std::vector<utl::TraceD> analogicalTraces){

    fLog("Applying FFT Integrator...", 3);

    /**************************************************************************
     *First add every 16 channels and transform the trace into Fourier space  *
     **************************************************************************/

    double binning = (maxTimePreFE-minTimePreFE) / fNDiscretization;
    const double numberOfChannelsToGroup = backEnd.GetNumberOfChannelsToGroup();
    int nroChannels = analogicalTraces.size();
    double nroOfAdders = nroChannels/numberOfChannelsToGroup;

    std::vector<TimeTrace> tracesPostFirstAdder;
    for (int adder = 0; adder < nroOfAdders; adder++) { //4

        int fromChannel = adder*numberOfChannelsToGroup;
        int toChannel = (adder+1)*numberOfChannelsToGroup;
        toChannel = std::min(toChannel, nroChannels);
        fLog("Adding pulses. From/To channel:", 3)(fromChannel)(toChannel);

        utl::FFTDataContainer<utl::Trace, TimeTrace::ValueType, FrequencyTrace::ValueType> fft;
        fLog("Time discretization binning ", 4)(binning / fUnits.GetTimeUnit())(fUnits.GetTimeName());
        fft.GetTimeSeries().SetBinning(binning);

        for (unsigned int i = 0; i < fNDiscretization; ++i) {
                double value = 0;
                //not efficient, but there are not that many iterations.
                for (int channel = fromChannel; channel < toChannel; channel++) { //16
                  value += analogicalTraces[channel][i];
                }
                fft.GetTimeSeries().PushBack(value);
                totalAnalogicalInput[i] += value;
        }

        ApplyTransferBlock(fft, backEnd, backEnd.eFirstAdder);

        fLog("Obtaining the amplified pulses in time-space...", 3);
        tracesPostFirstAdder.push_back(fft.GetTimeSeries());
    }

    utl::FFTDataContainer<utl::Trace, TimeTrace::ValueType, FrequencyTrace::ValueType> fft;
    fft.GetTimeSeries().SetBinning(binning);
    for(unsigned int i = 0; i<fNDiscretization; i++){
        double value = 0;
        for(unsigned int adder = 0; adder < tracesPostFirstAdder.size(); adder++){
                value += backEnd.ApplySaturation(tracesPostFirstAdder[adder][i],backEnd.eFirstAdder);
        }
        fft.GetTimeSeries().PushBack(value);
    }

    /********************************
     *Second Adder                      *
     ********************************/
    ApplyTransferBlock(fft, backEnd, backEnd.eSecondAdder);
    TimeTrace& tracePreLGChannel = fft.GetTimeSeries();

    utl::FFTDataContainer<utl::Trace, TimeTrace::ValueType, FrequencyTrace::ValueType> fftHGChannel;
    fftHGChannel.GetTimeSeries().SetBinning(binning);
    TimeTrace& tracePreHGChannel = fftHGChannel.GetTimeSeries();

    for(unsigned int i = 0; i<fNDiscretization; i++){

        tracePreLGChannel[i] = backEnd.ApplySaturation(tracePreLGChannel[i],backEnd.eSecondAdder);
        tracePreHGChannel.PushBack(tracePreLGChannel[i]);
    }

    /********************************************************
     *Amplification (low and high gain) and ADC         *
     *******************************************************/
    for(unsigned int i = 0; i<fNDiscretization; i++){
        /*fft.GetTimeSeries()[i] = backEnd.ApplySaturation(fft.GetTimeSeries()[i],backEnd.eLowGainAmplifier);
        fftHGChannel.GetTimeSeries()[i] = backEnd.ApplySaturation(fftHGChannel.GetTimeSeries()[i],backEnd.eHighGainAmplifier);*/
        traceAfterADCLowGain.PushBack(backEnd.ApplySaturation(fft.GetTimeSeries()[i],backEnd.eLowGainAmplifier));
        traceAfterADCHighGain.PushBack(backEnd.ApplySaturation(fftHGChannel.GetTimeSeries()[i],backEnd.eHighGainAmplifier));
    }

    ApplyTransferBlock(fft, backEnd, backEnd.eADC);
    ApplyTransferBlock(fftHGChannel, backEnd, backEnd.eADC);
    for(unsigned int i = 0; i<fNDiscretization; i++){
        traceAfterADCLowGain.PushBack(backEnd.ApplySaturation(fft.GetTimeSeries()[i],backEnd.eADC));
        traceAfterADCHighGain.PushBack(backEnd.ApplySaturation(fftHGChannel.GetTimeSeries()[i],backEnd.eADC));
    }
}

/*******************************************************************************
 *                                                             *
 *   Integrator BackEnd Transfer Apply transfer and saturation of all blocks   *
 *                                                                             *
 *******************************************************************************/

void
MdOptoElectronicSimulator::ApplyBackEndTransfer(const mdet::BackEndSiPM& backEnd, const double maxTimePreFE, const double minTimePreFE,
                TimeTrace& traceAfterADCLowGain, TimeTrace& traceAfterADCHighGain, TVectorD& totalAnalogicalInput, const std::vector<utl::TraceD> analogicalTraces){
        fLog("Applying FFT Integrator...", 3);

        /****************************************************************************
         *      First add every channel and transform the trace into Fourier space  *
         ****************************************************************************/

        double binning = (maxTimePreFE-minTimePreFE) / fNDiscretization;
        int nroChannels = analogicalTraces.size();

        utl::FFTDataContainer<utl::Trace, TimeTrace::ValueType, FrequencyTrace::ValueType> fft;
        fLog("Time discretization binning ", 4)(binning / fUnits.GetTimeUnit())(fUnits.GetTimeName());
        fft.GetTimeSeries().SetBinning(binning);

        for (unsigned int i = 0; i < fNDiscretization; ++i) {
                double value = 0;
                //not efficient, but there are not that many iterations.
                for (int channel = 0; channel < nroChannels; channel++) {
                         value += analogicalTraces[channel][i];
                }
                fft.GetTimeSeries().PushBack(value);
                totalAnalogicalInput[i] += value;
        }
        /********************************
         *      1st and 2nd adders      *
         ********************************/

        /*ApplyTransferBlock(fft, backEnd, backEnd.eFirstAdder);
        end_time = std::chrono::system_clock::to_time_t(std::chrono::system_clock::now());
        std::cout << "Apply second adder " << std::ctime(&end_time);
        ApplyTransferBlock(fft, backEnd, backEnd.eSecondAdder);*/

        utl::FFTDataContainer<utl::Trace, TimeTrace::ValueType, FrequencyTrace::ValueType> fftHG;
        ApplyTransferBlocks(fft, fftHG, backEnd);
        TimeTrace& tracePreLGChannel = fft.GetTimeSeries();
        TimeTrace& tracePreHGChannel = fftHG.GetTimeSeries();

        for(unsigned int i = 0; i<fNDiscretization; i++){
                traceAfterADCLowGain.PushBack(backEnd.ApplySaturation(tracePreLGChannel[i],backEnd.eLowGainAmplifier));
                traceAfterADCHighGain.PushBack(backEnd.ApplySaturation(tracePreHGChannel[i],backEnd.eHighGainAmplifier));
        }
}

void
MdOptoElectronicSimulator::ApplyTransferBlock(utl::FFTDataContainer<utl::Trace, TimeTrace::ValueType, FrequencyTrace::ValueType>& fft,
                const mdet::BackEndSiPM& backEnd, BackEndSiPM::TransferStep step){

        /********************************************************
        *       ApplyfirstAddertransferandchecksaturation*
        ********************************************************/

        fLog("Applying transfer function in frequency space.",3);
        FrequencyTrace& freqTrace = fft.GetFrequencySpectrum();
        const FrequencyTrace::SizeType nFreq = fft.GetConstFrequencySpectrum().GetSize();
        const double freqBinning = fft.GetConstFrequencySpectrum().GetBinning();

        fLog("Generating frequency response.",3);
        const unsigned int tabFFTLogLevel = 4;
        std::unique_ptr<utl::TabularStream> tabFFT;
        if (fLog.GetLevel() >= tabFFTLogLevel){
          Init(tabFFT,4);
          *tabFFT<<utl::hline
                 <<"Size"<<utl::endc//1
                 <<"Binning("<<fUnits.GetFrequencyName()<<")"<<utl::endc//2
                 <<"1st freq("<<fUnits.GetFrequencyName()<<")"<<utl::endc//3
                 <<"Last freq("<<fUnits.GetFrequencyName()<<")"<<utl::endr//4
                 <<nFreq<<utl::endc//1
                 <<fLog<<(freqBinning/fUnits.GetFrequencyUnit())<<utl::endc//2
                 <<fLog<<(fft.GetFrequencyOfBin(0)/fUnits.GetFrequencyUnit())<<utl::endc//3
                 <<fLog<<(fft.GetFrequencyOfBin(nFreq-1)/fUnits.GetFrequencyUnit())<<utl::endr//4
                 <<utl::hline;
          fLog(tabFFT->Str(),tabFFTLogLevel);
        }
        for (FrequencyTrace::SizeType freqBin=0; freqBin<nFreq; ++freqBin){
           double freq = fft.GetFrequencyOfBin(freqBin) / fUnits.GetFrequencyUnit();//To MHz
           //double preAmp= std::abs(freqTrace[freqBin])/fUnits.GetElectricCurrentUnit();
           std::complex<double> c = backEnd.ComputeTransfer(freq, step);
           //const std::complex<double> intervention(-5,0);
           freqTrace[freqBin] *= c;
        }
}

void
MdOptoElectronicSimulator::ApplyTransferBlocks(utl::FFTDataContainer<utl::Trace, TimeTrace::ValueType, FrequencyTrace::ValueType>& fft,
                utl::FFTDataContainer<utl::Trace, TimeTrace::ValueType, FrequencyTrace::ValueType>& fftHG,
                const mdet::BackEndSiPM& backEnd){

        /*******************************************************
        *       ApplyfirstAddertransferandchecksaturation      *
        ********************************************************/

        fLog("Applying transfer function in frequency space.",3);
        FrequencyTrace& freqTrace = fft.GetFrequencySpectrum();
        FrequencyTrace& freqTraceHG = fftHG.GetFrequencySpectrum();
        const FrequencyTrace::SizeType nFreq = fft.GetConstFrequencySpectrum().GetSize();
        const double freqBinning = fft.GetConstFrequencySpectrum().GetBinning();
        freqTraceHG.SetBinning(freqBinning);

        fLog("Generating frequency response.",3);
        const unsigned int tabFFTLogLevel = 4;
        std::unique_ptr<utl::TabularStream> tabFFT;
        if (fLog.GetLevel() >= tabFFTLogLevel){
          Init(tabFFT,4);
          *tabFFT<<utl::hline
                 <<"Size"<<utl::endc//1
                 <<"Binning("<<fUnits.GetFrequencyName()<<")"<<utl::endc//2
                 <<"1st freq("<<fUnits.GetFrequencyName()<<")"<<utl::endc//3
                 <<"Last freq("<<fUnits.GetFrequencyName()<<")"<<utl::endr//4
                 <<nFreq<<utl::endc//1
                 <<fLog<<(freqBinning/fUnits.GetFrequencyUnit())<<utl::endc//2
                 <<fLog<<(fft.GetFrequencyOfBin(0)/fUnits.GetFrequencyUnit())<<utl::endc//3
                 <<fLog<<(fft.GetFrequencyOfBin(nFreq-1)/fUnits.GetFrequencyUnit())<<utl::endr//4
                 <<utl::hline;
          fLog(tabFFT->Str(),tabFFTLogLevel);
        }
        for (FrequencyTrace::SizeType freqBin=0; freqBin<nFreq; ++freqBin){
           double freq = fft.GetFrequencyOfBin(freqBin) / fUnits.GetFrequencyUnit();//To MHz
           std::complex<double> cLG =  backEnd.ComputeTransfer(freq, backEnd.eSimplifiedLG);
           std::complex<double> cHG =  backEnd.ComputeTransfer(freq, backEnd.eSimplifiedHG);
           freqTraceHG.PushBack(freqTrace[freqBin]*cHG);
           freqTrace[freqBin] *= cLG;
        }
}

/********************************
 *                              *
 * ADC Sampling                 *
 *                              *
 ********************************/
void
MdOptoElectronicSimulator::SampleTraceADC(const double minTimePostFE, const double maxTimePostFE, const mdet::FrontEndSiPM& frontEnd,
                const utl::TimeInterval& traceStart, TimeTrace& traceAfterADCLowGain, TimeTrace& traceAfterADCHighGain,
                utl::TraceUSI& traceIntegratorA, utl::TraceUSI& traceIntegratorB){

      const double delayBinaryADC = frontEnd.GetDelayBinaryADC();
      const double sRate = frontEnd.GetSampleTimeADC();
      const unsigned int nBinsTrace = (int)((frontEnd.GetPreT1BufferLength() + frontEnd.GetPostT1BufferLength())*frontEnd.GetMeanSampleRatePeriod()/sRate);
      const double binning = (maxTimePostFE-minTimePostFE)/fNDiscretization;

      // sampleTimes.reserve(nBinsTrace);
      fLog("ADC analogical pulses (HG/LG).", 3);

      //Pre and Post buffer are given in terms of binary sampling rate
       const double startTime = traceStart.GetInterval() - frontEnd.GetMeanSampleRatePeriod() * frontEnd.GetPreT1BufferLength();
       const double stopTime  = traceStart.GetInterval() + frontEnd.GetMeanSampleRatePeriod() * frontEnd.GetPostT1BufferLength();

      //const double startTime  = traceStart.GetInterval();
      //const double stopTime   = startTime + sRate * nBinsTrace;

      if (startTime + delayBinaryADC < minTimePostFE) {
        fLog("No signal in the range", 4)
            (startTime / fUnits.GetTimeUnit())(fUnits.GetTimeName())
            ("to")
            ("(")(minTimePostFE / fUnits.GetTimeUnit())("+")(delayBinaryADC / fUnits.GetTimeUnit())(")=")((minTimePostFE+delayBinaryADC)/fUnits.GetTimeUnit())(fUnits.GetTimeName())
            (".");
      }
      fLog("ADC High/Low gain channels:", 4);

      //double chargeLG = 0.0;
      //double chargeHG = 0.0;
      unsigned short offsetLG = frontEnd.GetADCCounts(frontEnd.GetModule().GetBackEndSiPM().GetLowGainAmplifierOffset());
      unsigned short offsetHG = frontEnd.GetADCCounts(frontEnd.GetModule().GetBackEndSiPM().GetHighGainAmplifierOffset());

      for (double sampleTime = startTime; traceIntegratorA.GetSize() < nBinsTrace || fLog('\n', 4)(sampleTime);
        sampleTime += sRate) {

        unsigned short adcLG;
        unsigned short adcHG;
        //Introduce delay between counter and integrator
        if (sampleTime < minTimePostFE + delayBinaryADC || maxTimePostFE + delayBinaryADC < sampleTime) {
          adcLG = offsetLG;
          adcHG = offsetHG;
        } else {

          int bin = (int)((sampleTime-minTimePostFE-delayBinaryADC)/binning);
          double signal1 = traceAfterADCLowGain[bin];
          double signal2 = traceAfterADCHighGain[bin];
          adcLG = frontEnd.GetADCCounts(signal1);
          adcHG = frontEnd.GetADCCounts(signal2);
          //chargeLG += adcLG-offsetLG;
          //chargeHG += adcHG-offsetHG;

          // Only log the relevant part:
          fLog(adcLG, 4, false)('-')(adcHG)('@')(sampleTime / fUnits.GetTimeUnit())(fUnits.GetTimeName())('|');
        }
        if(fIncludeBaseLineFluctuationIntegrator){
        adcLG += frontEnd.GetADCBaseLineFluctuationLG();
        adcHG += frontEnd.GetADCBaseLineFluctuationHG();
        }
        traceIntegratorA.PushBack(adcLG);
        traceIntegratorB.PushBack(adcHG);
      }

      if (maxTimePostFE + delayBinaryADC < stopTime) {
        // See inner body of the previous loop (then condition).
        fLog("Samples false in the range ", 4)
            ("(")(maxTimePostFE / fUnits.GetTimeUnit())("+")(delayBinaryADC / fUnits.GetTimeUnit())(")=")((maxTimePostFE+delayBinaryADC)/fUnits.GetTimeUnit())(fUnits.GetTimeName())
            ("to")
            (stopTime / fUnits.GetTimeUnit())(fUnits.GetTimeName())
            ("as no signal in it.");
      }
}

void
MdOptoElectronicSimulator::PlotIntegrator(const double traceStartTime, const double minTimePreFE, const double binning, const mdet::FrontEndSiPM& frontEnd,
                TVectorD& traceAnalogical, TimeTrace& traceIntegratorAAmplifier, TimeTrace& traceIntegratorBAmplifier,
                utl::TraceUSI& traceIntegratorA, utl::TraceUSI& traceIntegratorB, const double delay){

          if (fPlotFileExtensions.empty()     &&
              !fGeneratePreFETotalPulseOutput &&
              !fGeneratePostFETotalPulseOutput)
            return;

          TVectorD traceIntegratorAVec(traceIntegratorA.GetSize());
          TVectorD traceIntegratorBVec(traceIntegratorB.GetSize());
          TVectorD traceIntegratorAVecAnalog(traceIntegratorAAmplifier.GetSize());
          TVectorD traceIntegratorBVecAnalog(traceIntegratorBAmplifier.GetSize());
          TVectorD analogicalTime(traceAnalogical.GetNoElements());
          TVectorD digitalTime(traceIntegratorA.GetSize());

          const double sRate = frontEnd.GetSampleTimeADC();

          for (int i=0; i<traceAnalogical.GetNoElements(); i++){
                  analogicalTime[i] = (i*binning + minTimePreFE + delay);
          double signal1 = (double)traceIntegratorAAmplifier[i];
          double signal2 = (double)traceIntegratorBAmplifier[i];
          unsigned short adcLG = frontEnd.GetADCCounts(signal1);
          unsigned short adcHG = frontEnd.GetADCCounts(signal2);
                  traceIntegratorAVecAnalog[i] = adcLG;
                  traceIntegratorBVecAnalog[i] = adcHG;
          }

          for (unsigned int i=0; i<traceIntegratorA.GetSize(); i++){
                  digitalTime[i] = (i*sRate + traceStartTime);
                  traceIntegratorAVec[i] = ((double)traceIntegratorA[i]);
                  traceIntegratorBVec[i] = ((double)traceIntegratorB[i]);
          }

          fLog("Generating plots.", 3);
          // The canvas...
          std::stringstream canvasName;
          canvasName << "c"
                     << frontEnd.GetModule().GetCounter().GetId() << "_"
                     << frontEnd.GetModule().GetId();

          TCanvas canvas(canvasName.str().c_str() ,"Integrator", 200, 10, 700, 500);
          canvas.SetTickx();
          canvas.SetTicky();
          canvas.SetBorderMode(0);
          canvas.SetFillColor(0);
          canvas.SetFrameBorderMode(0);
          //canvas.SetGridx();
          //canvas.SetGridy();

          TGraph traceIntegratorAPlot(digitalTime, traceIntegratorAVec);
          //traceIntegratorAPlot.SetLineWidth(3);
          //traceIntegratorAPlot.SetLineColor(kRed);
          traceIntegratorAPlot.SetMarkerStyle(kFullCircle);
          traceIntegratorAPlot.SetMarkerColor(kRed);
          traceIntegratorAPlot.SetMarkerSize(1);

          TGraph traceIntegratorBPlot(digitalTime, traceIntegratorBVec);
          //traceIntegratorBPlot.SetLineWidth(3);
          //traceIntegratorBPlot.SetLineColor(kOrange);
          traceIntegratorBPlot.SetMarkerStyle(kFullCircle);
          traceIntegratorBPlot.SetMarkerColor(kOrange);
          traceIntegratorBPlot.SetMarkerSize(1);
          traceIntegratorBPlot.SetTitle("");

          TGraph traceIntegratorAAnalogPlot(analogicalTime, traceIntegratorAVecAnalog);
          traceIntegratorAAnalogPlot.SetLineWidth(3);
          traceIntegratorAAnalogPlot.SetLineColor(kRed);
          //traceIntegratorAAnalogPlot.SetMarkerStyle(31);

          TGraph traceIntegratorBAnalogPlot(analogicalTime, traceIntegratorBVecAnalog);
          traceIntegratorBAnalogPlot.SetLineWidth(3);
          traceIntegratorBAnalogPlot.SetLineColor(kOrange);
          //traceIntegratorBAnalogPlot.SetMarkerStyle(31);

          traceIntegratorBPlot.GetXaxis()->SetTitle(("Time (" + fUnits.GetTimeName() + ")").c_str());
          traceIntegratorBPlot.GetYaxis()->SetTitle("Signal (ADC Counts)");
          canvas.Update();
          canvas.cd();

          traceIntegratorBPlot.Draw("AP");
          traceIntegratorAPlot.Draw("SAME P");
          traceIntegratorAAnalogPlot.Draw("SAME");
          traceIntegratorBAnalogPlot.Draw("SAME");

          traceIntegratorBPlot.GetXaxis()->SetRangeUser(minTimePreFE+delay-200,minTimePreFE+delay+200);

          TLegend leg;
          leg.SetNColumns( 2 );
          leg.SetEntrySeparation( 0.8 );
          leg.SetX1NDC( 0.10 );
          leg.SetX2NDC( 0.90 );
          leg.SetY1NDC( 1 - fStyle->GetPadTopMargin() );
          leg.SetY2NDC( 1 );
          leg.SetBorderSize( 0 );
          leg.SetTextFont( fStyle->GetTitleFont() );
          leg.SetTextSize( fStyle->GetTitleSize()+0.01 );
          leg.AddEntry(&traceIntegratorAAnalogPlot, "Analogue Low Gain", "L");
          leg.AddEntry(&traceIntegratorBAnalogPlot, "Analogue High Gain", "L");
          leg.AddEntry(&traceIntegratorAPlot, "FPGA Sample", "P");
          leg.AddEntry(&traceIntegratorBPlot, "FPGA Sample", "P");
          leg.SetFillColor(canvas.GetFillColor());
          leg.Draw();
          canvas.Update();

          if (fPlotOutcome)
            PrintPlots(canvas, fPlotFileExtensions, Tag("outcomeIntegrator", frontEnd)(fRunNumber));
}


VModule::ResultFlag
MdOptoElectronicSimulator::SimulateElectronics(
  mevt::Module& evtModule,
  const mdet::Module& module,
  const SignalsMap& sm,
  const utl::TimeStamp& eventTime)
{
  if ( !module.IsSiPM() ) {
    for (SignalsMap::const_iterator i = sm.begin(), fpe = sm.end(); i != fpe; ++i) {
           SignalsMap::key_type pixelId = i->first;
       /*
        * If not present, make it: originally a warning was emited, but because the timestart
        * is already set at this point, then that warning was emited.
        */
        const mdet::Channel& channel     = module.GetChannelFor(module.GetPixel(pixelId));
        if (! evtModule.HasChannel(channel.GetId()))
          evtModule.MakeChannel(channel.GetId());

        fLog("Simulating PMT-electronic channel:", 2)(channel.GetId())("...");
        mevt::Channel& evtChannel = evtModule.GetChannel(channel.GetId());
        if (!evtChannel.HasTrace())
          evtChannel.MakeTrace();
        utl::TraceB& trace = evtChannel.GetTrace();
        double span;
        // The last parameter defines the trace start: the interval is calculated between
        // the reference time of the event, and the registered time for the trace start.
        // Note that the particles' GetTime (which return an interval) is assumed to
        // return the interval with respect to the event time.
        const utl::TimeStamp& start = evtChannel.GetTraceStartTime();
        fLog("Trace start at", 3)(start)("...");
        ProcessPulses(channel, i->second, trace, span, start - eventTime /*, eventTime*/);
        fLog("Total analog particle pulse over-threshold time span", 3)
        (span / fUnits.GetTimeUnit())(fUnits.GetTimeName())('.');
      }//end for loop
    } else {
      /*
       * As the integrator add ups the 64 traces, It is not a good idea to have different binnings.
       * There are two ways to go. Either I first determine the binning for all traces and process
       * the counter and integrator with this binning. Or we process first the counter, and then the
       * integrator (which means re evaluating all the PE functions). I beleive the first option is less consuming.
       */
      std::vector<utl::TraceD> analogicalTraces;

      double minTimePreFE =  std::numeric_limits<double>::max();
      double maxTimePreFE = -std::numeric_limits<double>::max();


      utl::TimeStamp startTimeTrace;

      for (SignalsMap::const_iterator i = sm.begin(), fpe = sm.end(); i != fpe; ++i) {
         double minTimeSignal;
         double maxTimeSignal;

         GetPulseTimeSpan(i->second,  minTimeSignal , maxTimeSignal);
         minTimePreFE = std::min(minTimePreFE, minTimeSignal);
         maxTimePreFE = std::max(maxTimePreFE, maxTimeSignal);

      }

      for (SignalsMap::const_iterator i = sm.begin(), fpe = sm.end(); i != fpe; ++i) {
         SignalsMap::key_type pixelId = i->first;
         const mdet::ChannelSiPM& channel     = module.GetChannelSiPMFor(module.GetSiPM(pixelId));
         if (! evtModule.HasChannel(channel.GetId()))
           evtModule.MakeChannel(channel.GetId());

         fLog("Simulating SiPM-electronic channel:", 2)(channel.GetId())("...");
         mevt::Channel& evtChannel = evtModule.GetChannel(channel.GetId());
         if (!evtChannel.HasTrace())
           evtChannel.MakeTrace();

         utl::TraceB& trace = evtChannel.GetTrace();
         utl::TraceD traceAnalogical; //check if i don't need a pointer here

         double span;

         const utl::TimeStamp start = evtChannel.GetTraceStartTime();
         //All channels in a module have the same start time
         //keep one (needed to ADC startTime)
         startTimeTrace = start;

         fLog("Trace start at", 3)(start)("...");
         ProcessPulses(channel, i->second, trace, traceAnalogical, span, start - eventTime, minTimePreFE, maxTimePreFE);

         fLog("Total analog particle pulse over-threshold time span", 2)
        (span / fUnits.GetTimeUnit())(fUnits.GetTimeName())('.');

         analogicalTraces.push_back(traceAnalogical);
      }


      if (fInjectNoiseBinary)
        InjectDigitalNoise(module, evtModule);
      /*
       * The integrator analysis could be in the ProcessPulses method.
       * However, this method is already a mess, so, I am creating a new one for the sake of clarity (and mine).
       * ProcessPulse works at the level of channel, while ProcessPulsesIntegrator
       * works at the level of Module.
       * */
      if (analogicalTraces.size()>0){
        if (!evtModule.HasIntegratorATrace())
          evtModule.MakeIntegratorATrace();

        if (!evtModule.HasIntegratorBTrace())
          evtModule.MakeIntegratorBTrace();

        utl::TraceUSI& traceIntegratorA = evtModule.GetIntegratorATrace();
        utl::TraceUSI& traceIntegratorB = evtModule.GetIntegratorBTrace();

        fLog("Generating integrator traces", 2);
        ProcessPulsesIntegrator(module, analogicalTraces, traceIntegratorA, traceIntegratorB,  startTimeTrace - eventTime, minTimePreFE, maxTimePreFE);
      }
  }
  return eSuccess;
}

VModule::ResultFlag
MdOptoElectronicSimulator::Finish()
{
  return eSuccess;
}
/*
//Helper to set trace start time independently of Optical device
template<class FrontEndType>
void SetTraceStartTime( mevt::Counter::ModuleIterator mIt, const FrontEndType& frontEnd, const utl::TimeStamp& eventTime, const double & deltaLatch ) {

    const double preT1BufferLength  = frontEnd.GetPreT1BufferLength();
    const double sRate = frontEnd.GetMeanSampleRatePeriod();
    //shift the position of the T1 to the circular md-buffer (preT1BufferLength)
    utl::TimeInterval timeShiftpre (sRate*preT1BufferLength );

    for (typename FrontEndType::ChannelConstIterator ch = frontEnd.ChannelsBegin(), ec = frontEnd.ChannelsEnd(); ch != ec; ++ch) {
     if (!mIt->HasChannel(ch->GetId()))
       mIt->MakeChannel(ch->GetId());
     // Compute the trace start wrt the station trigger T1 GPS time
     const utl::TimeStamp traceStart = eventTime + utl::TimeInterval(deltaLatch) - timeShiftpre;
     mIt->GetChannel(ch->GetId()).SetTraceStartTime(traceStart);
    }// end for channel loop
    return;
}
*/
//Helper to set trace start time independently of Optical device
template<class FrontEndType>
void SetTraceStartTime( mevt::Counter::ModuleIterator mIt, const FrontEndType& frontEnd, const utl::TimeStamp& eventTime, const double & deltaLatch ) {

    //const double preT1BufferLength  = frontEnd.GetPreT1BufferLength();
    //const double sRate = frontEnd.GetMeanSampleRatePeriod();
    //shift the position of the T1 to the circular md-buffer (preT1BufferLength)
    //utl::TimeInterval timeShiftpre (sRate*preT1BufferLength );

    for (typename FrontEndType::ChannelConstIterator ch = frontEnd.ChannelsBegin(), ec = frontEnd.ChannelsEnd(); ch != ec; ++ch) {
     if (!mIt->HasChannel(ch->GetId()))
       mIt->MakeChannel(ch->GetId());
     // As simulated pulses are base on particle times, the signals in both detectors (SD and UMD) are align by definition
     const utl::TimeStamp traceStart = (eventTime + utl::TimeInterval(deltaLatch)) /*- timeShiftpre*/;
     mIt->GetChannel(ch->GetId()).SetTraceStartTime(traceStart);
    }// end for channel loop
    return;
}
int
MdOptoElectronicSimulator::GetTriggerTimeFromSD(evt::Event& event, const mdet::Counter& cmDet, mevt::MEvent::CounterIterator cIt , double& deltaLatch)
{
  if (!event.HasSEvent()) {
    ERROR("SEvent does not exist.");
    return 0;
  }
  sevt::SEvent& sEvent = event.GetSEvent();

  int partnerId = cmDet.GetAssociatedTankId();

  if (!sEvent.HasStation(partnerId)) {
    std::ostringstream message;
    message << "partner station with id " << partnerId << "was not found in SEvent\n";
    WARNING(message);
    return 0;
  }

  sevt::Station& sStation = sEvent.GetStation(partnerId);
  const sdet::Station& detStation = det::Detector::GetInstance().GetSDetector().GetStation(sStation);

  if (!sStation.HasSimData())
    return 0;

  int ret = 0;
  sevt::StationSimData& sSimData = sStation.GetSimData();
  utl::TimeStamp firstT1Sd;
  std::string algo;

  for (sevt::StationSimData::TriggerTimeIterator it = sSimData.TriggerTimesBegin(); it != sSimData.TriggerTimesEnd(); ++it) {

    utl::TimeStamp trigTime = *it;
    const sevt::StationTriggerData& trig = sSimData.GetTriggerData(trigTime);

    if (it == sSimData.TriggerTimesBegin()) {
      firstT1Sd = trigTime;
      //algo      = trig.GetAlgorithmName();
      if (trig.IsT1Threshold())
        algo = "T1Thr ";
      if (trig.IsTimeOverThreshold())
        algo = "ToT ";
      if (trig.IsTimeOverThresholdDeconvoluted())
        algo = "ToTd ";
      if (trig.IsMultiplicityOfPositiveSteps())
        algo = "MoPs ";
      if (trig.IsT2Threshold())
        algo = "T2thr ";
    }

    if (!(trig.IsT2() || trig.IsT1()))
      continue;

    ret = 1;
    if (trig.IsT2())
       ret = 2;

    if (trigTime < firstT1Sd) {
      //There is a previous trigger
      firstT1Sd = trigTime;
      //algo      = trig.GetAlgorithmName();
      if (trig.IsT1Threshold())
        algo = "T1Thr ";
      if (trig.IsTimeOverThreshold())
        algo = "ToT ";
      if (trig.IsTimeOverThresholdDeconvoluted())
        algo = "ToTd ";
      if (trig.IsMultiplicityOfPositiveSteps())
        algo = "MoPs ";
      if (trig.IsT2Threshold())
        algo = "T2thr ";

    }
    //std::cout << "Trigger time: " << trigTime                << std::endl;
    //std::cout << "AlgorithmName: "<< trig.GetAlgorithmName() << std::endl;
    //std::cout << "Algorithm: "    << trig.GetAlgorith()      << std::endl;

  } // end loop on TriggerTimeIterator

  if (!ret && !fForcedSDTrigger)
   return 0;

  //utl::TimeStamp eventTime = event.GetMEvent().GetHeader().GetTime();
  utl::TimeStamp eventTime = event.GetSEvent().GetHeader().GetTime();
  utl::TimeInterval shift;
  utl::TimeStamp TLatch;//SD GPS LatchTime

  //avoids break if forced WCD trigger (set T1 time stamp as event time + trace length)
  if ( fForcedSDTrigger )
    firstT1Sd = eventTime + detStation.GetFADCTraceLength() * detStation.GetFADCBinSize();

  shift = utl::TimeInterval((detStation.GetFADCTraceLength() - detStation.GetLatchBin())*detStation.GetFADCBinSize());

  TLatch = firstT1Sd - shift;//Absolute (GPS) time of the latch bin
  //Time interval from Eventime to TLatch
  deltaLatch = TLatch - eventTime;//Relative to T0==EventTime

  //Bins for the SD in relative to T0==EventTime (early stations < 0 and late stations > 0)
  int binTrgStop = int((firstT1Sd-eventTime)/detStation.GetFADCBinSize());
  int binTrgLatch= binTrgStop + detStation.GetLatchBin() - 1 - detStation.GetFADCTraceLength();
  int binTrgStart= binTrgStop - detStation.GetFADCTraceLength();

  utl::TabularStream tab("r|r|r|r|r|r|r|l");
  tab << "Sation" << utl::endc
      << "T1SD-EvtTime(ns)" << utl::endc
      << "StartBin" << utl::endc
      << "LatchBin" << utl::endc
      << "StopBin"  << utl::endc
      << "T1MD"     << utl::endc
      << "Type"     << utl::endc
      << "" << utl::endr;

  tab  << sStation.GetId() << utl::endc;
  tab  << firstT1Sd - eventTime   << utl::endc;
  tab  << binTrgStart << utl::endc;
  tab  << binTrgLatch << utl::endc;
  tab  << binTrgStop  << utl::endc;
  tab  << deltaLatch  << utl::endc;
  tab  << algo        << utl::endr;

  std::cout << tab << std::endl;

  for (mevt::Counter::ModuleIterator mIt = cIt->ModulesBegin(), em = cIt->ModulesEnd();
       mIt != em; ++mIt) {

    const int mId = mIt->GetId();
    const mdet::Module & mmDet = cmDet.GetModule(mId);

    if ( !mmDet.IsSiPM() ) { 
      const mdet::FrontEnd& frontEnd = mmDet.GetFrontEnd();
      //SetTraceStartTime <mdet::FrontEnd>  ( mIt, frontEnd, eventTime, deltaLatch );
      SetTraceStartTime <mdet::FrontEnd>  ( mIt, frontEnd, eventTime, (fForcedSDTrigger?0:binTrgLatch*detStation.GetFADCBinSize()) );
    }  else  {                
      const mdet::FrontEndSiPM& frontEnd = mmDet.GetFrontEndSiPM();
      //SetTraceStartTime <mdet::FrontEndSiPM>  ( mIt, frontEnd, eventTime, deltaLatch );
      SetTraceStartTime <mdet::FrontEndSiPM>  ( mIt, frontEnd, eventTime, (fForcedSDTrigger?0:binTrgLatch*detStation.GetFADCBinSize()) );
    }
  }//end loop on modules
  return ret;
}


void
MdOptoElectronicSimulator::Dump(const TF1& fun, const std::string& suffix)
{
  double min;
  double max;
  fun.GetRange(min, max);
  // Temporal step for the case when the whole pulse is sampled.
  const double stepWhole  = (max - min)        / fNPulseSamples;
  // Temporal step for the case when a fixed time window is sampled.
  const double stepWindow = fPulseSampleWindow / fNPulseSamples;
  for (unsigned int i = 0; i < fNPulseSamples; ++i) {
    double xWhole  = min + i * stepWhole;
    double xWindow = min + i * stepWindow;
    fPulseFile << fRunNumber        << '\t'
               << suffix            << '\t'
               << xWhole            << '\t'
               << fun.Eval(xWhole)  << '\t'
               << xWindow           << '\t'
               << fun.Eval(xWindow)
               << std::endl;
  }
  /*
   * Also separate each batch with a blank, which can be handy for plotting with something like
   * GNUPlot, which consideres white lines to mark the ending of data series.
   */
  fPulseFile << std::endl;
}

void
MdOptoElectronicSimulator::Init(std::unique_ptr<utl::TabularStream>& pt, unsigned int nCol)
  const
{
  std::stringstream fmt;
  for(unsigned int i = 0; i < nCol; ++i)
    fmt << "|r";
  fmt << "|";
  pt.reset(new utl::TabularStream(fmt.str()));
  /*
   * At last this call doesn't work as expected
   * because it only ends applying to the
   * current column.
   */
  fLog.ApplyConfigurationOn(*pt);
}

}