FdElectronicsSimulator.cc

Go to the documentation of this file.

#include <iostream>
#include <iomanip>
#include <string>

#include <utl/Reader.h>
#include <utl/ErrorLogger.h>
#include <utl/AugerUnits.h>
#include <utl/MathConstants.h>
#include <utl/Math.h>
#include <utl/PhysicalConstants.h>
#include <utl/RandomEngine.h>
#include <utl/TabulatedFunction.h>
#include <utl/TabulatedFunctionErrors.h>
#include <utl/Trace.h>
#include <utl/TraceAlgorithm.h>
#include <utl/MultiTabulatedFunction.h>
#include <utl/AugerException.h>
#include <utl/Md5Signature.h>

#include <fwk/CentralConfig.h>
#include <fwk/RandomEngineRegistry.h>

#include <det/Detector.h>
#include <fdet/FDetector.h>
#include <fdet/Eye.h>
#include <fdet/Telescope.h>
#include <fdet/Channel.h>
#include <fdet/Pixel.h>
#include <fdet/Mirror.h>
#include <fdet/Filter.h>
#include <fdet/Camera.h>

#include <evt/Event.h>

#include <fevt/FEvent.h>
#include <fevt/Eye.h>
#include <fevt/Telescope.h>
#include <fevt/TelescopeSimData.h>
#include <fevt/PixelSimData.h>
#include <fevt/Pixel.h>
#include <fevt/ChannelSimData.h>
#include <fevt/Channel.h>

#include <CLHEP/Random/RandPoisson.h>
#include <CLHEP/Random/RandGauss.h>

#include "FdElectronicsSimulator.h"

using namespace FdElectronicsSimulatorOG;
using namespace std;
using namespace utl;
using namespace evt;
using namespace fwk;
using namespace det;
using namespace fevt;
using namespace fdet;


using CLHEP::RandGauss;
using CLHEP::RandPoisson;


FdElectronicsSimulator::FdElectronicsSimulator() :
  fVerbosityLevel(0),
  fBaseline(0),
  fUseMonitoring(false),
  fThresholdBoxCarReducedInMon(true),
  fUseNewThreshold(false),
  fThresholdMode(false),
  fRandomEngine(0)
{
}

FdElectronicsSimulator::~FdElectronicsSimulator()
{
}

VModule::ResultFlag
FdElectronicsSimulator::Init()
{
  CentralConfig* cc = CentralConfig::GetInstance();
  Branch topB = cc->GetTopBranch("FdElectronicsSimulator");

  fRandomEngine =
    &RandomEngineRegistry::GetInstance().Get(RandomEngineRegistry::eDetector);

  fVerbosityLevel = 0;
  if (topB.GetChild("verbosityLevel"))
    topB.GetChild("verbosityLevel").GetData(fVerbosityLevel);

  fPeFluctuations = true;
  if (topB.GetChild("NoPeFluctuations")) {
    topB.GetChild("NoPeFluctuations").GetData(fPeFluctuations);
    fPeFluctuations = !fPeFluctuations;
  }

  fDoAntiAliasingFilter = true;
  if (topB.GetChild("NoAntiAliasingFilter")) {
    topB.GetChild("NoAntiAliasingFilter").GetData(fDoAntiAliasingFilter);
    fDoAntiAliasingFilter = !fDoAntiAliasingFilter;
  }

  if (topB.GetChild("thresholdMode")) {
    topB.GetChild("thresholdMode").GetData(fThresholdMode);
  }

  if (topB.GetChild("useNewThreshold"))
    topB.GetChild("useNewThreshold").GetData(fUseNewThreshold);

  if (fUseNewThreshold){
    Branch topBThreshold = cc->GetTopBranch("FSimulationThreshold");
    Branch it = topBThreshold.GetFirstChild();
    while (it){
      TabulatedFunction thresholdValues;
      it.GetData(thresholdValues);
      AttributeMap thresholdSignature = it.GetAttributes();
      string configSignature = thresholdSignature.begin()->second;
      fSimThresholdValues.insert(std::pair<string, TabulatedFunction>(configSignature, thresholdValues));
      ++it;
    }
  }

  topB.GetChild("baseline").GetData(fBaseline);
  topB.GetChild("useMonitoringInfo").GetData(fUseMonitoring);
  fThresholdBoxCarReducedInMon = true;
  if (topB.GetChild("thresholdBoxCarReducedInMon")) {
    topB.GetChild("thresholdBoxCarReducedInMon").GetData(fThresholdBoxCarReducedInMon);
  }

  // info output
  ostringstream info;
  info << " Version: "
       << GetVersionInfo(VModule::eRevisionNumber) << "\n"
          " Parameters:\n"
          "             baseline: " << fBaseline << "\n";
  if (fUseMonitoring)
    info << "  use monitoring info: " << fUseMonitoring << '\n';
  if (!fPeFluctuations)
    info << "      PE fluctuations: OFF\n";
  if (!fDoAntiAliasingFilter)
    info << "  Anti-aliasing filer: OFF\n";
  if (fVerbosityLevel)
    info << "      verbosity level: " << fVerbosityLevel << '\n';

  INFO(info);
  return eSuccess;
}// end of Init


