CMS 3D CMS Logo

List of all members | Public Member Functions | Private Member Functions | Private Attributes
Generator Class Reference

#include <Generator.h>

Public Member Functions

virtual const double eventWeight () const
 
 Generator (const edm::ParameterSet &p)
 
virtual const HepMC::GenEventgenEvent () const
 
virtual const math::XYZTLorentzVectorgenVertex () const
 
void HepMC2G4 (const HepMC::GenEvent *g, G4Event *e)
 
void nonCentralEvent2G4 (const HepMC::GenEvent *g, G4Event *e)
 
void setGenEvent (const HepMC::GenEvent *inpevt)
 
virtual ~Generator ()
 

Private Member Functions

bool isExotic (int pdgcode) const
 
bool isExoticNonDetectable (int pdgcode) const
 
bool IsInTheFilterList (int pdgcode) const
 
void particleAssignDaughters (G4PrimaryParticle *p, HepMC::GenParticle *hp, double length)
 
bool particlePassesPrimaryCuts (const G4ThreeVector &p) const
 
void setGenId (G4PrimaryParticle *p, int id) const
 

Private Attributes

HepMC::GenEventevt_
 
bool fEtaCuts
 
bool fFiductialCuts
 
LumiMonitorFilterfLumiFilter
 
bool fPCuts
 
bool fPDGFilter
 
bool fPhiCuts
 
bool fPtransCut
 
double maxZCentralCMS
 
std::vector< int > pdgFilter
 
bool pdgFilterSel
 
double theDecLenCut
 
double theDecRCut2
 
double theEtaCutForHector
 
double theMaxEtaCut
 
double theMaxPCut
 
double theMaxPhiCut
 
double theMinEtaCut
 
double theMinPCut
 
double theMinPhiCut
 
double theMinPtCut2
 
int verbose
 
math::XYZTLorentzVectorvtx_
 
double weight_
 
double Z_hector
 
double Z_lmax
 
double Z_lmin
 

Detailed Description

Definition at line 19 of file Generator.h.

Constructor & Destructor Documentation

◆ Generator()

Generator::Generator ( const edm::ParameterSet p)

Definition at line 23 of file Generator.cc.

References LumiMonitorFilter::Describe(), JetChargeProducer_cfi::exp, fEtaCuts, fFiductialCuts, fLumiFilter, fPCuts, fPDGFilter, fPhiCuts, fPtransCut, cuy::ii, BXlumiParameters_cfi::lumi, visualization-live-secondInstance_cfg::m, maxZCentralCMS, AlCaHLTBitMon_ParallelJobs::p, pdgFilter, pdgFilterSel, contentValuesCheck::ss, theDecLenCut, theDecRCut2, theEtaCutForHector, theMaxEtaCut, theMaxPCut, theMaxPhiCut, theMinEtaCut, theMinPCut, theMinPhiCut, theMinPtCut2, Z_hector, Z_lmax, and Z_lmin.

