GtRockGammaTool.cc

Go to the documentation of this file.

#include "GtRockGammaTool.h"

#include "GaudiKernel/IRndmGenSvc.h"
#include "GaudiKernel/IParticlePropertySvc.h"
#include "GaudiKernel/ParticleProperty.h"
#include "GaudiKernel/Vector3DTypes.h"
#include "GaudiKernel/PhysicalConstants.h"
#include "HepMC/GenEvent.h"
#include "HepMC/GenParticle.h"
#include "HepMC/GenVertex.h"

#include "CLHEP/Units/SystemOfUnits.h"
#include "CLHEP/Units/PhysicalConstants.h"

#include "DetHelpers/IPositionerTool.h"

using namespace CLHEP;

GtSourceRect::GtSourceRect(const CLHEP::Hep3Vector& origin,
                           const CLHEP::Hep3Vector& edgeA,
                           const CLHEP::Hep3Vector& edgeB)
  : m_origin(origin),
    m_edgeA(edgeA),
    m_edgeB(edgeB)
{
  // Pre-calculate area and normal vector
  m_area = edgeA.mag()*edgeB.mag();
  m_normal = edgeA.cross(edgeB).unit();
}

GtSourceRect::~GtSourceRect()
{
  ;
}

double GtSourceRect::area()
{
  return m_area;
}

CLHEP::Hep3Vector GtSourceRect::origin()
{
  return m_origin;
}

CLHEP::Hep3Vector GtSourceRect::normal()
{
  return m_normal;
}

CLHEP::Hep3Vector GtSourceRect::randomPoint(Rndm::Numbers& random)
{
  // Choose a random point on the rectangular surface
  double x=m_origin.x();
  double y=m_origin.y();
  double z=m_origin.z();
  double randA = random()-0.5;
  double randB = random()-0.5;
  x+=randA*m_edgeA.x();
  y+=randA*m_edgeA.y();
  z+=randA*m_edgeA.z();
  x+=randB*m_edgeB.x();
  y+=randB*m_edgeB.y();
  z+=randB*m_edgeB.z();
  return CLHEP::Hep3Vector(x,y,z);
}

GtSourceOct::GtSourceOct(const CLHEP::Hep3Vector& origin,
                         const CLHEP::Hep3Vector& edgeA,
                         const CLHEP::Hep3Vector& edgeB,
                         const double& bevelLength)
  : GtSourceRect(origin,edgeA,edgeB),
    m_bevelLength(bevelLength)
{
  // Pre-calculate area
  m_area = edgeA.mag()*edgeB.mag() - 2*bevelLength*bevelLength;
}

GtSourceOct::~GtSourceOct()
{
  ;
}

CLHEP::Hep3Vector GtSourceOct::randomPoint(Rndm::Numbers& random)
{
  // Choose a random point on the octagonal surface
  static const double sqrt2 = sqrt(2);
  double cutDist = ((m_edgeA.mag()+m_edgeB.mag())/(2*sqrt2)) 
    - m_bevelLength/sqrt2;
  CLHEP::Hep3Vector point;
  while(true){
    point = GtSourceRect::randomPoint(random);
    // Check if point is in beveled corners
    double pointDist = (fabs(point.x()) + fabs(point.y())) / sqrt2;
    if( pointDist < cutDist ) break;
  }
  return point;
}


GtRockGammaTool::GtRockGammaTool(const std::string& type,
                                 const std::string& name, 
                                 const IInterface* parent) 
  :GaudiTool(type,name,parent),
   m_totalSurfaceArea(0),
   m_walls()
{
  // Initialization
  m_pdfGammas.clear();   
  m_pdfGammas.push_back(1);
  m_edgesGammas.clear(); 
  m_edgesGammas.push_back(0.5*MeV); 
  m_edgesGammas.push_back(10*MeV);

  declareInterface<IHepMCEventMutator>(this);
  declareProperty("GammaSpectrum",m_pdfGammas,"User-Defined PDF for Gamma Spectrum");
  declareProperty("GammaBinEdges",m_edgesGammas,"Bin Edges of Gamma Spectrum");

}

GtRockGammaTool::~GtRockGammaTool() 
{
  for(unsigned int wallIdx=0;wallIdx<m_walls.size();wallIdx++){
    delete m_walls[wallIdx];
    m_walls[wallIdx] = 0;
  }
  m_walls.clear();
}

