Skip to content
Snippets Groups Projects
Commit d0b59711 authored by Andrii Verbytskyi's avatar Andrii Verbytskyi Committed by Adam Edward Barton
Browse files

Removing fortran from TruthIO package

parent 3ce49b4a
No related branches found
No related tags found
No related merge requests found
......@@ -22,7 +22,6 @@ find_package( HepPDT )
# Component(s) in the package:
atlas_add_component( TruthIO
src/*.cxx
src/*.F
src/components/*.cxx
INCLUDE_DIRS ${HEPPDT_INCLUDE_DIRS}
LINK_LIBRARIES ${HEPPDT_LIBRARIES} AtlasHepMCLib AtlasHepMCfioLib AthenaBaseComps GaudiKernel GeneratorModulesLib StoreGateLib SGtests EventInfo GeneratorObjects )
......
......@@ -21,5 +21,8 @@ private:
// Setable Properties:-
std::string m_key;
std::string m_input_file;
std::ifstream m_file;
bool read_hepevt_particle( int i);
bool read_hepevt_event_header();
};
/*
Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
*/
extern "C" {
const unsigned int hepevt_bytes_allocation_ATLAS =
sizeof(long int) * ( 2 + 6 * 10000 )
+ sizeof(double) * ( 9 * 10000 );
struct {
char data[hepevt_bytes_allocation_ATLAS];
} hepevt_;
}
\ No newline at end of file
......@@ -14,13 +14,63 @@
#include "StoreGate/StoreGateSvc.h"
extern "C" {
void openfile_(int & ifile, const char* filename, int);
void readevt_(int& ifile);
void closefile_(int& ifile);
}
bool ReadHepEvtFromAscii::read_hepevt_event_header()
{
const size_t max_e_buffer_size=512;
char buf_e[max_e_buffer_size];
bool eventline=false;
int im=0, pm=0;
while(!eventline)
{
m_file.getline(buf_e,max_e_buffer_size);
if( strlen(buf_e) == 0 ) return false;
std::stringstream st_e(buf_e);
char attr=' ';
eventline=false;
while (!eventline)
{
if (!(st_e>>attr)) break;
if (attr==' ') continue;
else eventline=false;
if (attr=='E')
{
eventline=static_cast<bool>(st_e>>im>>pm);
}
}
}
HepMC::HEPEVT_Wrapper::set_event_number(im);
HepMC::HEPEVT_Wrapper::set_number_entries(pm);
return eventline;
}
bool ReadHepEvtFromAscii::read_hepevt_particle( int i)
{
const size_t max_p_buffer_size=512;
const size_t max_v_buffer_size=512;
char buf_p[max_p_buffer_size];
char buf_v[max_v_buffer_size];
int intcodes[6];
double fltcodes1[5];
double fltcodes2[4];
m_file.getline(buf_p,max_p_buffer_size);
if( strlen(buf_p) == 0 ) return false;
m_file.getline(buf_v,max_v_buffer_size);
if( strlen(buf_v) == 0 ) return false;
std::stringstream st_p(buf_p);
std::stringstream st_v(buf_v);
if (!static_cast<bool>(st_p>>intcodes[0]>>intcodes[1]>>intcodes[2]>>intcodes[3]>>intcodes[4]>>intcodes[5]>>fltcodes1[0]>>fltcodes1[1]>>fltcodes1[2]>>fltcodes1[3]>>fltcodes1[4])) { return false;}
if (!static_cast<bool>(st_v>>fltcodes2[0]>>fltcodes2[1]>>fltcodes2[2]>>fltcodes2[3])) { return false;}
HepMC::HEPEVT_Wrapper::set_status(i,intcodes[0]);
HepMC::HEPEVT_Wrapper::set_id(i,intcodes[1]);
HepMC::HEPEVT_Wrapper::set_parents(i,intcodes[2],std::max(intcodes[2],intcodes[3]));/* Pythia writes second mother 0*/
HepMC::HEPEVT_Wrapper::set_children(i,intcodes[4],intcodes[5]);
HepMC::HEPEVT_Wrapper::set_momentum(i,fltcodes1[0],fltcodes1[1],fltcodes1[2],fltcodes1[3]);
HepMC::HEPEVT_Wrapper::set_mass(i,fltcodes1[4]);
HepMC::HEPEVT_Wrapper::set_position(i,fltcodes2[0],fltcodes2[1],fltcodes2[2],fltcodes2[3]);
return true;
}
ReadHepEvtFromAscii::ReadHepEvtFromAscii(const std::string& name, ISvcLocator* pSvcLocator)
: AthAlgorithm(name, pSvcLocator),
m_sgSvc(0)
......@@ -46,10 +96,8 @@ StatusCode ReadHepEvtFromAscii::initialize(){
HepMC::HEPEVT_Wrapper::set_max_number_entries(10000);
// Initialize input file
int ifile=5;
closefile_(ifile);
openfile_(ifile,m_input_file.c_str(),m_input_file.size());
m_file=std::ifstream(m_input_file.c_str());
if( !m_file.is_open() ) return StatusCode::FAILURE;
// Initialization terminated
return StatusCode::SUCCESS;
}
......@@ -73,8 +121,10 @@ StatusCode ReadHepEvtFromAscii::execute() {
}
HepMC::GenEvent* evt = new HepMC::GenEvent();
int ifile=5;
readevt_(ifile);
HepMC::HEPEVT_Wrapper::zero_everything();
bool fileok=read_hepevt_event_header();
for (int i=1; (i<=HepMC::HEPEVT_Wrapper::number_entries())&&fileok; i++) fileok=read_hepevt_particle(i);
if (!fileok) return StatusCode::FAILURE;
HepMC::IO_HEPEVT hepio;
hepio.set_print_inconsistency_errors(0);
......@@ -100,8 +150,7 @@ StatusCode ReadHepEvtFromAscii::finalize() {
msg(MSG::INFO) << ">>> ReadHepEvtFromAscii from finalize" << endmsg;
int ifile=5;
closefile_(ifile);
if(m_file.is_open()) m_file.close();
// End of finalization step
return StatusCode::SUCCESS;
......
subroutine openFile(ifile,filename)
integer ifile
character*(*) filename
logical ex
inquire(FILE=filename,EXIST=ex)
if(.not. ex) write(*,*) 'File does not exist:',filename
open(ifile,FILE=filename,status='old')
end
subroutine readEvt(ifile)
integer ifile
c...Hepevt common block
integer nmxhep,nevhep,nhep,isthep,idhep,jmohep,jdahep
double precision phep,vhep
PARAMETER (NMXHEP=10000)
COMMON/HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
& JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP)
integer i
read(ifile,*) nevhep,nhep
do i=1,nhep
read(ifile,*) dummy,isthep(i),idhep(i),jmohep(1,i),
& jmohep(2,i),jdahep(1,i),jdahep(2,i)
read(ifile,*) phep(1,i),phep(2,i),phep(3,i),phep(4,i),phep(5,i)
read(ifile,*) vhep(1,i),vhep(2,i),vhep(3,i),vhep(4,i)
enddo
end
subroutine closeFile(ifile)
integer ifile
close(ifile)
end
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment