File indexing completed on 2025-02-04 02:43:34
0001 #include "SimG4Core/Generators/interface/Generator3.h"
0002 #include "SimG4Core/Generators/interface/HepMCParticle.h"
0003 #include "SimG4Core/Generators/interface/LumiMonitorFilter.h"
0004
0005 #include "SimG4Core/Notification/interface/SimG4Exception.h"
0006
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008
0009 #include "HepPDT/ParticleID.hh"
0010
0011 #include "G4Event.hh"
0012 #include "G4HEPEvtParticle.hh"
0013 #include "G4Log.hh"
0014 #include "G4ParticleDefinition.hh"
0015 #include "G4PhysicalConstants.hh"
0016 #include <CLHEP/Units/SystemOfUnits.h>
0017 #include "G4UnitsTable.hh"
0018
0019 #include "HepMC3/Print.h"
0020
0021 #include <sstream>
0022
0023 using namespace edm;
0024
0025 Generator3::Generator3(const ParameterSet &p)
0026 : fPCuts(p.getParameter<bool>("ApplyPCuts")),
0027 fPtransCut(p.getParameter<bool>("ApplyPtransCut")),
0028 fEtaCuts(p.getParameter<bool>("ApplyEtaCuts")),
0029 fPhiCuts(p.getParameter<bool>("ApplyPhiCuts")),
0030 fSlepton(p.getParameter<bool>("IsSlepton")),
0031 theMinPhiCut(p.getParameter<double>("MinPhiCut")),
0032 theMaxPhiCut(p.getParameter<double>("MaxPhiCut")),
0033 theMinEtaCut(p.getParameter<double>("MinEtaCut")),
0034 theMaxEtaCut(p.getParameter<double>("MaxEtaCut")),
0035 theMinPCut(p.getParameter<double>("MinPCut")),
0036 theMaxPCut(p.getParameter<double>("MaxPCut")),
0037 theEtaCutForHector(p.getParameter<double>("EtaCutForHector")),
0038 verbose(p.getUntrackedParameter<int>("Verbosity", 0)),
0039 fLumiFilter(nullptr),
0040 evt_(nullptr),
0041 vtx_(nullptr),
0042 weight_(0),
0043 Z_lmin(0),
0044 Z_lmax(0),
0045 Z_hector(0),
0046 pdgFilterSel(false),
0047 fPDGFilter(false) {
0048 bool lumi = p.getParameter<bool>("ApplyLumiMonitorCuts");
0049 if (lumi) {
0050 fLumiFilter = new LumiMonitorFilter();
0051 }
0052
0053 double theRDecLenCut = p.getParameter<double>("RDecLenCut") * CLHEP::cm;
0054 theDecRCut2 = theRDecLenCut * theRDecLenCut;
0055
0056 theMinPtCut2 = theMinPCut * theMinPCut;
0057
0058 double theDecLenCut = p.getParameter<double>("LDecLenCut") * CLHEP::cm;
0059
0060 maxZCentralCMS = p.getParameter<double>("MaxZCentralCMS") * CLHEP::m;
0061
0062 fFiductialCuts = (fPCuts || fPtransCut || fEtaCuts || fPhiCuts);
0063
0064 pdgFilter.resize(0);
0065 if (p.exists("PDGselection")) {
0066 pdgFilterSel = (p.getParameter<edm::ParameterSet>("PDGselection")).getParameter<bool>("PDGfilterSel");
0067 pdgFilter = (p.getParameter<edm::ParameterSet>("PDGselection")).getParameter<std::vector<int>>("PDGfilter");
0068 if (!pdgFilter.empty()) {
0069 fPDGFilter = true;
0070 std::stringstream ss;
0071 ss << "SimG4Core/Generator3: ";
0072 if (pdgFilterSel) {
0073 ss << " Selecting only PDG ID = ";
0074 } else {
0075 ss << " Filtering out PDG ID = ";
0076 }
0077 for (unsigned int ii = 0; ii < pdgFilter.size(); ++ii) {
0078 ss << pdgFilter[ii] << " ";
0079 }
0080 edm::LogVerbatim("SimG4CoreGenerator3") << ss.str();
0081 }
0082 }
0083
0084 if (fEtaCuts) {
0085 Z_lmax = theRDecLenCut * ((1 - exp(-2 * theMaxEtaCut)) / (2 * exp(-theMaxEtaCut)));
0086 Z_lmin = theRDecLenCut * ((1 - exp(-2 * theMinEtaCut)) / (2 * exp(-theMinEtaCut)));
0087 }
0088
0089 Z_hector = theRDecLenCut * ((1 - exp(-2 * theEtaCutForHector)) / (2 * exp(-theEtaCutForHector)));
0090
0091 edm::LogVerbatim("SimG4CoreGenerator3")
0092 << "SimG4Core/Generator3: Rdecaycut= " << theRDecLenCut / CLHEP::cm
0093 << " cm; Zdecaycut= " << theDecLenCut / CLHEP::cm << "Z_min= " << Z_lmin / CLHEP::cm
0094 << " cm; Z_max= " << Z_lmax / CLHEP::cm << " cm;\n"
0095 << " MaxZCentralCMS = " << maxZCentralCMS / CLHEP::m << " m;"
0096 << " Z_hector = " << Z_hector / CLHEP::cm << " cm\n"
0097 << " ApplyCuts: " << fFiductialCuts << " PCuts: " << fPCuts << " PtransCut: " << fPtransCut
0098 << " EtaCut: " << fEtaCuts << " PhiCut: " << fPhiCuts << " LumiMonitorCut: " << lumi;
0099 if (fFiductialCuts) {
0100 edm::LogVerbatim("SimG4CoreGenerator3")
0101 << "SimG4Core/Generator3: Pmin(GeV)= " << theMinPCut << "; Pmax(GeV)= " << theMaxPCut
0102 << "; EtaMin= " << theMinEtaCut << "; EtaMax= " << theMaxEtaCut << "; PhiMin(rad)= " << theMinPhiCut
0103 << "; PhiMax(rad)= " << theMaxPhiCut;
0104 }
0105 if (lumi) {
0106 fLumiFilter->Describe();
0107 }
0108 }
0109
0110 Generator3::~Generator3() { delete fLumiFilter; }
0111
0112 void Generator3::HepMC2G4(const HepMC3::GenEvent *evt_orig, G4Event *g4evt) {
0113 HepMC3::GenEvent *evt = new HepMC3::GenEvent(*evt_orig);
0114
0115
0116 if ((evt->vertices()).empty()) {
0117 std::stringstream ss;
0118 ss << "SimG4Core/Generator3: in event " << g4evt->GetEventID() << " Corrupted Event - GenEvent with no vertex \n";
0119 throw SimG4Exception(ss.str());
0120 }
0121
0122 if (!evt->weights().empty()) {
0123 weight_ = evt->weights()[0];
0124 for (unsigned int iw = 1; iw < evt->weights().size(); ++iw) {
0125
0126 if (evt->weights()[iw] <= 0)
0127 break;
0128 weight_ *= evt->weights()[iw];
0129 }
0130 }
0131
0132 if (vtx_ != nullptr) {
0133 delete vtx_;
0134 }
0135
0136 for (const HepMC3::GenVertexPtr &v : evt->vertices()) {
0137 vtx_ =
0138 new math::XYZTLorentzVector((v->position()).x(), (v->position()).y(), (v->position()).z(), (v->position()).t());
0139 break;
0140 }
0141
0142 edm::LogVerbatim("SimG4CoreGenerator3")
0143 << "Generator3: primary Vertex = (" << vtx_->x() << ", " << vtx_->y() << ", " << vtx_->z() << ")";
0144
0145 unsigned int ng4vtx = 0;
0146 unsigned int ng4par = 0;
0147
0148 for (const HepMC3::GenVertexPtr &vitr : evt->vertices()) {
0149
0150
0151
0152
0153 G4bool qvtx = false;
0154 for (const HepMC3::GenParticlePtr &pitr : vitr->particles_out()) {
0155
0156
0157
0158
0159 int status = pitr->status();
0160 int pdg = pitr->pid();
0161 if (status > 3 && isExotic(pdg) && (!(isExoticNonDetectable(pdg)))) {
0162
0163
0164
0165
0166
0167 status = 2;
0168 }
0169 if (status == 2 && abs(pdg) == 9900015) {
0170 status = 3;
0171 }
0172
0173
0174 if (status == 1) {
0175
0176 if (fPDGFilter && !pdgFilterSel && IsInTheFilterList(pdg)) {
0177 continue;
0178 }
0179
0180 qvtx = true;
0181 if (verbose > 2)
0182 LogDebug("SimG4CoreGenerator3")
0183 << "GenVertex barcode = " << vitr->id() << " " << qvtx
0184 << " selected for GenParticle barcode = " << pitr->id();
0185 break;
0186 }
0187
0188
0189
0190 else if (status == 2) {
0191 if (pitr->end_vertex() != nullptr) {
0192 double xx = (pitr->end_vertex())->position().x();
0193 double yy = (pitr->end_vertex())->position().y();
0194 double r_dd = xx * xx + yy * yy;
0195 if (r_dd > theDecRCut2) {
0196 qvtx = true;
0197 if (verbose > 2)
0198 LogDebug("SimG4CoreGenerator3")
0199 << "GenVertex barcode = " << vitr->id() << " selected for GenParticle barcode = " << pitr->id()
0200 << " radius = " << std::sqrt(r_dd);
0201 break;
0202 }
0203 } else {
0204
0205
0206 qvtx = true;
0207 break;
0208 }
0209 }
0210 }
0211
0212
0213
0214 if (!qvtx) {
0215 continue;
0216 }
0217
0218 double x1 = vitr->position().x() * CLHEP::mm;
0219 double y1 = vitr->position().y() * CLHEP::mm;
0220 double z1 = vitr->position().z() * CLHEP::mm;
0221 double t1 = vitr->position().t() * CLHEP::mm / CLHEP::c_light;
0222
0223 G4PrimaryVertex *g4vtx = new G4PrimaryVertex(x1, y1, z1, t1);
0224
0225 for (const HepMC3::GenParticlePtr &pitr : vitr->particles_out()) {
0226 int status = pitr->status();
0227 int pdg = pitr->pid();
0228 bool hasDecayVertex = (nullptr != pitr->end_vertex());
0229
0230
0231 if (fPDGFilter) {
0232 bool isInTheList = IsInTheFilterList(pdg);
0233 if ((!pdgFilterSel && isInTheList) || (pdgFilterSel && !isInTheList)) {
0234 if (0 < verbose)
0235 edm::LogVerbatim("SimG4CoreGenerator3")
0236 << " Skiped GenParticle barcode= " << pitr->id() << " PDGid= " << pdg << " status= " << status
0237 << " isExotic: " << isExotic(pdg) << " isExoticNotDet: " << isExoticNonDetectable(pdg)
0238 << " isInTheList: " << isInTheList << " hasDecayVertex: " << hasDecayVertex;
0239 continue;
0240 }
0241 }
0242
0243 if (0 < verbose) {
0244 edm::LogVerbatim("SimG4CoreGenerator3")
0245 << "Generator3: pdg= " << pdg << " status= " << status << " hasPreDefinedDecay: " << hasDecayVertex
0246 << "\n isExotic: " << isExotic(pdg) << " isExoticNotDet: " << isExoticNonDetectable(pdg)
0247 << " isInTheList: " << IsInTheFilterList(pdg) << "\n"
0248 << " MaxZCentralCMS = " << maxZCentralCMS / CLHEP::m << " m; (x,y,z,t): (" << x1 << "," << y1 << "," << z1
0249 << "," << t1 << ")";
0250 }
0251 if (status > 3 && isExotic(pdg) && (!(isExoticNonDetectable(pdg)))) {
0252 status = hasDecayVertex ? 2 : 1;
0253 }
0254 if (status == 2 && abs(pdg) == 9900015) {
0255 status = 3;
0256 }
0257
0258
0259 if (2 == status && !hasDecayVertex) {
0260 edm::LogWarning("SimG4CoreGenerator3: in event ")
0261 << g4evt->GetEventID() << " a particle "
0262 << " pdgid= " << pdg << " has status=2 but has no decay vertex, so will be fully tracked by Geant4";
0263 status = 1;
0264 }
0265
0266 double x2 = x1;
0267 double y2 = y1;
0268 double z2 = z1;
0269 double decay_length = 0.0;
0270 if (2 == status) {
0271 x2 = pitr->end_vertex()->position().x();
0272 y2 = pitr->end_vertex()->position().y();
0273 z2 = pitr->end_vertex()->position().z();
0274 decay_length = std::sqrt((x2 - x1) * (x2 - x1) + (y2 - y1) * (y2 - y1) + (z2 - z1) * (z2 - z1));
0275 }
0276
0277 bool toBeAdded = !fFiductialCuts;
0278
0279 double px = pitr->momentum().px();
0280 double py = pitr->momentum().py();
0281 double pz = pitr->momentum().pz();
0282 double ptot = std::sqrt(px * px + py * py + pz * pz);
0283 math::XYZTLorentzVector p(px, py, pz, pitr->momentum().e());
0284
0285 double ximpact = x1;
0286 double yimpact = y1;
0287 double zimpact = z1;
0288
0289
0290
0291 const double minTan = 1.e-20;
0292 if (std::abs(z1) < Z_hector && std::abs(pz) >= minTan * ptot) {
0293 zimpact = (pz > 0.0) ? Z_hector : -Z_hector;
0294 double del = (zimpact - z1) / pz;
0295 ximpact += del * px;
0296 yimpact += del * py;
0297 }
0298 double rimpact2 = ximpact * ximpact + yimpact * yimpact;
0299
0300 if (verbose > 2)
0301 LogDebug("SimG4CoreGenerator3") << "Processing GenParticle barcode= " << pitr->id() << " pdg= " << pdg
0302 << " status= " << pitr->status() << " st= " << status
0303 << " rimpact(cm)= " << std::sqrt(rimpact2) / CLHEP::cm
0304 << " zimpact(cm)= " << zimpact / CLHEP::cm << " ptot(GeV)= " << ptot
0305 << " pz(GeV)= " << pz;
0306
0307
0308
0309 if (1 == status && std::abs(zimpact) >= Z_hector && rimpact2 <= theDecRCut2) {
0310
0311 toBeAdded = (2112 == std::abs(pdg) || 22 == pdg);
0312 if (verbose > 1) {
0313 edm::LogVerbatim("SimG4CoreGenerator3")
0314 << "GenParticle barcode = " << pitr->id() << " very forward; to be added: " << toBeAdded;
0315 }
0316 } else {
0317
0318 if (1 == status && (std::abs(zimpact) < Z_hector || rimpact2 > theDecRCut2)) {
0319
0320 if (fPCuts && (ptot < theMinPCut || ptot > theMaxPCut)) {
0321 continue;
0322 }
0323
0324 if (fPhiCuts) {
0325 double phi = p.phi();
0326 if (phi < theMinPhiCut || phi > theMaxPhiCut) {
0327 continue;
0328 }
0329 }
0330
0331
0332 if (fEtaCuts) {
0333
0334 double xi = x1;
0335 double yi = y1;
0336 double zi = z1;
0337
0338
0339 if (std::abs(pz) >= minTan * ptot) {
0340 if ((zi >= Z_lmax) & (pz < 0.0)) {
0341 zi = Z_lmax;
0342 } else if ((zi <= Z_lmin) & (pz > 0.0)) {
0343 zi = Z_lmin;
0344 } else {
0345 if (pz > 0) {
0346 zi = Z_lmax;
0347 } else {
0348 zi = Z_lmin;
0349 }
0350 }
0351 double del = (zi - z1) / pz;
0352 xi += del * px;
0353 yi += del * py;
0354 }
0355
0356 if ((zi >= Z_lmin) & (zi <= Z_lmax) & (xi * xi + yi * yi < theDecRCut2)) {
0357 continue;
0358 }
0359 }
0360 const HepMC3::GenParticle *ppointer = pitr.get();
0361 if (fLumiFilter && !fLumiFilter->isGoodForLumiMonitor(ppointer)) {
0362 continue;
0363 }
0364 toBeAdded = true;
0365 if (verbose > 1)
0366 edm::LogVerbatim("SimG4CoreGenerator3") << "GenParticle barcode = " << pitr->id() << " passed case 1";
0367
0368
0369
0370
0371 } else if (2 == status && x2 * x2 + y2 * y2 >= theDecRCut2 && (fSlepton || std::abs(z2) < Z_hector)) {
0372 toBeAdded = true;
0373 if (verbose > 1)
0374 edm::LogVerbatim("SimG4CoreGenerator3") << "GenParticle barcode = " << pitr->id() << " passed case 2"
0375 << " decay_length(cm)= " << decay_length / CLHEP::cm;
0376 }
0377 }
0378 if (toBeAdded) {
0379 G4PrimaryParticle *g4prim = new G4PrimaryParticle(pdg, px * CLHEP::GeV, py * CLHEP::GeV, pz * CLHEP::GeV);
0380
0381 if (g4prim->GetG4code() != nullptr) {
0382 g4prim->SetMass(g4prim->GetG4code()->GetPDGMass());
0383 double charge = g4prim->GetG4code()->GetPDGCharge();
0384
0385
0386 if (fPtransCut && 1 == status && 0.0 != charge && px * px + py * py < theMinPtCut2) {
0387 delete g4prim;
0388 continue;
0389 }
0390 g4prim->SetCharge(charge);
0391 }
0392
0393
0394
0395
0396 setGenId(g4prim, pitr->id());
0397
0398 if (2 == status) {
0399 particleAssignDaughters(g4prim, pitr.get(), decay_length);
0400 }
0401 if (verbose > 1)
0402 g4prim->Print();
0403
0404 ++ng4par;
0405 g4vtx->SetPrimary(g4prim);
0406 edm::LogVerbatim("SimG4CoreGenerator3") << " " << ng4par << ". new Geant4 particle pdg= " << pdg
0407 << " Ptot(GeV/c)= " << ptot << " Pt= " << std::sqrt(px * px + py * py)
0408 << " status= " << status << "; dir= " << g4prim->GetMomentumDirection();
0409 }
0410 }
0411
0412 if (verbose > 1)
0413 g4vtx->Print();
0414 g4evt->AddPrimaryVertex(g4vtx);
0415 ++ng4vtx;
0416 }
0417
0418
0419
0420 if (ng4vtx == 0) {
0421 G4PrimaryVertex *g4vtx = new G4PrimaryVertex(0.0, 0.0, 0.0, 0.0);
0422 if (verbose > 1)
0423 g4vtx->Print();
0424
0425 g4evt->AddPrimaryVertex(g4vtx);
0426 }
0427
0428 edm::LogVerbatim("SimG4CoreGenerator3")
0429 << "The list of Geant4 primaries includes " << ng4par << " particles in " << ng4vtx << " vertex";
0430
0431 delete evt;
0432 }
0433
0434 void Generator3::particleAssignDaughters(G4PrimaryParticle *g4p, HepMC3::GenParticle *vp, double decaylength) {
0435 if (verbose > 1) {
0436 LogDebug("SimG4CoreGenerator3") << "Special case of long decay length \n"
0437 << "Assign daughters to mother with decaylength=" << decaylength / CLHEP::cm
0438 << " cm";
0439 }
0440 math::XYZTLorentzVector p(vp->momentum().px(), vp->momentum().py(), vp->momentum().pz(), vp->momentum().e());
0441
0442
0443 double proper_time = decaylength / (p.Beta() * p.Gamma() * c_light);
0444 g4p->SetProperTime(proper_time);
0445
0446 if (verbose > 2) {
0447 LogDebug("SimG4CoreGenerator3") << " px= " << p.px() << " py= " << p.py() << " pz= " << p.pz() << " e= " << p.e()
0448 << " beta= " << p.Beta() << " gamma= " << p.Gamma()
0449 << " Proper time= " << proper_time / CLHEP::ns << " ns";
0450 }
0451
0452
0453
0454 double x1 = vp->end_vertex()->position().x();
0455 double y1 = vp->end_vertex()->position().y();
0456 double z1 = vp->end_vertex()->position().z();
0457
0458 for (const HepMC3::GenParticlePtr &vpdec : vp->end_vertex()->particles_out()) {
0459
0460 math::XYZTLorentzVector pdec(
0461 vpdec->momentum().px(), vpdec->momentum().py(), vpdec->momentum().pz(), vpdec->momentum().e());
0462
0463
0464 G4PrimaryParticle *g4daught =
0465 new G4PrimaryParticle(vpdec->pid(), pdec.x() * CLHEP::GeV, pdec.y() * CLHEP::GeV, pdec.z() * CLHEP::GeV);
0466
0467 if (g4daught->GetG4code() != nullptr) {
0468 g4daught->SetMass(g4daught->GetG4code()->GetPDGMass());
0469 g4daught->SetCharge(g4daught->GetG4code()->GetPDGCharge());
0470 }
0471
0472
0473
0474 setGenId(g4daught, vpdec->id());
0475
0476 int status = vpdec->status();
0477 if (verbose > 1)
0478 LogDebug("SimG4CoreGenerator3::::particleAssignDaughters")
0479 << "Assigning a " << vpdec->pid() << " as daughter of a " << vp->pid() << " status=" << status;
0480
0481 bool isInList;
0482 if (fSlepton) {
0483 std::vector<int> fParticleList = {1000011, 1000013, 1000015, 2000011, 2000013, 2000015};
0484 isInList = (std::find(fParticleList.begin(), fParticleList.end(), std::abs(vp->pdg_id())) != fParticleList.end());
0485 } else {
0486 isInList = std::abs(vp->pdg_id()) == 1000015;
0487 }
0488 if ((status == 2 || (status == 23 && isInList) || (status > 50 && status < 100)) &&
0489 vpdec->end_vertex() != nullptr) {
0490 double x2 = vpdec->end_vertex()->position().x();
0491 double y2 = vpdec->end_vertex()->position().y();
0492 double z2 = vpdec->end_vertex()->position().z();
0493 double dd = std::sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2) + (z1 - z2) * (z1 - z2));
0494 particleAssignDaughters(g4daught, vpdec.get(), dd);
0495 }
0496 vpdec->set_status(1000 + status);
0497 g4p->SetDaughter(g4daught);
0498
0499 if (verbose > 1)
0500 g4daught->Print();
0501 }
0502 }
0503
0504
0505 bool Generator3::particlePassesPrimaryCuts(const G4ThreeVector &p) const {
0506 bool flag = true;
0507 double ptot = p.mag();
0508 if (fPCuts && (ptot < theMinPCut * CLHEP::GeV || ptot > theMaxPCut * CLHEP::GeV)) {
0509 flag = false;
0510 }
0511 if (fEtaCuts && flag) {
0512 double pz = p.z();
0513 if (ptot < pz + 1.e-10) {
0514 flag = false;
0515
0516 } else {
0517 double eta = 0.5 * G4Log((ptot + pz) / (ptot - pz));
0518 if (eta < theMinEtaCut || eta > theMaxEtaCut) {
0519 flag = false;
0520 }
0521 }
0522 }
0523 if (fPhiCuts && flag) {
0524 double phi = p.phi();
0525 if (phi < theMinPhiCut || phi > theMaxPhiCut) {
0526 flag = false;
0527 }
0528 }
0529
0530 if (verbose > 2)
0531 LogDebug("SimG4CoreGenerator3") << "Generator ptot(GeV)= " << ptot / CLHEP::GeV << " eta= " << p.eta()
0532 << " phi= " << p.phi() << " Flag= " << flag;
0533
0534 return flag;
0535 }
0536
0537 bool Generator3::isExotic(int pdgcode) const {
0538 int pdgid = std::abs(pdgcode);
0539 return ((pdgid >= 1000000 && pdgid < 4000000 && pdgid != 3000022) ||
0540 pdgid == 17 ||
0541 pdgid == 34 ||
0542 pdgid == 37);
0543 }
0544
0545 bool Generator3::isExoticNonDetectable(int pdgcode) const {
0546 int pdgid = std::abs(pdgcode);
0547 HepPDT::ParticleID pid(pdgcode);
0548 int charge = pid.threeCharge();
0549 return ((charge == 0) && (pdgid >= 1000000 && pdgid < 1000040));
0550 }
0551
0552 bool Generator3::IsInTheFilterList(int pdgcode) const {
0553 int pdgid = std::abs(pdgcode);
0554 for (auto &pdg : pdgFilter) {
0555 if (std::abs(pdg) == pdgid) {
0556 return true;
0557 }
0558 }
0559 return false;
0560 }
0561
0562 void Generator3::nonCentralEvent2G4(const HepMC3::GenEvent *evt, G4Event *g4evt) {
0563 #if 0
0564
0565 int i = g4evt->GetNumberOfPrimaryVertex();
0566 for (HepMC3::GenEvent::particle_const_iterator it = evt->particles_begin(); it != evt->particles_end(); ++it) {
0567 ++i;
0568 HepMC3::GenParticle *gp = (*it);
0569
0570
0571 if (gp->status() != 1)
0572 continue;
0573
0574 int pdg = gp->pdg_id();
0575 G4PrimaryParticle *g4p = new G4PrimaryParticle(
0576 pdg, gp->momentum().px() * CLHEP::GeV, gp->momentum().py() * CLHEP::GeV, gp->momentum().pz() * CLHEP::GeV);
0577 if (g4p->GetG4code() != nullptr) {
0578 g4p->SetMass(g4p->GetG4code()->GetPDGMass());
0579 g4p->SetCharge(g4p->GetG4code()->GetPDGCharge());
0580 }
0581 setGenId(g4p, i);
0582 G4PrimaryVertex *v = new G4PrimaryVertex(gp->production_vertex()->position().x() * CLHEP::mm,
0583 gp->production_vertex()->position().y() * CLHEP::mm,
0584 gp->production_vertex()->position().z() * CLHEP::mm,
0585 gp->production_vertex()->position().t() * CLHEP::mm / CLHEP::c_light);
0586 v->SetPrimary(g4p);
0587 g4evt->AddPrimaryVertex(v);
0588 if (verbose > 0)
0589 v->Print();
0590 }
0591
0592 #endif
0593 }