1 | // HepMCInterface.cc is a part of the PYTHIA event generator. |
---|
2 | // Copyright (C) 2012 Mikhail Kirsanov, Torbjorn Sjostrand. |
---|
3 | // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details. |
---|
4 | // Please respect the MCnet Guidelines, see GUIDELINES for details. |
---|
5 | |
---|
6 | // Function definitions (not found in the header) for the I_Pythia8 class, |
---|
7 | // which converts a PYTHIA event record to the standard HepMC format. |
---|
8 | |
---|
9 | // Mikhail.Kirsanov@cern.ch |
---|
10 | |
---|
11 | // For GCC versions >= 4.6.0 can turn off shadow warnings |
---|
12 | #if ((__GNUC__*100)+__GNUC_MINOR__) >= 406 |
---|
13 | #pragma GCC diagnostic ignored "-Wshadow" |
---|
14 | #endif |
---|
15 | |
---|
16 | #include "HepMCInterface.h" |
---|
17 | #include "HepMC/GenEvent.h" |
---|
18 | |
---|
19 | // Switch shadow warnings back on |
---|
20 | #if ((__GNUC__*100)+__GNUC_MINOR__) >= 406 |
---|
21 | #pragma GCC diagnostic warning "-Wshadow" |
---|
22 | #endif |
---|
23 | |
---|
24 | namespace HepMC { |
---|
25 | |
---|
26 | //========================================================================== |
---|
27 | |
---|
28 | // Constructor and destructor. |
---|
29 | |
---|
30 | I_Pythia8::I_Pythia8(): |
---|
31 | m_trust_mothers_before_daughters(true), |
---|
32 | m_trust_both_mothers_and_daughters(false), |
---|
33 | m_print_inconsistency_errors(true), |
---|
34 | m_crash_on_problem(false), |
---|
35 | m_freepartonwarnings(true), |
---|
36 | m_convert_to_mev(false), |
---|
37 | m_mom_scale_factor(1.), |
---|
38 | m_internal_event_number(0), |
---|
39 | hasSetConvert(false), |
---|
40 | lengthScaleFactor(1.) {;} |
---|
41 | |
---|
42 | I_Pythia8::~I_Pythia8() {;} |
---|
43 | |
---|
44 | //-------------------------------------------------------------------------- |
---|
45 | |
---|
46 | // Main method for conversion from PYTHIA event to HepMC event. |
---|
47 | // Read one event from Pythia8 and fill GenEvent, |
---|
48 | // and return T/F = success/failure. |
---|
49 | |
---|
50 | bool I_Pythia8::fill_next_event( Pythia8::Event& pyev, GenEvent* evt, |
---|
51 | int ievnum ) { |
---|
52 | |
---|
53 | // 1. Error if no event passed. |
---|
54 | if ( !evt ) { |
---|
55 | std::cerr << "I_Pythia8::fill_next_event error - passed null event." |
---|
56 | << std::endl; |
---|
57 | return 0; |
---|
58 | } |
---|
59 | |
---|
60 | // Event number counter. |
---|
61 | if ( ievnum >= 0 ) { |
---|
62 | evt->set_event_number(ievnum); |
---|
63 | m_internal_event_number = ievnum; |
---|
64 | } else { |
---|
65 | evt->set_event_number(m_internal_event_number); |
---|
66 | m_internal_event_number++; |
---|
67 | } |
---|
68 | |
---|
69 | // Decide whether conversion from GeV to MeV is necessary. |
---|
70 | #ifdef HEPMC_HAS_UNITS |
---|
71 | m_mom_scale_factor = HepMC::Units::conversion_factor(HepMC::Units::GEV, |
---|
72 | evt->momentum_unit()); |
---|
73 | lengthScaleFactor = HepMC::Units::conversion_factor(HepMC::Units::MM, |
---|
74 | evt->length_unit()); |
---|
75 | #else |
---|
76 | if (!hasSetConvert) { |
---|
77 | m_convert_to_mev = true; |
---|
78 | m_mom_scale_factor = 1000.; |
---|
79 | } |
---|
80 | lengthScaleFactor = 1.; |
---|
81 | #endif |
---|
82 | |
---|
83 | // 2. Create a particle instance for each entry and fill a map, and |
---|
84 | // a vector which maps from the particle index to the GenParticle address. |
---|
85 | std::vector<GenParticle*> hepevt_particles( pyev.size() ); |
---|
86 | int i, istatus; |
---|
87 | for ( i = 1; i < pyev.size(); ++i ) { |
---|
88 | |
---|
89 | // The new HepMC status codes of February 2009 are now used by |
---|
90 | // default. If you want to recover the previous simpler behaviour |
---|
91 | // comment out the next line and uncomment the following one. |
---|
92 | istatus = pyev.statusHepMC(i); |
---|
93 | // istatus = (pyev[i].status() > 0) ? 1 : 2; |
---|
94 | |
---|
95 | // Fill the particle. |
---|
96 | hepevt_particles[i] = new GenParticle( |
---|
97 | FourVector( m_mom_scale_factor * pyev[i].p().px(), |
---|
98 | m_mom_scale_factor * pyev[i].p().py(), |
---|
99 | m_mom_scale_factor * pyev[i].p().pz(), |
---|
100 | m_mom_scale_factor * pyev[i].p().e() ), |
---|
101 | pyev[i].id(), istatus ); |
---|
102 | hepevt_particles[i]->suggest_barcode(i); |
---|
103 | hepevt_particles[i]->set_generated_mass( |
---|
104 | m_mom_scale_factor * pyev[i].m() ); |
---|
105 | |
---|
106 | // Colour flow uses index 1 and 2. |
---|
107 | int colType = pyev[i].colType(); |
---|
108 | if (colType == 1 || colType == 2) |
---|
109 | hepevt_particles[i]->set_flow(1, pyev[i].col()); |
---|
110 | if (colType == -1 || colType == 2) |
---|
111 | hepevt_particles[i]->set_flow(2, pyev[i].acol()); |
---|
112 | } |
---|
113 | |
---|
114 | // Here we assume that the first two particles in the list |
---|
115 | // are the incoming beam particles. |
---|
116 | evt->set_beam_particles( hepevt_particles[1], hepevt_particles[2] ); |
---|
117 | |
---|
118 | // 3 + 4. Loop over particles AGAIN, this time creating vertices. |
---|
119 | // We build EITHER the production or decay vertex for each entry in |
---|
120 | // hepevt, depending on the switch m_trust_mothers_before_daughters. |
---|
121 | // Note: the HEPEVT pointers are bi-directional, so sufficient to do one. |
---|
122 | for ( i = 1; i < pyev.size(); ++i ) { |
---|
123 | |
---|
124 | // 3. Build the production_vertex (if necessary). |
---|
125 | if ( m_trust_mothers_before_daughters || |
---|
126 | m_trust_both_mothers_and_daughters ) { |
---|
127 | GenParticle *p = hepevt_particles[i]; |
---|
128 | |
---|
129 | // 3a. Search to see if a production vertex already exists. |
---|
130 | std::vector<int> mothers = pyev.motherList(i); |
---|
131 | unsigned int imother = 0; |
---|
132 | int mother = -1; // note that in Pythia8 there is a particle number 0! |
---|
133 | if ( !mothers.empty() ) mother = mothers[imother]; |
---|
134 | GenVertex* prod_vtx = p->production_vertex(); |
---|
135 | while ( !prod_vtx && mother > 0 ) { |
---|
136 | prod_vtx = hepevt_particles[mother]->end_vertex(); |
---|
137 | if ( prod_vtx ) prod_vtx->add_particle_out( p ); |
---|
138 | imother++; |
---|
139 | if ( imother < mothers.size() ) mother = mothers[imother]; |
---|
140 | else mother = -1; |
---|
141 | } |
---|
142 | |
---|
143 | // 3b. If no suitable production vertex exists - and the particle has |
---|
144 | // at least one mother or position information to store - make one. |
---|
145 | FourVector prod_pos( lengthScaleFactor * pyev[i].xProd(), |
---|
146 | lengthScaleFactor * pyev[i].yProd(), |
---|
147 | lengthScaleFactor * pyev[i].zProd(), |
---|
148 | lengthScaleFactor * pyev[i].tProd() ); |
---|
149 | unsigned int nparents = mothers.size(); |
---|
150 | if ( !prod_vtx && ( nparents > 0 || prod_pos != FourVector() ) ) { |
---|
151 | prod_vtx = new GenVertex(); |
---|
152 | prod_vtx->add_particle_out( p ); |
---|
153 | evt->add_vertex( prod_vtx ); |
---|
154 | } |
---|
155 | |
---|
156 | // 3c. If prod_vtx doesn't already have position specified, fill it. |
---|
157 | if ( prod_vtx && prod_vtx->position() == FourVector() ) |
---|
158 | prod_vtx->set_position( prod_pos ); |
---|
159 | |
---|
160 | // 3d. loop over mothers to make sure their end_vertices are consistent. |
---|
161 | imother = 0; |
---|
162 | mother = -1; |
---|
163 | if ( !mothers.empty() ) mother = mothers[imother]; |
---|
164 | while ( prod_vtx && mother > 0 ) { |
---|
165 | |
---|
166 | // If end vertex of the mother isn't specified, do it now. |
---|
167 | if ( !hepevt_particles[mother]->end_vertex() ) { |
---|
168 | prod_vtx->add_particle_in( hepevt_particles[mother] ); |
---|
169 | |
---|
170 | // Problem scenario: the mother already has a decay vertex which |
---|
171 | // differs from the daughter's production vertex. This means there is |
---|
172 | // internal inconsistency in the HEPEVT event record. Print an error. |
---|
173 | // Note: we could provide a fix by joining the two vertices with a |
---|
174 | // dummy particle if the problem arises often. |
---|
175 | } else if (hepevt_particles[mother]->end_vertex() != prod_vtx ) { |
---|
176 | if ( m_print_inconsistency_errors ) std::cerr |
---|
177 | << "HepMC::I_Pythia8: inconsistent mother/daugher " |
---|
178 | << "information in Pythia8 event " << std::endl |
---|
179 | << "i= " << i << " mother = " << mother |
---|
180 | << "\n This warning can be turned off with the " |
---|
181 | << "I_Pythia8::print_inconsistency_errors switch." << std::endl; |
---|
182 | } |
---|
183 | |
---|
184 | // Variant with motherList. |
---|
185 | imother++; |
---|
186 | if ( imother < mothers.size() ) mother = mothers[imother]; |
---|
187 | else mother = -1; |
---|
188 | } |
---|
189 | |
---|
190 | // 4. Building from the mothers not implemented so far. |
---|
191 | } else { |
---|
192 | std::cerr << "trust_daughters_before_mothers not implemented" |
---|
193 | << std::endl; |
---|
194 | return 0; |
---|
195 | } |
---|
196 | } |
---|
197 | |
---|
198 | // 5. Check for particles which come from nowhere, i.e. are without |
---|
199 | // mothers or daughters. These need to be attached to a vertex, or else |
---|
200 | // they will never become part of the event. |
---|
201 | for ( i = 1; i < pyev.size(); ++i ) { |
---|
202 | if ( !hepevt_particles[i]->end_vertex() && |
---|
203 | !hepevt_particles[i]->production_vertex() ) { |
---|
204 | std::cerr << "hanging particle " << i << std::endl; |
---|
205 | GenVertex* prod_vtx = new GenVertex(); |
---|
206 | prod_vtx->add_particle_out( hepevt_particles[i] ); |
---|
207 | evt->add_vertex( prod_vtx ); |
---|
208 | } |
---|
209 | |
---|
210 | // Also check for free partons (= gluons and quarks; not diquarks?). |
---|
211 | if ( m_freepartonwarnings ) { |
---|
212 | if ( hepevt_particles[i]->pdg_id() == 21 && |
---|
213 | !hepevt_particles[i]->end_vertex() ) { |
---|
214 | std::cerr << "gluon without end vertex " << i << std::endl; |
---|
215 | if ( m_crash_on_problem ) exit(1); |
---|
216 | } |
---|
217 | if ( abs(hepevt_particles[i]->pdg_id()) <= 6 && |
---|
218 | !hepevt_particles[i]->end_vertex() ) { |
---|
219 | std::cerr << "quark without end vertex " << i << std::endl; |
---|
220 | if ( m_crash_on_problem ) exit(1); |
---|
221 | } |
---|
222 | } |
---|
223 | } |
---|
224 | |
---|
225 | // Done. |
---|
226 | return true; |
---|
227 | |
---|
228 | } |
---|
229 | |
---|
230 | //-------------------------------------------------------------------------- |
---|
231 | |
---|
232 | // Conversion from PYTHIA event to HepMC event, with PDF and |
---|
233 | // some other info included. Return T/F = success/failure. |
---|
234 | |
---|
235 | bool I_Pythia8::fill_next_event( Pythia8::Pythia& pythia, GenEvent* evt, |
---|
236 | int ievnum, bool convertGluonTo0 ) { |
---|
237 | |
---|
238 | // If hadronization is switched off then do not warn for free partons. |
---|
239 | bool doHadr = pythia.flag("HadronLevel:all") && |
---|
240 | pythia.flag("HadronLevel:Hadronize"); |
---|
241 | if (!doHadr) m_freepartonwarnings = false; |
---|
242 | |
---|
243 | // Let the method above convert the event record. Check that it worked. |
---|
244 | bool result = fill_next_event( pythia.event, evt, ievnum ); |
---|
245 | if ( result ) { |
---|
246 | |
---|
247 | // Store PDF information. |
---|
248 | put_pdf_info(evt, pythia, convertGluonTo0 ); |
---|
249 | |
---|
250 | // Store process code, scale, alpha_em, alpha_s. |
---|
251 | evt->set_signal_process_id(pythia.info.code()); |
---|
252 | evt->set_event_scale(pythia.info.pTHat()); |
---|
253 | if (evt->alphaQED() <= 0) evt->set_alphaQED( pythia.info.alphaEM() ); |
---|
254 | if (evt->alphaQCD() <= 0) evt->set_alphaQCD( pythia.info.alphaS() ); |
---|
255 | |
---|
256 | // Store event weight, which is in units of fb for Les Houches Event |
---|
257 | // strategies +-4. |
---|
258 | evt->weights().push_back( pythia.info.weight() ); |
---|
259 | |
---|
260 | // HepMC 2.05 supports cross-section information in pb: |
---|
261 | #ifdef HEPMC_HAS_CROSS_SECTION |
---|
262 | HepMC::GenCrossSection xsec; |
---|
263 | xsec.set_cross_section( pythia.info.sigmaGen() * 1e9, |
---|
264 | pythia.info.sigmaErr() * 1e9); |
---|
265 | evt->set_cross_section(xsec); |
---|
266 | #endif |
---|
267 | |
---|
268 | } |
---|
269 | |
---|
270 | // Done. |
---|
271 | return result; |
---|
272 | |
---|
273 | } |
---|
274 | |
---|
275 | //-------------------------------------------------------------------------- |
---|
276 | |
---|
277 | // Conversion of PDF info. |
---|
278 | |
---|
279 | void I_Pythia8::put_pdf_info( GenEvent* evt, Pythia8::Pythia& pythia, |
---|
280 | bool convertGluonTo0 ) { |
---|
281 | |
---|
282 | // Flavours of incoming partons. |
---|
283 | int id1pdf = pythia.info.id1pdf(); |
---|
284 | int id2pdf = pythia.info.id2pdf(); |
---|
285 | if ( convertGluonTo0 ) { |
---|
286 | if ( id1pdf == 21 ) id1pdf = 0; |
---|
287 | if ( id2pdf == 21 ) id2pdf = 0; |
---|
288 | } |
---|
289 | |
---|
290 | // x, Q and x*f(x) for incoming partons. |
---|
291 | double x1pdf = pythia.info.x1pdf(); |
---|
292 | double x2pdf = pythia.info.x2pdf(); |
---|
293 | double Qfac = pythia.info.QFac(); |
---|
294 | double pdf1 = pythia.info.pdf1(); |
---|
295 | double pdf2 = pythia.info.pdf2(); |
---|
296 | |
---|
297 | // Store PDF info and done. |
---|
298 | evt->set_pdf_info( PdfInfo( id1pdf, id2pdf, x1pdf, x2pdf, Qfac, pdf1, pdf2) ); |
---|
299 | |
---|
300 | } |
---|
301 | |
---|
302 | //========================================================================== |
---|
303 | |
---|
304 | } // end namespace HepMC |
---|