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