VModule::ResultFlag
FdElectronicsSimulator::Run(evt::Event& event)
{
  if (!event.HasFEvent()) {
    ERROR ("Event has no FEvent.");
    return eFailure;
  }

  fDrumCalibMode = !event.HasSimShower();

  if (fDrumCalibMode) {
    INFO ("Drum calibration mode!");
  }

  FEvent& fEvent = event.GetFEvent();
  for (fevt::FEvent::EyeIterator iEye =
         fEvent.EyesBegin(fevt::ComponentSelector::eInDAQ);
       iEye != fEvent.EyesEnd(fevt::ComponentSelector::eInDAQ);
       ++iEye) {

    if (iEye->GetStatus()==fevt::ComponentSelector::eDeSelected)
      continue;

    fevt::Eye& eyeEvent = *iEye;
    bool eyeHasAtLeastOneGoodTelescope = false;

    if (fVerbosityLevel>2) {
      ostringstream info;
      info << " Eye=" << eyeEvent.GetId();
      INFO(info);
    }
    
    for (fevt::Eye::TelescopeIterator iTel =
           eyeEvent.TelescopesBegin(fevt::ComponentSelector::eInDAQ);
         iTel != eyeEvent.TelescopesEnd(fevt::ComponentSelector::eInDAQ);
         ++iTel) {

      fevt::Telescope& telEvent = *iTel;

      if (fVerbosityLevel>2) {
        ostringstream info;
        info << " Eye=" << eyeEvent.GetId() << " Tel=" << telEvent.GetId();
        INFO(info);
      }

      for (fevt::Telescope::PixelIterator iPixel =
             telEvent.PixelsBegin(fevt::ComponentSelector::eInDAQ);
           iPixel != telEvent.PixelsEnd(fevt::ComponentSelector::eInDAQ);
           ++iPixel) {

        if (iPixel->GetStatus() == fevt::ComponentSelector::eDeSelected)
          continue;

        if (fVerbosityLevel>2) {
          ostringstream info;
          info << " Eye=" << eyeEvent.GetId()
               << " Tel=" << telEvent.GetId()
               << " Pix=" << iPixel->GetId();
          INFO(info);
        }

        if (iPixel->HasSimData()) {

          // Simulate electronics only for the cameras that contain photons
          if (InitCamera(telEvent)) {
            ElecSim(telEvent);
            eyeHasAtLeastOneGoodTelescope = true;
          } else
            telEvent.SetStatus(fevt::ComponentSelector::eDeSelected);

          break;
        }

      } // end loop over Pixels

    } // end loop over Telescopes

    if (!eyeHasAtLeastOneGoodTelescope)
      eyeEvent.SetStatus(fevt::ComponentSelector::eDeSelected);

  } // end loop over Eyes
  return eSuccess;
} // end of Run


VModule::ResultFlag
FdElectronicsSimulator::Finish()
{
  return eSuccess;
}// end of Finish


