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::GenEvent * genEvent () const
 
virtual const math::XYZTLorentzVectorgenVertex () const
 
void HepMC2G4 (const HepMC::GenEvent *g, G4Event *e)
 
void nonBeamEvent2G4 (const HepMC::GenEvent *g, G4Event *e)
 
void setGenEvent (const HepMC::GenEvent *inpevt)
 
virtual ~Generator ()
 

Private Member Functions

bool isExotic (HepMC::GenParticle *p) 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::GenEvent * evt_
 
bool fEtaCuts
 
LumiMonitorFilterfLumiFilter
 
bool fPCuts
 
bool fPDGFilter
 
bool fPhiCuts
 
bool fPtransCut
 
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 ( const edm::ParameterSet p)

Definition at line 22 of file Generator.cc.

References LumiMonitorFilter::Describe(), edm::ParameterSet::exists(), JetChargeProducer_cfi::exp, fEtaCuts, fLumiFilter, fPCuts, fPDGFilter, fPhiCuts, fPtransCut, edm::ParameterSet::getParameter(), cuy::ii, csvLumiCalc::lumi, pdgFilter, pdgFilterSel, theDecLenCut, theDecRCut2, theEtaCutForHector, theMaxEtaCut, theMinEtaCut, theMinPCut, theMinPtCut2, Z_hector, Z_lmax, and Z_lmin.

22  :
23  fPCuts(p.getParameter<bool>("ApplyPCuts")),
24  fPtransCut(p.getParameter<bool>("ApplyPtransCut")),
25  fEtaCuts(p.getParameter<bool>("ApplyEtaCuts")),
26  fPhiCuts(p.getParameter<bool>("ApplyPhiCuts")),
27  theMinPhiCut(p.getParameter<double>("MinPhiCut")), // in radians (CMS standard)
28  theMaxPhiCut(p.getParameter<double>("MaxPhiCut")),
29  theMinEtaCut(p.getParameter<double>("MinEtaCut")),
30  theMaxEtaCut(p.getParameter<double>("MaxEtaCut")),
31  theMinPCut(p.getParameter<double>("MinPCut")), // in GeV (CMS standard)
32  theMaxPCut(p.getParameter<double>("MaxPCut")),
33  theEtaCutForHector(p.getParameter<double>("EtaCutForHector")),
34  verbose(p.getUntrackedParameter<int>("Verbosity",0)),
35  fLumiFilter(nullptr),
36  evt_(nullptr),
37  vtx_(nullptr),
38  weight_(0),
39  Z_lmin(0),
40  Z_lmax(0),
41  Z_hector(0),
42  pdgFilterSel(false),
43  fPDGFilter(false)
44 {
45  bool lumi = p.getParameter<bool>("ApplyLumiMonitorCuts");
46  if(lumi) { fLumiFilter = new LumiMonitorFilter(); }
47 
48  double theRDecLenCut = p.getParameter<double>("RDecLenCut")*cm;
49  theDecRCut2 = theRDecLenCut*theRDecLenCut;
50 
52 
53  double theDecLenCut = p.getParameter<double>("LDecLenCut")*cm;
54 
55  pdgFilter.resize(0);
56  if ( p.exists("PDGselection") ) {
57  pdgFilterSel = (p.getParameter<edm::ParameterSet>("PDGselection")).
58  getParameter<bool>("PDGfilterSel");
59  pdgFilter = (p.getParameter<edm::ParameterSet>("PDGselection")).
60  getParameter<std::vector< int > >("PDGfilter");
61  if(0 < pdgFilter.size()) {
62  fPDGFilter = true;
63  for ( unsigned int ii = 0; ii < pdgFilter.size(); ++ii) {
64  if (pdgFilterSel) {
65  edm::LogWarning("SimG4CoreGenerator")
66  << " *** Selecting only PDG ID = " << pdgFilter[ii];
67  } else {
68  edm::LogWarning("SimG4CoreGenerator")
69  << " *** Filtering out PDG ID = " << pdgFilter[ii];
70  }
71  }
72  }
73  }
74 
75  if(fEtaCuts){
76  Z_lmax = theRDecLenCut*((1 - exp(-2*theMaxEtaCut) )/( 2*exp(-theMaxEtaCut)));
77  Z_lmin = theRDecLenCut*((1 - exp(-2*theMinEtaCut) )/( 2*exp(-theMinEtaCut)));
78  }
79 
80  Z_hector = theRDecLenCut*((1 - exp(-2*theEtaCutForHector))
81  / ( 2*exp(-theEtaCutForHector) ) );
82 
83  edm::LogInfo("SimG4CoreGenerator")
84  << " Rdecaycut= " << theRDecLenCut/cm
85  << " cm; Zdecaycut= " << theDecLenCut/cm
86  << "Z_min= " << Z_lmin/cm << " cm; Z_max= " << Z_lmax
87  << " cm; Z_hector = " << Z_hector << " cm\n"
88  << "ApplyPCuts: " << fPCuts << " ApplyPtransCut: " << fPtransCut
89  << " ApplyEtaCuts: " << fEtaCuts
90  << " ApplyPhiCuts: " << fPhiCuts
91  << " ApplyLumiMonitorCuts: " << lumi;
92  if(lumi) { fLumiFilter->Describe(); }
93 }
double theDecRCut2
Definition: Generator.h:54
bool fPhiCuts
Definition: Generator.h:46
T getParameter(std::string const &) const
T getUntrackedParameter(std::string const &, T const &) const
double theMinPtCut2
Definition: Generator.h:52
double theEtaCutForHector
Definition: Generator.h:55
double Z_lmax
Definition: Generator.h:62
double theMaxPhiCut
Definition: Generator.h:48
bool fPtransCut
Definition: Generator.h:44
bool exists(std::string const &parameterName) const
checks if a parameter exists
double theDecLenCut
Definition: Generator.h:56
double theMinEtaCut
Definition: Generator.h:49
bool pdgFilterSel
Definition: Generator.h:64
double theMinPCut
Definition: Generator.h:51
LumiMonitorFilter * fLumiFilter
Definition: Generator.h:58
bool fPDGFilter
Definition: Generator.h:65
double theMinPhiCut
Definition: Generator.h:47
double weight_
Definition: Generator.h:61
double theMaxEtaCut
Definition: Generator.h:50
ii
Definition: cuy.py:588
double theMaxPCut
Definition: Generator.h:53
void Describe() const
math::XYZTLorentzVector * vtx_
Definition: Generator.h:60
HepMC::GenEvent * evt_
Definition: Generator.h:59
double Z_lmin
Definition: Generator.h:62
int verbose
Definition: Generator.h:57
std::vector< int > pdgFilter
Definition: Generator.h:63
bool fEtaCuts
Definition: Generator.h:45
double Z_hector
Definition: Generator.h:62
bool fPCuts
Definition: Generator.h:43
Generator::~Generator ( )
virtual

