In This Package:

GtDiffuserBallTool.cc

Go to the documentation of this file.

#include "GtDiffuserBallTool.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 "CLHEP/Vector/ThreeVector.h"

#include "DetHelpers/IPositionerTool.h"

using namespace CLHEP;

GtDiffuserBallTool::GtDiffuserBallTool(const std::string& type,
                                       const std::string& name, 
                                       const IInterface* parent) 
  :GaudiTool(type,name,parent)
{
  // Initialization
  m_pid = 0;
  m_pdfWave.clear();   m_pdfWave.push_back(1);
  m_edgesWave.clear(); m_edgesWave.push_back(429*nm); m_edgesWave.push_back(431*nm);
  m_pdfTime.clear();   m_pdfTime.push_back(1);
  m_edgesTime.clear(); m_edgesTime.push_back(0*ns);   m_edgesTime.push_back(1*ns);

  declareInterface<IHepMCEventMutator>(this);
  declareProperty("PhotonsPerEvent",m_particlesPerEvent=3500,"Number of photons per event");
  declareProperty("PhotonsPerEventMode",m_particlesPerEventMode="Fixed","Type of distribution of number of photons per event");
  declareProperty("PhotonsPerEventSpread",m_particlesPerEventSpread=0,"Sigma/FWHM of number of photons per event");

  declareProperty("WavelengthMode",m_wavelengthMode="Fixed","Wavelength Mode");
  declareProperty("Wavelength", m_wavelength=430.0*nm,"Fixed (Mean) Photon Wavelength");
  declareProperty("WavelengthSpread",m_wavelengthSpread=5.0*nm,"Wavelength Spectrum Spread");
  declareProperty("PdfWave",m_pdfWave,"User-Defined PDF for Wavelength Spectrum");
  declareProperty("PdfEdgesWave",m_edgesWave,"Bin Edges of PdfWave");

  declareProperty("TimingMode",m_timingMode="Instant","Timing Mode");
  declareProperty("PdfTime",m_pdfTime,"User-Defined PDF for Timing Distribution");
  declareProperty("PdfEdgesTime",m_edgesTime,"Bin Edges of PdfTime"); 

  declareProperty("Radius",m_radius=0.9525*cm,"Radius of Diffuser Ball");
  declareProperty("Anisotropy",m_anisotropy=true,"Turn on/off anisotropy");
  declareProperty("BrightSpotTheta",m_bright_theta=0.0,"Polar angle of bright spot");
  declareProperty("BrightSpotPhi",m_bright_phi=0.0,"Azimuthal angle of bright spot");
  declareProperty("PhotonScaleWeight",m_photonScaleWeight = 3.74,
                  "Reweight photons based on maximum quantum efficiency");
  declareProperty("GeomPosTool",m_geomPosToolName="diffuserBallGeomPos",
          "Name of tool used to dynamically place the diffuser ball geometry");

}

GtDiffuserBallTool::~GtDiffuserBallTool() 
{
}