StatusCode GtRockGammaTool::initialize() 
{
  // Initialize the rock gamma tool
  // Random number service
  IRndmGenSvc *rgs = 0;
  if (service("RndmGenSvc",rgs,true).isFailure()) {
    fatal() << "Failed to get random service" << endreq;
    return StatusCode::FAILURE;        
  }
  
  StatusCode sc;
  sc = m_uni.initialize(rgs, Rndm::Flat(0,1));
  if (sc.isFailure()) {
    fatal() << "Failed to initialize uniform random numbers" << endreq;
    return StatusCode::FAILURE;
  }

  // check user input for user-defined PDFs
  if (m_pdfGammas.size()+1 != m_edgesGammas.size() ) {
    fatal() << "There should be 1 more number of bin edges than bins when defining pdf" << endreq;
    return StatusCode::FAILURE; 
  }

  // Fill pre-calculated Cumulative Distribution Functions
  m_cdfGammas.push_back( 0 );
  for (unsigned int bin=0; bin<m_pdfGammas.size(); bin++){
    m_cdfGammas.push_back( m_cdfGammas[bin] + m_pdfGammas[bin] );
  }
  // Normalize
  for (unsigned int bin=0; bin<m_cdfGammas.size(); bin++){
    m_cdfGammas[bin] /= m_cdfGammas.back();
  }

  // Pre-calculate some useful quantities 

  // FIXME: Current surfaces are defined as an octagon consistent with
  // the water pool dead volume, but the top surface of the octagon
  // has been shifted vertically to encompass most of the RPC.  It
  // would be better to define two volumes (Pool+RPC) to completely
  // encompass the RPC.

  // This defines 10 walls for the octagonal space:
  //        +Y
  // -X+Y /----\ +X+Y
  //  -X |      | +X    and (+Z, -Z)
  // -X-Y \----/ +X-Y
  //        -Y

  // FIXME: should extract these numbers from DetDesc Geometry
  double poolDeadSizeZ=10.0*m;
  double nearPoolDeadSizeX=16.0*m;
  double nearPoolDeadSizeY=10.0*m;
  //double farPoolDeadSizeX=16.0*m;
  //double farPoolDeadSizeY=16.0*m;
  // Bevel size is the distance that the octagonal corner walls remove
  // from the initial rectangular pool shape.  It is twice the
  // distance from the bevel wall to the corner that would be formed
  // by the X,Y walls.
  double poolDeadBevelSize=4.249*m;
  // RPC space above detector
  // Note: I have pre-rotated the RPC XY coords by 90deg to match the pool
  double rpcRoofSpace=0.50*m;
  //double rpcNearRoofSizeX=18.5*m;
  //double rpcNearRoofSizeY=12.5*m;
  double rpcNearRoofSizeZ=0.30*m;
  double rpcNearRoofHeightZ= rpcRoofSpace + rpcNearRoofSizeZ;
  //double rpcFarRoofSizeX=18.5*m;
  //double rpcFarRoofSizeY=18.5*m;
  //double rpcFarRoofSizeZ=0.30*m;
  //double rpcFarRoofHeightZ= rpcRoofSpace + rpcFarRoofSizeZ;

  // Calculate some helper parameters
  // FIXME: Hardcoded near side walls.  Must change for Far site.
  double poolAndRpcZ = poolDeadSizeZ + rpcNearRoofHeightZ;
  double sqrt2 = sqrt(2);
  // Extra RPC space shifts wall center higher by half the RPC thickness
  double zOrigin = rpcNearRoofHeightZ / 2.;
  // The bevelXYOffset is the distance along the X or Y axis removed by the bevel.
  double bevelXYOffset = poolDeadBevelSize/sqrt2;

  m_walls.clear();
  {
    // +X wall
    CLHEP::Hep3Vector origin(nearPoolDeadSizeX/2.,0,zOrigin);
    CLHEP::Hep3Vector edgeA(0,
                            0,
                            poolAndRpcZ);
    CLHEP::Hep3Vector edgeB(0,
                            (nearPoolDeadSizeY-2*bevelXYOffset),
                            0);
    m_walls.push_back( new GtSourceRect(origin,edgeA,edgeB) );
  }
  {
    // +X+Y wall
    CLHEP::Hep3Vector origin(nearPoolDeadSizeX/2.-bevelXYOffset/2.,
                             nearPoolDeadSizeY/2.-bevelXYOffset/2.,
                             zOrigin);
    CLHEP::Hep3Vector edgeA(0,
                            0,
                            poolAndRpcZ);
    CLHEP::Hep3Vector edgeB(-bevelXYOffset,
                            bevelXYOffset,
                            0);
    m_walls.push_back( new GtSourceRect(origin,edgeA,edgeB) );
  }
  {
    // +Y wall
    CLHEP::Hep3Vector origin(0,nearPoolDeadSizeY/2.,zOrigin);
    CLHEP::Hep3Vector edgeA(0,
                            0,
                            poolAndRpcZ);
    CLHEP::Hep3Vector edgeB(-(nearPoolDeadSizeX-2*bevelXYOffset),
                            0,
                            0);
    m_walls.push_back( new GtSourceRect(origin,edgeA,edgeB) );
  }
  {
    // -X+Y wall
    CLHEP::Hep3Vector origin(-(nearPoolDeadSizeX/2.-bevelXYOffset/2.),
                             nearPoolDeadSizeY/2.-bevelXYOffset/2.,
                             zOrigin);
    CLHEP::Hep3Vector edgeA(0,
                            0,
                            poolAndRpcZ);
    CLHEP::Hep3Vector edgeB(-bevelXYOffset,
                            -bevelXYOffset,
                            0);
    m_walls.push_back( new GtSourceRect(origin,edgeA,edgeB) );
  }
  {
    // -X wall
    CLHEP::Hep3Vector origin(-nearPoolDeadSizeX/2.,0,zOrigin);
    CLHEP::Hep3Vector edgeA(0,
                            0,
                            poolAndRpcZ);
    CLHEP::Hep3Vector edgeB(0,
                            -(nearPoolDeadSizeY-2*bevelXYOffset),
                            0);
    m_walls.push_back( new GtSourceRect(origin,edgeA,edgeB) );
  }
  {
    // -X-Y wall
    CLHEP::Hep3Vector origin(-(nearPoolDeadSizeX/2.-bevelXYOffset/2.),
                             -(nearPoolDeadSizeY/2.-bevelXYOffset/2.),
                             zOrigin);
    CLHEP::Hep3Vector edgeA(0,
                            0,
                            poolAndRpcZ);
    CLHEP::Hep3Vector edgeB(bevelXYOffset,
                            -bevelXYOffset,
                            0);
    m_walls.push_back( new GtSourceRect(origin,edgeA,edgeB) );
  }
  {
    // -Y wall
    CLHEP::Hep3Vector origin(0,-nearPoolDeadSizeY/2.,zOrigin);
    CLHEP::Hep3Vector edgeA(0,
                            0,
                            poolAndRpcZ);
    CLHEP::Hep3Vector edgeB(nearPoolDeadSizeX-2*bevelXYOffset,
                            0,
                            0);
    m_walls.push_back( new GtSourceRect(origin,edgeA,edgeB) );
  }
  {
    // +X-Y wall
    CLHEP::Hep3Vector origin(nearPoolDeadSizeX/2.-bevelXYOffset/2.,
                             -(nearPoolDeadSizeY/2.-bevelXYOffset/2.),
                             zOrigin);
    CLHEP::Hep3Vector edgeA(0,
                            0,
                            poolAndRpcZ);
    CLHEP::Hep3Vector edgeB(bevelXYOffset,
                            bevelXYOffset,
                            0);
    m_walls.push_back( new GtSourceRect(origin,edgeA,edgeB) );
  }
  {
    // +Z wall
    CLHEP::Hep3Vector origin(0,0,poolDeadSizeZ/2. + rpcNearRoofHeightZ);
    CLHEP::Hep3Vector edgeA(0,
                            nearPoolDeadSizeY,
                            0);
    CLHEP::Hep3Vector edgeB(nearPoolDeadSizeX,
                            0,
                            0);
    m_walls.push_back( new GtSourceOct(origin,edgeA,edgeB,bevelXYOffset) );
  }
  {
    // -Z wall
    CLHEP::Hep3Vector origin(0,0,-poolDeadSizeZ/2.);
    CLHEP::Hep3Vector edgeA(nearPoolDeadSizeX,
                            0,
                            0);
    CLHEP::Hep3Vector edgeB(0,
                            nearPoolDeadSizeY,
                            0);
    m_walls.push_back( new GtSourceOct(origin,edgeA,edgeB,bevelXYOffset) );
  }

  m_totalSurfaceArea = 0;
  for(unsigned int wallIdx=0;wallIdx<m_walls.size();wallIdx++){
    debug() << "Adding wall with area: " << m_walls[wallIdx]->area() << endreq;
    m_totalSurfaceArea += m_walls[wallIdx]->area();
  }
  debug() << "Total surface area: " << m_totalSurfaceArea << endreq;

  return this->GaudiTool::initialize();
}

