ArtificialShowerFileParticleIterator.cc

Go to the documentation of this file.

#include "ArtificialShowerFileParticleIterator.h"

#include <iostream>
#include <utl/PhysicalConstants.h>
#include <fwk/RandomEngineRegistry.h>
#include <utl/RandomEngine.h>
#include <CLHEP/Random/RandFlat.h>
#include <CLHEP/Random/RandPoisson.h>
#include <CLHEP/Random/RandGauss.h>
#include <fwk/CentralConfig.h>
#include <utl/Reader.h>
#include <utl/config.h>

using namespace std;
using namespace utl;
using namespace io;
using CLHEP::RandFlat;
using CLHEP::RandPoisson;
using namespace fwk;


ArtificialShowerFileParticleIterator::ArtificialShowerFileParticleIterator() :
  fCurrentParticle(utl::Particle::eUndefined,
                   utl::Particle::eShower,
                   Point(0,0,0, CoordinateSystem::GetRootCoordinateSystem()),
                   Vector(0,0,0, CoordinateSystem::GetRootCoordinateSystem()),
                   TimeInterval(0), 0),
  fMaxRadius(4000.0*m),
  fMuonEnergy(2.0*GeV),
  fPhotonEnergy(80.0*MeV),
  fElectronEnergy(80.0*MeV),
  fMuonDensity(0.1/m/m),
  fPhotonDensity(0.1/m/m),
  fElectronDensity(0.1/m/m),
  fParticleZenith(0.0*deg),
  fParticleAzimuth(0.0*deg),
  fParticleTime(2000.0*ns),
  fParticleWeight(100.0),
  fRandomAzimuth(false),
  fGeneratedParticles(0)
{
  ReadXMLParameters();

  fTotalParticles = (fMuonDensity + fPhotonDensity + fElectronDensity) * fMaxRadius * fMaxRadius * kPi;

}


ArtificialShowerFileParticleIterator::ArtificialShowerFileParticleIterator(double maxRadius,
                                                                           double muonEnergy, double photonEnergy, double electronEnergy,
                                                                           double muonDensity, double photonDensity, double electronDensity,
                                                                           double particleZenith, double particleAzimuth, double particleTime,
                                                                           double particleWeight, bool randomAzimuth) :
  fCurrentParticle(utl::Particle::eUndefined,
                   utl::Particle::eShower,
                   Point(0,0,0, CoordinateSystem::GetRootCoordinateSystem()),
                   Vector(0,0,0, CoordinateSystem::GetRootCoordinateSystem()),
                   TimeInterval(0), 0),
  fMaxRadius(maxRadius),
  fMuonEnergy(muonEnergy),
  fPhotonEnergy(photonEnergy),
  fElectronEnergy(electronEnergy),
  fMuonDensity(muonDensity),
  fPhotonDensity(photonDensity),
  fElectronDensity(electronDensity),
  fParticleZenith(particleZenith),
  fParticleAzimuth(particleAzimuth),
  fParticleTime(particleTime),
  fParticleWeight(particleWeight),
  fRandomAzimuth(randomAzimuth),
  fGeneratedParticles(0)
{
  fTotalParticles = (fMuonDensity + fPhotonDensity + fElectronDensity) * fMaxRadius * fMaxRadius * kPi;
}



ArtificialShowerFileParticleIterator::~ArtificialShowerFileParticleIterator()
{
}


void ArtificialShowerFileParticleIterator::SetShowerZenith(const double showerZenith)
{
  fShowerZenith = showerZenith;
}



void ArtificialShowerFileParticleIterator::SetShowerAzimuth(const double showerAzimuth)
{
  fShowerAzimuth = showerAzimuth;
}




