source: trunk/source/processes/electromagnetic/adjoint/src/G4AdjointCSManager.cc@ 1197

Last change on this file since 1197 was 1196, checked in by garnier, 16 years ago

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