StatusCode GtRockGammaTool::finalize()
{
  return this->GaudiTool::finalize();
}

StatusCode GtRockGammaTool::mutate(HepMC::GenEvent& event)
{ 
  if (this->oneVertex(event).isFailure()) 
    return StatusCode::FAILURE;
  return StatusCode::SUCCESS;
}

StatusCode GtRockGammaTool::oneVertex(HepMC::GenEvent& event)
{
  // generate vertex time (within single event)
  double vertex_t = 0;

  // generate vertex position in local pool coordinates
  // (random point around pool edges, top, or bottom)
  unsigned int wallIdx = GetRandomWall();
  // Choose a random point on the wall
  CLHEP::Hep3Vector point = m_walls[wallIdx]->randomPoint(m_uni);
  HepMC::GenVertex* vertex = new HepMC::GenVertex(HepMC::FourVector(point.x(),
                                                                    point.y(),
                                                                    point.z(),
                                                                    vertex_t));
  if( !event.add_vertex(vertex) ){
    error() << "Failed to add vertex to event" << endreq;
    return StatusCode::FAILURE;
  }

  // generate momentum direction
  CLHEP::Hep3Vector dir;
  double dir_costh = m_uni(); // cos(theta) range 0 to 1
  double dir_phi = m_uni()*twopi; // range from 0 to 2*pi radians
  
  dir.setRThetaPhi(1.0,acos(dir_costh),dir_phi);

  // rotate along Y-axis by vertex direction's theta
  dir.rotateY(m_walls[wallIdx]->normal().theta());
  // rotate along Z-axis by vertex direction's phi
  dir.rotateZ(m_walls[wallIdx]->normal().phi());
  // Normalize direction vector
  dir = dir.unit();

  // generate random energy
  double energy = ConvertCdfRand(m_uni(),m_cdfGammas,m_edgesGammas); 
  HepMC::FourVector fourmom(energy*dir.x(),
                            energy*dir.y(),
                            energy*dir.z(),
                            energy);
  
  HepMC::GenParticle* particle = new HepMC::GenParticle(fourmom,22,
                                                        1/*=status*/);
  // Set photon polarization
  double pol_phi = m_uni()*CLHEP::twopi;
  particle->set_polarization(HepMC::Polarization(0.5*pi,pol_phi));
  // Add gamma to vertex
  vertex->add_particle_out(particle);

  return StatusCode::SUCCESS;
}

double GtRockGammaTool::ConvertCdfRand(const double& rand, 
                                       const std::vector<double>& cdf, 
                                       const std::vector<double>& edges) {
  // Convert a uniform random number in [0,1] to a random sample from
  // the given cumulative distribution function.
  int left = 0;
  int right = cdf.size();
  const double* cdfStart = &cdf[0];
  const double* edgesStart = &edges[0];

  while( right-left > 1 ){
    int mid = (right+left)/2;
    if( rand < *(cdfStart+mid) ) right = mid;
    else left = mid;
  }
  double value = *(edgesStart + left);
  double slope = (*(edgesStart+right) - *(edgesStart+left))
                /(*(cdfStart+right) - *(cdfStart+left));
  value += slope * (rand - *(cdfStart+left));
  return value;
}

unsigned int GtRockGammaTool::GetRandomWall() {
  // Choose a random wall, based on surface area
  double rand = m_uni();
  double randArea = rand*m_totalSurfaceArea;
  unsigned int wallIdx = 0;
  while(wallIdx<m_walls.size()){
    if( randArea < m_walls[wallIdx]->area() ){
      // Found the correct wall
      break;
    }
    randArea -= m_walls[wallIdx]->area();
    wallIdx++;
  }
  return wallIdx;
}