In This Package:

GtPositionerTool.cc

Go to the documentation of this file.

#include "GtPositionerTool.h"

#include "GaudiKernel/IssueSeverity.h"
#include "GaudiKernel/IRndmGenSvc.h"
#include "GaudiKernel/DataObject.h"
#include "GaudiKernel/Point3DTypes.h"


#include "DetDesc/SolidBox.h"
#include "DetDesc/ILVolume.h"
#include "DetDesc/IPVolume.h"
#include "DetDesc/IDetectorElement.h"
#include "DetDesc/IGeometryInfo.h"

#include "HepMC/GenEvent.h"
#include "HepMC/GenVertex.h"
// This may need to change when move to HepMC v2.
#include "CLHEP/Vector/LorentzVector.h"

GtPositionerTool::GtPositionerTool(const std::string& type,
                                   const std::string& name,
                                   const IInterface* parent)
    : GaudiTool(type,name,parent)
    , m_position(3,0)
    , m_spread(0)
    , m_really_initialized(false)
    , m_min(3,0)
    , m_max(3,0)
    , m_detelem(0)
{
    declareInterface<IHepMCEventMutator>(this);

    declareProperty("Volume",m_volume="",
                    "Name (TDS path) of physical volume in which to position events");
    declareProperty("Strategy",m_strategy="FullVolume",
                    "Strategy for positioning events.");
    declareProperty("Mode",m_mode="Uniform","Mode for choosing positions.");
    declareProperty("Position",m_position,"Fixed/mean point");
    declareProperty("Spread",m_spread,"Spread for Gausian/Uniform smearing");
    declareProperty("FillVolumes",m_fillVolumes,
                    "Logical volume names to fill with generated events");
    declareProperty("FillMaterials",m_fillMaterials,
                    "Material names to fill with generated events");
}

GtPositionerTool::~GtPositionerTool()
{
}

const IDetectorElement* GtPositionerTool::detelem()
{
    if (m_detelem) return m_detelem;

    if ("" == m_volume) return 0;

    if (!existDet<IDetectorElement>(m_volume)) return 0;

    m_detelem = getDet<IDetectorElement>(m_volume);
    return m_detelem;
}

StatusCode GtPositionerTool::initialize()
{
    return StatusCode::SUCCESS;
}

StatusCode GtPositionerTool::really_initialize()
{
    m_really_initialized = true;

    const IDetectorElement* detelem = this->detelem();
    if (!detelem) {
        fatal() << "Could not get detector element: \"" << m_volume << "\"" << endreq;
        return StatusCode::FAILURE;
    }

    // Get volume's solid's bounding box
    const ILVolume* lvol = detelem->geometry()->lvolume();
    if (!lvol) {
        fatal() << "Could not get volume: \"" << m_volume << "\"" << endreq;
        return StatusCode::FAILURE;
    }

    const SolidBase* solid = dynamic_cast<const SolidBase*>(lvol->solid());
    m_min[0] = solid->xMin(); m_max[0] = solid->xMax();
    m_min[1] = solid->yMin(); m_max[1] = solid->yMax();
    m_min[2] = solid->zMin(); m_max[2] = solid->zMax();

    // sanity check position
    for (int ind=0; ind<3; ++ind) {
        if (m_min[ind] > m_position[ind] || m_max[ind] < m_position[ind]) {
            // Let this be only warning - maybe user wants to do this.
            warning() << "GtPositioner tool given position not in detector element "
                      << m_volume << "'s logical volume's " << *lvol
                      << " bounding box, index " << ind << ": " 
                      << m_min[ind] << " ? " << m_position[ind] << " ? " << m_max[ind]
                      << endreq;
            //return StatusCode::FAILURE;
        }
    }


    // Possibly shrink trial range if requested spread is less than
    // request uniform spread.
    if (m_mode == "Uniform") {
        for (int ind=0; ind<3; ++ind) {
            if (m_min[ind] < m_position[ind] - m_spread)
                m_min[ind] = m_position[ind] - m_spread;
            if (m_max[ind] > m_position[ind] + m_spread)
                m_max[ind] = m_position[ind] + m_spread;
        }
    }

    // Set up random numbers
    IRndmGenSvc *rgs = 0;
    if (service("RndmGenSvc",rgs,true).isFailure()) {
        fatal() << "Failed to get random service" << endreq;
        return StatusCode::FAILURE;        
    }
    if (m_uni.initialize(rgs, Rndm::Flat(0,1)).isFailure()) {
        fatal() << "Failed to initialize uniform random numbers" << endreq;
        return StatusCode::FAILURE;
    }
    if (m_gauss.initialize(rgs, Rndm::Gauss(0,1)).isFailure()) {
        fatal() << "Failed to initialize Gaussian random numbers" << endreq;
        return StatusCode::FAILURE;
    }

    return StatusCode::SUCCESS;
}

