BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
G4HepMCInterface.cpp
Go to the documentation of this file.
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26// ====================================================================
27//
28// G4HepMCInterface.cc
29// $Id: G4HepMCInterface.cpp,v 1.24 2022/04/02 02:46:19 fucd Exp $
30//
31// ====================================================================
32
33#ifndef WIN32 // Temporarly disabled on Windows, until CLHEP
34 // will support the HepMC module
35
36# include "G4HepMCInterface.h"
37
38// # include "EvtGenBase/EvtPDL.hh"
39# include "G4Event.hh"
40# include "G4LorentzVector.hh"
41# include "G4ParticleTable.hh"
42# include "G4PrimaryParticle.hh"
43# include "G4PrimaryVertex.hh"
44# include "HepPDT/ParticleData.hh"
45# include "HepPDT/ParticleDataTable.hh"
46# include "Randomize.hh"
47# include "globals.hh"
48
49# include "G4PhysicalConstants.hh"
50# include "G4Svc.h"
51# include "G4Svc/IG4Svc.h"
52# include "G4SystemOfUnits.hh"
53# include "GaudiKernel/Bootstrap.h"
54# include "GaudiKernel/IAlgManager.h"
55# include "GaudiKernel/IDataProviderSvc.h"
56# include "GaudiKernel/IPartPropSvc.h"
57# include "GaudiKernel/ISvcLocator.h"
58# include "GaudiKernel/MsgStream.h"
59# include "GaudiKernel/SmartDataPtr.h"
60# include <errno.h>
61using namespace CLHEP;
62
63////////////////////////////////////
64G4HepMCInterface::G4HepMCInterface()
65 : hepmcEvent( 0 )
66// , m_geantinoPdgId( 100 )
67// , m_chargedgeantinoPdgId( 101 )
68////////////////////////////////////
69{
70
71 ISvcLocator* svcLocator = Gaudi::svcLocator();
72 StatusCode status = svcLocator->service( "RealizationSvc", m_RealizationSvc );
73 if ( !status.isSuccess() )
74 {
75 std::cout << " Could not initialize Realization Service in G4HepMCInterface" << std::endl;
76 }
77
78 // int pdlId1 = -1, pdlId2 = -1;
79 // auto algMan = svcLocator->as<IAlgManager>();
80 // if ( algMan )
81 // {
82 // if ( algMan->existsAlgorithm( "EvtDecay" ) )
83 // { // BesEvtGen used
84 // EvtId evtId1 = EvtPDL::getId( "geantino" );
85 // pdlId1 = evtId1.getId();
86 // if ( pdlId1 != -1 )
87 // { // denote geantino has been defined in pdt.table
88 // m_geantinoPdgId = EvtPDL::getStdHep( evtId1 );
89 // }
90 // EvtId evtId2 = EvtPDL::getId( "chargedgeantino" );
91 // pdlId2 = evtId2.getId();
92 // if ( pdlId2 != -1 )
93 // { // denote chargedgeantino has been defined in pdt.table
94 // m_chargedgeantinoPdgId = EvtPDL::getStdHep( evtId2 );
95 // }
96 // }
97 // }
98
99 // if ( pdlId1 == -1 || pdlId2 == -1 )
100 // { // generator is not BesEvtGen, defined in PDGTABLE.MeV; if not, will use default
101 // IPartPropSvc* p_PartPropSvc;
102 // StatusCode PartPropStatus = svcLocator->service( "PartPropSvc", p_PartPropSvc );
103 // if ( !PartPropStatus.isSuccess() || 0 == p_PartPropSvc )
104 // {
105 // std::cout << " Could not initialize Particle Properties Service in G4HepMCInterface"
106 // << std::endl;
107 // }
108 // else
109 // {
110 // HepPDT::ParticleDataTable* pdt = p_PartPropSvc->PDT();
111 // if ( pdlId1 == -1 )
112 // {
113 // const HepPDT::ParticleData* pd = pdt->particle( "geantino" );
114 // if ( pd ) m_geantinoPdgId = pd->pid();
115 // }
116 // if ( pdlId2 == -1 )
117 // {
118 // const HepPDT::ParticleData* pd = pdt->particle( "chargedgeantino" );
119 // if ( pd ) m_chargedgeantinoPdgId = pd->pid();
120 // }
121 // }
122 // }
123}
124
125/////////////////////////////////////
126G4HepMCInterface::~G4HepMCInterface()
127/////////////////////////////////////
128{
129 delete hepmcEvent;
130}
131void G4HepMCInterface::Print( const HepMC::GenEvent* hepmcevt ) {
132 G4int i = 0;
133 for ( HepMC::GenEvent::vertex_const_iterator vitr = hepmcevt->vertices_begin();
134 vitr != hepmcevt->vertices_end(); ++vitr )
135 { // loop for vertex ...
136 G4cout << G4endl << "vertex:" << i << " barcode:" << ( *vitr )->barcode() << " ";
137 i++;
138 // if((*vitr)->mother())
139 // G4cout<<"mother particle: "<<(*vitr)->mother()-> pdg_id()<<" ";
140 // if((*vitr)->secondMother())
141 // G4cout<<" second mother: "<<(*vitr)->secondMother()-> pdg_id()<<G4endl;
142 G4LorentzVector xvtx;
143 xvtx.setX( ( *vitr )->position().x() );
144 xvtx.setY( ( *vitr )->position().y() );
145 xvtx.setZ( ( *vitr )->position().z() );
146 xvtx.setT( ( *vitr )->position().t() );
147 G4cout << "x:" << xvtx.x() << " y:" << xvtx.y() << " z:" << xvtx.z()
148 << " t:" << xvtx.t() * mm / c_light << G4endl;
149
150 G4int j = 0;
151 for ( HepMC::GenVertex::particle_iterator pitr =
152 ( *vitr )->particles_begin( HepMC::children );
153 pitr != ( *vitr )->particles_end( HepMC::children ); ++pitr )
154 {
155 G4int pdgcode = ( *pitr )->pdg_id();
156 G4LorentzVector p;
157 p.setX( ( *pitr )->momentum().x() );
158 p.setY( ( *pitr )->momentum().y() );
159 p.setZ( ( *pitr )->momentum().z() );
160 p.setT( ( *pitr )->momentum().t() );
161 G4cout << "particle:" << j << " pdgcode:" << pdgcode << " ";
162 G4cout << "p:" << p.x() << " " << p.y() << " " << p.z() << " ";
163 if ( ( *pitr )->end_vertex() )
164 G4cout << "endVtx:" << ( *pitr )->end_vertex()->barcode() << " ";
165 G4cout << "status:" << ( *pitr )->status() << G4endl;
166 j++;
167 }
168 }
169}
170
171void G4HepMCInterface::PrintHPlist() {
172 G4cout << G4endl;
173 G4cout << "particles sent to Geant4: " << G4endl;
174 for ( size_t ii = 0; ii < HPlist.size(); ii++ )
175 {
176 G4int pdgcode = HPlist[ii]->GetTheParticle()->GetPDGcode();
177 G4cout << pdgcode << " ";
178 }
179 G4cout << G4endl;
180}
181
182G4int G4HepMCInterface::CheckType( const HepMC::GenEvent* hepmcevt ) {
183 G4int flag = 0;
184 for ( HepMC::GenEvent::vertex_const_iterator vitr = hepmcevt->vertices_begin();
185 vitr != hepmcevt->vertices_end(); ++vitr )
186 { // loop for vertex ...
187 for ( HepMC::GenVertex::particle_iterator pitr =
188 ( *vitr )->particles_begin( HepMC::children );
189 pitr != ( *vitr )->particles_end( HepMC::children ); ++pitr )
190 {
191 if ( ( *pitr )->end_vertex() )
192 {
193 flag = 1;
194 break;
195 }
196 }
197 if ( flag ) break;
198 }
199 if ( m_logLevel <= 4 )
200 {
201 if ( flag == 0 ) G4cout << G4endl << "generator is GENBES!" << G4endl;
202 else G4cout << G4endl << "generator is not GENBES!" << G4endl;
203 }
204 return flag;
205}
206
207// Panmh 2006.12
208////////////////////////////////////////////////////////////////
209void G4HepMCInterface::Boost( HepMC::GenEvent* hepmcevt )
210////////////////////////////////////////////////////////////////
211{
212 // suppose relavtive vectoy is v(vx,vy,vz),position are (x,y,z,t)in CMS and(X,Y,Z,T)in
213 // Lab,the formulae for boost vertex position from Cms to Lab in natural units are following:
214 // v2 = vx*vx + vy*vy + vz*vz;
215 // gamma = 1.0 / sqrt(1.0 - v2);
216 // bp = vx*x + vy*y + vz*z;
217 // X=x + gamma2*bp*vx + gamma*vx*t;
218 // Y=y + gamma2*bp*vy + gamma*vy*t;
219 // Z=z + gamma2*bp*vz + gamma*vz*t;
220 // T=gamma*(t + bp);
221 // the function of these formulae is the same as xvtx.boost(ThreeVector v)
222 // For the E_cms,need to loop and splice the momentum of all the final state particles
223 G4LorentzVector ptot( 0, 0, 0, 0 );
224 double E_cms, p_Mag, e_Mass2, M_cms, theta;
225 for ( HepMC::GenEvent::vertex_const_iterator vitr = hepmcevt->vertices_begin();
226 vitr != hepmcevt->vertices_end(); ++vitr )
227 { // loop for vertex ...
228
229 for ( HepMC::GenVertex::particle_iterator vpitr =
230 ( *vitr )->particles_begin( HepMC::children );
231 vpitr != ( *vitr )->particles_end( HepMC::children ); ++vpitr )
232 {
233
234 if ( ( *vpitr )->status() != 1 ) continue; // Not final state particle
235
236 G4LorentzVector p;
237 p.setX( ( *vpitr )->momentum().x() );
238 p.setY( ( *vpitr )->momentum().y() );
239 p.setZ( ( *vpitr )->momentum().z() );
240 p.setT( ( *vpitr )->momentum().t() );
241 ptot = p + ptot;
242 }
243 }
244 E_cms = ptot.e() * GeV;
245
246 // get G4Svc, allow user to access the beam momentum shift in JobOptions
247 ISvcLocator* svcLocator = Gaudi::svcLocator();
248 IG4Svc* tmpSvc;
249 G4Svc* m_G4Svc;
250 StatusCode sc = svcLocator->service( "G4Svc", tmpSvc );
251 m_G4Svc = dynamic_cast<G4Svc*>( tmpSvc );
252
253 G4ThreeVector v( 0, 0, 0 ); // velocity of cms as seen from lab
254 // for cms velocity
255 theta = 11 * mrad;
256 // theta=(m_G4Svc->GetBeamAngle())*mrad;
257 // theta is half of the angle between the two beams,namely,is the angle between the beam and
258 // Z axis.
259 M_cms = E_cms; // for p_cms=0 in the CMS
260 e_Mass2 = electron_mass_c2 * electron_mass_c2; // square of electron mass
261 p_Mag = sqrt( ( E_cms * E_cms - 4 * e_Mass2 ) / ( 2 * ( 1 - cos( pi - 2 * theta ) ) ) );
262 G4ThreeVector p_Positron( p_Mag * sin( theta ), 0, p_Mag * cos( theta ) );
263 G4ThreeVector p_Electron( p_Mag * sin( pi - theta ), 0, p_Mag * cos( pi - theta ) );
264 G4ThreeVector vee = p_Positron + p_Electron;
265 G4ThreeVector beamshift( m_G4Svc->GetBeamShiftPx(), m_G4Svc->GetBeamShiftPy(),
266 m_G4Svc->GetBeamShiftPz() );
267 if ( m_G4Svc->GetSetBeamShift() == true && fabs( M_cms - 3686 ) < 50.0 ) { vee = beamshift; }
268
269 double pee = vee.r();
270 M_cms = sqrt( M_cms * M_cms + pee * pee );
271 v = vee / M_cms;
272
273 for ( HepMC::GenEvent::vertex_const_iterator vitr = hepmcevt->vertices_begin();
274 vitr != hepmcevt->vertices_end(); ++vitr )
275 { // loop for vertex ...
276
277 // Boost vertex position from cms to lab
278 G4LorentzVector xvtx;
279 xvtx.setX( ( *vitr )->position().x() );
280 xvtx.setY( ( *vitr )->position().y() );
281 xvtx.setZ( ( *vitr )->position().z() );
282 xvtx.setT( ( *vitr )->position().t() );
283 xvtx.boost( v );
284 ( *vitr )->set_position( xvtx );
285
286 // Boost the particles four momentum from cms to lab
287 // the transform formulae are similary to the position boost
288 for ( HepMC::GenVertex::particle_iterator vpitr =
289 ( *vitr )->particles_begin( HepMC::children );
290 vpitr != ( *vitr )->particles_end( HepMC::children ); ++vpitr )
291 {
292 G4LorentzVector p;
293 p.setX( ( *vpitr )->momentum().x() );
294 p.setY( ( *vpitr )->momentum().y() );
295 p.setZ( ( *vpitr )->momentum().z() );
296 p.setT( ( *vpitr )->momentum().t() );
297 p = p.boost( v );
298 ( *vpitr )->set_momentum( p );
299 }
300 }
301}
302
303////////////////////////////////////////////////////////////////
304void G4HepMCInterface::HepMC2G4( HepMC::GenEvent* hepmcevt, G4Event* g4event )
305////////////////////////////////////////////////////////////////
306{
307 // Print HepMC Event before boost
308 if ( m_logLevel <= 4 ) Print( hepmcevt );
309
310 // get G4Svc
311 ISvcLocator* svcLocator = Gaudi::svcLocator();
312 IG4Svc* tmpSvc;
313 G4Svc* m_G4Svc;
314 StatusCode sc = svcLocator->service( "G4Svc", tmpSvc );
315 m_G4Svc = dynamic_cast<G4Svc*>( tmpSvc );
316
317 // considering 22rmad beam clossing angle ,need to boost the Cms to Lab
318 // need to boost ?
319 if ( m_G4Svc->GetBoostLab() == true ) Boost( hepmcevt );
320
321 // Print HepMC Event after boost
322 if ( m_logLevel <= 2 && m_G4Svc->GetBoostLab() == true ) Print( hepmcevt );
323
324 //*********************send particles to HPlist*************************//
325 G4int flag = CheckType( hepmcevt );
326 if ( flag )
327 {
328 for ( HepMC::GenEvent::vertex_const_iterator vitr = hepmcevt->vertices_begin();
329 vitr != hepmcevt->vertices_end(); ++vitr )
330 { // loop for vertex ...
331 G4int vtxFlag = 0;
332 G4int pOut = ( *vitr )->particles_out_size();
333 HepMC::GenVertex::particle_iterator pitr = ( *vitr )->particles_begin( HepMC::children );
334 G4int pdgcode = ( *pitr )->pdg_id();
335 if ( pOut == 1 && abs( pdgcode ) == 11 ) vtxFlag = 1;
336 G4int barcodeVtx = ( *vitr )->barcode();
337
338 for ( HepMC::GenVertex::particle_iterator pitr =
339 ( *vitr )->particles_begin( HepMC::children );
340 pitr != ( *vitr )->particles_end( HepMC::children ); ++pitr )
341 { // loop for particle from this vertex
342 G4int pdgcode = ( *pitr )->pdg_id();
343 if ( vtxFlag == 0 )
344 {
345 if ( pdgcode == -22 ) pdgcode = 22;
346 G4LorentzVector p;
347 p.setX( ( *pitr )->momentum().x() );
348 p.setY( ( *pitr )->momentum().y() );
349 p.setZ( ( *pitr )->momentum().z() );
350 p.setT( ( *pitr )->momentum().t() );
351 G4PrimaryParticle* particle = 0;
352 // geantino and chargedgeantino can't create through pdgcode;
353 // their PdgID = 0 in Geant4.10.07.p02, and if particle->SetPDGcode(100), Geant4 will
354 // not know. So reset PdgID of input HepMC::GenParticle to G4ParticleDefinition's; if
355 // not, BesSensitiveManager::UpdatePrimaryTrack(const G4Track* track) will exit(1)
356 // because MatchError.
357 // if ( pdgcode == m_geantinoPdgId )
358 // {
359 // G4ParticleDefinition* pd =
360 // G4ParticleTable::GetParticleTable()->FindParticle( "geantino" );
361 // particle = new G4PrimaryParticle( pd, p.x() * GeV, p.y() * GeV, p.z() * GeV );
362 // ( *pitr )->set_pdg_id( particle->GetPDGcode() );
363 // }
364 // else if ( pdgcode == m_chargedgeantinoPdgId )
365 // {
366 // G4ParticleDefinition* pd =
367 // G4ParticleTable::GetParticleTable()->FindParticle( "chargedgeantino" );
368 // particle = new G4PrimaryParticle( pd, p.x() * GeV, p.y() * GeV, p.z() * GeV );
369 // ( *pitr )->set_pdg_id( particle->GetPDGcode() );
370 // }
371 // pdgcode==0 dot define both in Geant4 and pdt.table
372 // else if(pdgcode==0){
373 //}
374 // else
375 // {
376 particle = new G4PrimaryParticle( pdgcode, p.x() * GeV, p.y() * GeV, p.z() * GeV );
377 // }
378 // if ( !particle )
379 // { G4cout << "ERROR! Cannot create particle PDGCode=" << pdgcode << G4endl; }
380 // G4PrimaryParticle* particle = new G4PrimaryParticle(pdgcode);
381 // particle->Set4Momentum(p.x()*GeV,p.y()*GeV,p.z()*GeV,p.t()*GeV);
382 G4int ISTHEP = ( *pitr )->status();
383
384 G4int barcodeEndVtx = 0;
385 if ( ( *pitr )->end_vertex() ) barcodeEndVtx = ( *pitr )->end_vertex()->barcode();
386
387 G4HepMCParticle* hepmcParticle =
388 new G4HepMCParticle( particle, ISTHEP, barcodeEndVtx );
389 HPlist.push_back( hepmcParticle );
390
391 for ( size_t ii = 0; ii < HPlist.size(); ii++ )
392 {
393 if ( barcodeVtx == HPlist[ii]->GetBarcodeEndVtx() )
394 {
395 HPlist[ii]->GetTheParticle()->SetDaughter( particle );
396 hepmcParticle->Done();
397 break;
398 }
399 }
400 }
401 }
402 }
403 }
404 else
405 {
406 for ( HepMC::GenEvent::vertex_const_iterator vitr = hepmcevt->vertices_begin();
407 vitr != hepmcevt->vertices_end(); ++vitr )
408 {
409 for ( HepMC::GenVertex::particle_iterator pitr =
410 ( *vitr )->particles_begin( HepMC::children );
411 pitr != ( *vitr )->particles_end( HepMC::children ); ++pitr )
412 {
413 G4int ISTHEP = ( *pitr )->status();
414 G4LorentzVector p;
415 p.setX( ( *pitr )->momentum().x() );
416 p.setY( ( *pitr )->momentum().y() );
417 p.setZ( ( *pitr )->momentum().z() );
418 p.setT( ( *pitr )->momentum().t() );
419
420 G4int pdgcode = ( *pitr )->pdg_id();
421 G4int barcodeEndVtx = 0;
422 if ( ( *pitr )->end_vertex() ) barcodeEndVtx = ( *pitr )->end_vertex()->barcode();
423 G4PrimaryParticle* particle = 0;
424 // geantino and chargedgeantino can't create through pdgcode;
425 // their PdgID = 0 in Geant4.10.07.p02, and if particle->SetPDGcode(100), Geant4 will
426 // not know. So reset PdgID of input HepMC::GenParticle to G4ParticleDefinition's; if
427 // not, BesSensitiveManager::UpdatePrimaryTrack(const G4Track* track) will exit(1)
428 // because MatchError.
429 // if ( pdgcode == m_geantinoPdgId )
430 // {
431 // G4ParticleDefinition* pd =
432 // G4ParticleTable::GetParticleTable()->FindParticle( "geantino" );
433 // particle = new G4PrimaryParticle( pd, p.x() * GeV, p.y() * GeV, p.z() * GeV );
434 // ( *pitr )->set_pdg_id( particle->GetPDGcode() );
435 // }
436 // else if ( pdgcode == m_chargedgeantinoPdgId )
437 // {
438 // G4ParticleDefinition* pd =
439 // G4ParticleTable::GetParticleTable()->FindParticle( "chargedgeantino" );
440 // particle = new G4PrimaryParticle( pd, p.x() * GeV, p.y() * GeV, p.z() * GeV );
441 // ( *pitr )->set_pdg_id( particle->GetPDGcode() );
442 // }
443 // pdgcode==0 dot define both in Geant4 and pdt.table
444 // else if(pdgcode==0){
445 //}
446 // else
447 // {
448 particle = new G4PrimaryParticle( pdgcode, p.x() * GeV, p.y() * GeV, p.z() * GeV );
449 // }
450 // if ( !particle )
451 // { G4cout << "ERROR! Cannot create particle PDGCode=" << pdgcode << G4endl; }
452 // G4PrimaryParticle* particle = new G4PrimaryParticle(pdgcode);
453 // particle->Set4Momentum(p.x()*GeV,p.y()*GeV,p.z()*GeV,p.t()*GeV);
454
455 G4HepMCParticle* hepmcParticle =
456 new G4HepMCParticle( particle, ISTHEP, barcodeEndVtx );
457 HPlist.push_back( hepmcParticle );
458 if ( ISTHEP > 1 ) hepmcParticle->Done();
459 }
460 }
461 }
462
463 // print particles in HPlist
464 if ( m_logLevel <= 4 ) PrintHPlist();
465
466 // get time and position information from G4Svc
467 G4double pmPosX, pmPosY, pmPosZ, pmTime;
468 G4double tmpPosX, tmpPosY, tmpPosZ, tmpT;
469 G4double beamPosX = 0, beamPosY = 0, beamPosZ = 0, beamSizeX = 0, beamSizeY = 0,
470 beamSizeZ = 0;
471
472 if ( m_RealizationSvc->UseDBFlag() == false )
473 {
474 beamPosX = m_G4Svc->GetBeamPosX();
475 beamPosY = m_G4Svc->GetBeamPosY();
476 beamPosZ = m_G4Svc->GetBeamPosZ();
477
478 beamSizeX = m_G4Svc->GetBeamSizeX();
479 beamSizeY = m_G4Svc->GetBeamSizeY();
480 beamSizeZ = m_G4Svc->GetBeamSizeZ() / sqrt( 2 );
481 }
482 if ( m_RealizationSvc->UseDBFlag() == true && m_RealizationSvc->ifReadBunch() == true )
483 {
484 beamPosX = m_RealizationSvc->getBunchPosX();
485 beamPosY = m_RealizationSvc->getBunchPosY();
486 beamPosZ = m_RealizationSvc->getBunchPosZ();
487
488 beamSizeX = m_RealizationSvc->getBunchSizeX();
489 beamSizeY = m_RealizationSvc->getBunchSizeY();
490 beamSizeZ = m_RealizationSvc->getBunchSizeZ();
491 }
492 if ( m_RealizationSvc->UseDBFlag() == true && m_RealizationSvc->ifReadBunch() == false )
493 {
494 beamPosX = m_G4Svc->GetBeamPosX();
495 beamPosY = m_G4Svc->GetBeamPosY();
496 beamPosZ = m_G4Svc->GetBeamPosZ();
497
498 beamSizeX = m_G4Svc->GetBeamSizeX();
499 beamSizeY = m_G4Svc->GetBeamSizeY();
500 beamSizeZ = m_G4Svc->GetBeamSizeZ() / sqrt( 2 );
501 }
502
503 G4double gaussX = G4RandGauss::shoot();
504 G4double gaussY = G4RandGauss::shoot();
505 G4double gaussZ = G4RandGauss::shoot();
506 G4double gaussT = G4RandGauss::shoot();
507
508 G4double beamStartTime = m_G4Svc->GetBeamStartTime() * ns;
509 G4double beamDeltaTime = m_G4Svc->GetBeamDeltaTime() * ns;
510 G4double nBunch = m_G4Svc->GetNBunch();
511 G4double bunchTimeSigma = m_G4Svc->GetBunchTimeSigma() * ns;
512
513 G4double ran = G4UniformRand();
514 G4double beamTime = bunchTimeSigma * G4RandGauss::shoot() + beamStartTime +
515 beamDeltaTime * int( ran * nBunch );
516
517 tmpPosX = ( beamPosX + beamSizeX * gaussX ) * mm;
518 tmpPosY = ( beamPosY + beamSizeY * gaussY ) * mm;
519 tmpPosZ = ( beamPosZ + beamSizeZ * gaussZ ) * mm;
520 tmpT = ( beamSizeZ * gaussT ) * mm / c_light + beamTime;
521
522 G4LorentzVector tmpv( tmpPosX, tmpPosY, tmpPosZ, tmpT * c_light / mm );
523
524 HepMC::GenEvent::vertex_const_iterator vitr0 = hepmcevt->vertices_begin();
525 G4LorentzVector xvtx0;
526 xvtx0.setX( ( *vitr0 )->position().x() );
527 xvtx0.setY( ( *vitr0 )->position().y() );
528 xvtx0.setZ( ( *vitr0 )->position().z() );
529 xvtx0.setT( ( *vitr0 )->position().t() );
530 pmPosX = xvtx0.x() * mm + tmpPosX;
531 pmPosY = xvtx0.y() * mm + tmpPosY;
532 pmPosZ = xvtx0.z() * mm + tmpPosZ;
533 pmTime = xvtx0.t() * mm / c_light + tmpT;
534
535 if ( m_logLevel <= 4 )
536 {
537 G4cout << G4endl;
538 G4cout << xvtx0.x() << " " << xvtx0.y() << " " << xvtx0.z() << " " << beamSizeX * gaussX
539 << " " << beamSizeY * gaussY << " " << beamSizeZ * gaussZ << " " << G4endl;
540 G4cout << xvtx0.t() * mm / c_light << " " << beamSizeZ * gaussT / sqrt( 2 ) * mm / c_light
541 << " " << beamTime << G4endl;
542 }
543 for ( HepMC::GenEvent::vertex_const_iterator vitr = hepmcevt->vertices_begin();
544 vitr != hepmcevt->vertices_end(); ++vitr )
545 {
546 G4LorentzVector xvtx;
547 xvtx.setX( ( *vitr )->position().x() );
548 xvtx.setY( ( *vitr )->position().y() );
549 xvtx.setZ( ( *vitr )->position().z() );
550 xvtx.setT( ( *vitr )->position().t() );
551 ( *vitr )->set_position( xvtx + tmpv );
552 }
553 m_G4Svc->SetBeamTime( pmTime );
554
555 G4PrimaryVertex* g4vtx = new G4PrimaryVertex( pmPosX, pmPosY, pmPosZ, pmTime );
556
557 // put initial particles to the vertex
558 for ( size_t ii = 0; ii < HPlist.size(); ii++ )
559 {
560 if ( HPlist[ii]->GetISTHEP() > 0 ) // ISTHEP of daughters had been
561 // set to negative
562 {
563 G4PrimaryParticle* initialParticle = HPlist[ii]->GetTheParticle();
564 g4vtx->SetPrimary( initialParticle );
565 }
566 }
567
568 // clear G4HepMCInterface
569 for ( size_t iii = 0; iii < HPlist.size(); iii++ ) { delete HPlist[iii]; }
570 HPlist.clear();
571
572 g4event->AddPrimaryVertex( g4vtx );
573}
574
575///////////////////////////////////////////////////////
576HepMC::GenEvent* G4HepMCInterface::GenerateHepMCEvent()
577///////////////////////////////////////////////////////
578{
579 HepMC::GenEvent* aevent = new HepMC::GenEvent();
580 return aevent;
581}
582
583//////////////////////////////////////////////////////////////
584void G4HepMCInterface::GeneratePrimaryVertex( G4Event* anEvent )
585//////////////////////////////////////////////////////////////
586{
587 // delete previous event object
588 delete hepmcEvent;
589
590 // generate next event
591 hepmcEvent = GenerateHepMCEvent();
592 if ( !hepmcEvent )
593 {
594 G4Exception( "G4HepMCInterface", "GeneratePrimaryVertex", RunMustBeAborted,
595 "HepMCInterface: no generated particles. run terminated..." );
596 return;
597 }
598 HepMC2G4( hepmcEvent, anEvent );
599}
600
601#endif
pi
double pi
Definition EvtPhokharaDef.hh:17
sin
double sin(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:185
cos
double cos(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:187
v
**********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
Definition KarLud.h:35
G4HepMCInterface::hepmcEvent
HepMC::GenEvent * hepmcEvent
Definition G4HepMCInterface.h:64
G4HepMCInterface::HepMC2G4
void HepMC2G4(HepMC::GenEvent *hepmcevt, G4Event *g4event)
Definition G4HepMCInterface.cpp:304
G4HepMCInterface::Print
void Print(const HepMC::GenEvent *hepmcevt)
Definition G4HepMCInterface.cpp:131
G4HepMCInterface::Boost
void Boost(HepMC::GenEvent *hepmcevt)
Definition G4HepMCInterface.cpp:209
G4HepMCInterface::HPlist
std::vector< G4HepMCParticle * > HPlist
Definition G4HepMCInterface.h:82
G4HepMCInterface::GeneratePrimaryVertex
virtual void GeneratePrimaryVertex(G4Event *anEvent)
Definition G4HepMCInterface.cpp:584
G4HepMCInterface::GenerateHepMCEvent
virtual HepMC::GenEvent * GenerateHepMCEvent()
Definition G4HepMCInterface.cpp:576
G4HepMCInterface::CheckType
G4int CheckType(const HepMC::GenEvent *hepmcevt)
Definition G4HepMCInterface.cpp:182
G4Svc
Definition G4Svc.h:31
G4Svc::GetBeamSizeZ
double GetBeamSizeZ()
Definition G4Svc.h:79
G4Svc::GetBeamPosX
double GetBeamPosX()
Definition G4Svc.h:73
G4Svc::GetBeamPosZ
double GetBeamPosZ()
Definition G4Svc.h:75
G4Svc::GetBeamShiftPz
double GetBeamShiftPz()
Definition G4Svc.h:83
G4Svc::GetBoostLab
bool GetBoostLab()
Definition G4Svc.h:95
G4Svc::GetNBunch
double GetNBunch()
Definition G4Svc.h:87
G4Svc::GetBeamShiftPx
double GetBeamShiftPx()
Definition G4Svc.h:81
G4Svc::GetBeamDeltaTime
double GetBeamDeltaTime()
Definition G4Svc.h:86
G4Svc::GetBunchTimeSigma
double GetBunchTimeSigma()
Definition G4Svc.h:88
G4Svc::GetBeamShiftPy
double GetBeamShiftPy()
Definition G4Svc.h:82
G4Svc::GetBeamStartTime
double GetBeamStartTime()
Definition G4Svc.h:85
G4Svc::GetSetBeamShift
bool GetSetBeamShift()
Definition G4Svc.h:96
G4Svc::SetBeamTime
void SetBeamTime(double value)
Definition G4Svc.h:91
G4Svc::GetBeamSizeX
double GetBeamSizeX()
Definition G4Svc.h:77
G4Svc::GetBeamPosY
double GetBeamPosY()
Definition G4Svc.h:74
G4Svc::GetBeamSizeY
double GetBeamSizeY()
Definition G4Svc.h:78
ns
#define ns(x)
Definition xmltok.c:1355