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