source: trunk/source/processes/electromagnetic/adjoint/src/G4VEmAdjointModel.cc @ 1317

Last change on this file since 1317 was 1228, checked in by garnier, 15 years ago

update geant4.9.3 tag
File size: 27.2 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer                                           *
4// *                                                                  *
5// * The  Geant4 software  is  copyright of the Copyright Holders  of *
6// * the Geant4 Collaboration.  It is provided  under  the terms  and *
7// * conditions of the Geant4 Software License,  included in the file *
8// * LICENSE and available at  http://cern.ch/geant4/license .  These *
9// * include a list of copyright holders.                             *
10// *                                                                  *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work  make  any representation or  warranty, express or implied, *
14// * regarding  this  software system or assume any liability for its *
15// * use.  Please see the license in the file  LICENSE  and URL above *
16// * for the full disclaimer and the limitation of liability.         *
17// *                                                                  *
18// * This  code  implementation is the result of  the  scientific and *
19// * technical work of the GEANT4 collaboration.                      *
20// * By using,  copying,  modifying or  distributing the software (or *
21// * any work based  on the software)  you  agree  to acknowledge its *
22// * use  in  resulting  scientific  publications,  and indicate your *
23// * acceptance of all terms of the Geant4 Software license.          *
24// ********************************************************************
25//
26// $Id: G4VEmAdjointModel.cc,v 1.5 2009/12/16 17:50:09 gunter Exp $
27// GEANT4 tag $Name: geant4-09-03 $
28//
29#include "G4VEmAdjointModel.hh"
30#include "G4AdjointCSManager.hh"
31#include "G4Integrator.hh"
32#include "G4TrackStatus.hh"
33#include "G4ParticleChange.hh"
34#include "G4AdjointElectron.hh"
35#include "G4AdjointInterpolator.hh"
36#include "G4PhysicsTable.hh"
37
38////////////////////////////////////////////////////////////////////////////////
39//
40G4VEmAdjointModel::G4VEmAdjointModel(const G4String& nam):
41name(nam)
42// lowLimit(0.1*keV), highLimit(100.0*TeV), fluc(0), name(nam), pParticleChange(0)
43{ 
44  G4AdjointCSManager::GetAdjointCSManager()->RegisterEmAdjointModel(this);
45  second_part_of_same_type =false;
46  theDirectEMModel=0;
47}
48////////////////////////////////////////////////////////////////////////////////
49//
50G4VEmAdjointModel::~G4VEmAdjointModel()
51{;}
52////////////////////////////////////////////////////////////////////////////////
53//                             
54G4double G4VEmAdjointModel::AdjointCrossSection(const G4MaterialCutsCouple* aCouple,
55                                G4double primEnergy,
56                                G4bool IsScatProjToProjCase)
57{ 
58  DefineCurrentMaterial(aCouple);
59  preStepEnergy=primEnergy;
60 
61  std::vector<G4double>* CS_Vs_Element = &CS_Vs_ElementForProdToProjCase;
62  if (IsScatProjToProjCase)   CS_Vs_Element = &CS_Vs_ElementForScatProjToProjCase;
63  lastCS = G4AdjointCSManager::GetAdjointCSManager()->ComputeAdjointCS(currentMaterial,
64                                                                        this, 
65                                                                        primEnergy,
66                                                                        currentTcutForDirectSecond,
67                                                                        IsScatProjToProjCase,
68                                                                        *CS_Vs_Element);
69  if (IsScatProjToProjCase) lastAdjointCSForScatProjToProjCase = lastCS;
70  else lastAdjointCSForProdToProjCase =lastCS;                                                         
71       
72 
73 
74  return lastCS;
75                                                                       
76}                               
77////////////////////////////////////////////////////////////////////////////////
78//
79//General implementation correct for energy loss process, for the photoelectric and compton scattering the method should be redefine 
80G4double G4VEmAdjointModel::DiffCrossSectionPerAtomPrimToSecond(
81                                      G4double kinEnergyProj, 
82                                      G4double kinEnergyProd,
83                                      G4double Z, 
84                                      G4double A)
85{
86 G4double dSigmadEprod=0;
87 G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
88 G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
89 
90 
91 if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){ //the produced particle should have a kinetic energy smaller than the projectile
92        G4double Tmax=kinEnergyProj;
93        if (second_part_of_same_type) Tmax = kinEnergyProj/2.;
94       
95        G4double E1=kinEnergyProd;
96        G4double E2=kinEnergyProd*1.000001;
97        G4double dE=(E2-E1);
98        G4double sigma1=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E1,1.e20);
99        G4double sigma2=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E2,1.e20);
100       
101        dSigmadEprod=(sigma1-sigma2)/dE;
102 }
103 return dSigmadEprod;   
104 
105 
106 
107}
108//The implementation here is correct for energy loss process, for the photoelectric and compton scattering the method should be redefine
109////////////////////////////////////////////////////////////////////////////////
110//
111G4double G4VEmAdjointModel::DiffCrossSectionPerAtomPrimToScatPrim(
112                                      G4double kinEnergyProj, 
113                                      G4double kinEnergyScatProj,
114                                      G4double Z, 
115                                      G4double A)
116{ G4double kinEnergyProd = kinEnergyProj - kinEnergyScatProj;
117  G4double dSigmadEprod;
118  if (kinEnergyProd <=0) dSigmadEprod=0;
119  else dSigmadEprod=DiffCrossSectionPerAtomPrimToSecond(kinEnergyProj,kinEnergyProd,Z,A);
120  return dSigmadEprod; 
121
122}
123
124////////////////////////////////////////////////////////////////////////////////
125//
126//The implementation here is correct for energy loss process, for the photoelectric and compton scattering the method should be redefine 
127G4double G4VEmAdjointModel::DiffCrossSectionPerVolumePrimToSecond(
128                                      const G4Material* aMaterial,
129                                      G4double kinEnergyProj, 
130                                      G4double kinEnergyProd)
131{
132 G4double dSigmadEprod=0;
133 G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
134 G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
135 
136 
137 if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){ 
138        G4double Tmax=kinEnergyProj;
139        if (second_part_of_same_type) Tmax = kinEnergyProj/2.;
140        G4double E1=kinEnergyProd;
141        G4double E2=kinEnergyProd*1.0001;
142        G4double dE=(E2-E1);
143        G4double sigma1=theDirectEMModel->CrossSectionPerVolume(aMaterial,theDirectPrimaryPartDef,kinEnergyProj,E1,E2);
144        G4double sigma2=theDirectEMModel->CrossSectionPerVolume(aMaterial,theDirectPrimaryPartDef,kinEnergyProj,E2,1.e50);
145        dSigmadEprod=(sigma1-sigma2)/dE;
146 }
147 return dSigmadEprod;   
148 
149 
150 
151}
152//The implementation here is correct for energy loss process, for the photoelectric and compton scattering the method should be redefine
153////////////////////////////////////////////////////////////////////////////////
154//
155G4double G4VEmAdjointModel::DiffCrossSectionPerVolumePrimToScatPrim(
156                                      const G4Material* aMaterial,
157                                      G4double kinEnergyProj, 
158                                      G4double kinEnergyScatProj)
159{ G4double kinEnergyProd = kinEnergyProj - kinEnergyScatProj;
160  G4double dSigmadEprod;
161  if (kinEnergyProd <=0) dSigmadEprod=0;
162  else dSigmadEprod=DiffCrossSectionPerVolumePrimToSecond(aMaterial,kinEnergyProj,kinEnergyProd);
163  return dSigmadEprod; 
164
165}
166///////////////////////////////////////////////////////////////////////////////////////////////////////////
167//
168G4double G4VEmAdjointModel::DiffCrossSectionFunction1(G4double kinEnergyProj){
169 
170 
171  G4double bias_factor = CS_biasing_factor*kinEnergyProdForIntegration/kinEnergyProj;   
172
173
174  if (UseMatrixPerElement ) {
175        return DiffCrossSectionPerAtomPrimToSecond(kinEnergyProj,kinEnergyProdForIntegration,ZSelectedNucleus,ASelectedNucleus)*bias_factor;
176  }
177  else  {
178        return DiffCrossSectionPerVolumePrimToSecond(SelectedMaterial,kinEnergyProj,kinEnergyProdForIntegration)*bias_factor;
179  }     
180}
181
182////////////////////////////////////////////////////////////////////////////////
183//
184G4double G4VEmAdjointModel::DiffCrossSectionFunction2(G4double kinEnergyProj){
185 
186 G4double bias_factor =  CS_biasing_factor*kinEnergyScatProjForIntegration/kinEnergyProj;
187 if (UseMatrixPerElement ) {
188        return DiffCrossSectionPerAtomPrimToScatPrim(kinEnergyProj,kinEnergyScatProjForIntegration,ZSelectedNucleus,ASelectedNucleus)*bias_factor;
189 }     
190 else { 
191        return DiffCrossSectionPerVolumePrimToScatPrim(SelectedMaterial,kinEnergyProj,kinEnergyScatProjForIntegration)*bias_factor;
192 
193 }     
194               
195}
196////////////////////////////////////////////////////////////////////////////////
197//
198
199G4double G4VEmAdjointModel::DiffCrossSectionPerVolumeFunctionForIntegrationOverEkinProj(G4double kinEnergyProd)
200{
201  return DiffCrossSectionPerVolumePrimToSecond(SelectedMaterial,kinEnergyProjForIntegration,kinEnergyProd);
202}
203////////////////////////////////////////////////////////////////////////////////
204//
205std::vector< std::vector<G4double>* > G4VEmAdjointModel::ComputeAdjointCrossSectionVectorPerAtomForSecond(     
206                                G4double kinEnergyProd,
207                                G4double Z, 
208                                G4double A ,
209                                G4int nbin_pro_decade) //nb bins pro order of magnitude of energy
210{ 
211  G4Integrator<G4VEmAdjointModel, double(G4VEmAdjointModel::*)(double)> integral;
212  ASelectedNucleus= int(A);
213  ZSelectedNucleus=int(Z);
214  kinEnergyProdForIntegration = kinEnergyProd;
215 
216  //compute the vector of integrated cross sections
217  //-------------------
218 
219  G4double minEProj= GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
220  G4double maxEProj= GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
221  G4double E1=minEProj;
222  std::vector< double>*  log_ESec_vector = new  std::vector< double>();
223  std::vector< double>*  log_Prob_vector = new  std::vector< double>();
224  log_ESec_vector->clear();
225  log_Prob_vector->clear();
226  log_ESec_vector->push_back(std::log(E1));
227  log_Prob_vector->push_back(-50.);
228 
229  G4double E2=std::pow(10.,double( int(std::log10(minEProj)*nbin_pro_decade)+1)/nbin_pro_decade);
230  G4double fE=std::pow(10.,1./nbin_pro_decade);
231  G4double int_cross_section=0.;
232 
233  if (std::pow(fE,5.)>(maxEProj/minEProj)) fE = std::pow(maxEProj/minEProj,0.2);
234 
235  while (E1 <maxEProj*0.9999999){
236        //G4cout<<E1<<'\t'<<E2<<G4endl;
237       
238        int_cross_section +=integral.Simpson(this,
239        &G4VEmAdjointModel::DiffCrossSectionFunction1,E1,std::min(E2,maxEProj*0.99999999), 5);
240        log_ESec_vector->push_back(std::log(std::min(E2,maxEProj)));
241        log_Prob_vector->push_back(std::log(int_cross_section));       
242        E1=E2;
243        E2*=fE;
244 
245  }
246  std::vector< std::vector<G4double>* > res_mat;
247  res_mat.clear();
248  if (int_cross_section >0.) {
249        res_mat.push_back(log_ESec_vector);
250        res_mat.push_back(log_Prob_vector);
251  }     
252 
253  return res_mat;
254}
255
256/////////////////////////////////////////////////////////////////////////////////////
257//
258std::vector< std::vector<G4double>* > G4VEmAdjointModel::ComputeAdjointCrossSectionVectorPerAtomForScatProj(
259                                G4double kinEnergyScatProj,
260                                G4double Z, 
261                                G4double A ,
262                                G4int nbin_pro_decade) //nb bins pro order of magnitude of energy
263{ G4Integrator<G4VEmAdjointModel, double(G4VEmAdjointModel::*)(double)> integral;
264  ASelectedNucleus=int(A);
265  ZSelectedNucleus=int(Z);
266  kinEnergyScatProjForIntegration = kinEnergyScatProj;
267 
268  //compute the vector of integrated cross sections
269  //-------------------
270 
271  G4double minEProj= GetSecondAdjEnergyMinForScatProjToProjCase(kinEnergyScatProj); 
272  G4double maxEProj= GetSecondAdjEnergyMaxForScatProjToProjCase(kinEnergyScatProj);
273  G4double dEmax=maxEProj-kinEnergyScatProj;
274  G4double dEmin=GetLowEnergyLimit();
275  G4double dE1=dEmin;
276  G4double dE2=dEmin;
277 
278 
279  std::vector< double>*  log_ESec_vector = new std::vector< double>();
280  std::vector< double>*  log_Prob_vector = new std::vector< double>();
281  log_ESec_vector->push_back(std::log(dEmin));
282  log_Prob_vector->push_back(-50.);
283  G4int nbins=std::max( int(std::log10(dEmax/dEmin))*nbin_pro_decade,5);
284  G4double fE=std::pow(dEmax/dEmin,1./nbins);
285 
286 
287 
288 
289 
290  G4double int_cross_section=0.;
291 
292  while (dE1 <dEmax*0.9999999999999){
293        dE2=dE1*fE;
294        int_cross_section +=integral.Simpson(this,
295        &G4VEmAdjointModel::DiffCrossSectionFunction2,minEProj+dE1,std::min(minEProj+dE2,maxEProj), 5);
296        //G4cout<<"int_cross_section "<<minEProj+dE1<<'\t'<<int_cross_section<<G4endl;
297        log_ESec_vector->push_back(std::log(std::min(dE2,maxEProj-minEProj)));
298        log_Prob_vector->push_back(std::log(int_cross_section));       
299        dE1=dE2;
300
301  }
302 
303 
304  std::vector< std::vector<G4double> *> res_mat; 
305  res_mat.clear();
306  if (int_cross_section >0.) {
307        res_mat.push_back(log_ESec_vector);
308        res_mat.push_back(log_Prob_vector);
309  }     
310 
311  return res_mat;
312}
313////////////////////////////////////////////////////////////////////////////////
314//
315std::vector< std::vector<G4double>* > G4VEmAdjointModel::ComputeAdjointCrossSectionVectorPerVolumeForSecond(     
316                                G4Material* aMaterial,
317                                G4double kinEnergyProd,
318                                G4int nbin_pro_decade) //nb bins pro order of magnitude of energy
319{ G4Integrator<G4VEmAdjointModel, double(G4VEmAdjointModel::*)(double)> integral;
320  SelectedMaterial= aMaterial;
321  kinEnergyProdForIntegration = kinEnergyProd;
322   //compute the vector of integrated cross sections
323  //-------------------
324 
325  G4double minEProj= GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
326  G4double maxEProj= GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
327  G4double E1=minEProj;
328  std::vector< double>*  log_ESec_vector = new  std::vector< double>();
329  std::vector< double>*  log_Prob_vector = new  std::vector< double>();
330  log_ESec_vector->clear();
331  log_Prob_vector->clear();
332  log_ESec_vector->push_back(std::log(E1));
333  log_Prob_vector->push_back(-50.);
334 
335  G4double E2=std::pow(10.,double( int(std::log10(minEProj)*nbin_pro_decade)+1)/nbin_pro_decade);
336  G4double fE=std::pow(10.,1./nbin_pro_decade);
337  G4double int_cross_section=0.;
338 
339  if (std::pow(fE,5.)>(maxEProj/minEProj)) fE = std::pow(maxEProj/minEProj,0.2);
340 
341  while (E1 <maxEProj*0.9999999){
342       
343        int_cross_section +=integral.Simpson(this,
344        &G4VEmAdjointModel::DiffCrossSectionFunction1,E1,std::min(E2,maxEProj*0.99999999), 5);
345        log_ESec_vector->push_back(std::log(std::min(E2,maxEProj)));
346        log_Prob_vector->push_back(std::log(int_cross_section));
347        E1=E2;
348        E2*=fE;
349 
350  }
351  std::vector< std::vector<G4double>* > res_mat;
352  res_mat.clear();
353 
354  if (int_cross_section >0.) {
355        res_mat.push_back(log_ESec_vector);
356        res_mat.push_back(log_Prob_vector);
357  }
358 
359 
360 
361  return res_mat;
362}
363
364/////////////////////////////////////////////////////////////////////////////////////
365//
366std::vector< std::vector<G4double>* > G4VEmAdjointModel::ComputeAdjointCrossSectionVectorPerVolumeForScatProj(
367                                G4Material* aMaterial,
368                                G4double kinEnergyScatProj,
369                                G4int nbin_pro_decade) //nb bins pro order of magnitude of energy
370{ G4Integrator<G4VEmAdjointModel, double(G4VEmAdjointModel::*)(double)> integral;
371  SelectedMaterial= aMaterial;
372  kinEnergyScatProjForIntegration = kinEnergyScatProj;
373 
374  //compute the vector of integrated cross sections
375  //-------------------
376 
377  G4double minEProj= GetSecondAdjEnergyMinForScatProjToProjCase(kinEnergyScatProj);
378  G4double maxEProj= GetSecondAdjEnergyMaxForScatProjToProjCase(kinEnergyScatProj);
379 
380 
381  G4double dEmax=maxEProj-kinEnergyScatProj;
382  G4double dEmin=GetLowEnergyLimit();
383  G4double dE1=dEmin;
384  G4double dE2=dEmin;
385 
386 
387  std::vector< double>*  log_ESec_vector = new std::vector< double>();
388  std::vector< double>*  log_Prob_vector = new std::vector< double>();
389  log_ESec_vector->push_back(std::log(dEmin));
390  log_Prob_vector->push_back(-50.);
391  G4int nbins=std::max( int(std::log10(dEmax/dEmin))*nbin_pro_decade,5);
392  G4double fE=std::pow(dEmax/dEmin,1./nbins);
393 
394  G4double int_cross_section=0.;
395 
396  while (dE1 <dEmax*0.9999999999999){
397        dE2=dE1*fE;
398        int_cross_section +=integral.Simpson(this,
399        &G4VEmAdjointModel::DiffCrossSectionFunction2,minEProj+dE1,std::min(minEProj+dE2,maxEProj), 5);
400        log_ESec_vector->push_back(std::log(std::min(dE2,maxEProj-minEProj)));
401        log_Prob_vector->push_back(std::log(int_cross_section));
402        dE1=dE2;
403
404  }
405 
406 
407 
408 
409 
410  std::vector< std::vector<G4double> *> res_mat;
411  res_mat.clear();
412  if (int_cross_section >0.) {
413        res_mat.push_back(log_ESec_vector);
414        res_mat.push_back(log_Prob_vector);
415  }     
416 
417  return res_mat;
418}
419//////////////////////////////////////////////////////////////////////////////
420//                             
421G4double G4VEmAdjointModel::SampleAdjSecEnergyFromCSMatrix(size_t MatrixIndex,G4double aPrimEnergy,G4bool IsScatProjToProjCase)
422{ 
423 
424 
425  G4AdjointCSMatrix* theMatrix= (*pOnCSMatrixForProdToProjBackwardScattering)[MatrixIndex];
426  if (IsScatProjToProjCase) theMatrix= (*pOnCSMatrixForScatProjToProjBackwardScattering)[MatrixIndex];
427  std::vector< double>* theLogPrimEnergyVector = theMatrix->GetLogPrimEnergyVector();
428 
429  if (theLogPrimEnergyVector->size() ==0){
430        G4cout<<"No data are contained in the given AdjointCSMatrix!"<<G4endl;
431        G4cout<<"The sampling procedure will be stopped."<<G4endl;
432        return 0.;
433       
434  }
435 
436  G4AdjointInterpolator* theInterpolator=G4AdjointInterpolator::GetInstance();
437  G4double aLogPrimEnergy = std::log(aPrimEnergy);
438  size_t ind =theInterpolator->FindPositionForLogVector(aLogPrimEnergy,*theLogPrimEnergyVector);
439 
440 
441  G4double aLogPrimEnergy1,aLogPrimEnergy2;
442  G4double aLogCS1,aLogCS2;
443  G4double log01,log02;
444  std::vector< double>* aLogSecondEnergyVector1 =0;
445  std::vector< double>* aLogSecondEnergyVector2  =0;
446  std::vector< double>* aLogProbVector1=0;
447  std::vector< double>* aLogProbVector2=0; 
448  std::vector< size_t>* aLogProbVectorIndex1=0;
449  std::vector< size_t>* aLogProbVectorIndex2=0;
450                                                                     
451  theMatrix->GetData(ind, aLogPrimEnergy1,aLogCS1,log01, aLogSecondEnergyVector1,aLogProbVector1,aLogProbVectorIndex1);
452  theMatrix->GetData(ind+1, aLogPrimEnergy2,aLogCS2,log02, aLogSecondEnergyVector2,aLogProbVector2,aLogProbVectorIndex2);
453 
454  G4double rand_var = G4UniformRand();
455  G4double log_rand_var= std::log(rand_var);
456  G4double log_Tcut =std::log(currentTcutForDirectSecond);
457  G4double Esec=0;
458  G4double log_dE1,log_dE2;
459  G4double log_rand_var1,log_rand_var2;
460  G4double log_E1,log_E2;
461  log_rand_var1=log_rand_var;
462  log_rand_var2=log_rand_var;
463 
464  G4double Emin=0.;
465  G4double Emax=0.;
466  if (theMatrix->IsScatProjToProjCase()){ //case where Tcut plays a role
467        Emin=GetSecondAdjEnergyMinForScatProjToProjCase(aPrimEnergy,currentTcutForDirectSecond);
468        Emax=GetSecondAdjEnergyMaxForScatProjToProjCase(aPrimEnergy);
469        G4double dE=0;
470        if (Emin < Emax ){
471            if (ApplyCutInRange) {
472                if (second_part_of_same_type && currentTcutForDirectSecond>aPrimEnergy) return aPrimEnergy;
473               
474                log_rand_var1=log_rand_var+theInterpolator->InterpolateForLogVector(log_Tcut,*aLogSecondEnergyVector1,*aLogProbVector1);
475                log_rand_var2=log_rand_var+theInterpolator->InterpolateForLogVector(log_Tcut,*aLogSecondEnergyVector2,*aLogProbVector2);
476               
477            }   
478            log_dE1 = theInterpolator->Interpolate(log_rand_var1,*aLogProbVector1,*aLogSecondEnergyVector1,"Lin");
479            log_dE2 = theInterpolator->Interpolate(log_rand_var2,*aLogProbVector2,*aLogSecondEnergyVector2,"Lin");
480             dE=std::exp(theInterpolator->LinearInterpolation(aLogPrimEnergy,aLogPrimEnergy1,aLogPrimEnergy2,log_dE1,log_dE2));
481        }
482       
483        Esec = aPrimEnergy +dE;
484        Esec=std::max(Esec,Emin);
485        Esec=std::min(Esec,Emax);
486       
487  }
488  else { //Tcut condition is already full-filled
489       
490        log_E1 = theInterpolator->Interpolate(log_rand_var,*aLogProbVector1,*aLogSecondEnergyVector1,"Lin");
491        log_E2 = theInterpolator->Interpolate(log_rand_var,*aLogProbVector2,*aLogSecondEnergyVector2,"Lin");
492       
493        Esec = std::exp(theInterpolator->LinearInterpolation(aLogPrimEnergy,aLogPrimEnergy1,aLogPrimEnergy2,log_E1,log_E2));
494        Emin=GetSecondAdjEnergyMinForProdToProjCase(aPrimEnergy);
495        Emax=GetSecondAdjEnergyMaxForProdToProjCase(aPrimEnergy);
496        Esec=std::max(Esec,Emin);
497        Esec=std::min(Esec,Emax);
498       
499  }
500       
501  return Esec;
502 
503 
504 
505 
506                                                                                   
507}
508
509//////////////////////////////////////////////////////////////////////////////
510//                             
511G4double G4VEmAdjointModel::SampleAdjSecEnergyFromCSMatrix(G4double aPrimEnergy,G4bool IsScatProjToProjCase)
512{ SelectCSMatrix(IsScatProjToProjCase);
513  return SampleAdjSecEnergyFromCSMatrix(indexOfUsedCrossSectionMatrix, aPrimEnergy, IsScatProjToProjCase);
514}
515//////////////////////////////////////////////////////////////////////////////
516//
517void G4VEmAdjointModel::SelectCSMatrix(G4bool IsScatProjToProjCase)
518{ 
519  indexOfUsedCrossSectionMatrix=0;
520  if (!UseMatrixPerElement) indexOfUsedCrossSectionMatrix = currentMaterialIndex;
521  else if (!UseOnlyOneMatrixForAllElements) { //Select Material
522        std::vector<G4double>* CS_Vs_Element = &CS_Vs_ElementForScatProjToProjCase;
523        lastCS=lastAdjointCSForScatProjToProjCase;
524        if ( !IsScatProjToProjCase) {
525                CS_Vs_Element = &CS_Vs_ElementForProdToProjCase;
526                lastCS=lastAdjointCSForProdToProjCase;
527        }       
528        G4double rand_var= G4UniformRand();
529        G4double SumCS=0.;
530        size_t ind=0;
531        for (size_t i=0;i<CS_Vs_Element->size();i++){
532                SumCS+=(*CS_Vs_Element)[i];
533                if (rand_var<=SumCS/lastCS){
534                        ind=i;
535                        break;
536                }
537        }
538        indexOfUsedCrossSectionMatrix = currentMaterial->GetElement(ind)->GetIndex();
539  }
540}
541//////////////////////////////////////////////////////////////////////////////
542//                             
543G4double G4VEmAdjointModel::SampleAdjSecEnergyFromDiffCrossSectionPerAtom(G4double prim_energy,G4bool IsScatProjToProjCase)
544{ 
545  // here we try to use the rejection method
546  //-----------------------------------------
547 
548  G4double E=0;
549  G4double x,xmin,greject,q;
550  if ( IsScatProjToProjCase){
551        G4double Emax = GetSecondAdjEnergyMaxForScatProjToProjCase(prim_energy);
552        G4double Emin= prim_energy+currentTcutForDirectSecond;
553        xmin=Emin/Emax;
554        G4double grejmax = DiffCrossSectionPerAtomPrimToScatPrim(Emin,prim_energy,1)*prim_energy; 
555
556        do {
557                q = G4UniformRand();
558                x = 1./(q*(1./xmin -1.) +1.);
559                E=x*Emax;
560                greject = DiffCrossSectionPerAtomPrimToScatPrim( E,prim_energy ,1)*prim_energy; 
561               
562        }       
563       
564        while( greject < G4UniformRand()*grejmax );     
565       
566  }
567  else {
568        G4double Emax = GetSecondAdjEnergyMaxForProdToProjCase(prim_energy);
569        G4double Emin=  GetSecondAdjEnergyMinForProdToProjCase(prim_energy);;
570        xmin=Emin/Emax;
571        G4double grejmax = DiffCrossSectionPerAtomPrimToSecond(Emin,prim_energy,1); 
572        do {
573                q = G4UniformRand();
574                x = std::pow(xmin, q);
575                E=x*Emax;
576                greject = DiffCrossSectionPerAtomPrimToSecond( E,prim_energy ,1); 
577               
578        }       
579       
580        while( greject < G4UniformRand()*grejmax );
581 
582 
583 
584  }
585 
586  return E;
587}
588
589////////////////////////////////////////////////////////////////////////////////
590//
591void G4VEmAdjointModel::CorrectPostStepWeight(G4ParticleChange* fParticleChange, 
592                                                  G4double old_weight, 
593                                                  G4double adjointPrimKinEnergy, 
594                                                  G4double projectileKinEnergy,
595                                                  G4bool IsScatProjToProjCase) 
596{
597 G4double new_weight=old_weight;
598 G4double w_corr =1./CS_biasing_factor;
599 w_corr*=G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection();
600 
601 
602 lastCS=lastAdjointCSForScatProjToProjCase;
603 if ( !IsScatProjToProjCase) lastCS=lastAdjointCSForProdToProjCase;
604 if (adjointPrimKinEnergy !=preStepEnergy){ //Is that in all cases needed???
605        G4double post_stepCS=AdjointCrossSection(currentCouple, adjointPrimKinEnergy
606                                                 ,IsScatProjToProjCase );
607        w_corr*=post_stepCS/lastCS;                     
608 }
609       
610 new_weight*=w_corr;
611
612 //G4cout<<"Post step "<<new_weight<<'\t'<<w_corr<<'\t'<<old_weight<<G4endl;
613 new_weight*=projectileKinEnergy/adjointPrimKinEnergy;//This is needed due to the biasing of diff CS
614                                                        //by the factor adjointPrimKinEnergy/projectileKinEnergy
615   
616
617
618 fParticleChange->SetParentWeightByProcess(false);
619 fParticleChange->SetSecondaryWeightByProcess(false);
620 fParticleChange->ProposeParentWeight(new_weight);     
621}
622//////////////////////////////////////////////////////////////////////////////
623//                             
624G4double G4VEmAdjointModel::GetSecondAdjEnergyMaxForScatProjToProjCase(G4double kinEnergyScatProj)
625{ G4double maxEProj= HighEnergyLimit;
626  if (second_part_of_same_type)  maxEProj=std::min(kinEnergyScatProj*2.,HighEnergyLimit);
627  return maxEProj;
628}
629//////////////////////////////////////////////////////////////////////////////
630//
631G4double G4VEmAdjointModel::GetSecondAdjEnergyMinForScatProjToProjCase(G4double PrimAdjEnergy,G4double Tcut)
632{ G4double Emin=PrimAdjEnergy;
633  if (ApplyCutInRange) Emin=PrimAdjEnergy+Tcut;
634  return Emin;
635}
636//////////////////////////////////////////////////////////////////////////////
637//                             
638G4double G4VEmAdjointModel::GetSecondAdjEnergyMaxForProdToProjCase(G4double )
639{ return HighEnergyLimit;
640}
641//////////////////////////////////////////////////////////////////////////////
642//
643G4double G4VEmAdjointModel::GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy)
644{ G4double minEProj=PrimAdjEnergy;
645  if (second_part_of_same_type)  minEProj=PrimAdjEnergy*2.;
646  return minEProj;
647}
648////////////////////////////////////////////////////////////////////////////////////////////
649//
650void  G4VEmAdjointModel::DefineCurrentMaterial(const G4MaterialCutsCouple* couple)
651{ if(couple != currentCouple) {
652        currentCouple   = const_cast<G4MaterialCutsCouple*> (couple);
653        currentMaterial = const_cast<G4Material*> (couple->GetMaterial());
654        currentCoupleIndex = couple->GetIndex();
655        currentMaterialIndex = currentMaterial->GetIndex();
656        size_t idx=56;
657        currentTcutForDirectPrim =0.00000000001;
658        if (theAdjEquivOfDirectPrimPartDef) {
659                if (theAdjEquivOfDirectPrimPartDef->GetParticleName() == "adj_gamma") idx = 0;
660                else if (theAdjEquivOfDirectPrimPartDef->GetParticleName() == "adj_e-") idx = 1;
661                else if (theAdjEquivOfDirectPrimPartDef->GetParticleName() == "adj_e+") idx = 2;
662                if (idx <56){
663                        const std::vector<G4double>* aVec = G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(idx);
664                        currentTcutForDirectPrim=(*aVec)[currentCoupleIndex];
665                }       
666        }       
667   
668        currentTcutForDirectSecond =0.00000000001;
669        if (theAdjEquivOfDirectPrimPartDef == theAdjEquivOfDirectSecondPartDef) {
670                currentTcutForDirectSecond = currentTcutForDirectPrim;
671        }
672        else { 
673                if (theAdjEquivOfDirectSecondPartDef){
674                        if (theAdjEquivOfDirectSecondPartDef->GetParticleName() == "adj_gamma") idx = 0;
675                        else if (theAdjEquivOfDirectSecondPartDef->GetParticleName() == "adj_e-") idx = 1;
676                        else if (theAdjEquivOfDirectSecondPartDef->GetParticleName() == "adj_e+") idx = 2;
677                        const std::vector<G4double>* aVec = G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(idx);
678                        currentTcutForDirectSecond=(*aVec)[currentCoupleIndex];
679                        if (idx <56){
680                                const std::vector<G4double>* aVec = G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(idx);
681                                currentTcutForDirectPrim=(*aVec)[currentCoupleIndex];
682                        }
683                       
684                       
685                }
686        }
687  }
688}
689////////////////////////////////////////////////////////////////////////////////////////////
690//
691void G4VEmAdjointModel::SetHighEnergyLimit(G4double aVal)
692{ HighEnergyLimit=aVal;
693  if (theDirectEMModel) theDirectEMModel->SetHighEnergyLimit( aVal);
694}
695////////////////////////////////////////////////////////////////////////////////////////////
696//
697void G4VEmAdjointModel::SetLowEnergyLimit(G4double aVal)
698{
699  LowEnergyLimit=aVal;
700  if (theDirectEMModel) theDirectEMModel->SetLowEnergyLimit( aVal);
701}
702////////////////////////////////////////////////////////////////////////////////////////////
703//
704void G4VEmAdjointModel::SetAdjointEquivalentOfDirectPrimaryParticleDefinition(G4ParticleDefinition* aPart)
705{
706  theAdjEquivOfDirectPrimPartDef=aPart;
707  if (theAdjEquivOfDirectPrimPartDef->GetParticleName() =="adj_e-")
708                                        theDirectPrimaryPartDef=G4Electron::Electron();
709  if (theAdjEquivOfDirectPrimPartDef->GetParticleName() =="adj_gamma")
710                                        theDirectPrimaryPartDef=G4Gamma::Gamma();
711}
Note: See TracBrowser for help on using the repository browser.