Definition at line 95 of file Generator.cc.

References fLumiFilter.

96 {
97  delete fLumiFilter;
98 }
LumiMonitorFilter * fLumiFilter
Definition: Generator.h:58

Member Function Documentation

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

Definition at line 29 of file Generator.h.

References evt_.

Referenced by RunManagerMTWorker::produce(), and RunManager::produce().

29 { return evt_; }
HepMC::GenEvent * evt_
Definition: Generator.h:59
virtual const math::XYZTLorentzVector* Generator::genVertex ( ) const
inlinevirtual

Definition at line 30 of file Generator.h.

References vtx_.

Referenced by RunManagerMTWorker::produce(), and RunManager::produce().

30 { return vtx_; }
math::XYZTLorentzVector * vtx_
Definition: Generator.h:60
void Generator::HepMC2G4 ( const HepMC::GenEvent *  g,
G4Event *  e 
)

Definition at line 100 of file Generator.cc.

References ALCARECOTkAlJpsiMuMu_cff::charge, class-composition::children, fEtaCuts, spr::find(), fLumiFilter, fPCuts, fPDGFilter, fPhiCuts, fPtransCut, GenParticle::GenParticle, GeV, isExotic(), LumiMonitorFilter::isGoodForLumiMonitor(), LogDebug, nullptr, AlCaHLTBitMon_ParallelJobs::p, particleAssignDaughters(), HiggsValidation_cfi::pdg_id, pdgFilter, pdgFilterSel, phi, setGenId(), mathSSE::sqrt(), mps_update::status, theDecRCut2, theMaxPCut, theMaxPhiCut, theMinPtCut2, vtx_, weight_, hybridSuperClusters_cfi::xi, geometryCSVtoXML::xx, geometryCSVtoXML::yy, Z_hector, Z_lmax, and Z_lmin.

