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