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 }
GenEvent I/O serialization for structured text files.
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
int event_number() const
Get event number.
Definition: GenEvent.h:135
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
Parser for HepMC2 I/O files.
std::vector< ConstGenParticlePtr > beams() const
Vector of beam particles.
Definition: GenEvent.cc:420
Definition of class ReaderAsciiHepMC2.
Stores event-related information.
Definition: GenEvent.h:41
Definition of class ReaderAscii.
Definition of class WriterAsciiHepMC2.
Definition of event attribute class GenPdfInfo.
int main(int argc, char **argv)
const std::vector< double > & weights() const
Get event weight values as a vector.
Definition: GenEvent.h:86
Definition of class GenEvent.
void clear()
Remove contents of this event.
Definition: GenEvent.cc:608
const std::vector< ConstGenParticlePtr > & particles() const
Get list of particles (const)
Definition: GenEvent.cc:39
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