MdCounterSimulator.cc

Go to the documentation of this file.

#include "MdCounterSimulator.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() << '_'
                                << 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 MdCounterSimulatorAG {

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

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

MdCounterSimulator::MdCounterSimulator() :
  fRunNumber(1)
{
}

MdCounterSimulator::~MdCounterSimulator()
{
}

VModule::ResultFlag
MdCounterSimulator::Init()
{
  // Configuration reading.
  utl::Branch config = fwk::CentralConfig::GetInstance()->GetTopBranch("MdCounterSimulator");
  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, "simType", tag, kSimulationTypeTags[0]);
  fSimType = utl::ConsecutiveEnumFactory<SimulationType, eFromMEventSimulatedScint, kSimulationTypeTags>::Create(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, false);
  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:
  switch (fSimType) {
  case eFromMEvent:
  case eFromMEventSimulatedScint:
    fPlotChannelPulses         = true;
    fPlotOutcome               = true;
    fPlotInputOutput           = true;
    fPlotTransferAmpResponse   = true;
    fPlotTransferPhaseResponse = true;
    fPlotDelay                 = true;
    fPlotFftPriorAmp           = true;
    fPlotFftPostAmp            = true;
    break;
  }
  // Now apply the toggle:
  if (fTogglePlotchannelPulses)
    fPlotChannelPulses         = !fPlotFftPostAmp;
  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
MdCounterSimulator::Run(evt::Event& event)
{
  VModule::ResultFlag r = eFailure; // default.

  fStyle->cd();
  // While loading the type, it's checked to match one of the enums.
  switch (fSimType) {
    case eFromMEvent:
      r = RunFromMEvent(event);
      break;
    case eFromMEventSimulatedScint:
      r = RunFromMEventScintillatorSimulated(event);
      break;
  }
  ++fRunNumber;
  return r;
}

VModule::ResultFlag
MdCounterSimulator::RunFromMEvent(evt::Event& event)
{
  INFO("RunFromMEvent");
  if (event.HasMEvent()) {
    mevt::MEvent& me = event.GetMEvent();
    const utl::TimeStamp& mEventTime = me.GetHeader().GetTime();
    fLog("The event time for 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());
      for (mevt::Counter::ModuleIterator im = ic->ModulesBegin(), em = ic->ModulesEnd(); im != em; ++im) {
        const mdet::Module& module = counter.GetModule(im->GetId());
        // Index with pixel id, unique on a module level, so no need of multi indexes.
        // We keep these so as to have together all the pulses within a module
        // , as they are related due to cross talk.
        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());
            //passing a const and a nonconst version of ScintillatorSimData, should be reorganized
            // Should be nonconst if you are going the write in the object. No need to pass twice the same object (???)
            VModule::ResultFlag f = SimulatePulses(scint, is->GetSimData(), is->GetSimData(), signalsMap);
            // Check status and do early return if so...
            if (f != eSuccess)
              return f;
          }
        }
        // On a module level we simulate the electronics.
        VModule::ResultFlag f = SimulateElectronics(*im, module, signalsMap, mEventTime);

        if (f != eSuccess)
          return f;
      }
    }
  } else {
    INFO("Event without MEvent: nothing to be done.");
  }
  return eSuccess;
}


