RPCSimulator.cc

Go to the documentation of this file.

#include "RPCSimulator.h"
#include "RPCChargeGenerator.h"

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

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

#include <evt/Event.h>
#include <evt/ShowerSimData.h>

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

#include <cevt/CEvent.h>

#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/ErrorLogger.h>
#include <utl/Particle.h>
#include <utl/AugerCoordinateSystem.h>
#include <utl/UTMPoint.h>
#include <utl/MultiTimeDistribution.h>
#include <utl/GeometryUtilities.h>

#include <CLHEP/Random/RandGauss.h>

#include <TH1.h>
#include <TProfile.h>
#include <TFile.h>

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

using namespace RPCSimulatorLX;
using namespace det;
using namespace evt;
using namespace fwk;
using namespace cevt;
using namespace utl;
using namespace std;
using CLHEP::RandGauss;


RPCSimulator::~RPCSimulator()
{
  delete fChargeHisto;
  delete fChargeProfileHisto;
  delete fChargeGenerator;
}


VModule::ResultFlag RPCSimulator::Init()
{
  INFO("RPCSimulator::Init()");

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

  const AttributeMap useAtt({{ "use", "yes" }});
  Branch rootB = topB.GetChild("RootHistogramsFilename", useAtt);
  if (rootB) {
    fUseRootHistos = true;
    rootB.GetData(fRootHistoFilename);
    fChargeHisto = new TH1D("chargeDist", "charge distribution", 200, 0, 100);
    fChargeProfileHisto = new TProfile("AverageChargeVsAngle", "Average charge vs angle of incidence at RPC", 90, 0, 90);
  }

  Branch chargeB = topB.GetChild("ChargeGenerator", useAtt);
  if (chargeB) {
    fUseChargeDistribution = true;
    const double k0 = chargeB.GetChild("k_param").Get<double>();
    const double th0 = chargeB.GetChild("th_param").Get<double>();
    const double referenceAngle = chargeB.GetChild("ref_angle").Get<double>();
    fChargeGenerator = new RPCChargeGenerator(k0, th0, referenceAngle);
  } else {
    WARNING("PDF for charge generation in the RPC switched off! "
            "A fixed charge of 3.1415 pC will be set for each particle crossing the gas!");
  }

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

  return eSuccess;
}


VModule::ResultFlag
RPCSimulator::Finish()
{
  INFO("RPCSimulator::Finish()");

  if (fUseRootHistos) {
    TFile outFile(fRootHistoFilename.c_str(), "RECREATE");

    fChargeHisto->Write();
    fChargeProfileHisto->Write();
    outFile.Write();
    outFile.Close();

    delete fChargeHisto;
    delete fChargeProfileHisto;

    fChargeHisto = nullptr;
    fChargeProfileHisto = nullptr;
  }

  delete fChargeGenerator;
  fChargeGenerator = nullptr;

  return eSuccess;
}