24  : fPCuts(p.getParameter<bool>("ApplyPCuts")),
25  fPtransCut(p.getParameter<bool>("ApplyPtransCut")),
26  fEtaCuts(p.getParameter<bool>("ApplyEtaCuts")),
27  fPhiCuts(p.getParameter<bool>("ApplyPhiCuts")),
28  theMinPhiCut(p.getParameter<double>("MinPhiCut")), // in radians (CMS standard)
29  theMaxPhiCut(p.getParameter<double>("MaxPhiCut")),
30  theMinEtaCut(p.getParameter<double>("MinEtaCut")),
31  theMaxEtaCut(p.getParameter<double>("MaxEtaCut")),
32  theMinPCut(p.getParameter<double>("MinPCut")), // in GeV (CMS standard)
33  theMaxPCut(p.getParameter<double>("MaxPCut")),
34  theEtaCutForHector(p.getParameter<double>("EtaCutForHector")),
35  verbose(p.getUntrackedParameter<int>("Verbosity", 0)),
36  fLumiFilter(nullptr),
37  evt_(nullptr),
38  vtx_(nullptr),
39  weight_(0),
40  Z_lmin(0),
41  Z_lmax(0),
42  Z_hector(0),
43  pdgFilterSel(false),
44  fPDGFilter(false) {
45  bool lumi = p.getParameter<bool>("ApplyLumiMonitorCuts");
46  if (lumi) {
48  }
49 
50  double theRDecLenCut = p.getParameter<double>("RDecLenCut") * CLHEP::cm;
51  theDecRCut2 = theRDecLenCut * theRDecLenCut;
52 
54 
55  double theDecLenCut = p.getParameter<double>("LDecLenCut") * CLHEP::cm;
56 
57  maxZCentralCMS = p.getParameter<double>("MaxZCentralCMS") * CLHEP::m;
58 
60 
61  pdgFilter.resize(0);
62  if (p.exists("PDGselection")) {
63  pdgFilterSel = (p.getParameter<edm::ParameterSet>("PDGselection")).getParameter<bool>("PDGfilterSel");
64  pdgFilter = (p.getParameter<edm::ParameterSet>("PDGselection")).getParameter<std::vector<int>>("PDGfilter");
65  if (!pdgFilter.empty()) {
66  fPDGFilter = true;
67  std::stringstream ss;
68  ss << "SimG4Core/Generator: ";
69  if (pdgFilterSel) {
70  ss << " Selecting only PDG ID = ";
71  } else {
72  ss << " Filtering out PDG ID = ";
73  }
74  for (unsigned int ii = 0; ii < pdgFilter.size(); ++ii) {
75  ss << pdgFilter[ii] << " ";
76  }
77  edm::LogVerbatim("SimG4CoreGenerator") << ss.str();
78  }
79  }
80 
81  if (fEtaCuts) {
82  Z_lmax = theRDecLenCut * ((1 - exp(-2 * theMaxEtaCut)) / (2 * exp(-theMaxEtaCut)));
83  Z_lmin = theRDecLenCut * ((1 - exp(-2 * theMinEtaCut)) / (2 * exp(-theMinEtaCut)));
84  }
85 
86  Z_hector = theRDecLenCut * ((1 - exp(-2 * theEtaCutForHector)) / (2 * exp(-theEtaCutForHector)));
87 
88  edm::LogVerbatim("SimG4CoreGenerator") << "SimG4Core/Generator: Rdecaycut= " << theRDecLenCut / CLHEP::cm
89  << " cm; Zdecaycut= " << theDecLenCut / CLHEP::cm
90  << "Z_min= " << Z_lmin / CLHEP::cm << " cm; Z_max= " << Z_lmax / CLHEP::cm
91  << " cm;\n"
92  << " MaxZCentralCMS = " << maxZCentralCMS / CLHEP::m
93  << " m;"
94  << " Z_hector = " << Z_hector / CLHEP::cm << " cm\n"
95  << " ApplyCuts: " << fFiductialCuts
96  << " PCuts: " << fPCuts << " PtransCut: " << fPtransCut
97  << " EtaCut: " << fEtaCuts << " PhiCut: " << fPhiCuts
98  << " LumiMonitorCut: " << lumi;
99  if (fFiductialCuts) {
100  edm::LogVerbatim("SimG4CoreGenerator")
101  << "SimG4Core/Generator: Pmin(GeV)= " << theMinPCut << "; Pmax(GeV)= " << theMaxPCut
102  << "; EtaMin= " << theMinEtaCut << "; EtaMax= " << theMaxEtaCut << "; PhiMin(rad)= " << theMinPhiCut
103  << "; PhiMax(rad)= " << theMaxPhiCut;
104  }
105  if (lumi) {
107  }
108 }
double theDecRCut2
Definition: Generator.h:55
Log< level::Info, true > LogVerbatim
bool fPhiCuts
Definition: Generator.h:46
double theMinPtCut2
Definition: Generator.h:53
double maxZCentralCMS
Definition: Generator.h:58
double theEtaCutForHector
Definition: Generator.h:56
double Z_lmax
Definition: Generator.h:64
double theMaxPhiCut
Definition: Generator.h:49
bool fPtransCut
Definition: Generator.h:44
double theDecLenCut
Definition: Generator.h:57
double theMinEtaCut
Definition: Generator.h:50
bool pdgFilterSel
Definition: Generator.h:66
double theMinPCut
Definition: Generator.h:52
LumiMonitorFilter * fLumiFilter
Definition: Generator.h:60
bool fFiductialCuts
Definition: Generator.h:47
bool fPDGFilter
Definition: Generator.h:67
double theMinPhiCut
Definition: Generator.h:48
double weight_
Definition: Generator.h:63
double theMaxEtaCut
Definition: Generator.h:51
ii
Definition: cuy.py:589
double theMaxPCut
Definition: Generator.h:54
math::XYZTLorentzVector * vtx_
Definition: Generator.h:62
HepMC::GenEvent * evt_
Definition: Generator.h:61
double Z_lmin
Definition: Generator.h:64
int verbose
Definition: Generator.h:59
std::vector< int > pdgFilter
Definition: Generator.h:65
bool fEtaCuts
Definition: Generator.h:45
double Z_hector
Definition: Generator.h:64
bool fPCuts
Definition: Generator.h:43

◆ ~Generator()

Generator::~Generator ( )
virtual

Definition at line 110 of file Generator.cc.

References fLumiFilter.

110 { delete fLumiFilter; }
LumiMonitorFilter * fLumiFilter
Definition: Generator.h:60

Member Function Documentation

◆ eventWeight()

virtual const double Generator::eventWeight ( ) const
inlinevirtual

Definition at line 32 of file Generator.h.

References weight_.

Referenced by RunManagerMTWorker::produce().

32 { return weight_; }
double weight_
Definition: Generator.h:63

◆ genEvent()

virtual const HepMC::GenEvent* Generator::genEvent ( ) const
inlinevirtual

Definition at line 30 of file Generator.h.

References evt_.

Referenced by RunManagerMTWorker::produce().

30 { return evt_; }
HepMC::GenEvent * evt_
Definition: Generator.h:61

◆ genVertex()

virtual const math::XYZTLorentzVector* Generator::genVertex ( ) const
inlinevirtual

Definition at line 31 of file Generator.h.

References vtx_.

Referenced by RunManagerMTWorker::produce().

31 { return vtx_; }
math::XYZTLorentzVector * vtx_
Definition: Generator.h:62

◆ HepMC2G4()

void Generator::HepMC2G4 ( const HepMC::GenEvent g,
G4Event *  e 
)

Definition at line 112 of file Generator.cc.

