| 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 //
27
28 #include "G4AdjointCSManager.hh"
29
30 #include "G4AdjointCSMatrix.hh"
31 #include "G4AdjointElectron.hh"
32 #include "G4AdjointGamma.hh"
33 #include "G4AdjointInterpolator.hh"
34 #include "G4AdjointProton.hh"
35 #include "G4Electron.hh"
36 #include "G4Element.hh"
37 #include "G4ElementTable.hh"
38 #include "G4Gamma.hh"
39 #include "G4ParticleDefinition.hh"
40 #include "G4PhysicalConstants.hh"
41 #include "G4PhysicsLogVector.hh"
42 #include "G4PhysicsTable.hh"
43 #include "G4ProductionCutsTable.hh"
44 #include "G4Proton.hh"
45 #include "G4SystemOfUnits.hh"
46 #include "G4VEmAdjointModel.hh"
47 #include "G4VEmProcess.hh"
48 #include "G4VEnergyLossProcess.hh"
49
50 G4ThreadLocal G4AdjointCSManager* G4AdjointCSM
51
52 constexpr G4double G4AdjointCSManager::fTmin;
53 constexpr G4double G4AdjointCSManager::fTmax;
54 constexpr G4int G4AdjointCSManager::fNbins;
55
56 //////////////////////////////////////////////
57 G4AdjointCSManager* G4AdjointCSManager::GetAdj
58 {
59 if(fInstance == nullptr)
60 {
61 static G4ThreadLocalSingleton<G4AdjointCSM
62 fInstance = inst.Instance();
63 }
64 return fInstance;
65 }
66
67 //////////////////////////////////////////////
68 G4AdjointCSManager::G4AdjointCSManager()
69 {
70 RegisterAdjointParticle(G4AdjointElectron::A
71 RegisterAdjointParticle(G4AdjointGamma::Adjo
72 RegisterAdjointParticle(G4AdjointProton::Adj
73 }
74
75 //////////////////////////////////////////////
76 G4AdjointCSManager::~G4AdjointCSManager()
77 {
78 for (auto& p : fAdjointCSMatricesForProdToPr
79 for (auto p1 : p) {
80 if (p1) {
81 delete p1;
82 p1 = nullptr;
83 }
84 }
85 p.clear();
86 }
87 fAdjointCSMatricesForProdToProj.clear();
88
89 for (auto& p : fAdjointCSMatricesForScatProj
90 for (auto p1 : p) {
91 if (p1) {
92 delete p1;
93 p1 = nullptr;
94 }
95 }
96 p.clear();
97 }
98 fAdjointCSMatricesForScatProjToProj.clear();
99
100 for (auto p : fAdjointModels) {
101 if (p) {
102 delete p;
103 p = nullptr;
104 }
105 }
106 fAdjointModels.clear();
107
108 for (auto p : fTotalAdjSigmaTable) {
109 p->clearAndDestroy();
110 delete p;
111 p = nullptr;
112 }
113 fTotalAdjSigmaTable.clear();
114
115 for (auto p : fSigmaTableForAdjointModelScat
116 p->clearAndDestroy();
117 delete p;
118 p = nullptr;
119 }
120 fSigmaTableForAdjointModelScatProjToProj.cle
121
122 for (auto p : fSigmaTableForAdjointModelProd
123 p->clearAndDestroy();
124 delete p;
125 p = nullptr;
126 }
127 fSigmaTableForAdjointModelProdToProj.clear()
128
129 for (auto p : fTotalFwdSigmaTable) {
130 p->clearAndDestroy();
131 delete p;
132 p = nullptr;
133 }
134 fTotalFwdSigmaTable.clear();
135
136 for (auto p : fForwardProcesses) {
137 delete p;
138 p = nullptr;
139 }
140 fForwardProcesses.clear();
141
142 for (auto p : fForwardLossProcesses) {
143 delete p;
144 p = nullptr;
145 }
146 fForwardLossProcesses.clear();
147 }
148
149 //////////////////////////////////////////////
150 std::size_t G4AdjointCSManager::RegisterEmAdjo
151 {
152 fAdjointModels.push_back(aModel);
153 fSigmaTableForAdjointModelScatProjToProj.pus
154 fSigmaTableForAdjointModelProdToProj.push_ba
155 return fAdjointModels.size() - 1;
156 }
157
158 //////////////////////////////////////////////
159 void G4AdjointCSManager::RegisterEmProcess(G4V
160 G4P
161 {
162 G4ParticleDefinition* anAdjPartDef =
163 GetAdjointParticleEquivalent(aFwdPartDef);
164 if(anAdjPartDef && aProcess)
165 {
166 RegisterAdjointParticle(anAdjPartDef);
167
168 for(std::size_t i = 0; i < fAdjointParticl
169 {
170 if(anAdjPartDef->GetParticleName() ==
171 fAdjointParticlesInAction[i]->GetPart
172 fForwardProcesses[i]->push_back(aProce
173 }
174 }
175 }
176
177 //////////////////////////////////////////////
178 void G4AdjointCSManager::RegisterEnergyLossPro
179 G4VEnergyLossProcess* aProcess, G4ParticleDe
180 {
181 G4ParticleDefinition* anAdjPartDef =
182 GetAdjointParticleEquivalent(aFwdPartDef);
183 if(anAdjPartDef && aProcess)
184 {
185 RegisterAdjointParticle(anAdjPartDef);
186 for(std::size_t i = 0; i < fAdjointParticl
187 {
188 if(anAdjPartDef->GetParticleName() ==
189 fAdjointParticlesInAction[i]->GetPart
190 fForwardLossProc
191
192 }
193 }
194 }
195
196 //////////////////////////////////////////////
197 void G4AdjointCSManager::RegisterAdjointPartic
198 {
199 G4bool found = false;
200 for(auto p : fAdjointParticlesInAction)
201 {
202 if(p->GetParticleName() == aPartDef->GetPa
203 {
204 found = true;
205 }
206 }
207 if(!found)
208 {
209 fForwardLossProcesses.push_back(new std::v
210 fTotalFwdSigmaTable.push_back(new G4Physic
211 fTotalAdjSigmaTable.push_back(new G4Physic
212 fForwardProcesses.push_back(new std::vecto
213 fAdjointParticlesInAction.push_back(aPartD
214 fEminForFwdSigmaTables.push_back(std::vect
215 fEminForAdjSigmaTables.push_back(std::vect
216 fEkinofFwdSigmaMax.push_back(std::vector<G
217 fEkinofAdjSigmaMax.push_back(std::vector<G
218 }
219 }
220
221 //////////////////////////////////////////////
222 void G4AdjointCSManager::BuildCrossSectionMatr
223 {
224 if(fCSMatricesBuilt)
225 return;
226 // The Tcut and Tmax matrices will be comput
227 // When Tcut changes, some PhysicsTable will
228 // for each MaterialCutCouple but not all th
229 // The Tcut defines a lower limit in the ene
230 // scattering. In the Projectile to Scattere
231 // E_ScatProj<E_Proj-Tcut
232 // Therefore in the adjoint case we have
233 // Eproj> E_ScatProj+Tcut
234 // This implies that when computing the adjo
235 // Epro from E_ScatProj+Tcut to Emax
236 // In the Projectile to Secondary case Tcut
237 // Esecond should be greater than Tcut to ha
238 // adjoint process.
239 // To avoid recomputing matrices for all cha
240 // we propose to compute the matrices only o
241 // and then to interpolate the probability f
242 // G4VAdjointEmModel)
243
244 fAdjointCSMatricesForScatProjToProj.clear();
245 fAdjointCSMatricesForProdToProj.clear();
246 const G4ElementTable* theElementTable = G4
247 const G4MaterialTable* theMaterialTable = G4
248
249 G4cout << "========== Computation of cross s
250 "models =========="
251 << G4endl;
252 for(const auto& aModel : fAdjointModels)
253 {
254 G4cout << "Build adjoint cross section mat
255 << G4endl;
256 if(aModel->GetUseMatrix())
257 {
258 std::vector<G4AdjointCSMatrix*>* aListOf
259 new std::vector<G4AdjointCSMatrix*>();
260 std::vector<G4AdjointCSMatrix*>* aListOf
261 new std::vector<G4AdjointCSMatrix*>();
262 if(aModel->GetUseMatrixPerElement())
263 {
264 if(aModel->GetUseOnlyOneMatrixForAllEl
265 {
266 std::vector<G4AdjointCSMatrix*> two_
267 BuildCrossSectionsModelAndElement(
268 aListOfMat1->push_back(two_matrices[
269 aListOfMat2->push_back(two_matrices[
270 }
271 else
272 {
273 for(const auto& anElement : *theElem
274 {
275 G4int Z = G4lrint(anElement->GetZ(
276 G4int A = G4lrint(anElement->GetN(
277 std::vector<G4AdjointCSMatrix*> tw
278 BuildCrossSectionsModelAndElemen
279 aListOfMat1->push_back(two_matrice
280 aListOfMat2->push_back(two_matrice
281 }
282 }
283 }
284 else
285 { // Per material case
286 for(const auto& aMaterial : *theMateri
287 {
288 std::vector<G4AdjointCSMatrix*> two_
289 BuildCrossSectionsModelAndMaterial
290 aListOfMat1->push_back(two_matrices[
291 aListOfMat2->push_back(two_matrices[
292 }
293 }
294 fAdjointCSMatricesForProdToProj.push_bac
295 fAdjointCSMatricesForScatProjToProj.push
296 aModel->SetCSMatrices(aListOfMat1, aList
297 }
298 else
299 {
300 G4cout << "The model " << aModel->GetNam
301 << " does not use cross section m
302 std::vector<G4AdjointCSMatrix*> two_empt
303 fAdjointCSMatricesForProdToProj.push_bac
304 fAdjointCSMatricesForScatProjToProj.push
305 }
306 }
307 G4cout << " All adjoint cross s
308 << G4endl;
309 G4cout
310 << "======================================
311 << G4endl;
312
313 fCSMatricesBuilt = true;
314 }
315
316 //////////////////////////////////////////////
317 void G4AdjointCSManager::BuildTotalSigmaTables
318 {
319 if(fSigmaTableBuilt)
320 return;
321
322 const G4ProductionCutsTable* theCoupleTable
323 G4ProductionCutsTable::GetProductionCutsTa
324
325 // Prepare the Sigma table for all AdjointEM
326 for(std::size_t i = 0; i < fAdjointModels.si
327 {
328 fSigmaTableForAdjointModelScatProjToProj[i
329 fSigmaTableForAdjointModelProdToProj[i]->c
330 for(std::size_t j = 0; j < theCoupleTable-
331 {
332 fSigmaTableForAdjointModelScatProjToProj
333 new G4PhysicsLogVector(fTmin, fTmax, f
334 fSigmaTableForAdjointModelProdToProj[i]-
335 new G4PhysicsLogVector(fTmin, fTmax, f
336 }
337 }
338
339 for(std::size_t i = 0; i < fAdjointParticles
340 {
341 G4ParticleDefinition* thePartDef = fAdjoin
342 DefineCurrentParticle(thePartDef);
343 fTotalFwdSigmaTable[i]->clearAndDestroy();
344 fTotalAdjSigmaTable[i]->clearAndDestroy();
345 fEminForFwdSigmaTables[i].clear();
346 fEminForAdjSigmaTables[i].clear();
347 fEkinofFwdSigmaMax[i].clear();
348 fEkinofAdjSigmaMax[i].clear();
349
350 for(std::size_t j = 0; j < theCoupleTable-
351 {
352 const G4MaterialCutsCouple* couple =
353 theCoupleTable->GetMaterialCutsCouple(
354
355 // make first the total fwd CS table for
356 G4PhysicsVector* aVector = new G4Physics
357 G4bool Emin_found = false;
358 G4double sigma_max = 0.;
359 G4double e_sigma_max = 0.;
360 for(std::size_t l = 0; l < fNbins; ++l)
361 {
362 G4double totCS = 0.;
363 G4double e = aVector->Energy(l);
364 for(std::size_t k = 0; k < fForwardPro
365 {
366 totCS += (*fForwardProcesses[i])[k]-
367 }
368 for(std::size_t k = 0; k < fForwardLos
369 {
370 if(thePartDef == fAdjIon)
371 { // e is considered already as the
372 std::size_t mat_index = couple->Ge
373 G4VEmModel* currentModel =
374 (*fForwardLossProcesses[i])[k]->
375
376 G4double chargeSqRatio = currentMo
377 fFwdIon, couple->GetMaterial(),
378 (*fForwardLossProcesses[i])[k]->Se
379
380 }
381 G4double e1 = e / fMassRatio;
382 totCS += (*fForwardLossProcesses[i])
383 }
384 aVector->PutValue(l, totCS);
385 if(totCS > sigma_max)
386 {
387 sigma_max = totCS;
388 e_sigma_max = e;
389 }
390 if(totCS > 0 && !Emin_found)
391 {
392 fEminForFwdSigmaTables[i].push_back(
393 Emin_found = true;
394 }
395 }
396
397 fEkinofFwdSigmaMax[i].push_back(e_sigma_
398
399 if(!Emin_found)
400 fEminForFwdSigmaTables[i].push_back(fT
401
402 fTotalFwdSigmaTable[i]->push_back(aVecto
403
404 Emin_found = false;
405 sigma_max = 0;
406 e_sigma_max = 0.;
407 G4PhysicsVector* aVector1 = new G4Physic
408 for(std::size_t eindex = 0; eindex < fNb
409 {
410 G4double e = aVector1->Energy(eind
411 G4double totCS = ComputeTotalAdjointCS
412 couple, thePartDef,
413 e * 0.9999999 / fMassRatio); // fMa
414 aVector1->PutValue(eindex, totCS);
415 if(totCS > sigma_max)
416 {
417 sigma_max = totCS;
418 e_sigma_max = e;
419 }
420 if(totCS > 0 && !Emin_found)
421 {
422 fEminForAdjSigmaTables[i].push_back(
423 Emin_found = true;
424 }
425 }
426 fEkinofAdjSigmaMax[i].push_back(e_sigma_
427 if(!Emin_found)
428 fEminForAdjSigmaTables[i].push_back(fT
429
430 fTotalAdjSigmaTable[i]->push_back(aVecto
431 }
432 }
433 fSigmaTableBuilt = true;
434 }
435
436 //////////////////////////////////////////////
437 G4double G4AdjointCSManager::GetTotalAdjointCS
438 G4ParticleDefinition* aPartDef, G4double Eki
439 const G4MaterialCutsCouple* aCouple)
440 {
441 DefineCurrentMaterial(aCouple);
442 DefineCurrentParticle(aPartDef);
443 return (((*fTotalAdjSigmaTable[fCurrentParti
444 ->Value(Ekin * fMassRatio));
445 }
446
447 //////////////////////////////////////////////
448 G4double G4AdjointCSManager::GetTotalForwardCS
449 G4ParticleDefinition* aPartDef, G4double Eki
450 const G4MaterialCutsCouple* aCouple)
451 {
452 DefineCurrentMaterial(aCouple);
453 DefineCurrentParticle(aPartDef);
454 return (((*fTotalFwdSigmaTable[fCurrentParti
455 ->Value(Ekin * fMassRatio));
456 }
457
458 //////////////////////////////////////////////
459 G4double G4AdjointCSManager::GetAdjointSigma(
460 G4double Ekin_nuc, std::size_t index_model,
461 const G4MaterialCutsCouple* aCouple)
462 {
463 DefineCurrentMaterial(aCouple);
464 if(is_scat_proj_to_proj)
465 return (((*fSigmaTableForAdjointModelScatP
466 [fCurrentMatIndex])->Value(Ekin
467 else
468 return (
469 ((*fSigmaTableForAdjointModelProdToProj[
470 ->Value(Ekin_nuc));
471 }
472
473 //////////////////////////////////////////////
474 void G4AdjointCSManager::GetEminForTotalCS(G4P
475 con
476 G4d
477 G4d
478 {
479 DefineCurrentMaterial(aCouple);
480 DefineCurrentParticle(aPartDef);
481 emin_adj = fEminForAdjSigmaTables[fCurrentPa
482 fMassRatio;
483 emin_fwd = fEminForFwdSigmaTables[fCurrentPa
484 fMassRatio;
485 }
486
487 //////////////////////////////////////////////
488 void G4AdjointCSManager::GetMaxFwdTotalCS(G4Pa
489 cons
490 G4do
491 G4do
492 {
493 DefineCurrentMaterial(aCouple);
494 DefineCurrentParticle(aPartDef);
495 e_sigma_max = fEkinofFwdSigmaMax[fCurrentPar
496 sigma_max = ((*fTotalFwdSigmaTable[fCurrentP
497 ->Value(e_sigma_max);
498 e_sigma_max /= fMassRatio;
499 }
500
501 //////////////////////////////////////////////
502 void G4AdjointCSManager::GetMaxAdjTotalCS(G4Pa
503 cons
504 G4do
505 G4do
506 {
507 DefineCurrentMaterial(aCouple);
508 DefineCurrentParticle(aPartDef);
509 e_sigma_max = fEkinofAdjSigmaMax[fCurrentPar
510 sigma_max = ((*fTotalAdjSigmaTable[fCurrentP
511 ->Value(e_sigma_max);
512 e_sigma_max /= fMassRatio;
513 }
514
515 //////////////////////////////////////////////
516 G4double G4AdjointCSManager::GetCrossSectionCo
517 G4ParticleDefinition* aPartDef, G4double Pre
518 const G4MaterialCutsCouple* aCouple, G4bool&
519 {
520 static G4double lastEkin = 0.;
521 static G4ParticleDefinition* lastPartDef;
522
523 G4double corr_fac = 1.;
524 if(fForwardCSMode && aPartDef)
525 {
526 if(lastEkin != PreStepEkin || aPartDef !=
527 aCouple != fCurrentCouple)
528 {
529 DefineCurrentMaterial(aCouple);
530 G4double preadjCS = GetTotalAdjointCS(aP
531 G4double prefwdCS = GetTotalForwardCS(aP
532 lastEkin = PreStepEkin;
533 lastPartDef = aPartDef;
534 if(prefwdCS > 0. && preadjCS > 0.)
535 {
536 fForwardCSUsed = true;
537 fLastCSCorrectionFactor = prefwdCS / p
538 }
539 else
540 {
541 fForwardCSUsed = false;
542 fLastCSCorrectionFactor = 1.;
543 }
544 }
545 corr_fac = fLastCSCorrectionFactor;
546 }
547 else
548 {
549 fForwardCSUsed = false;
550 fLastCSCorrectionFactor = 1.;
551 }
552 fwd_is_used = fForwardCSUsed;
553 return corr_fac;
554 }
555
556 //////////////////////////////////////////////
557 G4double G4AdjointCSManager::GetContinuousWeig
558 G4ParticleDefinition* aPartDef, G4double Pre
559 const G4MaterialCutsCouple* aCouple, G4doubl
560 {
561 G4double corr_fac = 1.;
562 G4double after_fwdCS = GetTotalForwardCS(aPa
563 G4double pre_adjCS = GetTotalAdjointCS(aPa
564 if(!fForwardCSUsed || pre_adjCS == 0. || aft
565 {
566 G4double pre_fwdCS = GetTotalForwardCS(aPa
567 corr_fac *= std::exp((pre_adjCS - pre_fwdC
568 fLastCSCorrectionFactor = 1.;
569 }
570 else
571 {
572 fLastCSCorrectionFactor = after_fwdCS / pr
573 }
574 return corr_fac;
575 }
576
577 //////////////////////////////////////////////
578 G4double G4AdjointCSManager::GetPostStepWeight
579 {
580 return 1. / fLastCSCorrectionFactor;
581 }
582
583 //////////////////////////////////////////////
584 G4double G4AdjointCSManager::ComputeAdjointCS(
585 G4Material* aMaterial, G4VEmAdjointModel* aM
586 G4double Tcut, G4bool isScatProjToProj, std:
587 {
588 G4double EminSec = 0.;
589 G4double EmaxSec = 0.;
590
591 static G4double lastPrimaryEnergy = 0.;
592 static G4double lastTcut = 0.;
593 static G4Material* lastMaterial = nullptr;
594
595 if(isScatProjToProj)
596 {
597 EminSec = aModel->GetSecondAdjEnergyMinFor
598 EmaxSec = aModel->GetSecondAdjEnergyMaxFor
599 }
600 else if(PrimEnergy > Tcut || !aModel->GetApp
601 {
602 EminSec = aModel->GetSecondAdjEnergyMinFor
603 EmaxSec = aModel->GetSecondAdjEnergyMaxFor
604 }
605 if(EminSec >= EmaxSec)
606 return 0.;
607
608 G4bool need_to_compute = false;
609 if(aMaterial != lastMaterial || PrimEnergy !
610 Tcut != lastTcut)
611 {
612 lastMaterial = aMaterial;
613 lastPrimaryEnergy = PrimEnergy;
614 lastTcut = Tcut;
615 fIndexOfAdjointEMModelInAction.clear();
616 fIsScatProjToProj.clear();
617 fLastAdjointCSVsModelsAndElements.clear();
618 need_to_compute = true;
619 }
620
621 std::size_t ind = 0;
622 if(!need_to_compute)
623 {
624 need_to_compute = true;
625 for(std::size_t i = 0; i < fIndexOfAdjoint
626 {
627 std::size_t ind1 = fIndexOfAdjointEMMode
628 if(aModel == fAdjointModels[ind1] &&
629 isScatProjToProj == fIsScatProjToProj
630 {
631 need_to_compute = false;
632 CS_Vs_Element = fLastAdjointCSVsMode
633 }
634 ++ind;
635 }
636 }
637
638 if(need_to_compute)
639 {
640 std::size_t ind_model = 0;
641 for(std::size_t i = 0; i < fAdjointModels.
642 {
643 if(aModel == fAdjointModels[i])
644 {
645 ind_model = i;
646 break;
647 }
648 }
649 G4double Tlow = Tcut;
650 if(!fAdjointModels[ind_model]->GetApplyCut
651 Tlow = fAdjointModels[ind_model]->GetLow
652 fIndexOfAdjointEMModelInAction.push_back(i
653 fIsScatProjToProj.push_back(isScatProjToPr
654 CS_Vs_Element.clear();
655 if(!aModel->GetUseMatrix())
656 {
657 CS_Vs_Element.push_back(aModel->AdjointC
658 fCurrentCouple, PrimEnergy, isScatProj
659 }
660 else if(aModel->GetUseMatrixPerElement())
661 {
662 std::size_t n_el = aMaterial->GetNumberO
663 if(aModel->GetUseOnlyOneMatrixForAllElem
664 {
665 G4AdjointCSMatrix* theCSMatrix;
666 if(isScatProjToProj)
667 {
668 theCSMatrix = fAdjointCSMatricesForS
669 }
670 else
671 theCSMatrix = fAdjointCSMatricesForP
672 G4double CS = 0.;
673 if(PrimEnergy > Tlow)
674 CS = ComputeAdjointCS(PrimEnergy, th
675 G4double factor = 0.;
676 for(G4int i = 0; i < (G4int)n_el; ++i)
677 { // this could be computed only once
678 factor += aMaterial->GetElement(i)->
679 aMaterial->GetVecNbOfAtoms
680 }
681 CS *= factor;
682 CS_Vs_Element.push_back(CS);
683 }
684 else
685 {
686 for(G4int i = 0; i < (G4int)n_el; ++i)
687 {
688 std::size_t ind_el = aMaterial->GetE
689 G4AdjointCSMatrix* theCSMatrix;
690 if(isScatProjToProj)
691 {
692 theCSMatrix =
693 fAdjointCSMatricesForScatProjToP
694 }
695 else
696 theCSMatrix = fAdjointCSMatricesFo
697 G4double CS = 0.;
698 if(PrimEnergy > Tlow)
699 CS = ComputeAdjointCS(PrimEnergy,
700 CS_Vs_Element.push_back(CS *
701 (aMaterial->
702 }
703 }
704 }
705 else
706 {
707 std::size_t ind_mat = aMaterial->GetInde
708 G4AdjointCSMatrix* theCSMatrix;
709 if(isScatProjToProj)
710 {
711 theCSMatrix = fAdjointCSMatricesForSca
712 }
713 else
714 theCSMatrix = fAdjointCSMatricesForPro
715 G4double CS = 0.;
716 if(PrimEnergy > Tlow)
717 CS = ComputeAdjointCS(PrimEnergy, theC
718 CS_Vs_Element.push_back(CS);
719 }
720 fLastAdjointCSVsModelsAndElements.push_bac
721 }
722
723 G4double CS = 0.;
724 for(const auto& cs_vs_el : CS_Vs_Element)
725 {
726 // We could put the progressive sum of the
727 // element itself
728 CS += cs_vs_el;
729 }
730 return CS;
731 }
732
733 //////////////////////////////////////////////
734 G4Element* G4AdjointCSManager::SampleElementFr
735 G4Material* aMaterial, G4VEmAdjointModel* aM
736 G4double Tcut, G4bool isScatProjToProj)
737 {
738 std::vector<G4double> CS_Vs_Element;
739 G4double CS = ComputeAdjointCS(aMaterial,
740 isScatProjToP
741 G4double SumCS = 0.;
742 std::size_t ind = 0;
743 for(std::size_t i = 0; i < CS_Vs_Element.siz
744 {
745 SumCS += CS_Vs_Element[i];
746 if(G4UniformRand() <= SumCS / CS)
747 {
748 ind = i;
749 break;
750 }
751 }
752
753 return const_cast<G4Element*>(aMaterial->Get
754 }
755
756 //////////////////////////////////////////////
757 G4double G4AdjointCSManager::ComputeTotalAdjoi
758 const G4MaterialCutsCouple* aCouple, G4Parti
759 G4double Ekin)
760 {
761 G4double TotalCS = 0.;
762
763 DefineCurrentMaterial(aCouple);
764
765 std::vector<G4double> CS_Vs_Element;
766 G4double CS;
767 G4VEmAdjointModel* adjModel = nullptr;
768 for(std::size_t i = 0; i < fAdjointModels.si
769 {
770 G4double Tlow = 0.;
771 adjModel = fAdjointModels[i];
772 if(!adjModel->GetApplyCutInRange())
773 Tlow = adjModel->GetLowEnergyLimit();
774 else
775 {
776 G4ParticleDefinition* theDirSecondPartDe
777 adjModel->GetAdjointEquivalentOfDirect
778 std::size_t idx = 56;
779 if(theDirSecondPartDef->GetParticleName(
780 idx = 0;
781 else if(theDirSecondPartDef->GetParticle
782 idx = 1;
783 else if(theDirSecondPartDef->GetParticle
784 idx = 2;
785 if(idx < 56)
786 {
787 const std::vector<G4double>* aVec =
788 G4ProductionCutsTable::GetProduction
789 idx);
790 Tlow = (*aVec)[aCouple->GetIndex()];
791 }
792 }
793 if(Ekin <= adjModel->GetHighEnergyLimit()
794 Ekin >= adjModel->GetLowEnergyLimit())
795 {
796 if(aPartDef ==
797 adjModel->GetAdjointEquivalentOfDirec
798 {
799 CS = ComputeAdjointCS(fCurrentMaterial
800 CS_Vs_Element);
801 TotalCS += CS;
802 (*fSigmaTableForAdjointModelScatProjTo
803 ->PutValue(fNbins, CS);
804 }
805 if(aPartDef ==
806 adjModel->GetAdjointEquivalentOfDirec
807 {
808 CS = ComputeAdjointCS(fCurrentMaterial
809 CS_Vs_Element);
810 TotalCS += CS;
811 (*fSigmaTableForAdjointModelProdToProj
812 fNbins, CS);
813 }
814 }
815 else
816 {
817 (*fSigmaTableForAdjointModelScatProjToPr
818 ->PutValue(fNbins, 0.);
819 (*fSigmaTableForAdjointModelProdToProj[i
820 fNbins, 0.);
821 }
822 }
823 return TotalCS;
824 }
825
826 //////////////////////////////////////////////
827 std::vector<G4AdjointCSMatrix*>
828 G4AdjointCSManager::BuildCrossSectionsModelAnd
829
830
831 {
832 G4AdjointCSMatrix* theCSMatForProdToProjBack
833 new G4AdjointCSMatrix(false);
834 G4AdjointCSMatrix* theCSMatForScatProjToProj
835 new G4AdjointCSMatrix(true);
836
837 // make the vector of primary energy of the
838 G4double EkinMin = aModel->GetLowEner
839 G4double EkinMaxForScat = aModel->GetHighEne
840 G4double EkinMaxForProd = aModel->GetHighEne
841 if(aModel->GetSecondPartOfSameType())
842 EkinMaxForProd = EkinMaxForProd / 2.;
843
844 // Product to projectile backward scattering
845 G4double dE = std::pow(10., 1. / nbin_pro_de
846 G4double E2 =
847 std::pow(10., double(int(std::log10(EkinMi
848 nbin_pro_decade) / dE;
849 G4double E1 = EkinMin;
850 // Loop checking, 07-Aug-2015, Vladimir Ivan
851 while(E1 < EkinMaxForProd)
852 {
853 E1 = std::max(EkinMin, E2);
854 E1 = std::min(EkinMaxForProd, E1);
855 std::vector<std::vector<double>*> aMat =
856 aModel->ComputeAdjointCrossSectionVector
857
858 if(aMat.size() >= 2)
859 {
860 std::vector<double>* log_ESecVec = aMat[
861 std::vector<double>* log_CSVec = aMat[
862 G4double log_adjointCS = log_C
863 // normalise CSVec such that it becomes
864 for(std::size_t j = 0; j < log_CSVec->si
865 {
866 if(j == 0)
867 (*log_CSVec)[j] = 0.;
868 else
869 (*log_CSVec)[j] =
870 std::log(1. - std::exp((*log_CSVec
871 }
872 (*log_CSVec)[log_CSVec->size() - 1] =
873 (*log_CSVec)[log_CSVec->size() - 2] -
874 theCSMatForProdToProjBackwardScattering-
875 std::log(E1), log_adjointCS, log_ESecV
876 }
877 E1 = E2;
878 E2 *= dE;
879 }
880
881 // Scattered projectile to projectile backwa
882 E2 = std::pow(10., G4double(int(std::log10(E
883 nbin_pro_decade) / dE;
884 E1 = EkinMin;
885 // Loop checking, 07-Aug-2015, Vladimir Ivan
886 while(E1 < EkinMaxForScat)
887 {
888 E1 = std::max(EkinMin, E2);
889 E1 = std::min(EkinMaxForScat, E1);
890 std::vector<std::vector<G4double>*> aMat =
891 aModel->ComputeAdjointCrossSectionVector
892 E1, Z, A, nbin_pro_decade);
893 if(aMat.size() >= 2)
894 {
895 std::vector<G4double>* log_ESecVec = aMa
896 std::vector<G4double>* log_CSVec = aMa
897 G4double log_adjointCS = log
898 // normalise CSVec such that it becomes
899 for(std::size_t j = 0; j < log_CSVec->si
900 {
901 if(j == 0)
902 (*log_CSVec)[j] = 0.;
903 else
904 (*log_CSVec)[j] =
905 std::log(1. - std::exp((*log_CSVec
906 }
907 (*log_CSVec)[log_CSVec->size() - 1] =
908 (*log_CSVec)[log_CSVec->size() - 2] -
909 theCSMatForScatProjToProjBackwardScatter
910 std::log(E1), log_adjointCS, log_ESecV
911 }
912 E1 = E2;
913 E2 *= dE;
914 }
915
916 std::vector<G4AdjointCSMatrix*> res;
917 res.push_back(theCSMatForProdToProjBackwardS
918 res.push_back(theCSMatForScatProjToProjBackw
919
920 return res;
921 }
922
923 //////////////////////////////////////////////
924 std::vector<G4AdjointCSMatrix*>
925 G4AdjointCSManager::BuildCrossSectionsModelAnd
926 G4VEmAdjointModel* aModel, G4Material* aMate
927 {
928 G4AdjointCSMatrix* theCSMatForProdToProjBack
929 new G4AdjointCSMatrix(false);
930 G4AdjointCSMatrix* theCSMatForScatProjToProj
931 new G4AdjointCSMatrix(true);
932
933 G4double EkinMin = aModel->GetLowEner
934 G4double EkinMaxForScat = aModel->GetHighEne
935 G4double EkinMaxForProd = aModel->GetHighEne
936 if(aModel->GetSecondPartOfSameType())
937 EkinMaxForProd /= 2.;
938
939 // Product to projectile backward scattering
940 G4double dE = std::pow(10., 1. / nbin_pro_de
941 G4double E2 =
942 std::pow(10., double(int(std::log10(EkinMi
943 nbin_pro_decade) / dE;
944 G4double E1 = EkinMin;
945 // Loop checking, 07-Aug-2015, Vladimir Ivan
946 while(E1 < EkinMaxForProd)
947 {
948 E1 = std::max(EkinMin, E2);
949 E1 = std::min(EkinMaxForProd, E1);
950 std::vector<std::vector<G4double>*> aMat =
951 aModel->ComputeAdjointCrossSectionVector
952 aMaterial, E1, nbin_pro_decade);
953 if(aMat.size() >= 2)
954 {
955 std::vector<G4double>* log_ESecVec = aMa
956 std::vector<G4double>* log_CSVec = aMa
957 G4double log_adjointCS = log
958
959 // normalise CSVec such that it becomes
960 for(std::size_t j = 0; j < log_CSVec->si
961 {
962 if(j == 0)
963 (*log_CSVec)[j] = 0.;
964 else
965 (*log_CSVec)[j] =
966 std::log(1. - std::exp((*log_CSVec
967 }
968 (*log_CSVec)[log_CSVec->size() - 1] =
969 (*log_CSVec)[log_CSVec->size() - 2] -
970 theCSMatForProdToProjBackwardScattering-
971 std::log(E1), log_adjointCS, log_ESecV
972 }
973
974 E1 = E2;
975 E2 *= dE;
976 }
977
978 // Scattered projectile to projectile backwa
979 E2 =
980 std::pow(10., G4double(G4int(std::log10(Ek
981 nbin_pro_decade) /
982 dE;
983 E1 = EkinMin;
984 while(E1 < EkinMaxForScat)
985 {
986 E1 = std::max(EkinMin, E2);
987 E1 = std::min(EkinMaxForScat, E1);
988 std::vector<std::vector<G4double>*> aMat =
989 aModel->ComputeAdjointCrossSectionVector
990 aMaterial, E1, nbin_pro_decade);
991 if(aMat.size() >= 2)
992 {
993 std::vector<G4double>* log_ESecVec = aMa
994 std::vector<G4double>* log_CSVec = aMa
995 G4double log_adjointCS = log
996
997 for(std::size_t j = 0; j < log_CSVec->si
998 {
999 if(j == 0)
1000 (*log_CSVec)[j] = 0.;
1001 else
1002 (*log_CSVec)[j] =
1003 std::log(1. - std::exp((*log_CSVe
1004 }
1005 (*log_CSVec)[log_CSVec->size() - 1] =
1006 (*log_CSVec)[log_CSVec->size() - 2] -
1007
1008 theCSMatForScatProjToProjBackwardScatte
1009 std::log(E1), log_adjointCS, log_ESec
1010 }
1011 E1 = E2;
1012 E2 *= dE;
1013 }
1014
1015 std::vector<G4AdjointCSMatrix*> res;
1016 res.push_back(theCSMatForProdToProjBackward
1017 res.push_back(theCSMatForScatProjToProjBack
1018
1019 return res;
1020 }
1021
1022 /////////////////////////////////////////////
1023 G4ParticleDefinition* G4AdjointCSManager::Get
1024 G4ParticleDefinition* theFwdPartDef)
1025 {
1026 if(theFwdPartDef->GetParticleName() == "e-"
1027 return G4AdjointElectron::AdjointElectron
1028 else if(theFwdPartDef->GetParticleName() ==
1029 return G4AdjointGamma::AdjointGamma();
1030 else if(theFwdPartDef->GetParticleName() ==
1031 return G4AdjointProton::AdjointProton();
1032 else if(theFwdPartDef == fFwdIon)
1033 return fAdjIon;
1034 return nullptr;
1035 }
1036
1037 /////////////////////////////////////////////
1038 G4ParticleDefinition* G4AdjointCSManager::Get
1039 G4ParticleDefinition* theAdjPartDef)
1040 {
1041 if(theAdjPartDef->GetParticleName() == "adj
1042 return G4Electron::Electron();
1043 else if(theAdjPartDef->GetParticleName() ==
1044 return G4Gamma::Gamma();
1045 else if(theAdjPartDef->GetParticleName() ==
1046 return G4Proton::Proton();
1047 else if(theAdjPartDef == fAdjIon)
1048 return fFwdIon;
1049 return nullptr;
1050 }
1051
1052 /////////////////////////////////////////////
1053 void G4AdjointCSManager::DefineCurrentMateria
1054 const G4MaterialCutsCouple* couple)
1055 {
1056 if(couple != fCurrentCouple)
1057 {
1058 fCurrentCouple = const_cast<G4Ma
1059 fCurrentMaterial = const_cast<G4Ma
1060 fCurrentMatIndex = couple->GetInde
1061 fLastCSCorrectionFactor = 1.;
1062 }
1063 }
1064
1065 /////////////////////////////////////////////
1066 void G4AdjointCSManager::DefineCurrentParticl
1067 const G4ParticleDefinition* aPartDef)
1068 {
1069
1070 if(aPartDef != fCurrentParticleDef)
1071 {
1072 fCurrentParticleDef = const_cast<G4Partic
1073 fMassRatio = 1.;
1074 if(aPartDef == fAdjIon)
1075 fMassRatio = proton_mass_c2 / aPartDef-
1076 fCurrentParticleIndex = 1000000;
1077 for(std::size_t i = 0; i < fAdjointPartic
1078 {
1079 if(aPartDef == fAdjointParticlesInActio
1080 fCurrentParticleIndex = i;
1081 }
1082 }
1083 }
1084
1085 /////////////////////////////////////////////
1086 G4double G4AdjointCSManager::ComputeAdjointCS
1087 G4double aPrimEnergy, G4AdjointCSMatrix* an
1088 {
1089 std::vector<double>* theLogPrimEnergyVector
1090 anAdjointCSMatrix->GetLogPrimEnergyVector
1091 if(theLogPrimEnergyVector->empty())
1092 {
1093 G4cout << "No data are contained in the g
1094 return 0.;
1095 }
1096 G4double log_Tcut = std::log(Tcut);
1097 G4double log_E = std::log(aPrimEnergy);
1098
1099 if(aPrimEnergy <= Tcut || log_E > theLogPri
1100 return 0.;
1101
1102 G4AdjointInterpolator* theInterpolator = G4
1103
1104 std::size_t ind =
1105 theInterpolator->FindPositionForLogVector
1106 G4double aLogPrimEnergy1, aLogPrimEnergy2;
1107 G4double aLogCS1, aLogCS2;
1108 G4double log01, log02;
1109 std::vector<G4double>* aLogSecondEnergyVect
1110 std::vector<G4double>* aLogSecondEnergyVect
1111 std::vector<G4double>* aLogProbVector1
1112 std::vector<G4double>* aLogProbVector2
1113 std::vector<std::size_t>* aLogProbVectorInd
1114 std::vector<std::size_t>* aLogProbVectorInd
1115
1116 anAdjointCSMatrix->GetData((G4int)ind, aLog
1117 aLogSecondEnergy
1118 aLogProbVectorIn
1119 anAdjointCSMatrix->GetData(G4int(ind + 1),
1120 aLogSecondEnergy
1121 aLogProbVectorIn
1122 if (! (aLogProbVector1 && aLogProbVector2 &
1123 aLogSecondEnergyVector1 && aLogS
1124 return 0.;
1125 }
1126
1127 if(anAdjointCSMatrix->IsScatProjToProj())
1128 { // case where the Tcut plays a role
1129 G4double log_minimum_prob1, log_minimum_p
1130 log_minimum_prob1 = theInterpolator->Inte
1131 log_Tcut, *aLogSecondEnergyVector1, *aL
1132 log_minimum_prob2 = theInterpolator->Inte
1133 log_Tcut, *aLogSecondEnergyVector2, *aL
1134 aLogCS1 += log_minimum_prob1;
1135 aLogCS2 += log_minimum_prob2;
1136 }
1137
1138 G4double log_adjointCS = theInterpolator->L
1139 log_E, aLogPrimEnergy1, aLogPrimEnergy2,
1140 return std::exp(log_adjointCS);
1141 }
1142