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