/*********************************************************************************
 *                                                                                                                                                               *
 * 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
MdCounterSimulator::RunFromMEventScintillatorSimulated(evt::Event& event)
{
  INFO("RunFromMEventScintillatorSimulated");
  if (event.HasMEvent()) {
    mevt::MEvent& me = event.GetMEvent();
    double T1 = 0;
    const utl::TimeStamp& mEventTime = event.GetSEvent().GetHeader().GetTime();
    fLog("The event time for SiPM-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, T1);
      //Check trigger condition in WCD
      if (triggerFound) {
        std::cout << "SD trigger found for counter " << ic->GetId() << " at T1 " << T1/utl::ns << "ns " << std::endl;
        std::cout << "Trigger is T" << triggerFound << "\n";
        ic->SetT1(T1);
      } else if( fForcedSDTrigger) {
        std::cout << "XXX NO SD trigger found for counter " << ic->GetId() << " but used anyway" << std::endl;
        ic->SetT1(T1);
      } 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
      }
      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();

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

                const double delay = simData.GetPhotonTime(*photon);
                                int finalPixId;
                                bool isSiPM = module.IsSiPM();
                                if(!isSiPM){
                                        const mdet::Pixel& pixel = module.GetPixelFor(scint);
                                        mdet::Pixel::SPE spe = pixel.MakeSPEAt(delay);

                                        if (fIgnoreCrossTalk) {
                                                finalPixId = pixel.GetId();
                                        } else {
                                                const mdet::Pixel& dst = pixel.GetPMT().ComputePulseDestination(pixel);
                                                finalPixId = dst.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);
                                        nroPhotons++;
                                        finalPixId = sipm.GetId();

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

                                }
            }
            //cout << "Nro Photons in pixel " <<  module.GetSiPMFor(scint).GetId() << " " << nroPhotons << endl;
          }
        }

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

VModule::ResultFlag
MdCounterSimulator::SimulatePulses(const mdet::Scintillator& scint,
                                   const mevt::ScintillatorSimData& ssd,
                                   mevt::ScintillatorSimData& ssd_nonconst,
                                   SignalsMap& sm)
{

  const mdet::Module& module = scint.GetModule();
  bool isSiPM = module.IsSiPM();
  //fLog("Simulating pulses for", 1)(scint.GetIdMessage())(".");
  // Store only if there will be any signal, at least due to this source.

  if (ssd.HasAnalogicTrace()) {
        if(!isSiPM){
                const mdet::Pixel& pixel = module.GetPixelFor(scint);
                sm[ pixel.GetId() ].fScintSimData = &ssd;
        }else{
                const mdet::SiPM& sipm = module.GetSiPMFor(scint);
                sm[ sipm.GetId() ].fScintSimData = &ssd;
        }
  }
  // To show information at the particle level.
  const unsigned int particleTabLogLevel = 3;
  std::unique_ptr<utl::TabularStream> particleTab;
  if (fLog.GetLevel() >= particleTabLogLevel) {
    Init(particleTab, 8);
    *particleTab << utl::hline
                 << "# hit"        << utl::endc // 1
                 << "id hit scint" << utl::endc // 2
                 << "spe number"   << utl::endc // 3
                 << "depth ("            << fUnits.GetLengthName() << ")" << utl::endc // 4
                 << "track ("            << fUnits.GetLengthName() << ")" << utl::endc // 5
                 << "total path ("       << fUnits.GetLengthName() << ")" << utl::endc // 6
                 << "in-scint path ("    << fUnits.GetLengthName() << ")" << utl::endc // 7
                 << "on-manifold path (" << fUnits.GetLengthName() << ")" << utl::endr // 8
                 << utl::hline;
  }
  // To show information at the spe level.
  const unsigned int speTabLogLevel = 4;
  std::unique_ptr<utl::TabularStream> speTab;
  if (fLog.GetLevel() >= speTabLogLevel) {
    Init(speTab, 9);
    *speTab << utl::hline
            << "# hit"        << utl::endc // 1
            << "# spe"        << utl::endc // 2
            << "path delay ("     << fUnits.GetTimeName() << ")" << utl::endc // 3
            << "particle delay (" << fUnits.GetTimeName() << ")" << utl::endc // 4
            << "fiber delay ("    << fUnits.GetTimeName() << ")" << utl::endc // 5
            << "scint delay  ("   << fUnits.GetTimeName() << ")" << utl::endc // 6
            << "total delay ("    << fUnits.GetTimeName() << ")" << utl::endc // 7
            << "source pix" << utl::endc                                      // 8
            << "end pix"    << utl::endr                                      // 9
            << utl::hline;
  }
  unsigned int nHit = 0;
  for (mevt::ScintillatorSimData::ConstParticleIterator i = ssd.ParticlesBegin(), sse = ssd.ParticlesEnd();
    i != sse || fLog("", 2); ++i) {
    // ...the log puts the newline.
    const utl::Particle& particle = *i;

    /*
     * Filter through the allowed particles, see the declaration for fAllowedParticles
     * for more comments.
     * If the allowed particles set is empty, then we take that as a wildcard that means
     * to allow every particle.
     */
    if (!fAllowedParticleTypes.size() || fAllowedParticleTypes.count(particle.GetType())) {
      // Construct a line assuming rectilinear propagation of the particle.
      const utl::Line particlePath(particle.GetPosition(), particle.GetDirection());
      // If put in the simulation data for this particular scintillator, the
      // particle should intersect it; despite that it doesn't hurt to check
      // that it actually does.
      if (scint.IntersectedBy(particlePath)) {
        fLog(".", 2); // hit shown with a period...
        ++nHit;       // count the hit...
        const mdet::Fiber& fiber = module.GetFiberFor(scint);

        // Construct the track inside the scintillator:
        // the origin is the incidence point.
        const utl::Segment track = scint.ComputeIntersectionWith(particlePath);
        const utl::Point hitPoint(track.GetBegin());
        const unsigned int hitLogLevel = 33;
        if (fLog.GetLevel() >= hitLogLevel) {
          fLog("Track begin in scint", hitLogLevel)
            ("x")(track.GetBegin().GetX(scint.GetLocalCoordinateSystem()) / fUnits.GetLengthUnit())
              (fUnits.GetLengthName())
            ("y")(track.GetBegin().GetY(scint.GetLocalCoordinateSystem()) / fUnits.GetLengthUnit())
              (fUnits.GetLengthName())
            ("z")(track.GetBegin().GetZ(scint.GetLocalCoordinateSystem()) / fUnits.GetLengthUnit())
              (fUnits.GetLengthName());
          fLog("Track end in scint", hitLogLevel)
            ("x")(track.GetEnd().GetX(scint.GetLocalCoordinateSystem()) / fUnits.GetLengthUnit())
              (fUnits.GetLengthName())
            ("y")(track.GetEnd().GetY(scint.GetLocalCoordinateSystem()) / fUnits.GetLengthUnit())
              (fUnits.GetLengthName())
            ("z")(track.GetEnd().GetZ(scint.GetLocalCoordinateSystem()) / fUnits.GetLengthUnit())
              (fUnits.GetLengthName());
        }
        // Distance inside scintillator.
        const double trackLen = track.GetLength();
        /*
         * Distance to PMT:
         * We take as the impinging point the middle of the track vector,
         * and then compute this longitudinal distance along the fiber/scint.
         * We only care about 1-D distance to the start of the Fiber in the scint.
         *
         * The distance along the scint is in the y direction; and in the fiber is
         * in the x direction.
         */
        //So move the initial point half-way: this our reference point.
        const utl::Point trackMiddle = track.GetBegin() + (trackLen / 2) * track.GetDirection();
        // Now compute a plane perpendicular to the fiber.
        const utl::Plane targetPlane(utl::Point(0, 0, 0, fiber.GetLocalCoordinateSystem()),
                                     utl::Vector(1, 0, 0, fiber.GetLocalCoordinateSystem()));
        // Now the on-scintillator distance is the distance of the middle point, to the plane,
        // so we are neglecting the components perpendicular to the fiber (ie we only care of the
        // distance along the fiber).
        const double onScint = utl::Distance(trackMiddle, targetPlane);
        // We can take this length from the actual detector, or use the override value.
        const double onManifold = fiber.GetOnManifoldLength();
        // The on-scint len + on-manifold until reaching the PMT.
        const double pathLen = onScint + onManifold;
        /*
         * This particle generated this number of SPEs:
         * Note that all the energy is taken as parameter for the computation of the n spe, which
         * ends up meaning that all the energy is deposited in the scintillator (remember that
         * the fiber is used for the computation because it's the one related to the meaning of
         * attenuation (in the sense of propagation along a direction), tough also the scintillator
         * is relevant.
         */
        //#warning "ENERGY DEPOSITED BY THE PARTICLE IS NOW PASSED in the WEIGHT field (a la MARTA :-))" NO!
        //--> This attribute must be changed (or create a new one). There is no need whatsoever to do things like this, not funny.
        //General comment: when you get the option to do things right or wrong, always choose right.
        unsigned int nSPE = fiber.ComputeSPENumber(pathLen, trackLen, particle.GetWeight());

        if (fLog.GetLevel() >= particleTabLogLevel) {
          // Compute the depth where the particle hit the scintillator:
          const double hitDepth = module.GetDepth(hitPoint);
          *particleTab << nHit          << utl::endc
                       << scint.GetId() << utl::endc
                       << nSPE          << utl::endc
                       << fLog << hitDepth   / fUnits.GetLengthUnit() << utl::endc
                       << fLog << trackLen   / fUnits.GetLengthUnit() << utl::endc
                       << fLog << pathLen    / fUnits.GetLengthUnit() << utl::endc
                       << fLog << onScint    / fUnits.GetLengthUnit() << utl::endc
                       << fLog << onManifold / fUnits.GetLengthUnit() << utl::endr;
        }
        // These are the same for every spe:
        const double pathDelay  = pathLen * fiber.GetRefractionIndex() / utl::kSpeedOfLight;
        // Now the arrival time  (see also RunFromMEvent).
        const double particleDelay = particle.GetTime().GetInterval();

        for (unsigned int n = 0; n < nSPE; ++n) {
          // Geometrical propagation...
          double delay = pathDelay;
          // In-fiber decay delay...
          const double fiberDelay = fiber.ComputeDecayDelay();
          delay += fiberDelay;
          //  and in-scint decay delay...
          const double scintDelay = scint.ComputeDecayDelay();
          delay += scintDelay;
          // and the particle arrival time.
          delay += particleDelay;

          int finalPixId;
          if(!isSiPM){
                          const mdet::Pixel& pixel = module.GetPixelFor(scint);
                          mdet::Pixel::SPE spe = pixel.MakeSPEAt(delay);

                          /*
                           * Compute the final destination taking into account (eventually)
                           * the crosstalk from the PMT.
                           */
                          if (fIgnoreCrossTalk) {
                                finalPixId = pixel.GetId();
                          } else {
                                // Now determine the actual final pix, due to crosstalk.
                                const mdet::Pixel& dst = pixel.GetPMT().ComputePulseDestination(pixel);
                                finalPixId = dst.GetId();
                          }
                          //signal map is only used inside this Module
                          sm[ finalPixId ].fPulses.push_back(spe);
                          //save pulse information also in MEvent, duplicating of fPulses structure in MEvent and locally should be removed in the future..
                          ssd_nonconst.AddSPEPulse(spe.GetMu(), spe.GetSigma(), spe.GetAmplitude(), finalPixId);

          }else{
                          const mdet::SiPM& sipm = module.GetSiPMFor(scint);
                          mdet::SiPM::PE pe = sipm.MakePEAt(delay);
                          finalPixId = sipm.GetId();

                          //signal map is only used inside this Module
                          sm[ finalPixId ].fSiPMPulses.push_back(pe);
                          //save pulse information also in MEvent, duplicating of fPulses structure in MEvent and locally should be removed in the future..
                          ssd_nonconst.AddPEPulse(pe.GetStartTime(), pe.GetAmplitude1(), pe.GetAmplitude2(), pe.GetAmplitude2(),
                                          pe.GetRiseTime(), pe.GetFallTime1(), pe.GetFallTime2(), pe.GetFallTime3(), finalPixId);


          }

          if (fLog.GetLevel() >= speTabLogLevel){
            *speTab << nHit << utl::endc
                    << fLog << (n+1)                                << utl::endc
                    << fLog << pathDelay     / fUnits.GetTimeUnit() << utl::endc
                    << fLog << particleDelay / fUnits.GetTimeUnit() << utl::endc
                    << fLog << fiberDelay    / fUnits.GetTimeUnit() << utl::endc
                    << fLog << scintDelay    / fUnits.GetTimeUnit() << utl::endc
                    << fLog << delay         / fUnits.GetTimeUnit() << utl::endc
                   // << fLog << sipm.GetId()                        << utl::endc
                    << finalPixId                                   << utl::endr;
          }
        } // SPE loop
        // A line for each group.
        if (fLog.GetLevel() >= speTabLogLevel)
          (*speTab) << utl::hline;
      } else {
        // Intersection check: Warn if no intersection.
        fLog("Particle missing the scintillator after hit", 1)(nHit);
      }
    } // Particle type check
  }

  // Information for the particles...
  if (fLog.GetLevel() >= particleTabLogLevel) {
    *particleTab << utl::hline;
    fLog("Impinging particles information:", particleTabLogLevel);
    fLog(particleTab->Str(), particleTabLogLevel);
  }
  // Information for the particles...
  if (fLog.GetLevel() >= speTabLogLevel){
    // No need as of the hline from each group.
    // *speTab << utl::hline;
    fLog("Generated spe information:", speTabLogLevel);
    fLog(speTab->Str(), speTabLogLevel);
  }
  return eSuccess;
}

