Generators/src/Generator3.cc

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       theMinPhiCut(p.getParameter<double>("MinPhiCut")),  // in radians (CMS standard)
0031       theMaxPhiCut(p.getParameter<double>("MaxPhiCut")),
0032       theMinEtaCut(p.getParameter<double>("MinEtaCut")),
0033       theMaxEtaCut(p.getParameter<double>("MaxEtaCut")),
0034       theMinPCut(p.getParameter<double>("MinPCut")),  // in GeV (CMS standard)
0035       theMaxPCut(p.getParameter<double>("MaxPCut")),
0036       theEtaCutForHector(p.getParameter<double>("EtaCutForHector")),
0037       verbose(p.getUntrackedParameter<int>("Verbosity", 0)),
0038       fLumiFilter(nullptr),
0039       evt_(nullptr),
0040       vtx_(nullptr),
0041       weight_(0),
0042       Z_lmin(0),
0043       Z_lmax(0),
0044       Z_hector(0),
0045       pdgFilterSel(false),
0046       fPDGFilter(false) {
0047   bool lumi = p.getParameter<bool>("ApplyLumiMonitorCuts");
0048   if (lumi) {
0049     fLumiFilter = new LumiMonitorFilter();
0050   }
0051
0052   double theRDecLenCut = p.getParameter<double>("RDecLenCut") * CLHEP::cm;
0053   theDecRCut2 = theRDecLenCut * theRDecLenCut;
0054
0055   theMinPtCut2 = theMinPCut * theMinPCut;
0056
0057   double theDecLenCut = p.getParameter<double>("LDecLenCut") * CLHEP::cm;
0058
0059   maxZCentralCMS = p.getParameter<double>("MaxZCentralCMS") * CLHEP::m;
0060
0061   fFiductialCuts = (fPCuts || fPtransCut || fEtaCuts || fPhiCuts);
0062
0063   pdgFilter.resize(0);
0064   if (p.exists("PDGselection")) {
0065     pdgFilterSel = (p.getParameter<edm::ParameterSet>("PDGselection")).getParameter<bool>("PDGfilterSel");
0066     pdgFilter = (p.getParameter<edm::ParameterSet>("PDGselection")).getParameter<std::vector<int>>("PDGfilter");
0067     if (!pdgFilter.empty()) {
0068       fPDGFilter = true;
0069       std::stringstream ss;
0070       ss << "SimG4Core/Generator3: ";
0071       if (pdgFilterSel) {
0072         ss << " Selecting only PDG ID = ";
0073       } else {
0074         ss << " Filtering out PDG ID = ";
0075       }
0076       for (unsigned int ii = 0; ii < pdgFilter.size(); ++ii) {
0077         ss << pdgFilter[ii] << "  ";
0078       }
0079       edm::LogVerbatim("SimG4CoreGenerator3") << ss.str();
0080     }
0081   }
0082
0083   if (fEtaCuts) {
0084     Z_lmax = theRDecLenCut * ((1 - exp(-2 * theMaxEtaCut)) / (2 * exp(-theMaxEtaCut)));
0085     Z_lmin = theRDecLenCut * ((1 - exp(-2 * theMinEtaCut)) / (2 * exp(-theMinEtaCut)));
0086   }
0087
0088   Z_hector = theRDecLenCut * ((1 - exp(-2 * theEtaCutForHector)) / (2 * exp(-theEtaCutForHector)));
0089
0090   edm::LogVerbatim("SimG4CoreGenerator3")
0091       << "SimG4Core/Generator3: Rdecaycut= " << theRDecLenCut / CLHEP::cm
0092       << " cm;  Zdecaycut= " << theDecLenCut / CLHEP::cm << "Z_min= " << Z_lmin / CLHEP::cm
0093       << " cm; Z_max= " << Z_lmax / CLHEP::cm << " cm;\n"
0094       << "                     MaxZCentralCMS = " << maxZCentralCMS / CLHEP::m << " m;"
0095       << " Z_hector = " << Z_hector / CLHEP::cm << " cm\n"
0096       << "                     ApplyCuts: " << fFiductialCuts << "  PCuts: " << fPCuts << "  PtransCut: " << fPtransCut
0097       << "  EtaCut: " << fEtaCuts << "  PhiCut: " << fPhiCuts << "  LumiMonitorCut: " << lumi;
0098   if (fFiductialCuts) {
0099     edm::LogVerbatim("SimG4CoreGenerator3")
0100         << "SimG4Core/Generator3: Pmin(GeV)= " << theMinPCut << "; Pmax(GeV)= " << theMaxPCut
0101         << "; EtaMin= " << theMinEtaCut << "; EtaMax= " << theMaxEtaCut << "; PhiMin(rad)= " << theMinPhiCut
0102         << "; PhiMax(rad)= " << theMaxPhiCut;
0103   }
0104   if (lumi) {
0105     fLumiFilter->Describe();
0106   }
0107 }
0108
0109 Generator3::~Generator3() { delete fLumiFilter; }
0110
0111 void Generator3::HepMC2G4(const HepMC3::GenEvent *evt_orig, G4Event *g4evt) {
0112   HepMC3::GenEvent *evt = new HepMC3::GenEvent(*evt_orig);
0113   //HepMC3::Print::listing(*evt);
0114
0115   if ((evt->vertices()).empty()) {
0116     std::stringstream ss;
0117     ss << "SimG4Core/Generator3: in event " << g4evt->GetEventID() << " Corrupted Event - GenEvent with no vertex \n";
0118     throw SimG4Exception(ss.str());
0119   }
0120
0121   if (!evt->weights().empty()) {
0122     weight_ = evt->weights()[0];
0123     for (unsigned int iw = 1; iw < evt->weights().size(); ++iw) {
0124       // terminate if the vector of weights contains a zero-weight
0125       if (evt->weights()[iw] <= 0)
0126         break;
0127       weight_ *= evt->weights()[iw];
0128     }
0129   }
0130
0131   if (vtx_ != nullptr) {
0132     delete vtx_;
0133   }
0134
0135   for (HepMC3::GenVertexPtr v : evt->vertices()) {
0136     vtx_ =
0137         new math::XYZTLorentzVector((v->position()).x(), (v->position()).y(), (v->position()).z(), (v->position()).t());
0138     break;
0139   }
0140
0141   edm::LogVerbatim("SimG4CoreGenerator3")
0142       << "Generator3: primary Vertex = (" << vtx_->x() << ", " << vtx_->y() << ", " << vtx_->z() << ")";
0143
0144   unsigned int ng4vtx = 0;
0145   unsigned int ng4par = 0;
0146
0147   for (HepMC3::GenVertexPtr vitr : evt->vertices()) {
0148     // loop for vertex, is it a real vertex?
0149     // Set qvtx to true for any particles that should be propagated by GEANT,
0150     // i.e., status 1 particles or status 2 particles that decay outside the
0151     // beampipe.
0152     G4bool qvtx = false;
0153     for (HepMC3::GenParticlePtr pitr : vitr->particles_out()) {
0154       // For purposes of this function, the status is defined as follows:
0155       // 1:  particles are not decayed by generator
0156       // 2:  particles are decayed by generator but need to be propagated by GEANT
0157       // 3:  particles are decayed by generator and do not need to be propagated by GEANT
0158       int status = pitr->status();
0159       int pdg = pitr->pid();
0160       if (status > 3 && isExotic(pdg) && (!(isExoticNonDetectable(pdg)))) {
0161         // In Pythia 8, there are many status codes besides 1, 2, 3.
0162         // By setting the status to 2 for exotic particles, they will be
0163         // checked: if its decay vertex is outside the beampipe, it will be
0164         // propagated by GEANT. Some Standard Model particles, e.g., K0, cannot
0165         // be propagated by GEANT, so do not change their status code.
0166         status = 2;
0167       }
0168       if (status == 2 && abs(pdg) == 9900015) {  // Additional photon?
0169         status = 3;
0170       }
0171
0172       // Particles which are not decayed by generator
0173       if (status == 1) {
0174         // filter out unwanted particles and vertices
0175         if (fPDGFilter && !pdgFilterSel && IsInTheFilterList(pdg)) {
0176           continue;
0177         }
0178
0179         qvtx = true;
0180         if (verbose > 2)
0181           LogDebug("SimG4CoreGenerator3")
0182               << "GenVertex barcode = " << vitr->id() << " " << qvtx
0183               << " selected for GenParticle barcode = " << pitr->id();  // barcode is substituted by id
0184         break;
0185       }
0186       // The selection is made considering if the partcile with status = 2
0187       // have the end_vertex with a radius greater than the radius of beampipe
0188       // cylinder (no requirement on the Z of the vertex is applyed).
0189       else if (status == 2) {
0190         if (pitr->end_vertex() != nullptr) {
0191           double xx = (pitr->end_vertex())->position().x();
0192           double yy = (pitr->end_vertex())->position().y();
0193           double r_dd = xx * xx + yy * yy;
0194           if (r_dd > theDecRCut2) {
0195             qvtx = true;
0196             if (verbose > 2)
0197               LogDebug("SimG4CoreGenerator3")
0198                   << "GenVertex barcode = " << vitr->id() << " selected for GenParticle barcode = " << pitr->id()
0199                   << " radius = " << std::sqrt(r_dd);  // barcode is substituted by id
0200             break;
0201           }
0202         } else {
0203           // particles with status 2 without end_vertex are
0204           // equivalent to stable
0205           qvtx = true;
0206           break;
0207         }
0208       }
0209     }
0210
0211     // if this vertex is inside fiductial volume inside the beam pipe
0212     // and has no long-lived secondary the vertex is not saved
0213     if (!qvtx) {
0214       continue;
0215     }
0216
0217     double x1 = vitr->position().x() * CLHEP::mm;
0218     double y1 = vitr->position().y() * CLHEP::mm;
0219     double z1 = vitr->position().z() * CLHEP::mm;
0220     double t1 = vitr->position().t() * CLHEP::mm / CLHEP::c_light;
0221
0222     G4PrimaryVertex *g4vtx = new G4PrimaryVertex(x1, y1, z1, t1);
0223
0224     for (HepMC3::GenParticlePtr pitr : vitr->particles_out()) {
0225       int status = pitr->status();
0226       int pdg = pitr->pid();
0227       bool hasDecayVertex = (nullptr != pitr->end_vertex());
0228
0229       // Filter on allowed particle species if required
0230       if (fPDGFilter) {
0231         bool isInTheList = IsInTheFilterList(pdg);
0232         if ((!pdgFilterSel && isInTheList) || (pdgFilterSel && !isInTheList)) {
0233           if (0 < verbose)
0234             edm::LogVerbatim("SimG4CoreGenerator3")
0235                 << " Skiped GenParticle barcode= " << pitr->id() << " PDGid= " << pdg << " status= " << status
0236                 << " isExotic: " << isExotic(pdg) << " isExoticNotDet: " << isExoticNonDetectable(pdg)
0237                 << " isInTheList: " << isInTheList << " hasDecayVertex: " << hasDecayVertex;
0238           continue;
0239         }
0240       }
0241
0242       if (0 < verbose) {
0243         edm::LogVerbatim("SimG4CoreGenerator3")
0244             << "Generator3: pdg= " << pdg << " status= " << status << " hasPreDefinedDecay: " << hasDecayVertex
0245             << "\n           isExotic: " << isExotic(pdg) << " isExoticNotDet: " << isExoticNonDetectable(pdg)
0246             << " isInTheList: " << IsInTheFilterList(pdg) << "\n"
0247             << " MaxZCentralCMS = " << maxZCentralCMS / CLHEP::m << " m;  (x,y,z,t): (" << x1 << "," << y1 << "," << z1
0248             << "," << t1 << ")";
0249       }
0250       if (status > 3 && isExotic(pdg) && (!(isExoticNonDetectable(pdg)))) {
0251         status = hasDecayVertex ? 2 : 1;
0252       }
0253       if (status == 2 && abs(pdg) == 9900015) {  // Additional photon?
0254         status = 3;
0255       }
0256
0257       // this particle has predefined decay but has no vertex
0258       if (2 == status && !hasDecayVertex) {
0259         edm::LogWarning("SimG4CoreGenerator3: in event ")
0260             << g4evt->GetEventID() << " a particle "
0261             << " pdgid= " << pdg << " has status=2 but has no decay vertex, so will be fully tracked by Geant4";
0262         status = 1;
0263       }
0264
0265       double x2 = x1;
0266       double y2 = y1;
0267       double z2 = z1;
0268       double decay_length = 0.0;
0269       if (2 == status) {
0270         x2 = pitr->end_vertex()->position().x();
0271         y2 = pitr->end_vertex()->position().y();
0272         z2 = pitr->end_vertex()->position().z();
0273         decay_length = std::sqrt((x2 - x1) * (x2 - x1) + (y2 - y1) * (y2 - y1) + (z2 - z1) * (z2 - z1));
0274       }
0275
0276       bool toBeAdded = !fFiductialCuts;
0277
0278       double px = pitr->momentum().px();
0279       double py = pitr->momentum().py();
0280       double pz = pitr->momentum().pz();
0281       double ptot = std::sqrt(px * px + py * py + pz * pz);
0282       math::XYZTLorentzVector p(px, py, pz, pitr->momentum().e());
0283
0284       double ximpact = x1;
0285       double yimpact = y1;
0286       double zimpact = z1;
0287
0288       // protection against numerical problems for extremely low momenta
0289       // compute impact point at transition to Hector
0290       const double minTan = 1.e-20;
0291       if (std::abs(z1) < Z_hector && std::abs(pz) >= minTan * ptot) {
0292         zimpact = (pz > 0.0) ? Z_hector : -Z_hector;
0293         double del = (zimpact - z1) / pz;
0294         ximpact += del * px;
0295         yimpact += del * py;
0296       }
0297       double rimpact2 = ximpact * ximpact + yimpact * yimpact;
0298
0299       if (verbose > 2)
0300         LogDebug("SimG4CoreGenerator3") << "Processing GenParticle barcode= " << pitr->id() << " pdg= " << pdg
0301                                         << " status= " << pitr->status() << " st= " << status
0302                                         << " rimpact(cm)= " << std::sqrt(rimpact2) / CLHEP::cm
0303                                         << " zimpact(cm)= " << zimpact / CLHEP::cm << " ptot(GeV)= " << ptot
0304                                         << " pz(GeV)= " << pz;
0305
0306       // Particles of status 1 trasnported along the beam pipe
0307       // HECTOR transport of protons are done in corresponding PPS producer
0308       if (1 == status && std::abs(zimpact) >= Z_hector && rimpact2 <= theDecRCut2) {
0309         // very forward n, nbar, gamma are allowed
0310         toBeAdded = (2112 == std::abs(pdg) || 22 == pdg);
0311         if (verbose > 1) {
0312           edm::LogVerbatim("SimG4CoreGenerator3")
0313               << "GenParticle barcode = " << pitr->id() << " very forward; to be added: " << toBeAdded;
0314         }
0315       } else {
0316         // Standard case: particles not decayed by the generator and not forward
0317         if (1 == status && (std::abs(zimpact) < Z_hector || rimpact2 > theDecRCut2)) {
0318           // Ptot cut for all particles
0319           if (fPCuts && (ptot < theMinPCut || ptot > theMaxPCut)) {
0320             continue;
0321           }
0322           // phi cut is applied mainly for particle gun
0323           if (fPhiCuts) {
0324             double phi = p.phi();
0325             if (phi < theMinPhiCut || phi > theMaxPhiCut) {
0326               continue;
0327             }
0328           }
0329           // eta cut is applied if position of the decay
0330           // is within vacuum chamber and limited in Z
0331           if (fEtaCuts) {
0332             // eta cut
0333             double xi = x1;
0334             double yi = y1;
0335             double zi = z1;
0336
0337             // can be propagated along Z
0338             if (std::abs(pz) >= minTan * ptot) {
0339               if ((zi >= Z_lmax) & (pz < 0.0)) {
0340                 zi = Z_lmax;
0341               } else if ((zi <= Z_lmin) & (pz > 0.0)) {
0342                 zi = Z_lmin;
0343               } else {
0344                 if (pz > 0) {
0345                   zi = Z_lmax;
0346                 } else {
0347                   zi = Z_lmin;
0348                 }
0349               }
0350               double del = (zi - z1) / pz;
0351               xi += del * px;
0352               yi += del * py;
0353             }
0354             // check eta cut
0355             if ((zi >= Z_lmin) & (zi <= Z_lmax) & (xi * xi + yi * yi < theDecRCut2)) {
0356               continue;
0357             }
0358           }
0359           const HepMC3::GenParticle *ppointer = pitr.get();
0360           if (fLumiFilter && !fLumiFilter->isGoodForLumiMonitor(ppointer)) {  // MK: this function is always true
0361             continue;
0362           }
0363           toBeAdded = true;
0364           if (verbose > 1)
0365             edm::LogVerbatim("SimG4CoreGenerator3") << "GenParticle barcode = " << pitr->id() << " passed case 1";
0366
0367           // Decay chain outside the fiducial cylinder defined by theRDecLenCut
0368           // are used for Geant4 tracking with predefined decay channel
0369           // In the case of decay in vacuum particle is not tracked by Geant4
0370         } else if (2 == status && x2 * x2 + y2 * y2 >= theDecRCut2 && std::abs(z2) < Z_hector) {
0371           toBeAdded = true;
0372           if (verbose > 1)
0373             edm::LogVerbatim("SimG4CoreGenerator3") << "GenParticle barcode = " << pitr->id() << " passed case 2"
0374                                                     << " decay_length(cm)= " << decay_length / CLHEP::cm;
0375         }
0376       }
0377       if (toBeAdded) {
0378         G4PrimaryParticle *g4prim = new G4PrimaryParticle(pdg, px * CLHEP::GeV, py * CLHEP::GeV, pz * CLHEP::GeV);
0379
0380         if (g4prim->GetG4code() != nullptr) {
0381           g4prim->SetMass(g4prim->GetG4code()->GetPDGMass());
0382           double charge = g4prim->GetG4code()->GetPDGCharge();
0383
0384           // apply Pt cut
0385           if (fPtransCut && 1 == status && 0.0 != charge && px * px + py * py < theMinPtCut2) {
0386             delete g4prim;
0387             continue;
0388           }
0389           g4prim->SetCharge(charge);
0390         }
0391
0392         // V.I. do not use SetWeight but the same code
0393         // value of the code compute inside TrackWithHistory
0394         // g4prim->SetWeight( 10000*(*vpitr)->barcode() ) ;
0395         setGenId(g4prim, pitr->id());
0396
0397         if (2 == status) {
0398           particleAssignDaughters(g4prim, pitr.get(), decay_length);
0399         }
0400         if (verbose > 1)
0401           g4prim->Print();
0402
0403         ++ng4par;
0404         g4vtx->SetPrimary(g4prim);
0405         edm::LogVerbatim("SimG4CoreGenerator3") << "   " << ng4par << ". new Geant4 particle pdg= " << pdg
0406                                                 << " Ptot(GeV/c)= " << ptot << " Pt= " << std::sqrt(px * px + py * py)
0407                                                 << " status= " << status << "; dir= " << g4prim->GetMomentumDirection();
0408       }
0409     }
0410
0411     if (verbose > 1)
0412       g4vtx->Print();
0413     g4evt->AddPrimaryVertex(g4vtx);
0414     ++ng4vtx;
0415   }
0416
0417   // Add a protection for completely empty events (produced by LHCTransport):
0418   // add a dummy vertex with no particle attached to it
0419   if (ng4vtx == 0) {
0420     G4PrimaryVertex *g4vtx = new G4PrimaryVertex(0.0, 0.0, 0.0, 0.0);
0421     if (verbose > 1)
0422       g4vtx->Print();
0423
0424     g4evt->AddPrimaryVertex(g4vtx);
0425   }
0426
0427   edm::LogVerbatim("SimG4CoreGenerator3")
0428       << "The list of Geant4 primaries includes " << ng4par << " particles in " << ng4vtx << " vertex";
0429
0430   delete evt;
0431 }
0432
0433 void Generator3::particleAssignDaughters(G4PrimaryParticle *g4p, HepMC3::GenParticle *vp, double decaylength) {
0434   if (verbose > 1) {  // MK: Check message below
0435     LogDebug("SimG4CoreGenerator3") << "Special case of long decay length \n"
0436                                     << "Assign daughters to mother with decaylength=" << decaylength / CLHEP::cm
0437                                     << " cm";
0438   }
0439   math::XYZTLorentzVector p(vp->momentum().px(), vp->momentum().py(), vp->momentum().pz(), vp->momentum().e());
0440
0441   // defined preassigned decay time
0442   double proper_time = decaylength / (p.Beta() * p.Gamma() * c_light);
0443   g4p->SetProperTime(proper_time);
0444
0445   if (verbose > 2) {
0446     LogDebug("SimG4CoreGenerator3") << " px= " << p.px() << " py= " << p.py() << " pz= " << p.pz() << " e= " << p.e()
0447                                     << " beta= " << p.Beta() << " gamma= " << p.Gamma()
0448                                     << " Proper time= " << proper_time / CLHEP::ns << " ns";
0449   }
0450
0451   // the particle will decay after the same length if it
0452   // has not interacted before
0453   double x1 = vp->end_vertex()->position().x();
0454   double y1 = vp->end_vertex()->position().y();
0455   double z1 = vp->end_vertex()->position().z();
0456
0457   for (HepMC3::GenParticlePtr vpdec : vp->end_vertex()->particles_out()) {
0458     // transform decay products such that in the rest frame of mother
0459     math::XYZTLorentzVector pdec(
0460         vpdec->momentum().px(), vpdec->momentum().py(), vpdec->momentum().pz(), vpdec->momentum().e());
0461
0462     // children should only be taken into account once
0463     G4PrimaryParticle *g4daught =
0464         new G4PrimaryParticle(vpdec->pid(), 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     // V.I. do not use SetWeight but the same code
0472     // value of the code compute inside TrackWithHistory
0473     setGenId(g4daught, vpdec->id());
0474
0475     int status = vpdec->status();
0476     if (verbose > 1)
0477       LogDebug("SimG4CoreGenerator3::::particleAssignDaughters")
0478           << "Assigning a " << vpdec->pid() << " as daughter of a " << vp->pid() << " status=" << status;
0479
0480     if ((status == 2 || (status == 23 && std::abs(vp->pid()) == 1000015) || (status > 50 && status < 100)) &&
0481         vpdec->end_vertex() != nullptr) {
0482       double x2 = vpdec->end_vertex()->position().x();
0483       double y2 = vpdec->end_vertex()->position().y();
0484       double z2 = vpdec->end_vertex()->position().z();
0485       double dd = std::sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2) + (z1 - z2) * (z1 - z2));
0486       particleAssignDaughters(g4daught, vpdec.get(), dd);
0487     }
0488     vpdec->set_status(1000 + status);
0489     g4p->SetDaughter(g4daught);
0490
0491     if (verbose > 1)
0492       g4daught->Print();
0493   }
0494 }
0495
0496 // Used for non-beam particles
0497 bool Generator3::particlePassesPrimaryCuts(const G4ThreeVector &p) const {
0498   bool flag = true;
0499   double ptot = p.mag();
0500   if (fPCuts && (ptot < theMinPCut * CLHEP::GeV || ptot > theMaxPCut * CLHEP::GeV)) {
0501     flag = false;
0502   }
0503   if (fEtaCuts && flag) {
0504     double pz = p.z();
0505     if (ptot < pz + 1.e-10) {
0506       flag = false;
0507
0508     } else {
0509       double eta = 0.5 * G4Log((ptot + pz) / (ptot - pz));
0510       if (eta < theMinEtaCut || eta > theMaxEtaCut) {
0511         flag = false;
0512       }
0513     }
0514   }
0515   if (fPhiCuts && flag) {
0516     double phi = p.phi();
0517     if (phi < theMinPhiCut || phi > theMaxPhiCut) {
0518       flag = false;
0519     }
0520   }
0521
0522   if (verbose > 2)
0523     LogDebug("SimG4CoreGenerator3") << "Generator ptot(GeV)= " << ptot / CLHEP::GeV << " eta= " << p.eta()
0524                                     << "  phi= " << p.phi() << " Flag= " << flag;
0525
0526   return flag;
0527 }
0528
0529 bool Generator3::isExotic(int pdgcode) const {
0530   int pdgid = std::abs(pdgcode);
0531   return ((pdgid >= 1000000 && pdgid < 4000000 && pdgid != 3000022) ||  // SUSY, R-hadron, and technicolor particles
0532           pdgid == 17 ||                                                // 4th generation lepton
0533           pdgid == 34 ||                                                // W-prime
0534           pdgid == 37);                                                 // charged Higgs
0535 }
0536
0537 bool Generator3::isExoticNonDetectable(int pdgcode) const {
0538   int pdgid = std::abs(pdgcode);
0539   HepPDT::ParticleID pid(pdgcode);
0540   int charge = pid.threeCharge();
0541   return ((charge == 0) && (pdgid >= 1000000 && pdgid < 1000040));  // SUSY
0542 }
0543
0544 bool Generator3::IsInTheFilterList(int pdgcode) const {
0545   int pdgid = std::abs(pdgcode);
0546   for (auto &pdg : pdgFilter) {
0547     if (std::abs(pdg) == pdgid) {
0548       return true;
0549     }
0550   }
0551   return false;
0552 }
0553
0554 void Generator3::nonCentralEvent2G4(const HepMC3::GenEvent *evt, G4Event *g4evt) {
0555 #if 0
0556
0557   int i = g4evt->GetNumberOfPrimaryVertex();
0558   for (HepMC3::GenEvent::particle_const_iterator it = evt->particles_begin(); it != evt->particles_end(); ++it) {
0559     ++i;
0560     HepMC3::GenParticle *gp = (*it);
0561
0562     // storing only particle with status == 1
0563     if (gp->status() != 1)
0564       continue;
0565
0566     int pdg = gp->pdg_id();
0567     G4PrimaryParticle *g4p = new G4PrimaryParticle(
0568         pdg, gp->momentum().px() * CLHEP::GeV, gp->momentum().py() * CLHEP::GeV, gp->momentum().pz() * CLHEP::GeV);
0569     if (g4p->GetG4code() != nullptr) {
0570       g4p->SetMass(g4p->GetG4code()->GetPDGMass());
0571       g4p->SetCharge(g4p->GetG4code()->GetPDGCharge());
0572     }
0573     setGenId(g4p, i);
0574     G4PrimaryVertex *v = new G4PrimaryVertex(gp->production_vertex()->position().x() * CLHEP::mm,
0575                                              gp->production_vertex()->position().y() * CLHEP::mm,
0576                                              gp->production_vertex()->position().z() * CLHEP::mm,
0577                                              gp->production_vertex()->position().t() * CLHEP::mm / CLHEP::c_light);
0578     v->SetPrimary(g4p);
0579     g4evt->AddPrimaryVertex(v);
0580     if (verbose > 0)
0581       v->Print();
0582   }  // end loop on HepMC particles
0583
0584 #endif
0585 }