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

/// 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;
}

int main(int argc, char ** argv) {

  std::cout << " -- Starting program..." << std::endl;

  const unsigned nZ = 100;
  const unsigned z0 = 100;
  const unsigned step = 2;

  //TFile *fout = TFile::Open("PLOTS/output_hepmc_vtxorig_MB.root","RECREATE");
  //TFile *fout = TFile::Open("PLOTS/output_hepmc_vtxmodif.root","RECREATE");
  TFile *fout = TFile::Open("PLOTS/output_hepmc_hggmodif.root","RECREATE");
  fout->cd();
  TH1F *hvtx_x = new TH1F("hvtx_x",";x (mm);vertices",1000,-10,10);
  TH1F *hvtx_y = new TH1F("hvtx_y",";y (mm);vertices",1000,-10,10);
  TH1F *hvtx_z = new TH1F("hvtx_z",";z (mm);vertices",1000,-500,500);
  TH1F *hvtx_t = new TH1F("hvtx_t",";t (ns);vertices",20000,-100,100);
  TH2F *hvtx_tvsz = new TH2F("hvtx_tvsz",";z (mm); t (ps);vertices",1000,-200,200,1200,-600,600);
  TH1F *hvtx_t_z[nZ];
  for (unsigned i(0);i<nZ;++i){
    std::ostringstream label;
    int zmin = -1*z0+i*step;
    int zmax = -1*z0+(i+1)*step;
    label << "hvtx_t_z_" << zmin << "_" << zmax;
    hvtx_t_z[i] = new TH1F(label.str().c_str(),";t (ps);vertices",120,-600,600);
  }
  TH1F *hProton = new TH1F("hProton",";p (GeV);particles",1000,0,100);
  TH1F *hNeutron = new TH1F("hNeutron",";p (GeV);particles",1000,0,100);
  TH1F *hPipm = new TH1F("hPipm",";p (GeV);particles",1000,0,100);
  TH1F *hProtonLog = new TH1F("hProtonLog",";log(p) (log(GeV));particles",100,-2,2);
  TH1F *hNeutronLog = new TH1F("hNeutronLog",";log(p) (log(GeV));particles",100,-2,2);
  TH1F *hPipmLog = new TH1F("hPipmLog",";log(p) (log(GeV));particles",100,-2,2);

  const unsigned nFiles = 1;//6;

  for (unsigned i(0);i<nFiles;++i){

    // specify an input file
    std::ostringstream lname;
    //lname << "/afs/cern.ch/work/p/pdauncey/public/Pythia140305_";
    //if (i==0) lname << "000000";
    //else if (i<10) lname << "00000" << i;
    //else if (i<100) lname << "0000" << i;
    //else if (i<1000) lname << "000" << i;
    //else if (i<10000) lname << "00" << i;
    //else if (i<100000) lname << "0" << i;
    //else lname << i;
    //lname << ".dat";
    //lname << "/afs/cern.ch/work/a/amagnan/public/HepMCFiles/ggHgg_origVtx.dat";
    //lname << "/afs/cern.ch/work/a/amagnan/public/HepMCFiles/ggHgg_modifyVtx.dat";
    //lname << "/afs/cern.ch/work/a/amagnan/public/HepMCFiles/vertexHLLHC.dat";
    //lname << "/afs/cern.ch/work/p/pdauncey/public/Pythia140305_000000.dat";
    lname << "/afs/cern.ch/work/a/amagnan/public/HepMCFiles/ggHgg_1428658356.dat";
    HepMC::IO_GenEvent ascii_in(lname.str().c_str(),std::ios::in);
    
    // get the first event
    HepMC::GenEvent* evt = ascii_in.read_next_event();
    // loop until we run out of events
    while ( evt ) {
      unsigned ievt =  evt->event_number();
      if (ievt%10000==0) std::cout << "Processing Event Number "
				 << ievt
				 << std::endl;


      GenVertex * parent = 0;
      HepMC::GenEvent::vertex_const_iterator q = evt->vertices_begin();

      //for (; q != evt->vertices_end(); ++q ){

	//if ((*q)->position().x()!=0 || (*q)->position().y()!=0) continue;

      double z = (*q)->position().z();
      double t = (*q)->position().t();
      hvtx_x->Fill((*q)->position().x());
      hvtx_y->Fill((*q)->position().y());
      hvtx_z->Fill(z);
      hvtx_t->Fill(t);
      //hvtx_tvsz->Fill(z,t*1000);
      hvtx_tvsz->Fill(z,t);

      if (fabs(z)<z0){
	unsigned idx = static_cast<unsigned>((z+z0)*1./step);
	if (idx>(nZ-1)) continue;
	hvtx_t_z[idx]->Fill(t*1000);
      }

	/*std::cout << " -- vtx pos: " << (*q)->position().x() << " " << (*q)->position().y() << " " << (*q)->position().z() << " nParticles: in=" << (*q)->particles_in_size() << " " << (*q)->particles_out_size()
		  << std::endl;
	for ( HepMC::GenVertex::particles_in_const_iterator p
		= (*q)->particles_in_const_begin(); p != (*q)->particles_in_const_end(); ++p ){

	  std::cout << " ---- in particle " << (*p)->pdg_id() << " status " << (*p)->status()
		    << std::endl;
	  
	}
 	for ( HepMC::GenVertex::particles_out_const_iterator p
		= (*q)->particles_out_const_begin(); p != (*q)->particles_out_const_end(); ++p ){

	  std::cout << " ---- out particle " << (*p)->pdg_id() << " status " << (*p)->status()
		    << std::endl;
//.........这里部分代码省略.........

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

	// 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;
}