void
MdCounterSimulator::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.ClearTrace();
  //
  // 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);
  bool isAnalogicTrace = si.fScintSimData && si.fScintSimData->HasAnalogicTrace();
  if (si.fPulses.empty() && ! isAnalogicTrace) {
    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;
    }
  }
  // May not be used if !isAnalogicTrace, put it here to have greater scope.
  double analogicTraceStartTime = 0;
  double analogicTraceEndTime   = 0;
  if (isAnalogicTrace) {
    analogicTraceStartTime = (si.fScintSimData->GetAnalogicTraceStartTime() - eventTime).GetInterval();
    analogicTraceEndTime   = analogicTraceStartTime + si.fScintSimData->GetAnalogicTrace().GetSize() *
                                                si.fScintSimData->GetAnalogicTrace().GetBinning();
    minTimePreFE = std::min(minTimePreFE, analogicTraceStartTime);
    maxTimePreFE = std::max(maxTimePreFE, analogicTraceEndTime);
  }
  // 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);
    if (isAnalogicTrace) {
      // Interpolation logic.
      if (analogicTraceStartTime <= time && time <= analogicTraceEndTime) {
        // TODO More complex interpolation.
        double delta = time - analogicTraceStartTime;
        delta /= si.fScintSimData->GetAnalogicTrace().GetBinning();
        utl::TraceD::SizeType index = (utl::TraceD::SizeType) delta;
        value += si.fScintSimData->GetAnalogicTrace()[index];
      }
    }
    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;
  // sampleTimes.reserve(nBinsTrace);
  fLog("Sampling pulses.", 3);
  const double startTime  = traceStart.GetInterval();
  const double stopTime   = startTime +
    sRate * (frontEnd.GetPreT1BufferLength() + 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.");
  }
  fLog("Samples:", 4);
  // 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.
  for (double sampleTime = startTime; trace.GetSize() < nBinsTrace || fLog('\n', 4)(sampleTime);
    sampleTime += sRate) {
    //
    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())('|');
    }
    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(kFullSquare);
  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 );
  leg.SetX2NDC( 1 );
  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, "Analogical pulse","L");
  leg.AddEntry(&discrimFun, "Discriminator pulse", "L");
  leg.AddEntry(&samples, "Digital samples", "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.ClearTrace();
 *
 * 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
MdCounterSimulator::ProcessPulses(const mdet::ChannelSiPM& channel,
                                  const SignalInformation& si,
                                  utl::TraceB& trace,
                                                                  utl::TraceD& traceAnalogical,
                                  double& span,
                                  const utl::TimeInterval& traceStart,
                                                                  double& minTimePreFE,
                                                                  double& maxTimePreFE,
                                                                  double& analogicTraceStartTime,
                                                                  double& analogicTraceEndTime)
{

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

  const double sRate = frontEnd.GetMeanSampleRatePeriod();
  trace.ClearTrace();
  traceAnalogical.ClearTrace();
  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);

  bool isAnalogicTrace = si.fScintSimData && si.fScintSimData->HasAnalogicTrace();
  if (si.fSiPMPulses.empty() && ! isAnalogicTrace) {
    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, analogicTraceStartTime, analogicTraceEndTime,
                  totalPulsePostFrontEndTrace, totalPulsePostDiscriminatorTrace, traceAnalogical);


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

  PlotChannel(traceStart.GetInterval(), minTimePreFE, pulseTimeSpan / fNDiscretization, channel, traceAnalogical,
                  totalPulsePostFrontEndTrace, totalPulsePostDiscriminatorTrace, trace);

}

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

