Hello,
I’m trying to simulate electron ionization in a target volume (icsd example)…
How can i get the direction of secondary electrons from SteppingAction?
I know that i can get the initial direction of primary electrons from : preStep->GetMomentumDirection()
and the final direction of primary electrons from: postStep->GetMomentumDirection()
But how can i find the direction of secondary electrons?
I tried this but it didn’t work and the result is different from the direction set in the physics model class:
G4double E= preStep->GetKineticEnergy();
G4double Ep= postStep->GetKineticEnergy();
G4double totalMomentum = std::sqrt(E*(E + 2.*CLHEP::electron_mass_c2 ));
G4double totalMomentump = std::sqrt(Ep*(Ep + 2.*CLHEP::electron_mass_c2 ));
G4ThreeVector primaryDirection = preStep->GetMomentumDirection();
G4ThreeVector finalDirection = postStep->GetMomentumDirection();
G4double secPx = totalMomentum*primaryDirection.x() - totalMomentump*finalDirection.x();
G4double secPy = totalMomentum*primaryDirection.y() - totalMomentump*finalDirection.y();
G4double secPz = totalMomentum*primaryDirection.z() - totalMomentump*finalDirection.z();
G4double secMomentum = std::sqrt(secPx*secPx + secPy*secPy + secPz*secPz);
secPx /= secMomentum;
secPy /= secMomentum;
secPz /= secMomentum;
G4ThreeVector secDirection;
secDirection.set(secPx,secPy,secPz);
In the Model Class the final direction of primary electron is calculated as follows:
// final_direction_primary_electrons (scattered electrons)
// deltaDirection is the direction of ejected electron
G4ParticleMomentum primaryDirection = particle->GetMomentumDirection();
G4double particleMass = particle->GetDefinition()->GetPDGMass();
G4double totalEnergy = k + particleMass;
G4double pSquare = k * (totalEnergy + particleMass);
G4double totalMomentum = std::sqrt(pSquare);
G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*secondaryKinetic + 2.*electron_mass_c2 );
G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z();
G4double finalMomentum = std::sqrt(finalPx*finalPx + finalPy*finalPy + finalPz*finalPz);
finalPx /= finalMomentum;
finalPy /= finalMomentum;
finalPz /= finalMomentum;
G4ThreeVector direction;
direction.set(finalPx,finalPy,finalPz);