StatusCode GtDiffuserBallTool::initialize() 
{
  // Particle identification
  IParticlePropertySvc * ppSvc = 
    svc< IParticlePropertySvc >( "ParticlePropertySvc" , true ) ;
  ParticleProperty* particle = 0;

  m_particleName = "opticalphoton";
  particle = ppSvc->find(m_particleName);
  if(particle) {
    m_pid = particle->pdgID();
  }

  // if failed to find particle, return failure
  if (!particle) {
    fatal() << "Failed to find particle named \"" 
            << m_particleName << "\" and with ID " 
            << m_pid << endreq;
    return StatusCode::FAILURE;
  }

  // check modes
  if ((m_timingMode != "Instant") && (m_timingMode != "Defined")) {
    fatal() << "Not a recognized vertex timing distribution mode" << endreq;
    return StatusCode::FAILURE;
  }
  if ((m_wavelengthMode != "Fixed")   && (m_wavelengthMode != "Uniform") &&
      (m_wavelengthMode != "Smeared") && (m_wavelengthMode != "Defined")) {
    fatal() << "Not a recognized wavelength spectrum mode" << endreq;
    return StatusCode::FAILURE;
  }

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

  // random numbers for particle number
  if (m_particlesPerEventMode == "Gaus") {
    sc = m_randPartNo.initialize(rgs, Rndm::Gauss(0,1));
    if (sc.isFailure()) {
      fatal() << "Failed to initialize random numbers for particle number distribution" << endreq;
      return StatusCode::FAILURE;
    }
  }
  else if (m_particlesPerEventMode == "Lorentz"){
    sc = m_randPartNo.initialize(rgs, Rndm::BreitWigner(0,1));
    if (sc.isFailure()) {
      fatal() << "Failed to initialize random numbers for particle number distribution" << endreq;
      return StatusCode::FAILURE;
    }
  }
  
  // random numbers for timing
  if (m_timingMode == "Defined") {
    sc = m_randTime.initialize(rgs, Rndm::DefinedPdf(m_pdfTime,0));
    if (sc.isFailure()) {
      fatal() << "Failed to initialize random numbers for timing distribution" << endreq;
      return StatusCode::FAILURE;
    }
  }
  
  // random numbers for wavelength spectrum
  if (m_wavelengthMode == "Uniform") {
    sc = m_randWave.initialize(rgs, Rndm::Flat(0,1));
    if (sc.isFailure()) {
      fatal() << "Failed to initialize random numbers for wavelength spectrum" << endreq;
      return StatusCode::FAILURE;
    }
  }
  else if (m_wavelengthMode == "Smeared") {
    sc = m_randWave.initialize(rgs, Rndm::Gauss(0,1));
    if (sc.isFailure()) {
      fatal() << "Failed to initialize random numbers for wavelength spectrum" << endreq;
      return StatusCode::FAILURE;
    }
  }
  else if (m_wavelengthMode == "Defined") {
    sc = m_randWave.initialize(rgs, Rndm::DefinedPdf(m_pdfWave,0));
    if (sc.isFailure()) {
      fatal() << "Failed to initialize random numbers for wavelength spectrum" << endreq;
      return StatusCode::FAILURE;
    }
  }

  // check user input for user-defined PDFs
  if (m_pdfWave.size()+1 != m_edgesWave.size() || 
      m_pdfTime.size()+1 != m_edgesTime.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
  if(m_pdfTime.size()>0){
    m_cdfTime.push_back( 0 );
    for (unsigned int bin=0; bin<m_pdfTime.size(); bin++)
      m_cdfTime.push_back( m_cdfTime[bin] + m_pdfTime[bin] );
    // Normalize
    for (unsigned int bin=0; bin<m_cdfTime.size(); bin++)
      m_cdfTime[bin] /= m_cdfTime.back();
  }
  if(m_pdfWave.size()>0){
    m_cdfWave.push_back( 0 );
    for (unsigned int bin=0; bin<m_pdfWave.size(); bin++)
      m_cdfWave.push_back( m_cdfWave[bin] + m_pdfWave[bin] );
    // Normalize
    for (unsigned int bin=0; bin<m_cdfWave.size(); bin++)
      m_cdfWave[bin] /= m_cdfWave.back();
  }
  
  if( m_photonScaleWeight > 1.0 ){
    // Reduce number of photons, and scale each photon weight to compensate
    info() << "Scaling each photon weight by " << m_photonScaleWeight 
           << " and reducing the number of photons from "
           << m_particlesPerEvent << " to " 
           << int( m_particlesPerEvent / m_photonScaleWeight )
           << endreq;
    m_particlesPerEvent = int( m_particlesPerEvent / m_photonScaleWeight );
  }

  if( m_geomPosToolName == "" ){
    info() << "Placement of Diffuser Ball geometry not requested.  "
           << "Optical photons will still be produced at the requested "
           << "location." << endreq;
  }else{
    IPositionerTool* geomPosTool = 0;
    try{
      geomPosTool = tool<IPositionerTool>(m_geomPosToolName);
    }
    catch(const GaudiException& exg) {
      fatal() << "Failed to get geometry tool: \"" << m_geomPosToolName << "\"" 
              << endreq;
      return StatusCode::FAILURE;
    }
    info () << "Placing diffuser ball geometry using " << m_geomPosToolName 
            << endreq;
    sc = geomPosTool->placeVolume();
    if(sc.isFailure()) return sc;
    geomPosTool->release();
    geomPosTool = 0;
  }

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

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

StatusCode GtDiffuserBallTool::mutate(HepMC::GenEvent& event)
{ 
  verbose() << "Mode: " << m_particlesPerEventMode << endreq;
  int particleNumber = m_particlesPerEvent;
  debug() << " particleNumber " << particleNumber << " mode " << m_particlesPerEventMode << endreq;
  if(m_particlesPerEventMode != "Fixed"){
    // Smear particle number
    particleNumber = m_particlesPerEvent + int(m_particlesPerEventSpread/m_photonScaleWeight * m_randPartNo()); 
    debug() << " smeared # of particles " << particleNumber << endreq;
    // Make sure that at least one dummy photon is generated
    if ( particleNumber < 1 ) { 
      particleNumber = 1; 
      debug() << " after requiring >0 photons. particleNumber is " << particleNumber << endreq;
    }
    else {
      // Cap photon number
      if ( particleNumber > int(m_particlesPerEvent * 4)) { particleNumber = int(m_particlesPerEvent * 4); }   
      debug() << " after cap. particleNumber is " << particleNumber << endreq;
    }
    debug() << "Generate " << particleNumber << " photons" << endreq;
  }
  for (int ind=0; ind<particleNumber; ++ind) {
    if (this->oneVertex(event).isFailure()) 
      return StatusCode::FAILURE;
  }
  return StatusCode::SUCCESS;
}

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

  // generate vertex position (random point on ball from which photon is emitted)
  double vertex_costh;
  double vertex_phi = m_uni()*CLHEP::twopi;  // range from 0 to 2*pi radians uniformly
  
  if (m_anisotropy){
    // create bright spot at top of diffuser ball, theta=0
    // theta=0 is 50% brighter than theta=pi (brightness taken to be linear)
    // see DocDB-1212

    bool accept = false;
    double vertex_theta;
    while (!accept){
      vertex_theta = m_uni()*CLHEP::pi; 
      double test  = m_uni()*1.25; // function below has max of 1.25 
      if (test <= (-vertex_theta/twopi + 1.5)*sin(vertex_theta)) 
        accept = true;
    }

    // rotate vertex to random bright spot location
    CLHEP::Hep3Vector vertex_bright;
    vertex_bright.setRThetaPhi(1.0,vertex_theta,vertex_phi);
    vertex_bright.rotateY(m_bright_theta);
    vertex_bright.rotateZ(m_bright_phi);
    vertex_bright.unit();

    vertex_costh = vertex_bright.cosTheta();
    vertex_phi   = vertex_bright.phi();
  }else{
    vertex_costh = 2*m_uni()-1;  // cos(theta) ranges from -1 to 1 uniformly
  }

  double vertex_sinth = sqrt(1-vertex_costh*vertex_costh);
  double vertex_x = m_radius*vertex_sinth*cos(vertex_phi);
  double vertex_y = m_radius*vertex_sinth*sin(vertex_phi);
  double vertex_z = m_radius*vertex_costh;
  HepMC::GenVertex* vertex = new HepMC::GenVertex(HepMC::FourVector(vertex_x,
                                                                    vertex_y,
                                                                    vertex_z,
                                                                    vertex_t));
  // Catch scaling of photon weight
  if( m_photonScaleWeight > 1.0 ){
    vertex->weights().push_back( m_photonScaleWeight );
  }
  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
                              // do not want photon traveling into diffuser ball
  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(acos(vertex_costh));
  // rotate along Z-axis by vertex direction's phi
  dir.rotateZ(vertex_phi);
  
  dir = dir.unit();
  
  // generate wavelength/momentum/energy
  double momentum = CLHEP::twopi*CLHEP::hbarc / GetWave();
  double energy = momentum;
  HepMC::FourVector fourmom(momentum*dir.x(),
                            momentum*dir.y(),
                            momentum*dir.z(),
                            energy);
  
  HepMC::GenParticle* particle = new HepMC::GenParticle(fourmom,m_pid,1/*=status*/);
  double pol_phi = m_uni()*CLHEP::twopi;
  particle->set_polarization(HepMC::Polarization(0.5*pi,pol_phi));
    
  vertex->add_particle_out(particle);

  return StatusCode::SUCCESS;
}

double GtDiffuserBallTool::GetTime()
{
  double time = -1.0;
  double current_rand = m_randTime();

  if (m_timingMode == "Instant")      time = 0.0;
  else if (m_timingMode == "Defined") time = ConvertCdfRand(current_rand,m_cdfTime,m_edgesTime);

  return time;
}

double GtDiffuserBallTool::GetWave()
{
  double wave = 0.0;
  double current_rand = m_randWave();
  double mean = m_wavelength;
  double spread = m_wavelengthSpread;

  if (m_wavelengthMode == "Fixed")        wave = m_wavelength;
  else if (m_wavelengthMode == "Uniform") wave = (2*current_rand-1)*spread + mean;
  else if (m_wavelengthMode == "Smeared") wave = current_rand*spread + mean;
  else if (m_wavelengthMode == "Defined") wave = ConvertCdfRand(current_rand,m_cdfWave,m_edgesWave); 

  return wave;
}

double GtDiffuserBallTool::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;
}

---