double
MdCounterSimulator::GetPulseTimeSpan(const SignalInformation& si, const utl::TimeStamp& eventTime, double& minTimePreFE, double& maxTimePreFE,
                double& analogicTraceStartTime, double& analogicTraceEndTime){

        bool isAnalogicTrace = si.fScintSimData && si.fScintSimData->HasAnalogicTrace();

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

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

        if (isAnalogicTrace) {
          analogicTraceStartTime = (si.fScintSimData->GetAnalogicTraceStartTime() - eventTime).GetInterval();
          analogicTraceEndTime   = analogicTraceStartTime + si.fScintSimData->GetAnalogicTrace().GetSize() *
                                                                                                  si.fScintSimData->GetAnalogicTrace().GetBinning();
          minTimePreFE = std::min(minTimePreFE, analogicTraceStartTime);
          maxTimePreFE = std::max(maxTimePreFE, analogicTraceEndTime);
        }
        // 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
MdCounterSimulator::ApplyCITIROCTransfer(const mdet::ChannelSiPM& channel, const SignalInformation& si, const double pulseTimeSpan, double minTimePreFE,
                double analogicTraceStartTime, double analogicTraceEndTime, TimeTrace& totalPulsePostFrontEndTrace,
                TimeTrace& totalPulsePostDiscriminatorTrace, utl::TraceD& traceAnalogical){

        TotalPulseFunctor<SiPMPulseContainer> totalPulsePreFrontEnd(si.fSiPMPulses,
                                                                  fUnits.GetTimeUnit(),
                                                                  fUnits.GetElectricPotentialUnit());
        bool isAnalogicTrace = si.fScintSimData && si.fScintSimData->HasAnalogicTrace();

        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);
          if (isAnalogicTrace) {

                if (analogicTraceStartTime <= time && time <= analogicTraceEndTime) {

                  double delta = time - analogicTraceStartTime;
                  delta /= si.fScintSimData->GetAnalogicTrace().GetBinning();
                  utl::TraceD::SizeType index = (utl::TraceD::SizeType) delta;
                  value += si.fScintSimData->GetAnalogicTrace()[index];
                }
          }
          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
MdCounterSimulator::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();
          const double stopTime   = startTime +
            sRate * (frontEnd.GetPreT1BufferLength() + frontEnd.GetPostT1BufferLength());

          if (startTime < minTimePostFE) {
            fLog("Samples false in the range", 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);
            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 ", 4)
                (maxTimePostFE / fUnits.GetTimeUnit())(fUnits.GetTimeName())
                ("to")
                (stopTime / fUnits.GetTimeUnit())(fUnits.GetTimeName())
                ("as no signal in it.");
          }

}