References funct::abs(), ALCARECOTkAlJpsiMuMu_cff::charge, class-composition::children, fEtaCuts, fFiductialCuts, fLumiFilter, fPCuts, fPDGFilter, fPhiCuts, fPtransCut, GenParticle::GenParticle, isExotic(), isExoticNonDetectable(), LumiMonitorFilter::isGoodForLumiMonitor(), IsInTheFilterList(), LogDebug, visualization-live-secondInstance_cfg::m, maxZCentralCMS, AlCaHLTBitMon_ParallelJobs::p, particleAssignDaughters(), dileptonTrigSettings_cff::pdg, pdgFilterSel, phi, multPhiCorr_741_25nsDY_cfi::px, multPhiCorr_741_25nsDY_cfi::py, setGenId(), mathSSE::sqrt(), contentValuesCheck::ss, mps_update::status, RandomServiceHelper::t1, theDecRCut2, theMaxPCut, theMaxPhiCut, theMinPtCut2, verbose, vtx_, weight_, testProducerWithPsetDescEmpty_cfi::x1, testProducerWithPsetDescEmpty_cfi::x2, protons_cff::xi, geometryCSVtoXML::xx, testProducerWithPsetDescEmpty_cfi::y1, testProducerWithPsetDescEmpty_cfi::y2, geometryCSVtoXML::yy, testProducerWithPsetDescEmpty_cfi::z2, Z_hector, Z_lmax, and Z_lmin.

Referenced by RunManagerMTWorker::generateEvent().

