85 {
86
87 if (aStep->GetTrack()->GetParentID()) {
88 aStep->GetTrack()->SetTrackStatus(fStopAndKill);
89 return;
90 }
91
92
96 }
97
98
99 const G4TouchableHistory* touchHist = static_cast<const G4TouchableHistory*>(aStep->GetPreStepPoint()->GetTouchable());
100
101 const G4LogicalVolume *lv= touchHist ? touchHist->GetVolume()->GetLogicalVolume() : nullptr;
102 const G4Material *
mat = lv ? lv->GetMaterial() :
nullptr;
103
104 std::vector<unsigned char> elements;
105 std::vector<unsigned char> fractions;
106
107
108 if (mat &&
mat->GetRadlen() < 200000.) {
109
111
112
113 double steplength = aStep->GetStepLength();
114
115
116 G4ThreeVector
pos = aStep->GetPreStepPoint()->GetPosition();
117
118 double X0 =
mat->GetRadlen();
119 double L0 =
mat->GetNuclearInterLength();
122
123
124 double rho =
mat->GetDensity()*CLHEP::mm3/CLHEP::gram;
125
126
127 size_t elNumber =
mat->GetNumberOfElements();
128 const G4ElementVector* elVector =
mat->GetElementVector();
129 double totAtoms =
mat->GetTotNbOfAtomsPerVolume();
130
131
132 elements.reserve(elNumber);
133 fractions.reserve(elNumber);
134
135 if (1 == elNumber) {
136
137 A =
mat->GetA()/CLHEP::gram;
139
140 unsigned int Zint = (
unsigned int)Z;
141
142 elements.push_back(static_cast<unsigned char>(Z));
143 fractions.push_back(UCHAR_MAX);
144
146
147 Trk::Material elMat(X0,L0,A,Z,rho);
148 G4String elName = (*elVector)[0]->GetName();
149
152 }
153
154
155 } else {
156
157 const G4double* atVector =
mat->GetVecNbOfAtomsPerVolume();
158 double totalFrac = 0.;
159
160 for (size_t iel = 0; iel < elNumber; ++iel) {
161
162 const G4Element* currentEl = (*elVector)[iel];
163 double currentNum = atVector ? atVector[iel] : 1.;
164 double relNbAtoms = currentNum/totAtoms;
165
166 double currentZ = currentEl->GetZ();
167
168 A += relNbAtoms*currentEl->GetA()/CLHEP::gram;
169 Z += relNbAtoms*currentEl->GetZ();
170 unsigned int Zint = (
unsigned int)(currentEl->GetZ());
171
172
173 elements.push_back(int(currentZ));
174
175 unsigned int relNbAtomsChar = (
unsigned int)(relNbAtoms*(UCHAR_MAX));
176 relNbAtomsChar = relNbAtomsChar > UCHAR_MAX ? UCHAR_MAX : relNbAtomsChar;
177
178
179 totalFrac += relNbAtoms;
180 if (relNbAtomsChar) {
181 fractions.push_back(relNbAtomsChar);
182
184 double curA = currentEl->GetA()/CLHEP::gram;
185 double curZ = currentEl->GetZ();
186
187
188
189 double curX0 = 1432.8*curA/(curZ*(curZ+1)*(11.319-std::log(curZ)));
190 double curL0 = 0.;
191 double curRho =
rho*relNbAtoms;
192 Trk::Material elMat(curX0,curL0,curA,curZ,curRho);
193 const G4String& elName = currentEl->GetName();
194
197 }
198 }
199 }
200 if ((totalFrac-1.)*(totalFrac-1.) > 10e-4 )
201 ATH_MSG_DEBUG(
"Total fractions do not add up to one at INPUT (" << totalFrac <<
") !");
202 }
203
204
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));
208 else
210 }
211 }
212 }
1.14.0