File indexing completed on 2026-09-15 08:28:23
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029 #include "SteppingAction.hh"
0030
0031 #include "RunAction.hh"
0032
0033 #include "G4ParticleTypes.hh"
0034 #include "G4SteppingManager.hh"
0035 #include "G4VProcess.hh"
0036
0037
0038
0039 SteppingAction::SteppingAction(RunAction* RuAct) : G4UserSteppingAction(), fRunAction(RuAct)
0040 {
0041 fMuonMass = G4MuonPlus::MuonPlus()->GetPDGMass();
0042 }
0043
0044
0045
0046 SteppingAction::~SteppingAction() {}
0047
0048
0049
0050 void SteppingAction::UserSteppingAction(const G4Step* aStep)
0051 {
0052 const G4VProcess* process = aStep->GetPostStepPoint()->GetProcessDefinedStep();
0053 if (process == 0) return;
0054 G4String processName = process->GetProcessName();
0055 fRunAction->CountProcesses(processName);
0056
0057 if (processName != "GammaToMuPair") return;
0058
0059 G4StepPoint* PrePoint = aStep->GetPreStepPoint();
0060 G4double EGamma = PrePoint->GetTotalEnergy();
0061 G4ThreeVector PGamma = PrePoint->GetMomentum();
0062
0063 G4double Eplus(0), Eminus(0);
0064 G4ThreeVector Pplus, Pminus;
0065 const G4TrackVector* secondary = fpSteppingManager->GetSecondary();
0066 for (size_t lp = 0; lp < (*secondary).size(); lp++) {
0067 if ((*secondary)[lp]->GetDefinition() == G4MuonPlus::MuonPlusDefinition()) {
0068 Eplus = (*secondary)[lp]->GetTotalEnergy();
0069 Pplus = (*secondary)[lp]->GetMomentum();
0070 }
0071 else {
0072 Eminus = (*secondary)[lp]->GetTotalEnergy();
0073 Pminus = (*secondary)[lp]->GetMomentum();
0074 }
0075 }
0076
0077 G4double xPlus = Eplus / EGamma, xMinus = Eminus / EGamma;
0078 G4double thetaPlus = PGamma.angle(Pplus), thetaMinus = PGamma.angle(Pminus);
0079 G4double GammaPlus = Eplus / fMuonMass;
0080 G4double GammaMinus = Eminus / fMuonMass;
0081
0082 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
0083
0084 if (0.0 == thetaPlus || 0.0 == thetaMinus) {
0085 G4cout << "SteppingAction: "
0086 << "thetaPlus= " << thetaPlus << " thetaMinus= " << thetaMinus
0087 << " gamPlus= " << GammaPlus << " gamMinus= " << GammaMinus << " "
0088 << thetaPlus * GammaPlus - thetaMinus * GammaMinus << G4endl;
0089 return;
0090 }
0091 analysisManager->FillH1(1, 1. / (1. + std::pow(thetaPlus * GammaPlus, 2)));
0092 analysisManager->FillH1(2, std::log10(thetaPlus * GammaPlus));
0093
0094 analysisManager->FillH1(3, std::log10(thetaMinus * GammaMinus));
0095 analysisManager->FillH1(4,
0096 std::log10(std::fabs(thetaPlus * GammaPlus - thetaMinus * GammaMinus)));
0097
0098 analysisManager->FillH1(5, xPlus);
0099 analysisManager->FillH1(6, xMinus);
0100 }
0101
0102