MdSimulation/Deprecated/MdEdepSimulatorAG/PrimaryGenerator.cc

Go to the documentation of this file.

#include "PrimaryGenerator.h"
#include "MdEdepSimulator.h"

#include <det/Detector.h>
#include <mevt/CounterSimData.h>

#include <mdet/Counter.h>

#include <utl/CoordinateSystem.h>
#include <utl/CoordinateSystemPtr.h>
#include <utl/Point.h>
#include <utl/Particle.h>
#include <utl/ReferenceEllipsoid.h>
#include <utl/AugerCoordinateSystem.h>

#include <cstddef>
#include <iostream>
#include <sstream>

#include <G4Event.hh>
#include <G4ParticleGun.hh>
#include <G4ParticleDefinition.hh>
#include <G4ParticleTable.hh>

#include <utl/AugerUnits.h>

#include <fwk/CentralConfig.h>
#include <utl/ErrorLogger.h>
#include <utl/Reader.h>
#include <utl/MathConstants.h>

using namespace EdepSimulatorAG;
using namespace fwk;
using namespace utl;
using namespace det;
using namespace mevt;
using std::cout;
using std::endl;

PrimaryGenerator::PrimaryGenerator() : 
  fParticleGun(new G4ParticleGun(1)), 
  fParticleTable(G4ParticleTable::GetParticleTable()) 
{
  CentralConfig* theConfig = CentralConfig::GetInstance(); 
  Branch topB = theConfig->GetTopBranch("MdEdepSimulator");
  std::vector<G4double> vScintiSize;
  Branch scintillatorB = topB.GetChild("scintillator");
  scintillatorB.GetChild("size").GetData(vScintiSize);
  fScint_Side_Z = vScintiSize.at(2)/2;
}

PrimaryGenerator::~PrimaryGenerator() 
{    
  delete fParticleGun;
}

void
PrimaryGenerator::GeneratePrimaries(G4Event* theEvent) 
{
  const mdet::Scintillator* theCurrentDetectorScintillator         = EdepSimulator::GetCurrentDetectorScintillator();
  mevt::ScintillatorSimData::ParticleIterator theCurrentParticleIt = EdepSimulator::GetCurrentParticleIt();
 
  //Creates coordinate system locally at scint position
  Point scintPos    = theCurrentDetectorScintillator->GetPosition();

  //Scintillator local coordinate system (particles are loaded in this CS by the MdUnderGroundInjector)
  CoordinateSystemPtr localCS = theCurrentDetectorScintillator->GetLocalCoordinateSystem();
  //Retrives particle coordinate system
  CoordinateSystemPtr csDir   = theCurrentParticleIt->GetDirection().GetCoordinateSystem();
  
  // set properties and pass to particle gun
  G4ParticleDefinition* particleDef = fParticleTable->FindParticle(theCurrentParticleIt->GetName());

  if (!particleDef) {
    std::ostringstream msg;
    msg << "Undefined particle type: "
        << theCurrentParticleIt->GetName();
    WARNING(msg);
    return;
  }

  double x(0), y(0), z(0);
  boost::tie(x, y, z) = (theCurrentParticleIt->GetPosition() - scintPos).GetCoordinates(localCS);    
  //force z to be on scintillator top  
  z = fScint_Side_Z;

 G4cerr
  << "Scint ID = " << theCurrentDetectorScintillator->GetId()
  << ", part. type = " << theCurrentParticleIt->GetName() 
  << ", x = " << x/utl::cm << " (cm), " 
  << "y = " << y/utl::cm << " (cm), " 
  << "z = " << z/utl::cm << " (cm), " 
  << "t = " << theCurrentParticleIt->GetTime().GetInterval()/utl::ns << " (ns)"
  << G4endl;

  // Convert from Auger units to G4 units
  G4ThreeVector position(x * CLHEP::m / utl::m, y * CLHEP::m / utl::m, z * CLHEP::m / utl::m);
  
  G4ThreeVector direction(theCurrentParticleIt->GetDirection().GetX(localCS),  // JMF 01-Oct-2014 Directions should be loaded in scint CS
                          theCurrentParticleIt->GetDirection().GetY(localCS),  // JMF 01-Oct-2014
                          theCurrentParticleIt->GetDirection().GetZ(localCS)); // JMF 01-Oct-2014

  fParticleGun->SetParticleDefinition(particleDef);
  fParticleGun->SetParticlePosition(position);
  fParticleGun->SetParticleMomentumDirection(direction);

  // Convert from Auger units to G4 units
  fParticleGun->SetParticleEnergy(theCurrentParticleIt->GetKineticEnergy() * CLHEP::eV/utl::eV);
  fParticleGun->SetParticleTime(theCurrentParticleIt->GetTime().GetInterval() * CLHEP::ns/utl::ns);
    
  fParticleGun->GeneratePrimaryVertex(theEvent);
}