ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4PrimaryGeneratorAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4PrimaryGeneratorAction.cc
2 
3 #include "PHG4InEvent.h"
4 #include "PHG4Particle.h"
6 #include "PHG4VtxPoint.h"
7 
8 #include <phool/phool.h>
9 
10 #include <Geant4/G4Event.hh>
11 #include <Geant4/G4IonTable.hh>
12 #include <Geant4/G4ParticleDefinition.hh> // for G4ParticleDefinition
13 #include <Geant4/G4ParticleTable.hh>
14 #include <Geant4/G4PrimaryParticle.hh>
15 #include <Geant4/G4PrimaryVertex.hh>
16 #include <Geant4/G4String.hh> // for G4String
17 #include <Geant4/G4SystemOfUnits.hh>
18 #include <Geant4/G4ThreeVector.hh>
19 #include <Geant4/G4Types.hh> // for G4double
20 
21 #include <cmath> // for sqrt
22 #include <cstdlib>
23 #include <iostream>
24 #include <map>
25 #include <string> // for operator<<
26 #include <utility> // for pair
27 
28 using namespace std;
29 
31 {
32  if (!inEvent)
33  {
34  return;
35  }
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();
39 
40  for (vtxiter = vtxbegin_end.first; vtxiter != vtxbegin_end.second; ++vtxiter)
41  {
42  // cout << "vtx number: " << vtxiter->first << endl;
43  // (*vtxiter->second).identify();
44  // expected units are cm !
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)
49  {
50  // cout << "PHG4PrimaryGeneratorAction: dealing with" << endl;
51  // (particle_iter->second)->identify();
52 
53  // this is really ugly, and maybe it can be streamlined. Initially it was clear cut, if we only give a particle by its name,
54  // we find it here in the G4 particle table, find the
55  // PDG id and then hand it off with the momentum to G4PrimaryParticle
56  // We also have the capability to give a particle a PDG id and then we don't need this translation (the pdg/particle name lookup is
57  // done somewhere else, maybe this should be rethought)
58  // The problem is that geantinos have the pdg pid = 0 but handing this off to the G4PrimaryParticle ctor will just drop it. So
59  // after going through this pdg id lookup once, we have to go through it again in case it is still zero and treat the
60  // geantinos specially. Probably this can be combined with some thought, but rigth now I don't have time for this
61  if (!(*particle_iter->second).get_pid())
62  {
63  G4String particleName = (*particle_iter->second).get_name();
65  G4ParticleDefinition* particledef = particleTable->FindParticle(particleName);
66  if (particledef)
67  {
68  (*particle_iter->second).set_pid(particledef->GetPDGEncoding());
69  }
70  else
71  {
72  cout << PHWHERE << "Cannot get PDG value for particle " << particleName
73  << ", dropping it" << endl;
74  continue;
75  }
76  }
77  G4PrimaryParticle* g4part = nullptr;
78  if (!(*particle_iter->second).get_pid()) // deal with geantinos which have pid=0
79  {
81  G4ParticleDefinition* particle_definition = particleTable->FindParticle((*particle_iter->second).get_name());
82  if (particle_definition)
83  {
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());
89 
90  // expected momentum unit is GeV
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());
101  }
102  else
103  {
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;
107  exit(1);
108  }
109  }
110  else
111  {
112  // expected momentum unit is GeV
113  if ((*particle_iter->second).isIon())
114  {
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);
121  }
122  else if ((*particle_iter->second).get_pid() > 1000000000) // PDG encoding for ion, even without explicit ion tag in PHG4Particle
123  {
124  G4ParticleDefinition* ion = G4IonTable::GetIonTable()->GetIon((*particle_iter->second).get_pid());
125  if (ion)
126  {
127  g4part = new G4PrimaryParticle(ion);
128  // explicit set the ion to be fully ionized.
129  // if partically ionized atom is used in the future, here is the entry point to update it.
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);
134  }
135  else
136  {
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();
139  }
140  }
141  else
142  {
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);
147  }
148  }
149 
150  //if (inEvent->isEmbeded(particle_iter->second))
151  // Do this for all primaries, not just the embedded particle, so that
152  // we can carry the barcode information forward.
153 
154  if (g4part)
155  {
156  PHG4UserPrimaryParticleInformation* userdata = new PHG4UserPrimaryParticleInformation(inEvent->isEmbeded(particle_iter->second));
157  userdata->set_user_barcode((*particle_iter->second).get_barcode());
158  g4part->SetUserInformation(userdata);
159  vertex->SetPrimary(g4part);
160  }
161  }
162  // vertex->Print();
163  anEvent->AddPrimaryVertex(vertex);
164  }
165  return;
166 }