112  {
113  HepMC::GenEvent *evt = new HepMC::GenEvent(*evt_orig);
114 
115  if (*(evt->vertices_begin()) == nullptr) {
116  std::stringstream ss;
117  ss << "SimG4Core/Generator: in event " << g4evt->GetEventID() << " Corrupted Event - GenEvent with no vertex \n";
118  throw SimG4Exception(ss.str());
119  }
120 
121  if (!evt->weights().empty()) {
122  weight_ = evt->weights()[0];
123  for (unsigned int iw = 1; iw < evt->weights().size(); ++iw) {
124  // terminate if the vector of weights contains a zero-weight
125  if (evt->weights()[iw] <= 0)
126  break;
127  weight_ *= evt->weights()[iw];
128  }
129  }
130 
131  if (vtx_ != nullptr) {
132  delete vtx_;
133  }
134  vtx_ = new math::XYZTLorentzVector((*(evt->vertices_begin()))->position().x(),
135  (*(evt->vertices_begin()))->position().y(),
136  (*(evt->vertices_begin()))->position().z(),
137  (*(evt->vertices_begin()))->position().t());
138 
139  edm::LogVerbatim("SimG4CoreGenerator") << "Generator: primary Vertex = (" << vtx_->x() << ", " << vtx_->y() << ", "
140  << vtx_->z() << ")";
141 
142  unsigned int ng4vtx = 0;
143  unsigned int ng4par = 0;
144 
145  for (HepMC::GenEvent::vertex_const_iterator vitr = evt->vertices_begin(); vitr != evt->vertices_end(); ++vitr) {
146  // loop for vertex, is it a real vertex?
147  // Set qvtx to true for any particles that should be propagated by GEANT,
148  // i.e., status 1 particles or status 2 particles that decay outside the
149  // beampipe.
150  G4bool qvtx = false;
151  HepMC::GenVertex::particle_iterator pitr;
152  for (pitr = (*vitr)->particles_begin(HepMC::children); pitr != (*vitr)->particles_end(HepMC::children); ++pitr) {
153  // For purposes of this function, the status is defined as follows:
154  // 1: particles are not decayed by generator
155  // 2: particles are decayed by generator but need to be propagated by GEANT
156  // 3: particles are decayed by generator and do not need to be propagated by GEANT
157  int status = (*pitr)->status();
158  int pdg = (*pitr)->pdg_id();
159  if (status > 3 && isExotic(pdg) && (!(isExoticNonDetectable(pdg)))) {
160  // In Pythia 8, there are many status codes besides 1, 2, 3.
161  // By setting the status to 2 for exotic particles, they will be
162  // checked: if its decay vertex is outside the beampipe, it will be
163  // propagated by GEANT. Some Standard Model particles, e.g., K0, cannot
164  // be propagated by GEANT, so do not change their status code.
165  status = 2;
166  }
167  if (status == 2 && abs(pdg) == 9900015) {
168  status = 3;
169  }
170 
171  // Particles which are not decayed by generator
172  if (status == 1) {
173  // filter out unwanted particles and vertices
175  continue;
176  }
177 
178  qvtx = true;
179  if (verbose > 2)
180  LogDebug("SimG4CoreGenerator") << "GenVertex barcode = " << (*vitr)->barcode()
181  << " selected for GenParticle barcode = " << (*pitr)->barcode();
182  break;
183  }
184  // The selection is made considering if the partcile with status = 2
185  // have the end_vertex with a radius greater than the radius of beampipe
186  // cylinder (no requirement on the Z of the vertex is applyed).
187  else if (status == 2) {
188  if ((*pitr)->end_vertex() != nullptr) {
189  double xx = (*pitr)->end_vertex()->position().x();
190  double yy = (*pitr)->end_vertex()->position().y();
191  double r_dd = xx * xx + yy * yy;
192  if (r_dd > theDecRCut2) {
193  qvtx = true;
194  if (verbose > 2)
195  LogDebug("SimG4CoreGenerator")
196  << "GenVertex barcode = " << (*vitr)->barcode()
197  << " selected for GenParticle barcode = " << (*pitr)->barcode() << " radius = " << std::sqrt(r_dd);
198  break;
199  }
200  } else {
201  // particles with status 2 without end_vertex are
202  // equivalent to stable
203  qvtx = true;
204  break;
205  }
206  }
207  }
208 
209  // if this vertex is inside fiductial volume inside the beam pipe
210  // and has no long-lived secondary the vertex is not saved
211  if (!qvtx) {
212  continue;
213  }
214 
215  double x1 = (*vitr)->position().x() * CLHEP::mm;
216  double y1 = (*vitr)->position().y() * CLHEP::mm;
217  double z1 = (*vitr)->position().z() * CLHEP::mm;
218  double t1 = (*vitr)->position().t() * CLHEP::mm / CLHEP::c_light;
219 
220  G4PrimaryVertex *g4vtx = new G4PrimaryVertex(x1, y1, z1, t1);
221 
222  for (pitr = (*vitr)->particles_begin(HepMC::children); pitr != (*vitr)->particles_end(HepMC::children); ++pitr) {
223  int status = (*pitr)->status();
224  int pdg = (*pitr)->pdg_id();
225  bool hasDecayVertex = (nullptr != (*pitr)->end_vertex());
226 
227  // Filter on allowed particle species if required
228  if (fPDGFilter) {
229  bool isInTheList = IsInTheFilterList(pdg);
230  if ((!pdgFilterSel && isInTheList) || (pdgFilterSel && !isInTheList)) {
231  if (0 < verbose)
232  edm::LogVerbatim("SimG4CoreGenerator")
233  << " Skiped GenParticle barcode= " << (*pitr)->barcode() << " PDGid= " << pdg << " status= " << status
234  << " isExotic: " << isExotic(pdg) << " isExoticNotDet: " << isExoticNonDetectable(pdg)
235  << " isInTheList: " << isInTheList << " hasDecayVertex: " << hasDecayVertex;
236  continue;
237  }
238  }
239 
240  if (0 < verbose) {
241  edm::LogVerbatim("SimG4CoreGenerator")
242  << "Generator: pdg= " << pdg << " status= " << status << " hasPreDefinedDecay: " << hasDecayVertex
243  << "\n isExotic: " << isExotic(pdg) << " isExoticNotDet: " << isExoticNonDetectable(pdg)
244  << " isInTheList: " << IsInTheFilterList(pdg) << "\n"
245  << " MaxZCentralCMS = " << maxZCentralCMS / CLHEP::m << " m; (x,y,z,t): (" << x1 << "," << y1 << "," << z1
246  << "," << t1 << ")";
247  }
248  if (status > 3 && isExotic(pdg) && (!(isExoticNonDetectable(pdg)))) {
249  status = hasDecayVertex ? 2 : 1;
250  }
251  if (status == 2 && abs(pdg) == 9900015) {
252  status = 3;
253  }
254 
255  // this particle has predefined decay but has no vertex
256  if (2 == status && !hasDecayVertex) {
257  edm::LogWarning("SimG4CoreGenerator: in event ")
258  << g4evt->GetEventID() << " a particle "
259  << " pdgid= " << pdg << " has status=2 but has no decay vertex, so will be fully tracked by Geant4";
260  status = 1;
261  }
262 
263  double x2 = x1;
264  double y2 = y1;
265  double z2 = z1;
266  double decay_length = 0.0;
267  if (2 == status) {
268  x2 = (*pitr)->end_vertex()->position().x();
269  y2 = (*pitr)->end_vertex()->position().y();
270  z2 = (*pitr)->end_vertex()->position().z();
271  decay_length = std::sqrt((x2 - x1) * (x2 - x1) + (y2 - y1) * (y2 - y1) + (z2 - z1) * (z2 - z1));
272  }
273 
274  bool toBeAdded = !fFiductialCuts;
275 
276  double px = (*pitr)->momentum().px();
277  double py = (*pitr)->momentum().py();
278  double pz = (*pitr)->momentum().pz();
279  double ptot = std::sqrt(px * px + py * py + pz * pz);
280  math::XYZTLorentzVector p(px, py, pz, (*pitr)->momentum().e());
281 
282  double ximpact = x1;
283  double yimpact = y1;
284  double zimpact = z1;
285 
286  // protection against numerical problems for extremely low momenta
287  // compute impact point at transition to Hector
288  const double minTan = 1.e-20;
289  if (std::abs(z1) < Z_hector && std::abs(pz) >= minTan * ptot) {
290  zimpact = (pz > 0.0) ? Z_hector : -Z_hector;
291  double del = (zimpact - z1) / pz;
292  ximpact += del * px;
293  yimpact += del * py;
294  }
295  double rimpact2 = ximpact * ximpact + yimpact * yimpact;
296 
297  if (verbose > 2)
298  LogDebug("SimG4CoreGenerator") << "Processing GenParticle barcode= " << (*pitr)->barcode() << " pdg= " << pdg
299  << " status= " << (*pitr)->status() << " st= " << status
300  << " rimpact(cm)= " << std::sqrt(rimpact2) / cm
301  << " zimpact(cm)= " << zimpact / cm << " ptot(GeV)= " << ptot
302  << " pz(GeV)= " << pz;
303 
304  // Particles of status 1 trasnported along the beam pipe
305  // HECTOR transport of protons are done in corresponding PPS producer
306  if (1 == status && std::abs(zimpact) >= Z_hector && rimpact2 <= theDecRCut2) {
307  // very forward n, nbar, gamma are allowed
308  toBeAdded = (2112 == std::abs(pdg) || 22 == pdg);
309  if (verbose > 1) {
310  edm::LogVerbatim("SimG4CoreGenerator")
311  << "GenParticle barcode = " << (*pitr)->barcode() << " very forward; to be added: " << toBeAdded;
312  }
313  } else {
314  // Standard case: particles not decayed by the generator and not forward
315  if (1 == status && (std::abs(zimpact) < Z_hector || rimpact2 > theDecRCut2)) {
316  // Ptot cut for all particles
317  if (fPCuts && (ptot < theMinPCut || ptot > theMaxPCut)) {
318  continue;
319  }
320  // phi cut is applied mainly for particle gun
321  if (fPhiCuts) {
322  double phi = p.phi();
323  if (phi < theMinPhiCut || phi > theMaxPhiCut) {
324  continue;
325  }
326  }
327  // eta cut is applied if position of the decay
328  // is within vacuum chamber and limited in Z
329  if (fEtaCuts) {
330  // eta cut
331  double xi = x1;
332  double yi = y1;
333  double zi = z1;
334 
335  // can be propagated along Z
336  if (std::abs(pz) >= minTan * ptot) {
337  if ((zi >= Z_lmax) & (pz < 0.0)) {
338  zi = Z_lmax;
339  } else if ((zi <= Z_lmin) & (pz > 0.0)) {
340  zi = Z_lmin;
341  } else {
342  if (pz > 0) {
343  zi = Z_lmax;
344  } else {
345  zi = Z_lmin;
346  }
347  }
348  double del = (zi - z1) / pz;
349  xi += del * px;
350  yi += del * py;
351  }
352  // check eta cut
353  if ((zi >= Z_lmin) & (zi <= Z_lmax) & (xi * xi + yi * yi < theDecRCut2)) {
354  continue;
355  }
356  }
357  if (fLumiFilter && !fLumiFilter->isGoodForLumiMonitor(*pitr)) {
358  continue;
359  }
360  toBeAdded = true;
361  if (verbose > 1)
362  edm::LogVerbatim("SimG4CoreGenerator")
363  << "GenParticle barcode = " << (*pitr)->barcode() << " passed case 1";
364 
365  // Decay chain outside the fiducial cylinder defined by theRDecLenCut
366  // are used for Geant4 tracking with predefined decay channel
367  // In the case of decay in vacuum particle is not tracked by Geant4
368  } else if (2 == status && x2 * x2 + y2 * y2 >= theDecRCut2 && std::abs(z2) < Z_hector) {
369  toBeAdded = true;
370  if (verbose > 1)
371  edm::LogVerbatim("SimG4CoreGenerator") << "GenParticle barcode = " << (*pitr)->barcode() << " passed case 2"
372  << " decay_length(cm)= " << decay_length / CLHEP::cm;
373  }
374  }
375  if (toBeAdded) {
376  G4PrimaryParticle *g4prim = new G4PrimaryParticle(pdg, px * GeV, py * GeV, pz * GeV);
377 
378  if (g4prim->GetG4code() != nullptr) {
379  g4prim->SetMass(g4prim->GetG4code()->GetPDGMass());
380  double charge = g4prim->GetG4code()->GetPDGCharge();
381 
382  // apply Pt cut
383  if (fPtransCut && 1 == status && 0.0 != charge && px * px + py * py < theMinPtCut2) {
384  delete g4prim;
385  continue;
386  }
387  g4prim->SetCharge(charge);
388  }
389 
390  // V.I. do not use SetWeight but the same code
391  // value of the code compute inside TrackWithHistory
392  // g4prim->SetWeight( 10000*(*vpitr)->barcode() ) ;
393  setGenId(g4prim, (*pitr)->barcode());
394 
395  if (2 == status) {
396  particleAssignDaughters(g4prim, (HepMC::GenParticle *)*pitr, decay_length);
397  }
398  if (verbose > 1)
399  g4prim->Print();
400 
401  ++ng4par;
402  g4vtx->SetPrimary(g4prim);
403  edm::LogVerbatim("SimG4CoreGenerator") << " " << ng4par << ". new Geant4 particle pdg= " << pdg
404  << " Ptot(GeV/c)= " << ptot << " Pt= " << std::sqrt(px * px + py * py)
405  << " status= " << status << "; dir= " << g4prim->GetMomentumDirection();
406  }
407  }
408 
409  if (verbose > 1)
410  g4vtx->Print();
411  g4evt->AddPrimaryVertex(g4vtx);
412  ++ng4vtx;
413  }
414 
415  // Add a protection for completely empty events (produced by LHCTransport):
416  // add a dummy vertex with no particle attached to it
417  if (ng4vtx == 0) {
418  G4PrimaryVertex *g4vtx = new G4PrimaryVertex(0.0, 0.0, 0.0, 0.0);
419  if (verbose > 1)
420  g4vtx->Print();
421 
422  g4evt->AddPrimaryVertex(g4vtx);
423  }
424 
425  edm::LogVerbatim("SimG4CoreGenerator") << "The list of Geant4 primaries includes " << ng4par << " particles in "
426  << ng4vtx << " vertex";
427  delete evt;
428 }
double theDecRCut2
Definition: Generator.h:55
Log< level::Info, true > LogVerbatim
bool fPhiCuts
Definition: Generator.h:46
bool isExoticNonDetectable(int pdgcode) const
Definition: Generator.cc:535
double theMinPtCut2
Definition: Generator.h:53
double maxZCentralCMS
Definition: Generator.h:58
double Z_lmax
Definition: Generator.h:64
double theMaxPhiCut
Definition: Generator.h:49
bool isExotic(int pdgcode) const
Definition: Generator.cc:527
bool fPtransCut
Definition: Generator.h:44
bool isGoodForLumiMonitor(const HepMC::GenParticle *) const
bool pdgFilterSel
Definition: Generator.h:66
void particleAssignDaughters(G4PrimaryParticle *p, HepMC::GenParticle *hp, double length)
Definition: Generator.cc:430
XYZTLorentzVectorD XYZTLorentzVector
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:29
LumiMonitorFilter * fLumiFilter
Definition: Generator.h:60
bool fFiductialCuts
Definition: Generator.h:47
T sqrt(T t)
Definition: SSEVec.h:19
bool fPDGFilter
Definition: Generator.h:67
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
double weight_
Definition: Generator.h:63
double theMaxPCut
Definition: Generator.h:54
math::XYZTLorentzVector * vtx_
Definition: Generator.h:62
double Z_lmin
Definition: Generator.h:64
int verbose
Definition: Generator.h:59
void setGenId(G4PrimaryParticle *p, int id) const
Definition: Generator.h:40
Log< level::Warning, false > LogWarning
bool IsInTheFilterList(int pdgcode) const
Definition: Generator.cc:542
bool fEtaCuts
Definition: Generator.h:45
double Z_hector
Definition: Generator.h:64
#define LogDebug(id)
bool fPCuts
Definition: Generator.h:43

