g4simulation/g4main/PHG4PrimaryGeneratorAction.cc

0001 #include "PHG4PrimaryGeneratorAction.h"
0002
0003 #include "PHG4InEvent.h"
0004 #include "PHG4Particle.h"
0005 #include "PHG4UserPrimaryParticleInformation.h"
0006 #include "PHG4VtxPoint.h"
0007
0008 #include <phool/phool.h>
0009
0010 #include <Geant4/G4Event.hh>
0011 #include <Geant4/G4IonTable.hh>
0012 #include <Geant4/G4ParticleDefinition.hh>  // for G4ParticleDefinition
0013 #include <Geant4/G4ParticleTable.hh>
0014 #include <Geant4/G4PrimaryParticle.hh>
0015 #include <Geant4/G4PrimaryVertex.hh>
0016 #include <Geant4/G4String.hh>  // for G4String
0017 #include <Geant4/G4SystemOfUnits.hh>
0018 #include <Geant4/G4ThreeVector.hh>
0019 #include <Geant4/G4Types.hh>  // for G4double
0020
0021 #include <cmath>  // for sqrt
0022 #include <cstdlib>
0023 #include <iostream>
0024 #include <map>
0025 #include <string>   // for operator<<
0026 #include <utility>  // for pair
0027
0028 using namespace std;
0029
0030 void PHG4PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
0031 {
0032   if (!inEvent)
0033   {
0034     return;
0035   }
0036   map<int, PHG4VtxPoint*>::const_iterator vtxiter;
0037   multimap<int, PHG4Particle*>::const_iterator particle_iter;
0038   std::pair<std::map<int, PHG4VtxPoint*>::const_iterator, std::map<int, PHG4VtxPoint*>::const_iterator> vtxbegin_end = inEvent->GetVertices();
0039
0040   for (vtxiter = vtxbegin_end.first; vtxiter != vtxbegin_end.second; ++vtxiter)
0041   {
0042     //       cout << "vtx number: " << vtxiter->first << endl;
0043     //       (*vtxiter->second).identify();
0044     // expected units are cm !
0045     G4ThreeVector position((*vtxiter->second).get_x() * cm, (*vtxiter->second).get_y() * cm, (*vtxiter->second).get_z() * cm);
0046     G4PrimaryVertex* vertex = new G4PrimaryVertex(position, (*vtxiter->second).get_t() * nanosecond);
0047     pair<multimap<int, PHG4Particle*>::const_iterator, multimap<int, PHG4Particle*>::const_iterator> particlebegin_end = inEvent->GetParticles(vtxiter->first);
0048     for (particle_iter = particlebegin_end.first; particle_iter != particlebegin_end.second; ++particle_iter)
0049     {
0050       // cout << "PHG4PrimaryGeneratorAction: dealing with" << endl;
0051       //  (particle_iter->second)->identify();
0052
0053       // this is really ugly, and maybe it can be streamlined. Initially it was clear cut, if we only give a particle by its name,
0054       // we find it here in the G4 particle table, find the
0055       // PDG id and then hand it off with the momentum to G4PrimaryParticle
0056       // 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
0057       // done somewhere else, maybe this should be rethought)
0058       // The problem is that geantinos have the pdg pid = 0 but handing this off to the G4PrimaryParticle ctor will just drop it. So
0059       // after going through this pdg id lookup once, we have to go through it again in case it is still zero and treat the
0060       // geantinos specially. Probably this can be combined with some thought, but rigth now I don't have time for this
0061       if (!(*particle_iter->second).get_pid())
0062       {
0063         G4String particleName = (*particle_iter->second).get_name();
0064         G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
0065         G4ParticleDefinition* particledef = particleTable->FindParticle(particleName);
0066         if (particledef)
0067         {
0068           (*particle_iter->second).set_pid(particledef->GetPDGEncoding());
0069         }
0070         else
0071         {
0072           cout << PHWHERE << "Cannot get PDG value for particle " << particleName
0073                << ", dropping it" << endl;
0074           continue;
0075         }
0076       }
0077       G4PrimaryParticle* g4part = nullptr;
0078       if (!(*particle_iter->second).get_pid())  // deal with geantinos which have pid=0
0079       {
0080         G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
0081         G4ParticleDefinition* particle_definition = particleTable->FindParticle((*particle_iter->second).get_name());
0082         if (particle_definition)
0083         {
0084           G4double mass = particle_definition->GetPDGMass();
0085           g4part = new G4PrimaryParticle(particle_definition);
0086           double ekin = sqrt((*particle_iter->second).get_px() * (*particle_iter->second).get_px() +
0087                              (*particle_iter->second).get_py() * (*particle_iter->second).get_py() +
0088                              (*particle_iter->second).get_pz() * (*particle_iter->second).get_pz());
0089
0090           // expected momentum unit is GeV
0091           g4part->SetKineticEnergy(ekin * GeV);
0092           g4part->SetMass(mass);
0093           G4ThreeVector v((*particle_iter->second).get_px(), (*particle_iter->second).get_py(), (*particle_iter->second).get_pz());
0094           G4ThreeVector vunit = v.unit();
0095           g4part->SetMomentumDirection(vunit);
0096           g4part->SetCharge(particle_definition->GetPDGCharge());
0097           G4ThreeVector particle_polarization;
0098           g4part->SetPolarization(particle_polarization.x(),
0099                                   particle_polarization.y(),
0100                                   particle_polarization.z());
0101         }
0102         else
0103         {
0104           cout << PHWHERE << " cannot get G4 particle definition" << endl;
0105           cout << "you should have never gotten here, please check this in detail" << endl;
0106           cout << "exiting now" << endl;
0107           exit(1);
0108         }
0109       }
0110       else
0111       {
0112         // expected momentum unit is GeV
0113         if ((*particle_iter->second).isIon())
0114         {
0115           G4ParticleDefinition* ion = G4IonTable::GetIonTable()->GetIon((*particle_iter->second).get_Z(), (*particle_iter->second).get_A(), (*particle_iter->second).get_ExcitEnergy() * GeV);
0116           g4part = new G4PrimaryParticle(ion);
0117           g4part->SetCharge((*particle_iter->second).get_IonCharge());
0118           g4part->SetMomentum((*particle_iter->second).get_px() * GeV,
0119                               (*particle_iter->second).get_py() * GeV,
0120                               (*particle_iter->second).get_pz() * GeV);
0121         }
0122         else if ((*particle_iter->second).get_pid() > 1000000000)  // PDG encoding for ion, even without explicit ion tag in PHG4Particle
0123         {
0124           G4ParticleDefinition* ion = G4IonTable::GetIonTable()->GetIon((*particle_iter->second).get_pid());
0125           if (ion)
0126           {
0127             g4part = new G4PrimaryParticle(ion);
0128             // explicit set the ion to be fully ionized.
0129             // if partically ionized atom is used in the future, here is the entry point to update it.
0130             g4part->SetCharge(ion->GetPDGCharge());
0131             g4part->SetMomentum((*particle_iter->second).get_px() * GeV,
0132                                 (*particle_iter->second).get_py() * GeV,
0133                                 (*particle_iter->second).get_pz() * GeV);
0134           }
0135           else
0136           {
0137             cout << __PRETTY_FUNCTION__ << ": WARNING : PDG ID of " << (*particle_iter->second).get_pid() << " is not a valid ion! Therefore, this particle is ignored in processing :";
0138             (*particle_iter->second).identify();
0139           }
0140         }
0141         else
0142         {
0143           g4part = new G4PrimaryParticle((*particle_iter->second).get_pid(),
0144                                          (*particle_iter->second).get_px() * GeV,
0145                                          (*particle_iter->second).get_py() * GeV,
0146                                          (*particle_iter->second).get_pz() * GeV);
0147         }
0148       }
0149
0150       // if (inEvent->isEmbeded(particle_iter->second))
0151       //  Do this for all primaries, not just the embedded particle, so that
0152       //  we can carry the barcode information forward.
0153
0154       if (g4part)
0155       {
0156         PHG4UserPrimaryParticleInformation* userdata = new PHG4UserPrimaryParticleInformation(inEvent->isEmbeded(particle_iter->second));
0157         userdata->set_user_barcode((*particle_iter->second).get_barcode());
0158         g4part->SetUserInformation(userdata);
0159         vertex->SetPrimary(g4part);
0160       }
0161     }
0162     //      vertex->Print();
0163     anEvent->AddPrimaryVertex(vertex);
0164   }
0165   return;
0166 }