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>
17 #include <boost/version.hpp>
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"
26 #include <boost/property_tree/ptree.hpp>
29 #include <boost/operators.hpp>
30 #include <boost/property_tree/xml_parser.hpp>
32 #include <gsl/gsl_histogram.h>
42 using boost::property_tree::ptree;
50 read_xml(config_file, proptree);
53 std::string input = proptree.get(
"TEST.INPUT",
"test.dat");
56 std::ifstream istr(input.c_str());
59 std::cout << __PRETTY_FUNCTION__ <<
": input file \""
60 << input <<
"\" not found!" << std::endl;
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);
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);
80 HepMC::IO_GenEvent ascii_in(istr);
83 while (ascii_in >> evt)
85 HepMC::HeavyIon *hi = evt->heavy_ion();
86 HepMC::GenVertex *primary_vtx = evt->barcode_to_vertex(-1);
87 double phi0 = hi->event_plane_angle();
89 HepMC::GenVertexParticleRange
r(*primary_vtx, HepMC::children);
90 for (HepMC::GenVertex::particle_iterator
it = r.begin();
it != r.end();
it++)
92 HepMC::GenParticle *p1 = (*it);
93 if (p1->status() == 103)
97 double eta1 = p1->momentum().eta();
102 double pt = p1->momentum().perp();
103 if (pt < 0.2 || pt > 10.0)
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);
119 std::cout <<
"pt,v2,v2e" << std::endl;
120 gsl_histogram_div(v2h, v2hn);
122 for (
size_t i = 0; i <
N; i++)
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;