00001
00002
00003
00004
00005
00006
00007
00008 #include "SingleParticleGun/SingleParticleGun.h"
00009 #include "GeneratorModule/GeneratorName.h"
00010
00011
00012 #include "GaudiKernel/MsgStream.h"
00013
00014
00015 #include <math.h>
00016 #include "CLHEP/Units/PhysicalConstants.h"
00017
00018 #include "CLHEP/Random/RandFlat.h"
00019 #include "CLHEP/Random/RandGauss.h"
00020 #include "BesKernel/IBesRndmGenSvc.h"
00021
00022 #include "HepPDT/ParticleDataTable.hh"
00023 #include "HepPDT/ParticleData.hh"
00024 #ifndef ENABLE_BACKWARDS_COMPATIBILITY
00025 typedef double HepDouble;
00026 #endif
00027 using HepMC::GenVertex;
00028 using HepMC::GenParticle;
00029
00030
00031
00032
00033 SingleParticleGun::SingleParticleGun(const std::string& name, ISvcLocator* pSvcLocator)
00034 : GenModule(name,pSvcLocator)
00035 {
00036
00037 declareProperty("Mode",m_Emode = SingleEnergyMode::PtMode);
00038
00039
00040
00041
00042 declareProperty("Pt",m_requestedPt = 5.0);
00043 declareProperty("E",m_requestedE = 5.0);
00044 declareProperty("Phi",m_requestedPhi = 0.0);
00045 declareProperty("Theta",m_requestedTheta = 3.14/4.0);
00046
00047 declareProperty("VertexX",m_requestedX = 0.0);
00048 declareProperty("VertexY",m_requestedY = 0.0);
00049 declareProperty("VertexZ",m_requestedZ = 0.0);
00050 declareProperty("Time",m_requestedT = 0.0);
00051
00052 declareProperty("MinPt",m_minPt = 1.);
00053 declareProperty("MinE",m_minE = 1.);
00054 declareProperty("MinTheta",m_minTheta = 0.);
00055 declareProperty("MinPhi",m_minPhi = 0.);
00056
00057 declareProperty("MaxPt",m_maxPt = 100.);
00058 declareProperty("MaxE",m_maxE = 100.);
00059 declareProperty("MaxTheta",m_maxTheta=CLHEP::pi);
00060 declareProperty("MaxPhi",m_maxPhi = CLHEP::twopi);
00061
00062 declareProperty("SigmaPt",m_sigmaPt = 0.1);
00063 declareProperty("SigmaE",m_sigmaE = 0.1);
00064
00065 declareProperty("SigmaTheta",m_sigmaTheta= 0.1);
00066 declareProperty("SigmaPhi",m_sigmaPhi = 0.1);
00067
00068 declareProperty("ModePt",m_PtGenMode = SingleParticleGunGenMode::FixedMode);
00069 declareProperty("ModeE",m_EGenMode = SingleParticleGunGenMode::FixedMode);
00070 declareProperty("ModeTheta",m_ThetaGenMode=SingleParticleGunGenMode::FixedMode);
00071
00072 declareProperty("ModePhi",m_PhiGenMode=SingleParticleGunGenMode::FixedMode);
00073
00074
00075 declareProperty("PdgCode",m_pdgCode=211);
00076 }
00077
00078
00079 SingleParticleGun::~SingleParticleGun(){
00080
00081
00082
00083 }
00084
00085
00086 StatusCode SingleParticleGun::genInitialize() {
00087
00088
00089
00090
00091
00092 MsgStream log(messageService(), name());
00093
00094 static const bool CREATEIFNOTTHERE(true);
00095 StatusCode RndmStatus = service("BesRndmGenSvc", p_BesRndmGenSvc, CREATEIFNOTTHERE);
00096 if (!RndmStatus.isSuccess() || 0 == p_BesRndmGenSvc)
00097 {
00098 log << MSG::ERROR << " Could not initialize Random Number Service" << endreq;
00099 return RndmStatus;
00100 }
00101
00102
00103
00104 const HepPDT::ParticleData* particle = m_particleTable->particle(HepPDT::ParticleID(abs(m_pdgCode)));
00105 m_mass = particle->mass().value();
00106
00107
00108
00109
00110 switch (m_Emode){
00111 case SingleEnergyMode::PtMode:
00112 if(!m_PtGenMode && (m_minPt > m_requestedPt || m_maxPt < m_requestedPt)
00113 || m_maxPt < m_minPt ) {
00114 log << MSG:: ERROR << " Pt min and max out of range. \n" <<
00115 " Will set Pt mode to Fixed!!!" << endreq;
00116 m_PtGenMode = SingleParticleGunGenMode::FixedMode;
00117 }
00118 if(!m_ThetaGenMode && (m_minTheta > m_requestedTheta || m_maxTheta < m_requestedTheta)
00119 || m_maxTheta < m_minTheta ) {
00120 log << MSG:: ERROR << " Theta min and max out of range. \n" <<
00121 " Will set Theta mode to Fixed!!!" << endreq;
00122 m_ThetaGenMode = SingleParticleGunGenMode::FixedMode;
00123 }
00124 if(!m_PhiGenMode && (m_minPhi > m_requestedPhi || m_maxPhi < m_requestedPhi)
00125 || m_maxPhi < m_minPhi ) {
00126 log << MSG:: ERROR << " Phi min and max out of range. \n" <<
00127 " Will set Phi mode to Fixed!!!" << endreq;
00128 m_PhiGenMode = SingleParticleGunGenMode::FixedMode;
00129 }
00130 break;
00131 case SingleEnergyMode::EMode:
00132 if(!m_EGenMode && (m_minE > m_requestedE || m_maxE < m_requestedE)
00133 || m_maxE < m_minE ) {
00134 log << MSG:: ERROR << " E min and max out of range. \n" <<
00135 " Will set E mode to Fixed!!!" << endreq;
00136 m_EGenMode = SingleParticleGunGenMode::FixedMode;
00137 }
00138 if(!m_ThetaGenMode && (m_minTheta > m_requestedTheta || m_maxTheta < m_requestedTheta)
00139 || m_maxTheta < m_minTheta ) {
00140 log << MSG:: ERROR << " Theta min and max out of range. \n" <<
00141 " Will set Theta mode to Fixed!!!" << endreq;
00142 m_ThetaGenMode = SingleParticleGunGenMode::FixedMode;
00143 }
00144 if(!m_PhiGenMode && (m_minPhi > m_requestedPhi || m_maxPhi < m_requestedPhi)
00145 || m_maxPhi < m_minPhi ) {
00146 log << MSG:: ERROR << " Phi min and max out of range. \n" <<
00147 " Will set Phi mode to Fixed!!!" << endreq;
00148 m_PhiGenMode = SingleParticleGunGenMode::FixedMode;
00149 }break;
00150 }
00151 m_events = 0;
00152 return StatusCode::SUCCESS;
00153 }
00154
00155
00156 StatusCode SingleParticleGun::callGenerator() {
00157
00158
00159 ++m_events;
00160
00161 CLHEP::HepRandomEngine* engine = p_BesRndmGenSvc->GetEngine("SINGLE");
00162 const long* s = engine->getSeeds();
00163 m_seeds.clear();
00164 m_seeds.push_back(s[0]);
00165 m_seeds.push_back(s[1]);
00166
00167
00168
00169 double pt,phi,theta,E,px,py,pz,p;
00170 switch(m_Emode){
00171
00172 case SingleEnergyMode::PtMode:
00173 pt = generateValue(m_PtGenMode,m_requestedPt, m_sigmaPt,
00174 m_minPt, m_maxPt);
00175 theta = generateValue(m_ThetaGenMode,m_requestedTheta, m_sigmaTheta,
00176 m_minTheta, m_maxTheta);
00177 phi = generateValue(m_PhiGenMode,m_requestedPhi, m_sigmaPhi,
00178 m_minPhi, m_maxPhi);
00179
00180
00181
00182 px = pt*cos(phi);
00183 py = pt*sin(phi);
00184 pz = pt/tan(theta);
00185 m_fourMom.setVectM(CLHEP::Hep3Vector(px,py,pz),m_mass);
00186 m_fourPos.set(m_requestedX,m_requestedY,m_requestedZ,m_requestedT);
00187 return StatusCode::SUCCESS;
00188
00189 case SingleEnergyMode::EMode:
00190 E = generateValue(m_EGenMode,m_requestedE, m_sigmaE,
00191 m_minE, m_maxE);
00192 theta = generateValue(m_ThetaGenMode,m_requestedTheta, m_sigmaTheta,
00193 m_minTheta, m_maxTheta);
00194 phi = generateValue(m_PhiGenMode,m_requestedPhi, m_sigmaPhi,
00195 m_minPhi, m_maxPhi);
00196
00197
00198
00199 if(E*E-m_mass*m_mass < 0.){
00200 MsgStream log(messageService(), name());
00201 log << MSG::ERROR << "You have Generated a Tachyon!! Increase energy or change particle ID" << endreq;
00202 return StatusCode::FAILURE;
00203 }
00204 p=sqrt(E*E-m_mass*m_mass);
00205 px = p*sin(theta)*cos(phi);
00206 py = p*sin(theta)*sin(phi);
00207 pz = p*cos(theta);
00208
00209 m_fourMom.setVectM(CLHEP::Hep3Vector(px,py,pz),m_mass);
00210 m_fourPos.set(m_requestedX,m_requestedY,m_requestedZ,m_requestedT);
00211 return StatusCode::SUCCESS;
00212 }
00213 return StatusCode::SUCCESS;
00214 }
00215
00216
00217 StatusCode SingleParticleGun::genFinalize() {
00218
00219 return StatusCode::SUCCESS;
00220 }
00221
00222
00223 StatusCode SingleParticleGun::fillEvt(GenEvent* evt) {
00224
00225
00226
00227 evt->set_event_number(m_events);
00228 GenVertex* v1 = new GenVertex(m_fourPos);
00229
00230 std::cout << " SingleParticleGun::fillEvt --> " <<std::endl;
00231 std::cout << " Event number:: " << m_events << std::endl;
00232 std::cout << " vertex.x = " << m_fourPos.x()
00233 << " vertex.y = " << m_fourPos.y()
00234 << " vertex.z = " << m_fourPos.z()
00235 << " vertex.t = " << m_fourPos.t() << std::endl;
00236
00237 std::cout << " momentum.x = " << m_fourMom.x()
00238 << " momentum.y = " << m_fourMom.y()
00239 << " momentum.z = " << m_fourMom.z()
00240 << " momentum.t = " << m_fourMom.t() << std::endl;
00241
00242 evt->add_vertex( v1 );
00243
00244
00245 v1->add_particle_out( new GenParticle( m_fourMom, m_pdgCode,1) );
00246
00247 std::cout << " particles_size = " << evt->particles_size()
00248 << " vertices_size = " << evt->vertices_size()
00249 << std::endl;
00250
00251 evt->set_signal_process_id(SINGLE);
00252 evt->set_random_states(m_seeds);
00253
00254 return StatusCode::SUCCESS;
00255 }
00256
00257
00258 double SingleParticleGun::generateValue(int mode, double val, double sigma,
00259 double min, double max) {
00260
00261 double tmp;
00262 int i = 0;
00263 CLHEP::HepRandomEngine* engine = p_BesRndmGenSvc->GetEngine("SINGLE");
00264 const int maxtries = 100;
00265 switch (mode) {
00266 case SingleParticleGunGenMode::FixedMode:
00267 return val;
00268 case SingleParticleGunGenMode::GaussMode:
00269 tmp = max+1.0;
00270 i = 0;
00271 do{
00272
00273 tmp = CLHEP::RandGauss::shoot(engine, (HepDouble) val,(HepDouble) sigma);
00274 i++;
00275 } while ( (tmp<min) || (tmp > max) && (i < maxtries));
00276 if(i>maxtries) {
00277 MsgStream log(messageService(), name());
00278 log << MSG::ERROR << "Cant generate value in range (min, max) "
00279 << val << "\t" << min << "\t" << max << endreq;
00280 }
00281 return tmp;
00282 case SingleParticleGunGenMode::FlatMode:
00283
00284 tmp = CLHEP::RandFlat::shoot(engine, (HepDouble) min,(HepDouble) max);
00285 return tmp;
00286 default:
00287 MsgStream log(messageService(), name());
00288 log << MSG::ERROR << "Unknown Generation Mode" << endreq;
00289 return 0.;
00290 }
00291 }