| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // Author: F. Poignant, floriane.poignant@gmai
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.),fTot
41 { }
42
43 STCyclotronRun::~STCyclotronRun()
44 { }
45
46 void STCyclotronRun::Merge(const G4Run* aRun)
47 {
48
49 const STCyclotronRun* localRun = static_cast
50
51 //Merging cumulable variables
52 fTotalEnergyDepositTarget += localRun->fTota
53 fTotalEnergyDepositFoil += localRun->fTotalE
54 fParticleTarget += localRun->fParticleTarget
55
56 //Constant over the different runs
57 if(localRun->fTargetVolume!=0)fTargetVolume
58 if(localRun->fFoilVolume!=0)fFoilVolume = lo
59 if(localRun->fPrimariesPerEvent!=0)fPrimarie
60 if(localRun->fTimePerEvent!=0)fTimePerEvent
61 if(localRun->fTargetThickness!=0)fTargetThic
62 if(localRun->fTargetDiameter!=0)fTargetDiame
63 if(localRun->fFoilThickness!=0)fFoilThicknes
64 fBeamName = localRun->fBeamName;
65 if(localRun->fBeamCurrent!=0.)fBeamCurrent =
66 if(localRun->fBeamEnergy!=0.)fBeamEnergy = l
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> locPrimaryIsotopeCo
78 for (itSI = locPrimaryIsotopeCountTarget.beg
79 {
80 G4String name = itSI->first;
81 G4int count = itSI->second;
82 fPrimaryIsotopeCountTarget[name] += coun
83 }
84
85 std::map<G4String,G4double> locPrimaryIsotop
86 for (itSD = locPrimaryIsotopeTimeTarget.begi
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> fIsotopeIDTar
95 std::map<G4String,G4String> locDecayIsotopeC
96 for (itSS = locDecayIsotopeCountTarget.begin
97 {
98 G4String nameDaughter = itSS->first;
99 G4String mum = itSS->second;
100 fDecayIsotopeCountTarget[nameDaughter] =
101 }
102 std::map<G4String,G4double> locDecayIsotopeT
103 for (itSD = locDecayIsotopeTimeTarget.begin
104 {
105 G4String nameDaughter = itSD->first;
106 G4double time = itSD->second;
107 fDecayIsotopeTimeTarget[nameDaughter] =
108 }
109 std::map<G4String,G4String> locParticlePare
110 for (itSS = locParticleParent.begin(); itSS
111 {
112 G4String nameDaughter = itSS->first;
113 G4String parent = itSS->second;
114 fParticleParent[nameDaughter] = parent;
115 }
116 std::map<G4int,G4String> locIsotopeIDTarget
117 for (itIS = locIsotopeIDTarget.begin(); itI
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> locStableIsotopeCou
127 for (itSI = locStableIsotopeCountTarget.begi
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> locParticleCountTar
136 for (itSI = locParticleCountTarget.begin();
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 irradia
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::ofstr
159 fOutPut1.open("Output_ParentIsotopes.txt",st
160 fOutPut2.open("Output_DaughterIsotopes.txt",
161 fOutPut3.open("Output_OtherParticles.txt",st
162 fOutPut4.open("Output_StableIsotopes.txt",st
163
164 //------------------------------------------
165 // Calculates the equivalent time for a giv
166 //------------------------------------------
167 G4double timePerEvent = fTimePerEvent; //in
168 G4double timeForARun = nbEvents*timePerEvent
169 G4double minDecay = 0.0001; //in seconds
170 G4double maxDecay = 1000000.; //in seconds
171
172 //------------------------------------------
173 // Rescale the value of the beam current
174 // the loss of primary particles due to t
175 //------------------------------------------
176
177 G4int totalPrimaries = fPrimariesPerEvent*nb
178 G4double currentFactor;
179 if(fParticleTarget>0.) currentFactor =(fPart
180 else currentFactor = 0.;
181
182
183 fOutPut << "//------------------------------
184 fOutPut << "// Parameters of the simulati
185 fOutPut << "//------------------------------
186 fOutPut << "Beam parameters: " << G4endl;
187 fOutPut << fBeamName << " - Name of beam pri
188 fOutPut << fBeamEnergy << " - Energy of beam
189 fOutPut << fBeamCurrent << " - Beam current
190 fOutPut << irradiationTime << " - Irradiatio
191 fOutPut << currentFactor << " - Current fact
192 fOutPut << "//------------------------------
193 fOutPut << "Simulation parameters: " << G4en
194 fOutPut << timePerEvent << " - Equivalent ti
195 fOutPut << nbEvents << " - Number of events"
196 fOutPut << fPrimariesPerEvent << " - Primari
197 fOutPut << fPrimariesPerEvent*nbEvents << "
198 fOutPut << "//------------------------------
199 fOutPut << "Geometry parameters: " << G4endl
200 fOutPut << fTargetThickness << " - target th
201 fOutPut << fTargetDiameter << " - target dia
202 fOutPut << fFoilThickness << " - foil thickn
203 //Add particle type, particle energy, beam d
204
205 //target material???
206
207 ////////////////////////////////////////////
208 //////////Calculation of the number of isoto
209 ////////////////////////////////////////////
210
211
212 //Maps to fill
213 std::map<G4String,G4double> fPrimaryIsotopeE
214 std::map<G4String,G4double> fPrimaryActivity
215 std::map<G4String,G4double> fDecayIsotopeEOB
216 std::map<G4String,G4double> fDecayActivityTa
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 << "//-----------------------------
229 << "// Data for parent isotopes /
230 << "//-----------------------------------/
231
232 for (it = fPrimaryIsotopeCountTarget.begin()
233 {
234
235
236 G4String name = it->first;
237 G4double count = (it->second)*currentF
238 G4double halfLifeTime = fPrimaryIsotopeT
239 G4String process = fParticleParent[name]
240
241 //Only store isotopes with a life time b
242 G4bool store;
243 if(halfLifeTime > minDecay && halfLifeTi
244 else store = false;
245
246
247 //Calculation of the yield (s-1)
248 G4double decayConstant = 1/(fPrimaryIsot
249
250
251 //--------------------------------------
252 // Number of particles per second
253 //--------------------------------------
254
255 G4double particlesPerSecond = fPrimaryIs
256
257 //--------------------------------------
258 // Calculation yield EOB
259 //--------------------------------------
260
261 fPrimaryIsotopeEOBTarget[name] = particl
262
263 //--------------------------------------
264 // Calculation of the activity
265 // conversion factor Bq to mCi
266 //--------------------------------------
267
268 G4double conv = 2.7E-8;
269 fPrimaryActivityTarget[name]= fPrimaryIs
270
271 if(store)
272 {
273
274 //----------------------------------------
275 // Incrementation for total primary act
276 //----------------------------------------
277
278 primaryActivityTotal = primaryActivityTota
279 }
280
281
282 //---------------------------//
283 // Printing out results //
284 //---------------------------//
285
286 if(store)
287 {
288 fOutPut1 << name << " - name of parent iso
289 fOutPut1 << count/currentFactor << " - num
290 fOutPut1 << decayConstant << " - decay con
291 fOutPut1 << halfLifeTime << " - half life
292 fOutPut1 << process << " - creation proces
293 fOutPut1 << particlesPerSecond << " - isot
294 fOutPut1 << fPrimaryIsotopeEOBTarget[name]
295 fOutPut1 << fPrimaryActivityTarget[name] <
296 fOutPut1 << "------------------------" <<
297 }
298
299 }
300
301
302 //------------------------------------------
303 // CASE 2 : isotopes from primary isotopes
304 //------------------------------------------
305
306 fOutPut2 << "//-----------------------------
307 << "// Data for daughter isotopes /
308 << "//-----------------------------------/
309
310 std::map<G4String,G4String>::iterator it1;
311
312 for (it1 = fDecayIsotopeCountTarget.begin();
313 {
314
315
316 G4String nameDaughter = it1->first;
317 G4String nameMum = it1->second;
318
319 G4double halfLifeTimeMum = fDecayIsotope
320 G4double halfLifeTimeDaughter = fPrimary
321
322 G4bool store;
323 if(halfLifeTimeMum > minDecay && halfLif
324 halfLifeTimeDaughter > minDecay && halfLife
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 b
333 //--------------------------------------
334
335 G4double decayConstantMum = fPrimaryIsot
336 ? 1/(fPrimaryI
337 G4double decayConstantDaughter = fDecayI
338 ? 1/(fDec
339
340
341 //--------------------------------------
342 // Number of particles per secon
343 //--------------------------------------
344
345 G4double particlesPerSecond = fPrimaryIs
346
347 //--------------------------------------
348 // Number of particles at the EOB
349 //--------------------------------------
350
351 fDecayIsotopeEOBTarget[nameDaughter] = p
352
353
354 //--------------------------------------
355 // Calculation of activity
356 // conversion factor Bq to mCu
357 //--------------------------------------
358
359 G4double conv = 2.7E-8;
360 fDecayActivityTarget[nameDaughter]= fDec
361
362 if(store)
363 {
364 decayActivityTotal = decayActivityTotal +
365 }
366
367 if(store)
368 {
369 fOutPut2 << nameDaughter << " - name of da
370 fOutPut2 << nameMum << " - name of parent
371 fOutPut2 << decayConstantDaughter << " - d
372 fOutPut2 << decayConstantMum << " - decay
373 fOutPut2 << halfLifeTimeDaughter << " - ha
374 fOutPut2 << halfLifeTimeMum << " - half li
375 fOutPut2 << particlesPerSecond << " - iso
376 fOutPut2 << fDecayIsotopeEOBTarget[nameDau
377 fOutPut2 << fDecayActivityTarget[nameDaugh
378 fOutPut2 << "------------------------" <<
379 }
380
381 }
382
383 //------------------------------------------
384 // Particles created, other than nuclei
385 //------------------------------------------
386
387 fOutPut3 << "//-----------------------------
388 << "// Data for other particles /
389 << "//-----------------------------------/
390
391 std::map<G4String, G4int>::iterator it3;
392 for(it3=fParticleCountTarget.begin(); it3!=
393 {
394 G4String name = it3->first;
395 G4double number = it3->second;
396 fOutPut3 << name << " - name of the part
397 fOutPut3 << number << " - number of part
398 fOutPut3 << "------------------------" <
399 }
400
401
402
403
404 fOutPut4 << "//-----------------------------
405 << "// Data for stable isotopes /
406 << "//-----------------------------------/
407
408 std::map<G4String, G4int>::iterator it6;
409 for(it6=fStableIsotopeCountTarget.begin();it
410 {
411 G4String isotope = it6 ->first;
412 G4int number = it6 -> second;
413 fOutPut4 << isotope << " - name of the i
414 fOutPut4 << number << " - number of isot
415 fOutPut4 << "------------------------" <
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 = fTota
443 G4double totalEnergyDepositTargetPerSecond =
444 //Heat calculation in W/mm3
445 G4double heatTarget = totalEnergyDepositTarg
446 G4double heatFoil = fTotalEnergyDepositFoil
447
448
449 //Output data in a .txt file
450 fOutPut << "//------------------------------
451 << "// Heating, total activity and proc
452 << "//-------------------------------------
453
454 fOutPut << "Total heating in the target : "
455 << heatTarget << " W/mm3" << G4endl;
456 fOutPut << "The total heating during the irr
457 fOutPut << "Total heating in the foil : " <<
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
472
473 void STCyclotronRun::PrimaryIsotopeCountTarget
474 {
475 fPrimaryIsotopeCountTarget[name]++;
476 fPrimaryIsotopeTimeTarget[name]=time;
477 }
478 //------------------------------------
479 void STCyclotronRun::CountStableIsotopes(G4Str
480 {
481 fStableIsotopeCountTarget[name]++;
482 }
483 //------------------------------------
484 void STCyclotronRun::DecayIsotopeCountTarget(G
485 {
486 fDecayIsotopeCountTarget[nameDaughter]=mum;
487 fDecayIsotopeTimeTarget[nameDaughter]=time;
488 }
489 //------------------------------------
490 void STCyclotronRun::ParticleParent(G4String i
491 {
492 fParticleParent[isotope]=parent;
493 }
494 //
495 //-----> Count other particles
496 //------------------------------------
497 void STCyclotronRun::ParticleCountTarget(G4Str
498 {
499 fParticleCountTarget[name]++;
500 }
501
502
503
504 //--------------------------------------------
505 // Accumulation functions for maps used only d
506 //
507 //-----> Isotope ID to obtain the "mother isot
508 //------------------------------------
509 void STCyclotronRun::StoreIsotopeID(G4int ID,
510 {
511 fIsotopeIDTarget[ID]=name;
512 }
513 //
514 std::map<G4int,G4String> STCyclotronRun::GetIs
515 {
516 return fIsotopeIDTarget;
517 }
518
519
520 void STCyclotronRun::EnergyDepositionTarget(G4
521 {
522 fTotalEnergyDepositTarget += edep;
523 }
524
525 void STCyclotronRun::EnergyDepositionFoil(G4do
526 {
527 fTotalEnergyDepositFoil += edep;
528 }
529
530 void STCyclotronRun::CountParticlesTarget()
531 {
532 fParticleTarget++;
533 }
534
535 void STCyclotronRun::SetFoilVolume(G4double fo
536 {
537 fFoilVolume = foilVolume;
538 }
539
540 void STCyclotronRun::SetFoilThickness(G4double
541 {
542 fFoilThickness = foilThickness;
543 }
544
545 void STCyclotronRun::SetTargetVolume(G4double
546 {
547 fTargetVolume = targetVolume;
548 }
549
550 void STCyclotronRun::SetTargetThickness(G4doub
551 {
552 fTargetThickness = targetThickness;
553 }
554
555 void STCyclotronRun::SetTargetDiameter(G4doubl
556 {
557 fTargetDiameter = targetDiameter;
558 }
559
560 void STCyclotronRun::SetPrimariesPerEvent(G4in
561 {
562 fPrimariesPerEvent = primaries;
563 }
564
565 void STCyclotronRun::SetTimePerEvent(G4double
566 {
567 fTimePerEvent = timePerEvent;
568 }
569
570 void STCyclotronRun::SetBeamName(G4String beam
571 {
572 fBeamName = beamName;
573 }
574
575 void STCyclotronRun::SetBeamCurrent(G4double b
576 {
577 fBeamCurrent = beamCurrent;
578 }
579
580 void STCyclotronRun::SetBeamEnergy(G4double be
581 {
582 fBeamEnergy = beamEnergy;
583 }
584