void
MdCounterSimulator::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
MdCounterSimulator::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();

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

          for(unsigned int i=0; i<trace.GetSize(); i++){
                  digitalTime[i] = (i*sRate + traceStartTime);
                  traceFPGA[i] = ((double)trace[i]*discriminatorHiLevel);//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(29);
          pulsePostFPGAPlot.SetMarkerSize(1.2);

          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");

          TLegend leg;
          leg.SetNColumns( 3 );
          leg.SetEntrySeparation( 0.8 );
          //leg.SetX1NDC( fStyle->GetPadLeftMargin() );
          //leg.SetX2NDC( 1 - fStyle->GetPadRightMargin() );
          leg.SetX1NDC( 0 );
          leg.SetX2NDC( 1 );
          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 pulse","L");
          leg.AddEntry(&pulsePostFrontEndPlot, "Fast Shaper output", "L");
          leg.AddEntry(&pulsePostDiscriminatorPlot, "Discriminator output", "L");
          leg.AddEntry(&pulsePostFPGAPlot, "FPGA samples", "P");
          leg.SetFillColor(canvas.GetFillColor());
          leg.Draw();
          canvas.Update();

         /* std::stringstream outputName;
          outputName << "./output/output_"
                            << channel.GetFrontEnd().GetModule().GetCounter().GetId() << "_"
                            << channel.GetFrontEnd().GetModule().GetId() << "_"
                            << channel.GetId()
                                                << ".root";


          canvas.SaveAs(outputName.str().c_str());*/
          if (fPlotOutcome)
            PrintPlots(canvas, fPlotFileExtensions, Tag("outcome", channel)(fRunNumber));

}

