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