[u/mrichter/AliRoot.git] / TEvtGen / HepMC / HEPEVT_Wrapper.cc

//////////////////////////////////////////////////////////////////////////
// Matt.Dobbs@Cern.CH       June 30, 2000
// Generic Wrapper for the fortran HEPEVT common block
// 
// The static data member's initializations must be separate from .h file.
// 
//////////////////////////////////////////////////////////////////////////

#include "HepMC/HEPEVT_Wrapper.h"

namespace HepMC {

    ////////////////////////////////////////
    // static data member initializations //
    ////////////////////////////////////////

    unsigned int HEPEVT_Wrapper::s_sizeof_int = 4;

    unsigned int HEPEVT_Wrapper::s_sizeof_real = sizeof(double);

    unsigned int HEPEVT_Wrapper::s_max_number_entries = 4000;

    ///////////////////
    // Print Methods //
    ///////////////////

    void HEPEVT_Wrapper::print_hepevt( std::ostream& ostr ) 
    {
	/// dumps the content of this HEPEVT event to ostr  (Width is 80)
	ostr << "________________________________________"
	     << "________________________________________" << std::endl;
	ostr << "***** HEPEVT Common Event#: " 
	     << event_number()
	     << ", " << number_entries() << " particles (max "
	     << max_number_entries() << ") *****";
	if ( is_double_precision() ) {
	    ostr << " Double Precision" << std::endl;
	} else {
	    ostr << " Single Precision" << std::endl;
	}
	ostr << sizeof_int() << "-byte integers, " 
	     << sizeof_real() << "-byte floating point numbers, "
	     << max_number_entries() << "-allocated entries." 
	     << std::endl;
	print_legend(ostr);
	ostr << "________________________________________"
	     << "________________________________________" << std::endl;
	for ( int i=1; i <= number_entries(); ++i ) {
	    print_hepevt_particle( i, ostr );
	}
	ostr << "________________________________________"
	     << "________________________________________" << std::endl;
    }

    void HEPEVT_Wrapper::print_legend( std::ostream& ostr )
    {
	char outline[81];
	sprintf( outline,"%4s %4s %4s %5s   %10s, %9s, %9s, %9s, %10s",
		 "Indx","Stat","Par-","chil-",
		 "(  P_x","P_y","P_z","Energy","M ) ");
	ostr << outline << std::endl;
	sprintf( outline,"%9s %4s %4s    %10s, %9s, %9s, %9s) %9s",
		 "ID ","ents","dren",
		 "Prod (   X","Y","Z","cT", "[mm]");
	ostr << outline << std::endl;
    }	

    void HEPEVT_Wrapper::print_hepevt_particle( int i, std::ostream& ostr ) 
    {
	/// dumps the content HEPEVT particle entry i   (Width is 120)
	/// here i is the C array index (i.e. it starts at 0 ... whereas the
	/// fortran array index starts at 1) So if there's 100 particles, the
	/// last valid index is 100-1=99
	char outline[81];
	sprintf( outline,
		 "%4d %+4d %4d %4d    (%9.3g, %9.3g, %9.3g, %9.3g, %9.3g)"
		 ,i, status(i), first_parent(i), first_child(i),
		 px(i), py(i), pz(i), e(i), m(i) );
	ostr << outline << "\n";
	sprintf( outline,"%+9d %4d %4d    (%9.3g, %9.3g, %9.3g, %9.3g)",
		 // old version was:" (%+9.2e, %+9.2e, %+9.2e, %+9.2e)"
		 id(i), last_parent(i), last_child(i), 
		 x(i), y(i), z(i), t(i) );
	ostr << outline << std::endl;
    }


