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

Last change on this file since 1192 was 966, checked in by garnier, 17 years ago

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