294{
295
298
299
300 ISvcLocator* svcLocator = Gaudi::svcLocator();
303 StatusCode sc=svcLocator->service("G4Svc", tmpSvc);
304 m_G4Svc=
dynamic_cast<G4Svc *
>(tmpSvc);
305
306
307
310
311
314
315
318 {
319 for(HepMC::GenEvent::vertex_const_iterator vitr= hepmcevt->vertices_begin();
320 vitr != hepmcevt->vertices_end(); ++vitr )
321 {
322 G4int vtxFlag=0;
323 G4int pOut = (*vitr)->particles_out_size();
324 HepMC::GenVertex::particle_iterator pitr = (*vitr)->particles_begin(HepMC::children);
325 G4int pdgcode= (*pitr)-> pdg_id();
326 if(pOut == 1 &&
abs(pdgcode) == 11)
327 vtxFlag=1;
328 G4int barcodeVtx = (*vitr)->barcode();
329
330 for (HepMC::GenVertex::particle_iterator pitr= (*vitr)->particles_begin(HepMC::children);
331 pitr != (*vitr)->particles_end(HepMC::children); ++pitr)
332 {
333 G4int pdgcode = (*pitr)->pdg_id();
334 if(vtxFlag==0)
335 {
336 if(pdgcode==-22) pdgcode=22;
337 G4LorentzVector p;
342 G4PrimaryParticle* particle = 0;
343
344
345
346
347 if(pdgcode==m_geantinoPdgId){
348 G4ParticleDefinition* pd = G4ParticleTable::GetParticleTable()->FindParticle("geantino");
349 particle = new G4PrimaryParticle(pd, p.x()*GeV, p.y()*GeV, p.z()*GeV);
350 (*pitr)->set_pdg_id(particle->GetPDGcode());
351 }
352 else if(pdgcode==m_chargedgeantinoPdgId){
353 G4ParticleDefinition* pd = G4ParticleTable::GetParticleTable()->FindParticle("chargedgeantino");
354 particle = new G4PrimaryParticle(pd, p.x()*GeV, p.y()*GeV, p.z()*GeV);
355 (*pitr)->set_pdg_id(particle->GetPDGcode());
356 }
357
358
359
360 else{
361 particle = new G4PrimaryParticle(pdgcode, p.x()*GeV, p.y()*GeV, p.z()*GeV);
362 }
363 if(!particle){
364 G4cout << "ERROR! Cannot create particle PDGCode=" << pdgcode << G4endl;
365 }
366
367
368 G4int ISTHEP = (*pitr)->status();
369
370 G4int barcodeEndVtx = 0;
371 if((*pitr)->end_vertex())
372 barcodeEndVtx = (*pitr)->end_vertex()->barcode();
373
375 HPlist.push_back(hepmcParticle);
376
377 for(
size_t ii=0; ii<
HPlist.size(); ii++ )
378 {
379 if(barcodeVtx ==
HPlist[ii]->GetBarcodeEndVtx())
380 {
381 HPlist[ii]->GetTheParticle()->SetDaughter(particle);
382 hepmcParticle->
Done();
383 break;
384 }
385 }
386 }
387 }
388 }
389 }
390 else
391 {
392 for(HepMC::GenEvent::vertex_const_iterator vitr= hepmcevt->vertices_begin();
393 vitr != hepmcevt->vertices_end(); ++vitr )
394 {
395 for (HepMC::GenVertex::particle_iterator pitr= (*vitr)->particles_begin(HepMC::children);
396 pitr != (*vitr)->particles_end(HepMC::children); ++pitr)
397 {
398 G4int ISTHEP = (*pitr)->status();
399 G4LorentzVector p;
404
405 G4int pdgcode = (*pitr)->pdg_id();
406 G4int barcodeEndVtx = 0;
407 if((*pitr)->end_vertex())
408 barcodeEndVtx = (*pitr)->end_vertex()->barcode();
409 G4PrimaryParticle* particle = 0;
410
411
412
413
414 if(pdgcode==m_geantinoPdgId){
415 G4ParticleDefinition* pd = G4ParticleTable::GetParticleTable()->FindParticle("geantino");
416 particle = new G4PrimaryParticle(pd, p.x()*GeV, p.y()*GeV, p.z()*GeV);
417 (*pitr)->set_pdg_id(particle->GetPDGcode());
418 }
419 else if(pdgcode==m_chargedgeantinoPdgId){
420 G4ParticleDefinition* pd = G4ParticleTable::GetParticleTable()->FindParticle("chargedgeantino");
421 particle = new G4PrimaryParticle(pd, p.x()*GeV, p.y()*GeV, p.z()*GeV);
422 (*pitr)->set_pdg_id(particle->GetPDGcode());
423 }
424
425
426
427 else{
428 particle = new G4PrimaryParticle(pdgcode, p.x()*GeV, p.y()*GeV, p.z()*GeV);
429 }
430 if(!particle){
431 G4cout << "ERROR! Cannot create particle PDGCode=" << pdgcode << G4endl;
432 }
433
434
435
437 HPlist.push_back(hepmcParticle);
438 if(ISTHEP>1)
439 hepmcParticle->
Done();
440 }
441 }
442 }
443
444
447
448
449 G4double pmPosX,pmPosY,pmPosZ,pmTime;
450 G4double tmpPosX,tmpPosY,tmpPosZ,tmpT;
451 G4double beamPosX,beamPosY,beamPosZ,beamSizeX,beamSizeY,beamSizeZ;
452
453 if(m_RealizationSvc->
UseDBFlag() ==
false) {
457
461 }
466
470 }
475
479 }
480
481 G4double gaussX = G4RandGauss::shoot();
482 G4double gaussY = G4RandGauss::shoot();
483 G4double gaussZ = G4RandGauss::shoot();
484 G4double gaussT = G4RandGauss::shoot();
485
490
491 G4double ran=G4UniformRand();
492 G4double beamTime = bunchTimeSigma*G4RandGauss::shoot() + beamStartTime + beamDeltaTime*int(ran*nBunch);
493
494 tmpPosX = (beamPosX + beamSizeX*gaussX ) *mm;
495 tmpPosY = (beamPosY + beamSizeY*gaussY ) *mm;
496 tmpPosZ = (beamPosZ + beamSizeZ*gaussZ ) *mm;
497 tmpT = (beamSizeZ * gaussT ) * mm/c_light +beamTime;
498
499 G4LorentzVector tmpv(tmpPosX,tmpPosY,tmpPosZ,tmpT*c_light/mm);
500
501 HepMC::GenEvent::vertex_const_iterator vitr0= hepmcevt->vertices_begin();
502 G4LorentzVector xvtx0 ;
503 xvtx0.setX( (*vitr0)-> position().
x());
504 xvtx0.setY( (*vitr0)-> position().
y());
505 xvtx0.setZ( (*vitr0)-> position().z());
506 xvtx0.setT( (*vitr0)-> position().
t());
507 pmPosX = xvtx0.x()*mm + tmpPosX;
508 pmPosY = xvtx0.y()*mm + tmpPosY;
509 pmPosZ = xvtx0.z()*mm + tmpPosZ;
510 pmTime = xvtx0.t()*mm/c_light + tmpT;
511
513 {
514 G4cout<<G4endl;
515 G4cout<<xvtx0.x()<<" "<<xvtx0.y()<<" "<<xvtx0.z()<<" "
516 <<beamSizeX*gaussX<<" "
517 <<beamSizeY*gaussY<<" "
518 <<beamSizeZ*gaussZ<<" "<<G4endl;
519 G4cout<<xvtx0.t()* mm/c_light<<" "
520 <<beamSizeZ * gaussT/sqrt(2)*mm/c_light<<" "
521 <<beamTime<<G4endl;
522 }
523 for(HepMC::GenEvent::vertex_const_iterator vitr= hepmcevt->vertices_begin();
524 vitr != hepmcevt->vertices_end(); ++vitr )
525 {
526 G4LorentzVector xvtx;
527 xvtx.setX((*vitr)-> position().
x());
528 xvtx.setY((*vitr)-> position().
y());
529 xvtx.setZ((*vitr)-> position().z());
530 xvtx.setT((*vitr)-> position().
t());
531 (*vitr)->set_position(xvtx+tmpv);
532 }
534
535 G4PrimaryVertex* g4vtx= new G4PrimaryVertex(pmPosX,pmPosY,pmPosZ,pmTime);
536
537
538 for(
size_t ii=0; ii<
HPlist.size(); ii++ )
539 {
540 if(
HPlist[ii]->GetISTHEP() > 0 )
541
542 {
543 G4PrimaryParticle* initialParticle =
HPlist[ii]->GetTheParticle();
544 g4vtx->SetPrimary( initialParticle );
545 }
546 }
547
548
549 for(
size_t iii=0;iii<
HPlist.size();iii++)
552
553 g4event->AddPrimaryVertex(g4vtx);
554}
void Print(const HepMC::GenEvent *hepmcevt)
void Boost(HepMC::GenEvent *hepmcevt)
std::vector< G4HepMCParticle * > HPlist
G4int CheckType(const HepMC::GenEvent *hepmcevt)
double GetBeamDeltaTime()
double GetBunchTimeSigma()
double GetBeamStartTime()
void SetBeamTime(double value)