/*
 *
 *Simulate the Integrator electronics. 1 run per module
 *
 * */
void
MdCounterSimulator::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.ClearTrace();
  traceIntegratorB.ClearTrace();
  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;
        }
  }

//  const double pulseTimeSpan = maxTimePreFE-minTimePreFE;


  /**********************************************
   *   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);
  PlotIntegrator(traceStart.GetInterval(), minTimePreFE, (maxTimePreFE-minTimePreFE) / fNDiscretization,
                  frontEnd, totalAnalogicalInput, traceAfterADCLowGain, traceAfterADCHighGain, traceIntegratorA, traceIntegratorB);
  //PlotIntegrator(traceStart.GetInterval(), minTimePreFE, (maxTimePreFE-minTimePreFE) / fNDiscretization, frontEnd, totalAnalogicalInput, traceAfterADCLowGain, traceAfterADCHighGain);
}


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

void
MdCounterSimulator::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               *
         *******************************************************/

        //ApplyTransferBlock(fft, backEnd, backEnd.eLowGainAmplifier);
        //ApplyTransferBlock(fftHGChannel, backEnd, backEnd.eHighGainAmplifier);

        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
MdCounterSimulator::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
MdCounterSimulator::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
MdCounterSimulator::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
MdCounterSimulator::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("Sampling pulses.", 3);
          const double startTime  = traceStart.GetInterval();
          const double stopTime   = startTime + sRate * nBinsTrace;

          if (startTime + delayBinaryADC < minTimePostFE) {
            fLog("Samples false in the range", 4)
                (startTime / fUnits.GetTimeUnit())(fUnits.GetTimeName())
                ("to")
                (minTimePostFE / fUnits.GetTimeUnit())(fUnits.GetTimeName())
                ("as no signal in it.");
          }
          fLog("Samples:", 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)('@')(sampleTime / fUnits.GetTimeUnit())(fUnits.GetTimeName())('|');
              fLog(adcHG, 4, false)('@')(sampleTime / fUnits.GetTimeUnit())(fUnits.GetTimeName())('|');
            }
            if(fIncludeBaseLineFluctuationIntegrator){
                adcLG += frontEnd.GetADCBaseLineFluctuationLG();
                adcHG += frontEnd.GetADCBaseLineFluctuationHG();
            }
            traceIntegratorA.PushBack(adcLG);
            traceIntegratorB.PushBack(adcHG);
          }
          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.");
          }

}

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

          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);
                  traceIntegratorAVecAnalog[i] = (double)traceIntegratorAAmplifier[i];
                  traceIntegratorBVecAnalog[i] = (double)traceIntegratorBAmplifier[i];
          }

          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 inputSiPM(analogicalTime, traceAnalogical);
          inputSiPM.SetLineWidth(3);
          inputSiPM.SetLineColor(kBlue);

          TGraph traceIntegratorAPlot(digitalTime, traceIntegratorAVec);
          traceIntegratorAPlot.SetLineWidth(3);
          traceIntegratorAPlot.SetLineColor(kRed);
          traceIntegratorAPlot.SetMarkerStyle(31);

          TGraph traceIntegratorBPlot(digitalTime, traceIntegratorBVec);
          traceIntegratorBPlot.SetLineWidth(3);
          traceIntegratorBPlot.SetLineColor(kOrange);
          traceIntegratorBPlot.SetMarkerStyle(31);

          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("APL");
          traceIntegratorAPlot.Draw("SAMEPL");
          inputSiPM.Draw("SAME");
          traceIntegratorAAnalogPlot.Draw("SAME");
          traceIntegratorBAnalogPlot.Draw("SAME");

          TLegend leg;
          leg.SetNColumns( 3 );
          leg.SetEntrySeparation( 0.8 );
          //leg.SetX1NDC( fStyle->GetPadLeftMargin() );
          //leg.SetX2NDC( 1 - fStyle->GetPadRightMargin() );
          leg.SetX1NDC( 0 );
          leg.SetX2NDC( 1 );
          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(&inputSiPM, "SiPM input","L");
          leg.AddEntry(&traceIntegratorAPlot, "ADC LG", "L");
          leg.AddEntry(&traceIntegratorBPlot, "ADC HG", "L");
          leg.SetFillColor(canvas.GetFillColor());
          leg.Draw();
          canvas.Update();

         /* std::stringstream outputName;
          outputName << "./output/Integretor_"
                            << frontEnd.GetModule().GetCounter().GetId() << "_"
                            << frontEnd.GetModule().GetId() << ".root";

          canvas.SaveAs(outputName.str().c_str());*/
          if (fPlotOutcome)
            PrintPlots(canvas, fPlotFileExtensions, Tag("outcomeIntegrator", frontEnd)(fRunNumber));

}



