Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4AdjointCSManager.cc
Go to the documentation of this file.
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// $Id$
27//
28
29#include <fstream>
30#include <iomanip>
31
32#include "G4AdjointCSManager.hh"
33
34#include "G4PhysicalConstants.hh"
35#include "G4SystemOfUnits.hh"
36#include "G4AdjointCSMatrix.hh"
37#include "G4AdjointInterpolator.hh"
38#include "G4AdjointCSMatrix.hh"
39#include "G4VEmAdjointModel.hh"
40#include "G4ElementTable.hh"
41#include "G4Element.hh"
42#include "G4ParticleDefinition.hh"
43#include "G4Element.hh"
44#include "G4VEmProcess.hh"
45#include "G4VEnergyLossProcess.hh"
46#include "G4PhysicsTable.hh"
47#include "G4PhysicsLogVector.hh"
48#include "G4PhysicsTableHelper.hh"
49#include "G4Electron.hh"
50#include "G4Gamma.hh"
51#include "G4Proton.hh"
52#include "G4AdjointElectron.hh"
53#include "G4AdjointGamma.hh"
54#include "G4AdjointProton.hh"
55#include "G4ProductionCutsTable.hh"
56#include "G4ProductionCutsTable.hh"
57
58G4AdjointCSManager* G4AdjointCSManager::theInstance = 0;
59///////////////////////////////////////////////////////
60//
61G4AdjointCSManager* G4AdjointCSManager::GetAdjointCSManager()
62{ if(theInstance == 0) {
63 static G4AdjointCSManager ins;
64 theInstance = &ins;
65 }
66 return theInstance;
67}
68
69///////////////////////////////////////////////////////
70//
71G4AdjointCSManager::G4AdjointCSManager()
72{ CrossSectionMatrixesAreBuilt=false;
73 TotalSigmaTableAreBuilt=false;
74 theTotalForwardSigmaTableVector.clear();
75 theTotalAdjointSigmaTableVector.clear();
76 listOfForwardEmProcess.clear();
77 listOfForwardEnergyLossProcess.clear();
78 theListOfAdjointParticlesInAction.clear();
79 EminForFwdSigmaTables.clear();
80 EminForAdjSigmaTables.clear();
81 EkinofFwdSigmaMax.clear();
82 EkinofAdjSigmaMax.clear();
83 listSigmaTableForAdjointModelScatProjToProj.clear();
84 listSigmaTableForAdjointModelProdToProj.clear();
85 Tmin=0.1*keV;
86 Tmax=100.*TeV;
87 nbins=320; //probably this should be decrease, that was choosen to avoid error in the CS value closed to CS jump.(For example at Tcut)
88
89 RegisterAdjointParticle(G4AdjointElectron::AdjointElectron());
90 RegisterAdjointParticle(G4AdjointGamma::AdjointGamma());
91 RegisterAdjointParticle(G4AdjointProton::AdjointProton());
92
93 verbose = 1;
94
95 lastPartDefForCS =0;
96 LastEkinForCS =0;
97 LastCSCorrectionFactor =1.;
98
99 forward_CS_mode = true;
100
101 currentParticleDef = 0;
102 currentCouple =0;
103 currentMaterial=0;
104 lastMaterial=0;
105
106
107 theAdjIon = 0;
108 theFwdIon = 0;
109
110
111
112
113}
114///////////////////////////////////////////////////////
115//
116G4AdjointCSManager::~G4AdjointCSManager()
117{;
118}
119///////////////////////////////////////////////////////
120//
121size_t G4AdjointCSManager::RegisterEmAdjointModel(G4VEmAdjointModel* aModel)
122{listOfAdjointEMModel.push_back(aModel);
123 listSigmaTableForAdjointModelScatProjToProj.push_back(new G4PhysicsTable);
124 listSigmaTableForAdjointModelProdToProj.push_back(new G4PhysicsTable);
125 return listOfAdjointEMModel.size() -1;
126
127}
128///////////////////////////////////////////////////////
129//
130void G4AdjointCSManager::RegisterEmProcess(G4VEmProcess* aProcess, G4ParticleDefinition* aFwdPartDef)
131{
132 G4ParticleDefinition* anAdjPartDef = GetAdjointParticleEquivalent(aFwdPartDef);
133 if (anAdjPartDef && aProcess){
134 RegisterAdjointParticle(anAdjPartDef);
135 G4int index=-1;
136
137 for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){
138 if (anAdjPartDef->GetParticleName() == theListOfAdjointParticlesInAction[i]->GetParticleName()) index=i;
139 }
140 listOfForwardEmProcess[index]->push_back(aProcess);
141 }
142}
143///////////////////////////////////////////////////////
144//
145void G4AdjointCSManager::RegisterEnergyLossProcess(G4VEnergyLossProcess* aProcess, G4ParticleDefinition* aFwdPartDef)
146{
147 G4ParticleDefinition* anAdjPartDef = GetAdjointParticleEquivalent(aFwdPartDef);
148 if (anAdjPartDef && aProcess){
149 RegisterAdjointParticle(anAdjPartDef);
150 G4int index=-1;
151 for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){
152 if (anAdjPartDef->GetParticleName() == theListOfAdjointParticlesInAction[i]->GetParticleName()) index=i;
153 }
154 listOfForwardEnergyLossProcess[index]->push_back(aProcess);
155 }
156}
157///////////////////////////////////////////////////////
158//
159void G4AdjointCSManager::RegisterAdjointParticle(G4ParticleDefinition* aPartDef)
160{ G4int index=-1;
161 for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){
162 if (aPartDef->GetParticleName() == theListOfAdjointParticlesInAction[i]->GetParticleName()) index=i;
163 }
164
165 if (index ==-1){
166 listOfForwardEnergyLossProcess.push_back(new std::vector<G4VEnergyLossProcess*>());
167 theTotalForwardSigmaTableVector.push_back(new G4PhysicsTable);
168 theTotalAdjointSigmaTableVector.push_back(new G4PhysicsTable);
169 listOfForwardEmProcess.push_back(new std::vector<G4VEmProcess*>());
170 theListOfAdjointParticlesInAction.push_back(aPartDef);
171 EminForFwdSigmaTables.push_back(std::vector<G4double> ());
172 EminForAdjSigmaTables.push_back(std::vector<G4double> ());
173 EkinofFwdSigmaMax.push_back(std::vector<G4double> ());
174 EkinofAdjSigmaMax.push_back(std::vector<G4double> ());
175
176 }
177}
178///////////////////////////////////////////////////////
179//
180void G4AdjointCSManager::BuildCrossSectionMatrices()
181{
182 if (CrossSectionMatrixesAreBuilt) return;
183 //Tcut, Tmax
184 //The matrices will be computed probably just once
185 //When Tcut will change some PhysicsTable will be recomputed
186 // for each MaterialCutCouple but not all the matrices
187 //The Tcut defines a lower limit in the energy of the Projectile before the scattering
188 //In the Projectile to Scattered Projectile case we have
189 // E_ScatProj<E_Proj-Tcut
190 //Therefore in the adjoint case we have
191 // Eproj> E_ScatProj+Tcut
192 //This implies that when computing the adjoint CS we should integrate over Epro
193 // from E_ScatProj+Tcut to Emax
194 //In the Projectile to Secondary case Tcut plays a role only in the fact that
195 // Esecond should be greater than Tcut to have the possibility to have any adjoint
196 //process
197 //To avoid to recompute the matrices for all changes of MaterialCutCouple
198 //We propose to compute the matrices only once for the minimum possible Tcut and then
199 //to interpolate the probability for a new Tcut (implemented in G4VAdjointEmModel)
200
201
202 theAdjointCSMatricesForScatProjToProj.clear();
203 theAdjointCSMatricesForProdToProj.clear();
204 const G4ElementTable* theElementTable = G4Element::GetElementTable();
205 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
206
207 G4cout<<"========== Computation of cross section matrices for adjoint models =========="<<G4endl;
208 for (size_t i=0; i<listOfAdjointEMModel.size();i++){
209 G4VEmAdjointModel* aModel =listOfAdjointEMModel[i];
210 G4cout<<"Build adjoint cross section matrices for "<<aModel->GetName()<<G4endl;
211 if (aModel->GetUseMatrix()){
212 std::vector<G4AdjointCSMatrix*>* aListOfMat1 = new std::vector<G4AdjointCSMatrix*>();
213 std::vector<G4AdjointCSMatrix*>* aListOfMat2 = new std::vector<G4AdjointCSMatrix*>();
214 aListOfMat1->clear();
215 aListOfMat2->clear();
216 if (aModel->GetUseMatrixPerElement()){
217 if (aModel->GetUseOnlyOneMatrixForAllElements()){
218 std::vector<G4AdjointCSMatrix*>
219 two_matrices=BuildCrossSectionsMatricesForAGivenModelAndElement(aModel,1, 1, 80);
220 aListOfMat1->push_back(two_matrices[0]);
221 aListOfMat2->push_back(two_matrices[1]);
222 }
223 else {
224 for (size_t j=0; j<theElementTable->size();j++){
225 G4Element* anElement=(*theElementTable)[j];
226 G4int Z = int(anElement->GetZ());
227 G4int A = int(anElement->GetA());
228 std::vector<G4AdjointCSMatrix*>
229 two_matrices=BuildCrossSectionsMatricesForAGivenModelAndElement(aModel,Z, A, 40);
230 aListOfMat1->push_back(two_matrices[0]);
231 aListOfMat2->push_back(two_matrices[1]);
232 }
233 }
234 }
235 else { //Per material case
236 for (size_t j=0; j<theMaterialTable->size();j++){
237 G4Material* aMaterial=(*theMaterialTable)[j];
238 std::vector<G4AdjointCSMatrix*>
239 two_matrices=BuildCrossSectionsMatricesForAGivenModelAndMaterial(aModel,aMaterial, 40);
240 aListOfMat1->push_back(two_matrices[0]);
241 aListOfMat2->push_back(two_matrices[1]);
242 }
243
244 }
245 theAdjointCSMatricesForProdToProj.push_back(*aListOfMat1);
246 theAdjointCSMatricesForScatProjToProj.push_back(*aListOfMat2);
247 aModel->SetCSMatrices(aListOfMat1, aListOfMat2);
248 }
249 else { G4cout<<"The model "<<aModel->GetName()<<" does not use cross section matrices"<<G4endl;
250 std::vector<G4AdjointCSMatrix*> two_empty_matrices;
251 theAdjointCSMatricesForProdToProj.push_back(two_empty_matrices);
252 theAdjointCSMatricesForScatProjToProj.push_back(two_empty_matrices);
253
254 }
255 }
256 G4cout<<" All adjoint cross section matrices are computed!"<<G4endl;
257 G4cout<<"======================================================================"<<G4endl;
258
259 CrossSectionMatrixesAreBuilt = true;
260
261
262}
263
264
265///////////////////////////////////////////////////////
266//
267void G4AdjointCSManager::BuildTotalSigmaTables()
268{ if (TotalSigmaTableAreBuilt) return;
269
270
271 const G4ProductionCutsTable* theCoupleTable= G4ProductionCutsTable::GetProductionCutsTable();
272
273
274 //Prepare the Sigma table for all AdjointEMModel, will be filled later on
275 for (size_t i=0; i<listOfAdjointEMModel.size();i++){
276 listSigmaTableForAdjointModelScatProjToProj[i]->clearAndDestroy();
277 listSigmaTableForAdjointModelProdToProj[i]->clearAndDestroy();
278 for (size_t j=0;j<theCoupleTable->GetTableSize();j++){
279 listSigmaTableForAdjointModelScatProjToProj[i]->push_back(new G4PhysicsLogVector(Tmin, Tmax, nbins));
280 listSigmaTableForAdjointModelProdToProj[i]->push_back(new G4PhysicsLogVector(Tmin, Tmax, nbins));
281 }
282 }
283
284
285
286 for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){
287 G4ParticleDefinition* thePartDef = theListOfAdjointParticlesInAction[i];
288 DefineCurrentParticle(thePartDef);
289 theTotalForwardSigmaTableVector[i]->clearAndDestroy();
290 theTotalAdjointSigmaTableVector[i]->clearAndDestroy();
291 EminForFwdSigmaTables[i].clear();
292 EminForAdjSigmaTables[i].clear();
293 EkinofFwdSigmaMax[i].clear();
294 EkinofAdjSigmaMax[i].clear();
295 //G4cout<<thePartDef->GetParticleName();
296
297 for (size_t j=0;j<theCoupleTable->GetTableSize();j++){
298 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(j);
299
300 /*
301 G4String file_name1=couple->GetMaterial()->GetName()+"_"+thePartDef->GetParticleName()+"_adj_totCS.txt";
302 G4String file_name2=couple->GetMaterial()->GetName()+"_"+thePartDef->GetParticleName()+"_fwd_totCS.txt";
303
304 std::fstream FileOutputAdjCS(file_name1, std::ios::out);
305 std::fstream FileOutputFwdCS(file_name2, std::ios::out);
306
307
308
309 FileOutputAdjCS<<std::setiosflags(std::ios::scientific);
310 FileOutputAdjCS<<std::setprecision(6);
311 FileOutputFwdCS<<std::setiosflags(std::ios::scientific);
312 FileOutputFwdCS<<std::setprecision(6);
313 */
314
315
316 //make first the total fwd CS table for FwdProcess
317 G4PhysicsVector* aVector = new G4PhysicsLogVector(Tmin, Tmax, nbins);
318 G4bool Emin_found=false;
319 G4double sigma_max =0.;
320 G4double e_sigma_max =0.;
321 for(size_t l=0; l<aVector->GetVectorLength(); l++) {
322 G4double totCS=0.;
323 G4double e=aVector->GetLowEdgeEnergy(l);
324 for (size_t k=0; k<listOfForwardEmProcess[i]->size(); k++){
325 totCS+=(*listOfForwardEmProcess[i])[k]->GetLambda(e, couple);
326 }
327 for (size_t k=0; k<listOfForwardEnergyLossProcess[i]->size(); k++){
328 if (thePartDef == theAdjIon) { // e is considered already as the scaled energy
329 size_t mat_index = couple->GetIndex();
330 G4VEmModel* currentModel = (*listOfForwardEnergyLossProcess[i])[k]->SelectModelForMaterial(e,mat_index);
331 G4double chargeSqRatio = currentModel->GetChargeSquareRatio(theFwdIon,couple->GetMaterial(),e/massRatio);
332 (*listOfForwardEnergyLossProcess[i])[k]->SetDynamicMassCharge(massRatio,chargeSqRatio);
333 }
334 G4double e1=e/massRatio;
335 totCS+=(*listOfForwardEnergyLossProcess[i])[k]->GetLambda(e1, couple);
336 }
337 aVector->PutValue(l,totCS);
338 if (totCS>sigma_max){
339 sigma_max=totCS;
340 e_sigma_max = e;
341
342 }
343 //FileOutputFwdCS<<e<<'\t'<<totCS<<G4endl;
344
345 if (totCS>0 && !Emin_found) {
346 EminForFwdSigmaTables[i].push_back(e);
347 Emin_found=true;
348 }
349
350
351 }
352 //FileOutputFwdCS.close();
353
354 EkinofFwdSigmaMax[i].push_back(e_sigma_max);
355
356
357 if(!Emin_found) EminForFwdSigmaTables[i].push_back(Tmax);
358
359 theTotalForwardSigmaTableVector[i]->push_back(aVector);
360
361
362 Emin_found=false;
363 sigma_max=0;
364 e_sigma_max =0.;
365 G4PhysicsVector* aVector1 = new G4PhysicsLogVector(Tmin, Tmax, nbins);
366 for(eindex=0; eindex<aVector->GetVectorLength(); eindex++) {
367 G4double e=aVector->GetLowEdgeEnergy(eindex);
368 G4double totCS =ComputeTotalAdjointCS(couple,thePartDef,e*0.9999999/massRatio); //massRatio needed for ions
369 aVector1->PutValue(eindex,totCS);
370 if (totCS>sigma_max){
371 sigma_max=totCS;
372 e_sigma_max = e;
373
374 }
375 //FileOutputAdjCS<<e<<'\t'<<totCS<<G4endl;
376 if (totCS>0 && !Emin_found) {
377 EminForAdjSigmaTables[i].push_back(e);
378 Emin_found=true;
379 }
380
381 }
382 //FileOutputAdjCS.close();
383 EkinofAdjSigmaMax[i].push_back(e_sigma_max);
384 if(!Emin_found) EminForAdjSigmaTables[i].push_back(Tmax);
385
386 theTotalAdjointSigmaTableVector[i]->push_back(aVector1);
387
388 }
389 }
390 TotalSigmaTableAreBuilt =true;
391
392}
393///////////////////////////////////////////////////////
394//
395G4double G4AdjointCSManager::GetTotalAdjointCS(G4ParticleDefinition* aPartDef, G4double Ekin,
396 const G4MaterialCutsCouple* aCouple)
397{ DefineCurrentMaterial(aCouple);
398 DefineCurrentParticle(aPartDef);
399 G4bool b;
400 return (((*theTotalAdjointSigmaTableVector[currentParticleIndex])[currentMatIndex])->GetValue(Ekin*massRatio, b));
401
402
403
404}
405///////////////////////////////////////////////////////
406//
407G4double G4AdjointCSManager::GetTotalForwardCS(G4ParticleDefinition* aPartDef, G4double Ekin,
408 const G4MaterialCutsCouple* aCouple)
409{ DefineCurrentMaterial(aCouple);
410 DefineCurrentParticle(aPartDef);
411 G4bool b;
412 return (((*theTotalForwardSigmaTableVector[currentParticleIndex])[currentMatIndex])->GetValue(Ekin*massRatio, b));
413
414
415}
416///////////////////////////////////////////////////////
417//
418G4double G4AdjointCSManager::GetAdjointSigma(G4double Ekin_nuc, size_t index_model,G4bool is_scat_proj_to_proj,
419 const G4MaterialCutsCouple* aCouple)
420{ DefineCurrentMaterial(aCouple);
421 G4bool b;
422 if (is_scat_proj_to_proj) return (((*listSigmaTableForAdjointModelScatProjToProj[index_model])[currentMatIndex])->GetValue(Ekin_nuc, b));
423 else return (((*listSigmaTableForAdjointModelProdToProj[index_model])[currentMatIndex])->GetValue(Ekin_nuc, b));
424}
425///////////////////////////////////////////////////////
426//
427void G4AdjointCSManager::GetEminForTotalCS(G4ParticleDefinition* aPartDef,
428 const G4MaterialCutsCouple* aCouple, G4double& emin_adj, G4double& emin_fwd)
429{ DefineCurrentMaterial(aCouple);
430 DefineCurrentParticle(aPartDef);
431 emin_adj = EminForAdjSigmaTables[currentParticleIndex][currentMatIndex]/massRatio;
432 emin_fwd = EminForFwdSigmaTables[currentParticleIndex][currentMatIndex]/massRatio;
433
434
435
436}
437///////////////////////////////////////////////////////
438//
439void G4AdjointCSManager::GetMaxFwdTotalCS(G4ParticleDefinition* aPartDef,
440 const G4MaterialCutsCouple* aCouple, G4double& e_sigma_max, G4double& sigma_max)
441{ DefineCurrentMaterial(aCouple);
442 DefineCurrentParticle(aPartDef);
443 e_sigma_max = EkinofFwdSigmaMax[currentParticleIndex][currentMatIndex];
444 G4bool b;
445 sigma_max =((*theTotalForwardSigmaTableVector[currentParticleIndex])[currentMatIndex])->GetValue(e_sigma_max, b);
446 e_sigma_max/=massRatio;
447
448
449}
450///////////////////////////////////////////////////////
451//
452void G4AdjointCSManager::GetMaxAdjTotalCS(G4ParticleDefinition* aPartDef,
453 const G4MaterialCutsCouple* aCouple, G4double& e_sigma_max, G4double& sigma_max)
454{ DefineCurrentMaterial(aCouple);
455 DefineCurrentParticle(aPartDef);
456 e_sigma_max = EkinofAdjSigmaMax[currentParticleIndex][currentMatIndex];
457 G4bool b;
458 sigma_max =((*theTotalAdjointSigmaTableVector[currentParticleIndex])[currentMatIndex])->GetValue(e_sigma_max, b);
459 e_sigma_max/=massRatio;
460
461
462}
463///////////////////////////////////////////////////////
464//
465G4double G4AdjointCSManager::GetCrossSectionCorrection(G4ParticleDefinition* aPartDef,G4double PreStepEkin,const G4MaterialCutsCouple* aCouple, G4bool& fwd_is_used,
466 G4double& fwd_TotCS)
467{ G4double corr_fac = 1.;
468 if (forward_CS_mode) {
469 fwd_TotCS=PrefwdCS;
470 if (LastEkinForCS != PreStepEkin || aPartDef != lastPartDefForCS || aCouple!=currentCouple) {
471 DefineCurrentMaterial(aCouple);
472 PreadjCS = GetTotalAdjointCS(aPartDef, PreStepEkin,aCouple);
473 PrefwdCS = GetTotalForwardCS(aPartDef, PreStepEkin,aCouple);
474 LastEkinForCS = PreStepEkin;
475 lastPartDefForCS = aPartDef;
476 if (PrefwdCS >0. && PreadjCS >0.) {
477 forward_CS_is_used = true;
478 LastCSCorrectionFactor = PrefwdCS/PreadjCS;
479 }
480 else {
481 forward_CS_is_used = false;
482 LastCSCorrectionFactor = 1.;
483
484 }
485
486 }
487 corr_fac =LastCSCorrectionFactor;
488
489
490
491 }
492 else {
493 forward_CS_is_used = false;
494 LastCSCorrectionFactor = 1.;
495 }
496 fwd_TotCS=PrefwdCS;
497 fwd_is_used = forward_CS_is_used;
498 return corr_fac;
499}
500///////////////////////////////////////////////////////
501//
502G4double G4AdjointCSManager::GetContinuousWeightCorrection(G4ParticleDefinition* aPartDef, G4double PreStepEkin,G4double AfterStepEkin,
503 const G4MaterialCutsCouple* aCouple, G4double step_length)
504{ G4double corr_fac = 1.;
505 //return corr_fac;
506 //G4double after_adjCS = GetTotalAdjointCS(aPartDef, AfterStepEkin,aCouple);
507 G4double after_fwdCS = GetTotalForwardCS(aPartDef, AfterStepEkin,aCouple);
508 G4double pre_adjCS = GetTotalAdjointCS(aPartDef, PreStepEkin,aCouple);
509 if (!forward_CS_is_used || pre_adjCS ==0. || after_fwdCS==0.) {
510 forward_CS_is_used=false;
511 G4double pre_fwdCS = GetTotalForwardCS(aPartDef, PreStepEkin,aCouple);
512 corr_fac *=std::exp((pre_adjCS-pre_fwdCS)*step_length);
513 LastCSCorrectionFactor = 1.;
514 }
515 else {
516 LastCSCorrectionFactor = after_fwdCS/pre_adjCS;
517 }
518
519
520
521 return corr_fac;
522}
523///////////////////////////////////////////////////////
524//
525G4double G4AdjointCSManager::GetPostStepWeightCorrection( )
526{//return 1.;
527 return 1./LastCSCorrectionFactor;
528
529}
530///////////////////////////////////////////////////////
531//
532G4double G4AdjointCSManager::ComputeAdjointCS(G4Material* aMaterial,
533 G4VEmAdjointModel* aModel,
534 G4double PrimEnergy,
535 G4double Tcut,
536 G4bool IsScatProjToProjCase,
537 std::vector<G4double>& CS_Vs_Element)
538{
539
540 G4double EminSec=0;
541 G4double EmaxSec=0;
542
543 if (IsScatProjToProjCase){
544 EminSec= aModel->GetSecondAdjEnergyMinForScatProjToProjCase(PrimEnergy,Tcut);
545 EmaxSec= aModel->GetSecondAdjEnergyMaxForScatProjToProjCase(PrimEnergy);
546 }
547 else if (PrimEnergy > Tcut || !aModel->GetApplyCutInRange()) {
548 EminSec= aModel->GetSecondAdjEnergyMinForProdToProjCase(PrimEnergy);
549 EmaxSec= aModel->GetSecondAdjEnergyMaxForProdToProjCase(PrimEnergy);
550 }
551 if (EminSec >= EmaxSec) return 0.;
552
553
554 G4bool need_to_compute=false;
555 if ( aMaterial!= lastMaterial || PrimEnergy != lastPrimaryEnergy || Tcut != lastTcut){
556 lastMaterial =aMaterial;
557 lastPrimaryEnergy = PrimEnergy;
558 lastTcut=Tcut;
559 listOfIndexOfAdjointEMModelInAction.clear();
560 listOfIsScatProjToProjCase.clear();
561 lastAdjointCSVsModelsAndElements.clear();
562 need_to_compute=true;
563
564 }
565 size_t ind=0;
566 if (!need_to_compute){
567 need_to_compute=true;
568 for (size_t i=0;i<listOfIndexOfAdjointEMModelInAction.size();i++){
569 size_t ind1=listOfIndexOfAdjointEMModelInAction[i];
570 if (aModel == listOfAdjointEMModel[ind1] && IsScatProjToProjCase == listOfIsScatProjToProjCase[i]){
571 need_to_compute=false;
572 CS_Vs_Element = lastAdjointCSVsModelsAndElements[ind];
573 }
574 ind++;
575 }
576 }
577
578 if (need_to_compute){
579 size_t ind_model=0;
580 for (size_t i=0;i<listOfAdjointEMModel.size();i++){
581 if (aModel == listOfAdjointEMModel[i]){
582 ind_model=i;
583 break;
584 }
585 }
586 G4double Tlow=Tcut;
587 if (!listOfAdjointEMModel[ind_model]->GetApplyCutInRange()) Tlow =listOfAdjointEMModel[ind_model]->GetLowEnergyLimit();
588 listOfIndexOfAdjointEMModelInAction.push_back(ind_model);
589 listOfIsScatProjToProjCase.push_back(IsScatProjToProjCase);
590 CS_Vs_Element.clear();
591 if (!aModel->GetUseMatrix()){
592 CS_Vs_Element.push_back(aModel->AdjointCrossSection(currentCouple,PrimEnergy,IsScatProjToProjCase));
593
594
595 }
596 else if (aModel->GetUseMatrixPerElement()){
597 size_t n_el = aMaterial->GetNumberOfElements();
598 if (aModel->GetUseOnlyOneMatrixForAllElements()){
599 G4AdjointCSMatrix* theCSMatrix;
600 if (IsScatProjToProjCase){
601 theCSMatrix=theAdjointCSMatricesForScatProjToProj[ind_model][0];
602 }
603 else theCSMatrix=theAdjointCSMatricesForProdToProj[ind_model][0];
604 G4double CS =0.;
605 if (PrimEnergy > Tlow)
606 CS = ComputeAdjointCS(PrimEnergy,theCSMatrix,Tlow);
607 G4double factor=0.;
608 for (size_t i=0;i<n_el;i++){ //this could be computed only once
609 //size_t ind_el = aMaterial->GetElement(i)->GetIndex();
610 factor+=aMaterial->GetElement(i)->GetZ()*aMaterial->GetVecNbOfAtomsPerVolume()[i];
611 }
612 CS *=factor;
613 CS_Vs_Element.push_back(CS);
614
615 }
616 else {
617 for (size_t i=0;i<n_el;i++){
618 size_t ind_el = aMaterial->GetElement(i)->GetIndex();
619 //G4cout<<aMaterial->GetName()<<G4endl;
620 G4AdjointCSMatrix* theCSMatrix;
621 if (IsScatProjToProjCase){
622 theCSMatrix=theAdjointCSMatricesForScatProjToProj[ind_model][ind_el];
623 }
624 else theCSMatrix=theAdjointCSMatricesForProdToProj[ind_model][ind_el];
625 G4double CS =0.;
626 if (PrimEnergy > Tlow)
627 CS = ComputeAdjointCS(PrimEnergy,theCSMatrix,Tlow);
628 //G4cout<<CS<<G4endl;
629 CS_Vs_Element.push_back(CS*(aMaterial->GetVecNbOfAtomsPerVolume()[i]));
630 }
631 }
632
633 }
634 else {
635 size_t ind_mat = aMaterial->GetIndex();
636 G4AdjointCSMatrix* theCSMatrix;
637 if (IsScatProjToProjCase){
638 theCSMatrix=theAdjointCSMatricesForScatProjToProj[ind_model][ind_mat];
639 }
640 else theCSMatrix=theAdjointCSMatricesForProdToProj[ind_model][ind_mat];
641 G4double CS =0.;
642 if (PrimEnergy > Tlow)
643 CS = ComputeAdjointCS(PrimEnergy,theCSMatrix,Tlow);
644 CS_Vs_Element.push_back(CS);
645
646
647 }
648 lastAdjointCSVsModelsAndElements.push_back(CS_Vs_Element);
649
650 }
651
652
653 G4double CS=0;
654 for (size_t i=0;i<CS_Vs_Element.size();i++){
655 CS+=CS_Vs_Element[i]; //We could put the progressive sum of the CS instead of the CS of an element itself
656
657 }
658 return CS;
659}
660///////////////////////////////////////////////////////
661//
662G4Element* G4AdjointCSManager::SampleElementFromCSMatrices(G4Material* aMaterial,
663 G4VEmAdjointModel* aModel,
664 G4double PrimEnergy,
665 G4double Tcut,
666 G4bool IsScatProjToProjCase)
667{ std::vector<G4double> CS_Vs_Element;
668 G4double CS = ComputeAdjointCS(aMaterial,aModel,PrimEnergy,Tcut,IsScatProjToProjCase,CS_Vs_Element);
669 G4double rand_var= G4UniformRand();
670 G4double SumCS=0.;
671 size_t ind=0;
672 for (size_t i=0;i<CS_Vs_Element.size();i++){
673 SumCS+=CS_Vs_Element[i];
674 if (rand_var<=SumCS/CS){
675 ind=i;
676 break;
677 }
678 }
679
680 return const_cast<G4Element*>(aMaterial->GetElement(ind));
681
682
683
684}
685///////////////////////////////////////////////////////
686//
687G4double G4AdjointCSManager::ComputeTotalAdjointCS(const G4MaterialCutsCouple* aCouple,
688 G4ParticleDefinition* aPartDef,
689 G4double Ekin)
690{
691 G4double TotalCS=0.;
692
693 DefineCurrentMaterial(aCouple);
694
695
696 std::vector<G4double> CS_Vs_Element;
697 G4double CS;
698 for (size_t i=0; i<listOfAdjointEMModel.size();i++){
699
700 G4double Tlow=0;
701 if (!listOfAdjointEMModel[i]->GetApplyCutInRange()) Tlow =listOfAdjointEMModel[i]->GetLowEnergyLimit();
702 else {
703 G4ParticleDefinition* theDirSecondPartDef =
704 GetForwardParticleEquivalent(listOfAdjointEMModel[i]->GetAdjointEquivalentOfDirectSecondaryParticleDefinition());
705 size_t idx=56;
706 if (theDirSecondPartDef->GetParticleName() == "gamma") idx = 0;
707 else if (theDirSecondPartDef->GetParticleName() == "e-") idx = 1;
708 else if (theDirSecondPartDef->GetParticleName() == "e+") idx = 2;
709 if (idx <56) {
710 const std::vector<G4double>* aVec = G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(idx);
711 Tlow =(*aVec)[aCouple->GetIndex()];
712 }
713
714
715 }
716 if ( Ekin<=listOfAdjointEMModel[i]->GetHighEnergyLimit() && Ekin>=listOfAdjointEMModel[i]->GetLowEnergyLimit()){
717 if (aPartDef == listOfAdjointEMModel[i]->GetAdjointEquivalentOfDirectPrimaryParticleDefinition()){
718 CS=ComputeAdjointCS(currentMaterial,
719 listOfAdjointEMModel[i],
720 Ekin, Tlow,true,CS_Vs_Element);
721 TotalCS += CS;
722 (*listSigmaTableForAdjointModelScatProjToProj[i])[currentMatIndex]->PutValue(eindex,CS);
723 }
724 if (aPartDef == listOfAdjointEMModel[i]->GetAdjointEquivalentOfDirectSecondaryParticleDefinition()){
725 CS = ComputeAdjointCS(currentMaterial,
726 listOfAdjointEMModel[i],
727 Ekin, Tlow,false, CS_Vs_Element);
728 TotalCS += CS;
729 (*listSigmaTableForAdjointModelProdToProj[i])[currentMatIndex]->PutValue(eindex,CS);
730 }
731
732 }
733 else {
734 (*listSigmaTableForAdjointModelScatProjToProj[i])[currentMatIndex]->PutValue(eindex,0.);
735 (*listSigmaTableForAdjointModelProdToProj[i])[currentMatIndex]->PutValue(eindex,0.);
736
737 }
738 }
739 return TotalCS;
740
741
742}
743///////////////////////////////////////////////////////
744//
745std::vector<G4AdjointCSMatrix*>
746G4AdjointCSManager::BuildCrossSectionsMatricesForAGivenModelAndElement(G4VEmAdjointModel* aModel,G4int Z,G4int A,
747 G4int nbin_pro_decade)
748{
749 G4AdjointCSMatrix* theCSMatForProdToProjBackwardScattering = new G4AdjointCSMatrix(false);
750 G4AdjointCSMatrix* theCSMatForScatProjToProjBackwardScattering = new G4AdjointCSMatrix(true);
751
752
753 //make the vector of primary energy of the adjoint particle, could try to make this just once ?
754
755 G4double EkinMin =aModel->GetLowEnergyLimit();
756 G4double EkinMaxForScat =aModel->GetHighEnergyLimit()*0.999;
757 G4double EkinMaxForProd =aModel->GetHighEnergyLimit()*0.999;
758 if (aModel->GetSecondPartOfSameType() )EkinMaxForProd =EkinMaxForProd/2.;
759
760
761 //Product to projectile backward scattering
762 //-----------------------------------------
763 G4double fE=std::pow(10.,1./nbin_pro_decade);
764 G4double E2=std::pow(10.,double( int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
765 G4double E1=EkinMin;
766 while (E1 <EkinMaxForProd){
767 E1=std::max(EkinMin,E2);
768 E1=std::min(EkinMaxForProd,E1);
769 std::vector< std::vector< double>* > aMat= aModel->ComputeAdjointCrossSectionVectorPerAtomForSecond(E1,Z,A,nbin_pro_decade);
770 if (aMat.size()>=2) {
771 std::vector< double>* log_ESecVec=aMat[0];
772 std::vector< double>* log_CSVec=aMat[1];
773 G4double log_adjointCS=log_CSVec->back();
774 //normalise CSVec such that it becomes a probability vector
775 for (size_t j=0;j<log_CSVec->size();j++) {
776 if (j==0) (*log_CSVec)[j] = 0.;
777 else (*log_CSVec)[j]=std::log(1.-std::exp((*log_CSVec)[j]-log_adjointCS) +1e-50);
778 }
779 (*log_CSVec)[log_CSVec->size()-1]=(*log_CSVec)[log_CSVec->size()-2]-std::log(1000.);
780 theCSMatForProdToProjBackwardScattering->AddData(std::log(E1),log_adjointCS,log_ESecVec,log_CSVec,0);
781 }
782 E1=E2;
783 E2*=fE;
784 }
785
786 //Scattered projectile to projectile backward scattering
787 //-----------------------------------------
788
789 E2=std::pow(10.,double( int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
790 E1=EkinMin;
791 while (E1 <EkinMaxForScat){
792 E1=std::max(EkinMin,E2);
793 E1=std::min(EkinMaxForScat,E1);
794 std::vector< std::vector< double>* > aMat= aModel->ComputeAdjointCrossSectionVectorPerAtomForScatProj(E1,Z,A,nbin_pro_decade);
795 if (aMat.size()>=2) {
796 std::vector< double>* log_ESecVec=aMat[0];
797 std::vector< double>* log_CSVec=aMat[1];
798 G4double log_adjointCS=log_CSVec->back();
799 //normalise CSVec such that it becomes a probability vector
800 for (size_t j=0;j<log_CSVec->size();j++) {
801 if (j==0) (*log_CSVec)[j] = 0.;
802 else (*log_CSVec)[j]=std::log(1.-std::exp((*log_CSVec)[j]-log_adjointCS)+1e-50);
803 }
804 (*log_CSVec)[log_CSVec->size()-1]=(*log_CSVec)[log_CSVec->size()-2]-std::log(1000.);
805 theCSMatForScatProjToProjBackwardScattering->AddData(std::log(E1),log_adjointCS,log_ESecVec,log_CSVec,0);
806 }
807 E1=E2;
808 E2*=fE;
809 }
810
811
812 std::vector<G4AdjointCSMatrix*> res;
813 res.clear();
814 res.push_back(theCSMatForProdToProjBackwardScattering);
815 res.push_back(theCSMatForScatProjToProjBackwardScattering);
816
817
818/*
819 G4String file_name;
820 std::stringstream astream;
821 G4String str_Z;
822 astream<<Z;
823 astream>>str_Z;
824 theCSMatForProdToProjBackwardScattering->Write(aModel->GetName()+G4String("_CSMat_Z")+str_Z+"_ProdToProj.txt");
825 theCSMatForScatProjToProjBackwardScattering->Write(aModel->GetName()+G4String("_CSMat_Z")+str_Z+"_ScatProjToProj.txt");
826
827*/
828
829
830 return res;
831
832
833}
834///////////////////////////////////////////////////////
835//
836std::vector<G4AdjointCSMatrix*>
837G4AdjointCSManager::BuildCrossSectionsMatricesForAGivenModelAndMaterial(G4VEmAdjointModel* aModel,
838 G4Material* aMaterial,
839 G4int nbin_pro_decade)
840{
841 G4AdjointCSMatrix* theCSMatForProdToProjBackwardScattering = new G4AdjointCSMatrix(false);
842 G4AdjointCSMatrix* theCSMatForScatProjToProjBackwardScattering = new G4AdjointCSMatrix(true);
843
844
845 //make the vector of primary energy of the adjoint particle, could try to make this just once ?
846
847 G4double EkinMin =aModel->GetLowEnergyLimit();
848 G4double EkinMaxForScat =aModel->GetHighEnergyLimit()*0.999;
849 G4double EkinMaxForProd =aModel->GetHighEnergyLimit()*0.999;
850 if (aModel->GetSecondPartOfSameType() )EkinMaxForProd =EkinMaxForProd/2.;
851
852
853
854
855
856
857
858 //Product to projectile backward scattering
859 //-----------------------------------------
860 G4double fE=std::pow(10.,1./nbin_pro_decade);
861 G4double E2=std::pow(10.,double( int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
862 G4double E1=EkinMin;
863 while (E1 <EkinMaxForProd){
864 E1=std::max(EkinMin,E2);
865 E1=std::min(EkinMaxForProd,E1);
866 std::vector< std::vector< double>* > aMat= aModel->ComputeAdjointCrossSectionVectorPerVolumeForSecond(aMaterial,E1,nbin_pro_decade);
867 if (aMat.size()>=2) {
868 std::vector< double>* log_ESecVec=aMat[0];
869 std::vector< double>* log_CSVec=aMat[1];
870 G4double log_adjointCS=log_CSVec->back();
871
872 //normalise CSVec such that it becomes a probability vector
873 for (size_t j=0;j<log_CSVec->size();j++) {
874 //G4cout<<"CSMan1 "<<(*log_CSVec)[j]<<G4endl;
875 if (j==0) (*log_CSVec)[j] = 0.;
876 else (*log_CSVec)[j]=std::log(1.-std::exp((*log_CSVec)[j]-log_adjointCS));
877 //G4cout<<"CSMan2 "<<(*log_CSVec)[j]<<G4endl;
878 }
879 (*log_CSVec)[log_CSVec->size()-1]=(*log_CSVec)[log_CSVec->size()-2]-std::log(1000.);
880 theCSMatForProdToProjBackwardScattering->AddData(std::log(E1),log_adjointCS,log_ESecVec,log_CSVec,0);
881 }
882
883
884
885 E1=E2;
886 E2*=fE;
887 }
888
889 //Scattered projectile to projectile backward scattering
890 //-----------------------------------------
891
892 E2=std::pow(10.,double( int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
893 E1=EkinMin;
894 while (E1 <EkinMaxForScat){
895 E1=std::max(EkinMin,E2);
896 E1=std::min(EkinMaxForScat,E1);
897 std::vector< std::vector< double>* > aMat= aModel->ComputeAdjointCrossSectionVectorPerVolumeForScatProj(aMaterial,E1,nbin_pro_decade);
898 if (aMat.size()>=2) {
899 std::vector< double>* log_ESecVec=aMat[0];
900 std::vector< double>* log_CSVec=aMat[1];
901 G4double log_adjointCS=log_CSVec->back();
902
903 for (size_t j=0;j<log_CSVec->size();j++) {
904 //G4cout<<"CSMan1 "<<(*log_CSVec)[j]<<G4endl;
905 if (j==0) (*log_CSVec)[j] = 0.;
906 else (*log_CSVec)[j]=std::log(1.-std::exp((*log_CSVec)[j]-log_adjointCS));
907 //G4cout<<"CSMan2 "<<(*log_CSVec)[j]<<G4endl;if (theAdjPartDef->GetParticleName() == "adj_gamma") return G4Gamma::Gamma();
908
909 }
910 (*log_CSVec)[log_CSVec->size()-1]=(*log_CSVec)[log_CSVec->size()-2]-std::log(1000.);
911
912 theCSMatForScatProjToProjBackwardScattering->AddData(std::log(E1),log_adjointCS,log_ESecVec,log_CSVec,0);
913 }
914 E1=E2;
915 E2*=fE;
916 }
917
918
919
920
921
922
923
924 std::vector<G4AdjointCSMatrix*> res;
925 res.clear();
926
927 res.push_back(theCSMatForProdToProjBackwardScattering);
928 res.push_back(theCSMatForScatProjToProjBackwardScattering);
929
930 /*
931 theCSMatForProdToProjBackwardScattering->Write(aModel->GetName()+"_CSMat_"+aMaterial->GetName()+"_ProdToProj.txt");
932 theCSMatForScatProjToProjBackwardScattering->Write(aModel->GetName()+"_CSMat_"+aMaterial->GetName()+"_ScatProjToProj.txt");
933*/
934
935
936 return res;
937
938
939}
940
941///////////////////////////////////////////////////////
942//
943G4ParticleDefinition* G4AdjointCSManager::GetAdjointParticleEquivalent(G4ParticleDefinition* theFwdPartDef)
944{
945 if (theFwdPartDef->GetParticleName() == "e-") return G4AdjointElectron::AdjointElectron();
946 else if (theFwdPartDef->GetParticleName() == "gamma") return G4AdjointGamma::AdjointGamma();
947 else if (theFwdPartDef->GetParticleName() == "proton") return G4AdjointProton::AdjointProton();
948 else if (theFwdPartDef ==theFwdIon) return theAdjIon;
949
950 return 0;
951}
952///////////////////////////////////////////////////////
953//
954G4ParticleDefinition* G4AdjointCSManager::GetForwardParticleEquivalent(G4ParticleDefinition* theAdjPartDef)
955{
956 if (theAdjPartDef->GetParticleName() == "adj_e-") return G4Electron::Electron();
957 else if (theAdjPartDef->GetParticleName() == "adj_gamma") return G4Gamma::Gamma();
958 else if (theAdjPartDef->GetParticleName() == "adj_proton") return G4Proton::Proton();
959 else if (theAdjPartDef == theAdjIon) return theFwdIon;
960 return 0;
961}
962///////////////////////////////////////////////////////
963//
964void G4AdjointCSManager::DefineCurrentMaterial(const G4MaterialCutsCouple* couple)
965{
966 if(couple != currentCouple) {
967 currentCouple = const_cast<G4MaterialCutsCouple*> (couple);
968 currentMaterial = const_cast<G4Material*> (couple->GetMaterial());
969 currentMatIndex = couple->GetIndex();
970 lastPartDefForCS =0;
971 LastEkinForCS =0;
972 LastCSCorrectionFactor =1.;
973 }
974}
975
976///////////////////////////////////////////////////////
977//
978void G4AdjointCSManager::DefineCurrentParticle(const G4ParticleDefinition* aPartDef)
979{
980 if(aPartDef != currentParticleDef) {
981
982 currentParticleDef= const_cast< G4ParticleDefinition* > (aPartDef);
983 massRatio=1;
984 if (aPartDef == theAdjIon) massRatio = proton_mass_c2/aPartDef->GetPDGMass();
985 currentParticleIndex=1000000;
986 for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){
987 if (aPartDef == theListOfAdjointParticlesInAction[i]) currentParticleIndex=i;
988 }
989
990 }
991}
992
993
994
995/////////////////////////////////////////////////////////////////////////////////////////////////
996//
997G4double G4AdjointCSManager::ComputeAdjointCS(G4double aPrimEnergy,G4AdjointCSMatrix*
998 anAdjointCSMatrix,G4double Tcut)
999{
1000 std::vector< double> *theLogPrimEnergyVector = anAdjointCSMatrix->GetLogPrimEnergyVector();
1001 if (theLogPrimEnergyVector->size() ==0){
1002 G4cout<<"No data are contained in the given AdjointCSMatrix!"<<G4endl;
1003 G4cout<<"The s"<<G4endl;
1004 return 0.;
1005
1006 }
1007 G4double log_Tcut = std::log(Tcut);
1008 G4double log_E =std::log(aPrimEnergy);
1009
1010 if (aPrimEnergy <= Tcut || log_E > theLogPrimEnergyVector->back()) return 0.;
1011
1012
1013
1014 G4AdjointInterpolator* theInterpolator=G4AdjointInterpolator::GetInstance();
1015
1016 size_t ind =theInterpolator->FindPositionForLogVector(log_E,*theLogPrimEnergyVector);
1017 G4double aLogPrimEnergy1,aLogPrimEnergy2;
1018 G4double aLogCS1,aLogCS2;
1019 G4double log01,log02;
1020 std::vector< double>* aLogSecondEnergyVector1 =0;
1021 std::vector< double>* aLogSecondEnergyVector2 =0;
1022 std::vector< double>* aLogProbVector1=0;
1023 std::vector< double>* aLogProbVector2=0;
1024 std::vector< size_t>* aLogProbVectorIndex1=0;
1025 std::vector< size_t>* aLogProbVectorIndex2=0;
1026
1027
1028 anAdjointCSMatrix->GetData(ind, aLogPrimEnergy1,aLogCS1,log01, aLogSecondEnergyVector1,aLogProbVector1,aLogProbVectorIndex1);
1029 anAdjointCSMatrix->GetData(ind+1, aLogPrimEnergy2,aLogCS2,log02, aLogSecondEnergyVector2,aLogProbVector2,aLogProbVectorIndex2);
1030 if (anAdjointCSMatrix->IsScatProjToProjCase()){ //case where the Tcut plays a role
1031 G4double log_minimum_prob1, log_minimum_prob2;
1032 log_minimum_prob1=theInterpolator->InterpolateForLogVector(log_Tcut,*aLogSecondEnergyVector1,*aLogProbVector1);
1033 log_minimum_prob2=theInterpolator->InterpolateForLogVector(log_Tcut,*aLogSecondEnergyVector2,*aLogProbVector2);
1034 aLogCS1+= log_minimum_prob1;
1035 aLogCS2+= log_minimum_prob2;
1036 }
1037
1038 G4double log_adjointCS = theInterpolator->LinearInterpolation(log_E,aLogPrimEnergy1,aLogPrimEnergy2,aLogCS1,aLogCS2);
1039 return std::exp(log_adjointCS);
1040
1041
1042}
G4ElementTable
std::vector< G4Element * > G4ElementTable
Definition: G4ElementTable.hh:43
G4MaterialTable
std::vector< G4Material * > G4MaterialTable
Definition: G4MaterialTable.hh:43
G4double
double G4double
Definition: G4Types.hh:64
G4int
int G4int
Definition: G4Types.hh:66
G4bool
bool G4bool
Definition: G4Types.hh:67
G4endl
#define G4endl
Definition: G4ios.hh:52
G4cout
G4DLLIMPORT std::ostream G4cout
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:53
G4AdjointCSManager::RegisterEmAdjointModel
size_t RegisterEmAdjointModel(G4VEmAdjointModel *)
Definition: G4AdjointCSManager.cc:121
G4AdjointCSManager::RegisterAdjointParticle
void RegisterAdjointParticle(G4ParticleDefinition *aPartDef)
Definition: G4AdjointCSManager.cc:159
G4AdjointCSManager::GetForwardParticleEquivalent
G4ParticleDefinition * GetForwardParticleEquivalent(G4ParticleDefinition *theAdjPartDef)
Definition: G4AdjointCSManager.cc:954
G4AdjointCSManager::GetMaxAdjTotalCS
void GetMaxAdjTotalCS(G4ParticleDefinition *aPartDef, const G4MaterialCutsCouple *aCouple, G4double &e_sigma_max, G4double &sigma_max)
Definition: G4AdjointCSManager.cc:452
G4AdjointCSManager::GetCrossSectionCorrection
G4double GetCrossSectionCorrection(G4ParticleDefinition *aPartDef, G4double PreStepEkin, const G4MaterialCutsCouple *aCouple, G4bool &fwd_is_used, G4double &fwd_TotCS)
Definition: G4AdjointCSManager.cc:465
G4AdjointCSManager::GetTotalForwardCS
G4double GetTotalForwardCS(G4ParticleDefinition *aPartDef, G4double Ekin, const G4MaterialCutsCouple *aCouple)
Definition: G4AdjointCSManager.cc:407
G4AdjointCSManager::ComputeTotalAdjointCS
G4double ComputeTotalAdjointCS(const G4MaterialCutsCouple *aMatCutCouple, G4ParticleDefinition *aPart, G4double PrimEnergy)
Definition: G4AdjointCSManager.cc:687
G4AdjointCSManager::GetAdjointParticleEquivalent
G4ParticleDefinition * GetAdjointParticleEquivalent(G4ParticleDefinition *theFwdPartDef)
Definition: G4AdjointCSManager.cc:943
G4AdjointCSManager::SampleElementFromCSMatrices
G4Element * SampleElementFromCSMatrices(G4Material *aMaterial, G4VEmAdjointModel *aModel, G4double PrimEnergy, G4double Tcut, G4bool IsScatProjToProjCase)
Definition: G4AdjointCSManager.cc:662
G4AdjointCSManager::GetContinuousWeightCorrection
G4double GetContinuousWeightCorrection(G4ParticleDefinition *aPartDef, G4double PreStepEkin, G4double AfterStepEkin, const G4MaterialCutsCouple *aCouple, G4double step_length)
Definition: G4AdjointCSManager.cc:502
G4AdjointCSManager::GetPostStepWeightCorrection
G4double GetPostStepWeightCorrection()
Definition: G4AdjointCSManager.cc:525
G4AdjointCSManager::GetEminForTotalCS
void GetEminForTotalCS(G4ParticleDefinition *aPartDef, const G4MaterialCutsCouple *aCouple, G4double &emin_adj, G4double &emin_fwd)
Definition: G4AdjointCSManager.cc:427
G4AdjointCSManager::GetAdjointCSManager
static G4AdjointCSManager * GetAdjointCSManager()
Definition: G4AdjointCSManager.cc:61
G4AdjointCSManager::GetTotalAdjointCS
G4double GetTotalAdjointCS(G4ParticleDefinition *aPartDef, G4double Ekin, const G4MaterialCutsCouple *aCouple)
Definition: G4AdjointCSManager.cc:395
G4AdjointCSManager::GetAdjointSigma
G4double GetAdjointSigma(G4double Ekin_nuc, size_t index_model, G4bool is_scat_proj_to_proj, const G4MaterialCutsCouple *aCouple)
Definition: G4AdjointCSManager.cc:418
G4AdjointCSManager::RegisterEmProcess
void RegisterEmProcess(G4VEmProcess *aProcess, G4ParticleDefinition *aPartDef)
Definition: G4AdjointCSManager.cc:130
G4AdjointCSManager::RegisterEnergyLossProcess
void RegisterEnergyLossProcess(G4VEnergyLossProcess *aProcess, G4ParticleDefinition *aPartDef)
Definition: G4AdjointCSManager.cc:145
G4AdjointCSManager::GetMaxFwdTotalCS
void GetMaxFwdTotalCS(G4ParticleDefinition *aPartDef, const G4MaterialCutsCouple *aCouple, G4double &e_sigma_max, G4double &sigma_max)
Definition: G4AdjointCSManager.cc:439
G4AdjointCSManager::ComputeAdjointCS
G4double ComputeAdjointCS(G4Material *aMaterial, G4VEmAdjointModel *aModel, G4double PrimEnergy, G4double Tcut, G4bool IsScatProjToProjCase, std::vector< G4double > &AdjointCS_for_each_element)
Definition: G4AdjointCSManager.cc:532
G4AdjointCSMatrix::GetData
G4bool GetData(unsigned int i, G4double &aPrimEnergy, G4double &aCS, G4double &log0, std::vector< double > *&aLogSecondEnergyVector, std::vector< double > *&aLogProbVector, std::vector< size_t > *&aLogProbVectorIndex)
Definition: G4AdjointCSMatrix.cc:108
G4AdjointCSMatrix::GetLogPrimEnergyVector
std::vector< double > * GetLogPrimEnergyVector()
Definition: G4AdjointCSMatrix.hh:80
G4AdjointCSMatrix::AddData
void AddData(G4double aPrimEnergy, G4double aCS, std::vector< double > *aLogSecondEnergyVector, std::vector< double > *aLogProbVector, size_t n_pro_decade=0)
Definition: G4AdjointCSMatrix.cc:71
G4AdjointElectron::AdjointElectron
static G4AdjointElectron * AdjointElectron()
Definition: G4AdjointElectron.cc:100
G4AdjointGamma::AdjointGamma
static G4AdjointGamma * AdjointGamma()
Definition: G4AdjointGamma.cc:84
G4AdjointInterpolator::LinearInterpolation
G4double LinearInterpolation(G4double &x, G4double &x1, G4double &x2, G4double &y1, G4double &y2)
Definition: G4AdjointInterpolator.cc:63
G4AdjointInterpolator::GetInstance
static G4AdjointInterpolator * GetInstance()
Definition: G4AdjointInterpolator.cc:43
G4AdjointInterpolator::FindPositionForLogVector
size_t FindPositionForLogVector(G4double &x, std::vector< G4double > &x_vec)
Definition: G4AdjointInterpolator.cc:153
G4AdjointInterpolator::InterpolateForLogVector
G4double InterpolateForLogVector(G4double &x, std::vector< G4double > &x_vec, std::vector< G4double > &y_vec)
Definition: G4AdjointInterpolator.cc:208
G4AdjointProton::AdjointProton
static G4AdjointProton * AdjointProton()
Definition: G4AdjointProton.cc:84
G4Electron::Electron
static G4Electron * Electron()
Definition: G4Electron.cc:94
G4Element::GetZ
G4double GetZ() const
Definition: G4Element.hh:131
G4Element::GetA
G4double GetA() const
Definition: G4Element.hh:138
G4Element::GetIndex
size_t GetIndex() const
Definition: G4Element.hh:182
G4Element::GetElementTable
static const G4ElementTable * GetElementTable()
Definition: G4Element.cc:399
G4Gamma::Gamma
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
G4MaterialCutsCouple::GetMaterial
const G4Material * GetMaterial() const
Definition: G4MaterialCutsCouple.hh:150
G4Material::GetMaterialTable
static const G4MaterialTable * GetMaterialTable()
Definition: G4Material.cc:562
G4Material::GetElement
const G4Element * GetElement(G4int iel) const
Definition: G4Material.hh:201
G4Material::GetNumberOfElements
size_t GetNumberOfElements() const
Definition: G4Material.hh:185
G4Material::GetVecNbOfAtomsPerVolume
const G4double * GetVecNbOfAtomsPerVolume() const
Definition: G4Material.hh:205
G4Material::GetIndex
size_t GetIndex() const
Definition: G4Material.hh:261
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition: G4ParticleDefinition.hh:115
G4PhysicsVector::GetVectorLength
size_t GetVectorLength() const
G4PhysicsVector::GetLowEdgeEnergy
virtual G4double GetLowEdgeEnergy(size_t binNumber) const
Definition: G4PhysicsVector.cc:158
G4PhysicsVector::PutValue
void PutValue(size_t index, G4double theValue)
G4ProductionCutsTable::GetMaterialCutsCouple
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
Definition: G4ProductionCutsTable.hh:262
G4ProductionCutsTable::GetEnergyCutsVector
const std::vector< G4double > * GetEnergyCutsVector(size_t pcIdx) const
Definition: G4ProductionCutsTable.hh:250
G4ProductionCutsTable::GetProductionCutsTable
static G4ProductionCutsTable * GetProductionCutsTable()
Definition: G4ProductionCutsTable.cc:63
G4Proton::Proton
static G4Proton * Proton()
Definition: G4Proton.cc:93
G4VEmAdjointModel::GetSecondPartOfSameType
G4bool GetSecondPartOfSameType()
Definition: G4VEmAdjointModel.hh:203
G4VEmAdjointModel::SetCSMatrices
void SetCSMatrices(std::vector< G4AdjointCSMatrix * > *Vec1CSMatrix, std::vector< G4AdjointCSMatrix * > *Vec2CSMatrix)
Definition: G4VEmAdjointModel.hh:174
G4VEmAdjointModel::GetUseOnlyOneMatrixForAllElements
G4bool GetUseOnlyOneMatrixForAllElements()
Definition: G4VEmAdjointModel.hh:213
G4VEmAdjointModel::GetSecondAdjEnergyMaxForProdToProjCase
virtual G4double GetSecondAdjEnergyMaxForProdToProjCase(G4double PrimAdjEnergy)
Definition: G4VEmAdjointModel.cc:671
G4VEmAdjointModel::ComputeAdjointCrossSectionVectorPerAtomForSecond
std::vector< std::vector< double > * > ComputeAdjointCrossSectionVectorPerAtomForSecond(G4double kinEnergyProd, G4double Z, G4double A=0., G4int nbin_pro_decade=10)
Definition: G4VEmAdjointModel.cc:238
G4VEmAdjointModel::GetLowEnergyLimit
G4double GetLowEnergyLimit()
Definition: G4VEmAdjointModel.hh:187
G4VEmAdjointModel::GetSecondAdjEnergyMaxForScatProjToProjCase
virtual G4double GetSecondAdjEnergyMaxForScatProjToProjCase(G4double PrimAdjEnergy)
Definition: G4VEmAdjointModel.cc:657
G4VEmAdjointModel::ComputeAdjointCrossSectionVectorPerVolumeForSecond
std::vector< std::vector< double > * > ComputeAdjointCrossSectionVectorPerVolumeForSecond(G4Material *aMaterial, G4double kinEnergyProd, G4int nbin_pro_decade=10)
Definition: G4VEmAdjointModel.cc:348
G4VEmAdjointModel::GetSecondAdjEnergyMinForProdToProjCase
virtual G4double GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy)
Definition: G4VEmAdjointModel.cc:676
G4VEmAdjointModel::ComputeAdjointCrossSectionVectorPerVolumeForScatProj
std::vector< std::vector< double > * > ComputeAdjointCrossSectionVectorPerVolumeForScatProj(G4Material *aMaterial, G4double kinEnergyProd, G4int nbin_pro_decade=10)
Definition: G4VEmAdjointModel.cc:399
G4VEmAdjointModel::AdjointCrossSection
virtual G4double AdjointCrossSection(const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool IsScatProjToProjCase)
Definition: G4VEmAdjointModel.cc:58
G4VEmAdjointModel::ComputeAdjointCrossSectionVectorPerAtomForScatProj
std::vector< std::vector< double > * > ComputeAdjointCrossSectionVectorPerAtomForScatProj(G4double kinEnergyProd, G4double Z, G4double A=0., G4int nbin_pro_decade=10)
Definition: G4VEmAdjointModel.cc:291
G4VEmAdjointModel::GetHighEnergyLimit
G4double GetHighEnergyLimit()
Definition: G4VEmAdjointModel.hh:185
G4VEmAdjointModel::GetSecondAdjEnergyMinForScatProjToProjCase
virtual G4double GetSecondAdjEnergyMinForScatProjToProjCase(G4double PrimAdjEnergy, G4double Tcut=0)
Definition: G4VEmAdjointModel.cc:664
G4VEmModel::GetChargeSquareRatio
virtual G4double GetChargeSquareRatio(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
Definition: G4VEmModel.cc:262