Skip to content
Snippets Groups Projects
Commit 8ce39f9f authored by lintao@ihep.ac.cn's avatar lintao@ihep.ac.cn
Browse files

Support geantino/chargedgeantino in simulation framework.

parent fb8ee579
No related branches found
No related tags found
No related merge requests found
...@@ -86,10 +86,22 @@ GtGunTool::mutate(MyHepMC::GenEvent& event) { ...@@ -86,10 +86,22 @@ GtGunTool::mutate(MyHepMC::GenEvent& event) {
double charge = 0; double charge = 0;
TParticlePDG* particle = db_pdg->GetParticle(particle_name.c_str()); TParticlePDG* particle = db_pdg->GetParticle(particle_name.c_str());
info() << "According to the " << particle_name << " get " << particle << endmsg;
if (particle) { if (particle) {
pdgcode = particle->PdgCode(); pdgcode = particle->PdgCode();
mass = particle->Mass(); // GeV mass = particle->Mass(); // GeV
charge = particle->Charge()/3; // in ROOT, it is in units of |e|/3 charge = particle->Charge()/3; // in ROOT, it is in units of |e|/3
} else if (particle_name == "geantino") {
info() << "Create geantino ... " << endmsg;
pdgcode = 0;
mass = 0;
charge = 0;
} else if (particle_name == "chargedgeantino") {
info() << "Create chargedgeantino ... " << endmsg;
pdgcode = 0;
mass = 0;
charge = 1; // according to the definition in G4ChargedGeantino.cc
} else { } else {
// guess it is pdg code // guess it is pdg code
pdgcode = atol(particle_name.c_str()); pdgcode = atol(particle_name.c_str());
...@@ -99,7 +111,18 @@ GtGunTool::mutate(MyHepMC::GenEvent& event) { ...@@ -99,7 +111,18 @@ GtGunTool::mutate(MyHepMC::GenEvent& event) {
} }
} }
double energy = m_energymins.value()[i]==m_energymaxs.value()[i] ? m_energymins.value()[i] : CLHEP::RandFlat::shoot(m_energymins.value()[i], m_energymaxs.value()[i]); if (i>=m_energymins.value().size()) {
error() << "Mismatch EnergyMins" << endmsg;
return false;
}
if (i>=m_energymaxs.value().size()) {
error() << "Mismatch EnergyMaxs" << endmsg;
return false;
}
double energy_min = m_energymins.value()[i];
double energy_max = m_energymaxs.value()[i];
double energy = energy_min==energy_max ? energy_min : CLHEP::RandFlat::shoot(energy_min, energy_max);
// create the MC particle // create the MC particle
edm4hep::MCParticle mcp = event.m_mc_vec.create(); edm4hep::MCParticle mcp = event.m_mc_vec.create();
...@@ -126,6 +149,26 @@ GtGunTool::mutate(MyHepMC::GenEvent& event) { ...@@ -126,6 +149,26 @@ GtGunTool::mutate(MyHepMC::GenEvent& event) {
// direction // direction
// by default, randomize the direction // by default, randomize the direction
if (i>=m_thetamins.value().size()) {
error() << "Mismatch ThetaMins" << endmsg;
return false;
}
if (i>=m_thetamaxs.value().size()) {
error() << "Mismatch ThetaMaxs" << endmsg;
return false;
}
if (i>=m_phimins.value().size()) {
error() << "Mismatch PhiMins" << endmsg;
return false;
}
if (i>=m_phimaxs.value().size()) {
error() << "Mismatch PhiMaxs" << endmsg;
return false;
}
double theta = m_thetamins.value()[i]==m_thetamaxs.value()[i] ? m_thetamins.value()[i] : CLHEP::RandFlat::shoot(m_thetamins.value()[i], m_thetamaxs.value()[i]); double theta = m_thetamins.value()[i]==m_thetamaxs.value()[i] ? m_thetamins.value()[i] : CLHEP::RandFlat::shoot(m_thetamins.value()[i], m_thetamaxs.value()[i]);
double phi = m_phimins .value()[i]==m_phimaxs .value()[i] ? m_phimins .value()[i] : CLHEP::RandFlat::shoot(m_phimins .value()[i], m_phimaxs .value()[i]); double phi = m_phimins .value()[i]==m_phimaxs .value()[i] ? m_phimins .value()[i] : CLHEP::RandFlat::shoot(m_phimins .value()[i], m_phimaxs .value()[i]);
double costheta = cos(theta*acos(-1)/180); double costheta = cos(theta*acos(-1)/180);
......
...@@ -40,8 +40,18 @@ bool G4PrimaryCnvTool::mutate(G4Event* anEvent) { ...@@ -40,8 +40,18 @@ bool G4PrimaryCnvTool::mutate(G4Event* anEvent) {
// pdg/particle // pdg/particle
int pdgcode = p.getPDG(); int pdgcode = p.getPDG();
G4ParticleTable* particletbl = G4ParticleTable::GetParticleTable(); G4ParticleTable* particletbl = G4ParticleTable::GetParticleTable();
G4ParticleDefinition* particle_def = particletbl->FindParticle(pdgcode); G4ParticleDefinition* particle_def = nullptr;
// handle the several exceptions
if (pdgcode == 0 && p.getCharge() == 0) {
// this is geantino
particle_def = particletbl->FindParticle("geantino");
} else if (pdgcode == 0 && p.getCharge() == 1) {
// this is chargedgeantino
particle_def = particletbl->FindParticle("chargedgeantino");
} else {
particle_def = particletbl->FindParticle(pdgcode);
}
// momentum // momentum
const edm4hep::Vector3f& momentum = p.getMomentum(); const edm4hep::Vector3f& momentum = p.getMomentum();
G4PrimaryParticle* g4prim = new G4PrimaryParticle(particle_def, G4PrimaryParticle* g4prim = new G4PrimaryParticle(particle_def,
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment