hepmc - Blame information for rev 43

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