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

◆ ~Generator()

Generator::~Generator ( )
virtual

Definition at line 111 of file Generator.cc.

References fLumiFilter.

111 { delete fLumiFilter; }
LumiMonitorFilter * fLumiFilter
Definition: Generator.h:61

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:64

◆ 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:62

◆ 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:63

◆ HepMC2G4()

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

Definition at line 113 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().

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

◆ isExotic()

bool Generator::isExotic ( int  pdgcode) const
private

Definition at line 531 of file Generator.cc.

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

Referenced by HepMC2G4().

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

◆ isExoticNonDetectable()

bool Generator::isExoticNonDetectable ( int  pdgcode) const
private

Definition at line 539 of file Generator.cc.

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

Referenced by HepMC2G4().

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

◆ IsInTheFilterList()

bool Generator::IsInTheFilterList ( int  pdgcode) const
private

Definition at line 546 of file Generator.cc.

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

Referenced by HepMC2G4().

546  {
547  int pdgid = std::abs(pdgcode);
548  for (auto &pdg : pdgFilter) {
549  if (std::abs(pdg) == pdgid) {
550  return true;
551  }
552  }
553  return false;
554 }
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
std::vector< int > pdgFilter
Definition: Generator.h:66

◆ nonCentralEvent2G4()

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

Definition at line 556 of file Generator.cc.

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

Referenced by RunManagerMTWorker::generateEvent().

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

References funct::abs(), class-composition::children, createTree::dd, fFixG4Primary, 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().

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

499  {
500  bool flag = true;
501  double ptot = p.mag();
502  if (fPCuts && (ptot < theMinPCut * CLHEP::GeV || ptot > theMaxPCut * CLHEP::GeV)) {
503  flag = false;
504  }
505  if (fEtaCuts && flag) {
506  double pz = p.z();
507  if (ptot < pz + 1.e-10) {
508  flag = false;
509 
510  } else {
511  double eta = 0.5 * G4Log((ptot + pz) / (ptot - pz));
512  if (eta < theMinEtaCut || eta > theMaxEtaCut) {
513  flag = false;
514  }
515  }
516  }
517  if (fPhiCuts && flag) {
518  double phi = p.phi();
519  if (phi < theMinPhiCut || phi > theMaxPhiCut) {
520  flag = false;
521  }
522  }
523 
524  if (verbose > 2)
525  LogDebug("SimG4CoreGenerator") << "Generator ptot(GeV)= " << ptot / GeV << " eta= " << p.eta()
526  << " phi= " << p.phi() << " Flag= " << flag;
527 
528  return flag;
529 }
bool fPhiCuts
Definition: Generator.h:46
double theMaxPhiCut
Definition: Generator.h:50
double theMaxEtaCut
Definition: Generator.h:52
double theMaxPCut
Definition: Generator.h:55
int verbose
Definition: Generator.h:60
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:62

◆ setGenId()

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

Member Data Documentation

◆ evt_

HepMC::GenEvent* Generator::evt_
private

Definition at line 62 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 48 of file Generator.h.

Referenced by Generator(), and HepMC2G4().

◆ fFixG4Primary

bool Generator::fFixG4Primary
private

Definition at line 47 of file Generator.h.

Referenced by particleAssignDaughters().

◆ fLumiFilter

LumiMonitorFilter* Generator::fLumiFilter
private

Definition at line 61 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 68 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 59 of file Generator.h.

Referenced by Generator(), and HepMC2G4().

◆ pdgFilter

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

Definition at line 66 of file Generator.h.

Referenced by Generator(), and IsInTheFilterList().

◆ pdgFilterSel

bool Generator::pdgFilterSel
private

Definition at line 67 of file Generator.h.

Referenced by Generator(), and HepMC2G4().

◆ theDecLenCut

double Generator::theDecLenCut
private

Definition at line 58 of file Generator.h.

Referenced by Generator().

◆ theDecRCut2

double Generator::theDecRCut2
private

Definition at line 56 of file Generator.h.

Referenced by Generator(), and HepMC2G4().

◆ theEtaCutForHector

double Generator::theEtaCutForHector
private

Definition at line 57 of file Generator.h.

Referenced by Generator().

◆ theMaxEtaCut

double Generator::theMaxEtaCut
private

Definition at line 52 of file Generator.h.

Referenced by Generator(), and particlePassesPrimaryCuts().

◆ theMaxPCut

double Generator::theMaxPCut
private

Definition at line 55 of file Generator.h.

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

◆ theMaxPhiCut

double Generator::theMaxPhiCut
private

Definition at line 50 of file Generator.h.

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

◆ theMinEtaCut

double Generator::theMinEtaCut
private

Definition at line 51 of file Generator.h.

Referenced by Generator().

◆ theMinPCut

double Generator::theMinPCut
private

Definition at line 53 of file Generator.h.

Referenced by Generator().

◆ theMinPhiCut

double Generator::theMinPhiCut
private

Definition at line 49 of file Generator.h.

Referenced by Generator().

◆ theMinPtCut2

double Generator::theMinPtCut2
private

Definition at line 54 of file Generator.h.

Referenced by Generator(), and HepMC2G4().

◆ verbose

int Generator::verbose
private

◆ vtx_

math::XYZTLorentzVector* Generator::vtx_
private

Definition at line 63 of file Generator.h.

Referenced by genVertex(), and HepMC2G4().

◆ weight_

double Generator::weight_
private

Definition at line 64 of file Generator.h.

Referenced by eventWeight(), and HepMC2G4().

◆ Z_hector

double Generator::Z_hector
private

Definition at line 65 of file Generator.h.

Referenced by Generator(), and HepMC2G4().

◆ Z_lmax

double Generator::Z_lmax
private

Definition at line 65 of file Generator.h.

Referenced by Generator(), and HepMC2G4().

◆ Z_lmin

double Generator::Z_lmin
private

Definition at line 65 of file Generator.h.

Referenced by Generator(), and HepMC2G4().