Referenced by RunManagerMTWorker::generateEvent(), RunManager::generateEvent(), and setGenEvent().

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

Definition at line 514 of file Generator.cc.

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

Referenced by eventWeight(), and HepMC2G4().

515 {
516  int pdgid = abs(p->pdg_id());
517  if ((pdgid >= 1000000 && pdgid < 4000000) || // SUSY, R-hadron, and technicolor particles
518  pdgid == 17 || // 4th generation lepton
519  pdgid == 34 || // W-prime
520  pdgid == 37) // charged Higgs
521  {
522  return true;
523  }
524 
525  return false;
526 }
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
void Generator::nonBeamEvent2G4 ( const HepMC::GenEvent *  g,
G4Event *  e 
)

Definition at line 529 of file Generator.cc.

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

Referenced by RunManagerMTWorker::generateEvent(), RunManager::generateEvent(), and setGenEvent().

530 {
531  int i = 0;
532  for(HepMC::GenEvent::particle_const_iterator it = evt->particles_begin();
533  it != evt->particles_end(); ++it ) {
534  ++i;
535  HepMC::GenParticle * gp = (*it);
536  int g_status = gp->status();
537  // storing only particle with status == 1
538  if (g_status == 1) {
539  int g_id = gp->pdg_id();
540  G4PrimaryParticle * g4p =
541  new G4PrimaryParticle(g_id,gp->momentum().px()*GeV,
542  gp->momentum().py()*GeV,
543  gp->momentum().pz()*GeV);
544  if (g4p->GetG4code() != 0) {
545  g4p->SetMass(g4p->GetG4code()->GetPDGMass());
546  g4p->SetCharge(g4p->GetG4code()->GetPDGCharge());
547  }
548  setGenId(g4p,i);
549  if (particlePassesPrimaryCuts(g4p->GetMomentum())) {
550  G4PrimaryVertex * v =
551  new G4PrimaryVertex(gp->production_vertex()->position().x()*mm,
552  gp->production_vertex()->position().y()*mm,
553  gp->production_vertex()->position().z()*mm,
554  gp->production_vertex()->position().t()*mm/c_light);
555  v->SetPrimary(g4p);
556  g4evt->AddPrimaryVertex(v);
557  if(verbose > 0) { v->Print(); }
558 
559  } else {
560  delete g4p;
561  }
562  }
563  } // end loop on HepMC particles
564 }
void setGenId(G4PrimaryParticle *p, int id) const
Definition: Generator.h:39
const double GeV
Definition: MathUtil.h:16
bool particlePassesPrimaryCuts(const G4ThreeVector &p) const
Definition: Generator.cc:481
void Generator::particleAssignDaughters ( G4PrimaryParticle *  p,
HepMC::GenParticle *  hp,
double  length 
)
private

Definition at line 404 of file Generator.cc.

References class-composition::children, createTree::dd, GeV, LogDebug, AlCaHLTBitMon_ParallelJobs::p, setGenId(), and mathSSE::sqrt().

Referenced by eventWeight(), and HepMC2G4().