    bool HEPEVT_Wrapper::check_hepevt_consistency( std::ostream& os )
    {
	/// This method inspects the HEPEVT common block and looks for 
	/// inconsistencies in the mother/daughter pointers
	bool isConsistent=true;
        char header[81];
        sprintf( header,
                 "\n\n\t**** WARNINGInconsistent HEPEVT input, Event %10d ****"
		 , HEPEVT_Wrapper::event_number() );

	for ( int i = 1; i <= HEPEVT_Wrapper::number_entries(); ++i ) {
	    // 1. check its mothers
	    int moth1 = HEPEVT_Wrapper::first_parent( i );
	    int moth2 = HEPEVT_Wrapper::last_parent( i );
	    if ( moth2<moth1 ) {
		if ( isConsistent ) {
		    os << header << std::endl;
		    isConsistent = false;
		    print_legend(os);
		}
		os << "Inconsistent entry " << i 
		   << " first parent > last parent " << std::endl;
		HEPEVT_Wrapper::print_hepevt_particle( i, os );
	    }
	    for ( int m = moth1; m<=moth2 && m!=0; ++m ) {
		if ( m>HEPEVT_Wrapper::number_entries() || m < 0 ) {
		    if ( isConsistent ) {
			os << header << std::endl;
			isConsistent = false;
			print_legend(os);
		    }
		    os << "Inconsistent entry " << i 
		       << " mother points out of range " << std::endl;
		    HEPEVT_Wrapper::print_hepevt_particle( i, os );
		}
		int mChild1 = HEPEVT_Wrapper::first_child(m);
		int mChild2 = HEPEVT_Wrapper::last_child(m);
		// we don't consider null pointers as inconsistent
		if ( mChild1==0 && mChild2==0 ) continue;
		if ( i<mChild1 || i>mChild2 ) {
		    if ( isConsistent ) {
			os << header << std::endl;
			isConsistent = false;
			print_legend(os);
		    }
		    os << "Inconsistent mother-daughter relationship between "
		       << i << " & " << m 
		       << " (try !trust_mother)" << std::endl;
		    HEPEVT_Wrapper::print_hepevt_particle( i, os );
		    HEPEVT_Wrapper::print_hepevt_particle( m, os );
		}
	    }
	    // 2. check its daughters
	    int dau1 = HEPEVT_Wrapper::first_child( i );
	    int dau2 = HEPEVT_Wrapper::last_child( i );
	    if ( dau2<dau1 ) {
		if ( isConsistent ) {
		    os << header << std::endl;
		    isConsistent = false;
		    print_legend(os);
		}
		os << "Inconsistent entry " << i 
		   << " first child > last child " << std::endl;
		HEPEVT_Wrapper::print_hepevt_particle( i, os );
	    }
	    for ( int d = dau1; d<=dau2 && d!=0; ++d ) {
		if ( d>HEPEVT_Wrapper::number_entries() || d < 0 ) {
		    if ( isConsistent ) {
			os << header << std::endl;
			isConsistent = false;
			print_legend(os);
		    }
		    os << "Inconsistent entry " << i 
		       << " child points out of range " << std::endl;
		    HEPEVT_Wrapper::print_hepevt_particle( i, os );
		}
		int d_moth1 = HEPEVT_Wrapper::first_parent(d);
		int d_moth2 = HEPEVT_Wrapper::last_parent(d);
		// we don't consider null pointers as inconsistent
		if ( d_moth1==0 && d_moth2==0 ) continue;
		if ( i<d_moth1 || i>d_moth2 ) {
		    if ( isConsistent ) {
			os << header << std::endl;
			isConsistent = false;
			print_legend(os);
		    }
		    os << "Inconsistent mother-daughter relationship between "
		       << i << " & " << d 
		       << " (try trust_mothers)"<< std::endl;
		    HEPEVT_Wrapper::print_hepevt_particle( i, os );
		    HEPEVT_Wrapper::print_hepevt_particle( d, os );
		}
	    }
	}
	if (!isConsistent) {
	    os << "Above lists all the inconsistencies in the HEPEVT common "
	       << "\n block which has been provided as input to HepMC. "
	       << "\n HepMC WILL have trouble interpreting the mother-daughter"
	       << "\n relationships ... but all other information "
	       << "\n (4-vectors etc) will be correctly transferred."
	       << "\n In order for HepMC to be able to interpret the mother/"
	       << "\n daughter hierachy, it MUST be given consistent input."
	       << "\n This is one of the design criteria of HepMC: "
	       << "\n consistency is enforced by the code.";
	    os << "\nThere is a switch in IO_HEPEVT, set-able using "
	       << "\n IO_HEPEVT::set_trust_mothers_before_daughters( bool )"
	       << "\n which you may want to try.";
	    os << "\nNote: if HEPEVT common block has been filled by pythia"
	       << "\n pyhepc, then the switch MSTP(128)=2 should be used in"
	       << "\n pythia, which instructs pythia not to put multiple "
	       << "\n copies of resonances in the event record.\n";
	    os << "To obtain a file summarizing the inconsistency, you should:"
	       << "\n\t ofstream myFile(\"myInconsistentEvent.txt\"); "
	       << "\n\t HEPEVT_Wrapper::check_hepevt_consistency(myFile); "
	       << "\n\t HEPEVT_Wrapper::print_hepevt(myFile); "
	       << "\n[now write the event to HepMC using something like"
	       << "\n\t\t myIO_HEPEVT->write_event(myEvent); ]"
	       << "\n\t myEvent->print( myFile ); "
	       << "      // print event as HepMC sees it"
	       << "\n ------------------------- Thank-you. \n\n" << std::endl;
	}
	return isConsistent;
    }

