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"
29 using namespace HepMC3;
30 bool massInfo( const GenEvent&, std::ostream& );
31 
32 int 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 
144 bool 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:135
std::vector< ConstGenParticlePtr > beams() const
Vector of beam particles.
Definition: GenEvent.cc:420
void add_attribute(const std::string &name, const std::shared_ptr< Attribute > &att, const int &id=0)
Add event attribute to event.
Definition: GenEvent.h:208
void clear()
Remove contents of this event.
Definition: GenEvent.cc:608
const std::vector< double > & weights() const
Get event weight values as a vector.
Definition: GenEvent.h:86
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.
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
std::string version()
Get the HepMC library version string.
Definition: Version.h:20
int main(int argc, char **argv)