10 #include <Geant4/G4Event.hh>
11 #include <Geant4/G4IonTable.hh>
12 #include <Geant4/G4ParticleDefinition.hh>
13 #include <Geant4/G4ParticleTable.hh>
14 #include <Geant4/G4PrimaryParticle.hh>
15 #include <Geant4/G4PrimaryVertex.hh>
16 #include <Geant4/G4String.hh>
17 #include <Geant4/G4SystemOfUnits.hh>
18 #include <Geant4/G4ThreeVector.hh>
19 #include <Geant4/G4Types.hh>
36 map<int, PHG4VtxPoint*>::const_iterator vtxiter;
37 multimap<int, PHG4Particle*>::const_iterator particle_iter;
38 std::pair<std::map<int, PHG4VtxPoint*>::const_iterator, std::map<int, PHG4VtxPoint*>::const_iterator> vtxbegin_end =
inEvent->GetVertices();
40 for (vtxiter = vtxbegin_end.first; vtxiter != vtxbegin_end.second; ++vtxiter)
45 G4ThreeVector
position((*vtxiter->second).get_x() *
cm, (*vtxiter->second).get_y() *
cm, (*vtxiter->second).get_z() *
cm);
46 G4PrimaryVertex*
vertex =
new G4PrimaryVertex(
position, (*vtxiter->second).get_t() * nanosecond);
47 pair<multimap<int, PHG4Particle*>::const_iterator, multimap<int, PHG4Particle*>::const_iterator> particlebegin_end =
inEvent->GetParticles(vtxiter->first);
48 for (particle_iter = particlebegin_end.first; particle_iter != particlebegin_end.second; ++particle_iter)
61 if (!(*particle_iter->second).get_pid())
63 G4String particleName = (*particle_iter->second).get_name();
64 G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
65 G4ParticleDefinition* particledef = particleTable->FindParticle(particleName);
68 (*particle_iter->second).set_pid(particledef->GetPDGEncoding());
72 cout <<
PHWHERE <<
"Cannot get PDG value for particle " << particleName
73 <<
", dropping it" << endl;
77 G4PrimaryParticle* g4part =
nullptr;
78 if (!(*particle_iter->second).get_pid())
80 G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
81 G4ParticleDefinition* particle_definition = particleTable->FindParticle((*particle_iter->second).get_name());
82 if (particle_definition)
84 G4double
mass = particle_definition->GetPDGMass();
85 g4part =
new G4PrimaryParticle(particle_definition);
86 double ekin = sqrt((*particle_iter->second).get_px() * (*particle_iter->second).get_px() +
87 (*particle_iter->second).get_py() * (*particle_iter->second).get_py() +
88 (*particle_iter->second).get_pz() * (*particle_iter->second).get_pz());
91 g4part->SetKineticEnergy(ekin *
GeV);
92 g4part->SetMass(mass);
93 G4ThreeVector
v((*particle_iter->second).get_px(), (*particle_iter->second).get_py(), (*particle_iter->second).get_pz());
94 G4ThreeVector vunit =
v.unit();
95 g4part->SetMomentumDirection(vunit);
96 g4part->SetCharge(particle_definition->GetPDGCharge());
97 G4ThreeVector particle_polarization;
98 g4part->SetPolarization(particle_polarization.x(),
99 particle_polarization.y(),
100 particle_polarization.z());
104 cout <<
PHWHERE <<
" cannot get G4 particle definition" << endl;
105 cout <<
"you should have never gotten here, please check this in detail" << endl;
106 cout <<
"exiting now" << endl;
113 if ((*particle_iter->second).isIon())
115 G4ParticleDefinition* ion = G4IonTable::GetIonTable()->GetIon((*particle_iter->second).get_Z(), (*particle_iter->second).get_A(), (*particle_iter->second).get_ExcitEnergy() *
GeV);
116 g4part =
new G4PrimaryParticle(ion);
117 g4part->SetCharge((*particle_iter->second).get_IonCharge());
118 g4part->SetMomentum((*particle_iter->second).get_px() *
GeV,
119 (*particle_iter->second).get_py() *
GeV,
120 (*particle_iter->second).get_pz() *
GeV);
122 else if ((*particle_iter->second).get_pid() > 1000000000)
124 G4ParticleDefinition* ion = G4IonTable::GetIonTable()->GetIon((*particle_iter->second).get_pid());
127 g4part =
new G4PrimaryParticle(ion);
130 g4part->SetCharge(ion->GetPDGCharge());
131 g4part->SetMomentum((*particle_iter->second).get_px() *
GeV,
132 (*particle_iter->second).get_py() *
GeV,
133 (*particle_iter->second).get_pz() *
GeV);
137 cout << __PRETTY_FUNCTION__ <<
": WARNING : PDG ID of " << (*particle_iter->second).get_pid() <<
" is not a valid ion! Therefore, this particle is ignored in processing :";
138 (*particle_iter->second).
identify();
143 g4part =
new G4PrimaryParticle((*particle_iter->second).get_pid(),
144 (*particle_iter->second).get_px() *
GeV,
145 (*particle_iter->second).get_py() *
GeV,
146 (*particle_iter->second).get_pz() *
GeV);
158 g4part->SetUserInformation(userdata);
159 vertex->SetPrimary(g4part);
163 anEvent->AddPrimaryVertex(vertex);