bool
FdElectronicsSimulator::InitCamera(fevt::Telescope& tel)
{
  const fdet::FDetector& detFD = Detector::GetInstance().GetFDetector();
  const fdet::Telescope& detTel = detFD.GetTelescope(tel);
  const fdet::Camera& detCamera = detTel.GetCamera();

  const bool hasCorrector = detTel.HasCorrectorRing();
  bool telHasAtLeastOneValidPixel = false;

  const double normWavelength = detFD.GetReferenceLambda();

  const unsigned int nPixels = detTel.GetLastPixelId();
  fPeBgNoiseTable.clear();
  fCalibCorrection.clear();
  fQEffAtNormWavelength.clear();
  
  // first loop to finalize the config checksum
  map<int, double> gainOfPixel;
  for (unsigned int pixelId=1; pixelId<=nPixels; pixelId++) {

    if (!tel.HasPixel(pixelId))
      continue;

    fevt::Pixel& pix = tel.GetPixel(pixelId);
    const fdet::Pixel& detPixel = detFD.GetPixel(pix);
    const fdet::Channel& detChannel = detFD.GetChannel(pix); // mapped channel
    const TabulatedFunction& qEff = detPixel.GetQEfficiency();
    const double absGain = detChannel.GetElectronicsGain();
    
    fQEffAtNormWavelength[pixelId] = qEff.Y(normWavelength);
    gainOfPixel[pixelId] = absGain;
  }

  PrepareTimeConvolution(detTel);

  //GetThresholdConfigSignature
  const string thresholdConfigSignature = detTel.GetCamera().GetThresholdConfigSignature();

  // Get the calibration constants
  for (unsigned int pixelId=1; pixelId<=nPixels; pixelId++) {

    if (!tel.HasPixel(pixelId))
      continue;

    fevt::Pixel& pix = tel.GetPixel(pixelId);
    if (pix.GetStatus() == fevt::ComponentSelector::eDeSelected)
      continue;

    if (!pix.HasSimData())
      continue;

    // Get the calibration constant for this channel
    const fdet::Pixel& detPixel = detFD.GetPixel(pix);
    const fdet::Channel& detChannel = detFD.GetChannel(pix); // mapped channel

    // get calibration constant
    const TabulatedFunction& calibration = detPixel.GetEndToEndCalibrationConstant();
    if (!calibration.GetNPoints()) {
      ostringstream errorMsg;
      errorMsg<< "Could not find calibration for eye=" << pix.GetEyeId()
              << " tel=" << pix.GetTelescopeId() << " pix=" << pixelId;
      ERROR(errorMsg);
      continue;
    }
    
    // check calibration status
    if ( detPixel.GetStatus() != fdet::Pixel::eGood ) {
      ostringstream warn;
      warn << " bad calibration data for eye " << pix.GetEyeId()
           << " tel " << pix.GetTelescopeId() << " pix " << pixelId
           << " status = " << detPixel.GetStatus();
      warn << " setting fevt::PixelStatus to DeSelected..."<<endl;
      WARNING(warn.str());
      pix.SetStatus(fevt::ComponentSelector::eDeSelected);
      continue;
    }
    telHasAtLeastOneValidPixel = true;

    const double caliconst =
      detPixel.GetEndToEndCalibrationConstantAtReferenceWavelength();
    const double opticalEfficiencyCorrection =
      detPixel.GetOpticalEfficiencyCorrectionAtReferenceWavelength();

    if (caliconst==0) {
        ostringstream err;
        err << "\n     **************************************************************************\n"
            << "     **************************************************************************\n"
            << "      ERROR caliconst=0 of pixelId=" << pixelId
            << " ! Cannot use telId=" << pix.GetTelescopeId() << " of eyeId=" << pix.GetEyeId() << '\n'
            << "     **************************************************************************\n"
            << "     **************************************************************************";
        ERROR(err);
        return false;
    }
    
    // Pixel constants into arrays
    fevt::PixelSimData& pxsimdata = pix.GetSimData();
    const double backphotons =
      pxsimdata.HasPhotonTrace(fevt::FdConstants::eBackground) ?
      pxsimdata.GetPhotonTrace(fevt::FdConstants::eBackground)[0] :
      pxsimdata.GetMeanBgPhotonFlux();
    
    const double absGain = gainOfPixel[pixelId];
    
    if (fVerbosityLevel>2) {
      cout << " drum=" << fDrumCalibMode //<< " (" <<!event.HasSimShower()<<")"
           << " PeToADC=" << absGain
           << " hasCorrector=" << hasCorrector
           << " caliconst=" << caliconst
           << endl;
    }

    //correction for calib const in reconstruction
    const string& configSignature = tel.GetSimData().GetConfigSignature();
    
    // The simulation of threshold values is done in ADC therefore no scaling to calib constants
    // In all other cases: adapt sim gain to data
    const double calibCorrection = fThresholdMode ? 1.0 : detPixel.GetSimulatedEndToEndCalibration(configSignature) / caliconst;

    // npe is calculated with default gain in FdBackgroundSimulator
    // npe is not affected by optical efficiency loss correction
    double meanBgPe = backphotons * fQEffAtNormWavelength[pixelId];

    const double gainVariance = detChannel.GetGainVariance();
    const double elecNoiseVariance = detChannel.GetElectronicNoiseVariance();
    
    // correct Gaussian approximation
    if (!fThresholdMode)
      meanBgPe *= fBandWidthFactor;    
    
    // from FdBackgroundSimulator we have: bgPhotons -> (ideal det simulation) -> ADC counts
    // but in ElecSim we do: photons -> (ideal det simulation) -> (calibCorrection) -> ADC counts
    // so we should apply inverse calibCorrection to bgPhotons to be consistent with ElecSim
    // since this is used only for Poisson variance calculation, we need to correct in quadrature 

    if (!detChannel.IsVirtual()) {
      meanBgPe /= calibCorrection * calibCorrection;
    }

    fPeBgNoiseTable[pixelId] = meanBgPe;


    // Set pixel threshold & n. samples (trigger)

    const unsigned int fltBoxCarSize = detCamera.GetFLTBoxcarSumLength();

    double threshold = 0;
    double ADCVar = 0;
    double baseline = 0;

    if (fUseMonitoring) {
      baseline = detChannel.GetBaseline();
      if (fThresholdBoxCarReducedInMon) {
        threshold = (int) detChannel.GetThreshold()*fltBoxCarSize;
      } else {
        threshold = (int) detChannel.GetThreshold();
      }
      ADCVar = detChannel.GetADCVariance();
    }
    else {
      baseline = fBaseline;
      ADCVar  = ( meanBgPe * (1 + gainVariance) * absGain*absGain )
        + elecNoiseVariance;
      if (meanBgPe) {
        if (!fUseNewThreshold) {
          const double baseN  = baseline * fltBoxCarSize;
          const double sigmaN = sqrt(ADCVar * fltBoxCarSize);
          const double referenceADCBinSize = 100*ns;
          const double cameraADCBinSize = detCamera.GetFADCBinSize();
          // equivalent numbers for 100 ns binning
          // (for which EvalRedThreshold was made for)
          const unsigned int equivalentBoxCarSize =
            fltBoxCarSize * (cameraADCBinSize/referenceADCBinSize);
          const double equivalentMeanPe =
            meanBgPe;// * (referenceADCBinSize/cameraADCBinSize);
          threshold = sigmaN*EvalRedThreshold(equivalentMeanPe,
                                              equivalentBoxCarSize,
                                              gainVariance) + baseN;
        } else {
          //New threshold values are baseline subtracted as the baseline is a free parameter in the simulation
          threshold = fSimThresholdValues[thresholdConfigSignature].Y(meanBgPe);
          threshold += fltBoxCarSize * baseline;
        }
      }
      else { // meanBgPe==0
        // this is for noisefree simulations
        threshold = baseline*fltBoxCarSize + sqrt(ADCVar)*3 + 1;
      }
    }

    fCalibCorrection[pixelId] = calibCorrection;
    fOpticalEfficiencyCorrection[pixelId] = opticalEfficiencyCorrection;

    //we should not correct measured quantities
    if (!fUseMonitoring) {
      threshold = int(baseline*fltBoxCarSize + (threshold-baseline*fltBoxCarSize)*calibCorrection + .5);
      ADCVar *= calibCorrection * calibCorrection;
    }
    
    pxsimdata.SetMean(baseline);
    pxsimdata.SetRMS(sqrt(ADCVar));
    pxsimdata.SetThreshold(threshold);
    pxsimdata.SetNumSamples(fltBoxCarSize);


    if (fVerbosityLevel>2) {
      cout << " pixelId=" << pixelId
           << " gain=" << absGain
           << " var=" << ADCVar
           << " base=" << baseline
           << " thresh=" << threshold 
           << " ph_bg=" << backphotons
           << " meanPE=" << meanBgPe
           << " num=" << fltBoxCarSize
           << " caliCorr=" << calibCorrection
           << endl;
    }

  } // end of loop over pixels


  // all pixels have been DeSelected -> DeSelect whole camera
  if (!telHasAtLeastOneValidPixel)
    {
      ostringstream warn;
      warn << "********** All pixels of tel=" << tel.GetId() << " at eye=" << tel.GetEyeId()
           << " have invalid calibration constants -> switch OFF telescope **********";
      WARNING(warn);
      return false;
    }

  return true;

} // end of InitCamera


