hepmc - Blame information for rev 96

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