11 #include "GaudiKernel/IDataProviderSvc.h"
12 #include "GaudiKernel/INTupleSvc.h"
13 #include "GaudiKernel/NTuple.h"
14 #include "GaudiKernel/SmartDataPtr.h"
15 #include "GaudiKernel/MsgStream.h"
16 #include "GaudiKernel/Bootstrap.h"
17 #include "GaudiKernel/ISvcLocator.h"
18 #include "GaudiKernel/IMessageSvc.h"
22 #include "G4Material.hh"
23 #include "G4Element.hh"
24 #include "G4StepPoint.hh"
25 #include "G4TouchableHistory.hh"
26 #include "G4LogicalVolume.hh"
34 AthMessaging(
Gaudi::svcLocator()->service< IMessageSvc >(
"MessageSvc" ),
"MaterialStepRecorder"),
35 m_evtStore(
"StoreGateSvc/StoreGateSvc",
"MaterialStepRecorder"),
36 m_matStepCollection(nullptr),
37 m_matStepCollectionName(
"MaterialStepRecords"),
38 m_recordComposition(true),
41 m_elementTable(nullptr),
42 m_elementTableName(
"ElementTable"),
43 m_runElementTable(nullptr)
69 ATH_MSG_ERROR(
"EndOfEventAction : recording ElementTable in StoreGate FAILED");
87 if (aStep->GetTrack()->GetParentID()) {
88 aStep->GetTrack()->SetTrackStatus(fStopAndKill);
99 const G4TouchableHistory* touchHist =
static_cast<const G4TouchableHistory*
>(aStep->GetPreStepPoint()->GetTouchable());
101 const G4LogicalVolume *lv= touchHist ? touchHist->GetVolume()->GetLogicalVolume() :
nullptr;
102 const G4Material *
mat = lv ? lv->GetMaterial() :
nullptr;
104 std::vector<unsigned char> elements;
105 std::vector<unsigned char> fractions;
108 if (
mat &&
mat->GetRadlen() < 200000.) {
113 double steplength = aStep->GetStepLength();
116 G4ThreeVector
pos = aStep->GetPreStepPoint()->GetPosition();
118 double X0 =
mat->GetRadlen();
119 double L0 =
mat->GetNuclearInterLength();
127 size_t elNumber =
mat->GetNumberOfElements();
128 const G4ElementVector* elVector =
mat->GetElementVector();
129 double totAtoms =
mat->GetTotNbOfAtomsPerVolume();
132 elements.reserve(elNumber);
133 fractions.reserve(elNumber);
140 unsigned int Zint = (
unsigned int)Z;
142 elements.push_back(
static_cast<unsigned char>(Z));
143 fractions.push_back(UCHAR_MAX);
148 G4String elName = (*elVector)[0]->GetName();
157 const G4double* atVector =
mat->GetVecNbOfAtomsPerVolume();
158 double totalFrac = 0.;
160 for (
size_t iel = 0; iel < elNumber; ++iel) {
162 const G4Element* currentEl = (*elVector)[iel];
163 double currentNum = atVector ? atVector[iel] : 1.;
164 double relNbAtoms = currentNum/totAtoms;
166 double currentZ = currentEl->GetZ();
169 Z += relNbAtoms*currentEl->GetZ();
170 unsigned int Zint = (
unsigned int)(currentEl->GetZ());
173 elements.push_back(
int(currentZ));
175 unsigned int relNbAtomsChar = (
unsigned int)(relNbAtoms*(UCHAR_MAX));
176 relNbAtomsChar = relNbAtomsChar > UCHAR_MAX ? UCHAR_MAX : relNbAtomsChar;
179 totalFrac += relNbAtoms;
180 if (relNbAtomsChar) {
181 fractions.push_back(relNbAtomsChar);
185 double curZ = currentEl->GetZ();
189 double curX0 = 1432.8*curA/(curZ*(curZ+1)*(11.319-
std::log(curZ)));
191 double curRho =
rho*relNbAtoms;
193 const G4String& elName = currentEl->GetName();
200 if ((totalFrac-1.)*(totalFrac-1.) > 10
e-4 )
201 ATH_MSG_DEBUG(
"Total fractions do not add up to one at INPUT (" << totalFrac <<
") !");
205 if (aStep->GetTrack()->GetParticleDefinition()->GetPDGEncoding() == 0) {
207 m_matStepCollection->
push_back(
new Trk::MaterialStep(
pos.x(),
pos.y(),
pos.z(), steplength,
X0,
L0,
A, Z,
rho, elements, fractions));