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();
109 int numberofTracks =
m_event->getParticles()->size();
110 int numberOfVertices = 1;
113 <<
" with " << numberOfVertices <<
" vertices "
114 <<
" and " << numberofTracks <<
" tracks" );
116 <<
" with " << numberofTracks <<
" tracks" );
119 std::vector<starlightParticle>::const_iterator
part =
124 << starlightParticleCodes::jetsetToGeant((*part).getCharge() * (*part).getPdgCode()) <<
" "
125 << (*part).GetPx() <<
" " << (*part).GetPy() <<
" "<< (*part).GetPz()
126 <<
" " <<
m_events <<
" " << ipart <<
" " << 0 <<
" "
127 << (*part).getCharge() * (*part).getPdgCode() );
132 return StatusCode::SUCCESS;
140 return StatusCode::SUCCESS;
154 evt->add_vertex( v1 );
159 std::vector<starlightParticle>::const_iterator
part =
167 int pid = (*part).getPdgCode();
168 int charge = (*part).getCharge();
170 int pidsign =
pid/std::abs(
pid);
173 if( chsign != pidsign && chsign != 0)
pid = -
pid;
175 double px = (*part).GetPx();
176 double py = (*part).GetPy();
177 double pz = (*part).GetPz();
178 double e = (*part).GetE();
180 if(std::abs(
pid)==13) {
185 if(std::abs(
pid)==22) {
220 double e = sqrt(px_tot*px_tot + py_tot*py_tot + pz_tot*pz_tot +
mass*
mass);
221 v1->add_particle_out(
230 return StatusCode::SUCCESS;
237 std::string lheFilename =
"events.lhe";
238 std::ofstream lheStream;
239 lheStream.open(lheFilename.c_str(), std::ofstream::trunc);
245 lheStream <<
"<LesHouchesEvents version=\"1.0\">\n";
246 lheStream <<
"<!--\n";
247 lheStream <<
"File generated using Starlight \n";
248 lheStream <<
"-->\n";
256 lheStream <<
"<init>\n";
258 lheStream <<
" 1.000000e+00 0.000000e+00 1.000000e+00 9999\n";
259 lheStream <<
"</init>\n";
262 std::unique_ptr<upcEvent> uevent(
new upcEvent);
265 lheStream <<
"<event>\n";
268 CLHEP::HepLorentzVector photon_system(0);
270 std::vector<starlightParticle>::const_iterator
part = (uevent->getParticles())->
begin();
271 for (
part = uevent->getParticles()->begin();
part != uevent->getParticles()->end(); ++
part, ++ipart)
273 CLHEP::HepLorentzVector particle_sl((*part).GetPx(), (*part).GetPy(), (*part).GetPz(), (*part).GetE());
274 photon_system += particle_sl;
275 ptscale += std::sqrt((*part).GetPx()*(*part).GetPx() + (*part).GetPy()*(*part).GetPy());
279 ptscale /=
static_cast<float> (ipart);
280 lheStream <<
" 4 9999 1.000000e+00 "<<ptscale<<
" 7.297e-03 2.569093e-01\n";
283 lheStream <<
" -11 -1 0 0 0 0 0.0000000000e+00 0.0000000000e+00 "
284 << photon_system.m()/2.*
std::exp(photon_system.rapidity())<<
" "
285 <<photon_system.m()/2.*
std::exp(photon_system.rapidity())
286 <<
" 0.0000000000e+00 0. 9.\n";
287 lheStream <<
" 11 -1 0 0 0 0 0.0000000000e+00 0.0000000000e+00 "
288 << -photon_system.m()/2.*
std::exp(-photon_system.rapidity())<<
" "
289 <<photon_system.m()/2.*
std::exp(-photon_system.rapidity())
290 <<
" 0.0000000000e+00 0. 9.\n";
294 lheStream <<
" 22 -1 0 0 0 0 0.0000000000e+00 0.0000000000e+00 "
295 << photon_system.m()/2.*
std::exp(photon_system.rapidity())<<
" "
296 <<photon_system.m()/2.*
std::exp(photon_system.rapidity())
297 <<
" 0.0000000000e+00 0. 9.\n";
298 lheStream <<
" 22 -1 0 0 0 0 0.0000000000e+00 0.0000000000e+00 "
299 << -photon_system.m()/2.*
std::exp(-photon_system.rapidity())<<
" "
300 <<photon_system.m()/2.*
std::exp(-photon_system.rapidity())
301 <<
" 0.0000000000e+00 0. 9.\n";
304 for (
part = uevent->getParticles()->begin();
part != uevent->getParticles()->end(); ++
part, ++ipart)
306 int pid = (*part).getPdgCode();
307 int charge = (*part).getCharge();
309 int pidsign =
pid/std::abs(
pid);
311 if( chsign != pidsign )
pid = -
pid;
313 double px = (*part).GetPx();
314 double py = (*part).GetPy();
315 double pz = (*part).GetPz();
316 double e = (*part).GetE();
317 double mass = (*part).getMass();
322 lheStream <<
pid<<
" 1 1 2 0 0 "<<
px<<
" "<<
py<<
" "<<
pz<<
" "<<
e<<
" "<<
mass<<
" 0. 9.\n";
325 lheStream <<
"</event>\n";
329 lheStream <<
"</LesHouchesEvents>";
345 "problems initializing input parameters. cannot initialize starlight.";
352 ATH_MSG_WARNING(
"problems initializing input parameters. cannot initialize starlight. " );
368 std::string myparam = mystring.
piece<std::string>(1);
369 if (myparam ==
"beam1Z")
373 else if (myparam ==
"beam1A")
377 else if (myparam ==
"beam2Z")
381 else if (myparam ==
"beam2A")
385 else if (myparam ==
"beam1Gamma")
389 else if (myparam ==
"beam2Gamma")
393 else if (myparam ==
"maxW")
397 else if (myparam ==
"minW")
401 else if (myparam ==
"nmbWBins")
405 else if (myparam ==
"maxRapidity")
409 else if (myparam ==
"nmbRapidityBins")
413 else if (myparam ==
"accCutPt")
417 else if (myparam ==
"minPt")
421 else if (myparam ==
"maxPt")
425 else if (myparam ==
"accCutEta")
429 else if (myparam ==
"minEta")
433 else if (myparam ==
"maxEta")
437 else if (myparam ==
"productionMode")
441 else if (myparam ==
"axionMass")
445 else if (myparam ==
"nmbEventsTot")
449 else if (myparam ==
"prodParticleId")
453 else if (myparam ==
"randomSeed")
457 else if (myparam ==
"outputFormat")
461 else if (myparam ==
"beamBreakupMode")
465 else if (myparam ==
"interferenceEnabled")
469 else if (myparam ==
"interferenceStrength")
473 else if (myparam ==
"coherentProduction")
477 else if (myparam ==
"incoherentFactor")
481 else if (myparam ==
"maxPtInterference")
485 else if (myparam ==
"nmbPtBinsInterference")
489 else if (myparam ==
"xsecMethod")
493 else if (myparam ==
"nThreads")
497 else if (myparam ==
"pythFullRec")
503 ATH_MSG_ERROR(
" ERROR in STARLIGHT INITIALIZATION PARAMETERS "
504 << myparam <<
" is an invalid parameter !" );