| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 /// \file SteppingAction.cc
27 /// \brief Implementation of the SteppingAction class
28 //
29 //
30 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
32
33 #include "SteppingAction.hh"
34
35 #include "HistoManager.hh"
36 #include "Run.hh"
37
38 #include "G4HadronicProcess.hh"
39 #include "G4ParticleTypes.hh"
40 #include "G4RunManager.hh"
41
42 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
43
44 // SteppingAction::SteppingAction()
45 //: G4UserSteppingAction()
46 //{ }
47
48 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
49
50 void SteppingAction::UserSteppingAction(const G4Step* aStep)
51 {
52 // check trackID and stepNumber
53 G4int trackID = aStep->GetTrack()->GetTrackID();
54 G4int stepNb = aStep->GetTrack()->GetCurrentStepNumber();
55 if (trackID * stepNb != 1) return;
56 // ok, we are at first interaction of the primary particle
57
58 Run* run = static_cast<Run*>(G4RunManager::GetRunManager()->GetNonConstCurrentRun());
59
60 // count processes
61 //
62 const G4StepPoint* endPoint = aStep->GetPostStepPoint();
63 G4VProcess* process = const_cast<G4VProcess*>(endPoint->GetProcessDefinedStep());
64 run->CountProcesses(process);
65
66 // check that an real interaction occured (eg. not a transportation)
67 G4StepStatus stepStatus = endPoint->GetStepStatus();
68 G4bool transmit = (stepStatus == fGeomBoundary || stepStatus == fWorldBoundary);
69 if (transmit) return;
70
71 // real processes : sum track length
72 //
73 G4double stepLength = aStep->GetStepLength();
74 run->SumTrack(stepLength);
75
76 // energy-momentum balance initialisation
77 //
78 const G4StepPoint* prePoint = aStep->GetPreStepPoint();
79 G4double Q = -prePoint->GetKineticEnergy();
80 G4ThreeVector Pbalance = -prePoint->GetMomentum();
81
82 // initialisation of the nuclear channel identification
83 //
84 G4ParticleDefinition* particle = aStep->GetTrack()->GetDefinition();
85 G4String partName = particle->GetParticleName();
86 G4String nuclearChannel = partName;
87 G4HadronicProcess* hproc = dynamic_cast<G4HadronicProcess*>(process);
88 const G4Isotope* target = NULL;
89 if (hproc) target = hproc->GetTargetIsotope();
90 G4String targetName = "XXXX";
91 if (target) targetName = target->GetName();
92 nuclearChannel += " + " + targetName + " --> ";
93 if (targetName == "XXXX") run->SetTargetXXX(true);
94
95 // scattered primary particle (if any)
96 //
97 G4AnalysisManager* analysis = G4AnalysisManager::Instance();
98 G4int ih = 1;
99 if (aStep->GetTrack()->GetTrackStatus() == fAlive) {
100 G4double energy = endPoint->GetKineticEnergy();
101 analysis->FillH1(ih, energy);
102 //
103 G4ThreeVector momentum = endPoint->GetMomentum();
104 Q += energy;
105 Pbalance += momentum;
106 //
107 nuclearChannel += partName + " + ";
108 }
109
110 // secondaries
111 //
112 const std::vector<const G4Track*>* secondary = aStep->GetSecondaryInCurrentStep();
113 for (size_t lp = 0; lp < (*secondary).size(); lp++) {
114 particle = (*secondary)[lp]->GetDefinition();
115 G4String name = particle->GetParticleName();
116 G4String type = particle->GetParticleType();
117 G4double energy = (*secondary)[lp]->GetKineticEnergy();
118 run->ParticleCount(name, energy);
119 // energy spectrum
120 ih = 0;
121 if (particle == G4Gamma::Gamma())
122 ih = 2;
123 else if (particle == G4Electron::Electron())
124 ih = 3;
125 else if (particle == G4Neutron::Neutron())
126 ih = 4;
127 else if (particle == G4Proton::Proton())
128 ih = 5;
129 else if (particle == G4Deuteron::Deuteron())
130 ih = 6;
131 else if (particle == G4Alpha::Alpha())
132 ih = 7;
133 else if (type == "nucleus")
134 ih = 8;
135 else if (type == "meson")
136 ih = 9;
137 else if (type == "baryon")
138 ih = 10;
139 if (ih > 0) analysis->FillH1(ih, energy);
140 // atomic mass
141 if (type == "nucleus") {
142 G4int A = particle->GetAtomicMass();
143 analysis->FillH1(13, A);
144 }
145 // energy-momentum balance
146 G4ThreeVector momentum = (*secondary)[lp]->GetMomentum();
147 Q += energy;
148 Pbalance += momentum;
149 // count e- from internal conversion together with gamma
150 if (particle == G4Electron::Electron()) particle = G4Gamma::Gamma();
151 // particle flag
152 fParticleFlag[particle]++;
153 }
154
155 // energy-momentum balance
156 G4double Pbal = Pbalance.mag();
157 run->Balance(Pbal);
158 ih = 11;
159 analysis->FillH1(ih, Q);
160 ih = 12;
161 analysis->FillH1(ih, Pbal);
162
163 // nuclear channel
164 const G4int kMax = 16;
165 const G4String conver[] = {"0", "", "2 ", "3 ", "4 ", "5 ", "6 ", "7 ", "8 ",
166 "9 ", "10 ", "11 ", "12 ", "13 ", "14 ", "15 ", "16 "};
167 std::map<G4ParticleDefinition*, G4int>::iterator ip;
168 for (ip = fParticleFlag.begin(); ip != fParticleFlag.end(); ip++) {
169 particle = ip->first;
170 G4String name = particle->GetParticleName();
171 G4int nb = ip->second;
172 if (nb > kMax) nb = kMax;
173 G4String Nb = conver[nb];
174 if (particle == G4Gamma::Gamma()) {
175 run->CountGamma(nb);
176 Nb = "N ";
177 name = "gamma or e-";
178 }
179 if (ip != fParticleFlag.begin()) nuclearChannel += " + ";
180 nuclearChannel += Nb + name;
181 }
182
183 /// G4cout << "\n nuclear channel: " << nuclearChannel << G4endl;
184 run->CountNuclearChannel(nuclearChannel, Q);
185
186 fParticleFlag.clear();
187
188 // kill event after first interaction
189 //
190 G4RunManager::GetRunManager()->AbortEvent();
191 }
192
193 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
194