Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
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