VModule::ResultFlag
MdCounterSimulator::SimulateElectronics(
  mevt::Module& evtModule,
  const mdet::Module& module,
  const SignalsMap& sm,
  const utl::TimeStamp& eventTime)
{

  bool isSiPM = module.IsSiPM();
  if(!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())('.');
          }

    } 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;
      double maxTimePreFE;
      double analogicTraceStartTime;
      double analogicTraceEndTime;
      utl::TimeStamp startTimeTrace;

          for (SignalsMap::const_iterator i = sm.begin(), fpe = sm.end(); i != fpe; ++i) {
              double minTimeSignal;
              double maxTimeSignal;
              double analogicTraceStartTimePEs;
              double analogicTraceEndTimePEs;
              GetPulseTimeSpan(i->second, eventTime, minTimeSignal , maxTimeSignal , analogicTraceStartTimePEs , analogicTraceEndTimePEs);
                  minTimePreFE = std::min(minTimePreFE, minTimeSignal);
                  maxTimePreFE = std::max(maxTimePreFE, maxTimeSignal);
                  analogicTraceStartTime = std::min(analogicTraceStartTime, analogicTraceStartTimePEs);
                  analogicTraceEndTime = std::max(analogicTraceEndTime, analogicTraceEndTimePEs);
          }

          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();
                startTimeTrace = start;
                fLog("Trace start at", 3)(start)("...");
                ProcessPulses(channel, i->second, trace, traceAnalogical, span, start - eventTime, minTimePreFE, maxTimePreFE, analogicTraceStartTime, analogicTraceEndTime);

                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
