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_ = ( _p.targetBeamA() > 1 ) ? nuclear_pid : 2212; //Everything is a proton right now
0108
0109   return 0;
0110 }
0111
0112 //______________________________________________________________________________
0113 int eventFileWriter::writeEvent(eXEvent &event, int eventnumber)
0114 {
0115
0116     int numberoftracks = 0;
0117     if(_writeFullPythia)
0118     {
0119         numberoftracks = event.getParticles()->size();
0120     }
0121     else
0122     {
0123         for(unsigned int i = 0; i<event.getParticles()->size(); ++i)
0124         {
0125             if(event.getParticles()->at(i).getStatus() >= 0) numberoftracks++;
0126         }
0127     }
0128
0129
0130     // sometimes we don't have tracks due to wrongly picked W , check it
0131     if(numberoftracks){
0132       eventnumber++;
0133
0134       _fileStream << "EVENT: " << eventnumber << " " << numberoftracks << " " << 1 << std::endl;
0135
0136       _fileStream <<"VERTEX: "<<0.<<" "<<0.<<" "<<0.<<" "<<0.<<" "<<1<<" "<<0<<" "<<0<<" "<<numberoftracks<<std::endl;
0137
0138       for( uint igam = 0 ; igam < event.getGammaEnergies()->size(); ++igam){
0139     _fileStream <<"GAMMA: "<<event.getGammaEnergies()->at(igam)<<" "<<event.getGammaMasses()->at(igam)<<std::endl;
0140     lorentzVector gam = event.getGamma()->at(igam);
0141     // Write the photon 4-vector out to file. Might be used in later iterations, so I keep it here
0142     //_fileStream <<"GAMMA VECTOR: "<<gam.GetPx()<<" "<<gam.GetPy()<<" "<<gam.GetPz()<<" "<<gam.GetE()<<" "<<-temp<<std::endl;
0143       }
0144       for( uint itarget = 0 ; itarget < event.getTarget()->size(); ++itarget){
0145     lorentzVector target = event.getTarget()->at(itarget);
0146     _fileStream <<"t: "<<event.getVertext()->at(itarget)<<std::endl;
0147     _fileStream <<"TARGET: "<<target.GetPx()<<" "<<target.GetPy()<<" "<<target.GetPz()<<" "<<target.GetE()<<std::endl;
0148       }
0149       for( uint iel = 0 ; iel < event.getSources()->size(); ++iel){
0150     lorentzVector el = event.getSources()->at(iel);
0151     _fileStream <<"SOURCE: "<<el.GetPx()<<" "<<el.GetPy()<<" "<<el.GetPz()<<" "<<el.GetE()<<std::endl;
0152       }
0153
0154       int ipart = 0;
0155       std::vector<starlightParticle>::const_iterator part = (event.getParticles())->begin();
0156
0157       for (part = event.getParticles()->begin(); part != event.getParticles()->end(); part++, ipart++)
0158     {
0159           if(!_writeFullPythia)
0160           {
0161               if((*part).getStatus() < 0) continue;
0162           }
0163       _fileStream << "TRACK: " << " "<< starlightParticleCodes::jetsetToGeant((*part).getPdgCode()) <<" "<< (*part).GetPx() << " " << (*part).GetPy()
0164               << " "<< (*part).GetPz() << " " << eventnumber << " " << ipart << " " << 0 << " "
0165               << (*part).getPdgCode();
0166
0167       if(_writeFullPythia)
0168       {
0169         lorentzVector vtx = (*part).getVertex();
0170         _fileStream << " " << vtx.GetPx() << " " << vtx.GetPy() << " " << vtx.GetPz() << " " << vtx.GetE();
0171         _fileStream << " " << (*part).getFirstParent() << " " << (*part).getLastParent() << " " << (*part).getFirstDaughter() << " " << (*part).getLastDaughter() << " " << (*part).getStatus();
0172       }
0173
0174       _fileStream <<std::endl;
0175     }
0176     }
0177
0178 #ifdef HEPMC3_ON
0179     if ( _writeFullHepMC3 ){
0180       _hepmc3writer->writeEvent(event,eventnumber);
0181     }
0182 #endif
0183
0184     return 0;
0185 }
0186
0187
0188 //______________________________________________________________________________
0189 int eventFileWriter::writeEventLUND(eXEvent &event, int eventnumber)
0190 {
0191     _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;
0192     _fileStream<<"============================================"<<std::endl;
0193
0194     int particleIndex=1;
0195     uint lineNumberOfFirstElectronDaughter=3;
0196     uint lineNumberOfLastElectronDaughter=4;
0197     uint lineNumberOfFirstIonDaughter=6;
0198     uint lineNumberOfLastIonDaughter=lineNumberOfFirstIonDaughter+event.getParticles()->size();
0199     uint numberOfElectrons = event.getSources()->size();
0200     uint numberOfIons = event.getTarget()->size();
0201     uint numberOfGammas = event.getGammaEnergies()->size();
0202
0203     auto sp0 = std::setw(8);//column widths
0204     auto sp1 = std::setw(14);
0205
0206     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]);
0207     _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;
0208     particleIndex++;
0209
0210     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]);
0211     _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;
0212     particleIndex++;
0213
0214     //FOR ALEX
0215     for( uint iel = 0 ; iel < numberOfElectrons; ++iel){
0216       lorentzVector el = event.getSources()->at(iel);
0217       _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;
0218       particleIndex++;
0219     }
0220
0221     for( uint igam = 0 ; igam < numberOfGammas; ++igam){
0222       lorentzVector gam = event.getGamma()->at(igam);
0223       _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;
0224       particleIndex++;
0225     }
0226
0227
0228
0229     for( uint iel = 0 ; iel < numberOfElectrons; ++iel){
0230       lorentzVector el = event.getSources()->at(iel);
0231       _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;
0232       particleIndex++;
0233     }
0234
0235     //FIXME This is only set up for protons right now
0236     for( uint itarget = 0 ; itarget < numberOfIons; ++itarget){
0237       lorentzVector target = event.getTarget()->at(itarget);
0238       _fileStream <<particleIndex <<sp0<<1 <<sp0<<2212<<"       "<<2<<"       "<<lineNumberOfFirstIonDaughter<<"       "<<lineNumberOfLastIonDaughter<<"  " <<sp1<<target.GetPx()<<"  " <<sp1<<target.GetPy()<<"  " <<sp1<<target.GetPz()<<"  " <<sp1<<target.GetE()<<"   "<<sp1<<ionmass<<"       0       0       0"<<std::endl;
0239       particleIndex++;
0240     }
0241
0242
0243     int ipart = 0;
0244     std::vector<starlightParticle>::const_iterator part = (event.getParticles())->begin();
0245
0246     for (part = event.getParticles()->begin(); part != event.getParticles()->end(); part++, ipart++)
0247     {
0248       _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;
0249       particleIndex++;
0250     }
0251
0252
0253     _fileStream<<"=============== Event finished ==============="<<std::endl;
0254     return 0;
0255 }
0256
0257
0258 int eventFileWriter::close()
0259 {
0260
0261 #ifdef HEPMC3_ON
0262     if ( _writeFullHepMC3 ){
0263       _hepmc3writer->close();
0264     }
0265 #endif
0266
0267     try
0268     {
0269         _fileStream.close();
0270     }
0271     catch (const ios::failure & error)
0272     {
0273         cerr << "I/O exception: " << error.what() << endl;
0274         return EXIT_FAILURE;
0275     }
0276     return 0;
0277 }