Back to home page

Project CMSSW displayed by LXR

 
 

    


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")),  // in radians (CMS standard)
0030       theMaxPhiCut(p.getParameter<double>("MaxPhiCut")),
0031       theMinEtaCut(p.getParameter<double>("MinEtaCut")),
0032       theMaxEtaCut(p.getParameter<double>("MaxEtaCut")),
0033       theMinPCut(p.getParameter<double>("MinPCut")),  // in GeV (CMS standard)
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       // terminate if the vector of weights contains a zero-weight
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     // loop for vertex, is it a real vertex?
0148     // Set qvtx to true for any particles that should be propagated by GEANT,
0149     // i.e., status 1 particles or status 2 particles that decay outside the
0150     // beampipe.
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       // 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)->pdg_id();
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) {
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("SimG4CoreGenerator") << "GenVertex barcode = " << (*vitr)->barcode()
0182                                          << " selected for GenParticle barcode = " << (*pitr)->barcode();
0183         break;
0184       }
0185       // The selection is made considering if the partcile with status = 2
0186       // have the end_vertex with a radius greater than the radius of beampipe
0187       // cylinder (no requirement on the Z of the vertex is applyed).
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           // particles with status 2 without end_vertex are
0203           // equivalent to stable
0204           qvtx = true;
0205           break;
0206         }
0207       }
0208     }
0209 
0210     // if this vertex is inside fiductial volume inside the beam pipe
0211     // and has no long-lived secondary the vertex is not saved
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       // Filter on allowed particle species if required
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       // this particle has predefined decay but has no vertex
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       // protection against numerical problems for extremely low momenta
0288       // compute impact point at transition to Hector
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       // Particles of status 1 trasnported along the beam pipe
0306       // HECTOR transport of protons are done in corresponding PPS producer
0307       if (1 == status && std::abs(zimpact) >= Z_hector && rimpact2 <= theDecRCut2) {
0308         // very forward n, nbar, gamma are allowed
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         // Standard case: particles not decayed by the generator and not forward
0316         if (1 == status && (std::abs(zimpact) < Z_hector || rimpact2 > theDecRCut2)) {
0317           // Ptot cut for all particles
0318           if (fPCuts && (ptot < theMinPCut || ptot > theMaxPCut)) {
0319             continue;
0320           }
0321           // phi cut is applied mainly for particle gun
0322           if (fPhiCuts) {
0323             double phi = p.phi();
0324             if (phi < theMinPhiCut || phi > theMaxPhiCut) {
0325               continue;
0326             }
0327           }
0328           // eta cut is applied if position of the decay
0329           // is within vacuum chamber and limited in Z
0330           if (fEtaCuts) {
0331             // eta cut
0332             double xi = x1;
0333             double yi = y1;
0334             double zi = z1;
0335 
0336             // can be propagated along Z
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             // check eta cut
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           // Decay chain outside the fiducial cylinder defined by theRDecLenCut
0367           // are used for Geant4 tracking with predefined decay channel
0368           // In the case of decay in vacuum particle is not tracked by Geant4
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           // apply Pt cut
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         // V.I. do not use SetWeight but the same code
0392         // value of the code compute inside TrackWithHistory
0393         // g4prim->SetWeight( 10000*(*vpitr)->barcode() ) ;
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   // Add a protection for completely empty events (produced by LHCTransport):
0417   // add a dummy vertex with no particle attached to it
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   // defined preassigned decay time
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   // the particle will decay after the same length if it
0450   // has not interacted before
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     // 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)->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     // V.I. do not use SetWeight but the same code
0472     // value of the code compute inside TrackWithHistory
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 // Used for non-beam particles
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) ||  // SUSY, R-hadron, and technicolor particles
0539           pdgid == 17 ||                                                // 4th generation lepton
0540           pdgid == 34 ||                                                // W-prime
0541           pdgid == 37);                                                 // charged Higgs
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));  // SUSY
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     // storing only particle with status == 1
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   }  // end loop on HepMC particles
0588 }