benchmarks/zdc/zdc_neutrons_reader.cxx

0001 #include "HepMC3/GenEvent.h"
0002 #include "HepMC3/ReaderAscii.h"
0003 #include "HepMC3/WriterAscii.h"
0004 #include "HepMC3/Print.h"
0005
0006 #include "TH1F.h"
0007 #include <iostream>
0008 #include <cstdlib>
0009 using namespace HepMC3;
0010
0011 void zdc_neutrons_reader(){
0012
0013   //-------------------------------------
0014   ReaderAscii hepmc_input("data/neutrons_zdc.hepmc");
0015   if( hepmc_input.failed() ) {
0016     std::cerr << " could not find data file\n";
0017     std::quick_exit(1);
0018   }
0019
0020
0021   int        events_parsed = 0;
0022   GenEvent   evt(Units::GEV, Units::MM);
0023
0024   TH1F* h_neutron_energy = new TH1F("n energy","; E [GeV]",100,0,200);
0025
0026   // Read event from input file
0027   hepmc_input.read_event(evt);
0028   // If reading failed before loop then fail hard
0029   // because the data wasn't found.
0030   if( hepmc_input.failed() ) {
0031     std::cerr << " could not find first event\n";
0032     std::quick_exit(1);
0033   }
0034
0035   while(!hepmc_input.failed()) {
0036
0037
0038     for(const auto& v : evt.vertices() ) {
0039       for(const auto& p : v->particles_out() ) {
0040         if(p->pid() == 2112) {
0041           h_neutron_energy->Fill(p->momentum().e());
0042         }
0043       }
0044     }
0045     evt.clear();
0046     events_parsed++;
0047     hepmc_input.read_event(evt);
0048   }
0049   std::cout << "Events parsed and written: " << events_parsed << std::endl;
0050
0051   TCanvas* c = new TCanvas();
0052   h_neutron_energy->Draw();
0053   c->SaveAs("results/zdc_neutrons_reader.png");
0054   c->SaveAs("results/zdc_neutrons_reader.pdf");
0055 }
0056
0057