StatusCode GtPositionerTool::finalize()
{
    return StatusCode::SUCCESS;
}

StatusCode GtPositionerTool::mutate(HepMC::GenEvent& event)
{
    if (! m_really_initialized) {
        info() << "really_initialized with detector element: \"" 
               << m_volume << "\"" << endreq;
        StatusCode sc = this->really_initialize();
        if (sc.isFailure()) return sc;
    }


    // short circuit when in Fixed mode
    if ("Fixed" == m_mode or "Relative" == m_mode) {
        this->setPosition(m_position,event);
        return StatusCode::SUCCESS;
    }


    std::vector<double> local_pos(3,0);
    if (this->genPosition(local_pos)) {
        this->setPosition(local_pos,event);
        return StatusCode::SUCCESS;
    }
    
    fatal() << "Failed to generate position" << endreq;
    return StatusCode::FAILURE;
}

bool GtPositionerTool::getPoint(std::vector<double>& pos)
{
    // Keep drawing Gaussians until we get, for each dimension, one
    // that is inside the bounds.
    if ("Smeared" == m_mode) {
        for (int ind=0; ind<3; ++ind) {
            do {
                pos[ind] = m_gauss() * m_spread + m_position[ind];
            } while(pos[ind] < m_min[ind] || pos[ind] > m_max[ind]);
        }
        return true;
    }

    // Simply draw uniforms 
    if ("Uniform" == m_mode) {
        for (int ind=0; ind<3; ++ind)
            pos[ind] = m_uni() * (m_max[ind]-m_min[ind]) + m_min[ind];
        return true;
    }

    // Huh, why are we here?
    fatal() << "Unknown positioning mode: " << m_mode << endreq;
    return false;
}

