hepmc - Blame information for rev 128

Subversion Repositories:
Rev:
Rev Author Line No. Line
88 garren 1 //////////////////////////////////////////////////////////////////////////
2 // testHepMCIteration.cc.in
3 //
4 // garren@fnal.gov, May 2007
5 // Use Matt's example_EventSelection along with example_UsingIterators
6 // to check HepMC iteration.
7 // Apply an event selection to the events in testHepMC.input
8 // Events containing a photon of pT > 25 GeV pass the selection.
9 // Use iterators on these events.
10 //////////////////////////////////////////////////////////////////////////
11  
12 #include <list>
13  
128 garren 14 #include "HepMC/IO_GenEvent.h"
88 garren 15 #include "HepMC/IO_AsciiParticles.h"
16 #include "HepMC/GenEvent.h"
17  
18 // define methods and classes used by this test
19 #include "IsGoodEvent.h"
20 #include "testHepMCIteration.h"
21  
22 bool findW( HepMC::GenEvent* evt, std::ofstream& os);
23 bool simpleIter( HepMC::GenEvent* evt, std::ofstream& os );
24  
25 int main() {
26 // declare an input strategy to read the data produced with the
27 // example_MyPythia
128 garren 28 HepMC::IO_GenEvent ascii_in("@srcdir@/testIOGenEvent.input",std::ios::in);
88 garren 29 // declare an instance of the event selection predicate
30 IsGoodEvent is_good_event;
31 // define an output stream
32 std::ofstream os( "testHepMCIteration.dat" );
33 //........................................EVENT LOOP
34 int icount=0;
35 int num_good_events=0;
36 HepMC::GenEvent* evt = ascii_in.read_next_event();
96 garren 37 HepMC::GenEvent* evcopy;
88 garren 38 while ( evt ) {
39 icount++;
40 if ( icount%50==1 ) std::cout << "Processing Event Number " << icount
41 << " its # " << evt->event_number()
42 << std::endl;
43 // icount of 100 should be the last event
44 if ( icount==100 ) std::cout << "Processing Event Number " << icount
45 << " its # " << evt->event_number()
46 << std::endl;
96 garren 47 evcopy = evt;
48 if ( is_good_event(evcopy) ) {
49 os << "Event " << evcopy->event_number() << " is good " << std::endl;
88 garren 50 ++num_good_events;
51 // test iterators
96 garren 52 simpleIter( evcopy, os );
53 findW( evcopy, os );
88 garren 54 }
96 garren 55 evcopy->clear();
88 garren 56  
57 // clean up and get next event
58 delete evt;
59 evt = ascii_in.read_next_event();
60 }
61 //........................................PRINT RESULT
62 std::cout << num_good_events << " out of " << icount
63 << " processed events passed the cuts. Finished." << std::endl;
64 }
65  
66 bool simpleIter( HepMC::GenEvent* evt, std::ofstream& os )
67 {
68 // use GenEvent::vertex_iterator to fill a list of all
69 // vertices in the event
70 std::list<HepMC::GenVertex*> allvertices;
71 for ( HepMC::GenEvent::vertex_iterator v = evt->vertices_begin();
72 v != evt->vertices_end(); ++v ) {
73 allvertices.push_back(*v);
74 }
75  
76 // we could do the same thing with the STL algorithm copy
77 std::list<HepMC::GenVertex*> allvertices2;
78 copy( evt->vertices_begin(), evt->vertices_end(),
79 back_inserter(allvertices2) );
80  
81 // fill a list of all final state particles in the event, by requiring
82 // that each particle satisfyies the IsFinalState predicate
83 IsFinalState isfinal;
84 std::list<HepMC::GenParticle*> finalstateparticles;
85 for ( HepMC::GenEvent::particle_iterator p = evt->particles_begin();
86 p != evt->particles_end(); ++p ) {
87 if ( isfinal(*p) ) finalstateparticles.push_back(*p);
88 }
89  
90 // an STL-like algorithm called HepMC::copy_if is provided in the
91 // GenEvent.h header to do this sort of operation more easily,
92 // you could get the identical results as above by using:
93 std::list<HepMC::GenParticle*> finalstateparticles2;
94 HepMC::copy_if( evt->particles_begin(), evt->particles_end(),
95 back_inserter(finalstateparticles2), IsFinalState() );
96  
97 // print all photons in the event that satisfy the IsPhoton criteria
98 os << "photons in event " << evt->event_number() << ":" << std::endl;
99 for ( HepMC::GenEvent::particle_iterator p = evt->particles_begin();
100 p != evt->particles_end(); ++p ) {
101 if ( IsPhoton(*p) ) (*p)->print( os );
102 }
103 return true;
104 }
105  
106 bool findW( HepMC::GenEvent* evt, std::ofstream& os )
107 {
108 int num_W=0;
109 // use GenEvent::particle_iterator to find all W's in the event,
110 // then
111 // (1) for each W user the GenVertex::particle_iterator with a range of
112 // parents to return and print the immediate mothers of these W's.
113 // (2) for each W user the GenVertex::particle_iterator with a range of
114 // descendants to return and print all descendants of these W's.
115 for ( HepMC::GenEvent::particle_iterator p = evt->particles_begin();
116 p != evt->particles_end(); ++p ) {
117 if ( IsWBoson(*p) ) {
118 ++num_W;
119 os << "A W boson has been found in event: " << evt->event_number() << std::endl;
120 (*p)->print( os );
121 // return all parents
122 // we do this by pointing to the production vertex of the W
123 // particle and asking for all particle parents of that vertex
124 os << "\t Its parents are: " << std::endl;
125 if ( (*p)->production_vertex() ) {
126 for ( HepMC::GenVertex::particle_iterator mother
127 = (*p)->production_vertex()->
128 particles_begin(HepMC::parents);
129 mother != (*p)->production_vertex()->
130 particles_end(HepMC::parents);
131 ++mother ) {
132 os << "\t";
133 (*mother)->print( os );
134 }
135 }
136 // return all descendants
137 // we do this by pointing to the end vertex of the W
138 // particle and asking for all particle descendants of that vertex
139 os << "\t\t Its descendants are: " << std::endl;
140 if ( (*p)->end_vertex() ) {
141 for ( HepMC::GenVertex::particle_iterator des
142 =(*p)->end_vertex()->
143 particles_begin(HepMC::descendants);
144 des != (*p)->end_vertex()->
145 particles_end(HepMC::descendants);
146 ++des ) {
147 os << "\t\t";
148 (*des)->print( os );
149 }
150 }
151 }
152 }
153 return true;
154 }
155