void
FdElectronicsSimulator::ElecSim(fevt::Telescope& tel)
{
  const fdet::FDetector& detFD = Detector::GetInstance().GetFDetector();
  const fdet::Telescope& detTel = detFD.GetTelescope(tel);
  const fdet::Camera& detCamera = detTel.GetCamera();
  
  const double timeBinSize    = detCamera.GetFADCBinSize();
  double cutoffFrequency      = detCamera.GetCutoffFrequency();

  const unsigned int SLTbin = detCamera.GetSLTTriggerBin();
  const double adcBinsPer100ns = detCamera.GetFADCBinSize() / (100.*ns);
  const unsigned int adcTriggerBin = (unsigned int)(SLTbin / adcBinsPer100ns);

  const unsigned int fadcTraceLength = detCamera.GetFADCTraceLength();

  fevt::TelescopeSimData& telSim = tel.GetSimData();
  const unsigned int traceSize = telSim.GetNumberOfPhotonBins()
    + adcTriggerBin                    // reserve space in front of the photon signal
    + (fadcTraceLength-adcTriggerBin); // reserve space after the photon signal

  // for the time convolution, we need a little larger trace
  const unsigned int nSafety = (fDoAntiAliasingFilter ? fErrFuncFactors.size() - 1 : 0);
  const unsigned int workTraceSize = traceSize + 2*nSafety;

  const unsigned int binShift = adcTriggerBin + nSafety;

  
  /*
    RU Mi 9. Aug 18:13:12 CEST 2017

    first loop over pixels with photon content and calculate signal after PMT and before 
    ADC for all channels and virtual channels:
   */
  
  map<unsigned int, vector<double>> traceMap;

  // loop all pixels
  for (unsigned int pixelId=1; pixelId<=detTel.GetLastPixelId(); pixelId++) {

    if (tel.HasPixel(pixelId)) {

      fevt::Pixel& pixel = tel.GetPixel(pixelId);
   
      if (pixel.GetStatus() != fevt::ComponentSelector::eDeSelected &&
          pixel.HasSimData()) {
        
        fevt::PixelSimData& pxsimdata = pixel.GetSimData();
        
        const fdet::Channel& detChannel = detTel.GetChannel(pixel);
        const unsigned int channelId = detChannel.GetId();
        const double gainVariance = detChannel.GetGainVariance(); // of PMT (!)
        
        const double meanPeBgr = fPeBgNoiseTable[pixelId]; // background light PE
        
        // AND the corresponding virtual pixel (low gain)
        const unsigned int virtualChannelId = detChannel.GetVirtualChannelId();
        
        if (fVerbosityLevel>2) {
          cout << " elec: "
               << " eyeId=" << tel.GetEyeId()
               << " telId=" << tel.GetId()
               << " pixelId=" << pixelId
               << " channelId=" << channelId
               << " SLTbin=" << SLTbin
               << " FADClength=" << fadcTraceLength
               << " PhotonTrace=" << telSim.GetNumberOfPhotonBins()
               << " simElecTrace=" << traceSize
               << endl;
        }

        // the simulated pixel light to (nPE*gain) - trace
        {
          if (traceMap.count(channelId)!=0) {
            ERROR("*** normal camera channel is used twice !!!  ***");
          }
          vector<double> originalChannelTrace(workTraceSize, 0);
          traceMap[channelId] = originalChannelTrace;
          
          vector<double> originalChannelTraceVirtual(workTraceSize, 0);
          if (traceMap.count(virtualChannelId)==0) {
            traceMap[virtualChannelId] = originalChannelTraceVirtual;
          }
        }
        
        vector<double>& tracePix = traceMap[channelId];
        vector<double>& traceVPix = traceMap[virtualChannelId];
        
        const double quantEff = fQEffAtNormWavelength[pixelId];
        
        // now applying fluctuations from PMT
        vector<double> gaussRandomArray(workTraceSize);
        RandGauss::shootArray(&fRandomEngine->GetEngine(), gaussRandomArray.size(),
                              &gaussRandomArray.front(), 0., 1.);
        

        TraceD const * photonTraceOrg = 0;
        unsigned int photonSize = 0;
        if (pxsimdata.HasPhotonTrace(fevt::FdConstants::eTotal)) {
          photonTraceOrg = &pxsimdata.GetPhotonTrace(fevt::FdConstants::eTotal);
          photonSize = photonTraceOrg->GetSize();
          if (fVerbosityLevel>3) 
            cout << "photonSize=" << photonSize << endl;
          if (photonSize+binShift > workTraceSize) {
            ERROR (" ***** simulated photon trace will be truncated!!!! ****");            
          }
        }
        
        // use photon signal, add PE noise from photon background
        for (unsigned int ti=0; ti<workTraceSize; ti++) {
          
          double nPePh = 0;
          // the photon pulse is shifted by binShift
          if (photonTraceOrg && ti>=binShift) {
            const unsigned int ti_shift = ti-binShift;
            if (ti_shift<photonSize) {
              nPePh = (*photonTraceOrg)[ti_shift] * quantEff;
              //apply optical efficiency loss correction
              if (!detChannel.IsVirtual()) {
                nPePh /= fOpticalEfficiencyCorrection[pixelId];
              }
            }
          }
          
          double nPe = nPePh + meanPeBgr;
          
          if (fVerbosityLevel>3) {
            cout << "pe-noise: channelId=" << channelId
                 << " virt-channelId=" << virtualChannelId
                 << " ti=" << ti
                 << " nPe=" << nPe
                 << " nPePh=" << nPePh
                 << " meanPeBgr=" << meanPeBgr
                 << " fPeFluctuations=" << fPeFluctuations
                 << " gauss=" << gaussRandomArray[ti]
                 << " gainVariance=" << gainVariance
                 << endl;
          }
          
          // next: add PE poisson fluctuations
          if (fPeFluctuations)
            nPe = RandPoisson::shoot(&fRandomEngine->GetEngine(), nPe);
          // add gaussian fluctuation due to gain
          nPe += gaussRandomArray[ti] * sqrt(nPe * gainVariance);
          nPe -= meanPeBgr; // BG is added explicitly below
          
          tracePix[ti]   = nPe;
          traceVPix[ti] += nPe; // add this up
        } // end loop time bins
      } // pixel has sim data
    } // has pixel
  }// end loop over Pixels
  
  
  /*
    RU Mi 9. Aug 18:13:12 CEST 2017

    Loop over channel-map and apply electronics gain, variance and filter smoothing
    for both normal and virtual channels. Also set saturation bit!

    copy traces into fevt::ChannelSimData Objects
   */

  const unsigned int fadcDynamicRange = detCamera.GetADCDynamicRange();
  const int maxADC = (1<<fadcDynamicRange) - 1;
  
  for (auto& iChannelMap : traceMap) {
    
    const unsigned int channelId = iChannelMap.first;
    vector<double>& channelTrace = iChannelMap.second;
    const unsigned int binMax = channelTrace.size();
    
    const fdet::Channel& detChannel = detTel.GetChannel(channelId);
    const double fadcConvert = detChannel.GetElectronicsGain();
    
    const double elecNoiseVariance = detChannel.GetElectronicNoiseVariance();
    const double noiseEqBdw = 0.5*cutoffFrequency*sqrt(kPi/log(2.0));
    const double samplingFrequency = 1./timeBinSize;
    const double analogElectronicsFactor = 2. * noiseEqBdw / samplingFrequency;
    
    const double sqrtNoiseVar = sqrt(elecNoiseVariance / analogElectronicsFactor *
                                     fBandWidthFactor);
    

    if (fVerbosityLevel>2) {
      cout << "channelId=" << channelId
           << " sqrtNoiseVar=" << sqrtNoiseVar
           << " workTraceSize=" << workTraceSize
           << " nSafety=" << nSafety
           << " binMax=" << binMax
           << " fadcConvert=" << fadcConvert
           << endl;
    }
    
    // fluctuations from electronics    
    vector<double> elecNoiseArray(workTraceSize);
    RandGauss::shootArray(&fRandomEngine->GetEngine(), elecNoiseArray.size(),
                          &elecNoiseArray.front(), 0., sqrtNoiseVar);
    
    // the normal channel trace:    
    for (unsigned int ti=0; ti<binMax; ti++) {
      channelTrace[ti] *= fadcConvert;
      channelTrace[ti] += elecNoiseArray[ti];
    }
    
    // convolute trace with gaussian response of anti-aliasing filter
    vector<double> convolutedPMTTrace(traceSize, 0.);
    if (fDoAntiAliasingFilter) {
      DoTimeConvolution(channelTrace, convolutedPMTTrace, nSafety);
    }
    else {
      for ( unsigned int ti=0; ti<traceSize; ti++) {
        convolutedPMTTrace[ti] = channelTrace[ti+nSafety];
      }
    }
    
    // get the baseline and calib-correction factor per pixel
    double baseline = fBaseline; // for virtual channels just use this
    double calibCorrection = 1; // does not exist for virtual channels
    if (!detChannel.IsVirtual()) {
      const int pixelId = detChannel.GetPixelId();
      fevt::Pixel& evtPix = tel.GetPixel(pixelId);
      fevt::PixelSimData& pxsimdata = evtPix.GetSimData();
      baseline = pxsimdata.GetMean(); // from DB or xml
      calibCorrection = fCalibCorrection[pixelId];
      
      // CalibCorrection: ONLY for real pixels!
      for (unsigned int ti=0; ti<traceSize; ti++) {
        convolutedPMTTrace[ti] *= calibCorrection;
      }
    }
    
    
    // create final digital FADC trace
    // create new data structures
    if (!tel.HasChannel(channelId))
      tel.MakeChannel(channelId);
    fevt::Channel& channel = tel.GetChannel(channelId);
    
    if (!channel.HasSimData())
      channel.MakeSimData();
    fevt::ChannelSimData& channelSimData = channel.GetSimData();
    
    if (!channelSimData.HasFADCTrace(fevt::FdConstants::eTotal))
      channelSimData.MakeFADCTrace(traceSize, timeBinSize, fevt::FdConstants::eTotal);
    TraceI& digitTrace = channelSimData.GetFADCTrace(fevt::FdConstants::eTotal);
        
    // now convert temporary pmtTrace into final FADC trace
    for (unsigned int ti=0; ti<traceSize; ti++) {
      digitTrace[ti] = min(max(0, int(convolutedPMTTrace[ti] + baseline + 0.5)), maxADC);
    } // end loop bins    
  } // loop channels
}// end of ElecSim


