30 #include "CLHEP/Vector/LorentzVector.h"
32 #include "reportingUtils.h"
33 #include "starlightconstants.h"
34 #include "starlightparticlecodes.h"
37 static const std::string starlight_stream =
"STARLIGHT";
56 ATH_MSG_INFO(
"===> January 20 2011 STARLIGHT INTERFACE VERSION. \n" );
67 return StatusCode::FAILURE;
83 ATH_MSG_INFO(
"===> dumping starlight events to lhef format. \n" );
87 return StatusCode::SUCCESS;
97 const EventContext& ctx = Gaudi::Hive::currentContext();
113 int numberofTracks =
m_event->getParticles()->size();
114 int numberOfVertices = 1;
117 <<
" with " << numberOfVertices <<
" vertices "
118 <<
" and " << numberofTracks <<
" tracks" );
120 <<
" with " << numberofTracks <<
" tracks" );
123 std::vector<starlightParticle>::const_iterator
part =
128 << starlightParticleCodes::jetsetToGeant((*part).getCharge() * (*part).getPdgCode()) <<
" "
129 << (*part).GetPx() <<
" " << (*part).GetPy() <<
" "<< (*part).GetPz()
130 <<
" " <<
m_events <<
" " << ipart <<
" " << 0 <<
" "
131 << (*part).getCharge() * (*part).getPdgCode() );
136 return StatusCode::SUCCESS;
144 return StatusCode::SUCCESS;
158 evt->add_vertex( v1 );
163 std::vector<starlightParticle>::const_iterator
part =
171 int pid = (*part).getPdgCode();
172 int charge = (*part).getCharge();
174 int pidsign =
pid/std::abs(
pid);
177 if( chsign != pidsign && chsign != 0)
pid = -
pid;
179 double px = (*part).GetPx();
180 double py = (*part).GetPy();
181 double pz = (*part).GetPz();
182 double e = (*part).GetE();
184 if(std::abs(
pid)==13) {
189 if(std::abs(
pid)==22) {
224 double e = sqrt(px_tot*px_tot + py_tot*py_tot + pz_tot*pz_tot +
mass*
mass);
225 v1->add_particle_out(
234 return StatusCode::SUCCESS;
241 std::string lheFilename =
"events.lhe";
242 std::ofstream lheStream;
243 lheStream.open(lheFilename.c_str(), std::ofstream::trunc);
249 lheStream <<
"<LesHouchesEvents version=\"1.0\">\n";
250 lheStream <<
"<!--\n";
251 lheStream <<
"File generated using Starlight \n";
252 lheStream <<
"-->\n";
260 lheStream <<
"<init>\n";
262 lheStream <<
" 1.000000e+00 0.000000e+00 1.000000e+00 9999\n";
263 lheStream <<
"</init>\n";
266 std::unique_ptr<upcEvent> uevent(
new upcEvent);
269 lheStream <<
"<event>\n";
276 CLHEP::HepLorentzVector photon_system(0);
278 std::vector<starlightParticle>::const_iterator
part = (uevent->getParticles())->
begin();
279 for (
part = uevent->getParticles()->begin();
part != uevent->getParticles()->end(); ++
part, ++ipart)
281 CLHEP::HepLorentzVector particle_sl((*part).GetPx(), (*part).GetPy(), (*part).GetPz(), (*part).GetE());
282 photon_system += particle_sl;
283 ptscale += std::sqrt((*part).GetPx()*(*part).GetPx() + (*part).GetPy()*(*part).GetPy());
287 ptscale /=
static_cast<float> (ipart);
288 lheStream <<
" 4 9999 1.000000e+00 "<<ptscale<<
" 7.297e-03 2.569093e-01\n";
291 lheStream <<
" -11 -1 0 0 0 0 0.0000000000e+00 0.0000000000e+00 "
292 << photon_system.m()/2.*
std::exp(photon_system.rapidity())<<
" "
293 <<photon_system.m()/2.*
std::exp(photon_system.rapidity())
294 <<
" 0.0000000000e+00 0. 9.\n";
295 lheStream <<
" 11 -1 0 0 0 0 0.0000000000e+00 0.0000000000e+00 "
296 << -photon_system.m()/2.*
std::exp(-photon_system.rapidity())<<
" "
297 <<photon_system.m()/2.*
std::exp(-photon_system.rapidity())
298 <<
" 0.0000000000e+00 0. 9.\n";
302 lheStream <<
" 22 -1 0 0 0 0 0.0000000000e+00 0.0000000000e+00 "
303 << photon_system.m()/2.*
std::exp(photon_system.rapidity())<<
" "
304 <<photon_system.m()/2.*
std::exp(photon_system.rapidity())
305 <<
" 0.0000000000e+00 0. 9.\n";
306 lheStream <<
" 22 -1 0 0 0 0 0.0000000000e+00 0.0000000000e+00 "
307 << -photon_system.m()/2.*
std::exp(-photon_system.rapidity())<<
" "
308 <<photon_system.m()/2.*
std::exp(-photon_system.rapidity())
309 <<
" 0.0000000000e+00 0. 9.\n";
312 for (
part = uevent->getParticles()->begin();
part != uevent->getParticles()->end(); ++
part, ++ipart)
314 int pid = (*part).getPdgCode();
315 int charge = (*part).getCharge();
317 int pidsign =
pid/std::abs(
pid);
319 if( chsign != pidsign )
pid = -
pid;
321 double px = (*part).GetPx();
322 double py = (*part).GetPy();
323 double pz = (*part).GetPz();
324 double e = (*part).GetE();
325 double mass = (*part).getMass();
330 lheStream <<
pid<<
" 1 1 2 0 0 "<<
px<<
" "<<
py<<
" "<<
pz<<
" "<<
e<<
" "<<
mass<<
" 0. 9.\n";
333 lheStream <<
"</event>\n";
337 lheStream <<
"</LesHouchesEvents>";
353 "problems initializing input parameters. cannot initialize starlight.";
360 ATH_MSG_WARNING(
"problems initializing input parameters. cannot initialize starlight. " );
376 std::string myparam = mystring.
piece<std::string>(1);
377 if (myparam ==
"beam1Z")
381 else if (myparam ==
"beam1A")
385 else if (myparam ==
"beam2Z")
389 else if (myparam ==
"beam2A")
393 else if (myparam ==
"beam1Gamma")
397 else if (myparam ==
"beam2Gamma")
401 else if (myparam ==
"maxW")
405 else if (myparam ==
"minW")
409 else if (myparam ==
"nmbWBins")
413 else if (myparam ==
"maxRapidity")
417 else if (myparam ==
"nmbRapidityBins")
421 else if (myparam ==
"accCutPt")
425 else if (myparam ==
"minPt")
429 else if (myparam ==
"maxPt")
433 else if (myparam ==
"accCutEta")
437 else if (myparam ==
"minEta")
441 else if (myparam ==
"maxEta")
445 else if (myparam ==
"productionMode")
449 else if (myparam ==
"axionMass")
453 else if (myparam ==
"nmbEventsTot")
457 else if (myparam ==
"prodParticleId")
461 else if (myparam ==
"randomSeed")
465 else if (myparam ==
"outputFormat")
469 else if (myparam ==
"beamBreakupMode")
473 else if (myparam ==
"interferenceEnabled")
477 else if (myparam ==
"interferenceStrength")
481 else if (myparam ==
"coherentProduction")
485 else if (myparam ==
"incoherentFactor")
489 else if (myparam ==
"maxPtInterference")
493 else if (myparam ==
"nmbPtBinsInterference")
497 else if (myparam ==
"xsecMethod")
501 else if (myparam ==
"nThreads")
505 else if (myparam ==
"pythFullRec")
511 ATH_MSG_ERROR(
" ERROR in STARLIGHT INITIALIZATION PARAMETERS "
512 << myparam <<
" is an invalid parameter !" );