ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
test.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file test.cc
1 //
2 // Inspired by code from ATLAS. Thanks!
3 //
4 #include <HepMC/GenEvent.h>
5 #include <HepMC/GenParticle.h>
6 #include <HepMC/GenRanges.h>
7 #include <HepMC/GenVertex.h>
8 #include <HepMC/HeavyIon.h>
9 #include <HepMC/IO_BaseClass.h>
10 #include <HepMC/IO_GenEvent.h>
11 #include <HepMC/IteratorRange.h>
12 #include <HepMC/SimpleVector.h>
13 
14 // this is an ugly hack, the gcc optimizer has a bug which
15 // triggers the uninitialized variable warning which
16 // stops compilation because of our -Werror
17 #include <boost/version.hpp> // to get BOOST_VERSION
18 #if (__GNUC__ == 4 && __GNUC_MINOR__ == 8 && (BOOST_VERSION == 106000 || BOOST_VERSION == 106700 || BOOST_VERSION == 107000))
19 #pragma GCC diagnostic ignored "-Wshadow"
20 #pragma GCC diagnostic ignored "-Wunused-parameter"
21 #pragma message "ignoring bogus gcc warning in boost header ptree.hpp"
22 #include <boost/property_tree/ptree.hpp>
23 #pragma GCC diagnostic warning "-Wshadow"
24 #pragma GCC diagnostic warning "-Wunused-parameter"
25 #else
26 #include <boost/property_tree/ptree.hpp>
27 #endif
28 
29 #include <boost/operators.hpp>
30 #include <boost/property_tree/xml_parser.hpp>
31 
32 #include <gsl/gsl_histogram.h>
33 
34 #include <algorithm> // for max
35 #include <cmath> // for cos
36 #include <cstdlib>
37 #include <iostream>
38 #include <string>
39 
40 int main()
41 {
42  using boost::property_tree::ptree;
43  ptree proptree;
44 
45  std::ifstream config_file("test.xml");
46 
47  if (config_file)
48  {
49  // Read XML configuration file.
50  read_xml(config_file, proptree);
51  }
52 
53  std::string input = proptree.get("TEST.INPUT", "test.dat");
54 
55  // Try to open input file.
56  std::ifstream istr(input.c_str());
57  if (!istr)
58  {
59  std::cout << __PRETTY_FUNCTION__ << ": input file \""
60  << input << "\" not found!" << std::endl;
61  return 1;
62  }
63 
64  // Book a GSL histogram
65  size_t n = 32;
66  double a, b;
67  double w = M_PI / n;
68  a = -M_PI / 2.0 - w / 2.0;
69  b = M_PI / 2.0 - w / 2.0;
70  gsl_histogram *h = gsl_histogram_alloc(n);
71  gsl_histogram_set_ranges_uniform(h, a, b);
72 
73  size_t N = 20;
74  double A = 0.2, B = 5.0;
75  gsl_histogram *v2h = gsl_histogram_alloc(N);
76  gsl_histogram_set_ranges_uniform(v2h, A, B);
77  gsl_histogram *v2hn = gsl_histogram_alloc(N);
78  gsl_histogram_set_ranges_uniform(v2hn, A, B);
79 
80  HepMC::IO_GenEvent ascii_in(istr);
81  HepMC::GenEvent *evt;
82 
83  while (ascii_in >> evt)
84  {
85  HepMC::HeavyIon *hi = evt->heavy_ion();
86  HepMC::GenVertex *primary_vtx = evt->barcode_to_vertex(-1);
87  double phi0 = hi->event_plane_angle();
88 
89  HepMC::GenVertexParticleRange r(*primary_vtx, HepMC::children);
90  for (HepMC::GenVertex::particle_iterator it = r.begin(); it != r.end(); it++)
91  {
92  HepMC::GenParticle *p1 = (*it);
93  if (p1->status() == 103)
94  {
95  continue;
96  }
97  double eta1 = p1->momentum().eta();
98  if (fabs(eta1) > 1.5)
99  {
100  continue;
101  }
102  double pt = p1->momentum().perp();
103  if (pt < 0.2 || pt > 10.0)
104  {
105  continue;
106  }
107  double phi1 = p1->momentum().phi();
108  double dphi = phi1 - phi0;
109  dphi = atan2(sin(dphi), cos(dphi));
110  gsl_histogram_accumulate(v2h, pt, cos(2 * dphi));
111  gsl_histogram_increment(v2hn, pt);
112  gsl_histogram_increment(h, dphi);
113  }
114 
115  delete evt;
116  ascii_in >> evt;
117  }
118 
119  std::cout << "pt,v2,v2e" << std::endl;
120  gsl_histogram_div(v2h, v2hn);
121  double lower, upper;
122  for (size_t i = 0; i < N; i++)
123  {
124  gsl_histogram_get_range(v2h, i, &lower, &upper);
125  double pt = (lower + upper) / 2.0;
126  double counts = gsl_histogram_get(v2hn, i);
127  double val = gsl_histogram_get(v2h, i);
128  double err = val / sqrt(counts);
129  std::cout << pt << ", " << val << ", " << err << std::endl;
130  }
131 
132  return 0;
133 }