TEvtGen/HepMC/IO_HERWIG.cc

   1 //////////////////////////////////////////////////////////////////////////
   2 // Matt.Dobbs@Cern.CH, October 2002
   3 // Herwig 6.400 IO class
   4 //////////////////////////////////////////////////////////////////////////
   5
   6 #include "HepMC/IO_HERWIG.h"
   7 #include "HepMC/GenEvent.h"
   8 #include <cstdio>       // needed for formatted output using sprintf
   9
  10 namespace HepMC {
  11
  12     IO_HERWIG::IO_HERWIG() : m_trust_mothers_before_daughters(false),
  13                              m_trust_both_mothers_and_daughters(true),
  14                              m_print_inconsistency_errors(true),
  15                              m_no_gaps_in_barcodes(true),
  16                              m_herwig_to_pdg_id(100,0)
  17     {
  18         // These arrays are copied from Lynn Garren's stdhep 5.01-5.06.
  19         //   see http://cepa.fnal.gov/psm/stdhep/
  20         // Translation from HERWIG particle ID's to PDG particle ID's.
  21         m_herwig_to_pdg_id[1] =1;
  22         m_herwig_to_pdg_id[2] =2;
  23         m_herwig_to_pdg_id[3] =3;
  24         m_herwig_to_pdg_id[4] =4;
  25         m_herwig_to_pdg_id[5] =5;
  26         m_herwig_to_pdg_id[6] =6;
  27         m_herwig_to_pdg_id[7] =7;
  28         m_herwig_to_pdg_id[8] =8;
  29
  30         m_herwig_to_pdg_id[11] =11;
  31         m_herwig_to_pdg_id[12] =12;
  32         m_herwig_to_pdg_id[13] =13;
  33         m_herwig_to_pdg_id[14] =14;
  34         m_herwig_to_pdg_id[15] =15;
  35         m_herwig_to_pdg_id[16] =16;
  36
  37         m_herwig_to_pdg_id[21] =21;
  38         m_herwig_to_pdg_id[22] =22;
  39         m_herwig_to_pdg_id[23] =23;
  40         m_herwig_to_pdg_id[24] =24;
  41         m_herwig_to_pdg_id[25] =25;
  42         m_herwig_to_pdg_id[26] =51; // <-- H_L0 (redundant with h0(25))
  43
  44         m_herwig_to_pdg_id[32] =32;
  45         m_herwig_to_pdg_id[35] =35;
  46         m_herwig_to_pdg_id[36] =36;
  47         m_herwig_to_pdg_id[37] =37;
  48         m_herwig_to_pdg_id[39] =39;
  49
  50         m_herwig_to_pdg_id[40] =40; //Charybdis Black Hole
  51
  52         m_herwig_to_pdg_id[81] =81;
  53         m_herwig_to_pdg_id[82] =82;
  54         m_herwig_to_pdg_id[83] =83;
  55         m_herwig_to_pdg_id[84] =84;
  56         m_herwig_to_pdg_id[85] =85;
  57         m_herwig_to_pdg_id[86] =86;
  58         m_herwig_to_pdg_id[87] =87;
  59         m_herwig_to_pdg_id[88] =88;
  60         m_herwig_to_pdg_id[89] =89;
  61         m_herwig_to_pdg_id[90] =90;
  62
  63         m_herwig_to_pdg_id[91] =91;
  64         m_herwig_to_pdg_id[92] =92;
  65         m_herwig_to_pdg_id[93] =93;
  66         m_herwig_to_pdg_id[94] =94;
  67         m_herwig_to_pdg_id[95] =95;
  68         m_herwig_to_pdg_id[96] =96;
  69         m_herwig_to_pdg_id[97] =97;
  70         m_herwig_to_pdg_id[98] =9920022; // <-- remnant photon
  71         m_herwig_to_pdg_id[99] =9922212; // <-- remnant nucleon
  72
  73         // These particle ID's have no antiparticle, so aren't allowed.
  74         m_no_antiparticles.insert(-21);
  75         m_no_antiparticles.insert(-22);
  76         m_no_antiparticles.insert(-23);
  77         m_no_antiparticles.insert(-25);
  78         m_no_antiparticles.insert(-51);
  79         m_no_antiparticles.insert(-35);
  80         m_no_antiparticles.insert(-36);
  81     }
  82
  83     IO_HERWIG::~IO_HERWIG(){}
  84
  85     void IO_HERWIG::print( std::ostream& ostr ) const {
  86         ostr << "IO_HERWIG: reads an event from the FORTRAN Herwig HEPEVT "
  87              << "common block. \n"
  88              << " trust_mothers_before_daughters = "
  89              << m_trust_mothers_before_daughters
  90              << " trust_both_mothers_and_daughters = "
  91              << m_trust_both_mothers_and_daughters
  92              << " print_inconsistency_errors = "
  93              << m_print_inconsistency_errors << std::endl;
  94     }
  95
  96     bool IO_HERWIG::fill_next_event( GenEvent* evt ) {
  97         /// read one event from the Herwig HEPEVT common block and fill GenEvent
  98         /// return T/F =success/failure
  99         //
 100         // 0. Test that evt pointer is not null and set event number
 101         if ( !evt ) {
 102             std::cerr
 103                 << "IO_HERWIG::fill_next_event error - passed null event."
 104                 << std::endl;
 105             return false;
 106         }
 107
 108         // 1. First we have to fix the HEPEVT input, which is all mucked up for
 109         //    herwig.
 110         repair_hepevt();
 111
 112         evt->set_event_number( HEPEVT_Wrapper::event_number() );
 113         // Herwig units are GeV and mm
 114         // It would be nice to set the units right here,
 115         // but this could cause problems with existing code that
 116         // might convert GeV to MeV without calling the appropriate HepMC method
 117
 118         //
 119         // 2. create a particle instance for each HEPEVT entry and fill a map
 120         //    create a vector which maps from the HEPEVT particle index to the
 121         //    GenParticle address
 122         //    (+1 in size accounts for hepevt_particle[0] which is unfilled)
 123         std::vector<GenParticle*> hepevt_particle(
 124                                         HEPEVT_Wrapper::number_entries()+1 );
 125         hepevt_particle[0] = 0;
 126         for ( int i1 = 1; i1 <= HEPEVT_Wrapper::number_entries(); ++i1 ) {
 127             hepevt_particle[i1] = build_particle(i1);
 128         }
 129         std::set<GenVertex*> new_vertices;
 130         //
 131         // Here we assume that the first two particles in the list
 132         // are the incoming beam particles.
 133         // Best make sure this is done before any rearranging...
 134         evt->set_beam_particles( hepevt_particle[1], hepevt_particle[2] );
 135         //
 136         // 3. We need to take special care with the hard process
 137         // vertex.  The problem we are trying to avoid is when the
 138         // partons entering the hard process also have daughters from
 139         // the parton shower. When this happens, each one can get its
 140         // own decay vertex, making it difficult to join them
 141         // later. We handle it by joining them together first, then
 142         // the other daughters get added on later.
 143         // Find the partons entering the hard vertex (status codes 121, 122).
 144         int index_121 = 0;
 145         int index_122 = 0;
 146         for ( int i = 1; i <=HEPEVT_Wrapper::number_entries(); i++ ) {
 147             if ( HEPEVT_Wrapper::status(i)==121 ) index_121=i;
 148             if ( HEPEVT_Wrapper::status(i)==122 ) index_122=i;
 149             if ( index_121!=0 && index_122!=0 ) break;
 150         }
 151         if ( index_121 && index_122 ) {
 152             GenVertex* hard_vtx = new GenVertex();
 153             hard_vtx->add_particle_in( hepevt_particle[index_121] );
 154             hard_vtx->add_particle_in( hepevt_particle[index_122] );
 155             // evt->add_vertex( hard_vtx ); // not necessary, its done in
 156                                             // set_signal_process_vertex
 157             //BPK - Atlas -> index_hard retained if it is a boson
 158             int index_hard = 0;
 159             for ( int i = 1; i <=HEPEVT_Wrapper::number_entries(); i++ ) {
 160               if ( HEPEVT_Wrapper::status(i)==120 ) index_hard=i;
 161               if ( index_hard!=0 ) break;
 162             }
 163
 164             if ( index_hard!=0) {
 165               hard_vtx->add_particle_out( hepevt_particle[index_hard] );
 166               GenVertex* hard_vtx2 = new GenVertex();
 167               hard_vtx2->add_particle_in( hepevt_particle[index_hard] );
 168                   for ( int i = 1; i <= HEPEVT_Wrapper::number_entries(); ++i ) {
 169                 if (  HEPEVT_Wrapper::first_parent(i)==index_hard ) {
 170                   hard_vtx2->add_particle_out( hepevt_particle[i] );
 171                 }
 172               }
 173               evt->set_signal_process_vertex( hard_vtx );
 174               evt->set_signal_process_vertex( hard_vtx2 );
 175             }
 176             else {
 177               evt->set_signal_process_vertex( hard_vtx );
 178             }
 179             //BPK - Atlas -<
 180         }
 181         //
 182         // 4. loop over HEPEVT particles AGAIN, this time creating vertices
 183         for ( int i = 1; i <= HEPEVT_Wrapper::number_entries(); ++i ) {
 184             // We go through and build EITHER the production or decay
 185             // vertex for each entry in hepevt, depending on the switch
 186             // m_trust_mothers_before_daughters (new 2001-02-28)
 187             // Note: since the HEPEVT pointers are bi-directional, it is
 188             ///      sufficient to do one or the other.
 189             //
 190             // 3. Build the production_vertex (if necessary)
 191             if ( m_trust_mothers_before_daughters ||
 192                  m_trust_both_mothers_and_daughters ) {
 193                 build_production_vertex( i, hepevt_particle, evt );
 194             }
 195             //
 196             // 4. Build the end_vertex (if necessary)
 197             //    Identical steps as for production vertex
 198             if ( !m_trust_mothers_before_daughters ||
 199                  m_trust_both_mothers_and_daughters ) {
 200                 build_end_vertex( i, hepevt_particle, evt );
 201             }
 202         }
 203         // 5.             01.02.2000
 204         // handle the case of particles in HEPEVT which come from nowhere -
 205         //  i.e. particles without mothers or daughters.
 206         //  These particles need to be attached to a vertex, or else they
 207         //  will never become part of the event. check for this situation.
 208         for ( int i3 = 1; i3 <= HEPEVT_Wrapper::number_entries(); ++i3 ) {
 209             // Herwig also has some non-physical entries in HEPEVT
 210             // like CMS, HARD, and CONE. These are flagged by
 211             // repair_hepevt by making their status and id zero. We
 212             // delete those particles here.
 213             if ( hepevt_particle[i3] && !hepevt_particle[i3]->parent_event()
 214                  && !hepevt_particle[i3]->pdg_id()
 215                  && !hepevt_particle[i3]->status() ) {
 216                 //std::cout << "IO_HERWIG::fill_next_event is deleting null "
 217                 //        << "particle" << std::endl;
 218                 //hepevt_particle[i3]->print();
 219                 delete hepevt_particle[i3];
 220             } else if ( hepevt_particle[i3] &&
 221                         !hepevt_particle[i3]->end_vertex() &&
 222                         !hepevt_particle[i3]->production_vertex() ) {
 223                 GenVertex* prod_vtx = new GenVertex();
 224                 prod_vtx->add_particle_out( hepevt_particle[i3] );
 225                 evt->add_vertex( prod_vtx );
 226             }
 227         }
 228         return true;
 229     }
 230
 231     void IO_HERWIG::build_production_vertex(int i,
 232                                             std::vector<GenParticle*>&
 233                                             hepevt_particle,
 234                                             GenEvent* evt ) {
 235         ///
 236         /// for particle in HEPEVT with index i, build a production vertex
 237         /// if appropriate, and add that vertex to the event
 238         GenParticle* p = hepevt_particle[i];
 239         // a. search to see if a production vertex already exists
 240         int mother = HEPEVT_Wrapper::first_parent(i);
 241         GenVertex* prod_vtx = p->production_vertex();
 242         while ( !prod_vtx && mother > 0 ) {
 243             prod_vtx = hepevt_particle[mother]->end_vertex();
 244             if ( prod_vtx ) prod_vtx->add_particle_out( p );
 245             // increment mother for next iteration
 246             if ( ++mother > HEPEVT_Wrapper::last_parent(i) ) mother = 0;
 247         }
 248         // b. if no suitable production vertex exists - and the particle
 249         // has atleast one mother or position information to store -
 250         // make one
 251         FourVector prod_pos( HEPEVT_Wrapper::x(i), HEPEVT_Wrapper::y(i),
 252                                    HEPEVT_Wrapper::z(i), HEPEVT_Wrapper::t(i)
 253                                  );
 254         if ( !prod_vtx && (HEPEVT_Wrapper::number_parents(i)>0
 255                            || prod_pos!=FourVector(0,0,0,0)) )
 256         {
 257             prod_vtx = new GenVertex();
 258             prod_vtx->add_particle_out( p );
 259             evt->add_vertex( prod_vtx );
 260         }
 261         // c. if prod_vtx doesn't already have position specified, fill it
 262         if ( prod_vtx && prod_vtx->position()==FourVector(0,0,0,0) ) {
 263             prod_vtx->set_position( prod_pos );
 264         }
 265         // d. loop over mothers to make sure their end_vertices are
 266         //     consistent
 267         mother = HEPEVT_Wrapper::first_parent(i);
 268         while ( prod_vtx && mother > 0 ) {
 269             if ( !hepevt_particle[mother]->end_vertex() ) {
 270                 // if end vertex of the mother isn't specified, do it now
 271                 prod_vtx->add_particle_in( hepevt_particle[mother] );
 272             } else if (hepevt_particle[mother]->end_vertex() != prod_vtx ) {
 273                 // problem scenario --- the mother already has a decay
 274                 // vertex which differs from the daughter's produciton
 275                 // vertex. This means there is internal
 276                 // inconsistency in the HEPEVT event record. Print an
 277                 // error
 278                 // Note: we could provide a fix by joining the two
 279                 //       vertices with a dummy particle if the problem
 280                 //       arrises often with any particular generator.
 281                 if ( m_print_inconsistency_errors ) {
 282                   std::cerr
 283                     << "HepMC::IO_HERWIG: inconsistent mother/daugher "
 284                     << "information in HEPEVT event "
 285                     << HEPEVT_Wrapper::event_number()
 286                     << ". \n I recommend you try "
 287                     << "inspecting the event first with "
 288                     << "\n\tHEPEVT_Wrapper::check_hepevt_consistency()"
 289                     << "\n This warning can be turned off with the "
 290                     << "IO_HERWIG::print_inconsistency_errors switch."
 291                     << std::endl;
 292                   hepevt_particle[mother]->print(std::cerr);
 293                   std::cerr
 294                     << "problem vertices are: (prod_vtx, mother)" << std::endl;
 295                   if ( prod_vtx ) prod_vtx->print(std::cerr);
 296                   hepevt_particle[mother]->end_vertex()->print(std::cerr);
 297                 }
 298             }
 299             if ( ++mother > HEPEVT_Wrapper::last_parent(i) ) mother = 0;
 300         }
 301     }
 302
 303     void IO_HERWIG::build_end_vertex
 304     ( int i, std::vector<GenParticle*>& hepevt_particle, GenEvent* evt )
 305     {
 306         ///
 307         /// for particle in HEPEVT with index i, build an end vertex
 308         /// if appropriate, and add that vertex to the event
 309         //    Identical steps as for build_production_vertex
 310         GenParticle* p = hepevt_particle[i];
 311         // a.
 312         int daughter = HEPEVT_Wrapper::first_child(i);
 313         GenVertex* end_vtx = p->end_vertex();
 314         while ( !end_vtx && daughter > 0 ) {
 315             end_vtx = hepevt_particle[daughter]->production_vertex();
 316             if ( end_vtx ) end_vtx->add_particle_in( p );
 317             if ( ++daughter > HEPEVT_Wrapper::last_child(i) ) daughter = 0;
 318         }
 319         // b. (different from 3c. because HEPEVT particle can not know its
 320         //        decay position )
 321         if ( !end_vtx && HEPEVT_Wrapper::number_children(i)>0 ) {
 322             end_vtx = new GenVertex();
 323             end_vtx->add_particle_in( p );
 324             evt->add_vertex( end_vtx );
 325         }
 326         // c+d. loop over daughters to make sure their production vertices
 327         //    point back to the current vertex.
 328         //    We get the vertex position from the daughter as well.
 329         daughter = HEPEVT_Wrapper::first_child(i);
 330         while ( end_vtx && daughter > 0 ) {
 331             if ( !hepevt_particle[daughter]->production_vertex() ) {
 332                 // if end vertex of the mother isn't specified, do it now
 333                 end_vtx->add_particle_out( hepevt_particle[daughter] );
 334                 //
 335                 // 2001-03-29 M.Dobbs, fill vertex the position.
 336                 if ( end_vtx->position()==FourVector(0,0,0,0) ) {
 337                     FourVector prod_pos( HEPEVT_Wrapper::x(daughter),
 338                                                HEPEVT_Wrapper::y(daughter),
 339                                                HEPEVT_Wrapper::z(daughter),
 340                                                HEPEVT_Wrapper::t(daughter)
 341                         );
 342                     if ( prod_pos != FourVector(0,0,0,0) ) {
 343                         end_vtx->set_position( prod_pos );
 344                     }
 345                 }
 346             } else if (hepevt_particle[daughter]->production_vertex()
 347                        != end_vtx){
 348                 // problem scenario --- the daughter already has a prod
 349                 // vertex which differs from the mother's end
 350                 // vertex. This means there is internal
 351                 // inconsistency in the HEPEVT event record. Print an
 352                 // error
 353                 if ( m_print_inconsistency_errors ) std::cerr
 354                     << "HepMC::IO_HERWIG: inconsistent mother/daugher "
 355                     << "information in HEPEVT event "
 356                     << HEPEVT_Wrapper::event_number()
 357                     << ". \n I recommend you try "
 358                     << "inspecting the event first with "
 359                     << "\n\tHEPEVT_Wrapper::check_hepevt_consistency()"
 360                     << "\n This warning can be turned off with the "
 361                     << "IO_HERWIG::print_inconsistency_errors switch."
 362                     << std::endl;
 363             }
 364             if ( ++daughter > HEPEVT_Wrapper::last_child(i) ) daughter = 0;
 365         }
 366         if ( !p->end_vertex() && !p->production_vertex() ) {
 367             // Added 2001-11-04, to try and handle Isajet problems.
 368             build_production_vertex( i, hepevt_particle, evt );
 369         }
 370     }
 371
 372     GenParticle* IO_HERWIG::build_particle( int index ) {
 373         /// Builds a particle object corresponding to index in HEPEVT
 374         //
 375         GenParticle* p
 376             = new GenParticle( FourVector( HEPEVT_Wrapper::px(index),
 377                                                  HEPEVT_Wrapper::py(index),
 378                                                  HEPEVT_Wrapper::pz(index),
 379                                                  HEPEVT_Wrapper::e(index) ),
 380                                HEPEVT_Wrapper::id(index),
 381                                HEPEVT_Wrapper::status(index) );
 382         p->setGeneratedMass( HEPEVT_Wrapper::m(index) );
 383         p->suggest_barcode( index );
 384         return p;
 385     }
 386
 387     int IO_HERWIG::find_in_map( const std::map<GenParticle*,int>& m,
 388                                 GenParticle* p) const {
 389         std::map<GenParticle*,int>::const_iterator iter = m.find(p);
 390         if ( iter == m.end() ) return 0;
 391         return iter->second;
 392     }
 393
 394     void IO_HERWIG::repair_hepevt() const {
 395         ///  This routine takes the HEPEVT common block as used in HERWIG,
 396         ///  and converts it into the HEPEVT common block in the standard format
 397         ///
 398         ///  This means it:
 399         ///    - removes the color structure, which herwig overloads
 400         ///      into the mother/daughter fields
 401         ///    - zeros extra entries for hard subprocess, etc.
 402         ///
 403         ///
 404         /// Special HERWIG status codes
 405         ///   101,102   colliding beam particles
 406         ///   103       beam-beam collision CMS vector
 407         ///   120       hard subprocess CMS vector
 408         ///   121,122   hard subprocess colliding partons
 409         ///   123-129   hard subprocess outgoing particles
 410         ///   141-149   (ID=94) mirror image of hard subrpocess particles
 411         ///   100       (ID=0 cone)
 412         ///
 413         /// Special HERWIG particle id's
 414         ///   91 clusters
 415         ///   94 jets
 416         ///   0  others with no pdg code
 417
 418         // Make sure hepvt isn't empty.
 419         if ( HEPEVT_Wrapper::number_entries() <= 0 ) return;
 420
 421         // Find the index of the beam-beam collision and of the hard subprocess
 422         // Later we will assume that
 423         //              101 ---> 121 \.
 424         //                             X  Hard subprocess
 425         //              102 ---> 122 /
 426         //
 427         int index_collision = 0;
 428         int index_hard = 0;
 429         int index_101 = 0;
 430         int index_102 = 0;
 431         int index_121 = 0;
 432         int index_122 = 0;
 433
 434         for ( int i = 1; i <=HEPEVT_Wrapper::number_entries(); i++ ) {
 435             if ( HEPEVT_Wrapper::status(i)==101 ) index_101=i;
 436             if ( HEPEVT_Wrapper::status(i)==102 ) index_102=i;
 437             if ( HEPEVT_Wrapper::status(i)==103 ) index_collision=i;
 438             if ( HEPEVT_Wrapper::status(i)==120 ) index_hard=i;
 439             if ( HEPEVT_Wrapper::status(i)==121 ) index_121=i;
 440             if ( HEPEVT_Wrapper::status(i)==122 ) index_122=i;
 441             if ( index_collision!=0 && index_hard!=0 && index_101!=0 &&
 442                  index_102!=0 && index_121!=0 && index_122!=0 ) break;
 443         }
 444
 445         // The mother daughter information for the hard subprocess entry (120)
 446         // IS correct, whereas the information for the particles participating
 447         // in the hard subprocess contains instead the color flow relationships
 448         // Transfer the hard subprocess info onto the other particles
 449         // in the hard subprocess.
 450         //
 451         // We cannot specify daughters of the incoming hard process particles
 452         // because they have some daughters (their showered versions) which
 453         // are not adjacent in the particle record, so we cannot properly
 454         // set the daughter indices in hepevt.
 455         //
 456         if (index_121) HEPEVT_Wrapper::set_parents(index_121, index_101, 0 );
 457         if (index_121) HEPEVT_Wrapper::set_children( index_121, 0, 0 );
 458         if (index_122) HEPEVT_Wrapper::set_parents(index_122, index_102, 0 );
 459         if (index_122) HEPEVT_Wrapper::set_children( index_122, 0, 0 );
 460
 461         for ( int i = HEPEVT_Wrapper::first_child(index_hard);
 462               i <= HEPEVT_Wrapper::last_child(index_hard); i++ ) {
 463             //BPK - Atlas ->
 464             if (index_hard && HEPEVT_Wrapper::id(index_hard) == 0 ) {
 465               HEPEVT_Wrapper::set_parents(
 466                   i, HEPEVT_Wrapper::first_parent(index_hard),
 467                   HEPEVT_Wrapper::last_parent(index_hard) );
 468             //BPK -> inconsistency in HWHGUP, desc from hard vert should point to it.
 469             } else if (  HEPEVT_Wrapper::first_parent(i)!=index_hard) {
 470               HEPEVT_Wrapper::set_parents(i,index_hard,HEPEVT_Wrapper::last_parent(i) );
 471             }
 472             //BPK - Atlas -<
 473
 474             // When the direct descendants of the hard process are hadrons,
 475             // then the 2nd child contains color flow information, and so
 476             // we zero it.
 477             // However, if the direct descendant is status=195, then it is
 478             // a non-hadron, and so the 2nd child does contain real mother
 479             // daughter relationships. ( particularly relevant for H->WW,
 480             //                           April 18, 2003 )
 481             // BPK - part of the inconsistency in HWHGUP problem
 482             if ( HEPEVT_Wrapper::status(i) != 195 && HEPEVT_Wrapper::status(i) != 155 ) {
 483               HEPEVT_Wrapper::set_children(i,HEPEVT_Wrapper::first_child(i),0);
 484             }
 485         }
 486
 487         // now zero the collision and hard entries.
 488         //BPK - Atlas ->
 489         if (index_hard && HEPEVT_Wrapper::id(index_hard) == 0 ) zero_hepevt_entry(index_hard);
 490         if (index_hard && HEPEVT_Wrapper::id(index_collision) == 0  ) zero_hepevt_entry(index_collision);
 491         //BPK - Atlas -<
 492
 493         //     Loop over the particles individually and handle oddities
 494         for ( int i=1; i <=HEPEVT_Wrapper::number_entries(); i++ ) {
 495
 496             //       ----------- Fix ID codes ----------
 497             //       particles with ID=94 are mirror images of their mothers:
 498             if ( HEPEVT_Wrapper::id(i)==94 ) {
 499                 HEPEVT_Wrapper::set_id(
 500                     i, HEPEVT_Wrapper::id( HEPEVT_Wrapper::first_parent(i) ) );
 501             }
 502
 503             //     ----------- fix STATUS codes ------
 504             //     status=100 particles are "cones" which carry only color info
 505             //     throw them away
 506             if ( HEPEVT_Wrapper::status(i)==100 ) zero_hepevt_entry(i);
 507
 508
 509             // NOTE: status 101,102 particles are the beam particles.
 510             //       status 121,129 particles are the hard subprocess particles
 511             // we choose to allow the herwig particles to have herwig
 512             // specific codes, and so we don't bother to change these
 513             // to status =3.
 514
 515
 516
 517
 518             //  ----------- fix some MOTHER/DAUGHTER relationships
 519             //  Whenever the mother points to the hard process, it is referring
 520             //  to a color flow, so we zero it.
 521             if ( HEPEVT_Wrapper::last_parent(i)==index_hard ) {
 522                 HEPEVT_Wrapper::set_parents(
 523                     i, HEPEVT_Wrapper::first_parent(i), 0 );
 524             }
 525
 526             // It makes no sense to have a mother that is younger than you are!
 527
 528             if ( HEPEVT_Wrapper::first_parent(i) >= i ) {
 529                 HEPEVT_Wrapper::set_parents( i, 0, 0 );
 530             }
 531             if ( HEPEVT_Wrapper::last_parent(i) >= i ) {
 532                 HEPEVT_Wrapper::set_parents(
 533                     i, HEPEVT_Wrapper::first_parent(i), 0 );
 534             }
 535
 536             // Whenever the second mother/daughter has a lower index than the
 537             // first, it means the second mother/daughter contains color
 538             // info. Purge it.
 539             if ( HEPEVT_Wrapper::last_parent(i) <=
 540                  HEPEVT_Wrapper::first_parent(i) ) {
 541                 HEPEVT_Wrapper::set_parents(
 542                     i, HEPEVT_Wrapper::first_parent(i), 0 );
 543             }
 544
 545             if ( HEPEVT_Wrapper::last_child(i) <=
 546                  HEPEVT_Wrapper::first_child(i) ) {
 547                 HEPEVT_Wrapper::set_children(
 548                     i, HEPEVT_Wrapper::first_child(i), 0 );
 549             }
 550
 551             // The mothers & daughters of a soft centre of mass (stat=170) seem
 552             // to be correct, but they are out of sequence. The information is
 553             // elsewhere in the event record, so zero it.
 554             //
 555             if ( HEPEVT_Wrapper::status(i) == 170 ) {
 556                 HEPEVT_Wrapper::set_parents( i, 0, 0 );
 557                 HEPEVT_Wrapper::set_children( i, 0, 0 );
 558             }
 559
 560             // Recognise clusters.
 561             // Case 1: cluster has particle parents.
 562             // Clusters normally DO point to its two
 563             // correct mothers, but those 2 mothers are rarely adjacent in the
 564             // event record ... so the mother information might say something
 565             // like 123,48 where index123 and index48 really are the correct
 566             // mothers... however the hepevt standard states that the mother
 567             // pointers should give the index range. So we would have to
 568             // reorder the event record and add entries if we wanted to use
 569             // it. Instead we just zero the mothers, since all of that
 570             // information is contained in the daughter information of the
 571             // mothers.
 572             // Case 2: cluster has a soft process centre of mass (stat=170)
 573             // as parent. This is ok, keep it.
 574             //
 575             // Note if we were going directly to HepMC, then we could
 576             //  use this information properly!
 577
 578             if ( HEPEVT_Wrapper::id(i)==91 ) {
 579                 // if the cluster comes from a SOFT (id=0,stat=170)
 580                 if ( HEPEVT_Wrapper::status(HEPEVT_Wrapper::first_parent(i))
 581                      == 170 ) {
 582                     ; // In this case the mothers are ok
 583                 } else {
 584                     HEPEVT_Wrapper::set_parents( i, 0, 0 );
 585                 }
 586             }
 587         }
 588
 589         //     ---------- Loop over the particles individually and look
 590         //                for mother/daughter inconsistencies.
 591         // We consider a mother daughter relationship to be valid
 592         // ONLy when the mother points to the daughter AND the
 593         // daughter points back (true valid bidirectional
 594         // pointers) OR when a one thing points to the other, but
 595         // the other points to zero. If this isn't true, we zero
 596         // the offending relationship.
 597
 598         for ( int i=1; i <=HEPEVT_Wrapper::number_entries(); i++ ) {
 599             // loop over parents
 600             int ifirst = HEPEVT_Wrapper::first_parent(i);
 601             int ilast = HEPEVT_Wrapper::last_parent(i);
 602             if ( ilast == 0 ) ilast = HEPEVT_Wrapper::first_parent(i);
 603             bool first_is_acceptable = true;
 604             bool last_is_acceptable = true;
 605             // check for out of range.
 606             if ( ifirst>=i || ifirst<0 ) first_is_acceptable = false;
 607             if ( ilast>=i || ilast<ifirst || ilast<0 )last_is_acceptable=false;
 608             if ( first_is_acceptable ) {
 609                 for ( int j = ifirst; j<=ilast; j++ ) {
 610                     // these are the acceptable outcomes
 611                     if ( HEPEVT_Wrapper::first_child(j)==i ) {;}
 612                     // watch out
 613                     else if ( HEPEVT_Wrapper::first_child(j) <=i &&
 614                               HEPEVT_Wrapper::last_child(j) >=i ) {;}
 615                     else if ( HEPEVT_Wrapper::first_child(j) ==0 &&
 616                               HEPEVT_Wrapper::last_child(j) ==0 ) {;}
 617
 618                     // Error Condition:
 619                     // modified by MADobbs@lbl.gov April 21, 2003
 620                     // we distinguish between the first parent and all parents
 621                     //  being incorrect
 622                     else if (j==ifirst) { first_is_acceptable = false; break; }
 623                     else { last_is_acceptable = false; break; }
 624                 }
 625             }
 626             // if any one of the mothers gave a bad outcome, zero all mothers
 627             //BPK - Atlas ->
 628             // do not disconnect photons (most probably from photos)
 629             if ( HEPEVT_Wrapper::id(i) == 22 && HEPEVT_Wrapper::status(i) == 1 )
 630               { first_is_acceptable = true; }
 631             //BPK - Atlas -<
 632             if ( !first_is_acceptable ) {
 633               HEPEVT_Wrapper::set_parents( i, 0, 0 );
 634             } else if ( !last_is_acceptable ) {
 635               HEPEVT_Wrapper::set_parents(i,HEPEVT_Wrapper::first_parent(i),0);
 636             }
 637         }
 638         // Note: it's important to finish the mother loop, before
 639         // starting the daughter loop ... since many mother relations
 640         // will be zero'd which will validate the daughters.... i.e.,
 641         // we want relationships like:
 642         //      IHEP    ID      IDPDG IST MO1 MO2 DA1 DA2
 643         //        27 TQRK           6   3  26  26  30  30
 644         //        30 TQRK           6 155  26  11  31  32
 645         // to come out right.
 646
 647         for ( int i=1; i <=HEPEVT_Wrapper::number_entries(); i++ ) {
 648             // loop over daughters
 649             int ifirst = HEPEVT_Wrapper::first_child(i);
 650             int ilast = HEPEVT_Wrapper::last_child(i);
 651             if ( ilast==0 ) ilast = HEPEVT_Wrapper::first_child(i);
 652             bool is_acceptable = true;
 653             // check for out of range.
 654             if ( ifirst<=i || ifirst<0 ) is_acceptable = false;
 655             if ( ilast<=i || ilast<ifirst || ilast<0 ) is_acceptable = false;
 656             if ( is_acceptable ) {
 657                 for ( int j = ifirst; j<=ilast; j++ ) {
 658                     // these are the acceptable outcomes
 659                     if ( HEPEVT_Wrapper::first_parent(j)==i ) {;}
 660                     else if ( HEPEVT_Wrapper::first_parent(j) <=i &&
 661                               HEPEVT_Wrapper::last_parent(j) >=i ) {;}
 662                     else if ( HEPEVT_Wrapper::first_parent(j) ==0 &&
 663                               HEPEVT_Wrapper::last_parent(j) ==0 ) {;}
 664                     else { is_acceptable = false; } // error condition
 665                 }
 666             }
 667             // if any one of the children gave a bad outcome, zero all children
 668             if ( !is_acceptable ) HEPEVT_Wrapper::set_children( i, 0, 0 );
 669         }
 670
 671         // fixme
 672
 673         for ( int i=1; i <=HEPEVT_Wrapper::number_entries(); i++ ) {
 674             HEPEVT_Wrapper::set_id(
 675                 i, translate_herwig_to_pdg_id(HEPEVT_Wrapper::id(i)) );
 676         }
 677
 678
 679         if ( m_no_gaps_in_barcodes ) remove_gaps_in_hepevt();
 680     }
 681
 682     void IO_HERWIG::remove_gaps_in_hepevt() const {
 683         /// in this scenario, we do not allow there to be zero-ed
 684         /// entries in the HEPEVT common block, and so be reshuffle
 685         /// the common block, removing the zeero-ed entries as we
 686         /// go and making sure we keep the mother/daughter
 687         /// relationships appropriate
 688         std::vector<int> mymap(HEPEVT_Wrapper::number_entries()+1,0);
 689         int ilast = 0;
 690         for ( int i=1; i <=HEPEVT_Wrapper::number_entries(); i++ ) {
 691             if (HEPEVT_Wrapper::status(i)==0 && HEPEVT_Wrapper::id(i)==0) {
 692                 // we remove all entries for which stat=0, id=0
 693                 mymap[i]=0;
 694             } else {
 695                 ilast += 1;
 696                 if ( ilast != i ) {
 697                     HEPEVT_Wrapper::set_status(ilast,
 698                                                HEPEVT_Wrapper::status(i) );
 699                     HEPEVT_Wrapper::set_id(ilast, HEPEVT_Wrapper::id(i) );
 700                     HEPEVT_Wrapper::set_parents(
 701                         ilast,
 702                         HEPEVT_Wrapper::first_parent(i),
 703                         HEPEVT_Wrapper::last_parent(i) );
 704                     HEPEVT_Wrapper::set_children(
 705                         ilast,
 706                         HEPEVT_Wrapper::first_child(i),
 707                         HEPEVT_Wrapper::last_child(i) );
 708                     HEPEVT_Wrapper::set_momentum(
 709                         ilast,
 710                         HEPEVT_Wrapper::px(i), HEPEVT_Wrapper::py(i),
 711                         HEPEVT_Wrapper::pz(i), HEPEVT_Wrapper::e(i)  );
 712                     HEPEVT_Wrapper::set_mass(ilast, HEPEVT_Wrapper::m(i) );
 713                     HEPEVT_Wrapper::set_position(
 714                         ilast, HEPEVT_Wrapper::x(i),HEPEVT_Wrapper::y(i),
 715                         HEPEVT_Wrapper::z(i),HEPEVT_Wrapper::t(i) );
 716                 }
 717                 mymap[i]=ilast;
 718             }
 719         }
 720
 721         // M. Dobbs (from Borut) - April 26, to fix tauolo/herwig past
 722         // the end problem with daughter pointers:
 723         // HEPEVT_Wrapper::set_number_entries( ilast );
 724
 725         // Finally we need to re-map the mother/daughter pointers.
 726         for ( int i=1; i <=ilast; i++ ) {
 727
 728             HEPEVT_Wrapper::set_parents(
 729                 i,
 730                 mymap[HEPEVT_Wrapper::first_parent(i)],
 731                 mymap[HEPEVT_Wrapper::last_parent(i)] );
 732             HEPEVT_Wrapper::set_children(
 733                 i,
 734                 mymap[HEPEVT_Wrapper::first_child(i)],
 735                 mymap[HEPEVT_Wrapper::last_child(i)] );
 736         }
 737         // M. Dobbs (from Borut, part B) - April 26, to fix tauolo/herwig past
 738         // the end problem with daughter pointers:
 739         HEPEVT_Wrapper::set_number_entries( ilast );
 740     }
 741
 742     void IO_HERWIG::zero_hepevt_entry( int i ) const {
 743       if ( i <=0 || i > HepMC::HEPEVT_Wrapper::max_number_entries() ) return;
 744       HEPEVT_Wrapper::set_status( i, 0 );
 745       HEPEVT_Wrapper::set_id( i, 0 );
 746       HEPEVT_Wrapper::set_parents( i, 0, 0 );
 747       HEPEVT_Wrapper::set_children( i, 0, 0 );
 748       HEPEVT_Wrapper::set_momentum( i, 0, 0, 0, 0 );
 749       HEPEVT_Wrapper::set_mass( i, 0 );
 750       HEPEVT_Wrapper::set_position( i, 0, 0, 0, 0 );
 751     }
 752
 753     int IO_HERWIG::translate_herwig_to_pdg_id( int id ) const {
 754         /// This routine is copied from Lynn Garren's stdhep 5.01.
 755         ///   see http:///cepa.fnal.gov/psm/stdhep/
 756
 757                                                // example -9922212
 758         int hwtran = id;                       //         -9922212
 759         int ida    = abs(id);                  //          9922212
 760         int j1     = ida%10;                   //                2
 761         int i1     = (ida/10)%10;              //               1
 762         int i2     = (ida/100)%10;             //              2
 763         int i3     = (ida/1000)%10;            //             2
 764         //int i4     =(ida/10000)%10;          //            2
 765         //int i5     =(ida/100000)%10;         //           9
 766         //int k99    = (ida/100000)%100;       //          9
 767         int ksusy  = (ida/1000000)%10;         //         0
 768         //int ku     = (ida/10000000)%10;      //        0
 769         int kqn    = (ida/1000000000)%10;      //       0
 770
 771         if ( kqn==1 ) {
 772             //  ions not recognized
 773             hwtran=0;
 774             if ( m_print_inconsistency_errors ) {
 775                 std::cerr << "IO_HERWIG::translate_herwig_to_pdg_id " << id
 776                           << "nonallowed ion" << std::endl;
 777             }
 778         }
 779         else if (ida < 100) {
 780             // Higgs, etc.
 781             hwtran = m_herwig_to_pdg_id[ida];
 782             if ( id < 0 ) hwtran *= -1;
 783             // check for illegal antiparticles
 784             if ( id < 0 ) {
 785                 if ( hwtran>=-99 && hwtran<=-81) hwtran=0;
 786                 if ( m_no_antiparticles.count(hwtran) ) hwtran=0;
 787             }
 788         }
 789         else if ( ksusy==1 || ksusy==2 ) { ; }
 790         //  SUSY
 791         else if ( i1!=0 && i3!=0 && j1==2 ) {;}
 792         // spin 1/2 baryons
 793         else if ( i1!=0 && i3!=0 && j1==4 ) {;}
 794         // spin 3/2 baryons
 795         else if ( i1!=0 && i2!=0 && i3==0 ) {
 796             // mesons
 797             // check for illegal antiparticles
 798             if ( i1==i2 && id<0) hwtran=0;
 799         }
 800         else if ( i2!=0 && i3!=0 && i1==0 ) {;}
 801         // diquarks
 802         else {
 803             // undefined
 804             hwtran=0;
 805         }
 806
 807         // check for illegal anti KS, KL
 808         if ( id==-130 || id==-310 ) hwtran=0;
 809
 810         if ( hwtran==0 && ida!=0 && m_print_inconsistency_errors ) {
 811             std::cerr
 812                 << "IO_HERWIG::translate_herwig_to_pdg_id HERWIG particle "
 813                 << id << " translates to zero." << std::endl;
 814         }
 815
 816         return hwtran;
 817     }
 818
 819 } // HepMC
 820
 821
 822
 823