Line data Source code
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 0 : IO_HERWIG::IO_HERWIG() : m_trust_mothers_before_daughters(false),
13 0 : m_trust_both_mothers_and_daughters(true),
14 0 : m_print_inconsistency_errors(true),
15 0 : m_no_gaps_in_barcodes(true),
16 0 : m_herwig_to_pdg_id(100,0)
17 0 : {
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 0 : m_herwig_to_pdg_id[1] =1;
22 0 : m_herwig_to_pdg_id[2] =2;
23 0 : m_herwig_to_pdg_id[3] =3;
24 0 : m_herwig_to_pdg_id[4] =4;
25 0 : m_herwig_to_pdg_id[5] =5;
26 0 : m_herwig_to_pdg_id[6] =6;
27 0 : m_herwig_to_pdg_id[7] =7;
28 0 : m_herwig_to_pdg_id[8] =8;
29 :
30 0 : m_herwig_to_pdg_id[11] =11;
31 0 : m_herwig_to_pdg_id[12] =12;
32 0 : m_herwig_to_pdg_id[13] =13;
33 0 : m_herwig_to_pdg_id[14] =14;
34 0 : m_herwig_to_pdg_id[15] =15;
35 0 : m_herwig_to_pdg_id[16] =16;
36 :
37 0 : m_herwig_to_pdg_id[21] =21;
38 0 : m_herwig_to_pdg_id[22] =22;
39 0 : m_herwig_to_pdg_id[23] =23;
40 0 : m_herwig_to_pdg_id[24] =24;
41 0 : m_herwig_to_pdg_id[25] =25;
42 0 : m_herwig_to_pdg_id[26] =51; // <-- H_L0 (redundant with h0(25))
43 :
44 0 : m_herwig_to_pdg_id[32] =32;
45 0 : m_herwig_to_pdg_id[35] =35;
46 0 : m_herwig_to_pdg_id[36] =36;
47 0 : m_herwig_to_pdg_id[37] =37;
48 0 : m_herwig_to_pdg_id[39] =39;
49 :
50 0 : m_herwig_to_pdg_id[40] =40; //Charybdis Black Hole
51 :
52 0 : m_herwig_to_pdg_id[81] =81;
53 0 : m_herwig_to_pdg_id[82] =82;
54 0 : m_herwig_to_pdg_id[83] =83;
55 0 : m_herwig_to_pdg_id[84] =84;
56 0 : m_herwig_to_pdg_id[85] =85;
57 0 : m_herwig_to_pdg_id[86] =86;
58 0 : m_herwig_to_pdg_id[87] =87;
59 0 : m_herwig_to_pdg_id[88] =88;
60 0 : m_herwig_to_pdg_id[89] =89;
61 0 : m_herwig_to_pdg_id[90] =90;
62 :
63 0 : m_herwig_to_pdg_id[91] =91;
64 0 : m_herwig_to_pdg_id[92] =92;
65 0 : m_herwig_to_pdg_id[93] =93;
66 0 : m_herwig_to_pdg_id[94] =94;
67 0 : m_herwig_to_pdg_id[95] =95;
68 0 : m_herwig_to_pdg_id[96] =96;
69 0 : m_herwig_to_pdg_id[97] =97;
70 0 : m_herwig_to_pdg_id[98] =9920022; // <-- remnant photon
71 0 : m_herwig_to_pdg_id[99] =9922212; // <-- remnant nucleon
72 :
73 : // These particle ID's have no antiparticle, so aren't allowed.
74 0 : m_no_antiparticles.insert(-21);
75 0 : m_no_antiparticles.insert(-22);
76 0 : m_no_antiparticles.insert(-23);
77 0 : m_no_antiparticles.insert(-25);
78 0 : m_no_antiparticles.insert(-51);
79 0 : m_no_antiparticles.insert(-35);
80 0 : m_no_antiparticles.insert(-36);
81 0 : }
82 :
83 0 : IO_HERWIG::~IO_HERWIG(){}
84 :
85 : void IO_HERWIG::print( std::ostream& ostr ) const {
86 0 : ostr << "IO_HERWIG: reads an event from the FORTRAN Herwig HEPEVT "
87 0 : << "common block. \n"
88 0 : << " trust_mothers_before_daughters = "
89 0 : << m_trust_mothers_before_daughters
90 0 : << " trust_both_mothers_and_daughters = "
91 0 : << m_trust_both_mothers_and_daughters
92 0 : << " print_inconsistency_errors = "
93 0 : << m_print_inconsistency_errors << std::endl;
94 0 : }
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 0 : if ( !evt ) {
102 : std::cerr
103 0 : << "IO_HERWIG::fill_next_event error - passed null event."
104 0 : << std::endl;
105 0 : return false;
106 : }
107 :
108 : // 1. First we have to fix the HEPEVT input, which is all mucked up for
109 : // herwig.
110 0 : repair_hepevt();
111 :
112 0 : 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 0 : std::vector<GenParticle*> hepevt_particle(
124 0 : HEPEVT_Wrapper::number_entries()+1 );
125 0 : hepevt_particle[0] = 0;
126 0 : for ( int i1 = 1; i1 <= HEPEVT_Wrapper::number_entries(); ++i1 ) {
127 0 : hepevt_particle[i1] = build_particle(i1);
128 : }
129 0 : 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 0 : 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 0 : for ( int i = 1; i <=HEPEVT_Wrapper::number_entries(); i++ ) {
147 0 : if ( HEPEVT_Wrapper::status(i)==121 ) index_121=i;
148 0 : if ( HEPEVT_Wrapper::status(i)==122 ) index_122=i;
149 0 : if ( index_121!=0 && index_122!=0 ) break;
150 : }
151 0 : if ( index_121 && index_122 ) {
152 0 : GenVertex* hard_vtx = new GenVertex();
153 0 : hard_vtx->add_particle_in( hepevt_particle[index_121] );
154 0 : 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 0 : for ( int i = 1; i <=HEPEVT_Wrapper::number_entries(); i++ ) {
160 0 : if ( HEPEVT_Wrapper::status(i)==120 ) index_hard=i;
161 0 : if ( index_hard!=0 ) break;
162 : }
163 :
164 0 : if ( index_hard!=0) {
165 0 : hard_vtx->add_particle_out( hepevt_particle[index_hard] );
166 0 : GenVertex* hard_vtx2 = new GenVertex();
167 0 : hard_vtx2->add_particle_in( hepevt_particle[index_hard] );
168 0 : for ( int i = 1; i <= HEPEVT_Wrapper::number_entries(); ++i ) {
169 0 : if ( HEPEVT_Wrapper::first_parent(i)==index_hard ) {
170 0 : hard_vtx2->add_particle_out( hepevt_particle[i] );
171 : }
172 : }
173 0 : evt->set_signal_process_vertex( hard_vtx );
174 0 : evt->set_signal_process_vertex( hard_vtx2 );
175 0 : }
176 : else {
177 0 : evt->set_signal_process_vertex( hard_vtx );
178 : }
179 : //BPK - Atlas -<
180 0 : }
181 : //
182 : // 4. loop over HEPEVT particles AGAIN, this time creating vertices
183 0 : 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 0 : if ( m_trust_mothers_before_daughters ||
192 0 : m_trust_both_mothers_and_daughters ) {
193 0 : 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 0 : if ( !m_trust_mothers_before_daughters ||
199 0 : m_trust_both_mothers_and_daughters ) {
200 0 : 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 0 : 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 0 : if ( hepevt_particle[i3] && !hepevt_particle[i3]->parent_event()
214 0 : && !hepevt_particle[i3]->pdg_id()
215 0 : && !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 0 : delete hepevt_particle[i3];
220 0 : } else if ( hepevt_particle[i3] &&
221 0 : !hepevt_particle[i3]->end_vertex() &&
222 0 : !hepevt_particle[i3]->production_vertex() ) {
223 0 : GenVertex* prod_vtx = new GenVertex();
224 0 : prod_vtx->add_particle_out( hepevt_particle[i3] );
225 0 : evt->add_vertex( prod_vtx );
226 0 : }
227 : }
228 : return true;
229 0 : }
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 0 : GenParticle* p = hepevt_particle[i];
239 : // a. search to see if a production vertex already exists
240 0 : int mother = HEPEVT_Wrapper::first_parent(i);
241 0 : GenVertex* prod_vtx = p->production_vertex();
242 0 : while ( !prod_vtx && mother > 0 ) {
243 0 : prod_vtx = hepevt_particle[mother]->end_vertex();
244 0 : if ( prod_vtx ) prod_vtx->add_particle_out( p );
245 : // increment mother for next iteration
246 0 : 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 0 : FourVector prod_pos( HEPEVT_Wrapper::x(i), HEPEVT_Wrapper::y(i),
252 0 : HEPEVT_Wrapper::z(i), HEPEVT_Wrapper::t(i)
253 : );
254 0 : if ( !prod_vtx && (HEPEVT_Wrapper::number_parents(i)>0
255 0 : || prod_pos!=FourVector(0,0,0,0)) )
256 : {
257 0 : prod_vtx = new GenVertex();
258 0 : prod_vtx->add_particle_out( p );
259 0 : evt->add_vertex( prod_vtx );
260 0 : }
261 : // c. if prod_vtx doesn't already have position specified, fill it
262 0 : if ( prod_vtx && prod_vtx->position()==FourVector(0,0,0,0) ) {
263 0 : prod_vtx->set_position( prod_pos );
264 0 : }
265 : // d. loop over mothers to make sure their end_vertices are
266 : // consistent
267 0 : mother = HEPEVT_Wrapper::first_parent(i);
268 0 : while ( prod_vtx && mother > 0 ) {
269 0 : if ( !hepevt_particle[mother]->end_vertex() ) {
270 : // if end vertex of the mother isn't specified, do it now
271 0 : prod_vtx->add_particle_in( hepevt_particle[mother] );
272 0 : } 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 0 : if ( m_print_inconsistency_errors ) {
282 : std::cerr
283 0 : << "HepMC::IO_HERWIG: inconsistent mother/daugher "
284 0 : << "information in HEPEVT event "
285 0 : << HEPEVT_Wrapper::event_number()
286 0 : << ". \n I recommend you try "
287 0 : << "inspecting the event first with "
288 0 : << "\n\tHEPEVT_Wrapper::check_hepevt_consistency()"
289 0 : << "\n This warning can be turned off with the "
290 0 : << "IO_HERWIG::print_inconsistency_errors switch."
291 0 : << std::endl;
292 0 : hepevt_particle[mother]->print(std::cerr);
293 : std::cerr
294 0 : << "problem vertices are: (prod_vtx, mother)" << std::endl;
295 0 : if ( prod_vtx ) prod_vtx->print(std::cerr);
296 0 : hepevt_particle[mother]->end_vertex()->print(std::cerr);
297 0 : }
298 : }
299 0 : if ( ++mother > HEPEVT_Wrapper::last_parent(i) ) mother = 0;
300 : }
301 0 : }
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 0 : GenParticle* p = hepevt_particle[i];
311 : // a.
312 0 : int daughter = HEPEVT_Wrapper::first_child(i);
313 0 : GenVertex* end_vtx = p->end_vertex();
314 0 : while ( !end_vtx && daughter > 0 ) {
315 0 : end_vtx = hepevt_particle[daughter]->production_vertex();
316 0 : if ( end_vtx ) end_vtx->add_particle_in( p );
317 0 : 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 0 : if ( !end_vtx && HEPEVT_Wrapper::number_children(i)>0 ) {
322 0 : end_vtx = new GenVertex();
323 0 : end_vtx->add_particle_in( p );
324 0 : evt->add_vertex( end_vtx );
325 0 : }
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 0 : daughter = HEPEVT_Wrapper::first_child(i);
330 0 : while ( end_vtx && daughter > 0 ) {
331 0 : if ( !hepevt_particle[daughter]->production_vertex() ) {
332 : // if end vertex of the mother isn't specified, do it now
333 0 : end_vtx->add_particle_out( hepevt_particle[daughter] );
334 : //
335 : // 2001-03-29 M.Dobbs, fill vertex the position.
336 0 : if ( end_vtx->position()==FourVector(0,0,0,0) ) {
337 0 : FourVector prod_pos( HEPEVT_Wrapper::x(daughter),
338 0 : HEPEVT_Wrapper::y(daughter),
339 0 : HEPEVT_Wrapper::z(daughter),
340 0 : HEPEVT_Wrapper::t(daughter)
341 : );
342 0 : if ( prod_pos != FourVector(0,0,0,0) ) {
343 0 : end_vtx->set_position( prod_pos );
344 0 : }
345 0 : }
346 0 : } else if (hepevt_particle[daughter]->production_vertex()
347 0 : != 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 0 : if ( m_print_inconsistency_errors ) std::cerr
354 0 : << "HepMC::IO_HERWIG: inconsistent mother/daugher "
355 0 : << "information in HEPEVT event "
356 0 : << HEPEVT_Wrapper::event_number()
357 0 : << ". \n I recommend you try "
358 0 : << "inspecting the event first with "
359 0 : << "\n\tHEPEVT_Wrapper::check_hepevt_consistency()"
360 0 : << "\n This warning can be turned off with the "
361 0 : << "IO_HERWIG::print_inconsistency_errors switch."
362 0 : << std::endl;
363 : }
364 0 : if ( ++daughter > HEPEVT_Wrapper::last_child(i) ) daughter = 0;
365 : }
366 0 : if ( !p->end_vertex() && !p->production_vertex() ) {
367 : // Added 2001-11-04, to try and handle Isajet problems.
368 0 : build_production_vertex( i, hepevt_particle, evt );
369 0 : }
370 0 : }
371 :
372 : GenParticle* IO_HERWIG::build_particle( int index ) {
373 : /// Builds a particle object corresponding to index in HEPEVT
374 : //
375 : GenParticle* p
376 0 : = new GenParticle( FourVector( HEPEVT_Wrapper::px(index),
377 0 : HEPEVT_Wrapper::py(index),
378 0 : HEPEVT_Wrapper::pz(index),
379 0 : HEPEVT_Wrapper::e(index) ),
380 0 : HEPEVT_Wrapper::id(index),
381 0 : HEPEVT_Wrapper::status(index) );
382 0 : p->setGeneratedMass( HEPEVT_Wrapper::m(index) );
383 0 : p->suggest_barcode( index );
384 0 : return p;
385 0 : }
386 :
387 : int IO_HERWIG::find_in_map( const std::map<GenParticle*,int>& m,
388 : GenParticle* p) const {
389 0 : std::map<GenParticle*,int>::const_iterator iter = m.find(p);
390 0 : if ( iter == m.end() ) return 0;
391 0 : return iter->second;
392 0 : }
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 0 : 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 0 : for ( int i = 1; i <=HEPEVT_Wrapper::number_entries(); i++ ) {
435 0 : if ( HEPEVT_Wrapper::status(i)==101 ) index_101=i;
436 0 : if ( HEPEVT_Wrapper::status(i)==102 ) index_102=i;
437 0 : if ( HEPEVT_Wrapper::status(i)==103 ) index_collision=i;
438 0 : if ( HEPEVT_Wrapper::status(i)==120 ) index_hard=i;
439 0 : if ( HEPEVT_Wrapper::status(i)==121 ) index_121=i;
440 0 : if ( HEPEVT_Wrapper::status(i)==122 ) index_122=i;
441 0 : if ( index_collision!=0 && index_hard!=0 && index_101!=0 &&
442 0 : 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 0 : if (index_121) HEPEVT_Wrapper::set_parents(index_121, index_101, 0 );
457 0 : if (index_121) HEPEVT_Wrapper::set_children( index_121, 0, 0 );
458 0 : if (index_122) HEPEVT_Wrapper::set_parents(index_122, index_102, 0 );
459 0 : if (index_122) HEPEVT_Wrapper::set_children( index_122, 0, 0 );
460 :
461 0 : for ( int i = HEPEVT_Wrapper::first_child(index_hard);
462 0 : i <= HEPEVT_Wrapper::last_child(index_hard); i++ ) {
463 : //BPK - Atlas ->
464 0 : if (index_hard && HEPEVT_Wrapper::id(index_hard) == 0 ) {
465 0 : HEPEVT_Wrapper::set_parents(
466 0 : i, HEPEVT_Wrapper::first_parent(index_hard),
467 0 : HEPEVT_Wrapper::last_parent(index_hard) );
468 : //BPK -> inconsistency in HWHGUP, desc from hard vert should point to it.
469 0 : } else if ( HEPEVT_Wrapper::first_parent(i)!=index_hard) {
470 0 : HEPEVT_Wrapper::set_parents(i,index_hard,HEPEVT_Wrapper::last_parent(i) );
471 0 : }
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 0 : if ( HEPEVT_Wrapper::status(i) != 195 && HEPEVT_Wrapper::status(i) != 155 ) {
483 0 : HEPEVT_Wrapper::set_children(i,HEPEVT_Wrapper::first_child(i),0);
484 0 : }
485 : }
486 :
487 : // now zero the collision and hard entries.
488 : //BPK - Atlas ->
489 0 : if (index_hard && HEPEVT_Wrapper::id(index_hard) == 0 ) zero_hepevt_entry(index_hard);
490 0 : 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 0 : 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 0 : if ( HEPEVT_Wrapper::id(i)==94 ) {
499 0 : HEPEVT_Wrapper::set_id(
500 0 : i, HEPEVT_Wrapper::id( HEPEVT_Wrapper::first_parent(i) ) );
501 0 : }
502 :
503 : // ----------- fix STATUS codes ------
504 : // status=100 particles are "cones" which carry only color info
505 : // throw them away
506 0 : 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 0 : if ( HEPEVT_Wrapper::last_parent(i)==index_hard ) {
522 0 : HEPEVT_Wrapper::set_parents(
523 0 : i, HEPEVT_Wrapper::first_parent(i), 0 );
524 0 : }
525 :
526 : // It makes no sense to have a mother that is younger than you are!
527 :
528 0 : if ( HEPEVT_Wrapper::first_parent(i) >= i ) {
529 0 : HEPEVT_Wrapper::set_parents( i, 0, 0 );
530 0 : }
531 0 : if ( HEPEVT_Wrapper::last_parent(i) >= i ) {
532 0 : HEPEVT_Wrapper::set_parents(
533 0 : i, HEPEVT_Wrapper::first_parent(i), 0 );
534 0 : }
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 0 : if ( HEPEVT_Wrapper::last_parent(i) <=
540 0 : HEPEVT_Wrapper::first_parent(i) ) {
541 0 : HEPEVT_Wrapper::set_parents(
542 0 : i, HEPEVT_Wrapper::first_parent(i), 0 );
543 0 : }
544 :
545 0 : if ( HEPEVT_Wrapper::last_child(i) <=
546 0 : HEPEVT_Wrapper::first_child(i) ) {
547 0 : HEPEVT_Wrapper::set_children(
548 0 : i, HEPEVT_Wrapper::first_child(i), 0 );
549 0 : }
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 0 : if ( HEPEVT_Wrapper::status(i) == 170 ) {
556 0 : HEPEVT_Wrapper::set_parents( i, 0, 0 );
557 0 : HEPEVT_Wrapper::set_children( i, 0, 0 );
558 0 : }
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 0 : if ( HEPEVT_Wrapper::id(i)==91 ) {
579 : // if the cluster comes from a SOFT (id=0,stat=170)
580 0 : if ( HEPEVT_Wrapper::status(HEPEVT_Wrapper::first_parent(i))
581 0 : == 170 ) {
582 : ; // In this case the mothers are ok
583 : } else {
584 0 : 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 0 : for ( int i=1; i <=HEPEVT_Wrapper::number_entries(); i++ ) {
599 : // loop over parents
600 0 : int ifirst = HEPEVT_Wrapper::first_parent(i);
601 0 : int ilast = HEPEVT_Wrapper::last_parent(i);
602 0 : 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 0 : if ( ifirst>=i || ifirst<0 ) first_is_acceptable = false;
607 0 : if ( ilast>=i || ilast<ifirst || ilast<0 )last_is_acceptable=false;
608 0 : if ( first_is_acceptable ) {
609 0 : for ( int j = ifirst; j<=ilast; j++ ) {
610 : // these are the acceptable outcomes
611 0 : if ( HEPEVT_Wrapper::first_child(j)==i ) {;}
612 : // watch out
613 0 : else if ( HEPEVT_Wrapper::first_child(j) <=i &&
614 0 : HEPEVT_Wrapper::last_child(j) >=i ) {;}
615 0 : else if ( HEPEVT_Wrapper::first_child(j) ==0 &&
616 0 : 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 0 : else if (j==ifirst) { first_is_acceptable = false; break; }
623 0 : else { last_is_acceptable = false; break; }
624 : }
625 0 : }
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 0 : if ( HEPEVT_Wrapper::id(i) == 22 && HEPEVT_Wrapper::status(i) == 1 )
630 0 : { first_is_acceptable = true; }
631 : //BPK - Atlas -<
632 0 : if ( !first_is_acceptable ) {
633 0 : HEPEVT_Wrapper::set_parents( i, 0, 0 );
634 0 : } else if ( !last_is_acceptable ) {
635 0 : HEPEVT_Wrapper::set_parents(i,HEPEVT_Wrapper::first_parent(i),0);
636 0 : }
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 0 : for ( int i=1; i <=HEPEVT_Wrapper::number_entries(); i++ ) {
648 : // loop over daughters
649 0 : int ifirst = HEPEVT_Wrapper::first_child(i);
650 0 : int ilast = HEPEVT_Wrapper::last_child(i);
651 0 : if ( ilast==0 ) ilast = HEPEVT_Wrapper::first_child(i);
652 : bool is_acceptable = true;
653 : // check for out of range.
654 0 : if ( ifirst<=i || ifirst<0 ) is_acceptable = false;
655 0 : if ( ilast<=i || ilast<ifirst || ilast<0 ) is_acceptable = false;
656 0 : if ( is_acceptable ) {
657 0 : for ( int j = ifirst; j<=ilast; j++ ) {
658 : // these are the acceptable outcomes
659 0 : if ( HEPEVT_Wrapper::first_parent(j)==i ) {;}
660 0 : else if ( HEPEVT_Wrapper::first_parent(j) <=i &&
661 0 : HEPEVT_Wrapper::last_parent(j) >=i ) {;}
662 0 : else if ( HEPEVT_Wrapper::first_parent(j) ==0 &&
663 0 : HEPEVT_Wrapper::last_parent(j) ==0 ) {;}
664 : else { is_acceptable = false; } // error condition
665 : }
666 0 : }
667 : // if any one of the children gave a bad outcome, zero all children
668 0 : if ( !is_acceptable ) HEPEVT_Wrapper::set_children( i, 0, 0 );
669 : }
670 :
671 : // fixme
672 :
673 0 : for ( int i=1; i <=HEPEVT_Wrapper::number_entries(); i++ ) {
674 0 : HEPEVT_Wrapper::set_id(
675 0 : i, translate_herwig_to_pdg_id(HEPEVT_Wrapper::id(i)) );
676 : }
677 :
678 :
679 0 : if ( m_no_gaps_in_barcodes ) remove_gaps_in_hepevt();
680 0 : }
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 0 : std::vector<int> mymap(HEPEVT_Wrapper::number_entries()+1,0);
689 : int ilast = 0;
690 0 : for ( int i=1; i <=HEPEVT_Wrapper::number_entries(); i++ ) {
691 0 : if (HEPEVT_Wrapper::status(i)==0 && HEPEVT_Wrapper::id(i)==0) {
692 : // we remove all entries for which stat=0, id=0
693 0 : mymap[i]=0;
694 0 : } else {
695 0 : ilast += 1;
696 0 : if ( ilast != i ) {
697 0 : HEPEVT_Wrapper::set_status(ilast,
698 0 : HEPEVT_Wrapper::status(i) );
699 0 : HEPEVT_Wrapper::set_id(ilast, HEPEVT_Wrapper::id(i) );
700 0 : HEPEVT_Wrapper::set_parents(
701 : ilast,
702 0 : HEPEVT_Wrapper::first_parent(i),
703 0 : HEPEVT_Wrapper::last_parent(i) );
704 0 : HEPEVT_Wrapper::set_children(
705 : ilast,
706 0 : HEPEVT_Wrapper::first_child(i),
707 0 : HEPEVT_Wrapper::last_child(i) );
708 0 : HEPEVT_Wrapper::set_momentum(
709 : ilast,
710 0 : HEPEVT_Wrapper::px(i), HEPEVT_Wrapper::py(i),
711 0 : HEPEVT_Wrapper::pz(i), HEPEVT_Wrapper::e(i) );
712 0 : HEPEVT_Wrapper::set_mass(ilast, HEPEVT_Wrapper::m(i) );
713 0 : HEPEVT_Wrapper::set_position(
714 0 : ilast, HEPEVT_Wrapper::x(i),HEPEVT_Wrapper::y(i),
715 0 : HEPEVT_Wrapper::z(i),HEPEVT_Wrapper::t(i) );
716 : }
717 0 : 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 0 : for ( int i=1; i <=ilast; i++ ) {
727 :
728 0 : HEPEVT_Wrapper::set_parents(
729 : i,
730 0 : mymap[HEPEVT_Wrapper::first_parent(i)],
731 0 : mymap[HEPEVT_Wrapper::last_parent(i)] );
732 0 : HEPEVT_Wrapper::set_children(
733 : i,
734 0 : mymap[HEPEVT_Wrapper::first_child(i)],
735 0 : 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 0 : HEPEVT_Wrapper::set_number_entries( ilast );
740 0 : }
741 :
742 : void IO_HERWIG::zero_hepevt_entry( int i ) const {
743 0 : if ( i <=0 || i > HepMC::HEPEVT_Wrapper::max_number_entries() ) return;
744 0 : HEPEVT_Wrapper::set_status( i, 0 );
745 0 : HEPEVT_Wrapper::set_id( i, 0 );
746 0 : HEPEVT_Wrapper::set_parents( i, 0, 0 );
747 0 : HEPEVT_Wrapper::set_children( i, 0, 0 );
748 0 : HEPEVT_Wrapper::set_momentum( i, 0, 0, 0, 0 );
749 0 : HEPEVT_Wrapper::set_mass( i, 0 );
750 0 : HEPEVT_Wrapper::set_position( i, 0, 0, 0, 0 );
751 0 : }
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 0 : int hwtran = id; // -9922212
759 0 : int ida = abs(id); // 9922212
760 0 : int j1 = ida%10; // 2
761 0 : int i1 = (ida/10)%10; // 1
762 0 : int i2 = (ida/100)%10; // 2
763 0 : 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 0 : int ksusy = (ida/1000000)%10; // 0
768 : //int ku = (ida/10000000)%10; // 0
769 0 : int kqn = (ida/1000000000)%10; // 0
770 :
771 0 : if ( kqn==1 ) {
772 : // ions not recognized
773 0 : hwtran=0;
774 0 : if ( m_print_inconsistency_errors ) {
775 0 : std::cerr << "IO_HERWIG::translate_herwig_to_pdg_id " << id
776 0 : << "nonallowed ion" << std::endl;
777 0 : }
778 : }
779 0 : else if (ida < 100) {
780 : // Higgs, etc.
781 0 : hwtran = m_herwig_to_pdg_id[ida];
782 0 : if ( id < 0 ) hwtran *= -1;
783 : // check for illegal antiparticles
784 0 : if ( id < 0 ) {
785 0 : if ( hwtran>=-99 && hwtran<=-81) hwtran=0;
786 0 : if ( m_no_antiparticles.count(hwtran) ) hwtran=0;
787 : }
788 : }
789 0 : else if ( ksusy==1 || ksusy==2 ) { ; }
790 : // SUSY
791 0 : else if ( i1!=0 && i3!=0 && j1==2 ) {;}
792 : // spin 1/2 baryons
793 0 : else if ( i1!=0 && i3!=0 && j1==4 ) {;}
794 : // spin 3/2 baryons
795 0 : else if ( i1!=0 && i2!=0 && i3==0 ) {
796 : // mesons
797 : // check for illegal antiparticles
798 0 : if ( i1==i2 && id<0) hwtran=0;
799 : }
800 0 : else if ( i2!=0 && i3!=0 && i1==0 ) {;}
801 : // diquarks
802 : else {
803 : // undefined
804 0 : hwtran=0;
805 : }
806 :
807 : // check for illegal anti KS, KL
808 0 : if ( id==-130 || id==-310 ) hwtran=0;
809 :
810 0 : if ( hwtran==0 && ida!=0 && m_print_inconsistency_errors ) {
811 : std::cerr
812 0 : << "IO_HERWIG::translate_herwig_to_pdg_id HERWIG particle "
813 0 : << id << " translates to zero." << std::endl;
814 0 : }
815 :
816 0 : return hwtran;
817 0 : }
818 :
819 : } // HepMC
820 :
821 :
822 :
823 :
|