void ArtificialShowerFileParticleIterator::ReadXMLParameters()
{
  Branch topB = CentralConfig::GetInstance()->GetTopBranch("ArtificialShowerFile");
  Branch maxRadiusB = topB.GetChild("ShowerRadius");
  maxRadiusB.GetData(fMaxRadius);
  Branch muonEnergyB = topB.GetChild("MuonEnergy");
  muonEnergyB.GetData(fMuonEnergy);
  Branch photonEnergyB = topB.GetChild("PhotonEnergy");
  photonEnergyB.GetData(fPhotonEnergy);
  Branch electronEnergyB = topB.GetChild("ElectronEnergy");
  electronEnergyB.GetData(fElectronEnergy);
  Branch muonDensityB = topB.GetChild("MuonDensity");
  muonDensityB.GetData(fMuonDensity);
  Branch photonDensityB = topB.GetChild("PhotonDensity");
  photonDensityB.GetData(fPhotonDensity);
  Branch electronDensityB = topB.GetChild("ElectronDensity");
  electronDensityB.GetData(fElectronDensity);
  Branch particleZenithB = topB.GetChild("ParticleZenith");
  particleZenithB.GetData(fParticleZenith);
  Branch particleAzimuthB = topB.GetChild("ParticleAzimuth");
  particleAzimuthB.GetData(fParticleAzimuth);
  Branch particleTimeB = topB.GetChild("ParticleTimeWindow");
  particleTimeB.GetData(fParticleTime);
  Branch particleWeightB = topB.GetChild("ParticleWeight");
  particleWeightB.GetData(fParticleWeight);
  Branch randomAzimuthB = topB.GetChild("RandomParticleAzimuth");
  randomAzimuthB.GetData(fRandomAzimuth);
}



void ArtificialShowerFileParticleIterator::GenerateParticle(const CoordinateSystemPtr& groundCS)
{
  fRandomEngine = &RandomEngineRegistry::GetInstance().Get(RandomEngineRegistry::eDetector).GetEngine();


  const double particleType = RandFlat::shoot(fRandomEngine,0.0,1.0);
  const double totalDensity = fMuonDensity + fPhotonDensity + fElectronDensity;
  const double particleTime = RandFlat::shoot(fRandomEngine,0.0,1.0) * fParticleTime;
  double particleX = RandFlat::shoot(fRandomEngine, -fMaxRadius, fMaxRadius);
  double particleY = RandFlat::shoot(fRandomEngine, -fMaxRadius, fMaxRadius);

  while ( particleX*particleX + particleY*particleY > fMaxRadius * fMaxRadius )
    {
      particleX = RandFlat::shoot(fRandomEngine, -fMaxRadius, fMaxRadius);
      particleY = RandFlat::shoot(fRandomEngine, -fMaxRadius, fMaxRadius);
    }

  const double particleAz = fRandomAzimuth ? RandFlat::shoot(fRandomEngine, 0.0, 2.0*kPi) : fParticleAzimuth / rad;

  const double dirX = sin(fParticleZenith / rad) * cos(particleAz);
  const double dirY = sin(fParticleZenith / rad) * sin(particleAz);
  const double dirZ = cos(fParticleZenith / rad) * -1.0;

  double type;
  double energy;


  if (particleType < fMuonDensity / totalDensity)
    {
      type=utl::Particle::eMuon;
      energy=fMuonEnergy;
    }
  else if (particleType < (fPhotonDensity + fMuonDensity) / totalDensity)
    {
      type=utl::Particle::ePhoton;
      energy=fPhotonEnergy;
    }
  else
    {
      type=utl::Particle::eElectron;
      energy=fElectronEnergy;
    }

  fCurrentParticle = Particle(type, utl::Particle::eShower,
                              Point( particleX, particleY, 0.0, groundCS ),
                              Vector( dirX, dirY, dirZ, groundCS ),
                              TimeInterval(particleTime),
                              fParticleWeight,
                              energy);
}


Particle*
ArtificialShowerFileParticleIterator::GetOneParticle(const utl::CoordinateSystemPtr& groundCS)
{

  if (fGeneratedParticles > fTotalParticles)
    return NULL;


  GenerateParticle(groundCS);

  fGeneratedParticles++;
  
  return &fCurrentParticle;
}


void
ArtificialShowerFileParticleIterator::Rewind()
{
  fGeneratedParticles = 0;
}


/*
CoordinateSystemPtr
ArtificialShowerFileParticleIterator::ComputeExternalShowerCoordinateSystem(const utl::CoordinateSystemPtr& cs)
{
  return CoordinateSystem::RotationY(fShowerZenith, CoordinateSystem::RotationZ(fShowerAzimuth, cs));
}
*/