HepMC3 event record library
testMass.cc
1//-------------------------------------------------------------------
2// testMass.cc.in
3//
4// garren@fnal.gov, March 2006
5// Read events written by example_MyPythia.cc
6// Select events containing a photon of pT > 25 GeV
7// Add arbitrary PDF information to one of the good events
8// Add arbitrary HeavyIon information to one of the good events
9// Write the selected events and read them back in using an istream
10//-------------------------------------------------------------------
11
12#include <cmath> // for min()
13#include <ostream>
14
15#include "HepMC3/GenParticle.h"
16#include "HepMC3/GenEvent.h"
17#include "HepMC3/GenPdfInfo.h"
18#include "HepMC3/GenHeavyIon.h"
19
20
21#include "HepMC3/Version.h"
22#include "HepMC3/ReaderAscii.h"
23#include "HepMC3/WriterAscii.h"
26
27// define methods and classes used by this test
28#include "IsGoodEvent.h"
29using namespace HepMC3;
30bool massInfo( const GenEvent&, std::ostream& );
31
32int main()
33{
34 // declare an input strategy to read the data produced with the
35 ReaderAsciiHepMC2 ascii_in("inputMass.hepmc");
36 if (ascii_in.failed()) return 1;
37 // declare another IO_GenEvent for output
38 WriterAsciiHepMC2 ascii_out("testMass1.out");
39 // declare an instance of the event selection predicate
40 IsGoodEventDIS is_good_event;
41 // send version to output
43 //........................................EVENT LOOP
44 int icount=0;
45 int num_good_events=0;
46 double x1=0., x2=0., q=0., xf1=0., xf2=0.;
47 GenEvent evt;
48 while ( !ascii_in.failed() )
49 {
50 bool readOK=ascii_in.read_event(evt);
51 if (!readOK) return 1;
52 icount++;
53 if ( icount%50==1 ) std::cout << "Processing Event Number " << icount<< " its # " << evt.event_number() << std::endl;
54 if ( is_good_event(evt) )
55 {
56 if (num_good_events == 0 )
57 {
58 // add some arbitrary PDF information
59 x1 = std::min(0.8, 0.07 * icount);
60 x2 = 1-x1;
61 q = 1.69 * icount;
62 // use beam momentum
63 if( evt.beams().size()==2 )
64 {
65 GenParticlePtr bp1 = evt.beams().at(0);
66 GenParticlePtr bp2 = evt.beams().at(1);
67 xf1 = x1*bp1->momentum().p3mod();
68 xf2 = x2*bp1->momentum().p3mod();
69 }
70 else
71 {
72 xf1 = x1*0.34;
73 xf2 = x2*0.34;
74 }
75 // provide optional pdf set id numbers (two ints at the end of the constructor)
76 std::shared_ptr< GenPdfInfo> pdf = std::make_shared< GenPdfInfo>();
77 evt.add_attribute("GenPdfInfo",pdf);
78 pdf->set( 2, 3, x1, x2, q, xf1, xf2, 230, 230);
79 // add some arbitrary HeavyIon information
80 std::shared_ptr< GenHeavyIon> ion = std::make_shared< GenHeavyIon>();
81 evt.add_attribute("GenHeavyIon",ion);
82 ion->set(23,11,12,15,3,5,0,0,0,0.0145,0.0,0.0,0.0,0.23);
83 }
84 std::cout << "saving Event " << evt.event_number() << std::endl;
85 if( evt.weights().size() > 0 )
86 {
87 std::cout << "Weights: ";
88 for ( std::vector<double>::const_iterator w=evt.weights().begin(); w!=evt.weights().end(); ++w )
89 std::cout <<" "<<*w;
90 std::cout << std::endl;
91 }
92 ascii_out.write_event(evt);
93 ++num_good_events;
94 }
95 // clean up and get next event
96 evt.clear();
97 }
98 //........................................PRINT RESULT
99 std::cout << num_good_events << " out of " << icount
100 << " processed events passed the cuts. Finished." << std::endl;
101 ascii_in.close();
102 ascii_out.close();
103 // now read the file we just created
104 // declare an input strategy
105 std::ifstream istr( "testMass1.out" );
106 if( !istr )
107 {
108 std::cerr << "testMass: cannot open " << std::endl;
109 exit(1);
110 }
111 ReaderAsciiHepMC2 xin(istr);
112 if (xin.failed()) return 1;
113 // declare another IO_GenEvent for output
114 WriterAsciiHepMC2 xout("testMass2.out");
115 if (xout.failed()) return 1;
116 //........................................EVENT LOOP
117 int ixin=0;
118 while ( !xin.failed() )
119 {
120 bool readOK=xin.read_event(evt);
121 if (!readOK) return 1;
122 ixin++;
123 std::cout << "reading Event " << evt.event_number() << std::endl;
124 if( evt.weights().size() > 0 )
125 {
126 std::cout << "Weights: ";
127 for ( std::vector<double>::const_iterator w=evt.weights().begin(); w!=evt.weights().end(); ++w )
128 std::cout <<" "<<*w;
129 std::cout << std::endl;
130 }
131 xout.write_event(evt);
132 // look at mass info
133 if (! massInfo(evt,std::cout)) return 1;
134 // clean up and get next event
135 evt.clear();
136 }
137 //........................................PRINT RESULT
138 std::cout << ixin << " events in the second pass. Finished." << std::endl;
139 xin.close();
140 xout.close();
141 return 0;
142}
143
144bool massInfo( const GenEvent& e, std::ostream& os )
145{
146 for (ConstGenParticlePtr p: e.particles()) {
147 double gm = p->generated_mass();
148 double m = p->momentum().m();
149 double d = std::abs(m-gm);
150 if( d > 1.0e-4 && gm>1.0e-4)
151 {
152 os << "Event " << e.event_number()
153 << " Particle " << (p)->pdg_id()
154 << " generated mass " << gm
155 << " mass from momentum " << m
156 << " difference " << d << std::endl;
157 return false;
158 }
159 }
160 return true;
161}
Definition of class GenEvent.
Definition of attribute class GenHeavyIon.
Definition of class GenParticle.
Definition of event attribute class GenPdfInfo.
Definition of class ReaderAsciiHepMC2.
Definition of class ReaderAscii.
Definition of class WriterAsciiHepMC2.
Definition of class WriterAscii.
Stores event-related information.
Definition: GenEvent.h:41
int event_number() const
Get event number.
Definition: GenEvent.h:148
std::vector< ConstGenParticlePtr > beams() const
Vector of beam particles.
Definition: GenEvent.cc:416
void add_attribute(const std::string &name, const std::shared_ptr< Attribute > &att, const int &id=0)
Definition: GenEvent.cc:806
const std::vector< double > & weights() const
Get event weight values as a vector.
Definition: GenEvent.h:98
void clear()
Remove contents of this event.
Definition: GenEvent.cc:599
const std::vector< ConstGenParticlePtr > & particles() const
Get list of particles (const)
Definition: GenEvent.cc:39
Parser for HepMC2 I/O files.
GenEvent I/O serialization for structured text files.
HepMC3 main namespace.
std::string version()
Get the HepMC library version string.
Definition: Version.h:20