◆ isExotic()

bool Generator::isExotic ( int  pdgcode) const
private

Definition at line 527 of file Generator.cc.

References funct::abs(), and EgammaValidation_cff::pdgid.

Referenced by HepMC2G4().

527  {
528  int pdgid = std::abs(pdgcode);
529  return ((pdgid >= 1000000 && pdgid < 4000000 && pdgid != 3000022) || // SUSY, R-hadron, and technicolor particles
530  pdgid == 17 || // 4th generation lepton
531  pdgid == 34 || // W-prime
532  pdgid == 37); // charged Higgs
533 }
Abs< T >::type abs(const T &t)
Definition: Abs.h:22

◆ isExoticNonDetectable()

bool Generator::isExoticNonDetectable ( int  pdgcode) const
private

Definition at line 535 of file Generator.cc.

References funct::abs(), ALCARECOTkAlJpsiMuMu_cff::charge, LHEGenericFilter_cfi::ParticleID, and EgammaValidation_cff::pdgid.

Referenced by HepMC2G4().

535  {
536  int pdgid = std::abs(pdgcode);
537  HepPDT::ParticleID pid(pdgcode);
538  int charge = pid.threeCharge();
539  return ((charge == 0) && (pdgid >= 1000000 && pdgid < 1000040)); // SUSY
540 }
Abs< T >::type abs(const T &t)
Definition: Abs.h:22