// --------------------------------------------------------------------------------------
void
FdElectronicsSimulator::PrepareTimeConvolution(const fdet::Telescope& detTel)
{
  const fdet::Camera& detCamera = detTel.GetCamera();
  const double timeBinSize = detCamera.GetFADCBinSize();

  // see IEEE Trans.Nucl.Sci.50:1208-1213,2003
  const double cutoffFrequency = detCamera.GetCutoffFrequency();
  const double tau = sqrt(log(2.)) / (kTwoPi * cutoffFrequency);
  const double noiseEqBdw = 0.5 * cutoffFrequency * sqrt(kPi/log(2.));
  const double samplingFrequency = 1./timeBinSize;
  const double analogElectronicsFactor = 2 * noiseEqBdw / samplingFrequency;

  fErrFuncFactors.clear();

  const double epsilon = 1e-6; // precision
  const int maxIter = 50; // [bins]
  int iter = 0;

  double sumE = 0.;
  double sumESquare = 0.;

  // integrate the gaussian in timebins
  while ( sumE < 1 - epsilon &&
          iter < maxIter ) {

    const double t1 = (iter - 0.5) * timeBinSize;
    const double t2 = (iter + 0.5) * timeBinSize;

    const double E1 = 0.5 * erf(t1 / (sqrt(2.)*tau));
    const double E2 = 0.5 * erf(t2 / (sqrt(2.)*tau));
    const double E = E2 - E1;

    fErrFuncFactors.push_back(E);

    if ( iter == 0 ) {
      sumESquare += E*E;
      sumE += E;
    }
    else {
      sumESquare += 2*E*E;
      sumE += 2*E;
    }

    ++iter;
  }

  if ( iter+1 >= maxIter || fErrFuncFactors.empty() ) {
    ostringstream msg;
    msg << " Error filling fErrFuncFactors !!! iter="
        << iter << ",  sumE=" << sumE;
    ERROR(msg);
    throw AugerException(msg.str());
  }

  // correction factor for approximate formula used in
  // FdBackGroundSimulator
  fBandWidthFactor = analogElectronicsFactor / sumESquare;

  return;
}

