generators/flowAfterburner/test.cc

0001 //
0002 // Inspired by code from ATLAS.  Thanks!
0003 //
0004 #include <HepMC/GenEvent.h>
0005 #include <HepMC/GenParticle.h>
0006 #include <HepMC/GenRanges.h>
0007 #include <HepMC/GenVertex.h>
0008 #include <HepMC/HeavyIon.h>
0009 #include <HepMC/IO_BaseClass.h>
0010 #include <HepMC/IO_GenEvent.h>
0011 #include <HepMC/IteratorRange.h>
0012 #include <HepMC/SimpleVector.h>
0013
0014 // this is an ugly hack, the gcc optimizer has a bug which
0015 // triggers the uninitialized variable warning which
0016 // stops compilation because of our -Werror
0017 #include <boost/version.hpp>  // to get BOOST_VERSION
0018 #if (__GNUC__ == 4 && __GNUC_MINOR__ == 8 && (BOOST_VERSION == 106000 || BOOST_VERSION == 106700 || BOOST_VERSION == 107000))
0019 #pragma GCC diagnostic ignored "-Wshadow"
0020 #pragma GCC diagnostic ignored "-Wunused-parameter"
0021 #pragma message "ignoring bogus gcc warning in boost header ptree.hpp"
0022 #include <boost/property_tree/ptree.hpp>
0023 #pragma GCC diagnostic warning "-Wshadow"
0024 #pragma GCC diagnostic warning "-Wunused-parameter"
0025 #else
0026 #include <boost/property_tree/ptree.hpp>
0027 #endif
0028
0029 #include <boost/operators.hpp>
0030 #include <boost/property_tree/xml_parser.hpp>
0031
0032 #include <gsl/gsl_histogram.h>
0033
0034 #include <algorithm>  // for max
0035 #include <cmath>      // for cos
0036 #include <cstdlib>
0037 #include <iostream>
0038 #include <string>
0039
0040 // NOLINTNEXTLINE(bugprone-exception-escape)
0041 int main()
0042 {
0043   using boost::property_tree::ptree;
0044   ptree proptree;
0045
0046   std::ifstream config_file("test.xml");
0047
0048   if (config_file)
0049   {
0050     // Read XML configuration file.
0051     read_xml(config_file, proptree);
0052   }
0053
0054   std::string input = proptree.get("TEST.INPUT", "test.dat");
0055
0056   // Try to open input file.
0057   std::ifstream istr(input.c_str());
0058   if (!istr)
0059   {
0060     std::cout << __PRETTY_FUNCTION__ << ": input file \""
0061               << input << "\" not found!" << std::endl;
0062     return 1;
0063   }
0064
0065   // Book a GSL histogram
0066   size_t n = 32;
0067   double a;
0068   double b;
0069   double w = M_PI / n;
0070   a = -M_PI / 2.0 - w / 2.0;
0071   b = M_PI / 2.0 - w / 2.0;
0072   gsl_histogram *h = gsl_histogram_alloc(n);
0073   gsl_histogram_set_ranges_uniform(h, a, b);
0074
0075   size_t N = 20;
0076   double A = 0.2;
0077   double B = 5.0;
0078   gsl_histogram *v2h = gsl_histogram_alloc(N);
0079   gsl_histogram_set_ranges_uniform(v2h, A, B);
0080   gsl_histogram *v2hn = gsl_histogram_alloc(N);
0081   gsl_histogram_set_ranges_uniform(v2hn, A, B);
0082
0083   HepMC::IO_GenEvent ascii_in(istr);
0084   HepMC::GenEvent *evt;
0085
0086   while (ascii_in >> evt)
0087   {
0088     HepMC::HeavyIon *hi = evt->heavy_ion();
0089     HepMC::GenVertex *primary_vtx = evt->barcode_to_vertex(-1);
0090     double phi0 = hi->event_plane_angle();
0091
0092     HepMC::GenVertexParticleRange r(*primary_vtx, HepMC::children);
0093     for (HepMC::GenVertex::particle_iterator it = r.begin(); it != r.end(); it++)
0094     {
0095       HepMC::GenParticle *p1 = (*it);
0096       if (p1->status() == 103)
0097       {
0098         continue;
0099       }
0100       double eta1 = p1->momentum().eta();
0101       if (fabs(eta1) > 1.5)
0102       {
0103         continue;
0104       }
0105       double pt = p1->momentum().perp();
0106       if (pt < 0.2 || pt > 10.0)
0107       {
0108         continue;
0109       }
0110       double phi1 = p1->momentum().phi();
0111       double dphi = phi1 - phi0;
0112       dphi = atan2(sin(dphi), cos(dphi));
0113       gsl_histogram_accumulate(v2h, pt, cos(2 * dphi));
0114       gsl_histogram_increment(v2hn, pt);
0115       gsl_histogram_increment(h, dphi);
0116     }
0117
0118     delete evt;
0119     ascii_in >> evt;
0120   }
0121
0122   std::cout << "pt,v2,v2e" << std::endl;
0123   gsl_histogram_div(v2h, v2hn);
0124   double lower;
0125   double upper;
0126   for (size_t i = 0; i < N; i++)
0127   {
0128     gsl_histogram_get_range(v2h, i, &lower, &upper);
0129     double pt = (lower + upper) / 2.0;
0130     double counts = gsl_histogram_get(v2hn, i);
0131     double val = gsl_histogram_get(v2h, i);
0132     double err = val / sqrt(counts);
0133     std::cout << pt << ", " << val << ", " << err << std::endl;
0134   }
0135
0136   return 0;
0137 }