1 #include "AnalysisExample.h"
2 #include <limits>
3 namespace HepMC3
4 {
5 HEPMC3_DECLARE_WRITER_FILE(AnalysisExample)
HEPMC3_DECLARE_WRITER_STREAM(AnalysisExample)6 HEPMC3_DECLARE_WRITER_STREAM(AnalysisExample)
7
8 AnalysisExample::AnalysisExample(const std::string &filename,std::shared_ptr<GenRunInfo> /*run*/): m_file(filename),
9 m_stream(&m_file)
10 {
11 if ( !m_file.is_open() ) {
12 HEPMC3_ERROR( "AnalysisExample: could not open output file: "<<filename )
13 }
14 m_sum_of_weights=0;
15 m_sum_of_weights2=0;
16 m_bins["rapidity"]=std::vector<double> {-std::numeric_limits<double>::infinity(), -5.0,-4.5,-4.0,-3.5,-3.0,-2.5,-2.0,-1.5,-1.0,-0.5,0.0,0.5,1.0,1.5,2.0,2.5,3.0,3.5,4.0,4.5,5.0,std::numeric_limits<double>::infinity()};
17 m_vals["rapidity"]=std::vector<double>(m_bins.at("rapidity").size()-1,0.0);
18 m_errs["rapidity"]=std::vector<double>(m_bins.at("rapidity").size()-1,0.0);
19 }
20
AnalysisExample(std::ostream & stream,std::shared_ptr<GenRunInfo> run)21 AnalysisExample::AnalysisExample(std::ostream &stream, std::shared_ptr<GenRunInfo> run)
22 : m_file(),
23 m_stream(&stream)
24 {
25 m_sum_of_weights=0;
26 m_sum_of_weights2=0;
27 m_bins["rapidity"]=std::vector<double> {-std::numeric_limits<double>::infinity(), -5.0,-4.5,-4.0,-3.5,-3.0,-2.5,-2.0,-1.5,-1.0,-0.5,0.0,0.5,1.0,1.5,2.0,2.5,3.0,3.5,4.0,4.5,5.0,std::numeric_limits<double>::infinity()};
28 m_vals["rapidity"]=std::vector<double>(m_bins.at("rapidity").size()-1,0.0);
29 m_errs["rapidity"]=std::vector<double>(m_bins.at("rapidity").size()-1,0.0);
30
31 }
32
write_event(const GenEvent & evt)33 void AnalysisExample::write_event(const GenEvent &evt)
34 {
35 double w=evt.weight();
36 m_sum_of_weights+=w;
37 m_sum_of_weights2+=w*w;
38 for(auto p: evt.particles() )
39 {
40 if (p->status()!=1) continue;
41 double eta=p->momentum().eta();
42 int bin=std::distance(m_bins["rapidity"].begin(), lower_bound(m_bins["rapidity"].begin(),m_bins["rapidity"].end(),eta))-1;
43 if (bin<0) bin=0;
44 m_vals["rapidity"][bin]+=w;
45 m_errs["rapidity"][bin]+=w*w;
46 }
47 }
48
close()49 void AnalysisExample::close() {
50 if (!m_stream) return;
51 std::ofstream* ofs = dynamic_cast<std::ofstream*>(m_stream);
52 for (size_t i=1; i<m_vals["rapidity"].size()-1; i++)
53 {
54 double val=m_vals["rapidity"][i]/m_sum_of_weights/(m_bins["rapidity"][i+1]-m_bins["rapidity"][i]);
55 double err=sqrt(m_errs["rapidity"][i])/m_sum_of_weights/(m_bins["rapidity"][i+1]-m_bins["rapidity"][i]);
56 (*ofs)<< std::fixed << std::setprecision( 6 )<<m_bins["rapidity"][i]<<" "<<m_bins["rapidity"][i+1]<<" "<<val<<" "<<err<<std::endl;
57 }
58 if (ofs && !ofs->is_open()) return;
59 if (ofs) ofs->close();
60
61 }
62
63 } // namespace HepMC3
64