void
FdElectronicsSimulator::DoTimeConvolution(const vector<double>& originalTrace,
                                          vector<double>& convolutedTrace,
                                          unsigned int timeOffset)
{
  const int convolutedTraceSize = convolutedTrace.size();
  const int originialTraceSize  = originalTrace.size();
  const int nConv = (int) fErrFuncFactors.size();

  // check trace sizes...
  if ( originialTraceSize-convolutedTraceSize < 2*(nConv-1) ) {
    // this should never happend !!!!
    ostringstream msg;
    msg << " Trace size mismatch!!! "
        << convolutedTraceSize << " " << originialTraceSize << " " << fErrFuncFactors.size();
    ERROR(msg);
    throw AugerException(msg.str());
  }

  for ( int i=-nConv+1; i<nConv; ++i ) {
    const double w  = fErrFuncFactors[abs(i)];
    for ( int tConv = 0; tConv<convolutedTraceSize; ++tConv ) {
      const int tOrig = tConv + timeOffset + i;
      convolutedTrace[tConv] += w * originalTrace[tOrig];
    }
  }
}



//-----------  see  GAP-2003-083

// threshold should produce a 100Hz FLT trigger rate (in data):
// dT = 1000*100ns = 100us, nFLT=1sec/100us/10=1e3, n(100Hz)=1sec*100Hz = 1e2, pFLT=