VModule::ResultFlag
RPCSimulator::Run(Event& event)
{
  INFO("RPCSimulator::Run()");

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

  const sdet::SDetector& sDetector = det::Detector::GetInstance().GetSDetector();
  const cdet::CDetector& cDetector = det::Detector::GetInstance().GetCDetector();

  CEvent &cEvent = event.GetCEvent();

  for (CEvent::StationIterator cIt = cEvent.StationsBegin(); cIt != cEvent.StationsEnd(); ++cIt) {

    Station& station = *cIt;

    if (!station.HasSimData()) {
      cout << "No sim data for station " << station.GetId() << " !!!!";
      continue;
    }

    StationSimData &simData = station.GetSimData();

    if (simData.ParticlesBegin() == simData.ParticlesEnd())
      continue;

    const CoordinateSystemPtr cs = sDetector.GetStation(station.GetId()).GetLocalCoordinateSystem();

    const cdet::Station& dStation = cDetector.GetStation(station);

    for (StationSimData::ParticleIterator partIt = simData.ParticlesBegin();
         partIt != simData.ParticlesEnd(); ++partIt) {

      const Particle& part = *partIt;

      // Consider only partices with energy deposited in the gas of the RPC
      // Energy deposited is stored in the weight data member of class utl::Particle

      if (part.GetWeight() <= 0) // Consider only particles depositing energy in the gas
        continue;

      const Point& partPosition = part.GetPosition();
      const Vector& partMomentum = part.GetDirection();

      const double particleTime = part.GetTime().GetInterval();

      // CBT start by creating the trace of particles in each pad

      const unsigned padId = dStation.GetPadId(partPosition);

      if (!padId) {
        WARNING("Pad Id is zero !!!");
        continue;
      }

      if (!station.HasPad(padId))
        station.MakePad(padId);

      Pad& pad = station.GetPad(padId);
      if (!pad.HasSimData())
        pad.MakeSimData();

      PadSimData& padSimData = pad.GetSimData();

      // Fill time distributions of particles arriving at the gas of the RPC, with dE>0.0

      if (!padSimData.HasParticleTimeDistribution())
        padSimData.MakeParticleTimeDistribution();

      TimeDistributionI& particleDist1 = padSimData.GetParticleTimeDistribution();
      particleDist1.AddTime(particleTime);

      // Make time distribution for any electron at RPC
      if (part.GetType() == utl::Particle::eElectron ||
          part.GetType() == utl::Particle::ePositron) {

        if (!padSimData.HasParticleTimeDistribution(Station::eElectron))
          padSimData.MakeParticleTimeDistribution(Station::eElectron);

        padSimData.GetParticleTimeDistribution(Station::eElectron).AddTime(particleTime);
      }

      // Make time distribution for any photon at RPC
      if (part.GetType() == utl::Particle::ePhoton) {

        if (!padSimData.HasParticleTimeDistribution(Station::ePhoton))
          padSimData.MakeParticleTimeDistribution(Station::ePhoton);

        padSimData.GetParticleTimeDistribution(Station::ePhoton).AddTime(particleTime);
      }

      // Now for more detailed components
      Station::SignalComponent component = Station::GetSignalComponent(part);

      if (!padSimData.HasParticleTimeDistribution(component))
        padSimData.MakeParticleTimeDistribution(component);

      TimeDistributionI& particleDist2 = padSimData.GetParticleTimeDistribution(component);
      particleDist2.AddTime(particleTime);

      // its ptheta direction (momentum direction)
      const double theta = partMomentum.GetTheta(cs);

      const double pC = utl::pico*coulomb;

      // Generate charge according to particle length inside RPC (i.e. prop to sec(theta))
      const double charge = fUseChargeDistribution ? fChargeGenerator->GetCharge(theta) : kPi*pC;

      if (fUseRootHistos) {
        fChargeHisto->Fill(charge / pC);
        fChargeProfileHisto->Fill(theta / deg, charge / pC);
      }

      const double time = particleTime;

      // Fill charge distributions of particles arriving at the gas of the RPC, with dE>0.0

      if (!padSimData.HasChargeTimeDistribution())
        padSimData.MakeChargeTimeDistribution();

      TimeDistributionD& chargeDist1 = padSimData.GetChargeTimeDistribution();
      chargeDist1.AddTime(time, charge);

      // Make charge time distribution for any electron at RPC
      if (part.GetType() == utl::Particle::eElectron ||
          part.GetType() == utl::Particle::ePositron) {

        if (!padSimData.HasChargeTimeDistribution(Station::eElectron))
          padSimData.MakeChargeTimeDistribution(Station::eElectron);

        padSimData.GetChargeTimeDistribution(Station::eElectron).AddTime(time, charge);
      }

      // Make charge time distribution for any photon at RPC
      if (part.GetType() == utl::Particle::ePhoton) {

        if (!padSimData.HasChargeTimeDistribution(Station::ePhoton))
          padSimData.MakeChargeTimeDistribution(Station::ePhoton);

        padSimData.GetChargeTimeDistribution(Station::ePhoton).AddTime(time, charge);
      }

      // Now for more detailed components
      if (!padSimData.HasChargeTimeDistribution(component))
        padSimData.MakeChargeTimeDistribution(component);

      TimeDistributionD& chargeDist2 = padSimData.GetChargeTimeDistribution(component);
      chargeDist2.AddTime(time, charge);

    }
  }

  return eSuccess;
}