estarlight/src/eventfilewriter.cpp

0001
0002 //
0003 //    Copyright 2010
0004 //
0005 //    This file is part of starlight.
0006 //
0007 //    starlight is free software: you can redistribute it and/or modify
0008 //    it under the terms of the GNU General Public License as published by
0009 //    the Free Software Foundation, either version 3 of the License, or
0010 //    (at your option) any later version.
0011 //
0012 //    starlight is distributed in the hope that it will be useful,
0013 //    but WITHOUT ANY WARRANTY; without even the implied warranty of
0014 //    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
0015 //    GNU General Public License for more details.
0016 //
0017 //    You should have received a copy of the GNU General Public License
0018 //    along with starlight. If not, see <http://www.gnu.org/licenses/>.
0019 //
0020 ///////////////////////////////////////////////////////////////////////////
0021 //
0022 // File and Version Information:
0023 // $Rev:: 228                         $: revision of last commit
0024 // $Author:: butter                   $: author of last commit
0025 // $Date:: 2016-01-18 17:10:17 +0000 #$: date of last commit
0026 //
0027 // Description:
0028 //
0029 //
0030 //
0031 ///////////////////////////////////////////////////////////////////////////
0032
0033 #ifdef HEPMC3_ON
0034 #include "hepmc3writer.h"
0035 #endif
0036
0037 #include <iostream>
0038 #include <exception>
0039 #include <cstdlib>
0040
0041 #include "eventfilewriter.h"
0042 #include "starlightparticlecodes.h"
0043
0044 using namespace std;
0045
0046 //______________________________________________________________________________
0047 eventFileWriter::eventFileWriter()
0048 : fileWriter()
0049 ,_writeFullPythia(false)
0050 ,_writeFullHepMC3(false)
0051 { }
0052
0053 //______________________________________________________________________________
0054 eventFileWriter::eventFileWriter(std::string filename)
0055 : fileWriter(filename)
0056 { }
0057
0058 //______________________________________________________________________________
0059 int eventFileWriter::writeInit(inputParameters &_p)
0060 {
0061   _fileStream<<"CONFIG_OPT: "<<_p.productionMode()<<" "<<_p.prodParticleId()<<" "<<_p.nmbEvents()
0062          <<" "<<_p.quantumGlauber()<<" "<<_p.impulseVM()<<" "<<_p.randomSeed()<<std::endl;
0063   _fileStream<<"ELECTRON_BEAM: "<<" "<<_p.electronBeamLorentzGamma()<<std::endl;
0064   _fileStream<<"TARGET_BEAM: "<<_p.targetBeamZ()<<" "<<_p.targetBeamA()<<" "<<_p.targetBeamLorentzGamma()<<std::endl;
0065   _fileStream<<"PHOTON: "<<_p.nmbEnergyBins()<<" "<<_p.fixedQ2Range()<<" "<<_p.nmbGammaQ2Bins()
0066          <<" "<<_p.minGammaQ2()<<" "<<_p.maxGammaQ2()<<std::endl;
0067
0068   if ( _writeFullHepMC3 )
0069     {
0070 #ifdef HEPMC3_ON
0071       _hepmc3writer = new hepMC3Writer();
0072       _hepmc3writer->initWriter(_p);
0073 #else
0074       std::cout << "***HepMC3 file format not written --- eStarlight not compiled with HepMC3 Enabled***" << std::endl;
0075       _writeFullHepMC3 = false;
0076 #endif
0077     }
0078
0079   return 0;
0080 }
0081
0082 //______________________________________________________________________________
0083 int eventFileWriter::writeInitLUND(inputParameters &_p)
0084 {
0085   _fileStream<<"PYTHIA EVENT FILE"<<std::endl;
0086   _fileStream<<"============================================"<<std::endl;
0087   _fileStream<<"I, ievent, genevent, subprocess, nucleon,                 targetparton, xtargparton, beamparton, xbeamparton,               thetabeamprtn, truey, trueQ2, truex, trueW2, trueNu, leptonphi,   s_hat, t_hat, u_hat, pt2_hat, Q2_hat, F2, F1, R, sigma_rad,       SigRadCor, EBrems, photonflux, t-diff, nrTracks"<<std::endl;
0088   _fileStream<<"============================================"<<std::endl;
0089   _fileStream<<"I  K(I,1)  K(I,2)  K(I,3)  K(I,4)  K(I,5)  P(I,1)  P(I,2)  P(I,3)  P(I,4)  P(I,5)  V(I,1)  V(I,2)  V(I,3)"<<std::endl;
0090   _fileStream<<"============================================"<<std::endl;
0091
0092   //Proton and Electron Mass
0093   double p_mass = 0.93827;
0094   double e_mass = 0.000511;
0095
0096   //Calculate Four Vector
0097   double _electronBeam_E  = e_mass*_p.electronBeamLorentzGamma();
0098   double _targetBeam_E  = p_mass*_p.targetBeamLorentzGamma();
0099   double _electronBeam_pz = -1.0*sqrt(_electronBeam_E*_electronBeam_E - e_mass*e_mass);
0100   double _targetBeam_pz = sqrt(_targetBeam_E*_targetBeam_E - p_mass*p_mass);
0101
0102   //Define Four Vector
0103   _electronBeam_four_vector_={0,0,_electronBeam_pz,_electronBeam_E};
0104   _targetBeam_four_vector_ ={0,0,_targetBeam_pz,_targetBeam_E};
0105   _electronBeam_pdg_id_ = 11;
0106   int nuclear_pid = _p.targetBeamA()*10 + _p.targetBeamZ()*10000 + 1000000000;// Form 100ZZZAAAl; l=0
0107   _targetBeam_pdg_id_=nuclear_pid;
0108   // The following line was a temporary fix to make eSTARlight compatible with the ePIC simulations which could not handle the full nuclear_PIDs.  Removed Oct. 20. 2024
0109   //_targetBeam_pdg_id_ = ( _p.targetBeamA() > 1 ) ? nuclear_pid : 2212; //Everything is a proton right now
0110
0111   return 0;
0112 }
0113
0114 //______________________________________________________________________________
0115 int eventFileWriter::writeEvent(eXEvent &event, int eventnumber)
0116 {
0117
0118     int numberoftracks = 0;
0119     if(_writeFullPythia)
0120     {
0121         numberoftracks = event.getParticles()->size();
0122     }
0123     else
0124     {
0125         for(unsigned int i = 0; i<event.getParticles()->size(); ++i)
0126         {
0127             if(event.getParticles()->at(i).getStatus() >= 0) numberoftracks++;
0128         }
0129     }
0130
0131
0132     // sometimes we don't have tracks due to wrongly picked W , check it
0133     if(numberoftracks){
0134       eventnumber++;
0135
0136       _fileStream << "EVENT: " << eventnumber << " " << numberoftracks << " " << 1 << std::endl;
0137
0138       _fileStream <<"VERTEX: "<<0.<<" "<<0.<<" "<<0.<<" "<<0.<<" "<<1<<" "<<0<<" "<<0<<" "<<numberoftracks<<std::endl;
0139
0140       for( uint igam = 0 ; igam < event.getGammaEnergies()->size(); ++igam){
0141     _fileStream <<"GAMMA: "<<event.getGammaEnergies()->at(igam)<<" "<<event.getGammaMasses()->at(igam)<<std::endl;
0142     lorentzVector gam = event.getGamma()->at(igam);
0143     // Write the photon 4-vector out to file. Might be used in later iterations, so I keep it here
0144     //_fileStream <<"GAMMA VECTOR: "<<gam.GetPx()<<" "<<gam.GetPy()<<" "<<gam.GetPz()<<" "<<gam.GetE()<<" "<<-temp<<std::endl;
0145       }
0146       for( uint itarget = 0 ; itarget < event.getTarget()->size(); ++itarget){
0147     lorentzVector target = event.getTarget()->at(itarget);
0148     _fileStream <<"t: "<<event.getVertext()->at(itarget)<<std::endl;
0149     _fileStream <<"TARGET: "<<target.GetPx()<<" "<<target.GetPy()<<" "<<target.GetPz()<<" "<<target.GetE()<<std::endl;
0150       }
0151       for( uint iel = 0 ; iel < event.getSources()->size(); ++iel){
0152     lorentzVector el = event.getSources()->at(iel);
0153     _fileStream <<"SOURCE: "<<el.GetPx()<<" "<<el.GetPy()<<" "<<el.GetPz()<<" "<<el.GetE()<<std::endl;
0154       }
0155
0156       int ipart = 0;
0157       std::vector<starlightParticle>::const_iterator part = (event.getParticles())->begin();
0158
0159       for (part = event.getParticles()->begin(); part != event.getParticles()->end(); part++, ipart++)
0160     {
0161           if(!_writeFullPythia)
0162           {
0163               if((*part).getStatus() < 0) continue;
0164           }
0165       _fileStream << "TRACK: " << " "<< starlightParticleCodes::jetsetToGeant((*part).getPdgCode()) <<" "<< (*part).GetPx() << " " << (*part).GetPy()
0166               << " "<< (*part).GetPz() << " " << eventnumber << " " << ipart << " " << 0 << " "
0167               << (*part).getPdgCode();
0168
0169       if(_writeFullPythia)
0170       {
0171         lorentzVector vtx = (*part).getVertex();
0172         _fileStream << " " << vtx.GetPx() << " " << vtx.GetPy() << " " << vtx.GetPz() << " " << vtx.GetE();
0173         _fileStream << " " << (*part).getFirstParent() << " " << (*part).getLastParent() << " " << (*part).getFirstDaughter() << " " << (*part).getLastDaughter() << " " << (*part).getStatus();
0174       }
0175
0176       _fileStream <<std::endl;
0177     }
0178     }
0179
0180 #ifdef HEPMC3_ON
0181     if ( _writeFullHepMC3 ){
0182       _hepmc3writer->writeEvent(event,eventnumber);
0183     }
0184 #endif
0185
0186     return 0;
0187 }
0188
0189
0190 //______________________________________________________________________________
0191 int eventFileWriter::writeEventLUND(eXEvent &event, inputParameters &pin, int eventnumber)
0192 {
0193     _fileStream << "0         "<<eventnumber<<"         0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0"<<std::endl;
0194     _fileStream<<"============================================"<<std::endl;
0195
0196     int particleIndex=1;
0197     uint lineNumberOfFirstElectronDaughter=3;
0198     uint lineNumberOfLastElectronDaughter=4;
0199     uint lineNumberOfFirstIonDaughter=6;
0200     uint lineNumberOfLastIonDaughter=lineNumberOfFirstIonDaughter+event.getParticles()->size();
0201     uint numberOfElectrons = event.getSources()->size();
0202     uint numberOfIons = event.getTarget()->size();
0203     uint numberOfGammas = event.getGammaEnergies()->size();
0204
0205     auto sp0 = std::setw(12);//column widths.  Change from 8 to 12 on Oct.20, 2024 at Wanchang's request to accommodate heavy nuclei ID's.
0206     auto sp1 = std::setw(14);
0207
0208     double electronmass = sqrt(_electronBeam_four_vector_[3]*_electronBeam_four_vector_[3] - _electronBeam_four_vector_[2]*_electronBeam_four_vector_[2] - _electronBeam_four_vector_[1]*_electronBeam_four_vector_[1] - _electronBeam_four_vector_[0]*_electronBeam_four_vector_[0]);
0209     _fileStream <<particleIndex << sp0 <<21 << sp0 <<_electronBeam_pdg_id_<<"       "<<0<<"       "<<lineNumberOfFirstElectronDaughter<<"       "<<lineNumberOfLastElectronDaughter<<"  " << std::setprecision(6) << sp1 <<_electronBeam_four_vector_[0]<<"  " << sp1 <<_electronBeam_four_vector_[1]<<"  " << sp1 <<_electronBeam_four_vector_[2]<<"  " <<sp1<<_electronBeam_four_vector_[3]<<"   "<<sp1<<electronmass<<"       0       0       0"<<std::endl;
0210     particleIndex++;
0211
0212     double ionmass = sqrt(_targetBeam_four_vector_[3]*_targetBeam_four_vector_[3] - _targetBeam_four_vector_[2]*_targetBeam_four_vector_[2] - _targetBeam_four_vector_[1]*_targetBeam_four_vector_[1] - _targetBeam_four_vector_[0]*_targetBeam_four_vector_[0]);
0213     _fileStream <<particleIndex <<sp0<<21 << sp0 <<_targetBeam_pdg_id_<<"       "<<0<<"       "<<lineNumberOfFirstIonDaughter<<"       "<<lineNumberOfLastIonDaughter<<"  " << sp1 <<_targetBeam_four_vector_[0]<<"  " <<sp1<<_targetBeam_four_vector_[1]<<"  " <<sp1<<_targetBeam_four_vector_[2]<<"  " <<sp1<<_targetBeam_four_vector_[3]<<"   "<<sp1<<ionmass<<"       0       0       0"<<std::endl;
0214     particleIndex++;
0215
0216     //FOR ALEX
0217     for( uint iel = 0 ; iel < numberOfElectrons; ++iel){
0218       lorentzVector el = event.getSources()->at(iel);
0219       _fileStream <<particleIndex<<sp0<<21 <<sp0<<11<<"       "<<1<<"       "<<0<<"       "<<0<<"  " <<sp1<<el.GetPx()<<"  " <<sp1<<el.GetPy()<<"  " <<sp1<<el.GetPz()<<"  " <<sp1<<el.GetE()<<"   "<<sp1<<electronmass<<"       0       0       0"<<std::endl;
0220       particleIndex++;
0221     }
0222
0223     for( uint igam = 0 ; igam < numberOfGammas; ++igam){
0224       lorentzVector gam = event.getGamma()->at(igam);
0225       _fileStream <<particleIndex <<sp0<<21 <<sp0<<22<<"       "<<1<<"       "<<0<<"       "<<0<<"  " <<sp1<<gam.GetPx()<<"  " <<sp1<<gam.GetPy()<<"  " <<sp1<<gam.GetPz()<<"  " <<sp1<<event.getGammaEnergies()->at(igam)<<"   "<<sp1<<event.getGammaMasses()->at(igam)<<"       0       0       0"<<std::endl;
0226       particleIndex++;
0227     }
0228
0229
0230
0231     for( uint iel = 0 ; iel < numberOfElectrons; ++iel){
0232       lorentzVector el = event.getSources()->at(iel);
0233       _fileStream <<particleIndex<<sp0<<1 <<sp0<<11<<"       "<<3<<"       "<<0<<"       "<<0<<"  " <<sp1<<el.GetPx()<<"  " <<sp1<<el.GetPy()<<"  " <<sp1<<el.GetPz()<<"  " <<sp1<<el.GetE()<<"   "<<sp1<<electronmass<<"       0       0       0"<<std::endl;
0234       particleIndex++;
0235     }
0236
0237     //FIXME This is only set up for protons right now
0238      int nuclear_pid = pin.targetBeamA()*10 + pin.targetBeamZ()*10000 + 1000000000;// Form 100ZZZAAAl; l=0
0239     for( uint itarget = 0 ; itarget < numberOfIons; ++itarget){
0240       lorentzVector target = event.getTarget()->at(itarget);
0241       _fileStream <<particleIndex <<sp0<<1 <<sp0<<nuclear_pid <<"       "<<2<<"       "<<lineNumberOfFirstIonDaughter<<"       "<<lineNumberOfLastIonDaughter<<"  " <<sp1<<target.GetPx()<<"  " <<sp1<<target.GetPy()<<"  " <<sp1<<target.GetPz()<<"  " <<sp1<<target.GetE()<<"   "<<sp1<<ionmass<<"       0       0       0"<<std::endl;
0242       particleIndex++;
0243     }
0244
0245
0246     int ipart = 0;
0247     std::vector<starlightParticle>::const_iterator part = (event.getParticles())->begin();
0248
0249     for (part = event.getParticles()->begin(); part != event.getParticles()->end(); part++, ipart++)
0250     {
0251       _fileStream <<particleIndex <<sp0<<1 <<sp0<<(*part).getPdgCode()<<"       "<<2<<"       "<<0<<"       "<<0<<"  " <<sp1<<(*part).GetPx()<<"  " <<sp1<<(*part).GetPy()<<"  " <<sp1<<(*part).GetPz()<<"  " <<sp1<<(*part).GetE()<<"   "<<sp1<<(*part).M()<<"       0       0       0"<<std::endl;
0252       particleIndex++;
0253     }
0254
0255
0256     _fileStream<<"=============== Event finished ==============="<<std::endl;
0257     return 0;
0258 }
0259
0260
0261 int eventFileWriter::close()
0262 {
0263
0264 #ifdef HEPMC3_ON
0265     if ( _writeFullHepMC3 ){
0266       _hepmc3writer->close();
0267     }
0268 #endif
0269
0270     try
0271     {
0272         _fileStream.close();
0273     }
0274     catch (const ios::failure & error)
0275     {
0276         cerr << "I/O exception: " << error.what() << endl;
0277         return EXIT_FAILURE;
0278     }
0279     return 0;
0280 }