407 {
408  if ( verbose > 1 ) {
409  LogDebug("SimG4CoreGenerator")
410  << "Special case of long decay length \n"
411  << "Assign daughters with to mother with decaylength="
412  << decaylength/cm << " cm";
413  }
414  math::XYZTLorentzVector p(vp->momentum().px(), vp->momentum().py(),
415  vp->momentum().pz(), vp->momentum().e());
416 
417  // defined preassigned decay time
418  double proper_time = decaylength/(p.Beta()*p.Gamma()*c_light);
419  g4p->SetProperTime(proper_time);
420 
421  if( verbose > 2 ) {
422  LogDebug("SimG4CoreGenerator") <<" px= "<< p.px()
423  <<" py= "<< p.py()
424  <<" pz= "<< p.pz()
425  <<" e= " << p.e()
426  <<" beta= "<< p.Beta()
427  <<" gamma= " << p.Gamma()
428  <<" Proper time= " <<proper_time/ns <<" ns";
429  }
430 
431  // the particle will decay after the same length if it
432  // has not interacted before
433  double x1 = vp->end_vertex()->position().x();
434  double y1 = vp->end_vertex()->position().y();
435  double z1 = vp->end_vertex()->position().z();
436 
437  for (HepMC::GenVertex::particle_iterator
438  vpdec= vp->end_vertex()->particles_begin(HepMC::children);
439  vpdec != vp->end_vertex()->particles_end(HepMC::children); ++vpdec) {
440 
441  //transform decay products such that in the rest frame of mother
442  math::XYZTLorentzVector pdec((*vpdec)->momentum().px(),
443  (*vpdec)->momentum().py(),
444  (*vpdec)->momentum().pz(),
445  (*vpdec)->momentum().e());
446 
447  // children should only be taken into account once
448  G4PrimaryParticle * g4daught=
449  new G4PrimaryParticle((*vpdec)->pdg_id(), pdec.x()*GeV,
450  pdec.y()*GeV, pdec.z()*GeV);
451 
452  if ( g4daught->GetG4code() != 0 )
453  {
454  g4daught->SetMass( g4daught->GetG4code()->GetPDGMass() ) ;
455  g4daught->SetCharge( g4daught->GetG4code()->GetPDGCharge() ) ;
456  }
457 
458  // V.I. do not use SetWeight but the same code
459  // value of the code compute inside TrackWithHistory
460  setGenId( g4daught, (*vpdec)->barcode() );
461 
462  if ( verbose > 2 ) LogDebug("SimG4CoreGenerator")
463  <<"Assigning a "<< (*vpdec)->pdg_id()
464  <<" as daughter of a " << vp->pdg_id();
465  if ( (*vpdec)->status() == 2 && (*vpdec)->end_vertex() != nullptr)
466  {
467  double x2 = (*vpdec)->end_vertex()->position().x();
468  double y2 = (*vpdec)->end_vertex()->position().y();
469  double z2 = (*vpdec)->end_vertex()->position().z();
470  double dd = std::sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)+(z1-z2)*(z1-z2));
471  particleAssignDaughters(g4daught,*vpdec,dd);
472  }
473  (*vpdec)->set_status(1000+(*vpdec)->status());
474  g4p->SetDaughter(g4daught);
475 
476  if ( verbose > 1 ) g4daught->Print();
477  }
478 }
#define LogDebug(id)
void setGenId(G4PrimaryParticle *p, int id) const
Definition: Generator.h:39
const double GeV
Definition: MathUtil.h:16
void particleAssignDaughters(G4PrimaryParticle *p, HepMC::GenParticle *hp, double length)
Definition: Generator.cc:404
XYZTLorentzVectorD XYZTLorentzVector
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:29
T sqrt(T t)
Definition: SSEVec.h:18
bool Generator::particlePassesPrimaryCuts ( const G4ThreeVector &  p) const
private

Definition at line 481 of file Generator.cc.

References MillePedeFileConverter_cfg::e, stringResolutionProvider_cfi::eta, fEtaCuts, RemoveAddSevLevel::flag, fPCuts, fPhiCuts, GeV, LogDebug, phi, theMaxEtaCut, theMaxPCut, and theMaxPhiCut.

Referenced by eventWeight(), and nonBeamEvent2G4().

482 {
483  bool flag = true;
484  double ptot = p.mag();
485  if (fPCuts && (ptot < theMinPCut*GeV || ptot > theMaxPCut*GeV)) {
486  flag = false;
487  }
488  if (fEtaCuts && flag) {
489  double pz = p.z();
490  if(ptot < pz + 1.e-10) {
491  flag = false;
492 
493  } else {
494  double eta = 0.5*G4Log((ptot + pz)/(ptot - pz));
495  if(eta < theMinEtaCut || eta > theMaxEtaCut) {
496  flag = false;
497  }
498  }
499  }
500  if (fPhiCuts && flag) {
501  double phi = p.phi();
502  if(phi < theMinPhiCut || phi > theMaxPhiCut) {
503  flag = false;
504  }
505  }
506 
507  if ( verbose > 2 ) LogDebug("SimG4CoreGenerator")
508  << "Generator ptot(GeV)= " << ptot/GeV << " eta= " << p.eta()
509  << " phi= " << p.phi() << " Flag= " << flag;
510 
511  return flag;
512 }
#define LogDebug(id)
bool fPhiCuts
Definition: Generator.h:46
const double GeV
Definition: MathUtil.h:16
double theMaxPhiCut
Definition: Generator.h:48
double theMaxEtaCut
Definition: Generator.h:50
double theMaxPCut
Definition: Generator.h:53
bool fEtaCuts
Definition: Generator.h:45
bool fPCuts
Definition: Generator.h:43
void Generator::setGenEvent ( const HepMC::GenEvent *  inpevt)
inline

