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