| 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@gmai 26 // Author: F. Poignant, floriane.poignant@gmail.com
27 // 27 //
28 // file STCyclotronRun.cc 28 // file STCyclotronRun.cc
29 29
30 #include "STCyclotronRun.hh" 30 #include "STCyclotronRun.hh"
31 #include "G4RunManager.hh" 31 #include "G4RunManager.hh"
32 #include "G4Event.hh" 32 #include "G4Event.hh"
33 33
34 #include "G4SDManager.hh" 34 #include "G4SDManager.hh"
35 #include "G4HCofThisEvent.hh" 35 #include "G4HCofThisEvent.hh"
36 #include "G4THitsMap.hh" 36 #include "G4THitsMap.hh"
37 #include "G4SystemOfUnits.hh" 37 #include "G4SystemOfUnits.hh"
38 38
39 STCyclotronRun::STCyclotronRun() 39 STCyclotronRun::STCyclotronRun()
40 : G4Run(),fTotalEnergyDepositTarget(0.),fTot 40 : G4Run(),fTotalEnergyDepositTarget(0.),fTotalEnergyDepositFoil(0.),fParticleTarget(0),fTargetThickness(0.),fTargetDiameter(0.),fFoilThickness(0.),fTargetVolume(0.),fFoilVolume(0.),fPrimariesPerEvent(0),fTimePerEvent(0),fBeamName(""),fBeamCurrent(0.),fBeamEnergy(0.)
41 { } 41 { }
42 42
43 STCyclotronRun::~STCyclotronRun() 43 STCyclotronRun::~STCyclotronRun()
44 { } 44 { }
45 45
46 void STCyclotronRun::Merge(const G4Run* aRun) 46 void STCyclotronRun::Merge(const G4Run* aRun)
47 { 47 {
48 48
49 const STCyclotronRun* localRun = static_cast 49 const STCyclotronRun* localRun = static_cast<const STCyclotronRun*>(aRun);
50 50
51 //Merging cumulable variables 51 //Merging cumulable variables
52 fTotalEnergyDepositTarget += localRun->fTota 52 fTotalEnergyDepositTarget += localRun->fTotalEnergyDepositTarget;
53 fTotalEnergyDepositFoil += localRun->fTotalE 53 fTotalEnergyDepositFoil += localRun->fTotalEnergyDepositFoil;
54 fParticleTarget += localRun->fParticleTarget 54 fParticleTarget += localRun->fParticleTarget;
55 55
56 //Constant over the different runs 56 //Constant over the different runs
57 if(localRun->fTargetVolume!=0)fTargetVolume 57 if(localRun->fTargetVolume!=0)fTargetVolume = localRun->fTargetVolume;
58 if(localRun->fFoilVolume!=0)fFoilVolume = lo 58 if(localRun->fFoilVolume!=0)fFoilVolume = localRun->fFoilVolume;
59 if(localRun->fPrimariesPerEvent!=0)fPrimarie 59 if(localRun->fPrimariesPerEvent!=0)fPrimariesPerEvent = localRun->fPrimariesPerEvent;
60 if(localRun->fTimePerEvent!=0)fTimePerEvent 60 if(localRun->fTimePerEvent!=0)fTimePerEvent = localRun->fTimePerEvent;
61 if(localRun->fTargetThickness!=0)fTargetThic 61 if(localRun->fTargetThickness!=0)fTargetThickness = localRun->fTargetThickness;
62 if(localRun->fTargetDiameter!=0)fTargetDiame 62 if(localRun->fTargetDiameter!=0)fTargetDiameter = localRun->fTargetDiameter;
63 if(localRun->fFoilThickness!=0)fFoilThicknes 63 if(localRun->fFoilThickness!=0)fFoilThickness = localRun->fFoilThickness;
64 fBeamName = localRun->fBeamName; 64 fBeamName = localRun->fBeamName;
65 if(localRun->fBeamCurrent!=0.)fBeamCurrent = 65 if(localRun->fBeamCurrent!=0.)fBeamCurrent = localRun->fBeamCurrent;
66 if(localRun->fBeamEnergy!=0.)fBeamEnergy = l 66 if(localRun->fBeamEnergy!=0.)fBeamEnergy = localRun->fBeamEnergy;
67 67
68 68
69 69
70 //<<<----toMerge 70 //<<<----toMerge
71 std::map<G4String,G4int>::iterator itSI; 71 std::map<G4String,G4int>::iterator itSI;
72 std::map<G4String,G4double>::iterator itSD; 72 std::map<G4String,G4double>::iterator itSD;
73 std::map<G4String,G4String>::iterator itSS; 73 std::map<G4String,G4String>::iterator itSS;
74 std::map<G4int,G4String>::iterator itIS; 74 std::map<G4int,G4String>::iterator itIS;
75 75
76 //----Merging results for primary isotopes 76 //----Merging results for primary isotopes
77 std::map<G4String,G4int> locPrimaryIsotopeCo 77 std::map<G4String,G4int> locPrimaryIsotopeCountTarget = localRun->fPrimaryIsotopeCountTarget;
78 for (itSI = locPrimaryIsotopeCountTarget.beg 78 for (itSI = locPrimaryIsotopeCountTarget.begin(); itSI != locPrimaryIsotopeCountTarget.end(); itSI++)
79 { 79 {
80 G4String name = itSI->first; 80 G4String name = itSI->first;
81 G4int count = itSI->second; 81 G4int count = itSI->second;
82 fPrimaryIsotopeCountTarget[name] += coun 82 fPrimaryIsotopeCountTarget[name] += count;
83 } 83 }
84 84
85 std::map<G4String,G4double> locPrimaryIsotop 85 std::map<G4String,G4double> locPrimaryIsotopeTimeTarget = localRun->fPrimaryIsotopeTimeTarget;
86 for (itSD = locPrimaryIsotopeTimeTarget.begi 86 for (itSD = locPrimaryIsotopeTimeTarget.begin(); itSD != locPrimaryIsotopeTimeTarget.end(); itSD++)
87 { 87 {
88 G4String name = itSD->first; 88 G4String name = itSD->first;
89 G4double time = itSD->second; 89 G4double time = itSD->second;
90 fPrimaryIsotopeTimeTarget[name] = time; 90 fPrimaryIsotopeTimeTarget[name] = time;
91 } 91 }
92 92
93 //----Merging results for decay isotopes 93 //----Merging results for decay isotopes
94 // std::map<G4int,G4String> fIsotopeIDTar 94 // std::map<G4int,G4String> fIsotopeIDTarget;
95 std::map<G4String,G4String> locDecayIsotopeC 95 std::map<G4String,G4String> locDecayIsotopeCountTarget = localRun->fDecayIsotopeCountTarget;
96 for (itSS = locDecayIsotopeCountTarget.begin 96 for (itSS = locDecayIsotopeCountTarget.begin(); itSS != locDecayIsotopeCountTarget.end(); itSS++)
97 { 97 {
98 G4String nameDaughter = itSS->first; 98 G4String nameDaughter = itSS->first;
99 G4String mum = itSS->second; 99 G4String mum = itSS->second;
100 fDecayIsotopeCountTarget[nameDaughter] = 100 fDecayIsotopeCountTarget[nameDaughter] = mum;
101 } 101 }
102 std::map<G4String,G4double> locDecayIsotopeT 102 std::map<G4String,G4double> locDecayIsotopeTimeTarget = localRun->fDecayIsotopeTimeTarget;
103 for (itSD = locDecayIsotopeTimeTarget.begin 103 for (itSD = locDecayIsotopeTimeTarget.begin(); itSD != locDecayIsotopeTimeTarget.end(); itSD++)
104 { 104 {
105 G4String nameDaughter = itSD->first; 105 G4String nameDaughter = itSD->first;
106 G4double time = itSD->second; 106 G4double time = itSD->second;
107 fDecayIsotopeTimeTarget[nameDaughter] = 107 fDecayIsotopeTimeTarget[nameDaughter] = time;
108 } 108 }
109 std::map<G4String,G4String> locParticlePare 109 std::map<G4String,G4String> locParticleParent = localRun->fParticleParent;
110 for (itSS = locParticleParent.begin(); itSS 110 for (itSS = locParticleParent.begin(); itSS != locParticleParent.end(); itSS++)
111 { 111 {
112 G4String nameDaughter = itSS->first; 112 G4String nameDaughter = itSS->first;
113 G4String parent = itSS->second; 113 G4String parent = itSS->second;
114 fParticleParent[nameDaughter] = parent; 114 fParticleParent[nameDaughter] = parent;
115 } 115 }
116 std::map<G4int,G4String> locIsotopeIDTarget 116 std::map<G4int,G4String> locIsotopeIDTarget = localRun->fIsotopeIDTarget;
117 for (itIS = locIsotopeIDTarget.begin(); itI 117 for (itIS = locIsotopeIDTarget.begin(); itIS != locIsotopeIDTarget.end(); itIS++)
118 { 118 {
119 G4int ID = itIS->first; 119 G4int ID = itIS->first;
120 G4String name = itIS->second; 120 G4String name = itIS->second;
121 fIsotopeIDTarget[ID] = name; 121 fIsotopeIDTarget[ID] = name;
122 } 122 }
123 123
124 124
125 //----Merging results for stable isotopes 125 //----Merging results for stable isotopes
126 std::map<G4String,G4int> locStableIsotopeCou 126 std::map<G4String,G4int> locStableIsotopeCountTarget = localRun->fStableIsotopeCountTarget;
127 for (itSI = locStableIsotopeCountTarget.begi 127 for (itSI = locStableIsotopeCountTarget.begin(); itSI != locStableIsotopeCountTarget.end(); itSI++)
128 { 128 {
129 G4String name = itSI->first; 129 G4String name = itSI->first;
130 G4int count = itSI->second; 130 G4int count = itSI->second;
131 fStableIsotopeCountTarget[name] += count 131 fStableIsotopeCountTarget[name] += count;
132 } 132 }
133 133
134 //----Merging results for particles 134 //----Merging results for particles
135 std::map<G4String,G4int> locParticleCountTar 135 std::map<G4String,G4int> locParticleCountTarget = localRun->fParticleCountTarget;
136 for (itSI = locParticleCountTarget.begin(); 136 for (itSI = locParticleCountTarget.begin(); itSI != locParticleCountTarget.end(); itSI++)
137 { 137 {
138 G4String name = itSI->first; 138 G4String name = itSI->first;
139 G4int count = itSI->second; 139 G4int count = itSI->second;
140 fParticleCountTarget[name] += count; 140 fParticleCountTarget[name] += count;
141 } 141 }
142 142
143 143
144 G4Run::Merge(aRun); 144 G4Run::Merge(aRun);
145 145
146 } 146 }
147 147
148 void STCyclotronRun::EndOfRun(G4double irradia 148 void STCyclotronRun::EndOfRun(G4double irradiationTime)
149 { 149 {
150 150
151 G4int nbEvents = GetNumberOfEvent(); 151 G4int nbEvents = GetNumberOfEvent();
152 152
153 if (nbEvents == 0) return; 153 if (nbEvents == 0) return;
154 154
155 //------------------------------------------ 155 //------------------------------------------------------
156 // Opening the ASCII file 156 // Opening the ASCII file
157 //------------------------------------------ 157 //------------------------------------------------------
158 fOutPut.open("Output_General.txt",std::ofstr 158 fOutPut.open("Output_General.txt",std::ofstream::out);
159 fOutPut1.open("Output_ParentIsotopes.txt",st 159 fOutPut1.open("Output_ParentIsotopes.txt",std::ofstream::out);
160 fOutPut2.open("Output_DaughterIsotopes.txt", 160 fOutPut2.open("Output_DaughterIsotopes.txt",std::ofstream::out);
161 fOutPut3.open("Output_OtherParticles.txt",st 161 fOutPut3.open("Output_OtherParticles.txt",std::ofstream::out);
162 fOutPut4.open("Output_StableIsotopes.txt",st 162 fOutPut4.open("Output_StableIsotopes.txt",std::ofstream::out);
163 163
164 //------------------------------------------ 164 //------------------------------------------------------
165 // Calculates the equivalent time for a giv 165 // Calculates the equivalent time for a given run
166 //------------------------------------------ 166 //------------------------------------------------------
167 G4double timePerEvent = fTimePerEvent; //in 167 G4double timePerEvent = fTimePerEvent; //in seconds
168 G4double timeForARun = nbEvents*timePerEvent 168 G4double timeForARun = nbEvents*timePerEvent; //in seconds
169 G4double minDecay = 0.0001; //in seconds 169 G4double minDecay = 0.0001; //in seconds
170 G4double maxDecay = 1000000.; //in seconds 170 G4double maxDecay = 1000000.; //in seconds
171 171
172 //------------------------------------------ 172 //------------------------------------------------------
173 // Rescale the value of the beam current 173 // Rescale the value of the beam current to account for
174 // the loss of primary particles due to t 174 // the loss of primary particles due to the foil.
175 //------------------------------------------ 175 //------------------------------------------------------
176 176
177 G4int totalPrimaries = fPrimariesPerEvent*nb 177 G4int totalPrimaries = fPrimariesPerEvent*nbEvents;
178 G4double currentFactor; 178 G4double currentFactor;
179 if(fParticleTarget>0.) currentFactor =(fPart 179 if(fParticleTarget>0.) currentFactor =(fParticleTarget*1.)/(totalPrimaries*1.);
180 else currentFactor = 0.; 180 else currentFactor = 0.;
181 181
182 182
183 fOutPut << "//------------------------------ 183 fOutPut << "//-----------------------------------//" << G4endl;
184 fOutPut << "// Parameters of the simulati 184 fOutPut << "// Parameters of the simulation: //" << G4endl;
185 fOutPut << "//------------------------------ 185 fOutPut << "//-----------------------------------//" << G4endl;
186 fOutPut << "Beam parameters: " << G4endl; 186 fOutPut << "Beam parameters: " << G4endl;
187 fOutPut << fBeamName << " - Name of beam pri 187 fOutPut << fBeamName << " - Name of beam primary particles." << G4endl;
188 fOutPut << fBeamEnergy << " - Energy of beam 188 fOutPut << fBeamEnergy << " - Energy of beam primary particles (MeV)." << G4endl;
189 fOutPut << fBeamCurrent << " - Beam current 189 fOutPut << fBeamCurrent << " - Beam current (Ampere)." << G4endl;
190 fOutPut << irradiationTime << " - Irradiatio 190 fOutPut << irradiationTime << " - Irradiation time in hour(s)." << G4endl;
191 fOutPut << currentFactor << " - Current fact 191 fOutPut << currentFactor << " - Current factor." << G4endl;
192 fOutPut << "//------------------------------ 192 fOutPut << "//-----------------------------------//" << G4endl;
193 fOutPut << "Simulation parameters: " << G4en 193 fOutPut << "Simulation parameters: " << G4endl;
194 fOutPut << timePerEvent << " - Equivalent ti 194 fOutPut << timePerEvent << " - Equivalent time per event (s)." << G4endl;
195 fOutPut << nbEvents << " - Number of events" 195 fOutPut << nbEvents << " - Number of events" << G4endl;
196 fOutPut << fPrimariesPerEvent << " - Primari 196 fOutPut << fPrimariesPerEvent << " - Primaries per event" << G4endl;
197 fOutPut << fPrimariesPerEvent*nbEvents << " 197 fOutPut << fPrimariesPerEvent*nbEvents << " - Total number of particles sent." << G4endl;
198 fOutPut << "//------------------------------ 198 fOutPut << "//-----------------------------------//" << G4endl;
199 fOutPut << "Geometry parameters: " << G4endl 199 fOutPut << "Geometry parameters: " << G4endl;
200 fOutPut << fTargetThickness << " - target th 200 fOutPut << fTargetThickness << " - target thickness (mm)." << G4endl;
201 fOutPut << fTargetDiameter << " - target dia 201 fOutPut << fTargetDiameter << " - target diameter (mm)." << G4endl;
202 fOutPut << fFoilThickness << " - foil thickn 202 fOutPut << fFoilThickness << " - foil thickness (mm)." << G4endl;
203 //Add particle type, particle energy, beam d 203 //Add particle type, particle energy, beam diameter, beam current
204 204
205 //target material??? 205 //target material???
206 206
207 //////////////////////////////////////////// 207 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
208 //////////Calculation of the number of isoto 208 //////////Calculation of the number of isotopes at the end of the irradiation and the activity generated/////////
209 //////////////////////////////////////////// 209 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
210 210
211 211
212 //Maps to fill 212 //Maps to fill
213 std::map<G4String,G4double> fPrimaryIsotopeE 213 std::map<G4String,G4double> fPrimaryIsotopeEOBTarget;
214 std::map<G4String,G4double> fPrimaryActivity 214 std::map<G4String,G4double> fPrimaryActivityTarget;
215 std::map<G4String,G4double> fDecayIsotopeEOB 215 std::map<G4String,G4double> fDecayIsotopeEOBTarget;
216 std::map<G4String,G4double> fDecayActivityTa 216 std::map<G4String,G4double> fDecayActivityTarget;
217 217
218 218
219 //------------------------------------------ 219 //----------------------------------------------
220 // CASE 1 : Parent isotopes 220 // CASE 1 : Parent isotopes
221 //------------------------------------------ 221 //----------------------------------------------
222 222
223 223
224 std::map<G4String,G4int>::iterator it; 224 std::map<G4String,G4int>::iterator it;
225 G4double primaryActivityTotal = 0.; 225 G4double primaryActivityTotal = 0.;
226 G4double decayActivityTotal=0.; 226 G4double decayActivityTotal=0.;
227 227
228 fOutPut1 << "//----------------------------- 228 fOutPut1 << "//-----------------------------------//\n"
229 << "// Data for parent isotopes / 229 << "// Data for parent isotopes //\n"
230 << "//-----------------------------------/ 230 << "//-----------------------------------//\n" << G4endl;
231 231
232 for (it = fPrimaryIsotopeCountTarget.begin() 232 for (it = fPrimaryIsotopeCountTarget.begin(); it != fPrimaryIsotopeCountTarget.end(); it++)
233 { 233 {
234 234
235 235
236 G4String name = it->first; 236 G4String name = it->first;
237 G4double count = (it->second)*currentF 237 G4double count = (it->second)*currentFactor;
238 G4double halfLifeTime = fPrimaryIsotopeT 238 G4double halfLifeTime = fPrimaryIsotopeTimeTarget[name]*10E-10/3600.*std::log(2.);
239 G4String process = fParticleParent[name] 239 G4String process = fParticleParent[name];
240 240
241 //Only store isotopes with a life time b 241 //Only store isotopes with a life time between minDecay and maxDecay.
242 G4bool store; 242 G4bool store;
243 if(halfLifeTime > minDecay && halfLifeTi 243 if(halfLifeTime > minDecay && halfLifeTime < maxDecay) store = true;
244 else store = false; 244 else store = false;
245 245
246 246
247 //Calculation of the yield (s-1) 247 //Calculation of the yield (s-1)
248 G4double decayConstant = 1/(fPrimaryIsot 248 G4double decayConstant = 1/(fPrimaryIsotopeTimeTarget[name]*10E-10);
249 249
250 250
251 //-------------------------------------- 251 //----------------------------------------------
252 // Number of particles per second 252 // Number of particles per second
253 //-------------------------------------- 253 //----------------------------------------------
254 254
255 G4double particlesPerSecond = fPrimaryIs 255 G4double particlesPerSecond = fPrimaryIsotopeCountTarget[name]*currentFactor/timeForARun;
256 256
257 //-------------------------------------- 257 //----------------------------------------------
258 // Calculation yield EOB 258 // Calculation yield EOB
259 //-------------------------------------- 259 //----------------------------------------------
260 260
261 fPrimaryIsotopeEOBTarget[name] = particl 261 fPrimaryIsotopeEOBTarget[name] = particlesPerSecond/decayConstant * (1. - std::exp(-irradiationTime*3600*decayConstant));
262 262
263 //-------------------------------------- 263 //----------------------------------------------
264 // Calculation of the activity 264 // Calculation of the activity
265 // conversion factor Bq to mCi 265 // conversion factor Bq to mCi
266 //-------------------------------------- 266 //----------------------------------------------
267 267
268 G4double conv = 2.7E-8; 268 G4double conv = 2.7E-8;
269 fPrimaryActivityTarget[name]= fPrimaryIs 269 fPrimaryActivityTarget[name]= fPrimaryIsotopeEOBTarget[name]*decayConstant*conv;
270 270
271 if(store) 271 if(store)
272 { 272 {
273 273
274 //---------------------------------------- 274 //----------------------------------------------
275 // Incrementation for total primary act 275 // Incrementation for total primary activity
276 //---------------------------------------- 276 //----------------------------------------------
277 277
278 primaryActivityTotal = primaryActivityTota 278 primaryActivityTotal = primaryActivityTotal + fPrimaryActivityTarget[name];
279 } 279 }
280 280
281 281
282 //---------------------------// 282 //---------------------------//
283 // Printing out results // 283 // Printing out results //
284 //---------------------------// 284 //---------------------------//
285 285
286 if(store) 286 if(store)
287 { 287 {
288 fOutPut1 << name << " - name of parent iso 288 fOutPut1 << name << " - name of parent isotope." << G4endl;
289 fOutPut1 << count/currentFactor << " - num 289 fOutPut1 << count/currentFactor << " - number of isotopes created during the simulation." << G4endl;
290 fOutPut1 << decayConstant << " - decay con 290 fOutPut1 << decayConstant << " - decay constant in s-1." << G4endl;
291 fOutPut1 << halfLifeTime << " - half life 291 fOutPut1 << halfLifeTime << " - half life time in hour(s)." << G4endl;
292 fOutPut1 << process << " - creation proces 292 fOutPut1 << process << " - creation process." << G4endl;
293 fOutPut1 << particlesPerSecond << " - isot 293 fOutPut1 << particlesPerSecond << " - isotope per sec." << G4endl;
294 fOutPut1 << fPrimaryIsotopeEOBTarget[name] 294 fOutPut1 << fPrimaryIsotopeEOBTarget[name] << " - yield EOB." << G4endl;
295 fOutPut1 << fPrimaryActivityTarget[name] < 295 fOutPut1 << fPrimaryActivityTarget[name] << " - activity (mCi) at the EOB." << G4endl;
296 fOutPut1 << "------------------------" << 296 fOutPut1 << "------------------------" << G4endl;
297 } 297 }
298 298
299 } 299 }
300 300
301 301
302 //------------------------------------------ 302 //----------------------------------------------
303 // CASE 2 : isotopes from primary isotopes 303 // CASE 2 : isotopes from primary isotopes decay
304 //------------------------------------------ 304 //----------------------------------------------
305 305
306 fOutPut2 << "//----------------------------- 306 fOutPut2 << "//-----------------------------------//\n"
307 << "// Data for daughter isotopes / 307 << "// Data for daughter isotopes //\n"
308 << "//-----------------------------------/ 308 << "//-----------------------------------//\n" << G4endl;
309 309
310 std::map<G4String,G4String>::iterator it1; 310 std::map<G4String,G4String>::iterator it1;
311 311
312 for (it1 = fDecayIsotopeCountTarget.begin(); 312 for (it1 = fDecayIsotopeCountTarget.begin(); it1 != fDecayIsotopeCountTarget.end(); it1++)
313 { 313 {
314 314
315 315
316 G4String nameDaughter = it1->first; 316 G4String nameDaughter = it1->first;
317 G4String nameMum = it1->second; 317 G4String nameMum = it1->second;
318 318
319 G4double halfLifeTimeMum = fDecayIsotope 319 G4double halfLifeTimeMum = fDecayIsotopeTimeTarget[nameDaughter]*10E10/3600;
320 G4double halfLifeTimeDaughter = fPrimary 320 G4double halfLifeTimeDaughter = fPrimaryIsotopeTimeTarget[nameMum]*10E10/3600;
321 321
322 G4bool store; 322 G4bool store;
323 if(halfLifeTimeMum > minDecay && halfLif 323 if(halfLifeTimeMum > minDecay && halfLifeTimeMum < maxDecay &&
324 halfLifeTimeDaughter > minDecay && halfLife 324 halfLifeTimeDaughter > minDecay && halfLifeTimeDaughter < maxDecay){store=true;}
325 else{store=false;} 325 else{store=false;}
326 326
327 327
328 328
329 //-------------------------------------- 329 //----------------------------------------------
330 // Calculation of the yield 330 // Calculation of the yield
331 // fParticleTime[name] is the time 331 // fParticleTime[name] is the time
332 // life of the particle, divided b 332 // life of the particle, divided by ln(2), in nS
333 //-------------------------------------- 333 //----------------------------------------------
334 334
335 G4double decayConstantMum = fPrimaryIsot 335 G4double decayConstantMum = fPrimaryIsotopeTimeTarget[nameMum]
336 ? 1/(fPrimaryI 336 ? 1/(fPrimaryIsotopeTimeTarget[nameMum]*10.E-10) : 0.;
337 G4double decayConstantDaughter = fDecayI 337 G4double decayConstantDaughter = fDecayIsotopeTimeTarget[nameDaughter]
338 ? 1/(fDec 338 ? 1/(fDecayIsotopeTimeTarget[nameDaughter]*10.E-10) : 0.;
339 339
340 340
341 //-------------------------------------- 341 //----------------------------------------------
342 // Number of particles per secon 342 // Number of particles per second
343 //-------------------------------------- 343 //----------------------------------------------
344 344
345 G4double particlesPerSecond = fPrimaryIs 345 G4double particlesPerSecond = fPrimaryIsotopeCountTarget[nameMum]*currentFactor/timeForARun;
346 346
347 //-------------------------------------- 347 //----------------------------------------------
348 // Number of particles at the EOB 348 // Number of particles at the EOB
349 //-------------------------------------- 349 //----------------------------------------------
350 350
351 fDecayIsotopeEOBTarget[nameDaughter] = p 351 fDecayIsotopeEOBTarget[nameDaughter] = particlesPerSecond*((1 - std::exp(-irradiationTime*3600*decayConstantDaughter))/decayConstantDaughter + (std::exp(-irradiationTime*3600*decayConstantDaughter) - std::exp(-irradiationTime*3600*decayConstantMum))/(decayConstantDaughter-decayConstantMum));
352 352
353 353
354 //-------------------------------------- 354 //----------------------------------------------
355 // Calculation of activity 355 // Calculation of activity
356 // conversion factor Bq to mCu 356 // conversion factor Bq to mCu
357 //-------------------------------------- 357 //----------------------------------------------
358 358
359 G4double conv = 2.7E-8; 359 G4double conv = 2.7E-8;
360 fDecayActivityTarget[nameDaughter]= fDec 360 fDecayActivityTarget[nameDaughter]= fDecayIsotopeEOBTarget[nameDaughter]*decayConstantDaughter*conv;
361 361
362 if(store) 362 if(store)
363 { 363 {
364 decayActivityTotal = decayActivityTotal + 364 decayActivityTotal = decayActivityTotal + fDecayActivityTarget[nameDaughter];
365 } 365 }
366 366
367 if(store) 367 if(store)
368 { 368 {
369 fOutPut2 << nameDaughter << " - name of da 369 fOutPut2 << nameDaughter << " - name of daughter isotope." << G4endl;
370 fOutPut2 << nameMum << " - name of parent 370 fOutPut2 << nameMum << " - name of parent isotope." << G4endl;
371 fOutPut2 << decayConstantDaughter << " - d 371 fOutPut2 << decayConstantDaughter << " - decay constant of daughter in s-1." << G4endl;
372 fOutPut2 << decayConstantMum << " - decay 372 fOutPut2 << decayConstantMum << " - decay constant of mum in s-1." << G4endl;
373 fOutPut2 << halfLifeTimeDaughter << " - ha 373 fOutPut2 << halfLifeTimeDaughter << " - half life time of daughter in hour(s)." << G4endl;
374 fOutPut2 << halfLifeTimeMum << " - half li 374 fOutPut2 << halfLifeTimeMum << " - half life time of mum in hour(s)." << G4endl;
375 fOutPut2 << particlesPerSecond << " - iso 375 fOutPut2 << particlesPerSecond << " - isotope per sec." << G4endl;
376 fOutPut2 << fDecayIsotopeEOBTarget[nameDau 376 fOutPut2 << fDecayIsotopeEOBTarget[nameDaughter] << " - yield at the EOB." << G4endl;
377 fOutPut2 << fDecayActivityTarget[nameDaugh 377 fOutPut2 << fDecayActivityTarget[nameDaughter] << " - activity (mCi) at the EOB." << G4endl;
378 fOutPut2 << "------------------------" << 378 fOutPut2 << "------------------------" << G4endl;
379 } 379 }
380 380
381 } 381 }
382 382
383 //------------------------------------------ 383 //----------------------------------------------
384 // Particles created, other than nuclei 384 // Particles created, other than nuclei
385 //------------------------------------------ 385 //----------------------------------------------
386 386
387 fOutPut3 << "//----------------------------- 387 fOutPut3 << "//-----------------------------------//\n"
388 << "// Data for other particles / 388 << "// Data for other particles //\n"
389 << "//-----------------------------------/ 389 << "//-----------------------------------//" << G4endl;
390 390
391 std::map<G4String, G4int>::iterator it3; 391 std::map<G4String, G4int>::iterator it3;
392 for(it3=fParticleCountTarget.begin(); it3!= 392 for(it3=fParticleCountTarget.begin(); it3!= fParticleCountTarget.end(); it3++)
393 { 393 {
394 G4String name = it3->first; 394 G4String name = it3->first;
395 G4double number = it3->second; 395 G4double number = it3->second;
396 fOutPut3 << name << " - name of the part 396 fOutPut3 << name << " - name of the particle" << G4endl;
397 fOutPut3 << number << " - number of part 397 fOutPut3 << number << " - number of particles" << G4endl;
398 fOutPut3 << "------------------------" < 398 fOutPut3 << "------------------------" << G4endl;
399 } 399 }
400 400
401 401
402 402
403 403
404 fOutPut4 << "//----------------------------- 404 fOutPut4 << "//-----------------------------------//\n"
405 << "// Data for stable isotopes / 405 << "// Data for stable isotopes //\n"
406 << "//-----------------------------------/ 406 << "//-----------------------------------//\n" << G4endl;
407 407
408 std::map<G4String, G4int>::iterator it6; 408 std::map<G4String, G4int>::iterator it6;
409 for(it6=fStableIsotopeCountTarget.begin();it 409 for(it6=fStableIsotopeCountTarget.begin();it6!=fStableIsotopeCountTarget.end();it6++)
410 { 410 {
411 G4String isotope = it6 ->first; 411 G4String isotope = it6 ->first;
412 G4int number = it6 -> second; 412 G4int number = it6 -> second;
413 fOutPut4 << isotope << " - name of the i 413 fOutPut4 << isotope << " - name of the isotope" << G4endl;
414 fOutPut4 << number << " - number of isot 414 fOutPut4 << number << " - number of isotopes" << G4endl;
415 fOutPut4 << "------------------------" < 415 fOutPut4 << "------------------------" << G4endl;
416 } 416 }
417 417
418 418
419 419
420 //Clear the maps 420 //Clear the maps
421 fPrimaryIsotopeEOBTarget.clear(); 421 fPrimaryIsotopeEOBTarget.clear();
422 fPrimaryActivityTarget.clear(); 422 fPrimaryActivityTarget.clear();
423 fDecayIsotopeEOBTarget.clear(); 423 fDecayIsotopeEOBTarget.clear();
424 fDecayActivityTarget.clear(); 424 fDecayActivityTarget.clear();
425 425
426 426
427 //Clear the maps 427 //Clear the maps
428 fPrimaryIsotopeCountTarget.clear(); 428 fPrimaryIsotopeCountTarget.clear();
429 fPrimaryIsotopeTimeTarget.clear(); 429 fPrimaryIsotopeTimeTarget.clear();
430 fDecayIsotopeCountTarget.clear(); 430 fDecayIsotopeCountTarget.clear();
431 fDecayIsotopeTimeTarget.clear(); 431 fDecayIsotopeTimeTarget.clear();
432 fParticleParent.clear(); 432 fParticleParent.clear();
433 fParticleCountTarget.clear(); 433 fParticleCountTarget.clear();
434 fStableIsotopeCountTarget.clear(); 434 fStableIsotopeCountTarget.clear();
435 fIsotopeIDTarget.clear(); 435 fIsotopeIDTarget.clear();
436 436
437 437
438 //----------------------------- 438 //-----------------------------
439 // Calculation of heat 439 // Calculation of heat
440 //----------------------------- 440 //-----------------------------
441 441
442 G4double totalEnergyDepositTargetEOB = fTota 442 G4double totalEnergyDepositTargetEOB = fTotalEnergyDepositTarget/timeForARun * irradiationTime * 3600.;
443 G4double totalEnergyDepositTargetPerSecond = 443 G4double totalEnergyDepositTargetPerSecond = fTotalEnergyDepositTarget/timeForARun;
444 //Heat calculation in W/mm3 444 //Heat calculation in W/mm3
445 G4double heatTarget = totalEnergyDepositTarg 445 G4double heatTarget = totalEnergyDepositTargetPerSecond/fTargetVolume * 1.60E-13;
446 G4double heatFoil = fTotalEnergyDepositFoil 446 G4double heatFoil = fTotalEnergyDepositFoil / fFoilVolume * 1.60E-13;
447 447
448 448
449 //Output data in a .txt file 449 //Output data in a .txt file
450 fOutPut << "//------------------------------ 450 fOutPut << "//-------------------------------------------------//\n"
451 << "// Heating, total activity and proc 451 << "// Heating, total activity and process data //\n"
452 << "//------------------------------------- 452 << "//-------------------------------------------------//" << G4endl;
453 453
454 fOutPut << "Total heating in the target : " 454 fOutPut << "Total heating in the target : "
455 << heatTarget << " W/mm3" << G4endl; 455 << heatTarget << " W/mm3" << G4endl;
456 fOutPut << "The total heating during the irr 456 fOutPut << "The total heating during the irradiation is " << totalEnergyDepositTargetEOB << "J/mm3" << G4endl;
457 fOutPut << "Total heating in the foil : " << 457 fOutPut << "Total heating in the foil : " << heatFoil << " W/mm3" << G4endl;
458 458
459 459
460 fOutPut.close(); 460 fOutPut.close();
461 fOutPut1.close(); 461 fOutPut1.close();
462 fOutPut2.close(); 462 fOutPut2.close();
463 fOutPut3.close(); 463 fOutPut3.close();
464 fOutPut4.close(); 464 fOutPut4.close();
465 465
466 } 466 }
467 467
468 468
469 469
470 470
471 // Accumulation functions for maps used at the 471 // Accumulation functions for maps used at the end of run action
472 472
473 void STCyclotronRun::PrimaryIsotopeCountTarget 473 void STCyclotronRun::PrimaryIsotopeCountTarget(G4String name,G4double time)
474 { 474 {
475 fPrimaryIsotopeCountTarget[name]++; 475 fPrimaryIsotopeCountTarget[name]++;
476 fPrimaryIsotopeTimeTarget[name]=time; 476 fPrimaryIsotopeTimeTarget[name]=time;
477 } 477 }
478 //------------------------------------ 478 //------------------------------------
479 void STCyclotronRun::CountStableIsotopes(G4Str 479 void STCyclotronRun::CountStableIsotopes(G4String name)
480 { 480 {
481 fStableIsotopeCountTarget[name]++; 481 fStableIsotopeCountTarget[name]++;
482 } 482 }
483 //------------------------------------ 483 //------------------------------------
484 void STCyclotronRun::DecayIsotopeCountTarget(G 484 void STCyclotronRun::DecayIsotopeCountTarget(G4String nameDaughter,G4String mum, G4double time)
485 { 485 {
486 fDecayIsotopeCountTarget[nameDaughter]=mum; 486 fDecayIsotopeCountTarget[nameDaughter]=mum;
487 fDecayIsotopeTimeTarget[nameDaughter]=time; 487 fDecayIsotopeTimeTarget[nameDaughter]=time;
488 } 488 }
489 //------------------------------------ 489 //------------------------------------
490 void STCyclotronRun::ParticleParent(G4String i 490 void STCyclotronRun::ParticleParent(G4String isotope, G4String parent)
491 { 491 {
492 fParticleParent[isotope]=parent; 492 fParticleParent[isotope]=parent;
493 } 493 }
494 // 494 //
495 //-----> Count other particles 495 //-----> Count other particles
496 //------------------------------------ 496 //------------------------------------
497 void STCyclotronRun::ParticleCountTarget(G4Str 497 void STCyclotronRun::ParticleCountTarget(G4String name)
498 { 498 {
499 fParticleCountTarget[name]++; 499 fParticleCountTarget[name]++;
500 } 500 }
501 501
502 502
503 503
504 //-------------------------------------------- 504 //-------------------------------------------------------------------------------------------------------------
505 // Accumulation functions for maps used only d 505 // Accumulation functions for maps used only during the run
506 // 506 //
507 //-----> Isotope ID to obtain the "mother isot 507 //-----> Isotope ID to obtain the "mother isotope" in SensitiveTarget()
508 //------------------------------------ 508 //------------------------------------
509 void STCyclotronRun::StoreIsotopeID(G4int ID, 509 void STCyclotronRun::StoreIsotopeID(G4int ID, G4String name)
510 { 510 {
511 fIsotopeIDTarget[ID]=name; 511 fIsotopeIDTarget[ID]=name;
512 } 512 }
513 // 513 //
514 std::map<G4int,G4String> STCyclotronRun::GetIs 514 std::map<G4int,G4String> STCyclotronRun::GetIsotopeID()
515 { 515 {
516 return fIsotopeIDTarget; 516 return fIsotopeIDTarget;
517 } 517 }
518 518
519 519
520 void STCyclotronRun::EnergyDepositionTarget(G4 520 void STCyclotronRun::EnergyDepositionTarget(G4double edep)
521 { 521 {
522 fTotalEnergyDepositTarget += edep; 522 fTotalEnergyDepositTarget += edep;
523 } 523 }
524 524
525 void STCyclotronRun::EnergyDepositionFoil(G4do 525 void STCyclotronRun::EnergyDepositionFoil(G4double edep)
526 { 526 {
527 fTotalEnergyDepositFoil += edep; 527 fTotalEnergyDepositFoil += edep;
528 } 528 }
529 529
530 void STCyclotronRun::CountParticlesTarget() 530 void STCyclotronRun::CountParticlesTarget()
531 { 531 {
532 fParticleTarget++; 532 fParticleTarget++;
533 } 533 }
534 534
535 void STCyclotronRun::SetFoilVolume(G4double fo 535 void STCyclotronRun::SetFoilVolume(G4double foilVolume)
536 { 536 {
537 fFoilVolume = foilVolume; 537 fFoilVolume = foilVolume;
538 } 538 }
539 539
540 void STCyclotronRun::SetFoilThickness(G4double 540 void STCyclotronRun::SetFoilThickness(G4double foilThickness)
541 { 541 {
542 fFoilThickness = foilThickness; 542 fFoilThickness = foilThickness;
543 } 543 }
544 544
545 void STCyclotronRun::SetTargetVolume(G4double 545 void STCyclotronRun::SetTargetVolume(G4double targetVolume)
546 { 546 {
547 fTargetVolume = targetVolume; 547 fTargetVolume = targetVolume;
548 } 548 }
549 549
550 void STCyclotronRun::SetTargetThickness(G4doub 550 void STCyclotronRun::SetTargetThickness(G4double targetThickness)
551 { 551 {
552 fTargetThickness = targetThickness; 552 fTargetThickness = targetThickness;
553 } 553 }
554 554
555 void STCyclotronRun::SetTargetDiameter(G4doubl 555 void STCyclotronRun::SetTargetDiameter(G4double targetDiameter)
556 { 556 {
557 fTargetDiameter = targetDiameter; 557 fTargetDiameter = targetDiameter;
558 } 558 }
559 559
560 void STCyclotronRun::SetPrimariesPerEvent(G4in 560 void STCyclotronRun::SetPrimariesPerEvent(G4int primaries)
561 { 561 {
562 fPrimariesPerEvent = primaries; 562 fPrimariesPerEvent = primaries;
563 } 563 }
564 564
565 void STCyclotronRun::SetTimePerEvent(G4double 565 void STCyclotronRun::SetTimePerEvent(G4double timePerEvent)
566 { 566 {
567 fTimePerEvent = timePerEvent; 567 fTimePerEvent = timePerEvent;
568 } 568 }
569 569
570 void STCyclotronRun::SetBeamName(G4String beam 570 void STCyclotronRun::SetBeamName(G4String beamName)
571 { 571 {
572 fBeamName = beamName; 572 fBeamName = beamName;
573 } 573 }
574 574
575 void STCyclotronRun::SetBeamCurrent(G4double b 575 void STCyclotronRun::SetBeamCurrent(G4double beamCurrent)
576 { 576 {
577 fBeamCurrent = beamCurrent; 577 fBeamCurrent = beamCurrent;
578 } 578 }
579 579
580 void STCyclotronRun::SetBeamEnergy(G4double be 580 void STCyclotronRun::SetBeamEnergy(G4double beamEnergy)
581 { 581 {
582 fBeamEnergy = beamEnergy; 582 fBeamEnergy = beamEnergy;
583 } 583 }
584 584