12 #include "GaudiKernel/ServiceHandle.h"
15 #include "GaudiKernel/ITHistSvc.h"
16 #include "Gaudi/Property.h"
90 return StatusCode::SUCCESS;
103 return StatusCode::SUCCESS;
126 return StatusCode::FAILURE;
138 return StatusCode::SUCCESS;
142 for (
const HepMC::GenEvent* itr : *mcColl){
143 HepMC::GenEvent
evt = (*itr);
171 return StatusCode::SUCCESS;
180 const HepMC::FourVector fv(
p->momentum().px() / 1000.,
181 p->momentum().py() / 1000.,
182 p->momentum().pz() / 1000.,
183 p->momentum().e() / 1000.);
194 const HepMC::FourVector fv(
p->momentum().px() * 1000.,
195 p->momentum().py() * 1000.,
196 p->momentum().pz() * 1000.,
197 p->momentum().e() * 1000.);
228 m_FileBeam1 << evt_number <<
"\t" << std::setprecision(17)
229 <<
m_PosRP1.
x() << std::setprecision(17) <<
"\t"
230 <<
m_PosRP1.
y() <<
"\t" << std::setprecision(17)
231 <<
m_PosRP1.
z() <<
"\t" << std::setprecision(17)
232 <<
m_MomRP1.
x() << std::setprecision(17) <<
"\t"
233 <<
m_MomRP1.
y() <<
"\t" << std::setprecision(17)
234 <<
m_MomRP1.
z() <<
"\t" << std::setprecision(17)
236 m_FileBeam2 << evt_number <<
"\t" << std::setprecision(17)
237 <<
m_PosRP3.
x() << std::setprecision(17) <<
"\t"
238 <<
m_PosRP3.
y() <<
"\t" << std::setprecision(17)
239 <<
m_PosRP3.
z() <<
"\t" << std::setprecision(17)
240 <<
m_MomRP3.
x() << std::setprecision(17) <<
"\t"
241 <<
m_MomRP3.
y() <<
"\t" << std::setprecision(17)
242 <<
m_MomRP3.
z() <<
"\t" << std::setprecision(17)
253 std::vector<FPTracker::Point> PosAtRP1;
254 std::vector<FPTracker::Point> PosAtRP3;
255 std::vector<FPTracker::Point> MomAtPR1;
256 std::vector<FPTracker::Point> MomAtPR3;
257 std::vector<double> EnergyRP1;
258 std::vector<double> EnergyRP3;
273 theta = std::acos(std::abs(
p->momentum().pz()) /
mom);
278 int pid =
p->pdg_id();
282 HepMC::FourVector Position =
p->production_vertex()->position();
283 HepMC::FourVector Momentum =
p->momentum();
291 if (
p->momentum().pz() > 0. &&
pid == 2212) {
296 p1->production_vertex()->position().x() / 1000.,
297 p1->production_vertex()->position().y() / 1000.,
298 p1->production_vertex()->position().z() / 1000.,
301 p1->momentum().pz());
317 m_FileBeam1 << evt_number <<
"\t" << std::setprecision(17)
318 <<
m_PosRP1.
x() << std::setprecision(17) <<
"\t"
319 <<
m_PosRP1.
y() <<
"\t" << std::setprecision(17)
320 <<
m_PosRP1.
z() <<
"\t" << std::setprecision(17)
321 <<
m_MomRP1.
x() << std::setprecision(17) <<
"\t"
322 <<
m_MomRP1.
y() <<
"\t" << std::setprecision(17)
323 <<
m_MomRP1.
z() <<
"\t" << std::setprecision(17)
327 }
else if (
p->momentum().pz() < 0. &&
pid == 2212) {
331 p2->production_vertex()->position().x() / 1000.,
332 p2->production_vertex()->position().y() / 1000.,
333 p2->production_vertex()->position().z() / 1000.,
336 p2->momentum().pz());
350 m_FileBeam2 << evt_number <<
"\t" << std::setprecision(17)
351 <<
m_PosRP3.
x() << std::setprecision(17) <<
"\t"
352 <<
m_PosRP3.
y() <<
"\t" << std::setprecision(17)
353 <<
m_PosRP3.
z() <<
"\t" << std::setprecision(17)
354 <<
m_MomRP3.
x() << std::setprecision(17) <<
"\t"
355 <<
m_MomRP3.
y() <<
"\t" << std::setprecision(17)
356 <<
m_MomRP3.
z() <<
"\t" << std::setprecision(17)
367 ATH_MSG_ERROR(
"Strange: More than two protons in this event!");
376 ATH_MSG_INFO(
"Add transproted particle into HepMC event record Beam 1");
379 for(
int i=0;
i<(
int)PosAtRP1.size();
i++){
381 HepMC::FourVector PositionVectorRP1 = HepMC::FourVector(PosAtRP1.at(
i).x()*1000.,PosAtRP1.at(
i).y()*1000.,PosAtRP1.at(
i).z()*1000.,0.*1000.);
382 HepMC::FourVector MomentumVectorRP1 = HepMC::FourVector(MomAtPR1.at(
i).x(),MomAtPR1.at(
i).y(),MomAtPR1.at(
i).z(),EnergyRP1.at(
i));
387 VertexRP1->add_particle_out(std::move(ParticleRP1));
390 evt.add_vertex(std::move(VertexRP1));
393 ATH_MSG_INFO (
"Add transproted particle into HepMC event record Beam 2" );
394 for(
int i=0;
i<(
int)PosAtRP3.size();
i++){
395 HepMC::FourVector PositionVectorRP3 = HepMC::FourVector(PosAtRP3.at(
i).x()*1000.,PosAtRP3.at(
i).y()*1000.,PosAtRP3.at(
i).z()*1000.,0.*1000.);
396 HepMC::FourVector MomentumVectorRP3 = HepMC::FourVector(MomAtPR3.at(
i).x(),MomAtPR3.at(
i).y(),MomAtPR3.at(
i).z(),EnergyRP3.at(
i));
399 VertexRP3->add_particle_out(std::move(ParticleRP3));
401 evt.add_vertex(std::move(VertexRP3));