202 {
203
204
205
206 const G4VProcess *
process = aStep->GetPostStepPoint()->GetProcessDefinedStep();
207
208
209 G4StepPoint *preStep=aStep->GetPreStepPoint();
210 G4StepPoint *postStep=aStep->GetPostStepPoint();
211
212 G4ThreeVector
mom = preStep->GetMomentum();
213 const G4ThreeVector&
pos = preStep->GetPosition();
214
215
216 G4Track *
track = aStep->GetTrack();
217
218 int trackID=
track->GetTrackID();
219 const bool trackIsAlive = (
track->GetTrackStatus() == fAlive ||
track->GetTrackStatus() == fStopButAlive);
220
222
223 m_pdg =
track->GetDefinition()->GetPDGEncoding();
224 const G4VProcess* creation =
track->GetCreatorProcess();
225 m_scIn = creation? creation->GetProcessSubType() : -1;
226
227 VTrackInformation * trackInfo=
static_cast<VTrackInformation*
>(
track->GetUserInformation());
229 HepMC::GenVertexPtr vtx = currentGenParticle ? currentGenParticle->production_vertex() :
nullptr;
230 m_gen = currentGenParticle? 0 : -1;
231
232 if (currentGenParticle) {
233 while (currentGenParticle && vtx ) {
234 int pdgID=currentGenParticle->pdg_id();
235 const HepMC::GenParticlePtr genmom = vtx->particles_in().size()>0 ? vtx->particles_in().front() :
nullptr;
236 if ( genmom && pdgID!=genmom->pdg_id() )
m_gen++;
237 else if (vtx->particles_out().size()>0 && currentGenParticle!=vtx->particles_out().front())
m_gen++;
238 vtx = genmom ? genmom->production_vertex() : nullptr;
239 currentGenParticle = genmom;
240 }
241 } else {
242
243 int parentID=
track->GetParentID();
244 std::map<int, int>::iterator genIt =
m_trackGenMap.find(parentID);
246 }
247
249
251
253 }
254
261
262
263
264 double stepLength = aStep->GetStepLength();
265 double radLengthInX0 = preStep->GetMaterial()->GetRadlen();
266 double l0 = preStep->GetMaterial()->GetNuclearInterLength();
267 float stepInX0 = stepLength/radLengthInX0;
268
269 if (stepInX0>1.e-06) {
270
273
274 if (l0>0.) {
276
277
278 double totNA = preStep->GetMaterial()->GetTotNbOfAtomsPerVolume();
279 const G4ElementVector* eVec = preStep->GetMaterial()->GetElementVector();
280 const G4double* atVector = preStep->GetMaterial() ->GetVecNbOfAtomsPerVolume();
281 float mFactor =stepInX0* preStep->GetMaterial()->GetDensity();
282
283 float zOverA = 0.; float frSum = 0.;
284 for (
unsigned int i=0;
i<eVec->size();
i++) {
285 float fEl = atVector ? atVector[
i]/totNA : 0.;
286 m_wZ += stepInX0*fEl*((*eVec)[
i]->GetZ());
287
288
289 zOverA += fEl*((*eVec)[
i]->GetZ())/((*eVec)[
i]->GetA());
290 frSum += fEl;
291 }
292 if (fabs(frSum-1.)>0.01)
ATH_MSG_DEBUG(
"G4 material description inconsistent, sum of element fractions:"<< frSum);
294 }
295
296 }
297
298
299
301
302 float eloss = postStep->GetMomentum().mag()-preStep->GetMomentum().mag();
303
306
307 ::iGeant4::Geant4TruthIncident truth( aStep, geoID);
308 unsigned int nSec = truth.numberOfChildren();
309 if (nSec>0 || !trackIsAlive ) {
310
314 G4ThreeVector
mom = preStep->GetMomentum();
316
319 m_vtx_phi = postStep->GetPosition().phi();
320
321 int iPrimSurv = !trackIsAlive ? 0 : 1;
323
324 G4ThreeVector pbal(mom);
325
326 if (iPrimSurv>0) {
330 pbal -= postStep->GetMomentum();
331 }
332
334 for (unsigned int isec=0; isec< nSec; isec++) {
335 G4ThreeVector secMom = truth.childP(isec);
336 if (isec<nSecMax) {
337 m_pdg_child[isec+iPrimSurv] = truth.childPdgCode(isec);
340 }
341 pbal -= secMom;
342 }
343
347
349
350
352 }
353 }
354
355
356 G4VPhysicalVolume *preVol=preStep->GetPhysicalVolume();
357 G4VPhysicalVolume *postVol=postStep->GetPhysicalVolume();
358
359 if (postVol==0) {
360
364
365 m_pth =
track->GetVertexMomentumDirection().theta();
366 m_pph =
track->GetVertexMomentumDirection().phi();
367 m_p =
track->GetVertexKineticEnergy();
369
373
374 m_thEnd=postStep->GetPosition().theta();
375 m_phEnd=postStep->GetPosition().phi();
376 m_dEnd=postStep->GetPosition().mag();
377
383
386
387 return;
388 }
389
390
391 if ( !trackIsAlive ) {
395
396 m_pth =
track->GetVertexMomentumDirection().theta();
397 m_pph =
track->GetVertexMomentumDirection().phi();
398 m_p =
track->GetVertexKineticEnergy();
400
404
405 m_thEnd=postStep->GetPosition().theta();
406 m_phEnd=postStep->GetPosition().phi();
407 m_dEnd=postStep->GetPosition().mag();
408
417 }
418
419 if ( preVol==postVol ) return;
420
421
422 const G4ThreeVector &postPos = postStep->GetPosition();
423 const G4ThreeVector &postMom = postStep->GetMomentum();
425 postPos.y(),
426 postPos.z(),
427 postMom.x(),
428 postMom.y(),
429 postMom.z() );
430
431
432 if ( nextGeoID == geoID) {
434 } else {
435 ATH_MSG_DEBUG(
"track moves from "<<geoID<<
" to "<<nextGeoID);
436
437
438
439 if (aStep->GetStepLength() == 0.) {
440 return;
441 }
442
443
445
446
447
448 const G4ThreeVector& g4pos =
track->GetPosition();
449
450 const HepGeom::Point3D<double>
position(g4pos.x(),g4pos.y(),g4pos.z());
451
452 G4ThreeVector g4mom =
track->GetMomentum();
453 const HepGeom::Vector3D<double>
momentum(g4mom.x(),g4mom.y(),g4mom.z());
454
455 bool dead=false;
456 if (
track->GetTrackStatus()==fStopAndKill) {
457 dead=true;
458 }
459
460 if (!dead) {
461
462
463
466
467 m_pth =
track->GetVertexMomentumDirection().theta();
468 m_pph =
track->GetVertexMomentumDirection().phi();
469 m_p =
track->GetVertexKineticEnergy();
471
475
476 m_thEnd=postStep->GetPosition().theta();
477 m_phEnd=postStep->GetPosition().phi();
478 m_dEnd=postStep->GetPosition().mag();
479
487 }
489
490 }
491
492 }
const std::string process
AtlasRegion
A simple enum of ATLAS regions and sub-detectors.
HepMC3::GenParticlePtr GenParticlePtr
HepMC3::GenVertexPtr GenVertexPtr
const Amg::Vector3D & position() const
Method to retrieve the position of the Intersection.