◆ IsInTheFilterList()

bool Generator::IsInTheFilterList ( int  pdgcode) const
private

Definition at line 542 of file Generator.cc.

References funct::abs(), pdgFilter, and EgammaValidation_cff::pdgid.

Referenced by HepMC2G4().

542  {
543  int pdgid = std::abs(pdgcode);
544  for (auto &pdg : pdgFilter) {
545  if (std::abs(pdg) == pdgid) {
546  return true;
547  }
548  }
549  return false;
550 }
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
std::vector< int > pdgFilter
Definition: Generator.h:65

◆ nonCentralEvent2G4()

void Generator::nonCentralEvent2G4 ( const HepMC::GenEvent g,
G4Event *  e 
)

Definition at line 552 of file Generator.cc.

References GenParticle::GenParticle, runTauDisplay::gp, mps_fire::i, ALPAKA_ACCELERATOR_NAMESPACE::vertexFinder::it, setGenId(), findQualityFiles::v, and verbose.

Referenced by RunManagerMTWorker::generateEvent().

552  {
553  int i = g4evt->GetNumberOfPrimaryVertex();
554  for (HepMC::GenEvent::particle_const_iterator it = evt->particles_begin(); it != evt->particles_end(); ++it) {
555  ++i;
556  HepMC::GenParticle *gp = (*it);
557 
558  // storing only particle with status == 1
559  if (gp->status() != 1)
560  continue;
561 
562  int pdg = gp->pdg_id();
563  G4PrimaryParticle *g4p = new G4PrimaryParticle(
564  pdg, gp->momentum().px() * CLHEP::GeV, gp->momentum().py() * CLHEP::GeV, gp->momentum().pz() * CLHEP::GeV);
565  if (g4p->GetG4code() != nullptr) {
566  g4p->SetMass(g4p->GetG4code()->GetPDGMass());
567  g4p->SetCharge(g4p->GetG4code()->GetPDGCharge());
568  }
569  setGenId(g4p, i);
570  G4PrimaryVertex *v = new G4PrimaryVertex(gp->production_vertex()->position().x() * CLHEP::mm,
571  gp->production_vertex()->position().y() * CLHEP::mm,
572  gp->production_vertex()->position().z() * CLHEP::mm,
573  gp->production_vertex()->position().t() * CLHEP::mm / CLHEP::c_light);
574  v->SetPrimary(g4p);
575  g4evt->AddPrimaryVertex(v);
576  if (verbose > 0)
577  v->Print();
578  } // end loop on HepMC particles
579 }
int verbose
Definition: Generator.h:59
void setGenId(G4PrimaryParticle *p, int id) const
Definition: Generator.h:40

◆ particleAssignDaughters()

void Generator::particleAssignDaughters ( G4PrimaryParticle *  p,
HepMC::GenParticle *  hp,
double  length 
)
private

Definition at line 430 of file Generator.cc.

