| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 /// \file SteppingAction.cc
27 /// \brief Implementation of the SteppingActio
28 //
29 //
30 //....oooOO0OOooo........oooOO0OOooo........oo
31 //....oooOO0OOooo........oooOO0OOooo........oo
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........oo
43
44 // SteppingAction::SteppingAction()
45 //: G4UserSteppingAction()
46 //{ }
47
48 //....oooOO0OOooo........oooOO0OOooo........oo
49
50 void SteppingAction::UserSteppingAction(const
51 {
52 // check trackID and stepNumber
53 G4int trackID = aStep->GetTrack()->GetTrackI
54 G4int stepNb = aStep->GetTrack()->GetCurrent
55 if (trackID * stepNb != 1) return;
56 // ok, we are at first interaction of the pr
57
58 Run* run = static_cast<Run*>(G4RunManager::G
59
60 // count processes
61 //
62 const G4StepPoint* endPoint = aStep->GetPost
63 G4VProcess* process = const_cast<G4VProcess*
64 run->CountProcesses(process);
65
66 // check that an real interaction occured (e
67 G4StepStatus stepStatus = endPoint->GetStepS
68 G4bool transmit = (stepStatus == fGeomBounda
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->GetPreS
79 G4double Q = -prePoint->GetKineticEnergy();
80 G4ThreeVector Pbalance = -prePoint->GetMomen
81
82 // initialisation of the nuclear channel ide
83 //
84 G4ParticleDefinition* particle = aStep->GetT
85 G4String partName = particle->GetParticleNam
86 G4String nuclearChannel = partName;
87 G4HadronicProcess* hproc = dynamic_cast<G4Ha
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(
94
95 // scattered primary particle (if any)
96 //
97 G4AnalysisManager* analysis = G4AnalysisMana
98 G4int ih = 1;
99 if (aStep->GetTrack()->GetTrackStatus() == f
100 G4double energy = endPoint->GetKineticEner
101 analysis->FillH1(ih, energy);
102 //
103 G4ThreeVector momentum = endPoint->GetMome
104 Q += energy;
105 Pbalance += momentum;
106 //
107 nuclearChannel += partName + " + ";
108 }
109
110 // secondaries
111 //
112 const std::vector<const G4Track*>* secondary
113 for (size_t lp = 0; lp < (*secondary).size()
114 particle = (*secondary)[lp]->GetDefinition
115 G4String name = particle->GetParticleName(
116 G4String type = particle->GetParticleType(
117 G4double energy = (*secondary)[lp]->GetKin
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]-
147 Q += energy;
148 Pbalance += momentum;
149 // count e- from internal conversion toget
150 if (particle == G4Electron::Electron()) pa
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 "
166 "9 ", "10 ", "11
167 std::map<G4ParticleDefinition*, G4int>::iter
168 for (ip = fParticleFlag.begin(); ip != fPart
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()) nuclearCh
180 nuclearChannel += Nb + name;
181 }
182
183 /// G4cout << "\n nuclear channel: " << nucl
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........oo
194