RPCElectronicsSimulator.cc

Go to the documentation of this file.

#include "RPCElectronicsSimulator.h"
using namespace RPCElectronicsSimulatorLX;

#include <iostream>
#include <fstream>
#include <iterator>

#include <det/Detector.h>
#include <cdet/CDetector.h>
#include <cdet/Station.h>

#include <evt/Event.h>
#include <evt/ShowerSimData.h>
using namespace evt;

#include <fwk/CentralConfig.h>
#include <fwk/RunController.h>
#include <fwk/LocalCoordinateSystem.h>
#include <fwk/CoordinateSystemRegistry.h>
#include <fwk/RandomEngineRegistry.h>

using namespace fwk;

#include <cevt/CEvent.h>
using namespace cevt;

#include <utl/config.h>
#include <utl/Trace.h>
#include <utl/AugerUnits.h>
#include <utl/Math.h>
#include <utl/MathConstants.h>
#include <utl/PhysicalConstants.h>
#include <utl/TimeStamp.h>
#include <utl/TabulatedFunction.h>
#include <utl/TabularStream.h>
#include <utl/ErrorLogger.h>
#include <utl/Particle.h>
#include <utl/AugerCoordinateSystem.h>
#include <utl/UTMPoint.h>
#include <utl/MultiTimeDistribution.h>
#include <utl/GeometryUtilities.h>
#include <utl/RandomEngine.h>

using namespace utl;

#include <CLHEP/Random/RandGauss.h>

#include <sstream>
#include <string>
#include <vector>

using CLHEP::RandGauss;

RPCElectronicsSimulator::RPCElectronicsSimulator() :
  fVerbose(0),
  fChargeThreshold(0.0),
  //  fDeadTime(0.0),
  fTimeResolution(0.0), 
  fGPSTimeJitter(0.0), 
  fRandomEngine(0)
{
}

RPCElectronicsSimulator::~RPCElectronicsSimulator() 
{;}

VModule::ResultFlag RPCElectronicsSimulator::Init(void) {

  std::ostringstream str;

  str << "RPCElectronicsSimulator::Init() \n";

  Branch topB =
    CentralConfig::GetInstance()->GetTopBranch("RPCElectronicsSimulator");

  topB.GetChild("Verbosity").GetData(fVerbose);
  topB.GetChild("ChargeThreshold").GetData(fChargeThreshold);
  //  topB.GetChild("DeadTime").GetData(fDeadTime);
  topB.GetChild("TimeResolution").GetData(fTimeResolution);
  topB.GetChild("GPSTimeJitter").GetData(fGPSTimeJitter);

  fRandomEngine = 
    &RandomEngineRegistry::GetInstance().Get(RandomEngineRegistry::eDetector);
 
  str << "ChargeThreshold = "    << fChargeThreshold/(pico*coulomb) << " pC \n"
      << "TimeResolution =  "    << fTimeResolution/ns              << " ns \n" 
      << "GPSTimeJitter =  " << fGPSTimeJitter/ns           << " ns \n" 
    //      << "DeadTime =        "    << fDeadTime/ns                    << " ns " 
      << "\n";
  
  INFO(str);
  
  return eSuccess;
}


VModule::ResultFlag RPCElectronicsSimulator::Finish(void) {
  
  INFO("RPCElectronicsSimulator::Finish()");
   
  return eSuccess;

}


VModule::ResultFlag RPCElectronicsSimulator::Run(Event &event){

  INFO("RPCElectronicsSimulator::Run()");

  if (!event.HasCEvent()) {
      ERROR("RPCElectronicsSimulator::No CEvent!");
      return eFailure;
  }


  const cdet::CDetector & detector = det::Detector::GetInstance().GetCDetector();

  CEvent &theCEvent = event.GetCEvent();

  TabularStream tab("|c|c|.|.|.|.|");

  if (fVerbose > 0) {
    tab << HLine('=');
    tab << "station" << endc 
        << "pad" << endc 
        << "time" << endc 
        << "charge" << endc 
        << "hits. time" << endc 
        << "time alive" << endr;

    tab << endc 
        << endc
        << "[ns]" << endc
        << "[pC]" << endc
        << "[ns]" << endc
        << "[ns]" << endr
        << HLine('=');
  }

  for (CEvent::StationIterator cIt = theCEvent.StationsBegin();
       cIt != theCEvent.StationsEnd(); ++cIt) {
    
    Station& rpcStation = *cIt;    
    
    if (rpcStation.GetNumberActivePads() == 0)
      continue;

    if (fVerbose > 0) {
      std::cout << "\n \n"
           << "Number of muons/active pads in station " 
           << rpcStation.GetId() << " "
           << rpcStation.GetSimData().GetNumberOfMuons() << "/" 
           << rpcStation.GetNumberActivePads() 
           << std::endl;
    }
    
    const float deadTime = detector.GetStation(rpcStation).GetDeadTime();

    // Smearing in the definition of the time of the station due to the GPS resolution
    const double timeGPS = 
      RandGauss::shoot(&fRandomEngine->GetEngine(),0.0,fGPSTimeJitter);

    // Loop over pads with signal for this station
    for(Station::PadIterator pIt = rpcStation.PadsBegin(); 
        pIt != rpcStation.PadsEnd(); ++pIt) {
      
      Pad& pad = *pIt;
      PadSimData& padSimData = pad.GetSimData();

      // Loop over components of the charge time distributions
      for (PadSimData::ChargeTimeDistributionsIterator tIt =
             padSimData.ChargeTimeDistributionsBegin();
           tIt != padSimData.ChargeTimeDistributionsEnd(); ++tIt) {
        
        const Station::SignalComponent component = 
          static_cast<Station::SignalComponent>(tIt->GetLabel());
        
        // Simulate signal digitization 
        TimeDistributionD& chargeDist = tIt->GetTimeDistribution();
        TimeDistributionD::SparseIterator qIt = chargeDist.SparseBegin();

        const double timeBin = chargeDist.GetBinning();

        double timeVeto = qIt->get<0>()*timeBin;

        for (; qIt != chargeDist.SparseEnd(); ++qIt) {

          const double time0 =  qIt->get<0>() * timeBin;
          const double charge = qIt->get<1>();

          if (fVerbose > 0 && component == Station::eTotal) {
            tab << ' ' << rpcStation.GetId() << ' ' << endc
                << ' ' << pad.GetId() << ' ' << endc 
                << time0/nanosecond << endc 
                << charge/(pico*coulomb)  << ' '  << endc;
          }

          if (charge > fChargeThreshold && time0 >= timeVeto) { 

            // Smear time
            const double time = 
              RandGauss::shoot(&fRandomEngine->GetEngine(),time0,fTimeResolution);        
            if (!padSimData.HasHitsTimeDistribution(component))
              padSimData.MakeHitsTimeDistribution(component);

            padSimData.GetHitsTimeDistribution(component).AddTime(time+timeGPS);
        
            timeVeto = time + deadTime;

            if (fVerbose > 0 && component == Station::eTotal) 
              tab << (time+timeGPS)/nanosecond << endc
                  << timeVeto/nanosecond << endr;
          }

          else {
            // No trigger
            if (fVerbose > 0 && component == Station::eTotal) 
              tab << " -- " << endc << timeVeto/nanosecond << endr;
          }
        }
        
      } // end loop over signal components
    } // end loop over station pads

    if (fVerbose > 0) 
      tab << hline;

  } // end loop over stations


  if (fVerbose > 0) 
    std::cout << tab << std::endl;

  return eSuccess;
}