EIC-SynRad/SynRadAna/ReadSynRadFiles.cc

0001 #include "ReadSynRadFiles.h"
0002
0003 #include <fun4all/Fun4AllReturnCodes.h>
0004
0005 #include <HepMC/GenEvent.h>
0006 #include <HepMC/GenParticle.h>  // for GenParticle
0007 #include <HepMC/GenVertex.h>
0008 #include <HepMC/PdfInfo.h>       // for PdfInfo
0009 #include <HepMC/SimpleVector.h>  // for FourVector
0010 #include <HepMC/Units.h>         // for GEV, MM
0011
0012 // eicsmear classes
0013 #include <eicsmear/erhic/EventMC.h>
0014 #include <eicsmear/erhic/ParticleMC.h>  // for ParticleMC
0015 #include <eicsmear/erhic/Pid.h>         // for Pid
0016
0017 // General Root and C++ classes
0018 #include <TBranch.h>  // for TBranch
0019 #include <TChain.h>
0020 #include <TVector3.h>  // for TVector3
0021
0022 #include <algorithm>  // copy
0023 #include <cassert>
0024 #include <iostream>  // for operator<<, endl, basic_ostream
0025 #include <iterator>  // ostream_operator
0026 #include <string>
0027 #include <vector>  // for vector
0028 #include <vector>
0029
0030 #include <boost/tokenizer.hpp>
0031
0032 class PHCompositeNode;
0033 class PHHepMCGenEvent;
0034
0035 using namespace std;
0036
0037 ///////////////////////////////////////////////////////////////////
0038
0039 ReadSynRadFiles::ReadSynRadFiles(const string &name)
0040   : SubsysReco(name)
0041   , filename("")
0042   , nEntries(1)
0043   , entry(1)
0044   , GenEvent(nullptr)
0045   , _node_name("PHHepMCGenEvent")
0046 {
0047   hepmc_helper.set_embedding_id(1);  // default embedding ID to 1
0048   return;
0049 }
0050
0051 ///////////////////////////////////////////////////////////////////
0052
0053 ReadSynRadFiles::~ReadSynRadFiles()
0054 {
0055 }
0056
0057 ///////////////////////////////////////////////////////////////////
0058
0059 bool ReadSynRadFiles::OpenInputFile(const string &name)
0060 {
0061   filename = name;
0062   m_csv_input.open(filename);
0063   assert(m_csv_input.is_open());
0064
0065   // skip header
0066   string line;
0067   std::getline(m_csv_input, line);
0068
0069   return true;
0070 }
0071
0072 int ReadSynRadFiles::Init(PHCompositeNode *topNode)
0073 {
0074   /* Create node tree */
0075   CreateNodeTree(topNode);
0076   return 0;
0077 }
0078
0079 int ReadSynRadFiles::process_event(PHCompositeNode *topNode)
0080 {
0081   /* Create GenEvent */
0082   HepMC::GenEvent *evt = new HepMC::GenEvent();
0083
0084   /* define the units (Pythia uses GeV and mm) */
0085   evt->use_units(HepMC::Units::GEV, HepMC::Units::CM);
0086
0087   /* Loop over all particles for this event in input file and fill
0088    * vector with HepMC particles */
0089   vector<HepMC::GenParticle *> hepmc_particles;
0090   vector<unsigned> origin_index;
0091   //  /* add HepMC particles to Hep MC vertices; skip first two particles
0092   //   * in loop, assuming that they are the beam particles */
0093   vector<HepMC::GenVertex *> hepmc_vertices;
0094
0095   typedef boost::tokenizer<boost::escaped_list_separator<char> > Tokenizer;
0096
0097   vector<string> vec;
0098   string line;
0099
0100   int SumPhoton(0);
0101   double SumFlux(0);
0102
0103   if (Verbosity())
0104   {
0105     cout << "ReadSynRadFiles::process_event - reading " << nEntries << " lines from " << filename << endl;
0106   }
0107   for (int ii = 1; ii <= nEntries; ii++)
0108   {
0109     if (not std::getline(m_csv_input, line))
0110     {
0111       cout << "ReadSynRadFiles::process_event - "
0112            << "input file end reached" << endl;
0113
0114       return Fun4AllReturnCodes::ABORTRUN;
0115     }
0116
0117     Tokenizer tok(line);
0118     vec.assign(tok.begin(), tok.end());
0119
0120     if (vec.size() != 14 and vec.size() != 15)
0121     {
0122       cout << "ReadSynRadFiles::process_event - "
0123            << "invalid input :" << line << endl;
0124
0125       return Fun4AllReturnCodes::ABORTRUN;
0126     }
0127
0128     //    Pos_X_[cm]  Pos_Y_[cm]  Pos_Z_[cm]  Pos_u Pos_v Dir_X Dir_Y Dir_Z Dir_theta_[rad] Dir_phi_[rad] LowFluxRatio  Energy_[eV] Flux_[photon/s] Power_[W]
0129
0130     int id = 0;
0131     const double Pos_X_cm = stod(vec[id++]);
0132     const double Pos_Y_cm = stod(vec[id++]);
0133     const double Pos_Z_cm = stod(vec[id++]);
0134     const double Pos_u = stod(vec[id++]);
0135     const double Pos_v = stod(vec[id++]);
0136     const double Dir_X = stod(vec[id++]);
0137     const double Dir_Y = stod(vec[id++]);
0138     const double Dir_Z = stod(vec[id++]);
0139     const double Dir_theta_rad = stod(vec[id++]);
0140     const double Dir_phi_rad = stod(vec[id++]);
0141     const double LowFluxRatio = stod(vec[id++]);
0142     const double Energy_eV = stod(vec[id++]);
0143     const double Flux_photon_s = stod(vec[id++]);
0144     const double Power_W = stod(vec[id++]);
0145     assert(id == 14);
0146
0147     if (Verbosity())
0148     {
0149       cout << "ReadSynRadFiles::process_event - " << line << " -> " << endl
0150            << Pos_X_cm << ", "
0151            << Pos_Y_cm << ", "
0152            << Pos_Z_cm << ", "
0153            << Pos_u << ", "
0154            << Pos_v << ", "
0155            << Dir_X << ", "
0156            << Dir_Y << ", "
0157            << Dir_Z << ", "
0158            << Dir_theta_rad << ", "
0159            << Dir_phi_rad << ", "
0160            << LowFluxRatio << ", "
0161            << Energy_eV << ", "
0162            << Flux_photon_s << ", "
0163            << Power_W << ". "
0164            << endl;
0165     }
0166
0167     double xz_sign = +1;
0168     //assert(m_reverseXZ);
0169     if (m_reverseXZ)
0170     {
0171       xz_sign = -1;
0172
0173       static bool once = true;
0174
0175       if (once)
0176       {
0177         cout << "ReadSynRadFiles::process_event - reverse x z axis direction for input photons" << endl;
0178         once = false;
0179       }
0180     }
0181
0182     const double E_GeV = Energy_eV / 1e9;
0183     //    const double E_GeV = Energy_eV / 1e3;
0184     const double px = E_GeV * Dir_X;
0185     const double py = E_GeV * Dir_Y;
0186     const double pz = E_GeV * Dir_Z;
0187
0188     /* Create HepMC particle record */
0189     HepMC::GenParticle *hepmcpart = new HepMC::GenParticle(HepMC::FourVector(px * xz_sign, py, pz * xz_sign, E_GeV), 22);
0190
0191     hepmcpart->set_status(1);
0192
0193     /* add particle information */
0194     hepmcpart->setGeneratedMass(0);
0195
0196     /* append particle to vector */
0197     hepmc_particles.push_back(hepmcpart);
0198     origin_index.push_back(ii);
0199
0200     HepMC::GenVertex *hepmcvtx = new HepMC::GenVertex(HepMC::FourVector(Pos_X_cm * xz_sign,
0201                                                                         Pos_Y_cm,
0202                                                                         Pos_Z_cm * xz_sign,
0203                                                                         0));
0204     hepmc_vertices.push_back(hepmcvtx);
0205     hepmcvtx->add_particle_out(hepmcpart);
0206
0207     ++SumPhoton;
0208     SumFlux += Flux_photon_s;
0209   }
0210
0211   /* Check if hepmc_particles and origin_index vectors are the same size */
0212   if (hepmc_particles.size() != origin_index.size())
0213   {
0214     cout << "ReadSynRadFiles::process_event - Lengths of HepMC particles and Origin index vectors do not match!" << endl;
0215
0216     delete evt;
0217     return Fun4AllReturnCodes::ABORTRUN;
0218   }
0219
0220   /* Add HepMC vertices to event */
0221   for (unsigned v = 0; v < hepmc_vertices.size(); v++)
0222     evt->add_vertex(hepmc_vertices.at(v));
0223
0224   // save weights
0225   auto &weightcontainer = evt->weights();
0226   weightcontainer.push_back(SumFlux);
0227   weightcontainer.push_back(SumPhoton);
0228
0229   /* pass HepMC to PHNode*/
0230   PHHepMCGenEvent *success = hepmc_helper.insert_event(evt);
0231   if (!success)
0232   {
0233     cout << "ReadSynRadFiles::process_event - Failed to add event to HepMC record!" << endl;
0234     return Fun4AllReturnCodes::ABORTRUN;
0235   }
0236
0237   entry++;
0238   /* Done */
0239   return 0;
0240 }
0241
0242 int ReadSynRadFiles::CreateNodeTree(PHCompositeNode *topNode)
0243 {
0244   hepmc_helper.create_node_tree(topNode);
0245
0246   return Fun4AllReturnCodes::EVENT_OK;
0247 }