    void HEPEVT_Wrapper::zero_everything()
    {
	set_event_number( 0 );
        set_number_entries( 0 );
	for ( int i = 1; i<=max_number_entries(); ++i ) {
	    set_status( i, 0 );
	    set_id( i, 0 );
	    set_parents( i, 0, 0 );
	    set_children( i, 0, 0 );
	    set_momentum( i, 0, 0, 0, 0 );
	    set_mass( i, 0 );
	    set_position( i, 0, 0, 0, 0 );
	}
    }

} // HepMC
Commit	Line	Data
0ca57c2f	1	//////////////////////////////////////////////////////////////////////////
	2	// Matt.Dobbs@Cern.CH June 30, 2000
	3	// Generic Wrapper for the fortran HEPEVT common block
	4	//
	5	// The static data member's initializations must be separate from .h file.
	6	//
	7	//////////////////////////////////////////////////////////////////////////
	8
	9	#include "HepMC/HEPEVT_Wrapper.h"
	10
	11	namespace HepMC {
	12
	13	////////////////////////////////////////
	14	// static data member initializations //
	15	////////////////////////////////////////
	16
	17	unsigned int HEPEVT_Wrapper::s_sizeof_int = 4;
	18
	19	unsigned int HEPEVT_Wrapper::s_sizeof_real = sizeof(double);
	20
	21	unsigned int HEPEVT_Wrapper::s_max_number_entries = 4000;
	22
	23	///////////////////
	24	// Print Methods //
	25	///////////////////
	26
	27	void HEPEVT_Wrapper::print_hepevt( std::ostream& ostr )
	28	{
	29	/// dumps the content of this HEPEVT event to ostr (Width is 80)
	30	ostr << "________________________________________"
	31	<< "________________________________________" << std::endl;
	32	ostr << "***** HEPEVT Common Event#: "
	33	<< event_number()
	34	<< ", " << number_entries() << " particles (max "
	35	<< max_number_entries() << ") *****";
	36	if ( is_double_precision() ) {
	37	ostr << " Double Precision" << std::endl;
	38	} else {
	39	ostr << " Single Precision" << std::endl;
	40	}
	41	ostr << sizeof_int() << "-byte integers, "
	42	<< sizeof_real() << "-byte floating point numbers, "
	43	<< max_number_entries() << "-allocated entries."
	44	<< std::endl;
	45	print_legend(ostr);
	46	ostr << "________________________________________"
	47	<< "________________________________________" << std::endl;
	48	for ( int i=1; i <= number_entries(); ++i ) {
	49	print_hepevt_particle( i, ostr );
	50	}
	51	ostr << "________________________________________"
	52	<< "________________________________________" << std::endl;
	53	}
	54
	55	void HEPEVT_Wrapper::print_legend( std::ostream& ostr )
	56	{
	57	char outline[81];
	58	sprintf( outline,"%4s %4s %4s %5s %10s, %9s, %9s, %9s, %10s",
	59	"Indx","Stat","Par-","chil-",
	60	"( P_x","P_y","P_z","Energy","M ) ");
	61	ostr << outline << std::endl;
	62	sprintf( outline,"%9s %4s %4s %10s, %9s, %9s, %9s) %9s",
	63	"ID ","ents","dren",
	64	"Prod ( X","Y","Z","cT", "[mm]");
65	ostr << outline << std::endl;
66	}
67
68	void HEPEVT_Wrapper::print_hepevt_particle( int i, std::ostream& ostr )
69	{
70	/// dumps the content HEPEVT particle entry i (Width is 120)
71	/// here i is the C array index (i.e. it starts at 0 ... whereas the
72	/// fortran array index starts at 1) So if there's 100 particles, the
73	/// last valid index is 100-1=99
74	char outline[81];
75	sprintf( outline,
76	"%4d %+4d %4d %4d (%9.3g, %9.3g, %9.3g, %9.3g, %9.3g)"
77	,i, status(i), first_parent(i), first_child(i),
78	px(i), py(i), pz(i), e(i), m(i) );
79	ostr << outline << "\n";
80	sprintf( outline,"%+9d %4d %4d (%9.3g, %9.3g, %9.3g, %9.3g)",
81	// old version was:" (%+9.2e, %+9.2e, %+9.2e, %+9.2e)"
82	id(i), last_parent(i), last_child(i),
83	x(i), y(i), z(i), t(i) );
84	ostr << outline << std::endl;
85	}
86
87
88	bool HEPEVT_Wrapper::check_hepevt_consistency( std::ostream& os )
89	{
90	/// This method inspects the HEPEVT common block and looks for
91	/// inconsistencies in the mother/daughter pointers
92	bool isConsistent=true;
93	char header[81];
94	sprintf( header,
95	"\n\n\t** WARNINGInconsistent HEPEVT input, Event %10d **"
96	, HEPEVT_Wrapper::event_number() );
97
98	for ( int i = 1; i <= HEPEVT_Wrapper::number_entries(); ++i ) {
99	// 1. check its mothers
100	int moth1 = HEPEVT_Wrapper::first_parent( i );
101	int moth2 = HEPEVT_Wrapper::last_parent( i );
102	if ( moth2<moth1 ) {
103	if ( isConsistent ) {
104	os << header << std::endl;
105	isConsistent = false;
106	print_legend(os);
107	}
108	os << "Inconsistent entry " << i
109	<< " first parent > last parent " << std::endl;
110	HEPEVT_Wrapper::print_hepevt_particle( i, os );
111	}
112	for ( int m = moth1; m<=moth2 && m!=0; ++m ) {
113	if ( m>HEPEVT_Wrapper::number_entries() \|\| m < 0 ) {
114	if ( isConsistent ) {
115	os << header << std::endl;
116	isConsistent = false;
117	print_legend(os);
118	}
119	os << "Inconsistent entry " << i
120	<< " mother points out of range " << std::endl;
121	HEPEVT_Wrapper::print_hepevt_particle( i, os );
122	}
123	int mChild1 = HEPEVT_Wrapper::first_child(m);
124	int mChild2 = HEPEVT_Wrapper::last_child(m);
125	// we don't consider null pointers as inconsistent
126	if ( mChild1==0 && mChild2==0 ) continue;
127	if ( i<mChild1 \|\| i>mChild2 ) {
128	if ( isConsistent ) {
129	os << header << std::endl;
130	isConsistent = false;
131	print_legend(os);
132	}
133	os << "Inconsistent mother-daughter relationship between "
134	<< i << " & " << m
135	<< " (try !trust_mother)" << std::endl;
136	HEPEVT_Wrapper::print_hepevt_particle( i, os );
137	HEPEVT_Wrapper::print_hepevt_particle( m, os );
138	}
139	}
140	// 2. check its daughters
141	int dau1 = HEPEVT_Wrapper::first_child( i );
142	int dau2 = HEPEVT_Wrapper::last_child( i );
143	if ( dau2<dau1 ) {
144	if ( isConsistent ) {
145	os << header << std::endl;
146	isConsistent = false;
147	print_legend(os);
148	}
149	os << "Inconsistent entry " << i
150	<< " first child > last child " << std::endl;
151	HEPEVT_Wrapper::print_hepevt_particle( i, os );
152	}
153	for ( int d = dau1; d<=dau2 && d!=0; ++d ) {
154	if ( d>HEPEVT_Wrapper::number_entries() \|\| d < 0 ) {
155	if ( isConsistent ) {
156	os << header << std::endl;
157	isConsistent = false;
158	print_legend(os);
159	}
160	os << "Inconsistent entry " << i
161	<< " child points out of range " << std::endl;
162	HEPEVT_Wrapper::print_hepevt_particle( i, os );
163	}
164	int d_moth1 = HEPEVT_Wrapper::first_parent(d);
165	int d_moth2 = HEPEVT_Wrapper::last_parent(d);
166	// we don't consider null pointers as inconsistent
167	if ( d_moth1==0 && d_moth2==0 ) continue;
168	if ( i<d_moth1 \|\| i>d_moth2 ) {
169	if ( isConsistent ) {
170	os << header << std::endl;
171	isConsistent = false;
172	print_legend(os);
173	}
174	os << "Inconsistent mother-daughter relationship between "
175	<< i << " & " << d
176	<< " (try trust_mothers)"<< std::endl;
177	HEPEVT_Wrapper::print_hepevt_particle( i, os );
178	HEPEVT_Wrapper::print_hepevt_particle( d, os );
179	}
180	}
181	}
182	if (!isConsistent) {
183	os << "Above lists all the inconsistencies in the HEPEVT common "
184	<< "\n block which has been provided as input to HepMC. "
185	<< "\n HepMC WILL have trouble interpreting the mother-daughter"
186	<< "\n relationships ... but all other information "
187	<< "\n (4-vectors etc) will be correctly transferred."
188	<< "\n In order for HepMC to be able to interpret the mother/"
189	<< "\n daughter hierachy, it MUST be given consistent input."
190	<< "\n This is one of the design criteria of HepMC: "
191	<< "\n consistency is enforced by the code.";
192	os << "\nThere is a switch in IO_HEPEVT, set-able using "
193	<< "\n IO_HEPEVT::set_trust_mothers_before_daughters( bool )"
194	<< "\n which you may want to try.";
195	os << "\nNote: if HEPEVT common block has been filled by pythia"
196	<< "\n pyhepc, then the switch MSTP(128)=2 should be used in"
197	<< "\n pythia, which instructs pythia not to put multiple "
198	<< "\n copies of resonances in the event record.\n";
199	os << "To obtain a file summarizing the inconsistency, you should:"
200	<< "\n\t ofstream myFile(\"myInconsistentEvent.txt\"); "
201	<< "\n\t HEPEVT_Wrapper::check_hepevt_consistency(myFile); "
202	<< "\n\t HEPEVT_Wrapper::print_hepevt(myFile); "
203	<< "\n[now write the event to HepMC using something like"
204	<< "\n\t\t myIO_HEPEVT->write_event(myEvent); ]"
205	<< "\n\t myEvent->print( myFile ); "
206	<< " // print event as HepMC sees it"
207	<< "\n ------------------------- Thank-you. \n\n" << std::endl;
208	}
209	return isConsistent;
210	}
211
212	void HEPEVT_Wrapper::zero_everything()
213	{
214	set_event_number( 0 );
215	set_number_entries( 0 );
216	for ( int i = 1; i<=max_number_entries(); ++i ) {
217	set_status( i, 0 );
218	set_id( i, 0 );
219	set_parents( i, 0, 0 );
220	set_children( i, 0, 0 );
221	set_momentum( i, 0, 0, 0, 0 );
222	set_mass( i, 0 );
223	set_position( i, 0, 0, 0, 0 );
224	}
225	}
226
227	} // HepMC
228