File indexing completed on 2025-02-23 09:20:26
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030 #include "STCyclotronRun.hh"
0031 #include "G4RunManager.hh"
0032 #include "G4Event.hh"
0033
0034 #include "G4SDManager.hh"
0035 #include "G4HCofThisEvent.hh"
0036 #include "G4THitsMap.hh"
0037 #include "G4SystemOfUnits.hh"
0038
0039 STCyclotronRun::STCyclotronRun()
0040 : 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.)
0041 { }
0042
0043 STCyclotronRun::~STCyclotronRun()
0044 { }
0045
0046 void STCyclotronRun::Merge(const G4Run* aRun)
0047 {
0048
0049 const STCyclotronRun* localRun = static_cast<const STCyclotronRun*>(aRun);
0050
0051
0052 fTotalEnergyDepositTarget += localRun->fTotalEnergyDepositTarget;
0053 fTotalEnergyDepositFoil += localRun->fTotalEnergyDepositFoil;
0054 fParticleTarget += localRun->fParticleTarget;
0055
0056
0057 if(localRun->fTargetVolume!=0)fTargetVolume = localRun->fTargetVolume;
0058 if(localRun->fFoilVolume!=0)fFoilVolume = localRun->fFoilVolume;
0059 if(localRun->fPrimariesPerEvent!=0)fPrimariesPerEvent = localRun->fPrimariesPerEvent;
0060 if(localRun->fTimePerEvent!=0)fTimePerEvent = localRun->fTimePerEvent;
0061 if(localRun->fTargetThickness!=0)fTargetThickness = localRun->fTargetThickness;
0062 if(localRun->fTargetDiameter!=0)fTargetDiameter = localRun->fTargetDiameter;
0063 if(localRun->fFoilThickness!=0)fFoilThickness = localRun->fFoilThickness;
0064 fBeamName = localRun->fBeamName;
0065 if(localRun->fBeamCurrent!=0.)fBeamCurrent = localRun->fBeamCurrent;
0066 if(localRun->fBeamEnergy!=0.)fBeamEnergy = localRun->fBeamEnergy;
0067
0068
0069
0070
0071 std::map<G4String,G4int>::iterator itSI;
0072 std::map<G4String,G4double>::iterator itSD;
0073 std::map<G4String,G4String>::iterator itSS;
0074 std::map<G4int,G4String>::iterator itIS;
0075
0076
0077 std::map<G4String,G4int> locPrimaryIsotopeCountTarget = localRun->fPrimaryIsotopeCountTarget;
0078 for (itSI = locPrimaryIsotopeCountTarget.begin(); itSI != locPrimaryIsotopeCountTarget.end(); itSI++)
0079 {
0080 G4String name = itSI->first;
0081 G4int count = itSI->second;
0082 fPrimaryIsotopeCountTarget[name] += count;
0083 }
0084
0085 std::map<G4String,G4double> locPrimaryIsotopeTimeTarget = localRun->fPrimaryIsotopeTimeTarget;
0086 for (itSD = locPrimaryIsotopeTimeTarget.begin(); itSD != locPrimaryIsotopeTimeTarget.end(); itSD++)
0087 {
0088 G4String name = itSD->first;
0089 G4double time = itSD->second;
0090 fPrimaryIsotopeTimeTarget[name] = time;
0091 }
0092
0093
0094
0095 std::map<G4String,G4String> locDecayIsotopeCountTarget = localRun->fDecayIsotopeCountTarget;
0096 for (itSS = locDecayIsotopeCountTarget.begin(); itSS != locDecayIsotopeCountTarget.end(); itSS++)
0097 {
0098 G4String nameDaughter = itSS->first;
0099 G4String mum = itSS->second;
0100 fDecayIsotopeCountTarget[nameDaughter] = mum;
0101 }
0102 std::map<G4String,G4double> locDecayIsotopeTimeTarget = localRun->fDecayIsotopeTimeTarget;
0103 for (itSD = locDecayIsotopeTimeTarget.begin(); itSD != locDecayIsotopeTimeTarget.end(); itSD++)
0104 {
0105 G4String nameDaughter = itSD->first;
0106 G4double time = itSD->second;
0107 fDecayIsotopeTimeTarget[nameDaughter] = time;
0108 }
0109 std::map<G4String,G4String> locParticleParent = localRun->fParticleParent;
0110 for (itSS = locParticleParent.begin(); itSS != locParticleParent.end(); itSS++)
0111 {
0112 G4String nameDaughter = itSS->first;
0113 G4String parent = itSS->second;
0114 fParticleParent[nameDaughter] = parent;
0115 }
0116 std::map<G4int,G4String> locIsotopeIDTarget = localRun->fIsotopeIDTarget;
0117 for (itIS = locIsotopeIDTarget.begin(); itIS != locIsotopeIDTarget.end(); itIS++)
0118 {
0119 G4int ID = itIS->first;
0120 G4String name = itIS->second;
0121 fIsotopeIDTarget[ID] = name;
0122 }
0123
0124
0125
0126 std::map<G4String,G4int> locStableIsotopeCountTarget = localRun->fStableIsotopeCountTarget;
0127 for (itSI = locStableIsotopeCountTarget.begin(); itSI != locStableIsotopeCountTarget.end(); itSI++)
0128 {
0129 G4String name = itSI->first;
0130 G4int count = itSI->second;
0131 fStableIsotopeCountTarget[name] += count;
0132 }
0133
0134
0135 std::map<G4String,G4int> locParticleCountTarget = localRun->fParticleCountTarget;
0136 for (itSI = locParticleCountTarget.begin(); itSI != locParticleCountTarget.end(); itSI++)
0137 {
0138 G4String name = itSI->first;
0139 G4int count = itSI->second;
0140 fParticleCountTarget[name] += count;
0141 }
0142
0143
0144 G4Run::Merge(aRun);
0145
0146 }
0147
0148 void STCyclotronRun::EndOfRun(G4double irradiationTime)
0149 {
0150
0151 G4int nbEvents = GetNumberOfEvent();
0152
0153 if (nbEvents == 0) return;
0154
0155
0156
0157
0158 fOutPut.open("Output_General.txt",std::ofstream::out);
0159 fOutPut1.open("Output_ParentIsotopes.txt",std::ofstream::out);
0160 fOutPut2.open("Output_DaughterIsotopes.txt",std::ofstream::out);
0161 fOutPut3.open("Output_OtherParticles.txt",std::ofstream::out);
0162 fOutPut4.open("Output_StableIsotopes.txt",std::ofstream::out);
0163
0164
0165
0166
0167 G4double timePerEvent = fTimePerEvent;
0168 G4double timeForARun = nbEvents*timePerEvent;
0169 G4double minDecay = 0.0001;
0170 G4double maxDecay = 1000000.;
0171
0172
0173
0174
0175
0176
0177 G4int totalPrimaries = fPrimariesPerEvent*nbEvents;
0178 G4double currentFactor;
0179 if(fParticleTarget>0.) currentFactor =(fParticleTarget*1.)/(totalPrimaries*1.);
0180 else currentFactor = 0.;
0181
0182
0183 fOutPut << "//-----------------------------------//" << G4endl;
0184 fOutPut << "// Parameters of the simulation: //" << G4endl;
0185 fOutPut << "//-----------------------------------//" << G4endl;
0186 fOutPut << "Beam parameters: " << G4endl;
0187 fOutPut << fBeamName << " - Name of beam primary particles." << G4endl;
0188 fOutPut << fBeamEnergy << " - Energy of beam primary particles (MeV)." << G4endl;
0189 fOutPut << fBeamCurrent << " - Beam current (Ampere)." << G4endl;
0190 fOutPut << irradiationTime << " - Irradiation time in hour(s)." << G4endl;
0191 fOutPut << currentFactor << " - Current factor." << G4endl;
0192 fOutPut << "//-----------------------------------//" << G4endl;
0193 fOutPut << "Simulation parameters: " << G4endl;
0194 fOutPut << timePerEvent << " - Equivalent time per event (s)." << G4endl;
0195 fOutPut << nbEvents << " - Number of events" << G4endl;
0196 fOutPut << fPrimariesPerEvent << " - Primaries per event" << G4endl;
0197 fOutPut << fPrimariesPerEvent*nbEvents << " - Total number of particles sent." << G4endl;
0198 fOutPut << "//-----------------------------------//" << G4endl;
0199 fOutPut << "Geometry parameters: " << G4endl;
0200 fOutPut << fTargetThickness << " - target thickness (mm)." << G4endl;
0201 fOutPut << fTargetDiameter << " - target diameter (mm)." << G4endl;
0202 fOutPut << fFoilThickness << " - foil thickness (mm)." << G4endl;
0203
0204
0205
0206
0207
0208
0209
0210
0211
0212
0213 std::map<G4String,G4double> fPrimaryIsotopeEOBTarget;
0214 std::map<G4String,G4double> fPrimaryActivityTarget;
0215 std::map<G4String,G4double> fDecayIsotopeEOBTarget;
0216 std::map<G4String,G4double> fDecayActivityTarget;
0217
0218
0219
0220
0221
0222
0223
0224 std::map<G4String,G4int>::iterator it;
0225 G4double primaryActivityTotal = 0.;
0226 G4double decayActivityTotal=0.;
0227
0228 fOutPut1 << "//-----------------------------------//\n"
0229 << "// Data for parent isotopes //\n"
0230 << "//-----------------------------------//\n" << G4endl;
0231
0232 for (it = fPrimaryIsotopeCountTarget.begin(); it != fPrimaryIsotopeCountTarget.end(); it++)
0233 {
0234
0235
0236 G4String name = it->first;
0237 G4double count = (it->second)*currentFactor;
0238 G4double halfLifeTime = fPrimaryIsotopeTimeTarget[name]*10E-10/3600.*std::log(2.);
0239 G4String process = fParticleParent[name];
0240
0241
0242 G4bool store;
0243 if(halfLifeTime > minDecay && halfLifeTime < maxDecay) store = true;
0244 else store = false;
0245
0246
0247
0248 G4double decayConstant = 1/(fPrimaryIsotopeTimeTarget[name]*10E-10);
0249
0250
0251
0252
0253
0254
0255 G4double particlesPerSecond = fPrimaryIsotopeCountTarget[name]*currentFactor/timeForARun;
0256
0257
0258
0259
0260
0261 fPrimaryIsotopeEOBTarget[name] = particlesPerSecond/decayConstant * (1. - std::exp(-irradiationTime*3600*decayConstant));
0262
0263
0264
0265
0266
0267
0268 G4double conv = 2.7E-8;
0269 fPrimaryActivityTarget[name]= fPrimaryIsotopeEOBTarget[name]*decayConstant*conv;
0270
0271 if(store)
0272 {
0273
0274
0275
0276
0277
0278 primaryActivityTotal = primaryActivityTotal + fPrimaryActivityTarget[name];
0279 }
0280
0281
0282
0283
0284
0285
0286 if(store)
0287 {
0288 fOutPut1 << name << " - name of parent isotope." << G4endl;
0289 fOutPut1 << count/currentFactor << " - number of isotopes created during the simulation." << G4endl;
0290 fOutPut1 << decayConstant << " - decay constant in s-1." << G4endl;
0291 fOutPut1 << halfLifeTime << " - half life time in hour(s)." << G4endl;
0292 fOutPut1 << process << " - creation process." << G4endl;
0293 fOutPut1 << particlesPerSecond << " - isotope per sec." << G4endl;
0294 fOutPut1 << fPrimaryIsotopeEOBTarget[name] << " - yield EOB." << G4endl;
0295 fOutPut1 << fPrimaryActivityTarget[name] << " - activity (mCi) at the EOB." << G4endl;
0296 fOutPut1 << "------------------------" << G4endl;
0297 }
0298
0299 }
0300
0301
0302
0303
0304
0305
0306 fOutPut2 << "//-----------------------------------//\n"
0307 << "// Data for daughter isotopes //\n"
0308 << "//-----------------------------------//\n" << G4endl;
0309
0310 std::map<G4String,G4String>::iterator it1;
0311
0312 for (it1 = fDecayIsotopeCountTarget.begin(); it1 != fDecayIsotopeCountTarget.end(); it1++)
0313 {
0314
0315
0316 G4String nameDaughter = it1->first;
0317 G4String nameMum = it1->second;
0318
0319 G4double halfLifeTimeMum = fDecayIsotopeTimeTarget[nameDaughter]*10E10/3600;
0320 G4double halfLifeTimeDaughter = fPrimaryIsotopeTimeTarget[nameMum]*10E10/3600;
0321
0322 G4bool store;
0323 if(halfLifeTimeMum > minDecay && halfLifeTimeMum < maxDecay &&
0324 halfLifeTimeDaughter > minDecay && halfLifeTimeDaughter < maxDecay){store=true;}
0325 else{store=false;}
0326
0327
0328
0329
0330
0331
0332
0333
0334
0335 G4double decayConstantMum = fPrimaryIsotopeTimeTarget[nameMum]
0336 ? 1/(fPrimaryIsotopeTimeTarget[nameMum]*10.E-10) : 0.;
0337 G4double decayConstantDaughter = fDecayIsotopeTimeTarget[nameDaughter]
0338 ? 1/(fDecayIsotopeTimeTarget[nameDaughter]*10.E-10) : 0.;
0339
0340
0341
0342
0343
0344
0345 G4double particlesPerSecond = fPrimaryIsotopeCountTarget[nameMum]*currentFactor/timeForARun;
0346
0347
0348
0349
0350
0351 fDecayIsotopeEOBTarget[nameDaughter] = particlesPerSecond*((1 - std::exp(-irradiationTime*3600*decayConstantDaughter))/decayConstantDaughter + (std::exp(-irradiationTime*3600*decayConstantDaughter) - std::exp(-irradiationTime*3600*decayConstantMum))/(decayConstantDaughter-decayConstantMum));
0352
0353
0354
0355
0356
0357
0358
0359 G4double conv = 2.7E-8;
0360 fDecayActivityTarget[nameDaughter]= fDecayIsotopeEOBTarget[nameDaughter]*decayConstantDaughter*conv;
0361
0362 if(store)
0363 {
0364 decayActivityTotal = decayActivityTotal + fDecayActivityTarget[nameDaughter];
0365 }
0366
0367 if(store)
0368 {
0369 fOutPut2 << nameDaughter << " - name of daughter isotope." << G4endl;
0370 fOutPut2 << nameMum << " - name of parent isotope." << G4endl;
0371 fOutPut2 << decayConstantDaughter << " - decay constant of daughter in s-1." << G4endl;
0372 fOutPut2 << decayConstantMum << " - decay constant of mum in s-1." << G4endl;
0373 fOutPut2 << halfLifeTimeDaughter << " - half life time of daughter in hour(s)." << G4endl;
0374 fOutPut2 << halfLifeTimeMum << " - half life time of mum in hour(s)." << G4endl;
0375 fOutPut2 << particlesPerSecond << " - isotope per sec." << G4endl;
0376 fOutPut2 << fDecayIsotopeEOBTarget[nameDaughter] << " - yield at the EOB." << G4endl;
0377 fOutPut2 << fDecayActivityTarget[nameDaughter] << " - activity (mCi) at the EOB." << G4endl;
0378 fOutPut2 << "------------------------" << G4endl;
0379 }
0380
0381 }
0382
0383
0384
0385
0386
0387 fOutPut3 << "//-----------------------------------//\n"
0388 << "// Data for other particles //\n"
0389 << "//-----------------------------------//" << G4endl;
0390
0391 std::map<G4String, G4int>::iterator it3;
0392 for(it3=fParticleCountTarget.begin(); it3!= fParticleCountTarget.end(); it3++)
0393 {
0394 G4String name = it3->first;
0395 G4double number = it3->second;
0396 fOutPut3 << name << " - name of the particle" << G4endl;
0397 fOutPut3 << number << " - number of particles" << G4endl;
0398 fOutPut3 << "------------------------" << G4endl;
0399 }
0400
0401
0402
0403
0404 fOutPut4 << "//-----------------------------------//\n"
0405 << "// Data for stable isotopes //\n"
0406 << "//-----------------------------------//\n" << G4endl;
0407
0408 std::map<G4String, G4int>::iterator it6;
0409 for(it6=fStableIsotopeCountTarget.begin();it6!=fStableIsotopeCountTarget.end();it6++)
0410 {
0411 G4String isotope = it6 ->first;
0412 G4int number = it6 -> second;
0413 fOutPut4 << isotope << " - name of the isotope" << G4endl;
0414 fOutPut4 << number << " - number of isotopes" << G4endl;
0415 fOutPut4 << "------------------------" << G4endl;
0416 }
0417
0418
0419
0420
0421 fPrimaryIsotopeEOBTarget.clear();
0422 fPrimaryActivityTarget.clear();
0423 fDecayIsotopeEOBTarget.clear();
0424 fDecayActivityTarget.clear();
0425
0426
0427
0428 fPrimaryIsotopeCountTarget.clear();
0429 fPrimaryIsotopeTimeTarget.clear();
0430 fDecayIsotopeCountTarget.clear();
0431 fDecayIsotopeTimeTarget.clear();
0432 fParticleParent.clear();
0433 fParticleCountTarget.clear();
0434 fStableIsotopeCountTarget.clear();
0435 fIsotopeIDTarget.clear();
0436
0437
0438
0439
0440
0441
0442 G4double totalEnergyDepositTargetEOB = fTotalEnergyDepositTarget/timeForARun * irradiationTime * 3600.;
0443 G4double totalEnergyDepositTargetPerSecond = fTotalEnergyDepositTarget/timeForARun;
0444
0445 G4double heatTarget = totalEnergyDepositTargetPerSecond/fTargetVolume * 1.60E-13;
0446 G4double heatFoil = fTotalEnergyDepositFoil / fFoilVolume * 1.60E-13;
0447
0448
0449
0450 fOutPut << "//-------------------------------------------------//\n"
0451 << "// Heating, total activity and process data //\n"
0452 << "//-------------------------------------------------//" << G4endl;
0453
0454 fOutPut << "Total heating in the target : "
0455 << heatTarget << " W/mm3" << G4endl;
0456 fOutPut << "The total heating during the irradiation is " << totalEnergyDepositTargetEOB << "J/mm3" << G4endl;
0457 fOutPut << "Total heating in the foil : " << heatFoil << " W/mm3" << G4endl;
0458
0459
0460 fOutPut.close();
0461 fOutPut1.close();
0462 fOutPut2.close();
0463 fOutPut3.close();
0464 fOutPut4.close();
0465
0466 }
0467
0468
0469
0470
0471
0472
0473 void STCyclotronRun::PrimaryIsotopeCountTarget(G4String name,G4double time)
0474 {
0475 fPrimaryIsotopeCountTarget[name]++;
0476 fPrimaryIsotopeTimeTarget[name]=time;
0477 }
0478
0479 void STCyclotronRun::CountStableIsotopes(G4String name)
0480 {
0481 fStableIsotopeCountTarget[name]++;
0482 }
0483
0484 void STCyclotronRun::DecayIsotopeCountTarget(G4String nameDaughter,G4String mum, G4double time)
0485 {
0486 fDecayIsotopeCountTarget[nameDaughter]=mum;
0487 fDecayIsotopeTimeTarget[nameDaughter]=time;
0488 }
0489
0490 void STCyclotronRun::ParticleParent(G4String isotope, G4String parent)
0491 {
0492 fParticleParent[isotope]=parent;
0493 }
0494
0495
0496
0497 void STCyclotronRun::ParticleCountTarget(G4String name)
0498 {
0499 fParticleCountTarget[name]++;
0500 }
0501
0502
0503
0504
0505
0506
0507
0508
0509 void STCyclotronRun::StoreIsotopeID(G4int ID, G4String name)
0510 {
0511 fIsotopeIDTarget[ID]=name;
0512 }
0513
0514 std::map<G4int,G4String> STCyclotronRun::GetIsotopeID()
0515 {
0516 return fIsotopeIDTarget;
0517 }
0518
0519
0520 void STCyclotronRun::EnergyDepositionTarget(G4double edep)
0521 {
0522 fTotalEnergyDepositTarget += edep;
0523 }
0524
0525 void STCyclotronRun::EnergyDepositionFoil(G4double edep)
0526 {
0527 fTotalEnergyDepositFoil += edep;
0528 }
0529
0530 void STCyclotronRun::CountParticlesTarget()
0531 {
0532 fParticleTarget++;
0533 }
0534
0535 void STCyclotronRun::SetFoilVolume(G4double foilVolume)
0536 {
0537 fFoilVolume = foilVolume;
0538 }
0539
0540 void STCyclotronRun::SetFoilThickness(G4double foilThickness)
0541 {
0542 fFoilThickness = foilThickness;
0543 }
0544
0545 void STCyclotronRun::SetTargetVolume(G4double targetVolume)
0546 {
0547 fTargetVolume = targetVolume;
0548 }
0549
0550 void STCyclotronRun::SetTargetThickness(G4double targetThickness)
0551 {
0552 fTargetThickness = targetThickness;
0553 }
0554
0555 void STCyclotronRun::SetTargetDiameter(G4double targetDiameter)
0556 {
0557 fTargetDiameter = targetDiameter;
0558 }
0559
0560 void STCyclotronRun::SetPrimariesPerEvent(G4int primaries)
0561 {
0562 fPrimariesPerEvent = primaries;
0563 }
0564
0565 void STCyclotronRun::SetTimePerEvent(G4double timePerEvent)
0566 {
0567 fTimePerEvent = timePerEvent;
0568 }
0569
0570 void STCyclotronRun::SetBeamName(G4String beamName)
0571 {
0572 fBeamName = beamName;
0573 }
0574
0575 void STCyclotronRun::SetBeamCurrent(G4double beamCurrent)
0576 {
0577 fBeamCurrent = beamCurrent;
0578 }
0579
0580 void STCyclotronRun::SetBeamEnergy(G4double beamEnergy)
0581 {
0582 fBeamEnergy = beamEnergy;
0583 }