Generators/src/Generator.cc

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