// Generate a local position in the logical volume
bool GtPositionerTool::genPosition(std::vector<double>& local_pos)
{
    if ("Surface" == m_strategy) {
        //generate point on surface, 
        // munge event
        // return
        //...
        fatal() << "Surface strategy not yet supported" << endreq;
        return false;
    }


    if (!this->getPoint(local_pos)) {
        fatal () << "Failed to get local point" << endreq;
        return false;
    }

    Gaudi::XYZPoint local_point(local_pos[0],local_pos[1],local_pos[2]);

    // Check if inside logical volume
    const ILVolume* lvol = m_detelem->geometry()->lvolume();
    if (!lvol) {
        fatal() << "Could not get volume: \"" << m_volume << "\"" << endreq;
        return false;
    }

    // not inside, try again
    if (!lvol->isInside(local_point)) 
        return this->genPosition(local_pos);

    if ("AvoidDaughters" == m_strategy) {
        // Place in top-most volume only
        ILVolume::ReplicaPath replicaPath;
        if (lvol->belongsTo(local_point,1,replicaPath).isFailure()){
          fatal() << "Invalid point in " << m_volume << ": " << local_point
                  << endreq;
        }
        // Check if point is inside a daughter volume
        if (replicaPath.size() == 0)
          return true;
        else
          return this->genPosition(local_pos);
    } else if("VolumeType" == m_strategy || "Material" == m_strategy) {
      // Fill by physical volume or material type
      //info() << "Starting volume search" << endreq;
      ILVolume::PVolumePath pVolumePath;
      //info() << "   Local point: " << local_point << endreq;
      if (lvol->belongsTo(local_point,-1,pVolumePath).isFailure()){
        fatal() << "Invalid point in " << m_volume << ": " << local_point
                << endreq;
      }
      // Check if point is inside a daughter volume
      //info() << "Depth of PVolumes: " <<  pVolumePath.size() << endreq;
      //for(unsigned int idx = 0; idx < pVolumePath.size(); idx++){
      //        info() << "    PVolume: " << pVolumePath[idx]->name() << endreq;
      //}
      if (pVolumePath.size() == 0){
        // No daughter volumes, so check this volume
        if("Material" == m_strategy){
          if(std::find(m_fillMaterials.begin(), m_fillMaterials.end(),
                       lvol->materialName()) != m_fillMaterials.end()){
            debug() << "Placing in " <<  lvol->name() << endreq;
            return true;
          }
        }
        if("VolumeType" == m_strategy){
          if(std::find(m_fillVolumes.begin(), m_fillVolumes.end(),
                       baseName(lvol->name())) == m_fillVolumes.end()){
            //info() << "    Not the volume we want" << endreq;
            return this->genPosition(local_pos);
          }
        }
        return this->genPosition(local_pos);
      }else{
        const IPVolume* lastPV = pVolumePath[pVolumePath.size() - 1];
        //info() << "    Last PVolume" << endreq;
        //info() << "      Name: " << lastPV->name() << endreq;
        //info() << "      Material: " << lastPV->lvolume()->materialName() 
        //       << endreq;
        if("VolumeType" == m_strategy){
          if(std::find(m_fillVolumes.begin(), m_fillVolumes.end(),
                       baseName(lastPV->lvolumeName())) == m_fillVolumes.end()){
            //info() << "    Not the volume we want" << endreq;
            return this->genPosition(local_pos);
          }
        }
        if("Material" == m_strategy){
          std::string currentMaterial = 
            baseName(lastPV->lvolume()->materialName());
          //info() << "    Current Material: " << currentMaterial << endreq;
          if(std::find(m_fillMaterials.begin(), m_fillMaterials.end(),
                       currentMaterial) == m_fillMaterials.end()){
            //info() << "    Not the material we want" << endreq;
            return this->genPosition(local_pos);
          }
        }
        debug() << "Placing in " <<  lastPV->name() << endreq;
        return true;
      }
    } else if ("FullVolume" == m_strategy) {
      return true;
    }
    
    fatal() << "Failed to generate volume using strategy " << m_strategy << endreq;
    return false;
}

std::string GtPositionerTool::baseName(const std::string& path){
  size_t pos = path.rfind('/');
  if( pos == std::string::npos ) return path;
  return path.substr(pos+1);
}

#include "CLHEP/Units/SystemOfUnits.h"

void GtPositionerTool::setPosition(std::vector<double>& vpos, HepMC::GenEvent& event)
{
    HepMC::FourVector position(vpos[0],vpos[1],vpos[2],0);
    debug() << "Setting position: " 
            <<position.x()<<","<<position.y()<<","<<position.z()
            << endreq;
    HepMC::GenEvent::vertex_iterator vtx, done = event.vertices_end();
    for (vtx = event.vertices_begin(); vtx != done; ++vtx) {
        if ("Relative" == m_mode) {
            HepMC::FourVector relative_pos = (*vtx)->position();
            (*vtx)->set_position(HepMC::FourVector(position.x()+relative_pos.x(),
                                                   position.y()+relative_pos.y(),
                                                   position.z()+relative_pos.z(),
                                                   position.t()+relative_pos.t()));
        }
        else{
            HepMC::FourVector position4d = position;
            position4d.setT((*vtx)->position().t());
            (*vtx)->set_position(position4d);
            debug() << "\tvertex @ " << "("
                    << (*vtx)->position().x()/CLHEP::cm << ","
                    << (*vtx)->position().y()/CLHEP::cm << ","
                    << (*vtx)->position().z()/CLHEP::cm << ","
                    << (*vtx)->position().t()/CLHEP::second << ") [cm,cm,cm,second]" << endreq;

        }
    }
}

---