References funct::abs(), class-composition::children, createTree::dd, LogDebug, AlCaHLTBitMon_ParallelJobs::p, setGenId(), mathSSE::sqrt(), mps_update::status, verbose, testProducerWithPsetDescEmpty_cfi::x1, testProducerWithPsetDescEmpty_cfi::x2, testProducerWithPsetDescEmpty_cfi::y1, testProducerWithPsetDescEmpty_cfi::y2, and testProducerWithPsetDescEmpty_cfi::z2.

Referenced by HepMC2G4().

430  {
431  if (verbose > 1) {
432  LogDebug("SimG4CoreGenerator") << "Special case of long decay length \n"
433  << "Assign daughters with to mother with decaylength=" << decaylength / cm << " cm";
434  }
435  math::XYZTLorentzVector p(vp->momentum().px(), vp->momentum().py(), vp->momentum().pz(), vp->momentum().e());
436 
437  // defined preassigned decay time
438  double proper_time = decaylength / (p.Beta() * p.Gamma() * c_light);
439  g4p->SetProperTime(proper_time);
440 
441  if (verbose > 2) {
442  LogDebug("SimG4CoreGenerator") << " px= " << p.px() << " py= " << p.py() << " pz= " << p.pz() << " e= " << p.e()
443  << " beta= " << p.Beta() << " gamma= " << p.Gamma()
444  << " Proper time= " << proper_time / ns << " ns";
445  }
446 
447  // the particle will decay after the same length if it
448  // has not interacted before
449  double x1 = vp->end_vertex()->position().x();
450  double y1 = vp->end_vertex()->position().y();
451  double z1 = vp->end_vertex()->position().z();
452 
453  for (HepMC::GenVertex::particle_iterator vpdec = vp->end_vertex()->particles_begin(HepMC::children);
454  vpdec != vp->end_vertex()->particles_end(HepMC::children);
455  ++vpdec) {
456  // transform decay products such that in the rest frame of mother
458  (*vpdec)->momentum().px(), (*vpdec)->momentum().py(), (*vpdec)->momentum().pz(), (*vpdec)->momentum().e());
459 
460  // children should only be taken into account once
461  G4PrimaryParticle *g4daught =
462  new G4PrimaryParticle((*vpdec)->pdg_id(), pdec.x() * CLHEP::GeV, pdec.y() * CLHEP::GeV, pdec.z() * CLHEP::GeV);
463 
464  if (g4daught->GetG4code() != nullptr) {
465  g4daught->SetMass(g4daught->GetG4code()->GetPDGMass());
466  g4daught->SetCharge(g4daught->GetG4code()->GetPDGCharge());
467  }
468 
469  // V.I. do not use SetWeight but the same code
470  // value of the code compute inside TrackWithHistory
471  setGenId(g4daught, (*vpdec)->barcode());
472 
473  int status = (*vpdec)->status();
474  if (verbose > 1)
475  LogDebug("SimG4CoreGenerator::::particleAssignDaughters")
476  << "Assigning a " << (*vpdec)->pdg_id() << " as daughter of a " << vp->pdg_id() << " status=" << status;
477 
478  if ((status == 2 || (status == 23 && std::abs(vp->pdg_id()) == 1000015) || (status > 50 && status < 100)) &&
479  (*vpdec)->end_vertex() != nullptr) {
480  double x2 = (*vpdec)->end_vertex()->position().x();
481  double y2 = (*vpdec)->end_vertex()->position().y();
482  double z2 = (*vpdec)->end_vertex()->position().z();
483  double dd = std::sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2) + (z1 - z2) * (z1 - z2));
484  particleAssignDaughters(g4daught, *vpdec, dd);
485  }
486  (*vpdec)->set_status(1000 + status);
487  g4p->SetDaughter(g4daught);
488 
489  if (verbose > 1)
490  g4daught->Print();
491  }
492 }
string dd
Definition: createTree.py:154
void particleAssignDaughters(G4PrimaryParticle *p, HepMC::GenParticle *hp, double length)
Definition: Generator.cc:430
XYZTLorentzVectorD XYZTLorentzVector
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:29
T sqrt(T t)
Definition: SSEVec.h:19
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
int verbose
Definition: Generator.h:59
void setGenId(G4PrimaryParticle *p, int id) const
Definition: Generator.h:40
#define LogDebug(id)

◆ particlePassesPrimaryCuts()

bool Generator::particlePassesPrimaryCuts ( const G4ThreeVector &  p) const
private

Definition at line 495 of file Generator.cc.

References MillePedeFileConverter_cfg::e, PVValHelper::eta, fEtaCuts, RemoveAddSevLevel::flag, fPCuts, fPhiCuts, LogDebug, AlCaHLTBitMon_ParallelJobs::p, phi, theMaxEtaCut, theMaxPCut, theMaxPhiCut, and verbose.

495  {
496  bool flag = true;
497  double ptot = p.mag();
498  if (fPCuts && (ptot < theMinPCut * CLHEP::GeV || ptot > theMaxPCut * CLHEP::GeV)) {
499  flag = false;
500  }
501  if (fEtaCuts && flag) {
502  double pz = p.z();
503  if (ptot < pz + 1.e-10) {
504  flag = false;
505 
506  } else {
507  double eta = 0.5 * G4Log((ptot + pz) / (ptot - pz));
508  if (eta < theMinEtaCut || eta > theMaxEtaCut) {
509  flag = false;
510  }
511  }
512  }
513  if (fPhiCuts && flag) {
514  double phi = p.phi();
515  if (phi < theMinPhiCut || phi > theMaxPhiCut) {
516  flag = false;
517  }
518  }
519 
520  if (verbose > 2)
521  LogDebug("SimG4CoreGenerator") << "Generator ptot(GeV)= " << ptot / GeV << " eta= " << p.eta()
522  << " phi= " << p.phi() << " Flag= " << flag;
523 
524  return flag;
525 }
bool fPhiCuts
Definition: Generator.h:46
double theMaxPhiCut
Definition: Generator.h:49
double theMaxEtaCut
Definition: Generator.h:51
double theMaxPCut
Definition: Generator.h:54
int verbose
Definition: Generator.h:59
bool fEtaCuts
Definition: Generator.h:45
#define LogDebug(id)
bool fPCuts
Definition: Generator.h:43

