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