| 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 //
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* G4AdjointCSManager::fInstance = nullptr;
51
52 constexpr G4double G4AdjointCSManager::fTmin;
53 constexpr G4double G4AdjointCSManager::fTmax;
54 constexpr G4int G4AdjointCSManager::fNbins;
55
56 ///////////////////////////////////////////////////////
57 G4AdjointCSManager* G4AdjointCSManager::GetAdjointCSManager()
58 {
59 if(fInstance == nullptr)
60 {
61 static G4ThreadLocalSingleton<G4AdjointCSManager> inst;
62 fInstance = inst.Instance();
63 }
64 return fInstance;
65 }
66
67 ///////////////////////////////////////////////////////
68 G4AdjointCSManager::G4AdjointCSManager()
69 {
70 RegisterAdjointParticle(G4AdjointElectron::AdjointElectron());
71 RegisterAdjointParticle(G4AdjointGamma::AdjointGamma());
72 RegisterAdjointParticle(G4AdjointProton::AdjointProton());
73 }
74
75 ///////////////////////////////////////////////////////
76 G4AdjointCSManager::~G4AdjointCSManager()
77 {
78 for (auto& p : fAdjointCSMatricesForProdToProj) {
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 : fAdjointCSMatricesForScatProjToProj) {
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 : fSigmaTableForAdjointModelScatProjToProj) {
116 p->clearAndDestroy();
117 delete p;
118 p = nullptr;
119 }
120 fSigmaTableForAdjointModelScatProjToProj.clear();
121
122 for (auto p : fSigmaTableForAdjointModelProdToProj) {
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::RegisterEmAdjointModel(G4VEmAdjointModel* aModel)
151 {
152 fAdjointModels.push_back(aModel);
153 fSigmaTableForAdjointModelScatProjToProj.push_back(new G4PhysicsTable);
154 fSigmaTableForAdjointModelProdToProj.push_back(new G4PhysicsTable);
155 return fAdjointModels.size() - 1;
156 }
157
158 ///////////////////////////////////////////////////////
159 void G4AdjointCSManager::RegisterEmProcess(G4VEmProcess* aProcess,
160 G4ParticleDefinition* aFwdPartDef)
161 {
162 G4ParticleDefinition* anAdjPartDef =
163 GetAdjointParticleEquivalent(aFwdPartDef);
164 if(anAdjPartDef && aProcess)
165 {
166 RegisterAdjointParticle(anAdjPartDef);
167
168 for(std::size_t i = 0; i < fAdjointParticlesInAction.size(); ++i)
169 {
170 if(anAdjPartDef->GetParticleName() ==
171 fAdjointParticlesInAction[i]->GetParticleName())
172 fForwardProcesses[i]->push_back(aProcess);
173 }
174 }
175 }
176
177 ///////////////////////////////////////////////////////
178 void G4AdjointCSManager::RegisterEnergyLossProcess(
179 G4VEnergyLossProcess* aProcess, G4ParticleDefinition* aFwdPartDef)
180 {
181 G4ParticleDefinition* anAdjPartDef =
182 GetAdjointParticleEquivalent(aFwdPartDef);
183 if(anAdjPartDef && aProcess)
184 {
185 RegisterAdjointParticle(anAdjPartDef);
186 for(std::size_t i = 0; i < fAdjointParticlesInAction.size(); ++i)
187 {
188 if(anAdjPartDef->GetParticleName() ==
189 fAdjointParticlesInAction[i]->GetParticleName())
190 fForwardLossProcesses[i]->push_back(aProcess);
191
192 }
193 }
194 }
195
196 ///////////////////////////////////////////////////////
197 void G4AdjointCSManager::RegisterAdjointParticle(G4ParticleDefinition* aPartDef)
198 {
199 G4bool found = false;
200 for(auto p : fAdjointParticlesInAction)
201 {
202 if(p->GetParticleName() == aPartDef->GetParticleName())
203 {
204 found = true;
205 }
206 }
207 if(!found)
208 {
209 fForwardLossProcesses.push_back(new std::vector<G4VEnergyLossProcess*>());
210 fTotalFwdSigmaTable.push_back(new G4PhysicsTable);
211 fTotalAdjSigmaTable.push_back(new G4PhysicsTable);
212 fForwardProcesses.push_back(new std::vector<G4VEmProcess*>());
213 fAdjointParticlesInAction.push_back(aPartDef);
214 fEminForFwdSigmaTables.push_back(std::vector<G4double>());
215 fEminForAdjSigmaTables.push_back(std::vector<G4double>());
216 fEkinofFwdSigmaMax.push_back(std::vector<G4double>());
217 fEkinofAdjSigmaMax.push_back(std::vector<G4double>());
218 }
219 }
220
221 ///////////////////////////////////////////////////////
222 void G4AdjointCSManager::BuildCrossSectionMatrices()
223 {
224 if(fCSMatricesBuilt)
225 return;
226 // The Tcut and Tmax matrices will be computed probably just once.
227 // When Tcut changes, some PhysicsTable will be recomputed
228 // for each MaterialCutCouple but not all the matrices.
229 // The Tcut defines a lower limit in the energy of the projectile before
230 // scattering. In the Projectile to Scattered Projectile case we have
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 adjoint CS we should integrate over
235 // Epro from E_ScatProj+Tcut to Emax
236 // In the Projectile to Secondary case Tcut plays a role only in the fact that
237 // Esecond should be greater than Tcut to have the possibility to have any
238 // adjoint process.
239 // To avoid recomputing matrices for all changes of MaterialCutCouple,
240 // we propose to compute the matrices only once for the minimum possible Tcut
241 // and then to interpolate the probability for a new Tcut (implemented in
242 // G4VAdjointEmModel)
243
244 fAdjointCSMatricesForScatProjToProj.clear();
245 fAdjointCSMatricesForProdToProj.clear();
246 const G4ElementTable* theElementTable = G4Element::GetElementTable();
247 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
248
249 G4cout << "========== Computation of cross section matrices for adjoint "
250 "models =========="
251 << G4endl;
252 for(const auto& aModel : fAdjointModels)
253 {
254 G4cout << "Build adjoint cross section matrices for " << aModel->GetName()
255 << G4endl;
256 if(aModel->GetUseMatrix())
257 {
258 std::vector<G4AdjointCSMatrix*>* aListOfMat1 =
259 new std::vector<G4AdjointCSMatrix*>();
260 std::vector<G4AdjointCSMatrix*>* aListOfMat2 =
261 new std::vector<G4AdjointCSMatrix*>();
262 if(aModel->GetUseMatrixPerElement())
263 {
264 if(aModel->GetUseOnlyOneMatrixForAllElements())
265 {
266 std::vector<G4AdjointCSMatrix*> two_matrices =
267 BuildCrossSectionsModelAndElement(aModel, 1, 1, 80);
268 aListOfMat1->push_back(two_matrices[0]);
269 aListOfMat2->push_back(two_matrices[1]);
270 }
271 else
272 {
273 for(const auto& anElement : *theElementTable)
274 {
275 G4int Z = G4lrint(anElement->GetZ());
276 G4int A = G4lrint(anElement->GetN());
277 std::vector<G4AdjointCSMatrix*> two_matrices =
278 BuildCrossSectionsModelAndElement(aModel, Z, A, 40);
279 aListOfMat1->push_back(two_matrices[0]);
280 aListOfMat2->push_back(two_matrices[1]);
281 }
282 }
283 }
284 else
285 { // Per material case
286 for(const auto& aMaterial : *theMaterialTable)
287 {
288 std::vector<G4AdjointCSMatrix*> two_matrices =
289 BuildCrossSectionsModelAndMaterial(aModel, aMaterial, 40);
290 aListOfMat1->push_back(two_matrices[0]);
291 aListOfMat2->push_back(two_matrices[1]);
292 }
293 }
294 fAdjointCSMatricesForProdToProj.push_back(*aListOfMat1);
295 fAdjointCSMatricesForScatProjToProj.push_back(*aListOfMat2);
296 aModel->SetCSMatrices(aListOfMat1, aListOfMat2);
297 }
298 else
299 {
300 G4cout << "The model " << aModel->GetName()
301 << " does not use cross section matrices" << G4endl;
302 std::vector<G4AdjointCSMatrix*> two_empty_matrices;
303 fAdjointCSMatricesForProdToProj.push_back(two_empty_matrices);
304 fAdjointCSMatricesForScatProjToProj.push_back(std::move(two_empty_matrices));
305 }
306 }
307 G4cout << " All adjoint cross section matrices are computed!"
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::GetProductionCutsTable();
324
325 // Prepare the Sigma table for all AdjointEMModel, will be filled later on
326 for(std::size_t i = 0; i < fAdjointModels.size(); ++i)
327 {
328 fSigmaTableForAdjointModelScatProjToProj[i]->clearAndDestroy();
329 fSigmaTableForAdjointModelProdToProj[i]->clearAndDestroy();
330 for(std::size_t j = 0; j < theCoupleTable->GetTableSize(); ++j)
331 {
332 fSigmaTableForAdjointModelScatProjToProj[i]->push_back(
333 new G4PhysicsLogVector(fTmin, fTmax, fNbins));
334 fSigmaTableForAdjointModelProdToProj[i]->push_back(
335 new G4PhysicsLogVector(fTmin, fTmax, fNbins));
336 }
337 }
338
339 for(std::size_t i = 0; i < fAdjointParticlesInAction.size(); ++i)
340 {
341 G4ParticleDefinition* thePartDef = fAdjointParticlesInAction[i];
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->GetTableSize(); ++j)
351 {
352 const G4MaterialCutsCouple* couple =
353 theCoupleTable->GetMaterialCutsCouple((G4int)j);
354
355 // make first the total fwd CS table for FwdProcess
356 G4PhysicsVector* aVector = new G4PhysicsLogVector(fTmin, fTmax, fNbins);
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 < fForwardProcesses[i]->size(); ++k)
365 {
366 totCS += (*fForwardProcesses[i])[k]->GetCrossSection(e, couple);
367 }
368 for(std::size_t k = 0; k < fForwardLossProcesses[i]->size(); ++k)
369 {
370 if(thePartDef == fAdjIon)
371 { // e is considered already as the scaled energy
372 std::size_t mat_index = couple->GetIndex();
373 G4VEmModel* currentModel =
374 (*fForwardLossProcesses[i])[k]->SelectModelForMaterial(e,
375 mat_index);
376 G4double chargeSqRatio = currentModel->GetChargeSquareRatio(
377 fFwdIon, couple->GetMaterial(), e / fMassRatio);
378 (*fForwardLossProcesses[i])[k]->SetDynamicMassCharge(fMassRatio,
379 chargeSqRatio);
380 }
381 G4double e1 = e / fMassRatio;
382 totCS += (*fForwardLossProcesses[i])[k]->GetLambda(e1, couple);
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(e);
393 Emin_found = true;
394 }
395 }
396
397 fEkinofFwdSigmaMax[i].push_back(e_sigma_max);
398
399 if(!Emin_found)
400 fEminForFwdSigmaTables[i].push_back(fTmax);
401
402 fTotalFwdSigmaTable[i]->push_back(aVector);
403
404 Emin_found = false;
405 sigma_max = 0;
406 e_sigma_max = 0.;
407 G4PhysicsVector* aVector1 = new G4PhysicsLogVector(fTmin, fTmax, fNbins);
408 for(std::size_t eindex = 0; eindex < fNbins; ++eindex)
409 {
410 G4double e = aVector1->Energy(eindex);
411 G4double totCS = ComputeTotalAdjointCS(
412 couple, thePartDef,
413 e * 0.9999999 / fMassRatio); // fMassRatio needed for ions
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(e);
423 Emin_found = true;
424 }
425 }
426 fEkinofAdjSigmaMax[i].push_back(e_sigma_max);
427 if(!Emin_found)
428 fEminForAdjSigmaTables[i].push_back(fTmax);
429
430 fTotalAdjSigmaTable[i]->push_back(aVector1);
431 }
432 }
433 fSigmaTableBuilt = true;
434 }
435
436 ///////////////////////////////////////////////////////
437 G4double G4AdjointCSManager::GetTotalAdjointCS(
438 G4ParticleDefinition* aPartDef, G4double Ekin,
439 const G4MaterialCutsCouple* aCouple)
440 {
441 DefineCurrentMaterial(aCouple);
442 DefineCurrentParticle(aPartDef);
443 return (((*fTotalAdjSigmaTable[fCurrentParticleIndex])[fCurrentMatIndex])
444 ->Value(Ekin * fMassRatio));
445 }
446
447 ///////////////////////////////////////////////////////
448 G4double G4AdjointCSManager::GetTotalForwardCS(
449 G4ParticleDefinition* aPartDef, G4double Ekin,
450 const G4MaterialCutsCouple* aCouple)
451 {
452 DefineCurrentMaterial(aCouple);
453 DefineCurrentParticle(aPartDef);
454 return (((*fTotalFwdSigmaTable[fCurrentParticleIndex])[fCurrentMatIndex])
455 ->Value(Ekin * fMassRatio));
456 }
457
458 ///////////////////////////////////////////////////////
459 G4double G4AdjointCSManager::GetAdjointSigma(
460 G4double Ekin_nuc, std::size_t index_model, G4bool is_scat_proj_to_proj,
461 const G4MaterialCutsCouple* aCouple)
462 {
463 DefineCurrentMaterial(aCouple);
464 if(is_scat_proj_to_proj)
465 return (((*fSigmaTableForAdjointModelScatProjToProj[index_model])
466 [fCurrentMatIndex])->Value(Ekin_nuc));
467 else
468 return (
469 ((*fSigmaTableForAdjointModelProdToProj[index_model])[fCurrentMatIndex])
470 ->Value(Ekin_nuc));
471 }
472
473 ///////////////////////////////////////////////////////
474 void G4AdjointCSManager::GetEminForTotalCS(G4ParticleDefinition* aPartDef,
475 const G4MaterialCutsCouple* aCouple,
476 G4double& emin_adj,
477 G4double& emin_fwd)
478 {
479 DefineCurrentMaterial(aCouple);
480 DefineCurrentParticle(aPartDef);
481 emin_adj = fEminForAdjSigmaTables[fCurrentParticleIndex][fCurrentMatIndex] /
482 fMassRatio;
483 emin_fwd = fEminForFwdSigmaTables[fCurrentParticleIndex][fCurrentMatIndex] /
484 fMassRatio;
485 }
486
487 ///////////////////////////////////////////////////////
488 void G4AdjointCSManager::GetMaxFwdTotalCS(G4ParticleDefinition* aPartDef,
489 const G4MaterialCutsCouple* aCouple,
490 G4double& e_sigma_max,
491 G4double& sigma_max)
492 {
493 DefineCurrentMaterial(aCouple);
494 DefineCurrentParticle(aPartDef);
495 e_sigma_max = fEkinofFwdSigmaMax[fCurrentParticleIndex][fCurrentMatIndex];
496 sigma_max = ((*fTotalFwdSigmaTable[fCurrentParticleIndex])[fCurrentMatIndex])
497 ->Value(e_sigma_max);
498 e_sigma_max /= fMassRatio;
499 }
500
501 ///////////////////////////////////////////////////////
502 void G4AdjointCSManager::GetMaxAdjTotalCS(G4ParticleDefinition* aPartDef,
503 const G4MaterialCutsCouple* aCouple,
504 G4double& e_sigma_max,
505 G4double& sigma_max)
506 {
507 DefineCurrentMaterial(aCouple);
508 DefineCurrentParticle(aPartDef);
509 e_sigma_max = fEkinofAdjSigmaMax[fCurrentParticleIndex][fCurrentMatIndex];
510 sigma_max = ((*fTotalAdjSigmaTable[fCurrentParticleIndex])[fCurrentMatIndex])
511 ->Value(e_sigma_max);
512 e_sigma_max /= fMassRatio;
513 }
514
515 ///////////////////////////////////////////////////////
516 G4double G4AdjointCSManager::GetCrossSectionCorrection(
517 G4ParticleDefinition* aPartDef, G4double PreStepEkin,
518 const G4MaterialCutsCouple* aCouple, G4bool& fwd_is_used)
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 != lastPartDef ||
527 aCouple != fCurrentCouple)
528 {
529 DefineCurrentMaterial(aCouple);
530 G4double preadjCS = GetTotalAdjointCS(aPartDef, PreStepEkin, aCouple);
531 G4double prefwdCS = GetTotalForwardCS(aPartDef, PreStepEkin, aCouple);
532 lastEkin = PreStepEkin;
533 lastPartDef = aPartDef;
534 if(prefwdCS > 0. && preadjCS > 0.)
535 {
536 fForwardCSUsed = true;
537 fLastCSCorrectionFactor = prefwdCS / preadjCS;
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::GetContinuousWeightCorrection(
558 G4ParticleDefinition* aPartDef, G4double PreStepEkin, G4double AfterStepEkin,
559 const G4MaterialCutsCouple* aCouple, G4double step_length)
560 {
561 G4double corr_fac = 1.;
562 G4double after_fwdCS = GetTotalForwardCS(aPartDef, AfterStepEkin, aCouple);
563 G4double pre_adjCS = GetTotalAdjointCS(aPartDef, PreStepEkin, aCouple);
564 if(!fForwardCSUsed || pre_adjCS == 0. || after_fwdCS == 0.)
565 {
566 G4double pre_fwdCS = GetTotalForwardCS(aPartDef, PreStepEkin, aCouple);
567 corr_fac *= std::exp((pre_adjCS - pre_fwdCS) * step_length);
568 fLastCSCorrectionFactor = 1.;
569 }
570 else
571 {
572 fLastCSCorrectionFactor = after_fwdCS / pre_adjCS;
573 }
574 return corr_fac;
575 }
576
577 ///////////////////////////////////////////////////////
578 G4double G4AdjointCSManager::GetPostStepWeightCorrection()
579 {
580 return 1. / fLastCSCorrectionFactor;
581 }
582
583 ///////////////////////////////////////////////////////
584 G4double G4AdjointCSManager::ComputeAdjointCS(
585 G4Material* aMaterial, G4VEmAdjointModel* aModel, G4double PrimEnergy,
586 G4double Tcut, G4bool isScatProjToProj, std::vector<G4double>& CS_Vs_Element)
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->GetSecondAdjEnergyMinForScatProjToProj(PrimEnergy, Tcut);
598 EmaxSec = aModel->GetSecondAdjEnergyMaxForScatProjToProj(PrimEnergy);
599 }
600 else if(PrimEnergy > Tcut || !aModel->GetApplyCutInRange())
601 {
602 EminSec = aModel->GetSecondAdjEnergyMinForProdToProj(PrimEnergy);
603 EmaxSec = aModel->GetSecondAdjEnergyMaxForProdToProj(PrimEnergy);
604 }
605 if(EminSec >= EmaxSec)
606 return 0.;
607
608 G4bool need_to_compute = false;
609 if(aMaterial != lastMaterial || PrimEnergy != lastPrimaryEnergy ||
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 < fIndexOfAdjointEMModelInAction.size(); ++i)
626 {
627 std::size_t ind1 = fIndexOfAdjointEMModelInAction[i];
628 if(aModel == fAdjointModels[ind1] &&
629 isScatProjToProj == fIsScatProjToProj[i])
630 {
631 need_to_compute = false;
632 CS_Vs_Element = fLastAdjointCSVsModelsAndElements[ind];
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.size(); ++i)
642 {
643 if(aModel == fAdjointModels[i])
644 {
645 ind_model = i;
646 break;
647 }
648 }
649 G4double Tlow = Tcut;
650 if(!fAdjointModels[ind_model]->GetApplyCutInRange())
651 Tlow = fAdjointModels[ind_model]->GetLowEnergyLimit();
652 fIndexOfAdjointEMModelInAction.push_back(ind_model);
653 fIsScatProjToProj.push_back(isScatProjToProj);
654 CS_Vs_Element.clear();
655 if(!aModel->GetUseMatrix())
656 {
657 CS_Vs_Element.push_back(aModel->AdjointCrossSection(
658 fCurrentCouple, PrimEnergy, isScatProjToProj));
659 }
660 else if(aModel->GetUseMatrixPerElement())
661 {
662 std::size_t n_el = aMaterial->GetNumberOfElements();
663 if(aModel->GetUseOnlyOneMatrixForAllElements())
664 {
665 G4AdjointCSMatrix* theCSMatrix;
666 if(isScatProjToProj)
667 {
668 theCSMatrix = fAdjointCSMatricesForScatProjToProj[ind_model][0];
669 }
670 else
671 theCSMatrix = fAdjointCSMatricesForProdToProj[ind_model][0];
672 G4double CS = 0.;
673 if(PrimEnergy > Tlow)
674 CS = ComputeAdjointCS(PrimEnergy, theCSMatrix, Tlow);
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)->GetZ() *
679 aMaterial->GetVecNbOfAtomsPerVolume()[i];
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->GetElement(i)->GetIndex();
689 G4AdjointCSMatrix* theCSMatrix;
690 if(isScatProjToProj)
691 {
692 theCSMatrix =
693 fAdjointCSMatricesForScatProjToProj[ind_model][ind_el];
694 }
695 else
696 theCSMatrix = fAdjointCSMatricesForProdToProj[ind_model][ind_el];
697 G4double CS = 0.;
698 if(PrimEnergy > Tlow)
699 CS = ComputeAdjointCS(PrimEnergy, theCSMatrix, Tlow);
700 CS_Vs_Element.push_back(CS *
701 (aMaterial->GetVecNbOfAtomsPerVolume()[i]));
702 }
703 }
704 }
705 else
706 {
707 std::size_t ind_mat = aMaterial->GetIndex();
708 G4AdjointCSMatrix* theCSMatrix;
709 if(isScatProjToProj)
710 {
711 theCSMatrix = fAdjointCSMatricesForScatProjToProj[ind_model][ind_mat];
712 }
713 else
714 theCSMatrix = fAdjointCSMatricesForProdToProj[ind_model][ind_mat];
715 G4double CS = 0.;
716 if(PrimEnergy > Tlow)
717 CS = ComputeAdjointCS(PrimEnergy, theCSMatrix, Tlow);
718 CS_Vs_Element.push_back(CS);
719 }
720 fLastAdjointCSVsModelsAndElements.push_back(CS_Vs_Element);
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 CS instead of the CS of an
727 // element itself
728 CS += cs_vs_el;
729 }
730 return CS;
731 }
732
733 ///////////////////////////////////////////////////////
734 G4Element* G4AdjointCSManager::SampleElementFromCSMatrices(
735 G4Material* aMaterial, G4VEmAdjointModel* aModel, G4double PrimEnergy,
736 G4double Tcut, G4bool isScatProjToProj)
737 {
738 std::vector<G4double> CS_Vs_Element;
739 G4double CS = ComputeAdjointCS(aMaterial, aModel, PrimEnergy, Tcut,
740 isScatProjToProj, CS_Vs_Element);
741 G4double SumCS = 0.;
742 std::size_t ind = 0;
743 for(std::size_t i = 0; i < CS_Vs_Element.size(); ++i)
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->GetElement((G4int)ind));
754 }
755
756 ///////////////////////////////////////////////////////
757 G4double G4AdjointCSManager::ComputeTotalAdjointCS(
758 const G4MaterialCutsCouple* aCouple, G4ParticleDefinition* aPartDef,
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.size(); ++i)
769 {
770 G4double Tlow = 0.;
771 adjModel = fAdjointModels[i];
772 if(!adjModel->GetApplyCutInRange())
773 Tlow = adjModel->GetLowEnergyLimit();
774 else
775 {
776 G4ParticleDefinition* theDirSecondPartDef = GetForwardParticleEquivalent(
777 adjModel->GetAdjointEquivalentOfDirectSecondaryParticleDefinition());
778 std::size_t idx = 56;
779 if(theDirSecondPartDef->GetParticleName() == "gamma")
780 idx = 0;
781 else if(theDirSecondPartDef->GetParticleName() == "e-")
782 idx = 1;
783 else if(theDirSecondPartDef->GetParticleName() == "e+")
784 idx = 2;
785 if(idx < 56)
786 {
787 const std::vector<G4double>* aVec =
788 G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(
789 idx);
790 Tlow = (*aVec)[aCouple->GetIndex()];
791 }
792 }
793 if(Ekin <= adjModel->GetHighEnergyLimit() &&
794 Ekin >= adjModel->GetLowEnergyLimit())
795 {
796 if(aPartDef ==
797 adjModel->GetAdjointEquivalentOfDirectPrimaryParticleDefinition())
798 {
799 CS = ComputeAdjointCS(fCurrentMaterial, adjModel, Ekin, Tlow, true,
800 CS_Vs_Element);
801 TotalCS += CS;
802 (*fSigmaTableForAdjointModelScatProjToProj[i])[fCurrentMatIndex]
803 ->PutValue(fNbins, CS);
804 }
805 if(aPartDef ==
806 adjModel->GetAdjointEquivalentOfDirectSecondaryParticleDefinition())
807 {
808 CS = ComputeAdjointCS(fCurrentMaterial, adjModel, Ekin, Tlow, false,
809 CS_Vs_Element);
810 TotalCS += CS;
811 (*fSigmaTableForAdjointModelProdToProj[i])[fCurrentMatIndex]->PutValue(
812 fNbins, CS);
813 }
814 }
815 else
816 {
817 (*fSigmaTableForAdjointModelScatProjToProj[i])[fCurrentMatIndex]
818 ->PutValue(fNbins, 0.);
819 (*fSigmaTableForAdjointModelProdToProj[i])[fCurrentMatIndex]->PutValue(
820 fNbins, 0.);
821 }
822 }
823 return TotalCS;
824 }
825
826 ///////////////////////////////////////////////////////
827 std::vector<G4AdjointCSMatrix*>
828 G4AdjointCSManager::BuildCrossSectionsModelAndElement(G4VEmAdjointModel* aModel,
829 G4int Z, G4int A,
830 G4int nbin_pro_decade)
831 {
832 G4AdjointCSMatrix* theCSMatForProdToProjBackwardScattering =
833 new G4AdjointCSMatrix(false);
834 G4AdjointCSMatrix* theCSMatForScatProjToProjBackwardScattering =
835 new G4AdjointCSMatrix(true);
836
837 // make the vector of primary energy of the adjoint particle.
838 G4double EkinMin = aModel->GetLowEnergyLimit();
839 G4double EkinMaxForScat = aModel->GetHighEnergyLimit() * 0.999;
840 G4double EkinMaxForProd = aModel->GetHighEnergyLimit() * 0.999;
841 if(aModel->GetSecondPartOfSameType())
842 EkinMaxForProd = EkinMaxForProd / 2.;
843
844 // Product to projectile backward scattering
845 G4double dE = std::pow(10., 1. / nbin_pro_decade);
846 G4double E2 =
847 std::pow(10., double(int(std::log10(EkinMin) * nbin_pro_decade) + 1) /
848 nbin_pro_decade) / dE;
849 G4double E1 = EkinMin;
850 // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
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->ComputeAdjointCrossSectionVectorPerAtomForSecond(E1, Z, A,
857 nbin_pro_decade);
858 if(aMat.size() >= 2)
859 {
860 std::vector<double>* log_ESecVec = aMat[0];
861 std::vector<double>* log_CSVec = aMat[1];
862 G4double log_adjointCS = log_CSVec->back();
863 // normalise CSVec such that it becomes a probability vector
864 for(std::size_t j = 0; j < log_CSVec->size(); ++j)
865 {
866 if(j == 0)
867 (*log_CSVec)[j] = 0.;
868 else
869 (*log_CSVec)[j] =
870 std::log(1. - std::exp((*log_CSVec)[j] - log_adjointCS) + 1e-50);
871 }
872 (*log_CSVec)[log_CSVec->size() - 1] =
873 (*log_CSVec)[log_CSVec->size() - 2] - std::log(1000.);
874 theCSMatForProdToProjBackwardScattering->AddData(
875 std::log(E1), log_adjointCS, log_ESecVec, log_CSVec, 0);
876 }
877 E1 = E2;
878 E2 *= dE;
879 }
880
881 // Scattered projectile to projectile backward scattering
882 E2 = std::pow(10., G4double(int(std::log10(EkinMin) * nbin_pro_decade) + 1) /
883 nbin_pro_decade) / dE;
884 E1 = EkinMin;
885 // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
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->ComputeAdjointCrossSectionVectorPerAtomForScatProj(
892 E1, Z, A, nbin_pro_decade);
893 if(aMat.size() >= 2)
894 {
895 std::vector<G4double>* log_ESecVec = aMat[0];
896 std::vector<G4double>* log_CSVec = aMat[1];
897 G4double log_adjointCS = log_CSVec->back();
898 // normalise CSVec such that it becomes a probability vector
899 for(std::size_t j = 0; j < log_CSVec->size(); ++j)
900 {
901 if(j == 0)
902 (*log_CSVec)[j] = 0.;
903 else
904 (*log_CSVec)[j] =
905 std::log(1. - std::exp((*log_CSVec)[j] - log_adjointCS) + 1e-50);
906 }
907 (*log_CSVec)[log_CSVec->size() - 1] =
908 (*log_CSVec)[log_CSVec->size() - 2] - std::log(1000.);
909 theCSMatForScatProjToProjBackwardScattering->AddData(
910 std::log(E1), log_adjointCS, log_ESecVec, log_CSVec, 0);
911 }
912 E1 = E2;
913 E2 *= dE;
914 }
915
916 std::vector<G4AdjointCSMatrix*> res;
917 res.push_back(theCSMatForProdToProjBackwardScattering);
918 res.push_back(theCSMatForScatProjToProjBackwardScattering);
919
920 return res;
921 }
922
923 ///////////////////////////////////////////////////////
924 std::vector<G4AdjointCSMatrix*>
925 G4AdjointCSManager::BuildCrossSectionsModelAndMaterial(
926 G4VEmAdjointModel* aModel, G4Material* aMaterial, G4int nbin_pro_decade)
927 {
928 G4AdjointCSMatrix* theCSMatForProdToProjBackwardScattering =
929 new G4AdjointCSMatrix(false);
930 G4AdjointCSMatrix* theCSMatForScatProjToProjBackwardScattering =
931 new G4AdjointCSMatrix(true);
932
933 G4double EkinMin = aModel->GetLowEnergyLimit();
934 G4double EkinMaxForScat = aModel->GetHighEnergyLimit() * 0.999;
935 G4double EkinMaxForProd = aModel->GetHighEnergyLimit() * 0.999;
936 if(aModel->GetSecondPartOfSameType())
937 EkinMaxForProd /= 2.;
938
939 // Product to projectile backward scattering
940 G4double dE = std::pow(10., 1. / nbin_pro_decade);
941 G4double E2 =
942 std::pow(10., double(int(std::log10(EkinMin) * nbin_pro_decade) + 1) /
943 nbin_pro_decade) / dE;
944 G4double E1 = EkinMin;
945 // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
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->ComputeAdjointCrossSectionVectorPerVolumeForSecond(
952 aMaterial, E1, nbin_pro_decade);
953 if(aMat.size() >= 2)
954 {
955 std::vector<G4double>* log_ESecVec = aMat[0];
956 std::vector<G4double>* log_CSVec = aMat[1];
957 G4double log_adjointCS = log_CSVec->back();
958
959 // normalise CSVec such that it becomes a probability vector
960 for(std::size_t j = 0; j < log_CSVec->size(); ++j)
961 {
962 if(j == 0)
963 (*log_CSVec)[j] = 0.;
964 else
965 (*log_CSVec)[j] =
966 std::log(1. - std::exp((*log_CSVec)[j] - log_adjointCS));
967 }
968 (*log_CSVec)[log_CSVec->size() - 1] =
969 (*log_CSVec)[log_CSVec->size() - 2] - std::log(1000.);
970 theCSMatForProdToProjBackwardScattering->AddData(
971 std::log(E1), log_adjointCS, log_ESecVec, log_CSVec, 0);
972 }
973
974 E1 = E2;
975 E2 *= dE;
976 }
977
978 // Scattered projectile to projectile backward scattering
979 E2 =
980 std::pow(10., G4double(G4int(std::log10(EkinMin) * nbin_pro_decade) + 1) /
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->ComputeAdjointCrossSectionVectorPerVolumeForScatProj(
990 aMaterial, E1, nbin_pro_decade);
991 if(aMat.size() >= 2)
992 {
993 std::vector<G4double>* log_ESecVec = aMat[0];
994 std::vector<G4double>* log_CSVec = aMat[1];
995 G4double log_adjointCS = log_CSVec->back();
996
997 for(std::size_t j = 0; j < log_CSVec->size(); ++j)
998 {
999 if(j == 0)
1000 (*log_CSVec)[j] = 0.;
1001 else
1002 (*log_CSVec)[j] =
1003 std::log(1. - std::exp((*log_CSVec)[j] - log_adjointCS));
1004 }
1005 (*log_CSVec)[log_CSVec->size() - 1] =
1006 (*log_CSVec)[log_CSVec->size() - 2] - std::log(1000.);
1007
1008 theCSMatForScatProjToProjBackwardScattering->AddData(
1009 std::log(E1), log_adjointCS, log_ESecVec, log_CSVec, 0);
1010 }
1011 E1 = E2;
1012 E2 *= dE;
1013 }
1014
1015 std::vector<G4AdjointCSMatrix*> res;
1016 res.push_back(theCSMatForProdToProjBackwardScattering);
1017 res.push_back(theCSMatForScatProjToProjBackwardScattering);
1018
1019 return res;
1020 }
1021
1022 ///////////////////////////////////////////////////////
1023 G4ParticleDefinition* G4AdjointCSManager::GetAdjointParticleEquivalent(
1024 G4ParticleDefinition* theFwdPartDef)
1025 {
1026 if(theFwdPartDef->GetParticleName() == "e-")
1027 return G4AdjointElectron::AdjointElectron();
1028 else if(theFwdPartDef->GetParticleName() == "gamma")
1029 return G4AdjointGamma::AdjointGamma();
1030 else if(theFwdPartDef->GetParticleName() == "proton")
1031 return G4AdjointProton::AdjointProton();
1032 else if(theFwdPartDef == fFwdIon)
1033 return fAdjIon;
1034 return nullptr;
1035 }
1036
1037 ///////////////////////////////////////////////////////
1038 G4ParticleDefinition* G4AdjointCSManager::GetForwardParticleEquivalent(
1039 G4ParticleDefinition* theAdjPartDef)
1040 {
1041 if(theAdjPartDef->GetParticleName() == "adj_e-")
1042 return G4Electron::Electron();
1043 else if(theAdjPartDef->GetParticleName() == "adj_gamma")
1044 return G4Gamma::Gamma();
1045 else if(theAdjPartDef->GetParticleName() == "adj_proton")
1046 return G4Proton::Proton();
1047 else if(theAdjPartDef == fAdjIon)
1048 return fFwdIon;
1049 return nullptr;
1050 }
1051
1052 ///////////////////////////////////////////////////////
1053 void G4AdjointCSManager::DefineCurrentMaterial(
1054 const G4MaterialCutsCouple* couple)
1055 {
1056 if(couple != fCurrentCouple)
1057 {
1058 fCurrentCouple = const_cast<G4MaterialCutsCouple*>(couple);
1059 fCurrentMaterial = const_cast<G4Material*>(couple->GetMaterial());
1060 fCurrentMatIndex = couple->GetIndex();
1061 fLastCSCorrectionFactor = 1.;
1062 }
1063 }
1064
1065 ///////////////////////////////////////////////////////
1066 void G4AdjointCSManager::DefineCurrentParticle(
1067 const G4ParticleDefinition* aPartDef)
1068 {
1069
1070 if(aPartDef != fCurrentParticleDef)
1071 {
1072 fCurrentParticleDef = const_cast<G4ParticleDefinition*>(aPartDef);
1073 fMassRatio = 1.;
1074 if(aPartDef == fAdjIon)
1075 fMassRatio = proton_mass_c2 / aPartDef->GetPDGMass();
1076 fCurrentParticleIndex = 1000000;
1077 for(std::size_t i = 0; i < fAdjointParticlesInAction.size(); ++i)
1078 {
1079 if(aPartDef == fAdjointParticlesInAction[i])
1080 fCurrentParticleIndex = i;
1081 }
1082 }
1083 }
1084
1085 ////////////////////////////////////////////////////////////////////////////////
1086 G4double G4AdjointCSManager::ComputeAdjointCS(
1087 G4double aPrimEnergy, G4AdjointCSMatrix* anAdjointCSMatrix, G4double Tcut)
1088 {
1089 std::vector<double>* theLogPrimEnergyVector =
1090 anAdjointCSMatrix->GetLogPrimEnergyVector();
1091 if(theLogPrimEnergyVector->empty())
1092 {
1093 G4cout << "No data are contained in the given AdjointCSMatrix!" << G4endl;
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 > theLogPrimEnergyVector->back())
1100 return 0.;
1101
1102 G4AdjointInterpolator* theInterpolator = G4AdjointInterpolator::GetInstance();
1103
1104 std::size_t ind =
1105 theInterpolator->FindPositionForLogVector(log_E, *theLogPrimEnergyVector);
1106 G4double aLogPrimEnergy1, aLogPrimEnergy2;
1107 G4double aLogCS1, aLogCS2;
1108 G4double log01, log02;
1109 std::vector<G4double>* aLogSecondEnergyVector1 = nullptr;
1110 std::vector<G4double>* aLogSecondEnergyVector2 = nullptr;
1111 std::vector<G4double>* aLogProbVector1 = nullptr;
1112 std::vector<G4double>* aLogProbVector2 = nullptr;
1113 std::vector<std::size_t>* aLogProbVectorIndex1 = nullptr;
1114 std::vector<std::size_t>* aLogProbVectorIndex2 = nullptr;
1115
1116 anAdjointCSMatrix->GetData((G4int)ind, aLogPrimEnergy1, aLogCS1, log01,
1117 aLogSecondEnergyVector1, aLogProbVector1,
1118 aLogProbVectorIndex1);
1119 anAdjointCSMatrix->GetData(G4int(ind + 1), aLogPrimEnergy2, aLogCS2, log02,
1120 aLogSecondEnergyVector2, aLogProbVector2,
1121 aLogProbVectorIndex2);
1122 if (! (aLogProbVector1 && aLogProbVector2 &&
1123 aLogSecondEnergyVector1 && aLogSecondEnergyVector2)){
1124 return 0.;
1125 }
1126
1127 if(anAdjointCSMatrix->IsScatProjToProj())
1128 { // case where the Tcut plays a role
1129 G4double log_minimum_prob1, log_minimum_prob2;
1130 log_minimum_prob1 = theInterpolator->InterpolateForLogVector(
1131 log_Tcut, *aLogSecondEnergyVector1, *aLogProbVector1);
1132 log_minimum_prob2 = theInterpolator->InterpolateForLogVector(
1133 log_Tcut, *aLogSecondEnergyVector2, *aLogProbVector2);
1134 aLogCS1 += log_minimum_prob1;
1135 aLogCS2 += log_minimum_prob2;
1136 }
1137
1138 G4double log_adjointCS = theInterpolator->LinearInterpolation(
1139 log_E, aLogPrimEnergy1, aLogPrimEnergy2, aLogCS1, aLogCS2);
1140 return std::exp(log_adjointCS);
1141 }
1142