hepmc - Blame information for rev 306

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