C++ GenEvent::vertices_end方法代码示例

/// Smearing function
StatusCode GaussSmearVertex::smearVertex(HepMC::GenEvent& theEvent) {

  double dx = m_gaussDist() * sqrt(m_xsig) + m_xmean;
  double dy = m_gaussDist() * sqrt(m_ysig) + m_ymean;
  double dz = m_gaussDist() * sqrt(m_zsig) + m_zmean;
  double dt = m_gaussDist() * sqrt(m_tsig) + m_tmean;

  Gaudi::LorentzVector dpos(dx, dy, dz, dt);

  debug() << "Smearing vertices by " << dpos << endmsg;

  for (auto vit = theEvent.vertices_begin(); vit != theEvent.vertices_end(); ++vit) {
    Gaudi::LorentzVector pos((*vit)->position());
    pos += dpos;
    (*vit)->set_position(HepMC::FourVector(pos.x(), pos.y(), pos.z(), pos.t()));
  }

  return StatusCode::SUCCESS;
}

//.........这里部分代码省略.........

	// Interface for conversion from Pythia8::Event to HepMC event.
	HepMC::Pythia8ToHepMC ToHepMC;

	std::size_t counter = 0; // number of "interesting" (that satisfy all the requirements) events generated so far
	std::size_t total = 0; // total number of events generated so far

	if(verbose) {
		std::cout << "Starting to generate events" << std::endl;
	}

	std::map<std::size_t, std::size_t> stable_ptcs_count;

	while(counter < nevents) {
		if(pythia.next()) {
			++total;

			// creating HepMC event storage
			HepMC::GenEvent * hepmcevt = new HepMC::GenEvent(HepMC::Units::GEV, HepMC::Units::MM);

			// converting generated event to HepMC format
			ToHepMC.fill_next_event(pythia, hepmcevt);

			auto nstable = std::count_if(hepmcevt->particles_begin(), hepmcevt->particles_end(), [](HepMC::GenParticle const * const ptc_ptr) {return ptc_ptr->status() == 1;});

			if(nstable <= 7) {
				stable_ptcs_count[nstable]++;
				++counter;

				if(verbose && counter % 100 == 0) {
					std::cout << counter << " events with with 7 or less particles in the final state have been generated (" << total << " total). " << std::chrono::duration<double>(std::chrono::system_clock::now() - last_timestamp).count() / 100 << "events / sec" << std::endl;
					last_timestamp = std::chrono::system_clock::now();
				}

				// filling event info
				auto evinfo = fcc::EventInfo();
				evinfo.Number(counter); // Number takes int as its parameter, so here's a narrowing conversion (std::size_t to int). Should be safe unless we get 2^32 events or more. Then undefined behaviour
				evinfocoll.push_back(evinfo);

				// filling vertices
				std::unordered_map<HepMC::GenVertex *, fcc::GenVertex> vtx_map;
				for(auto iv = hepmcevt->vertices_begin(), endv = hepmcevt->vertices_end(); iv != endv; ++iv) {
					auto vtx = fcc::GenVertex();
					vtx.Position().X = (*iv)->position().x();
					vtx.Position().Y = (*iv)->position().y();
					vtx.Position().Z = (*iv)->position().z();
					vtx.Ctau((*iv)->position().t());
					vtx_map.emplace(*iv, vtx);

					vcoll.push_back(vtx);
				}

				// filling particles
				for(auto ip = hepmcevt->particles_begin(), endp = hepmcevt->particles_end(); ip != endp; ++ip) {
					auto ptc = fcc::MCParticle();
					auto & core = ptc.Core();
					core.Type = (*ip)->pdg_id();
					core.Status = (*ip)->status();

					core.Charge = pythia.particleData.charge(core.Type); // PYTHIA returns charge as a double value (in case it's quark), so here's a narrowing conversion (double to int), but here it's safe
					core.P4.Mass = (*ip)->momentum().m();
					core.P4.Px = (*ip)->momentum().px();
					core.P4.Py = (*ip)->momentum().py();
					core.P4.Pz = (*ip)->momentum().pz();

					auto prodvtx = vtx_map.find((*ip)->production_vertex());
					if(prodvtx != vtx_map.end()) {
						ptc.StartVertex(prodvtx->second);
					}
					auto endvtx = vtx_map.find((*ip)->end_vertex());
					if(endvtx != vtx_map.end()) {
						ptc.EndVertex(endvtx->second);
					}

					pcoll.push_back(ptc);
				}

				writer.writeEvent();
				store.clearCollections();
			}

			// freeing resources
			if(hepmcevt) {
				delete hepmcevt;
				hepmcevt = nullptr;
			}
		}
	}

	writer.finish();

	std::cout << counter << " events with 7 or less particles in the final state have been generated (" << total << " total)." << std::endl;
	for(auto const & nv : stable_ptcs_count) {
		std::cout << std::setw(4) << std::right << nv.first << std::setw(4) << std::right << nv.second << "(" << static_cast<long double>(nv.second) * static_cast<long double>(100) / static_cast<long double>(total) << "%)" << std::endl;
	}
	auto elapsed_seconds = std::chrono::duration<double>(std::chrono::system_clock::now() - start_time).count();
	std::cout << "Elapsed time: " << elapsed_seconds << " s (" << static_cast<long double>(counter) / static_cast<long double>(elapsed_seconds) << " events / s)" << std::endl;

	return EXIT_SUCCESS;
}