Geant4.10
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
58 G4ThreadLocal G4AdjointCSManager* G4AdjointCSManager::theInstance = 0;
59 ///////////////////////////////////////////////////////
60 //
61 G4AdjointCSManager* G4AdjointCSManager::GetAdjointCSManager()
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 
69 ///////////////////////////////////////////////////////
70 //
71 G4AdjointCSManager::G4AdjointCSManager()
72 { CrossSectionMatrixesAreBuilt=false;
73  TotalSigmaTableAreBuilt=false;
74  theTotalForwardSigmaTableVector.clear();
75  theTotalAdjointSigmaTableVector.clear();
76  listOfForwardEmProcess.clear();
77  listOfForwardEnergyLossProcess.clear();
78  theListOfAdjointParticlesInAction.clear();
79  EminForFwdSigmaTables.clear();
80  EminForAdjSigmaTables.clear();
81  EkinofFwdSigmaMax.clear();
82  EkinofAdjSigmaMax.clear();
83  listSigmaTableForAdjointModelScatProjToProj.clear();
84  listSigmaTableForAdjointModelProdToProj.clear();
85  Tmin=0.1*keV;
86  Tmax=100.*TeV;
87  nbins=320; //probably this should be decrease, that was choosen to avoid error in the CS value closed to CS jump.(For example at Tcut)
88 
89  RegisterAdjointParticle(G4AdjointElectron::AdjointElectron());
90  RegisterAdjointParticle(G4AdjointGamma::AdjointGamma());
91  RegisterAdjointParticle(G4AdjointProton::AdjointProton());
92 
93  verbose = 1;
94 
95  lastPartDefForCS =0;
96  LastEkinForCS =0;
97  LastCSCorrectionFactor =1.;
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 }
114 ///////////////////////////////////////////////////////
115 //
116 G4AdjointCSManager::~G4AdjointCSManager()
117 {;
118 }
119 ///////////////////////////////////////////////////////
120 //
121 size_t G4AdjointCSManager::RegisterEmAdjointModel(G4VEmAdjointModel* aModel)
122 {listOfAdjointEMModel.push_back(aModel);
123  listSigmaTableForAdjointModelScatProjToProj.push_back(new G4PhysicsTable);
124  listSigmaTableForAdjointModelProdToProj.push_back(new G4PhysicsTable);
125  return listOfAdjointEMModel.size() -1;
126 
127 }
128 ///////////////////////////////////////////////////////
129 //
130 void G4AdjointCSManager::RegisterEmProcess(G4VEmProcess* aProcess, G4ParticleDefinition* aFwdPartDef)
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 }
143 ///////////////////////////////////////////////////////
144 //
145 void G4AdjointCSManager::RegisterEnergyLossProcess(G4VEnergyLossProcess* aProcess, G4ParticleDefinition* aFwdPartDef)
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 }
157 ///////////////////////////////////////////////////////
158 //
159 void G4AdjointCSManager::RegisterAdjointParticle(G4ParticleDefinition* aPartDef)
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*>());
167  theTotalForwardSigmaTableVector.push_back(new G4PhysicsTable);
168  theTotalAdjointSigmaTableVector.push_back(new G4PhysicsTable);
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 }
178 ///////////////////////////////////////////////////////
179 //
180 void G4AdjointCSManager::BuildCrossSectionMatrices()
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 
202  theAdjointCSMatricesForScatProjToProj.clear();
203  theAdjointCSMatricesForProdToProj.clear();
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++){
209  G4VEmAdjointModel* aModel =listOfAdjointEMModel[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 
259  CrossSectionMatrixesAreBuilt = true;
260 
261 
262 }
263 
264 
265 ///////////////////////////////////////////////////////
266 //
267 void G4AdjointCSManager::BuildTotalSigmaTables()
268 { if (TotalSigmaTableAreBuilt) return;
269 
270 
271  const G4ProductionCutsTable* theCoupleTable= G4ProductionCutsTable::GetProductionCutsTable();
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++){
279  listSigmaTableForAdjointModelScatProjToProj[i]->push_back(new G4PhysicsLogVector(Tmin, Tmax, nbins));
280  listSigmaTableForAdjointModelProdToProj[i]->push_back(new G4PhysicsLogVector(Tmin, Tmax, nbins));
281  }
282  }
283 
284 
285 
286  for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){
287  G4ParticleDefinition* thePartDef = theListOfAdjointParticlesInAction[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
317  G4PhysicsVector* aVector = new G4PhysicsLogVector(Tmin, Tmax, nbins);
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  }
334  G4double e1=e/massRatio;
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  }
390  TotalSigmaTableAreBuilt =true;
391 
392 }
393 ///////////////////////////////////////////////////////
394 //
395 G4double G4AdjointCSManager::GetTotalAdjointCS(G4ParticleDefinition* aPartDef, G4double Ekin,
396  const G4MaterialCutsCouple* aCouple)
397 { DefineCurrentMaterial(aCouple);
398  DefineCurrentParticle(aPartDef);
399  G4bool b;
400  return (((*theTotalAdjointSigmaTableVector[currentParticleIndex])[currentMatIndex])->GetValue(Ekin*massRatio, b));
401 
402 
403 
404 }
405 ///////////////////////////////////////////////////////
406 //
407 G4double G4AdjointCSManager::GetTotalForwardCS(G4ParticleDefinition* aPartDef, G4double Ekin,
408  const G4MaterialCutsCouple* aCouple)
409 { DefineCurrentMaterial(aCouple);
410  DefineCurrentParticle(aPartDef);
411  G4bool b;
412  return (((*theTotalForwardSigmaTableVector[currentParticleIndex])[currentMatIndex])->GetValue(Ekin*massRatio, b));
413 
414 
415 }
416 ///////////////////////////////////////////////////////
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 }
425 ///////////////////////////////////////////////////////
426 //
427 void G4AdjointCSManager::GetEminForTotalCS(G4ParticleDefinition* aPartDef,
428  const G4MaterialCutsCouple* aCouple, G4double& emin_adj, G4double& emin_fwd)
429 { DefineCurrentMaterial(aCouple);
430  DefineCurrentParticle(aPartDef);
431  emin_adj = EminForAdjSigmaTables[currentParticleIndex][currentMatIndex]/massRatio;
432  emin_fwd = EminForFwdSigmaTables[currentParticleIndex][currentMatIndex]/massRatio;
433 
434 
435 
436 }
437 ///////////////////////////////////////////////////////
438 //
439 void G4AdjointCSManager::GetMaxFwdTotalCS(G4ParticleDefinition* aPartDef,
440  const G4MaterialCutsCouple* aCouple, G4double& e_sigma_max, G4double& sigma_max)
441 { DefineCurrentMaterial(aCouple);
442  DefineCurrentParticle(aPartDef);
443  e_sigma_max = EkinofFwdSigmaMax[currentParticleIndex][currentMatIndex];
444  G4bool b;
445  sigma_max =((*theTotalForwardSigmaTableVector[currentParticleIndex])[currentMatIndex])->GetValue(e_sigma_max, b);
446  e_sigma_max/=massRatio;
447 
448 
449 }
450 ///////////////////////////////////////////////////////
451 //
452 void G4AdjointCSManager::GetMaxAdjTotalCS(G4ParticleDefinition* aPartDef,
453  const G4MaterialCutsCouple* aCouple, G4double& e_sigma_max, G4double& sigma_max)
454 { DefineCurrentMaterial(aCouple);
455  DefineCurrentParticle(aPartDef);
456  e_sigma_max = EkinofAdjSigmaMax[currentParticleIndex][currentMatIndex];
457  G4bool b;
458  sigma_max =((*theTotalAdjointSigmaTableVector[currentParticleIndex])[currentMatIndex])->GetValue(e_sigma_max, b);
459  e_sigma_max/=massRatio;
460 
461 
462 }
463 ///////////////////////////////////////////////////////
464 //
465 G4double G4AdjointCSManager::GetCrossSectionCorrection(G4ParticleDefinition* aPartDef,G4double PreStepEkin,const G4MaterialCutsCouple* aCouple, G4bool& fwd_is_used,
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;
478  LastCSCorrectionFactor = PrefwdCS/PreadjCS;
479  }
480  else {
481  forward_CS_is_used = false;
482  LastCSCorrectionFactor = 1.;
483 
484  }
485 
486  }
487  corr_fac =LastCSCorrectionFactor;
488 
489 
490 
491  }
492  else {
493  forward_CS_is_used = false;
494  LastCSCorrectionFactor = 1.;
495  }
496  fwd_TotCS=PrefwdCS;
497  fwd_is_used = forward_CS_is_used;
498  return corr_fac;
499 }
500 ///////////////////////////////////////////////////////
501 //
502 G4double G4AdjointCSManager::GetContinuousWeightCorrection(G4ParticleDefinition* aPartDef, G4double PreStepEkin,G4double AfterStepEkin,
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);
513  LastCSCorrectionFactor = 1.;
514  }
515  else {
516  LastCSCorrectionFactor = after_fwdCS/pre_adjCS;
517  }
518 
519 
520 
521  return corr_fac;
522 }
523 ///////////////////////////////////////////////////////
524 //
525 G4double G4AdjointCSManager::GetPostStepWeightCorrection( )
526 {//return 1.;
527  return 1./LastCSCorrectionFactor;
528 
529 }
530 ///////////////////////////////////////////////////////
531 //
532 G4double G4AdjointCSManager::ComputeAdjointCS(G4Material* aMaterial,
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;
559  listOfIndexOfAdjointEMModelInAction.clear();
560  listOfIsScatProjToProjCase.clear();
561  lastAdjointCSVsModelsAndElements.clear();
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++){
569  size_t ind1=listOfIndexOfAdjointEMModelInAction[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 }
660 ///////////////////////////////////////////////////////
661 //
662 G4Element* G4AdjointCSManager::SampleElementFromCSMatrices(G4Material* aMaterial,
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 }
685 ///////////////////////////////////////////////////////
686 //
687 G4double G4AdjointCSManager::ComputeTotalAdjointCS(const G4MaterialCutsCouple* aCouple,
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()){
718  CS=ComputeAdjointCS(currentMaterial,
719  listOfAdjointEMModel[i],
720  Ekin, Tlow,true,CS_Vs_Element);
721  TotalCS += CS;
722  (*listSigmaTableForAdjointModelScatProjToProj[i])[currentMatIndex]->PutValue(eindex,CS);
723  }
724  if (aPartDef == listOfAdjointEMModel[i]->GetAdjointEquivalentOfDirectSecondaryParticleDefinition()){
725  CS = ComputeAdjointCS(currentMaterial,
726  listOfAdjointEMModel[i],
727  Ekin, Tlow,false, CS_Vs_Element);
728  TotalCS += CS;
729  (*listSigmaTableForAdjointModelProdToProj[i])[currentMatIndex]->PutValue(eindex,CS);
730  }
731 
732  }
733  else {
734  (*listSigmaTableForAdjointModelScatProjToProj[i])[currentMatIndex]->PutValue(eindex,0.);
735  (*listSigmaTableForAdjointModelProdToProj[i])[currentMatIndex]->PutValue(eindex,0.);
736 
737  }
738  }
739  return TotalCS;
740 
741 
742 }
743 ///////////////////////////////////////////////////////
744 //
745 std::vector<G4AdjointCSMatrix*>
746 G4AdjointCSManager::BuildCrossSectionsMatricesForAGivenModelAndElement(G4VEmAdjointModel* aModel,G4int Z,G4int A,
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 }
834 ///////////////////////////////////////////////////////
835 //
836 std::vector<G4AdjointCSMatrix*>
837 G4AdjointCSManager::BuildCrossSectionsMatricesForAGivenModelAndMaterial(G4VEmAdjointModel* aModel,
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 
941 ///////////////////////////////////////////////////////
942 //
943 G4ParticleDefinition* G4AdjointCSManager::GetAdjointParticleEquivalent(G4ParticleDefinition* theFwdPartDef)
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 }
952 ///////////////////////////////////////////////////////
953 //
954 G4ParticleDefinition* G4AdjointCSManager::GetForwardParticleEquivalent(G4ParticleDefinition* theAdjPartDef)
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 }
962 ///////////////////////////////////////////////////////
963 //
964 void G4AdjointCSManager::DefineCurrentMaterial(const G4MaterialCutsCouple* couple)
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;
972  LastCSCorrectionFactor =1.;
973  }
974 }
975 
976 ///////////////////////////////////////////////////////
977 //
978 void G4AdjointCSManager::DefineCurrentParticle(const G4ParticleDefinition* aPartDef)
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++){
987  if (aPartDef == theListOfAdjointParticlesInAction[i]) currentParticleIndex=i;
988  }
989 
990  }
991 }
992 
993 
994 
995 /////////////////////////////////////////////////////////////////////////////////////////////////
996 //
997 G4double G4AdjointCSManager::ComputeAdjointCS(G4double aPrimEnergy,G4AdjointCSMatrix*
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 
1014  G4AdjointInterpolator* theInterpolator=G4AdjointInterpolator::GetInstance();
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 }
G4AdjointGamma::AdjointGamma
static G4AdjointGamma * AdjointGamma()
Definition: G4AdjointGamma.cc:84
G4AdjointCSManager::GetTotalAdjointCS
G4double GetTotalAdjointCS(G4ParticleDefinition *aPartDef, G4double Ekin, const G4MaterialCutsCouple *aCouple)
Definition: G4AdjointCSManager.cc:395
G4AdjointCSManager::GetCrossSectionCorrection
G4double GetCrossSectionCorrection(G4ParticleDefinition *aPartDef, G4double PreStepEkin, const G4MaterialCutsCouple *aCouple, G4bool &fwd_is_used, G4double &fwd_TotCS)
Definition: G4AdjointCSManager.cc:465
G4AdjointCSManager::GetAdjointSigma
G4double GetAdjointSigma(G4double Ekin_nuc, size_t index_model, G4bool is_scat_proj_to_proj, const G4MaterialCutsCouple *aCouple)
Definition: G4AdjointCSManager.cc:418
G4ProductionCutsTable::GetEnergyCutsVector
const std::vector< G4double > * GetEnergyCutsVector(size_t pcIdx) const
Definition: G4ProductionCutsTable.hh:250
G4AdjointCSManager::ComputeTotalAdjointCS
G4double ComputeTotalAdjointCS(const G4MaterialCutsCouple *aMatCutCouple, G4ParticleDefinition *aPart, G4double PrimEnergy)
Definition: G4AdjointCSManager.cc:687
G4VEmAdjointModel::GetSecondAdjEnergyMaxForProdToProjCase
virtual G4double GetSecondAdjEnergyMaxForProdToProjCase(G4double PrimAdjEnergy)
Definition: G4VEmAdjointModel.cc:671
G4AdjointCSMatrix::GetData
G4bool GetData(unsigned int i, G4double &aPrimEnergy, G4double &aCS, G4double &log0, std::vector< double > *&aLogSecondEnergyVector, std::vector< double > *&aLogProbVector, std::vector< size_t > *&aLogProbVectorIndex)
Definition: G4AdjointCSMatrix.cc:108
G4AdjointCSManager::GetPostStepWeightCorrection
G4double GetPostStepWeightCorrection()
Definition: G4AdjointCSManager.cc:525
G4Element.hh
int
typedef int(XMLCALL *XML_NotStandaloneHandler)(void *userData)
G4Material::GetIndex
size_t GetIndex() const
Definition: G4Material.hh:260
G4AdjointCSManager::GetMaxFwdTotalCS
void GetMaxFwdTotalCS(G4ParticleDefinition *aPartDef, const G4MaterialCutsCouple *aCouple, G4double &e_sigma_max, G4double &sigma_max)
Definition: G4AdjointCSManager.cc:439
G4AdjointCSManager::RegisterEmAdjointModel
size_t RegisterEmAdjointModel(G4VEmAdjointModel *)
Definition: G4AdjointCSManager.cc:121
G4Element::GetZ
G4double GetZ() const
Definition: G4Element.hh:131
G4Material::GetMaterialTable
static G4MaterialTable * GetMaterialTable()
Definition: G4Material.cc:564
G4MaterialTable
std::vector< G4Material * > G4MaterialTable
Definition: G4MaterialTable.hh:41
G4PhysicsLogVector.hh
G4AdjointElectron::AdjointElectron
static G4AdjointElectron * AdjointElectron()
Definition: G4AdjointElectron.cc:100
G4VEmAdjointModel::GetUseOnlyOneMatrixForAllElements
G4bool GetUseOnlyOneMatrixForAllElements()
Definition: G4VEmAdjointModel.hh:213
G4AdjointCSManager::RegisterEnergyLossProcess
void RegisterEnergyLossProcess(G4VEnergyLossProcess *aProcess, G4ParticleDefinition *aPartDef)
Definition: G4AdjointCSManager.cc:145
G4PhysicsVector::GetVectorLength
size_t GetVectorLength() const
G4PhysicsVector::GetLowEdgeEnergy
G4double GetLowEdgeEnergy(size_t binNumber) const
Definition: G4PhysicsVector.cc:157
G4ThreadLocal
#define G4ThreadLocal
Definition: tls.hh:52
G4Material::GetElement
const G4Element * GetElement(G4int iel) const
Definition: G4Material.hh:200
G4Element::GetA
G4double GetA() const
Definition: G4Element.hh:138
G4int
int G4int
Definition: G4Types.hh:78
G4AdjointCSMatrix::AddData
void AddData(G4double aPrimEnergy, G4double aCS, std::vector< double > *aLogSecondEnergyVector, std::vector< double > *aLogProbVector, size_t n_pro_decade=0)
Definition: G4AdjointCSMatrix.cc:71
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition: G4ParticleDefinition.hh:159
G4VEmAdjointModel::GetSecondAdjEnergyMaxForScatProjToProjCase
virtual G4double GetSecondAdjEnergyMaxForScatProjToProjCase(G4double PrimAdjEnergy)
Definition: G4VEmAdjointModel.cc:657
G4AdjointInterpolator::InterpolateForLogVector
G4double InterpolateForLogVector(G4double &x, std::vector< G4double > &x_vec, std::vector< G4double > &y_vec)
Definition: G4AdjointInterpolator.cc:218
G4AdjointCSManager::RegisterAdjointParticle
void RegisterAdjointParticle(G4ParticleDefinition *aPartDef)
Definition: G4AdjointCSManager.cc:159
G4AdjointInterpolator.hh
G4Material::GetVecNbOfAtomsPerVolume
const G4double * GetVecNbOfAtomsPerVolume() const
Definition: G4Material.hh:204
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:87
G4cout
G4GLOB_DLL std::ostream G4cout
G4ElementTable.hh
G4VEmAdjointModel::GetSecondAdjEnergyMinForProdToProjCase
virtual G4double GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy)
Definition: G4VEmAdjointModel.cc:676
G4AdjointCSManager::GetAdjointParticleEquivalent
G4ParticleDefinition * GetAdjointParticleEquivalent(G4ParticleDefinition *theFwdPartDef)
Definition: G4AdjointCSManager.cc:943
G4PhysicsTable.hh
G4bool
bool G4bool
Definition: G4Types.hh:79
G4AdjointCSMatrix::GetLogPrimEnergyVector
std::vector< double > * GetLogPrimEnergyVector()
Definition: G4AdjointCSMatrix.hh:80
G4VEmProcess.hh
G4Proton.hh
G4PhysicsVector::PutValue
void PutValue(size_t index, G4double theValue)
G4AdjointInterpolator::FindPositionForLogVector
size_t FindPositionForLogVector(G4double &x, std::vector< G4double > &x_vec)
Definition: G4AdjointInterpolator.cc:162
G4Element::GetIndex
size_t GetIndex() const
Definition: G4Element.hh:181
G4Proton::Proton
static G4Proton * Proton()
Definition: G4Proton.cc:93
python.hepunit.proton_mass_c2
float proton_mass_c2
Definition: hepunit.py:275
G4ParticleDefinition.hh
G4VEmAdjointModel::GetLowEnergyLimit
G4double GetLowEnergyLimit()
Definition: G4VEmAdjointModel.hh:187
G4AdjointProton.hh
G4VEmAdjointModel::SetCSMatrices
void SetCSMatrices(std::vector< G4AdjointCSMatrix * > *Vec1CSMatrix, std::vector< G4AdjointCSMatrix * > *Vec2CSMatrix)
Definition: G4VEmAdjointModel.hh:174
G4Gamma::Gamma
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
G4AdjointCSManager::GetTotalForwardCS
G4double GetTotalForwardCS(G4ParticleDefinition *aPartDef, G4double Ekin, const G4MaterialCutsCouple *aCouple)
Definition: G4AdjointCSManager.cc:407
G4VEnergyLossProcess.hh
G4PhysicalConstants.hh
G4AdjointCSManager::GetEminForTotalCS
void GetEminForTotalCS(G4ParticleDefinition *aPartDef, const G4MaterialCutsCouple *aCouple, G4double &emin_adj, G4double &emin_fwd)
Definition: G4AdjointCSManager.cc:427
G4VEmAdjointModel::ComputeAdjointCrossSectionVectorPerAtomForScatProj
std::vector< std::vector< double > * > ComputeAdjointCrossSectionVectorPerAtomForScatProj(G4double kinEnergyProd, G4double Z, G4double A=0., G4int nbin_pro_decade=10)
Definition: G4VEmAdjointModel.cc:291
G4AdjointGamma.hh
G4AdjointCSManager::GetForwardParticleEquivalent
G4ParticleDefinition * GetForwardParticleEquivalent(G4ParticleDefinition *theAdjPartDef)
Definition: G4AdjointCSManager.cc:954
G4VEmModel::GetChargeSquareRatio
virtual G4double GetChargeSquareRatio(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
Definition: G4VEmModel.cc:322
G4VEmAdjointModel::AdjointCrossSection
virtual G4double AdjointCrossSection(const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool IsScatProjToProjCase)
Definition: G4VEmAdjointModel.cc:58
G4ProductionCutsTable::GetProductionCutsTable
static G4ProductionCutsTable * GetProductionCutsTable()
Definition: G4ProductionCutsTable.cc:63
G4AdjointCSManager::GetContinuousWeightCorrection
G4double GetContinuousWeightCorrection(G4ParticleDefinition *aPartDef, G4double PreStepEkin, G4double AfterStepEkin, const G4MaterialCutsCouple *aCouple, G4double step_length)
Definition: G4AdjointCSManager.cc:502
G4Electron.hh
G4AdjointCSManager::RegisterEmProcess
void RegisterEmProcess(G4VEmProcess *aProcess, G4ParticleDefinition *aPartDef)
Definition: G4AdjointCSManager.cc:130
G4AdjointCSManager.hh
G4ProductionCutsTable.hh
G4ParticleDefinition::GetPDGMass
G4double GetPDGMass() const
Definition: G4ParticleDefinition.hh:161
G4ProductionCutsTable::GetMaterialCutsCouple
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
Definition: G4ProductionCutsTable.hh:262
G4INCL::Math::max
T max(const T t1, const T t2)
brief Return the largest of the two arguments
Definition: G4INCLGlobals.hh:116
G4AdjointProton::AdjointProton
static G4AdjointProton * AdjointProton()
Definition: G4AdjointProton.cc:84
G4VEmAdjointModel::ComputeAdjointCrossSectionVectorPerAtomForSecond
std::vector< std::vector< double > * > ComputeAdjointCrossSectionVectorPerAtomForSecond(G4double kinEnergyProd, G4double Z, G4double A=0., G4int nbin_pro_decade=10)
Definition: G4VEmAdjointModel.cc:238
G4VEmAdjointModel::GetSecondAdjEnergyMinForScatProjToProjCase
virtual G4double GetSecondAdjEnergyMinForScatProjToProjCase(G4double PrimAdjEnergy, G4double Tcut=0)
Definition: G4VEmAdjointModel.cc:664
G4VEmAdjointModel::GetHighEnergyLimit
G4double GetHighEnergyLimit()
Definition: G4VEmAdjointModel.hh:185
G4AdjointInterpolator::LinearInterpolation
G4double LinearInterpolation(G4double &x, G4double &x1, G4double &x2, G4double &y1, G4double &y2)
Definition: G4AdjointInterpolator.cc:64
G4AdjointCSManager::GetMaxAdjTotalCS
void GetMaxAdjTotalCS(G4ParticleDefinition *aPartDef, const G4MaterialCutsCouple *aCouple, G4double &e_sigma_max, G4double &sigma_max)
Definition: G4AdjointCSManager.cc:452
G4INCL::Math::min
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
Definition: G4INCLGlobals.hh:121
G4VEmAdjointModel::ComputeAdjointCrossSectionVectorPerVolumeForSecond
std::vector< std::vector< double > * > ComputeAdjointCrossSectionVectorPerVolumeForSecond(G4Material *aMaterial, G4double kinEnergyProd, G4int nbin_pro_decade=10)
Definition: G4VEmAdjointModel.cc:348
G4AdjointElectron.hh
G4AdjointInterpolator::GetInstance
static G4AdjointInterpolator * GetInstance()
Definition: G4AdjointInterpolator.cc:41
G4VEmAdjointModel::GetSecondPartOfSameType
G4bool GetSecondPartOfSameType()
Definition: G4VEmAdjointModel.hh:203
G4Electron::Electron
static G4Electron * Electron()
Definition: G4Electron.cc:94
G4endl
#define G4endl
Definition: G4ios.hh:61
G4Material::GetNumberOfElements
size_t GetNumberOfElements() const
Definition: G4Material.hh:184
G4VEmAdjointModel::ComputeAdjointCrossSectionVectorPerVolumeForScatProj
std::vector< std::vector< double > * > ComputeAdjointCrossSectionVectorPerVolumeForScatProj(G4Material *aMaterial, G4double kinEnergyProd, G4int nbin_pro_decade=10)
Definition: G4VEmAdjointModel.cc:399
G4ElementTable
std::vector< G4Element * > G4ElementTable
Definition: G4ElementTable.hh:41
G4double
double G4double
Definition: G4Types.hh:76
G4SystemOfUnits.hh
G4AdjointCSManager::SampleElementFromCSMatrices
G4Element * SampleElementFromCSMatrices(G4Material *aMaterial, G4VEmAdjointModel *aModel, G4double PrimEnergy, G4double Tcut, G4bool IsScatProjToProjCase)
Definition: G4AdjointCSManager.cc:662
G4Element::GetElementTable
static G4ElementTable * GetElementTable()
Definition: G4Element.cc:395
G4PhysicsTableHelper.hh
G4AdjointCSMatrix.hh
G4AdjointCSManager::GetAdjointCSManager
static G4AdjointCSManager * GetAdjointCSManager()
Definition: G4AdjointCSManager.cc:61
G4VEmAdjointModel.hh
G4Gamma.hh
G4MaterialCutsCouple::GetMaterial
const G4Material * GetMaterial() const
Definition: G4MaterialCutsCouple.hh:150
G4AdjointCSManager::ComputeAdjointCS
G4double ComputeAdjointCS(G4Material *aMaterial, G4VEmAdjointModel *aModel, G4double PrimEnergy, G4double Tcut, G4bool IsScatProjToProjCase, std::vector< G4double > &AdjointCS_for_each_element)
Definition: G4AdjointCSManager.cc:532
Generated on Wed Apr 30 2014 15:54:36 for Geant4.10 by   doxygen 1.8.7