/search.css" rel="stylesheet" type="text/css"/> /search.js">
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

In This Package:

GtGunGenTool.cc
Go to the documentation of this file.
00001 
00015 #include "GtGunGenTool.h"
00016 
00017 #include "GaudiKernel/IRndmGenSvc.h"
00018 #include "GaudiKernel/IParticlePropertySvc.h"
00019 #include "GaudiKernel/ParticleProperty.h"
00020 #include "GaudiKernel/Vector3DTypes.h"
00021 
00022 #include "HepMC/GenEvent.h"
00023 #include "HepMC/GenParticle.h"
00024 #include "HepMC/GenVertex.h"
00025 
00026 #include "CLHEP/Units/SystemOfUnits.h"
00027 #include "CLHEP/Vector/ThreeVector.h"
00028 
00029 using namespace CLHEP;
00030 
00031 GtGunGenTool::GtGunGenTool(const std::string& type,
00032                            const std::string& name, 
00033                            const IInterface* parent)
00034     : GaudiTool(type,name,parent)
00035     , m_direction(3,0)
00036     , m_pid(0)
00037     , m_mass(0)
00038 {
00039     declareInterface<IHepMCEventMutator>(this);
00040 
00041     declareProperty("ParticleName", m_particleName = "", 
00042                     "PDG particle name");
00043     declareProperty("ParticlesPerEvent",m_particlesPerEvent = 1,
00044                     "Number of particles per event to generate.");
00045     declareProperty("MomentumMode", m_momentumMode = "Fixed","Momentum mode");
00046     declareProperty("Momentum", m_momentum = 1.0*MeV,"Particle momentum");
00047     m_momentum.verifier().setLower(0.*eV);
00048     declareProperty("MomentumSpread", m_momentumSpread = 0,"Momentum spread");
00049     declareProperty("MomentumInterpretation",m_momentumInterpretation = "Momentum","Valid interpretation values are Momentum(default), KineticEnergy, TotalEnergy") ;
00050     declareProperty("DirectionMode", m_directionMode = "Fixed","Direction mode");
00051     m_direction[2] = 1.0;
00052     declareProperty("Direction", m_direction,"Direction");
00053     declareProperty("DirectionSpread",  m_directionSpread = 0,"Direction spread for cos(zenith angle)");
00054     declareProperty("PolarizeMode", m_polarizeMode = "None",
00055                     "Type of polarization of the particle.");
00056 }
00057 
00058 GtGunGenTool::~GtGunGenTool()
00059 {
00060 }
00061 
00062 StatusCode GtGunGenTool::initialize()
00063 {
00064     this->GaudiTool::initialize();
00065 
00066     // Get mass of particle for later calculation of energy
00067     IParticlePropertySvc * ppSvc = 
00068         svc< IParticlePropertySvc >( "ParticlePropertySvc" , true ) ;
00069     ParticleProperty* particle = 0;
00070 
00071     // First try to find via PID
00072     if (m_pid) {
00073         particle = ppSvc->findByStdHepID(m_pid);
00074         if (particle) {
00075             m_mass = particle->mass();
00076             m_particleName = particle->particle();
00077         }
00078     }
00079     // If failed, then try via name
00080     if (!particle) {
00081         particle = ppSvc->find(m_particleName);
00082         if (particle) {
00083             m_mass = particle->mass();
00084             m_pid = particle->pdgID();
00085         }
00086     }
00087     // If still fail, complain
00088     if (!particle) {
00089         fatal() << "Failed to find particle named \"" 
00090                 << m_particleName << "\" and with ID " 
00091                 << m_pid << endreq;
00092         return StatusCode::FAILURE;
00093     }
00094 
00095     IRndmGenSvc *rgs = 0;
00096     if (service("RndmGenSvc",rgs,true).isFailure()) {
00097         fatal() << "Failed to get random service" << endreq;
00098         return StatusCode::FAILURE;        
00099     }
00100 
00101     StatusCode sc;
00102 
00103     sc = m_gauss.initialize(rgs, Rndm::Gauss(0,1));
00104     if (sc.isFailure()) {
00105         fatal() << "Failed to initialize gaussian random numbers" << endreq;
00106         return StatusCode::FAILURE;
00107     }
00108     sc = m_uni.initialize(rgs, Rndm::Flat(0,1));
00109     if (sc.isFailure()) {
00110         fatal() << "Failed to initialize uniform random numbers" << endreq;
00111         return StatusCode::FAILURE;
00112     }
00113 
00114     // force normalization
00115     Hep3Vector unit(m_direction[0],m_direction[1],m_direction[2]);
00116     unit = unit.unit();
00117     for (int ind=0; ind<3; ++ind) m_direction[ind] = unit[ind];
00118  
00119     return StatusCode::SUCCESS;
00120 }
00121 
00122 StatusCode GtGunGenTool::finalize()
00123 {
00124 
00125     return this->GaudiTool::finalize();
00126 }
00127 
00128 
00129 StatusCode GtGunGenTool::mutate(HepMC::GenEvent& event)
00130 {
00131     for (int ind=0; ind<m_particlesPerEvent; ++ind) {
00132         if (this->oneVertex(event).isFailure()) 
00133             return StatusCode::FAILURE;
00134     }
00135     return StatusCode::SUCCESS;
00136 }
00137 StatusCode GtGunGenTool::oneVertex(HepMC::GenEvent& event)
00138 {
00139     CLHEP::Hep3Vector dir;
00140 
00141     if (m_directionMode == "Fixed") {
00142         dir.set(m_direction[0],m_direction[1],m_direction[2]);
00143     }
00144     else {
00145         // cos(angle) from mean direction
00146         double costh=0;
00147         do {
00148             costh = generateNumber(m_directionMode,1,m_directionSpread);
00149         } while (costh < -1.0||costh > 1.0);
00150         
00151         double sinth = sqrt(1-costh*costh);
00152 
00153         // azimuth angle around mean direction
00154         double phi = generateNumber("Uniform",0.0,360*CLHEP::degree);
00155         double cosphi = cos(phi);
00156         double sinphi = sin(phi);
00157         
00158         dir.set(sinth*cosphi,sinth*sinphi,costh);
00159 
00160         // rotate along Y-axis by the mean direction's theta
00161         dir.rotateY(acos(m_direction[2]));
00162         // rotate along Z-axis by the mean direction's phi
00163         dir.rotateZ(atan2(m_direction[1],m_direction[0]));
00164 
00165     }
00166 
00167     dir = dir.unit();
00168 
00169     double momentum;
00170     do{
00171       momentum = generateNumber(m_momentumMode,m_momentum,m_momentumSpread);
00172     } while (momentum < 0.0);
00173 
00174     // What is the interpretation of momentum?
00175     double energy = sqrt(momentum*momentum + m_mass*m_mass); 
00176     if (m_momentumInterpretation == "TotalEnergy") {
00177       energy = momentum ; 
00178       if ( energy < m_mass ) {
00179         fatal() << "Cannot set total energy(MeV) = " << energy/CLHEP::MeV 
00180                 << " of " << m_particleName 
00181                 << " to be less than particle mass(MeV) = " << m_mass/CLHEP::MeV << endreq;
00182         return StatusCode::FAILURE;
00183       }
00184       momentum = sqrt( (energy+m_mass) * (energy-m_mass) );
00185     }
00186 
00187     if (m_momentumInterpretation == "KineticEnergy") {
00188       energy = momentum + m_mass;
00189       momentum = sqrt( (energy+m_mass) * (energy-m_mass) );
00190     }
00191     
00192     HepMC::ThreeVector p(momentum*dir);
00193     HepMC::FourVector fourmom(p.x(),p.y(),p.z(),energy);
00194 
00195     debug() << m_particleName << " (" << m_pid << ") generated with 4momentum = " 
00196             << momentum
00197             << " dir = " << dir
00198             << " momentum(MeV) = " << momentum/CLHEP::MeV
00199             << " energy(MeV) = " << energy/CLHEP::MeV
00200             << " kinetic energy(MeV) = " << (energy-m_mass)/CLHEP::MeV
00201             << endreq;
00202 
00203     HepMC::GenParticle* particle = new HepMC::GenParticle(fourmom,m_pid,1/*=status*/);
00204 
00205     if (m_polarizeMode == "Random") {
00206         double phi = generateNumber("Uniform",0.0,360*CLHEP::degree);
00207         particle->set_polarization(HepMC::Polarization(0.5*HepMC::HepMC_pi,phi));
00208     }
00209 
00210     HepMC::GenVertex* vertex = event.signal_process_vertex();
00211     if (!vertex) {
00212         vertex = new HepMC::GenVertex(HepMC::FourVector(0,0,0,0));
00213         event.set_signal_process_vertex(vertex);
00214     }
00215 
00216     vertex->add_particle_out(particle);
00217 
00218     return StatusCode::SUCCESS;
00219 }
00220 
00221 double GtGunGenTool::generateNumber(const std::string& mode, double mean, double spread)
00222 {
00223     if (mode == "Smeared") {
00224         return m_gauss() * spread + mean;
00225     }
00226     if (mode == "Uniform") {
00227         return (m_uni()*2-1) * spread +mean ;
00228     }
00229     return mean;
00230 }
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

Generated on Fri May 16 2014 10:18:50 for GenTools by doxygen 1.7.4