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