36#include "G4Svc/G4HepMCInterface.h"
39#include "G4LorentzVector.hh"
41#include "G4PrimaryParticle.hh"
42#include "G4PrimaryVertex.hh"
43#include "Randomize.hh"
45#include "G4Svc/IG4Svc.h"
46#include "G4Svc/G4Svc.h"
47#include "GaudiKernel/IDataProviderSvc.h"
48#include "GaudiKernel/AlgFactory.h"
49#include "GaudiKernel/MsgStream.h"
50#include "GaudiKernel/SvcFactory.h"
51#include "GaudiKernel/ISvcLocator.h"
52#include "GaudiKernel/SmartDataPtr.h"
53#include "GaudiKernel/Bootstrap.h"
61 ISvcLocator* svcLocator = Gaudi::svcLocator();
62 StatusCode status = svcLocator->service(
"RealizationSvc",tmpReal);
63 if (!status.isSuccess())
65 std::cout <<
" Could not initialize Realization Service in G4HepMCInterface" << std::endl;
81 for(HepMC::GenEvent::vertex_const_iterator vitr= hepmcevt->vertices_begin();
82 vitr != hepmcevt->vertices_end(); ++vitr ) {
83 G4cout<<G4endl<<
"vertex:"<<i<<
" barcode:"<<(*vitr)->barcode()<<
" ";
90 xvtx.setX( (*vitr)-> position().
x());
91 xvtx.setY( (*vitr)-> position().y());
92 xvtx.setZ( (*vitr)-> position().z());
93 xvtx.setT( (*vitr)-> position().
t());
94 G4cout<<
"x:"<<xvtx.x()<<
" y:"<<xvtx.y()<<
" z:"<<xvtx.z()
95 <<
" t:"<<xvtx.t()*mm/c_light<<G4endl;
98 for (HepMC::GenVertex::particle_iterator
99 pitr= (*vitr)->particles_begin(HepMC::children);
100 pitr != (*vitr)->particles_end(HepMC::children); ++pitr) {
101 G4int pdgcode= (*pitr)-> pdg_id();
107 G4cout<<
"particle:"<<j<<
" pdgcode:"<<pdgcode<<
" ";
108 G4cout<<
"p:"<<p.x()<<
" "<<p.y()<<
" "<<p.z()<<
" ";
109 if((*pitr)->end_vertex())
110 G4cout<<
"endVtx:"<<(*pitr)->end_vertex()->barcode()<<
" ";
111 G4cout<<
"status:"<<(*pitr)->status()<<G4endl;
120 G4cout<<
"particles sent to Geant4: "<<G4endl;
121 for(
size_t ii=0; ii<
HPlist.size(); ii++ )
123 G4int pdgcode =
HPlist[ii]->GetTheParticle()->GetPDGcode();
124 G4cout<<pdgcode<<
" ";
132 for(HepMC::GenEvent::vertex_const_iterator vitr= hepmcevt->vertices_begin();
133 vitr != hepmcevt->vertices_end(); ++vitr ) {
134 for (HepMC::GenVertex::particle_iterator
135 pitr= (*vitr)->particles_begin(HepMC::children);
136 pitr != (*vitr)->particles_end(HepMC::children); ++pitr) {
137 if((*pitr)->end_vertex())
145 G4cout<<G4endl<<
"generator is GENBES!"<<G4endl;
147 G4cout<<G4endl<<
"generator is not GENBES!"<<G4endl;
167 G4LorentzVector ptot=0;
168 double E_cms,p_Mag,e_Mass2,M_cms,theta;
169 for(HepMC::GenEvent::vertex_const_iterator vitr= hepmcevt->vertices_begin();
170 vitr != hepmcevt->vertices_end(); ++vitr ) {
172 for(HepMC::GenVertex::particle_iterator
173 vpitr = (*vitr)->particles_begin(HepMC::children);
174 vpitr != (*vitr)->particles_end(HepMC::children); ++vpitr) {
176 if((*vpitr)->status() !=1)
continue;
189 ISvcLocator* svcLocator = Gaudi::svcLocator();
192 StatusCode sc=svcLocator->service(
"G4Svc", tmpSvc);
193 m_G4Svc=
dynamic_cast<G4Svc *
>(tmpSvc);
202 e_Mass2=electron_mass_c2*electron_mass_c2;
203 p_Mag=sqrt((E_cms*E_cms-4*e_Mass2)/(2*(1-
cos(
pi-2*theta))));
204 G4ThreeVector p_Positron(p_Mag*
sin(theta),0,p_Mag*
cos(theta));
205 G4ThreeVector p_Electron(p_Mag*
sin(
pi-theta),0,p_Mag*
cos(
pi-theta));
206 G4ThreeVector vee=p_Positron+p_Electron;
208 if(m_G4Svc-> GetSetBeamShift()==
true && fabs(M_cms-3686)<50.0) {vee = beamshift;}
211 M_cms = sqrt(M_cms*M_cms + pee*pee);
215 for(HepMC::GenEvent::vertex_const_iterator vitr= hepmcevt->vertices_begin();
216 vitr != hepmcevt->vertices_end(); ++vitr ) {
219 G4LorentzVector xvtx;
220 xvtx.setX( (*vitr)-> position().
x());
221 xvtx.setY( (*vitr)-> position().y());
222 xvtx.setZ( (*vitr)-> position().z());
223 xvtx.setT( (*vitr)-> position().
t());
225 (*vitr)->set_position(xvtx);
229 for (HepMC::GenVertex::particle_iterator
230 vpitr = (*vitr)->particles_begin(HepMC::children);
231 vpitr != (*vitr)->particles_end(HepMC::children); ++vpitr) {
238 (*vpitr)->set_momentum(p);
253 ISvcLocator* svcLocator = Gaudi::svcLocator();
256 StatusCode sc=svcLocator->service(
"G4Svc", tmpSvc);
257 m_G4Svc=
dynamic_cast<G4Svc *
>(tmpSvc);
265 if(m_logLevel <=2 && m_G4Svc->GetBoostLab()==
true)
272 for(HepMC::GenEvent::vertex_const_iterator vitr= hepmcevt->vertices_begin();
273 vitr != hepmcevt->vertices_end(); ++vitr )
276 G4int pOut = (*vitr)->particles_out_size();
277 HepMC::GenVertex::particle_iterator pitr = (*vitr)->particles_begin(HepMC::children);
278 G4int pdgcode= (*pitr)-> pdg_id();
279 if(pOut == 1 && pdgcode == -11)
281 G4int barcodeVtx = (*vitr)->barcode();
283 for (HepMC::GenVertex::particle_iterator pitr= (*vitr)->particles_begin(HepMC::children);
284 pitr != (*vitr)->particles_end(HepMC::children); ++pitr)
286 G4int pdgcode = (*pitr)->pdg_id();
289 if(pdgcode==-22) pdgcode=22;
295 G4PrimaryParticle* particle =
new G4PrimaryParticle(pdgcode, p.x()*GeV, p.y()*GeV, p.z()*GeV);
298 G4int ISTHEP = (*pitr)->status();
300 G4int barcodeEndVtx = 0;
301 if((*pitr)->end_vertex())
302 barcodeEndVtx = (*pitr)->end_vertex()->barcode();
305 HPlist.push_back(hepmcParticle);
307 for(
size_t ii=0; ii<
HPlist.size(); ii++ )
309 if(barcodeVtx ==
HPlist[ii]->GetBarcodeEndVtx())
311 HPlist[ii]->GetTheParticle()->SetDaughter(particle);
312 hepmcParticle->
Done();
322 for(HepMC::GenEvent::vertex_const_iterator vitr= hepmcevt->vertices_begin();
323 vitr != hepmcevt->vertices_end(); ++vitr )
325 for (HepMC::GenVertex::particle_iterator pitr= (*vitr)->particles_begin(HepMC::children);
326 pitr != (*vitr)->particles_end(HepMC::children); ++pitr)
328 G4int ISTHEP = (*pitr)->status();
335 G4int pdgcode = (*pitr)->pdg_id();
336 G4int barcodeEndVtx = 0;
337 if((*pitr)->end_vertex())
338 barcodeEndVtx = (*pitr)->end_vertex()->barcode();
339 G4PrimaryParticle* particle =
new G4PrimaryParticle(pdgcode, p.x()*GeV, p.y()*GeV, p.z()*GeV);
344 HPlist.push_back(hepmcParticle);
346 hepmcParticle->
Done();
356 G4double pmPosX,pmPosY,pmPosZ,pmTime;
357 G4double tmpPosX,tmpPosY,tmpPosZ,tmpT;
358 G4double beamPosX,beamPosY,beamPosZ,beamSizeX,beamSizeY,beamSizeZ;
360 if(m_RealizationSvc->
UseDBFlag() ==
false) {
388 G4double gaussX = G4RandGauss::shoot();
389 G4double gaussY = G4RandGauss::shoot();
390 G4double gaussZ = G4RandGauss::shoot();
391 G4double gaussT = G4RandGauss::shoot();
398 G4double ran=G4UniformRand();
399 G4double beamTime = bunchTimeSigma*G4RandGauss::shoot() + beamStartTime + beamDeltaTime*int(ran*nBunch);
401 tmpPosX = (beamPosX + beamSizeX*gaussX ) *mm;
402 tmpPosY = (beamPosY + beamSizeY*gaussY ) *mm;
403 tmpPosZ = (beamPosZ + beamSizeZ*gaussZ ) *mm;
404 tmpT = (beamSizeZ * gaussT ) * mm/c_light +beamTime;
406 G4LorentzVector tmpv(tmpPosX,tmpPosY,tmpPosZ,tmpT*c_light/mm);
408 HepMC::GenEvent::vertex_const_iterator vitr0= hepmcevt->vertices_begin();
409 G4LorentzVector xvtx0 ;
410 xvtx0.setX( (*vitr0)-> position().
x());
411 xvtx0.setY( (*vitr0)-> position().y());
412 xvtx0.setZ( (*vitr0)-> position().z());
413 xvtx0.setT( (*vitr0)-> position().
t());
414 pmPosX = xvtx0.x()*mm + tmpPosX;
415 pmPosY = xvtx0.y()*mm + tmpPosY;
416 pmPosZ = xvtx0.z()*mm + tmpPosZ;
417 pmTime = xvtx0.t()*mm/c_light + tmpT;
422 G4cout<<xvtx0.x()<<
" "<<xvtx0.y()<<
" "<<xvtx0.z()<<
" "
423 <<beamSizeX*gaussX<<
" "
424 <<beamSizeY*gaussY<<
" "
425 <<beamSizeZ*gaussZ<<
" "<<G4endl;
426 G4cout<<xvtx0.t()* mm/c_light<<
" "
427 <<beamSizeZ * gaussT/sqrt(2)*mm/c_light<<
" "
430 for(HepMC::GenEvent::vertex_const_iterator vitr= hepmcevt->vertices_begin();
431 vitr != hepmcevt->vertices_end(); ++vitr )
433 G4LorentzVector xvtx;
434 xvtx.setX((*vitr)-> position().
x());
435 xvtx.setY((*vitr)-> position().y());
436 xvtx.setZ((*vitr)-> position().z());
437 xvtx.setT((*vitr)-> position().
t());
438 (*vitr)->set_position(xvtx+tmpv);
442 G4PrimaryVertex* g4vtx=
new G4PrimaryVertex(pmPosX,pmPosY,pmPosZ,pmTime);
445 for(
size_t ii=0; ii<
HPlist.size(); ii++ )
447 if(
HPlist[ii]->GetISTHEP() > 0 )
450 G4PrimaryParticle* initialParticle =
HPlist[ii]->GetTheParticle();
451 g4vtx->SetPrimary( initialParticle );
456 for(
size_t iii=0;iii<
HPlist.size();iii++)
460 g4event->AddPrimaryVertex(g4vtx);
468 HepMC::GenEvent* aevent=
new HepMC::GenEvent();
482 G4Exception(
"G4HepMCInterface",
"GeneratePrimaryVertex",RunMustBeAborted,
483 "HepMCInterface: no generated particles. run terminated..." );
**********INTEGER nmxhep !maximum number of particles DOUBLE PRECISION vhep INTEGER jdahep COMMON hepevt $ !serial number $ !number of particles $ !status code $ !particle ident KF $ !parent particles $ !childreen particles $ !four momentum
double sin(const BesAngle a)
double cos(const BesAngle a)
**********Class see also m_nmax DOUBLE PRECISION m_amel DOUBLE PRECISION m_x2 DOUBLE PRECISION m_alfinv DOUBLE PRECISION m_Xenph INTEGER m_KeyWtm INTEGER m_idyfs DOUBLE PRECISION m_zini DOUBLE PRECISION m_q2 DOUBLE PRECISION m_Wt_KF DOUBLE PRECISION m_WtCut INTEGER m_KFfin *COMMON c_KarLud $ !Input CMS energy[GeV] $ !CMS energy after beam spread beam strahlung[GeV] $ !Beam energy spread[GeV] $ !z boost due to beam spread $ !electron beam mass *ff pair spectrum $ !minimum v
virtual ~G4HepMCInterface()
std::vector< G4HepMCParticle * > HPlist
void HepMC2G4(HepMC::GenEvent *hepmcevt, G4Event *g4event)
HepMC::GenEvent * hepmcEvent
void Print(const HepMC::GenEvent *hepmcevt)
void Boost(HepMC::GenEvent *hepmcevt)
virtual void GeneratePrimaryVertex(G4Event *anEvent)
virtual HepMC::GenEvent * GenerateHepMCEvent()
G4int CheckType(const HepMC::GenEvent *hepmcevt)
double GetBeamDeltaTime()
double GetBunchTimeSigma()
double GetBeamStartTime()
void SetBeamTime(double value)