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"
27 using namespace HepMC3;
28 bool massInfo( const GenEvent&, std::ostream& );
29 
30 int 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 
142 bool 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 }
int event_number() const
Get event number.
Definition: GenEvent.h:136
GenEvent I/O serialization for structured text files.
HepMC3 main namespace.
Definition: WriterDOT.h:19
Definition of class GenParticle.
Definition of attribute class GenHeavyIon.
Definition of class WriterAscii.
std::string version()
Get the HepMC library version string.
Definition: Version.h:20
Parser for HepMC2 I/O files.
const std::vector< double > & weights() const
Get event weight values as a vector.
Definition: GenEvent.h:87
Definition of class ReaderAsciiHepMC2.
Stores event-related information.
Definition: GenEvent.h:42
Definition of class ReaderAscii.
const std::vector< ConstGenParticlePtr > & particles() const
Get list of particles (const)
Definition: GenEvent.cc:40
Definition of class WriterAsciiHepMC2.
void add_attribute(const string &name, const shared_ptr< Attribute > &att, const int &id=0)
Add event attribute to event.
Definition: GenEvent.h:202
Definition of event attribute class GenPdfInfo.
int main(int argc, char **argv)
std::vector< ConstGenParticlePtr > beams() const
Vector of beam particles.
Definition: GenEvent.cc:419
Definition of class GenEvent.
void clear()
Remove contents of this event.
Definition: GenEvent.cc:438
Feature< Feature_type > abs(const Feature< Feature_type > &input)
Obtain the absolute value of a Feature. This works as you&#39;d expect. If foo is a valid Feature...
Definition: Feature.h:316