Skip to content
Snippets Groups Projects
Commit 2bc9beed authored by Adam Edward Barton's avatar Adam Edward Barton
Browse files

Merge branch 'truthio_cleanup' into 'master'

Removing fortran from TruthIO package

See merge request atlas/athena!32990
parents a0fa1267 d0b59711
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