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