Line data Source code
1 : //////////////////////////////////////////////////////////////////////////
2 : // Matt.Dobbs@Cern.CH, January 2000
3 : // HEPEVT IO class
4 : //////////////////////////////////////////////////////////////////////////
5 :
6 : #include "HepMC/IO_HEPEVT.h"
7 : #include "HepMC/GenEvent.h"
8 : #include <cstdio> // needed for formatted output using sprintf
9 :
10 : namespace HepMC {
11 :
12 0 : IO_HEPEVT::IO_HEPEVT() : m_trust_mothers_before_daughters(1),
13 0 : m_trust_both_mothers_and_daughters(0),
14 0 : m_print_inconsistency_errors(1),
15 0 : m_trust_beam_particles(true)
16 0 : {}
17 :
18 0 : IO_HEPEVT::~IO_HEPEVT(){}
19 :
20 : void IO_HEPEVT::print( std::ostream& ostr ) const {
21 0 : ostr << "IO_HEPEVT: reads an event from the FORTRAN HEPEVT "
22 0 : << "common block. \n"
23 0 : << " trust_mothers_before_daughters = "
24 0 : << m_trust_mothers_before_daughters
25 0 : << " trust_both_mothers_and_daughters = "
26 0 : << m_trust_both_mothers_and_daughters
27 0 : << ", print_inconsistency_errors = "
28 0 : << m_print_inconsistency_errors << std::endl;
29 0 : }
30 :
31 : bool IO_HEPEVT::fill_next_event( GenEvent* evt ) {
32 : /// read one event from the HEPEVT common block and fill GenEvent
33 : /// return T/F =success/failure
34 : ///
35 : /// For HEPEVT commons built with the luhepc routine of Pythia 5.7
36 : /// the children pointers are not always correct (i.e. there is
37 : /// oftentimes an internal inconsistency between the parents and
38 : /// children pointers). The parent pointers always seem to be correct.
39 : /// Thus the switch trust_mothers_before_daughters=1 is appropriate for
40 : /// pythia. NOTE: you should also set the switch MSTP(128) = 2 in
41 : /// pythia (not the default!), so that pythia doesn't
42 : /// store two copies of resonances in the event record.
43 : /// The situation is opposite for the HEPEVT which comes from Isajet
44 : /// via stdhep, so then use the switch trust_mothers_before_daughters=0
45 : //
46 : // 1. test that evt pointer is not null and set event number
47 0 : if ( !evt ) {
48 : std::cerr
49 0 : << "IO_HEPEVT::fill_next_event error - passed null event."
50 0 : << std::endl;
51 0 : return false;
52 : }
53 0 : evt->set_event_number( HEPEVT_Wrapper::event_number() );
54 : //
55 : // 2. create a particle instance for each HEPEVT entry and fill a map
56 : // create a vector which maps from the HEPEVT particle index to the
57 : // GenParticle address
58 : // (+1 in size accounts for hepevt_particle[0] which is unfilled)
59 0 : std::vector<GenParticle*> hepevt_particle(
60 0 : HEPEVT_Wrapper::number_entries()+1 );
61 0 : hepevt_particle[0] = 0;
62 0 : for ( int i1 = 1; i1 <= HEPEVT_Wrapper::number_entries(); ++i1 ) {
63 0 : hepevt_particle[i1] = build_particle(i1);
64 : }
65 0 : std::set<GenVertex*> new_vertices;
66 : //
67 : // Here we assume that the first two particles in the list
68 : // are the incoming beam particles.
69 0 : if( trust_beam_particles() ) {
70 0 : evt->set_beam_particles( hepevt_particle[1], hepevt_particle[2] );
71 : }
72 : //
73 : // 3.+4. loop over HEPEVT particles AGAIN, this time creating vertices
74 0 : for ( int i = 1; i <= HEPEVT_Wrapper::number_entries(); ++i ) {
75 : // We go through and build EITHER the production or decay
76 : // vertex for each entry in hepevt, depending on the switch
77 : // m_trust_mothers_before_daughters (new 2001-02-28)
78 : // Note: since the HEPEVT pointers are bi-directional, it is
79 : /// sufficient to do one or the other.
80 : //
81 : // 3. Build the production_vertex (if necessary)
82 0 : if ( m_trust_mothers_before_daughters ||
83 0 : m_trust_both_mothers_and_daughters ) {
84 0 : build_production_vertex( i, hepevt_particle, evt );
85 : }
86 : //
87 : // 4. Build the end_vertex (if necessary)
88 : // Identical steps as for production vertex
89 0 : if ( !m_trust_mothers_before_daughters ||
90 0 : m_trust_both_mothers_and_daughters ) {
91 0 : build_end_vertex( i, hepevt_particle, evt );
92 : }
93 : }
94 : // 5. 01.02.2000
95 : // handle the case of particles in HEPEVT which come from nowhere -
96 : // i.e. particles without mothers or daughters.
97 : // These particles need to be attached to a vertex, or else they
98 : // will never become part of the event. check for this situation
99 0 : for ( int i3 = 1; i3 <= HEPEVT_Wrapper::number_entries(); ++i3 ) {
100 0 : if ( !hepevt_particle[i3]->end_vertex() &&
101 0 : !hepevt_particle[i3]->production_vertex() ) {
102 0 : GenVertex* prod_vtx = new GenVertex();
103 0 : prod_vtx->add_particle_out( hepevt_particle[i3] );
104 0 : evt->add_vertex( prod_vtx );
105 0 : }
106 : }
107 : return true;
108 0 : }
109 :
110 : void IO_HEPEVT::write_event( const GenEvent* evt ) {
111 : /// This writes an event out to the HEPEVT common block. The daughters
112 : /// field is NOT filled, because it is possible to contruct graphs
113 : /// for which the mothers and daughters cannot both be make sequential.
114 : /// This is consistent with how pythia fills HEPEVT (daughters are not
115 : /// necessarily filled properly) and how IO_HEPEVT reads HEPEVT.
116 : //
117 0 : if ( !evt ) return;
118 : //
119 : // map all particles onto a unique index
120 0 : std::vector<GenParticle*> index_to_particle(
121 0 : HEPEVT_Wrapper::max_number_entries()+1 );
122 0 : index_to_particle[0]=0;
123 0 : std::map<GenParticle*,int> particle_to_index;
124 : int particle_counter=0;
125 0 : for ( GenEvent::vertex_const_iterator v = evt->vertices_begin();
126 0 : v != evt->vertices_end(); ++v ) {
127 : // all "mothers" or particles_in are kept adjacent in the list
128 : // so that the mother indices in hepevt can be filled properly
129 0 : for ( GenVertex::particles_in_const_iterator p1
130 0 : = (*v)->particles_in_const_begin();
131 0 : p1 != (*v)->particles_in_const_end(); ++p1 ) {
132 0 : ++particle_counter;
133 0 : if ( particle_counter >
134 0 : HEPEVT_Wrapper::max_number_entries() ) break;
135 0 : index_to_particle[particle_counter] = *p1;
136 0 : particle_to_index[*p1] = particle_counter;
137 : }
138 : // daughters are entered only if they aren't a mother of
139 : // another vtx
140 0 : for ( GenVertex::particles_out_const_iterator p2
141 0 : = (*v)->particles_out_const_begin();
142 0 : p2 != (*v)->particles_out_const_end(); ++p2 ) {
143 0 : if ( !(*p2)->end_vertex() ) {
144 0 : ++particle_counter;
145 0 : if ( particle_counter >
146 0 : HEPEVT_Wrapper::max_number_entries() ) {
147 0 : break;
148 : }
149 0 : index_to_particle[particle_counter] = *p2;
150 0 : particle_to_index[*p2] = particle_counter;
151 0 : }
152 : }
153 : }
154 0 : if ( particle_counter > HEPEVT_Wrapper::max_number_entries() ) {
155 0 : particle_counter = HEPEVT_Wrapper::max_number_entries();
156 0 : }
157 : //
158 : // fill the HEPEVT event record
159 0 : HEPEVT_Wrapper::set_event_number( evt->event_number() );
160 0 : HEPEVT_Wrapper::set_number_entries( particle_counter );
161 0 : for ( int i = 1; i <= particle_counter; ++i ) {
162 0 : HEPEVT_Wrapper::set_status( i, index_to_particle[i]->status() );
163 0 : HEPEVT_Wrapper::set_id( i, index_to_particle[i]->pdg_id() );
164 0 : FourVector m = index_to_particle[i]->momentum();
165 0 : HEPEVT_Wrapper::set_momentum( i, m.px(), m.py(), m.pz(), m.e() );
166 0 : HEPEVT_Wrapper::set_mass( i, index_to_particle[i]->generatedMass() );
167 : // there should ALWAYS be particles in any vertex, but some generators
168 : // are making non-kosher HepMC events
169 0 : if ( index_to_particle[i]->production_vertex() &&
170 0 : index_to_particle[i]->production_vertex()->particles_in_size()) {
171 0 : FourVector p = index_to_particle[i]->
172 : production_vertex()->position();
173 0 : HEPEVT_Wrapper::set_position( i, p.x(), p.y(), p.z(), p.t() );
174 0 : int num_mothers = index_to_particle[i]->production_vertex()->
175 : particles_in_size();
176 0 : int first_mother = find_in_map( particle_to_index,
177 0 : *(index_to_particle[i]->
178 : production_vertex()->
179 : particles_in_const_begin()));
180 0 : int last_mother = first_mother + num_mothers - 1;
181 0 : if ( first_mother == 0 ) last_mother = 0;
182 0 : HEPEVT_Wrapper::set_parents( i, first_mother, last_mother );
183 0 : } else {
184 0 : HEPEVT_Wrapper::set_position( i, 0, 0, 0, 0 );
185 0 : HEPEVT_Wrapper::set_parents( i, 0, 0 );
186 : }
187 0 : HEPEVT_Wrapper::set_children( i, 0, 0 );
188 0 : }
189 0 : }
190 :
191 : void IO_HEPEVT::build_production_vertex(int i,
192 : std::vector<HepMC::GenParticle*>&
193 : hepevt_particle,
194 : GenEvent* evt ) {
195 : ///
196 : /// for particle in HEPEVT with index i, build a production vertex
197 : /// if appropriate, and add that vertex to the event
198 0 : GenParticle* p = hepevt_particle[i];
199 : // a. search to see if a production vertex already exists
200 0 : int mother = HEPEVT_Wrapper::first_parent(i);
201 0 : GenVertex* prod_vtx = p->production_vertex();
202 0 : while ( !prod_vtx && mother > 0 ) {
203 0 : prod_vtx = hepevt_particle[mother]->end_vertex();
204 0 : if ( prod_vtx ) prod_vtx->add_particle_out( p );
205 : // increment mother for next iteration
206 0 : if ( ++mother > HEPEVT_Wrapper::last_parent(i) ) mother = 0;
207 : }
208 : // b. if no suitable production vertex exists - and the particle
209 : // has atleast one mother or position information to store -
210 : // make one
211 0 : FourVector prod_pos( HEPEVT_Wrapper::x(i), HEPEVT_Wrapper::y(i),
212 0 : HEPEVT_Wrapper::z(i), HEPEVT_Wrapper::t(i)
213 : );
214 0 : if ( !prod_vtx && (HEPEVT_Wrapper::number_parents(i)>0
215 0 : || prod_pos!=FourVector(0,0,0,0)) )
216 : {
217 0 : prod_vtx = new GenVertex();
218 0 : prod_vtx->add_particle_out( p );
219 0 : evt->add_vertex( prod_vtx );
220 0 : }
221 : // c. if prod_vtx doesn't already have position specified, fill it
222 0 : if ( prod_vtx && prod_vtx->position()==FourVector(0,0,0,0) ) {
223 0 : prod_vtx->set_position( prod_pos );
224 0 : }
225 : // d. loop over mothers to make sure their end_vertices are
226 : // consistent
227 0 : mother = HEPEVT_Wrapper::first_parent(i);
228 0 : while ( prod_vtx && mother > 0 ) {
229 0 : if ( !hepevt_particle[mother]->end_vertex() ) {
230 : // if end vertex of the mother isn't specified, do it now
231 0 : prod_vtx->add_particle_in( hepevt_particle[mother] );
232 0 : } else if (hepevt_particle[mother]->end_vertex() != prod_vtx ) {
233 : // problem scenario --- the mother already has a decay
234 : // vertex which differs from the daughter's produciton
235 : // vertex. This means there is internal
236 : // inconsistency in the HEPEVT event record. Print an
237 : // error
238 : // Note: we could provide a fix by joining the two
239 : // vertices with a dummy particle if the problem
240 : // arrises often with any particular generator.
241 0 : if ( m_print_inconsistency_errors ) std::cerr
242 0 : << "HepMC::IO_HEPEVT: inconsistent mother/daugher "
243 0 : << "information in HEPEVT event "
244 0 : << HEPEVT_Wrapper::event_number()
245 0 : << ". \n I recommend you try "
246 0 : << "inspecting the event first with "
247 0 : << "\n\tHEPEVT_Wrapper::check_hepevt_consistency()"
248 0 : << "\n This warning can be turned off with the "
249 0 : << "IO_HEPEVT::print_inconsistency_errors switch."
250 0 : << std::endl;
251 : }
252 0 : if ( ++mother > HEPEVT_Wrapper::last_parent(i) ) mother = 0;
253 : }
254 0 : }
255 :
256 : void IO_HEPEVT::build_end_vertex
257 : ( int i, std::vector<HepMC::GenParticle*>& hepevt_particle, GenEvent* evt )
258 : {
259 : ///
260 : /// for particle in HEPEVT with index i, build an end vertex
261 : /// if appropriate, and add that vertex to the event
262 : // Identical steps as for build_production_vertex
263 0 : GenParticle* p = hepevt_particle[i];
264 : // a.
265 0 : int daughter = HEPEVT_Wrapper::first_child(i);
266 0 : GenVertex* end_vtx = p->end_vertex();
267 0 : while ( !end_vtx && daughter > 0 ) {
268 0 : end_vtx = hepevt_particle[daughter]->production_vertex();
269 0 : if ( end_vtx ) end_vtx->add_particle_in( p );
270 0 : if ( ++daughter > HEPEVT_Wrapper::last_child(i) ) daughter = 0;
271 : }
272 : // b. (different from 3c. because HEPEVT particle can not know its
273 : // decay position )
274 0 : if ( !end_vtx && HEPEVT_Wrapper::number_children(i)>0 ) {
275 0 : end_vtx = new GenVertex();
276 0 : end_vtx->add_particle_in( p );
277 0 : evt->add_vertex( end_vtx );
278 0 : }
279 : // c+d. loop over daughters to make sure their production vertices
280 : // point back to the current vertex.
281 : // We get the vertex position from the daughter as well.
282 0 : daughter = HEPEVT_Wrapper::first_child(i);
283 0 : while ( end_vtx && daughter > 0 ) {
284 0 : if ( !hepevt_particle[daughter]->production_vertex() ) {
285 : // if end vertex of the mother isn't specified, do it now
286 0 : end_vtx->add_particle_out( hepevt_particle[daughter] );
287 : //
288 : // 2001-03-29 M.Dobbs, fill vertex the position.
289 0 : if ( end_vtx->position()==FourVector(0,0,0,0) ) {
290 0 : FourVector prod_pos( HEPEVT_Wrapper::x(daughter),
291 0 : HEPEVT_Wrapper::y(daughter),
292 0 : HEPEVT_Wrapper::z(daughter),
293 0 : HEPEVT_Wrapper::t(daughter)
294 : );
295 0 : if ( prod_pos != FourVector(0,0,0,0) ) {
296 0 : end_vtx->set_position( prod_pos );
297 0 : }
298 0 : }
299 0 : } else if (hepevt_particle[daughter]->production_vertex()
300 0 : != end_vtx){
301 : // problem scenario --- the daughter already has a prod
302 : // vertex which differs from the mother's end
303 : // vertex. This means there is internal
304 : // inconsistency in the HEPEVT event record. Print an
305 : // error
306 0 : if ( m_print_inconsistency_errors ) std::cerr
307 0 : << "HepMC::IO_HEPEVT: inconsistent mother/daugher "
308 0 : << "information in HEPEVT event "
309 0 : << HEPEVT_Wrapper::event_number()
310 0 : << ". \n I recommend you try "
311 0 : << "inspecting the event first with "
312 0 : << "\n\tHEPEVT_Wrapper::check_hepevt_consistency()"
313 0 : << "\n This warning can be turned off with the "
314 0 : << "IO_HEPEVT::print_inconsistency_errors switch."
315 0 : << std::endl;
316 : }
317 0 : if ( ++daughter > HEPEVT_Wrapper::last_child(i) ) daughter = 0;
318 : }
319 0 : if ( !p->end_vertex() && !p->production_vertex() ) {
320 : // Added 2001-11-04, to try and handle Isajet problems.
321 0 : build_production_vertex( i, hepevt_particle, evt );
322 0 : }
323 0 : }
324 :
325 : GenParticle* IO_HEPEVT::build_particle( int index ) {
326 : /// Builds a particle object corresponding to index in HEPEVT
327 : //
328 : GenParticle* p
329 0 : = new GenParticle( FourVector( HEPEVT_Wrapper::px(index),
330 0 : HEPEVT_Wrapper::py(index),
331 0 : HEPEVT_Wrapper::pz(index),
332 0 : HEPEVT_Wrapper::e(index) ),
333 0 : HEPEVT_Wrapper::id(index),
334 0 : HEPEVT_Wrapper::status(index) );
335 0 : p->setGeneratedMass( HEPEVT_Wrapper::m(index) );
336 0 : p->suggest_barcode( index );
337 0 : return p;
338 0 : }
339 :
340 : int IO_HEPEVT::find_in_map( const std::map<HepMC::GenParticle*,int>& m,
341 : GenParticle* p) const {
342 0 : std::map<GenParticle*,int>::const_iterator iter = m.find(p);
343 0 : if ( iter == m.end() ) return 0;
344 0 : return iter->second;
345 0 : }
346 :
347 : } // HepMC
348 :
349 :
350 :
|