Definition at line 25 of file Generator.h.

References MillePedeFileConverter_cfg::e, evt_, g, HepMC2G4(), nonBeamEvent2G4(), and reco::return().

Referenced by RunManagerMTWorker::generateEvent(), and RunManager::generateEvent().

26  { evt_ = (HepMC::GenEvent*)inpevt; return ; }
HepMC::GenEvent * evt_
Definition: Generator.h:59
return(e1-e2)*(e1-e2)+dp *dp
void Generator::setGenId ( G4PrimaryParticle *  p,
int  id 
) const
inlineprivate

Definition at line 39 of file Generator.h.

Referenced by HepMC2G4(), nonBeamEvent2G4(), and particleAssignDaughters().

Member Data Documentation

HepMC::GenEvent* Generator::evt_
private

Definition at line 59 of file Generator.h.

Referenced by genEvent(), and setGenEvent().

bool Generator::fEtaCuts
private

Definition at line 45 of file Generator.h.

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

LumiMonitorFilter* Generator::fLumiFilter
private

Definition at line 58 of file Generator.h.

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

bool Generator::fPCuts
private

Definition at line 43 of file Generator.h.

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

bool Generator::fPDGFilter
private

Definition at line 65 of file Generator.h.

Referenced by Generator(), and HepMC2G4().

bool Generator::fPhiCuts
private

Definition at line 46 of file Generator.h.

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

bool Generator::fPtransCut
private

Definition at line 44 of file Generator.h.

Referenced by Generator(), and HepMC2G4().

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

Definition at line 63 of file Generator.h.

Referenced by Generator(), and HepMC2G4().

bool Generator::pdgFilterSel
private

Definition at line 64 of file Generator.h.

Referenced by Generator(), and HepMC2G4().

double Generator::theDecLenCut
private

Definition at line 56 of file Generator.h.

Referenced by Generator().

double Generator::theDecRCut2
private

Definition at line 54 of file Generator.h.

Referenced by Generator(), and HepMC2G4().

double Generator::theEtaCutForHector
private

Definition at line 55 of file Generator.h.

Referenced by Generator().

double Generator::theMaxEtaCut
private

Definition at line 50 of file Generator.h.

Referenced by Generator(), and particlePassesPrimaryCuts().

double Generator::theMaxPCut
private

Definition at line 53 of file Generator.h.

Referenced by HepMC2G4(), and particlePassesPrimaryCuts().

double Generator::theMaxPhiCut
private

Definition at line 48 of file Generator.h.

Referenced by HepMC2G4(), and particlePassesPrimaryCuts().

double Generator::theMinEtaCut
private

Definition at line 49 of file Generator.h.

Referenced by Generator().

double Generator::theMinPCut
private

Definition at line 51 of file Generator.h.

Referenced by Generator().

double Generator::theMinPhiCut
private

Definition at line 47 of file Generator.h.

double Generator::theMinPtCut2
private

Definition at line 52 of file Generator.h.

Referenced by Generator(), and HepMC2G4().

int Generator::verbose
private

Definition at line 57 of file Generator.h.

math::XYZTLorentzVector* Generator::vtx_
private

Definition at line 60 of file Generator.h.

Referenced by genVertex(), and HepMC2G4().

double Generator::weight_
private

Definition at line 61 of file Generator.h.

Referenced by eventWeight(), and HepMC2G4().

double Generator::Z_hector
private

Definition at line 62 of file Generator.h.

Referenced by Generator(), and HepMC2G4().

double Generator::Z_lmax
private

Definition at line 62 of file Generator.h.

Referenced by Generator(), and HepMC2G4().

double Generator::Z_lmin
private

Definition at line 62 of file Generator.h.

Referenced by Generator(), and HepMC2G4().