TEvtGen/HepMC/IO_HEPEVT.cc

   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                              m_trust_beam_particles(true)
  16     {}
  17
  18     IO_HEPEVT::~IO_HEPEVT(){}
  19
  20     void IO_HEPEVT::print( std::ostream& ostr ) const {
  21         ostr << "IO_HEPEVT: reads an event from the FORTRAN HEPEVT "
  22              << "common block. \n"
  23              << " trust_mothers_before_daughters = "
  24              << m_trust_mothers_before_daughters
  25              << " trust_both_mothers_and_daughters = "
  26              << m_trust_both_mothers_and_daughters
  27              << ", print_inconsistency_errors = "
  28              << m_print_inconsistency_errors << std::endl;
  29     }
  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         if ( !evt ) {
  48             std::cerr
  49                 << "IO_HEPEVT::fill_next_event error - passed null event."
  50                 << std::endl;
  51             return false;
  52         }
  53         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         std::vector<GenParticle*> hepevt_particle(
  60                                         HEPEVT_Wrapper::number_entries()+1 );
  61         hepevt_particle[0] = 0;
  62         for ( int i1 = 1; i1 <= HEPEVT_Wrapper::number_entries(); ++i1 ) {
  63             hepevt_particle[i1] = build_particle(i1);
  64         }
  65         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         if( trust_beam_particles() ) {
  70         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         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             if ( m_trust_mothers_before_daughters ||
  83                  m_trust_both_mothers_and_daughters ) {
  84                 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             if ( !m_trust_mothers_before_daughters ||
  90                  m_trust_both_mothers_and_daughters ) {
  91                 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         for ( int i3 = 1; i3 <= HEPEVT_Wrapper::number_entries(); ++i3 ) {
 100             if ( !hepevt_particle[i3]->end_vertex() &&
 101                         !hepevt_particle[i3]->production_vertex() ) {
 102                 GenVertex* prod_vtx = new GenVertex();
 103                 prod_vtx->add_particle_out( hepevt_particle[i3] );
 104                 evt->add_vertex( prod_vtx );
 105             }
 106         }
 107         return true;
 108     }
 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         if ( !evt ) return;
 118         //
 119         // map all particles onto a unique index
 120         std::vector<GenParticle*> index_to_particle(
 121             HEPEVT_Wrapper::max_number_entries()+1 );
 122         index_to_particle[0]=0;
 123         std::map<GenParticle*,int> particle_to_index;
 124         int particle_counter=0;
 125         for ( GenEvent::vertex_const_iterator v = evt->vertices_begin();
 126               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             for ( GenVertex::particles_in_const_iterator p1
 130                       = (*v)->particles_in_const_begin();
 131                   p1 != (*v)->particles_in_const_end(); ++p1 ) {
 132                 ++particle_counter;
 133                 if ( particle_counter >
 134                      HEPEVT_Wrapper::max_number_entries() ) break;
 135                 index_to_particle[particle_counter] = *p1;
 136                 particle_to_index[*p1] = particle_counter;
 137             }
 138             // daughters are entered only if they aren't a mother of
 139             // another vtx
 140             for ( GenVertex::particles_out_const_iterator p2
 141                       = (*v)->particles_out_const_begin();
 142                   p2 != (*v)->particles_out_const_end(); ++p2 ) {
 143                 if ( !(*p2)->end_vertex() ) {
 144                     ++particle_counter;
 145                     if ( particle_counter >
 146                          HEPEVT_Wrapper::max_number_entries() ) {
 147                         break;
 148                     }
 149                     index_to_particle[particle_counter] = *p2;
 150                     particle_to_index[*p2] = particle_counter;
 151                 }
 152             }
 153         }
 154         if ( particle_counter > HEPEVT_Wrapper::max_number_entries() ) {
 155             particle_counter = HEPEVT_Wrapper::max_number_entries();
 156         }
 157         //
 158         // fill the HEPEVT event record
 159         HEPEVT_Wrapper::set_event_number( evt->event_number() );
 160         HEPEVT_Wrapper::set_number_entries( particle_counter );
 161         for ( int i = 1; i <= particle_counter; ++i ) {
 162             HEPEVT_Wrapper::set_status( i, index_to_particle[i]->status() );
 163             HEPEVT_Wrapper::set_id( i, index_to_particle[i]->pdg_id() );
 164             FourVector m = index_to_particle[i]->momentum();
 165             HEPEVT_Wrapper::set_momentum( i, m.px(), m.py(), m.pz(), m.e() );
 166             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             if ( index_to_particle[i]->production_vertex() &&
 170                  index_to_particle[i]->production_vertex()->particles_in_size()) {
 171                 FourVector p = index_to_particle[i]->
 172                                      production_vertex()->position();
 173                 HEPEVT_Wrapper::set_position( i, p.x(), p.y(), p.z(), p.t() );
 174                 int num_mothers = index_to_particle[i]->production_vertex()->
 175                                   particles_in_size();
 176                 int first_mother = find_in_map( particle_to_index,
 177                                                 *(index_to_particle[i]->
 178                                                   production_vertex()->
 179                                                   particles_in_const_begin()));
 180                 int last_mother = first_mother + num_mothers - 1;
 181                 if ( first_mother == 0 ) last_mother = 0;
 182                 HEPEVT_Wrapper::set_parents( i, first_mother, last_mother );
 183             } else {
 184                 HEPEVT_Wrapper::set_position( i, 0, 0, 0, 0 );
 185                 HEPEVT_Wrapper::set_parents( i, 0, 0 );
 186             }
 187             HEPEVT_Wrapper::set_children( i, 0, 0 );
 188         }
 189     }
 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         GenParticle* p = hepevt_particle[i];
 199         // a. search to see if a production vertex already exists
 200         int mother = HEPEVT_Wrapper::first_parent(i);
 201         GenVertex* prod_vtx = p->production_vertex();
 202         while ( !prod_vtx && mother > 0 ) {
 203             prod_vtx = hepevt_particle[mother]->end_vertex();
 204             if ( prod_vtx ) prod_vtx->add_particle_out( p );
 205             // increment mother for next iteration
 206             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         FourVector prod_pos( HEPEVT_Wrapper::x(i), HEPEVT_Wrapper::y(i),
 212                                    HEPEVT_Wrapper::z(i), HEPEVT_Wrapper::t(i)
 213                                  );
 214         if ( !prod_vtx && (HEPEVT_Wrapper::number_parents(i)>0
 215                            || prod_pos!=FourVector(0,0,0,0)) )
 216         {
 217             prod_vtx = new GenVertex();
 218             prod_vtx->add_particle_out( p );
 219             evt->add_vertex( prod_vtx );
 220         }
 221         // c. if prod_vtx doesn't already have position specified, fill it
 222         if ( prod_vtx && prod_vtx->position()==FourVector(0,0,0,0) ) {
 223             prod_vtx->set_position( prod_pos );
 224         }
 225         // d. loop over mothers to make sure their end_vertices are
 226         //     consistent
 227         mother = HEPEVT_Wrapper::first_parent(i);
 228         while ( prod_vtx && mother > 0 ) {
 229             if ( !hepevt_particle[mother]->end_vertex() ) {
 230                 // if end vertex of the mother isn't specified, do it now
 231                 prod_vtx->add_particle_in( hepevt_particle[mother] );
 232             } 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                 if ( m_print_inconsistency_errors ) std::cerr
 242                     << "HepMC::IO_HEPEVT: inconsistent mother/daugher "
 243                     << "information in HEPEVT event "
 244                     << HEPEVT_Wrapper::event_number()
 245                     << ". \n I recommend you try "
 246                     << "inspecting the event first with "
 247                     << "\n\tHEPEVT_Wrapper::check_hepevt_consistency()"
 248                     << "\n This warning can be turned off with the "
 249                     << "IO_HEPEVT::print_inconsistency_errors switch."
 250                     << std::endl;
 251             }
 252             if ( ++mother > HEPEVT_Wrapper::last_parent(i) ) mother = 0;
 253         }
 254     }
 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         GenParticle* p = hepevt_particle[i];
 264         // a.
 265         int daughter = HEPEVT_Wrapper::first_child(i);
 266         GenVertex* end_vtx = p->end_vertex();
 267         while ( !end_vtx && daughter > 0 ) {
 268             end_vtx = hepevt_particle[daughter]->production_vertex();
 269             if ( end_vtx ) end_vtx->add_particle_in( p );
 270             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         if ( !end_vtx && HEPEVT_Wrapper::number_children(i)>0 ) {
 275             end_vtx = new GenVertex();
 276             end_vtx->add_particle_in( p );
 277             evt->add_vertex( end_vtx );
 278         }
 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         daughter = HEPEVT_Wrapper::first_child(i);
 283         while ( end_vtx && daughter > 0 ) {
 284             if ( !hepevt_particle[daughter]->production_vertex() ) {
 285                 // if end vertex of the mother isn't specified, do it now
 286                 end_vtx->add_particle_out( hepevt_particle[daughter] );
 287                 //
 288                 // 2001-03-29 M.Dobbs, fill vertex the position.
 289                 if ( end_vtx->position()==FourVector(0,0,0,0) ) {
 290                     FourVector prod_pos( HEPEVT_Wrapper::x(daughter),
 291                                                HEPEVT_Wrapper::y(daughter),
 292                                                HEPEVT_Wrapper::z(daughter),
 293                                                HEPEVT_Wrapper::t(daughter)
 294                         );
 295                     if ( prod_pos != FourVector(0,0,0,0) ) {
 296                         end_vtx->set_position( prod_pos );
 297                     }
 298                 }
 299             } else if (hepevt_particle[daughter]->production_vertex()
 300                        != 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                 if ( m_print_inconsistency_errors ) std::cerr
 307                     << "HepMC::IO_HEPEVT: inconsistent mother/daugher "
 308                     << "information in HEPEVT event "
 309                     << HEPEVT_Wrapper::event_number()
 310                     << ". \n I recommend you try "
 311                     << "inspecting the event first with "
 312                     << "\n\tHEPEVT_Wrapper::check_hepevt_consistency()"
 313                     << "\n This warning can be turned off with the "
 314                     << "IO_HEPEVT::print_inconsistency_errors switch."
 315                     << std::endl;
 316             }
 317             if ( ++daughter > HEPEVT_Wrapper::last_child(i) ) daughter = 0;
 318         }
 319         if ( !p->end_vertex() && !p->production_vertex() ) {
 320             // Added 2001-11-04, to try and handle Isajet problems.
 321             build_production_vertex( i, hepevt_particle, evt );
 322         }
 323     }
 324
 325     GenParticle* IO_HEPEVT::build_particle( int index ) {
 326         /// Builds a particle object corresponding to index in HEPEVT
 327         //
 328         GenParticle* p
 329             = new GenParticle( FourVector( HEPEVT_Wrapper::px(index),
 330                                                  HEPEVT_Wrapper::py(index),
 331                                                  HEPEVT_Wrapper::pz(index),
 332                                                  HEPEVT_Wrapper::e(index) ),
 333                                HEPEVT_Wrapper::id(index),
 334                                HEPEVT_Wrapper::status(index) );
 335         p->setGeneratedMass( HEPEVT_Wrapper::m(index) );
 336         p->suggest_barcode( index );
 337         return p;
 338     }
 339
 340     int IO_HEPEVT::find_in_map( const std::map<HepMC::GenParticle*,int>& m,
 341                                 GenParticle* p) const {
 342         std::map<GenParticle*,int>::const_iterator iter = m.find(p);
 343         if ( iter == m.end() ) return 0;
 344         return iter->second;
 345     }
 346
 347 } // HepMC
 348
 349
 350