MdSimulation/Deprecated/MdGroundPropagatorAG/PrimaryGenerator.cc

Go to the documentation of this file.

#include "PrimaryGenerator.h"
#include "MdGroundPropagator.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>

using namespace GroundPropagatorAG;
using namespace fwk;
using namespace det;
using namespace mevt;
using utl::CoordinateSystem;
using utl::CoordinateSystemPtr;
using utl::Point;
using utl::Particle;
using utl::ReferenceEllipsoid;
using utl::AugerCoordinateSystem;
using std::cout;
using std::endl;


PrimaryGenerator::PrimaryGenerator() :
  fParticleGun(new G4ParticleGun(1)),
  fParticleTable(G4ParticleTable::GetParticleTable())
{ }


PrimaryGenerator::PrimaryGenerator(const std::string&
#ifdef AMIGA_DEBUG
  statf
#endif
) :
  fParticleGun(new G4ParticleGun(1)),
  fParticleTable(G4ParticleTable::GetParticleTable())
{
#ifdef AMIGA_DEBUG
  fStatFile.open(statf.c_str(), std::ios_base::app);
  fStatFile << "#Id PartName X Y Z(km,SiteCoordSys) x y z(m,CounterCoordSys) Time(ns) Parent\n";
#endif
}


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


void
PrimaryGenerator::GeneratePrimaries(G4Event* const theEvent)
{
  const mdet::Counter* theCurrentDetectorCounter                 = GroundPropagator::GetCurrentDetectorCounter();
  mevt::CounterSimData::GrdParticleIterator theCurrentParticleIt = GroundPropagator::GetCurrentParticleIt();

  //Creates coordinate system locally at counter position
  Point counterPos    = theCurrentDetectorCounter->GetPosition();
  //const int counterId = theCurrentDetectorCounter->GetId();
  const utl::CoordinateSystemPtr localCS = theCurrentDetectorCounter->GetLocalCoordinateSystem();
  CoordinateSystemPtr siteCS  = det::Detector::GetInstance().GetSiteCoordinateSystem();
/*
  const ReferenceEllipsoid& e = ReferenceEllipsoid::GetWGS84();
  CoordinateSystemPtr localCS = AugerCoordinateSystem::Create(counterPos, e, siteCS);
*/
  //Calculates the centroid of the modules that build up the counter to fix the origin there
  double xc(0), yc(0), areaTot(0);
  areaTot = 0;
  xc = yc = 0;
  typedef mdet::Counter::ModuleConstIterator ModuleIt;
  //typedef mdet::Module::ScintillatorConstIterator ScintiIt;
  for( ModuleIt mIt = theCurrentDetectorCounter->ModulesBegin(); mIt != theCurrentDetectorCounter->ModulesEnd(); ++mIt ) {
    Point mpos = mIt->GetPosition();

    double area = mIt->GetArea();

    areaTot += area;

    xc += mpos.GetX(theCurrentDetectorCounter->GetLocalCoordinateSystem()) * area;
    yc += mpos.GetY(theCurrentDetectorCounter->GetLocalCoordinateSystem()) * area;
    //G4cerr << "Module position ";
    //G4cerr << "x " << mpos.GetX(theCurrentDetectorCounter->GetLocalCoordinateSystem()) / utl::m << "m, ";
    //G4cerr << "y " << mpos.GetY(theCurrentDetectorCounter->GetLocalCoordinateSystem()) / utl::m << "m\n";
    //G4cerr << "area " << area;
  }
  xc /= areaTot;
  yc /= areaTot;
  //G4cerr << "Centroid position ";
  //G4cerr << "xc " << xc / utl::m << "m, ";
  //G4cerr << "yc " << yc / utl::m << "m\n";

  //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() - counterPos).GetCoordinates(localCS);
  x += xc;
  y += yc;
  //force z to be on the ground
  z = 0;

  double xp(0), yp(0), zp(0);
  boost::tie(xp, yp, zp) = (theCurrentParticleIt->GetPosition()).GetCoordinates(siteCS);
  //force z to be on the ground
  zp = 0;

#ifdef AMIGA_DEBUG
  fStatFile
    << counterId << " "
    << theCurrentParticleIt->GetName()
    << " " << xp/ utl::km
    << " " << yp/ utl::km
    << " " << zp/ utl::km
    << " " << x/ utl::m
    << " " << y/ utl::m
    << " " << z/ utl::m
    << " " << theCurrentParticleIt->GetTime().GetInterval()/utl::ns
    << " " << theCurrentParticleIt->GetParent().GetType()
    << G4endl;
 #endif

  // 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(csDir),
                          theCurrentParticleIt->GetDirection().GetY(csDir),
                          theCurrentParticleIt->GetDirection().GetZ(csDir));

  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);
}