double
FdElectronicsSimulator::EvalRedThreshold (double npe, int nSamp,
                                          double gainVariance)
{
  // mean values of photoelectrons at photocathode
  const int peTableSize = 12;
  const double maxNpe = 10000.0;
  const double vecNpe[peTableSize] = { 2.0, 3.0, 5.0, 10.0, 20.0, 30.0, 50.0, 100.0, 200.0, 500.0, 2000.0, 10000.0};

  // cutoff values [P (>cutoff) = 1.e-5] for pure Poisson (npe)
  const int maxSamples = 20;
  const double fgThrFactPG00[maxSamples][peTableSize] = {
   {5.62553, 5.47466, 5.22473, 4.95916, 4.76080, 4.67745, 4.58724, 4.49212, 4.42888, 4.36924, 4.31726, 4.27359},
    {5.33020, 5.14876, 4.95916, 4.76080, 4.62016, 4.55890, 4.49212, 4.42888, 4.38149, 4.33885, 4.30203, 4.26984},
    {5.14876, 4.97170, 4.83641, 4.67745, 4.55890, 4.50640, 4.45387, 4.39864, 4.36009, 4.32552, 4.29530, 4.26818},
    {5.01369, 4.88375, 4.76080, 4.62016, 4.52161, 4.47194, 4.42888, 4.38149, 4.34717, 4.31726, 4.29122, 4.26718},
    {4.95916, 4.83641, 4.71542, 4.58724, 4.49212, 4.45387, 4.41152, 4.36924, 4.33885, 4.31192, 4.28845, 4.26650},
    {4.88375, 4.78669, 4.67745, 4.55890, 4.47194, 4.43772, 4.39864, 4.36009, 4.33231, 4.30772, 4.28642, 4.26600},
    {4.83906, 4.75062, 4.64813, 4.53423, 4.45954, 4.42361, 4.38857, 4.35303, 4.32738, 4.30466, 4.2848, 4.26562},
    {4.81387, 4.71102, 4.62016, 4.52161, 4.44715, 4.41481, 4.38149, 4.34717, 4.32337, 4.30203, 4.28353, 4.26531},
    {4.78669, 4.69794, 4.60003, 4.50640, 4.43772, 4.40640, 4.37459, 4.34290, 4.32022, 4.30000, 4.28248, 4.26505},
    {4.76080, 4.67745, 4.58724, 4.49212, 4.42889, 4.39864, 4.36924, 4.33885, 4.31726, 4.29813, 4.28158, 4.26483},
    {4.73603, 4.65687, 4.56868, 4.48066, 4.42147, 4.39268, 4.36448, 4.33546, 4.31500, 4.29665, 4.28079, 4.26465},
    {4.71102, 4.63861, 4.55890, 4.47194, 4.41482, 4.38738, 4.36009, 4.33231, 4.31285, 4.29530, 4.28011, 4.26449},
    {4.69514, 4.62323, 4.54854, 4.46522, 4.40773, 4.38301, 4.35605, 4.32974, 4.31094, 4.29409, 4.27952, 4.26435},
    {4.68729, 4.61054, 4.53423, 4.45954, 4.40377, 4.37855, 4.35303, 4.32738, 4.30918, 4.29303, 4.27898, 4.26422},
    {4.67745, 4.60003, 4.52760, 4.45387, 4.39865, 4.37460, 4.35001, 4.32552, 4.30772, 4.29210, 4.27851, 4.26410},
    {4.66370, 4.59111, 4.52161, 4.44715, 4.39506, 4.37123, 4.34717, 4.32337, 4.30649, 4.29122, 4.27809, 4.26400},
    {4.64398, 4.58317, 4.51415, 4.44087, 4.39013, 4.36829, 4.34518, 4.32175, 4.30512, 4.29046, 4.27769, 4.26391},
    {4.63861, 4.57557, 4.50640, 4.43772, 4.38739, 4.36546, 4.34290, 4.32022, 4.30409, 4.28970, 4.27734, 4.26382},
    {4.63335, 4.56769, 4.49896, 4.43135, 4.38454, 4.36238, 4.34077, 4.31873, 4.30305, 4.28904, 4.27700, 4.26374},
    {4.62016, 4.55890, 4.49212, 4.42888, 4.38149, 4.36009, 4.33885, 4.31726, 4.30203, 4.28845, 4.27670, 4.26367}
  };

  // cutoff values [P (>cutoff) = 1.e-5] for Poisson (npe) convoluted with
  // Gauss (0,0.41*npe)
  const double fgThrFactPG41[maxSamples][peTableSize] = {
    {5.47757, 5.18121, 4.90285, 4.69526, 4.55630, 4.50512, 4.45118, 4.39478, 4.35642, 4.32268, 4.29289, 4.27359},
    {5.16834, 4.92291, 4.74034, 4.56101, 4.47433, 4.43379, 4.39478, 4.35647, 4.32962, 4.30570, 4.28452, 4.26984},
    {4.99797, 4.83452, 4.64616, 4.50969, 4.43384, 4.40344, 4.37184, 4.33990, 4.31792, 4.29814, 4.28077, 4.26818},
    {4.92332, 4.74995, 4.60145, 4.47870, 4.41033, 4.38180, 4.35648, 4.32966, 4.31063, 4.29349, 4.27855, 4.26718},
    {4.87796, 4.71279, 4.56867, 4.45550, 4.39485, 4.37185, 4.34742, 4.32287, 4.30584, 4.29056, 4.27699, 4.26650},
    {4.80947, 4.67208, 4.54720, 4.43799, 4.38185, 4.36244, 4.33991, 4.31796, 4.30203, 4.28824, 4.27589, 4.26600},
    {4.79474, 4.65794, 4.53000, 4.42349, 4.37572, 4.35394, 4.33437, 4.31401, 4.29925, 4.28651, 4.27502, 4.26562},
    {4.75745, 4.63742, 4.51394, 4.41430, 4.36841, 4.34784, 4.32967, 4.31068, 4.29710, 4.28511, 4.27431, 4.26531},
    {4.72165, 4.61372, 4.50083, 4.40739, 4.36249, 4.34402, 4.32560, 4.30817, 4.29538, 4.28396, 4.27372, 4.26505},
    {4.71184, 4.60127, 4.48892, 4.39867, 4.35654, 4.33992, 4.32288, 4.30589, 4.29363, 4.28298, 4.27324, 4.26483},
    {4.69720, 4.58849, 4.47548, 4.39364, 4.35292, 4.33668, 4.32026, 4.30401, 4.29233, 4.28209, 4.27282, 4.26465},
    {4.68197, 4.56997, 4.46975, 4.38555, 4.34789, 4.33354, 4.31797, 4.30207, 4.29126, 4.28137, 4.27244, 4.26449},
    {4.66760, 4.56238, 4.46344, 4.38306, 4.34539, 4.33087, 4.31584, 4.30070, 4.29019, 4.28071, 4.27212, 4.26435},
    {4.65437, 4.55571, 4.45383, 4.37930, 4.34268, 4.32835, 4.31402, 4.29930, 4.28931, 4.28013, 4.27182, 4.26422},
    {4.64203, 4.54710, 4.45082, 4.37543, 4.33997, 4.32561, 4.31219, 4.29833, 4.28839, 4.27959, 4.27157, 4.26410},
    {4.63008, 4.53821, 4.44349, 4.37190, 4.33721, 4.32432, 4.31069, 4.29714, 4.28761, 4.27915, 4.27134, 4.26400},
    {4.61792, 4.52983, 4.43929, 4.36876, 4.33495, 4.32208, 4.30926, 4.29630, 4.28699, 4.27872, 4.27112, 4.26391},
    {4.60487, 4.52228, 4.43568, 4.36588, 4.33359, 4.32067, 4.30817, 4.29542, 4.28640, 4.27829, 4.27092, 4.26382},
    {4.59857, 4.51564, 4.42911, 4.36302, 4.33156, 4.31938, 4.30699, 4.29457, 4.28581, 4.27794, 4.27074, 4.26374},
    {4.59471, 4.50988, 4.42631, 4.35986, 4.32972, 4.31798, 4.30589, 4.29368, 4.28526, 4.27759, 4.27058, 4.26367}
  };

  if ( nSamp > maxSamples ) {
    ostringstream msg;
    msg << " number of samples to large "
        << nSamp << " max=" << maxSamples;
    ERROR(msg);
    throw OutOfBoundException(msg.str());
  }

  if ( npe > maxNpe ) {
    ostringstream msg;
    msg << " number of photo electrons to high "
        << npe << " max=" << maxNpe;
    ERROR(msg);
    throw OutOfBoundException(msg.str());
  }

  int i2 = 0;

  for (int i=0; i<peTableSize; i++)
    if (npe > vecNpe[i])
      i2 = i + 1;

  if (i2 == 0)
    i2 = 1;
  if (i2 == peTableSize)
    i2 = peTableSize-1;
  int i1 = i2 - 1;

  double fac00 = (fgThrFactPG00[nSamp-1][i2]
                  - fgThrFactPG00[nSamp-1][i1])/(log(vecNpe[i2])
                         - log(vecNpe[i1])) * (log(npe) - log(vecNpe[i1]));

  fac00 += fgThrFactPG00[nSamp-1][i1];

  double fac41 = (fgThrFactPG41[nSamp-1][i2] - fgThrFactPG41[nSamp-1][i1])
    /(log(vecNpe[i2]) - log(vecNpe[i1]))*(log(npe) - log(vecNpe[i1]));

  fac41 += fgThrFactPG41[nSamp-1][i1];

  const double fac = fac00 + (fac41 - fac00)*gainVariance/0.41;

  return fac;

}// end of EvalThreshold


// Configure (x)emacs for this file ...
// Local Variables:
// mode:c++
// compile-command: "cd $AUGER_BASE && make"
// End: