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