G4ASCIIAction.cc

Go to the documentation of this file.

#include "G4ASCIIAction.h"
#include "G4TankEventAction.h"
#include "G4TankSimulator.h"
#include "G4TankSteppingAction.h"
#include "G4TankTrackingAction.h"
#include "G4TankTrackInformation.h"
#include "ParticleCases.h"

#include <sevt/PMT.h>
#include <sevt/PMTConstants.h>
#include <sevt/PMTSimData.h>
#include <sevt/SEvent.h>
#include <sevt/Station.h>
#include <sevt/StationSimData.h>

#include <utl/ErrorLogger.h>
#include <utl/Particle.h>
#include <utl/Trace.h>
#include <utl/PhysicalConstants.h>

#include <fwk/RandomEngineRegistry.h>
#include <CLHEP/Random/RandFlat.h>
#include <CLHEP/Random/RandPoisson.h>
#include <CLHEP/Random/RandExponential.h>
#include <fwk/RunController.h>

#include <G4Event.hh>
#include <G4HCofThisEvent.hh>
#include <G4SDManager.hh>
#include <G4Step.hh>
#include <G4TouchableHistory.hh>
#include <G4Track.hh>
#include <G4Trajectory.hh>
#include <G4TrajectoryContainer.hh>
#include <G4VHitsCollection.hh>

#include <iostream>

using namespace std;
using namespace G4TankSimulatorASCII;
using namespace sevt;
using namespace utl;
using namespace fwk;


using CLHEP::RandFlat;
using CLHEP::RandPoisson;
using CLHEP::RandExponential;



G4ASCIIAction::G4ASCIIAction(const G4String name) :
  G4VSensitiveDetector(name)
{

  fG4TankSimulator =
    dynamic_cast<G4TankSimulator*>(&RunController::GetInstance().GetModule("G4TankSimulatorASCII"));

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


G4bool
G4ASCIIAction::ProcessHits(G4Step* theStep, G4TouchableHistory* /*R0Hist*/)
{
  // Fill ASCII  : use PMT4 to save the single photoelectron distribution

  if ( theStep->GetTrack()->GetDefinition()->GetPDGCharge() == 0 ) return true;

  if ( theStep->GetTrack()->GetVolume()->GetName() !="scin" )
    {

      return true;
    }

  const bool IsFromMuonDecay=theStep->GetTrack()->GetUserInformation();
  unsigned int fBarIndex=theStep->GetTrack()->GetVolume()->GetCopyNo();
  const unsigned int fPMTIndex=5; // Temporary.
  const double time = theStep->GetPreStepPoint()->GetGlobalTime();
  const double Edep=theStep->GetTotalEnergyDeposit();
  const double xPos=theStep->GetTrack()->GetPosition().x();
  if ( Edep <=0. ) return true;

  const int currentParticle = G4TankSteppingAction::GetCurrentParticleId();

  if (time >= 1*utl::second )
    return true;

  //---------------------------------------------------------------------------


  if (  theStep->GetTrack()->GetTrackID() ==1 )
    {
      G4TankTrackingAction::ASCII_Track_p += (theStep->GetStepLength());
      G4TankTrackingAction::ASCII_Edep_p += Edep;
    }

  G4TankTrackingAction::ASCII_Track += (theStep->GetStepLength());
  G4TankTrackingAction::ASCII_Edep  += (Edep);
  G4TankTrackingAction::ASCII_Edep_i[fBarIndex]  += Edep;
  if ( IsFromMuonDecay ) G4TankTrackingAction::ASCII_Edep_Michel  += (Edep);

  //---------------------------------------------------------------------------

  sevt::StationConstants::SignalComponent component;
  sevt::StationConstants::SignalComponent totcomponent=sevt::StationConstants::eTotal;
  const int particle = G4TankEventAction::GetPrimaryId();
  sevt::Station& station = *fG4TankSimulator->GetCurrentEventStationIt();

  if (!station.HasPMT(fPMTIndex)) {
    // if the station doesn't have a scintillator PMT, no PE traces can be stored
    return true;
  }

  if (!station.GetPMT(fPMTIndex).HasSimData())
    station.GetPMT(fPMTIndex).MakeSimData();

  PMTSimData& pmtSim = station.GetPMT(fPMTIndex).GetSimData();

  ostringstream msg;

  switch (particle) {
  case OFFLINE_PHOTON:
    component = sevt::StationConstants::ePhoton;
    break;
  case OFFLINE_ELECTRONS:
    component = sevt::StationConstants::eElectron;
    break;
  case OFFLINE_MUONS:
    component = sevt::StationConstants::eMuon;
    break;
  default:
    component=sevt::StationConstants::eHadron;
    break;
  }

  if (!pmtSim.HasPETimeDistribution(totcomponent)) pmtSim.MakePETimeDistribution(totcomponent);
  if (!pmtSim.HasPETimeDistribution(component)) pmtSim.MakePETimeDistribution(component);

  TimeDistributionI &totalDist  =  pmtSim.GetPETimeDistribution(totcomponent);
  TimeDistributionI &compDist   =  pmtSim.GetPETimeDistribution(component);

  const double attenuationLength=308.*CLHEP::cm,attenuationAmplitude=24.4;
  const double attenuationReferenceEnergy=2.25*CLHEP::MeV;


  double X=xPos+900.*CLHEP::mm+20.*CLHEP::cm; // X from the PMT
  double lambda=attenuationAmplitude*exp(-X/attenuationLength);

  double nPe_mean=lambda*Edep/attenuationReferenceEnergy;
  unsigned int nPe=RandPoisson::shoot(&fRandomEngine->GetEngine(), nPe_mean);

  const double refractionIndex=1.42;
  const double decayTime[2]={ 3.7, 3.5 } ;// fiber/scin [ns]
  const double pathDelay  = (xPos/CLHEP::m) * refractionIndex/ utl::kSpeedOfLight;


  // Loop over SPE
  for (unsigned int ipe=0;ipe<nPe;ipe++)
    {
      double delay=pathDelay;
      const double fiberDelay = RandExponential::shoot(&fRandomEngine->GetEngine(), decayTime[0]);
      const double scintDelay = RandExponential::shoot(&fRandomEngine->GetEngine(), decayTime[1]);
      delay += ( fiberDelay+scintDelay);

      totalDist.AddTime(time+delay);
      compDist.AddTime(time+delay);
      if ( IsFromMuonDecay )
        {
          G4TankTrackingAction::MichelElectronNPE[0]+=1;
          G4TankTrackingAction::MichelElectronNPE_time[0]+=(time+delay);
        }
    }

  return true;

}