MdCounterSimulator::Finish()
{
  return eSuccess;
}

int
MdCounterSimulator::GetTriggerTimeFromSD(evt::Event& event, const mdet::Counter& cmDet,mevt::MEvent::CounterIterator cIt , double& T1)
{
  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;

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

    if (it == sSimData.TriggerTimesBegin())
      firstT1Sd = *it;

    utl::TimeStamp trigTime = *it;
    const sevt::StationTriggerData& trig = sSimData.GetTriggerData(trigTime);
    std::cout << "IsT1 " << trig.IsT1() << " - " << "IsT2 " << trig.IsT2() << std::endl;

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

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

    if (trigTime < firstT1Sd) {
      std::cout << "There is a previous trigger\n";
      firstT1Sd = trigTime;
    }
    std::cout << "Trigger time: " << trigTime << std::endl;
    std::cout << "Algorithm: "    << trig.GetAlgorithmName() << 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 TstartTrace;

  //avoids break if forced WCD trigger (set T1 time stamp as event time)
  if ( fForcedSDTrigger )
    firstT1Sd = eventTime;
   else
    //see sdet::StationTriggerAlgorithm.h and TankTriggerSimulator.cc
    shift = utl::TimeInterval((detStation.GetFADCTraceLength() - detStation.GetLatchBin())*detStation.GetFADCBinSize());

  TstartTrace = firstT1Sd - shift;
  //Time interval from Eventime to TstartTrace
  T1 = TstartTrace - eventTime;
  std::cout << "eventTime "   << eventTime   << std::endl;
  std::cout << "T1 "          << T1          << std::endl;
  std::cout << "firstT1Sd "   << firstT1Sd   << std::endl;
  std::cout << "TstartTrace " << TstartTrace << std::endl;
  std::cout << "shift "       << shift       << 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 SiPM prepares to iterate over module channels (mdet)
    typedef mdet::FrontEndSiPM::ChannelConstIterator chIt;
    const mdet::FrontEndSiPM& frontEnd = mmDet.GetFrontEndSiPM();
    const double sRate = frontEnd.GetMeanSampleRatePeriod();
    const double preT1BufferLength  = frontEnd.GetPreT1BufferLength();
    const double postT1BufferLength = frontEnd.GetPostT1BufferLength();

    //shift the position of the T1 to the circular md-buffer (preT1BufferLength)
    utl::TimeInterval timeShiftpre (sRate*preT1BufferLength );
    utl::TimeInterval timeShiftpost(sRate*postT1BufferLength);

    for (chIt 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(T1) - timeShiftpre;
     mIt->GetChannel(ch->GetId()).SetTraceStartTime(traceStart);
   } // end for channel loop
  }//end loop on modules
  return ret;
}

void
MdCounterSimulator::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
MdCounterSimulator::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);
}

}