Newer
Older
// $Id: Geant4Converter.cpp 603 2013-06-13 21:15:14Z markus.frank $
//====================================================================
// AIDA Detector description implementation for LCD
//--------------------------------------------------------------------
//
// @author P.Kostka (main author)
// @author M.Frank (code reshuffeling into new DDG4 scheme)
//
//====================================================================
// Framework include files
#include "DD4hep/Plugins.h"
#include "DDG4/Geant4Primary.h"
#include "DDG4/Geant4Context.h"
#include "DDG4/Geant4InputAction.h"
#include "G4Event.hh"
using namespace std;
using namespace DD4hep::Simulation;
typedef DD4hep::ReferenceBitMask<int> PropertyMask;
/// Initializing constructor
Geant4EventReader::Geant4EventReader(const std::string& nam)
: m_name(nam), m_directAccess(false), m_currEvent(0)
{
}
/// Default destructor
Geant4EventReader::~Geant4EventReader() {
}
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
/// Skip event. To be implemented for sequential sources
Geant4EventReader::EventReaderStatus Geant4EventReader::skipEvent() {
if ( hasDirectAccess() ) {
++m_currEvent;
return EVENT_READER_OK;
}
std::vector<Particle*> particles;
++m_currEvent;
EventReaderStatus sc = readParticles(m_currEvent,particles);
for_each(particles.begin(),particles.end(),deleteObject<Particle>);
return sc;
}
/// Move to the indicated event number.
Geant4EventReader::EventReaderStatus
Geant4EventReader::moveToEvent(int event_number) {
if ( m_currEvent == event_number-1 ) {
return EVENT_READER_OK;
}
else if ( hasDirectAccess() ) {
m_currEvent = event_number-1;
return EVENT_READER_OK;
}
else if ( event_number<m_currEvent ) {
return EVENT_READER_ERROR;
}
else {
for(int i=m_currEvent; i<event_number;++i)
skipEvent();
m_currEvent = event_number-1;
return EVENT_READER_OK;
}
return EVENT_READER_ERROR;
}
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
/// Standard constructor
Geant4InputAction::Geant4InputAction(Geant4Context* context, const string& name)
: Geant4GeneratorAction(context,name), m_reader(0)
{
declareProperty("Input", m_input);
declareProperty("Sync", m_firstEvent=0);
declareProperty("Mask", m_mask = 0);
declareProperty("MomentumScale", m_momScale = 1.0);
m_needsControl = true;
}
/// Default destructor
Geant4InputAction::~Geant4InputAction() {
}
/// helper to report Geant4 exceptions
string Geant4InputAction::issue(int i) const {
stringstream str;
str << "Geant4InputAction[" << name() << "]: Event " << i << " ";
return str.str();
}
/// Read an event and return a LCCollection of MCParticles.
int Geant4InputAction::readParticles(int evt_number, std::vector<Particle*>& particles) {
int evid = evt_number + m_firstEvent;
if ( 0 == m_reader ) {
if ( m_input.empty() ) {
throw runtime_error("InputAction: No inoput file declared!");
}
string err;
TypeName tn = TypeName::split(m_input,"|");
try {
m_reader = PluginService::Create<Geant4EventReader*>(tn.first,tn.second);
if ( 0 == m_reader ) {
PluginDebug dbg();
m_reader = PluginService::Create<Geant4EventReader*>(tn.first,tn.second);
abortRun(issue(evid)+"Error creating reader plugin.",
"Failed to create file reader of type %s. Cannot open dataset %s",
tn.first.c_str(),tn.second.c_str());
return Geant4EventReader::EVENT_READER_NO_FACTORY;
}
}
catch(const exception& e) {
err = e.what();
}
if ( !err.empty() ) {
abortRun(issue(evid)+err,"Error when creating reader for file %s",m_input.c_str());
return Geant4EventReader::EVENT_READER_NO_FACTORY;
}
}
int status = m_reader->moveToEvent(evid);
if ( Geant4EventReader::EVENT_READER_OK != status ) {
abortRun(issue(evid)+"Error when moving to event - may be end of file.",
"Error when reading file %s",m_input.c_str());
}
status = m_reader->readParticles(evid,particles);
if ( Geant4EventReader::EVENT_READER_OK != status ) {
abortRun(issue(evid)+"Error when reading file - may be end of file.",
"Error when reading file %s",m_input.c_str());
}
return status;
}
/// Callback to generate primary particles
void Geant4InputAction::operator()(G4Event* event) {
vector<Particle*> primaries;
Geant4Event& evt = context()->event();
Geant4PrimaryEvent* prim = evt.extension<Geant4PrimaryEvent>();
Geant4PrimaryInteraction* inter = new Geant4PrimaryInteraction();
int result = readParticles(event->GetEventID(),primaries);
if ( result != Geant4EventReader::EVENT_READER_OK ) { // handle I/O error, but how?
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
return;
}
prim->add(m_mask, inter);
// check if there is at least one particle
if ( primaries.empty() ) return;
print("+++ Particle interaction with %d generator particles ++++++++++++++++++++++++",
int(primaries.size()));
Geant4Vertex* vtx = new Geant4Vertex();
vtx->x = 0;
vtx->y = 0;
vtx->z = 0;
vtx->time = 0;
inter->vertices.insert(make_pair(m_mask,vtx));
// build collection of MCParticles
for(size_t i=0; i<primaries.size(); ++i ) {
Geant4ParticleHandle p(primaries[i]);
const double mom_scale = m_momScale;
PropertyMask status(p->status);
p->psx = mom_scale*p->psx;
p->psy = mom_scale*p->psy;
p->psz = mom_scale*p->psz;
if ( p->parents.size() == 0 ) {
if ( status.isSet(G4PARTICLE_GEN_EMPTY) || status.isSet(G4PARTICLE_GEN_DOCUMENTATION) )
vtx->in.insert(p->id); // Beam particles and primary quarks etc.
else
vtx->out.insert(p->id); // Stuff, to be given to Geant4 together with daughters
}
inter->particles.insert(make_pair(p->id,p));
p.dump3(outputLevel()-1,name(),"+->");
}
}