| 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 // Author: F. Poignant, floriane.poignant@gma 26 // Author: F. Poignant, floriane.poignant@gmail.com
27 // 27 //
28 // file STCyclotronSensitiveTarget.cc 28 // file STCyclotronSensitiveTarget.cc
29 // 29 //
30 #include "STCyclotronRun.hh" 30 #include "STCyclotronRun.hh"
31 #include "STCyclotronSensitiveTarget.hh" 31 #include "STCyclotronSensitiveTarget.hh"
32 32
33 #include "G4AnalysisManager.hh" 33 #include "G4AnalysisManager.hh"
34 #include "G4RunManager.hh" 34 #include "G4RunManager.hh"
35 #include "G4HCofThisEvent.hh" 35 #include "G4HCofThisEvent.hh"
36 #include "G4UnitsTable.hh" 36 #include "G4UnitsTable.hh"
37 #include "G4Step.hh" 37 #include "G4Step.hh"
38 #include "G4SteppingManager.hh" 38 #include "G4SteppingManager.hh"
39 #include "G4ThreeVector.hh" 39 #include "G4ThreeVector.hh"
40 #include "G4SDManager.hh" 40 #include "G4SDManager.hh"
41 #include "G4ios.hh" 41 #include "G4ios.hh"
42 #include "G4SystemOfUnits.hh" 42 #include "G4SystemOfUnits.hh"
43 #include "G4ThreeVector.hh" 43 #include "G4ThreeVector.hh"
44 #include "G4Track.hh" 44 #include "G4Track.hh"
45 #include "G4ParticleDefinition.hh" 45 #include "G4ParticleDefinition.hh"
46 #include "G4DecayTable.hh" 46 #include "G4DecayTable.hh"
47 #include "G4VDecayChannel.hh" 47 #include "G4VDecayChannel.hh"
48 #include "G4TrackVector.hh" 48 #include "G4TrackVector.hh"
49 #include "G4VProcess.hh" 49 #include "G4VProcess.hh"
50 #include "G4Tubs.hh" 50 #include "G4Tubs.hh"
51 #include <map> 51 #include <map>
52 52
53 STCyclotronSensitiveTarget::STCyclotronSensiti 53 STCyclotronSensitiveTarget::STCyclotronSensitiveTarget(G4String name,
54 STCyclotronDetectorConstruc 54 STCyclotronDetectorConstruction* det)
55 : G4VSensitiveDetector(name), 55 : G4VSensitiveDetector(name),
56 fDet(det) 56 fDet(det)
57 { 57 {
58 fTempTrack = 0; 58 fTempTrack = 0;
59 fTempTrack1 = 0; 59 fTempTrack1 = 0;
60 fTempEnergy = 0.; 60 fTempEnergy = 0.;
61 fTempVector = G4ThreeVector(0.,0.,0.); 61 fTempVector = G4ThreeVector(0.,0.,0.);
62 fTrack=0; 62 fTrack=0;
63 } 63 }
64 64
65 STCyclotronSensitiveTarget::~STCyclotronSensit 65 STCyclotronSensitiveTarget::~STCyclotronSensitiveTarget()
66 { 66 {
67 delete fTrack; 67 delete fTrack;
68 } 68 }
69 69
70 G4bool STCyclotronSensitiveTarget::ProcessHits 70 G4bool STCyclotronSensitiveTarget::ProcessHits(G4Step* aStep, G4TouchableHistory*)
71 { 71 {
72 72
73 STCyclotronRun* fRun = static_cast<STCyclotr 73 STCyclotronRun* fRun = static_cast<STCyclotronRun*>(G4RunManager::GetRunManager()->GetNonConstCurrentRun());
74 fTrack = aStep->GetTrack(); 74 fTrack = aStep->GetTrack();
75 75
76 auto analysisManager = G4AnalysisManager::In 76 auto analysisManager = G4AnalysisManager::Instance();
77 77
78 //------------------------------------------ 78 //----------------------------------------------
79 // Volume info 79 // Volume info
80 //------------------------------------------ 80 //----------------------------------------------
81 G4double targetHalfDiameter= (fDet->GetTarge 81 G4double targetHalfDiameter= (fDet->GetTargetDiameter())/2.;
82 82
83 //------------------------------------------ 83 //----------------------------------------------
84 // Step information 84 // Step information
85 //------------------------------------------ 85 //----------------------------------------------
86 86
87 G4double edep = aStep->GetTotalEnergyDeposit 87 G4double edep = aStep->GetTotalEnergyDeposit();
88 G4double energy = aStep->GetPreStepPoint()-> 88 G4double energy = aStep->GetPreStepPoint()->GetKineticEnergy();
89 G4ThreeVector momentumDirection = aStep->Get 89 G4ThreeVector momentumDirection = aStep->GetPreStepPoint()->GetMomentumDirection();
90 G4ThreeVector vectorPosition = aStep->GetPre 90 G4ThreeVector vectorPosition = aStep->GetPreStepPoint()->GetPosition();
91 91
92 //------------------------------------------ 92 //----------------------------------------------
93 // Track 93 // Track
94 //------------------------------------------ 94 //----------------------------------------------
95 95
96 G4ParticleDefinition* thePartDef = fTrack->G 96 G4ParticleDefinition* thePartDef = fTrack->GetDefinition();
97 G4String partType= fTrack->GetDefinition()-> 97 G4String partType= fTrack->GetDefinition()->GetParticleType();
98 G4String name = fTrack->GetDefinition()->Get 98 G4String name = fTrack->GetDefinition()->GetParticleName();
99 G4double timeLife = fTrack->GetDefinition()- 99 G4double timeLife = fTrack->GetDefinition()->GetPDGLifeTime(); //<---
100 const G4VProcess* process = fTrack->GetCreat 100 const G4VProcess* process = fTrack->GetCreatorProcess();
101 101
102 //------------------------------------------ 102 //----------------------------------------------
103 // Collect general information concerning 103 // Collect general information concerning all of the particles
104 // Collect energy deposition ; separe dec 104 // Collect energy deposition ; separe decay case to beam case
105 //------------------------------------------ 105 //----------------------------------------------
106 106
107 fRun->EnergyDepositionTarget(edep); 107 fRun->EnergyDepositionTarget(edep);
108 108
109 109
110 //------------------------------------------ 110 //----------------------------------------------
111 //Collect information about protons and deut 111 //Collect information about protons and deuterons
112 //------------------------------------------ 112 //----------------------------------------------
113 113
114 if(name == "proton" || name == "deuteron") 114 if(name == "proton" || name == "deuteron")
115 { 115 {
116 116
117 if(fTrack->GetTrackID()!=fTempTrack && ( 117 if(fTrack->GetTrackID()!=fTempTrack && (momentumDirection.getZ()>0.) &&
118 vectorPosition.getX()<targetHalfDiameter && 118 vectorPosition.getX()<targetHalfDiameter &&
119 vectorPosition.getX()>-targetHalfDiameter & 119 vectorPosition.getX()>-targetHalfDiameter &&
120 vectorPosition.getY()<targetHalfDiameter && 120 vectorPosition.getY()<targetHalfDiameter &&
121 vectorPosition.getY()>-targetHalfDiameter) 121 vectorPosition.getY()>-targetHalfDiameter)
122 { 122 {
123 analysisManager->FillH2(0,vectorPosition.g 123 analysisManager->FillH2(0,vectorPosition.getX(),vectorPosition.getY());
124 analysisManager->FillH1(0,energy); 124 analysisManager->FillH1(0,energy);
125 fRun->CountParticlesTarget(); 125 fRun->CountParticlesTarget();
126 fTempTrack = fTrack->GetTrackID(); 126 fTempTrack = fTrack->GetTrackID();
127 } 127 }
128 128
129 if(fTempTrack1 == 0) 129 if(fTempTrack1 == 0)
130 { 130 {
131 fTempTrack1 = fTrack->GetTrackID(); 131 fTempTrack1 = fTrack->GetTrackID();
132 } 132 }
133 133
134 if(fTrack->GetTrackID()!=fTempTrack1 && 134 if(fTrack->GetTrackID()!=fTempTrack1 && (momentumDirection.getZ()>0.) &&
135 fTempVector.getX()<targetHalfDiameter && 135 fTempVector.getX()<targetHalfDiameter &&
136 fTempVector.getX()>-targetHalfDiameter && 136 fTempVector.getX()>-targetHalfDiameter &&
137 fTempVector.getY()<targetHalfDiameter && 137 fTempVector.getY()<targetHalfDiameter &&
138 fTempVector.getY()>-targetHalfDiameter ) 138 fTempVector.getY()>-targetHalfDiameter )
139 { 139 {
140 140
141 analysisManager->FillH2(4,fTempVector.getX 141 analysisManager->FillH2(4,fTempVector.getX(),fTempVector.getY());
142 analysisManager->FillH1(2,fTempEnergy); 142 analysisManager->FillH1(2,fTempEnergy);
143 fTempTrack1 = fTrack->GetTrackID(); 143 fTempTrack1 = fTrack->GetTrackID();
144 } 144 }
145 145
146 fTempVector = aStep->GetPostStepPoint() 146 fTempVector = aStep->GetPostStepPoint()->GetPosition(); //vectorPosition;
147 fTempEnergy = aStep->GetPostStepPoint()- 147 fTempEnergy = aStep->GetPostStepPoint()->GetKineticEnergy(); //energy;
148 148
149 analysisManager->FillH2(3,vectorPosition 149 analysisManager->FillH2(3,vectorPosition.getZ(),energy);
150 150
151 } 151 }
152 152
153 //------------------------------------------ 153 //----------------------------------------------
154 // Store ID for particles that are 154 // Store ID for particles that are
155 // not protons/electrons or deuterons 155 // not protons/electrons or deuterons
156 //------------------------------------------ 156 //----------------------------------------------
157 157
158 if((name != "proton") && (name != "e-") && ( 158 if((name != "proton") && (name != "e-") && (name != "deuteron"))
159 { 159 {
160 fRun->StoreIsotopeID(fTrack->GetTrackID( 160 fRun->StoreIsotopeID(fTrack->GetTrackID(),name);
161 } 161 }
162 162
163 //------------------------------------------ 163 //----------------------------------------------
164 // Collect of information for unstable isot 164 // Collect of information for unstable isotopes
165 // generated from an interaction with the t 165 // generated from an interaction with the target
166 //------------------------------------------ 166 //----------------------------------------------
167 167
168 if (name!="deuteron") 168 if (name!="deuteron")
169 { 169 {
170 if (( partType == "nucleus") && !(thePa 170 if (( partType == "nucleus") && !(thePartDef->GetPDGStable()) && (fTrack->GetCurrentStepNumber()==1) && timeLife!=0.)
171 { 171 {
172 //G4cout << "Saving unstable particles ... 172 //G4cout << "Saving unstable particles ..." << G4endl;
173 G4int Z=thePartDef->GetAtomicNumber(); 173 G4int Z=thePartDef->GetAtomicNumber();
174 G4int A=thePartDef->GetAtomicMass(); 174 G4int A=thePartDef->GetAtomicMass();
175 analysisManager->FillH2(2,Z,A); 175 analysisManager->FillH2(2,Z,A);
176 176
177 //---------------------------------------- 177 //----------------------------------------------
178 // isotopes count 178 // isotopes count
179 //---------------------------------------- 179 //----------------------------------------------
180 180
181 fRun->PrimaryIsotopeCountTarget(name,timeL 181 fRun->PrimaryIsotopeCountTarget(name,timeLife);
182 analysisManager->FillH1(4,fTrack->GetPosit 182 analysisManager->FillH1(4,fTrack->GetPosition().getZ());
183 //particle that created the nucleus 183 //particle that created the nucleus
184 std::map<G4int,G4String> parentID = fRun-> 184 std::map<G4int,G4String> parentID = fRun->GetIsotopeID();
185 G4String nameParent = parentID[fTrack->Get 185 G4String nameParent = parentID[fTrack->GetParentID()];
186 fRun->ParticleParent(name, process->GetPro 186 fRun->ParticleParent(name, process->GetProcessName());
187 187
188 //G4cout << name << " : " << process->GetP 188 //G4cout << name << " : " << process->GetProcessName() << " with track ID " << fTrack->GetTrackID() << " and step ID " << fTrack->GetCurrentStepNumber() << G4endl;
189 189
190 } 190 }
191 } 191 }
192 192
193 //------------------------------------------ 193 //----------------------------------------------
194 // Collect of information for stable isoto 194 // Collect of information for stable isotopes
195 // generated from an interaction with the 195 // generated from an interaction with the target
196 //------------------------------------------ 196 //----------------------------------------------
197 197
198 if (name!="deuteron") 198 if (name!="deuteron")
199 { 199 {
200 if (( partType == "nucleus") && (thePar 200 if (( partType == "nucleus") && (thePartDef->GetPDGStable()) && (process->GetProcessName() != "RadioactiveDecay") && (fTrack->GetCurrentStepNumber()==1) )
201 { 201 {
202 //---------------------------------------- 202 //----------------------------------------------
203 // isotopes count 203 // isotopes count
204 //---------------------------------------- 204 //----------------------------------------------
205 fRun->CountStableIsotopes(name); 205 fRun->CountStableIsotopes(name);
206 } 206 }
207 } 207 }
208 208
209 209
210 //------------------------------------------ 210 //----------------------------------------------
211 // Collect unstable isotopes from decay 211 // Collect unstable isotopes from decay
212 //------------------------------------------ 212 //----------------------------------------------
213 213
214 if (( partType == "nucleus") && !(thePartDe 214 if (( partType == "nucleus") && !(thePartDef->GetPDGStable()) && (process->GetProcessName() == "RadioactiveDecay") && (fTrack->GetCurrentStepNumber()==1) && timeLife!=0)
215 { 215 {
216 std::map<G4int,G4String>::iterator itbis 216 std::map<G4int,G4String>::iterator itbis;
217 std::map<G4int,G4String> parentID = fRun 217 std::map<G4int,G4String> parentID = fRun->GetIsotopeID();
218 G4String nameParent = parentID[fTrack->G 218 G4String nameParent = parentID[fTrack->GetParentID()];
219 fRun->DecayIsotopeCountTarget(name,nameP 219 fRun->DecayIsotopeCountTarget(name,nameParent,timeLife);
220 } 220 }
221 221
222 //------------------------------------------ 222 //----------------------------------------------
223 // Collect any other particles emitted 223 // Collect any other particles emitted
224 //------------------------------------------ 224 //----------------------------------------------
225 225
226 if((partType!="nucleus")&&(name!="proton")&& 226 if((partType!="nucleus")&&(name!="proton")&&(name!="deuteron"))
227 { 227 {
228 fRun->ParticleCountTarget(name); 228 fRun->ParticleCountTarget(name);
229 //Condition so the particle will be coun 229 //Condition so the particle will be counted for only one step
230 if((fTrack->GetCurrentStepNumber()==1)) 230 if((fTrack->GetCurrentStepNumber()==1))
231 { 231 {
232 if(process->GetProcessName() != "Radioacti 232 if(process->GetProcessName() != "RadioactiveDecay")
233 { 233 {
234 if(name=="e+"){ 234 if(name=="e+"){
235 analysisManager->FillH1(5,energy); 235 analysisManager->FillH1(5,energy);
236 } 236 }
237 if(name=="e-"){ 237 if(name=="e-"){
238 analysisManager->FillH1(6,energy); 238 analysisManager->FillH1(6,energy);
239 } 239 }
240 if(name=="gamma"){ 240 if(name=="gamma"){
241 analysisManager->FillH1(7,energy); 241 analysisManager->FillH1(7,energy);
242 } 242 }
243 if(name=="neutron"){ 243 if(name=="neutron"){
244 analysisManager->FillH1(8,energy); 244 analysisManager->FillH1(8,energy);
245 } 245 }
246 } 246 }
247 247
248 if(process->GetProcessName() == "Radioacti 248 if(process->GetProcessName() == "RadioactiveDecay")
249 { 249 {
250 if(name=="e+"){ 250 if(name=="e+"){
251 analysisManager->FillH1(9,energy); 251 analysisManager->FillH1(9,energy);
252 } 252 }
253 if(name=="e-"){ 253 if(name=="e-"){
254 analysisManager->FillH1(10,energy); 254 analysisManager->FillH1(10,energy);
255 } 255 }
256 if(name=="gamma"){ 256 if(name=="gamma"){
257 analysisManager->FillH1(11,energy); 257 analysisManager->FillH1(11,energy);
258 } 258 }
259 if(name=="neutron"){ 259 if(name=="neutron"){
260 analysisManager->FillH1(12,energy); 260 analysisManager->FillH1(12,energy);
261 } 261 }
262 if(name=="nu_e"){ 262 if(name=="nu_e"){
263 analysisManager->FillH1(13,energy); 263 analysisManager->FillH1(13,energy);
264 } 264 }
265 if(name=="anti_nu_e"){ 265 if(name=="anti_nu_e"){
266 analysisManager->FillH1(14,energy); 266 analysisManager->FillH1(14,energy);
267 } 267 }
268 } 268 }
269 } 269 }
270 } 270 }
271 271
272 272
273 fRun->SetTargetVolume(fDet->GetTargetVolume( 273 fRun->SetTargetVolume(fDet->GetTargetVolume());
274 fRun->SetTargetThickness(fDet->GetTargetThic 274 fRun->SetTargetThickness(fDet->GetTargetThickness());
275 fRun->SetTargetDiameter(fDet->GetTargetDiame 275 fRun->SetTargetDiameter(fDet->GetTargetDiameter());
276 276
277 return true; 277 return true;
278 278
279 } 279 }
280 280
281 281
282 282