◆ setGenEvent()

void Generator::setGenEvent ( const HepMC::GenEvent inpevt)
inline

Definition at line 24 of file Generator.h.

References evt_.

Referenced by RunManagerMTWorker::generateEvent().

24  {
25  evt_ = (HepMC::GenEvent *)inpevt;
26  return;
27  }
HepMC::GenEvent * evt_
Definition: Generator.h:61

◆ setGenId()

void Generator::setGenId ( G4PrimaryParticle *  p,
int  id 
) const
inlineprivate

Member Data Documentation

◆ evt_

HepMC::GenEvent* Generator::evt_
private

Definition at line 61 of file Generator.h.

Referenced by genEvent(), and setGenEvent().

◆ fEtaCuts

bool Generator::fEtaCuts
private

Definition at line 45 of file Generator.h.

Referenced by Generator(), HepMC2G4(), and particlePassesPrimaryCuts().

◆ fFiductialCuts

bool Generator::fFiductialCuts
private

Definition at line 47 of file Generator.h.

Referenced by Generator(), and HepMC2G4().

◆ fLumiFilter

LumiMonitorFilter* Generator::fLumiFilter
private

Definition at line 60 of file Generator.h.

Referenced by Generator(), HepMC2G4(), and ~Generator().

◆ fPCuts

bool Generator::fPCuts
private

Definition at line 43 of file Generator.h.

Referenced by Generator(), HepMC2G4(), and particlePassesPrimaryCuts().

◆ fPDGFilter

bool Generator::fPDGFilter
private

Definition at line 67 of file Generator.h.

Referenced by Generator(), and HepMC2G4().

◆ fPhiCuts

bool Generator::fPhiCuts
private

Definition at line 46 of file Generator.h.

Referenced by Generator(), HepMC2G4(), and particlePassesPrimaryCuts().

◆ fPtransCut

bool Generator::fPtransCut
private

Definition at line 44 of file Generator.h.

Referenced by Generator(), and HepMC2G4().

◆ maxZCentralCMS

double Generator::maxZCentralCMS
private

Definition at line 58 of file Generator.h.

Referenced by Generator(), and HepMC2G4().

◆ pdgFilter

std::vector<int> Generator::pdgFilter
private

Definition at line 65 of file Generator.h.

Referenced by Generator(), and IsInTheFilterList().

◆ pdgFilterSel

bool Generator::pdgFilterSel
private

Definition at line 66 of file Generator.h.

Referenced by Generator(), and HepMC2G4().

◆ theDecLenCut

double Generator::theDecLenCut
private

Definition at line 57 of file Generator.h.

Referenced by Generator().

◆ theDecRCut2

double Generator::theDecRCut2
private

Definition at line 55 of file Generator.h.

Referenced by Generator(), and HepMC2G4().

◆ theEtaCutForHector

double Generator::theEtaCutForHector
private

Definition at line 56 of file Generator.h.

Referenced by Generator().

◆ theMaxEtaCut

double Generator::theMaxEtaCut
private

Definition at line 51 of file Generator.h.

Referenced by Generator(), and particlePassesPrimaryCuts().

◆ theMaxPCut

double Generator::theMaxPCut
private

Definition at line 54 of file Generator.h.

Referenced by Generator(), HepMC2G4(), and particlePassesPrimaryCuts().

◆ theMaxPhiCut

double Generator::theMaxPhiCut
private

Definition at line 49 of file Generator.h.

Referenced by Generator(), HepMC2G4(), and particlePassesPrimaryCuts().

◆ theMinEtaCut

double Generator::theMinEtaCut
private

Definition at line 50 of file Generator.h.

Referenced by Generator().

◆ theMinPCut

double Generator::theMinPCut
private

Definition at line 52 of file Generator.h.

Referenced by Generator().

◆ theMinPhiCut

double Generator::theMinPhiCut
private

Definition at line 48 of file Generator.h.

Referenced by Generator().

◆ theMinPtCut2

double Generator::theMinPtCut2
private

Definition at line 53 of file Generator.h.

Referenced by Generator(), and HepMC2G4().

◆ verbose

int Generator::verbose
private

◆ vtx_

math::XYZTLorentzVector* Generator::vtx_
private

Definition at line 62 of file Generator.h.

Referenced by genVertex(), and HepMC2G4().

◆ weight_

double Generator::weight_
private

Definition at line 63 of file Generator.h.

Referenced by eventWeight(), and HepMC2G4().

◆ Z_hector

double Generator::Z_hector
private

Definition at line 64 of file Generator.h.

Referenced by Generator(), and HepMC2G4().

◆ Z_lmax

double Generator::Z_lmax
private

Definition at line 64 of file Generator.h.

Referenced by Generator(), and HepMC2G4().

◆ Z_lmin

double Generator::Z_lmin
private

Definition at line 64 of file Generator.h.

Referenced by Generator(), and HepMC2G4().