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

Last change on this file since 1168 was 966, checked in by garnier, 15 years ago

fichier ajoutes
File size: 33.5 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer                                           *
4// *                                                                  *
5// * The  Geant4 software  is  copyright of the Copyright Holders  of *
6// * the Geant4 Collaboration.  It is provided  under  the terms  and *
7// * conditions of the Geant4 Software License,  included in the file *
8// * LICENSE and available at  http://cern.ch/geant4/license .  These *
9// * include a list of copyright holders.                             *
10// *                                                                  *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work  make  any representation or  warranty, express or implied, *
14// * regarding  this  software system or assume any liability for its *
15// * use.  Please see the license in the file  LICENSE  and URL above *
16// * for the full disclaimer and the limitation of liability.         *
17// *                                                                  *
18// * This  code  implementation is the result of  the  scientific and *
19// * technical work of the GEANT4 collaboration.                      *
20// * By using,  copying,  modifying or  distributing the software (or *
21// * any work based  on the software)  you  agree  to acknowledge its *
22// * use  in  resulting  scientific  publications,  and indicate your *
23// * acceptance of all terms of the Geant4 Software license.          *
24// ********************************************************************
25//
26#include "G4VEmAdjointModel.hh"
27#include "G4AdjointCSManager.hh"
28
29
30#include "G4Integrator.hh"
31#include "G4TrackStatus.hh"
32#include "G4ParticleChange.hh"
33#include "G4AdjointElectron.hh"
34#include "G4AdjointInterpolator.hh"
35
36////////////////////////////////////////////////////////////////////////////////
37//
38G4VEmAdjointModel::G4VEmAdjointModel(const G4String& nam):
39name(nam)
40// lowLimit(0.1*keV), highLimit(100.0*TeV), fluc(0), name(nam), pParticleChange(0)
41{ G4AdjointCSManager::GetAdjointCSManager()->RegisterEmAdjointModel(this);
42  CorrectWeightMode =true;
43  UseMatrix =true;
44  UseMatrixPerElement = true;
45  ApplyCutInRange = true;
46  ApplyBiasing = true;
47  UseOnlyOneMatrixForAllElements = true;
48  IsIonisation =true; 
49  CS_biasing_factor =1.;
50  //ApplyBiasing = false;
51}
52////////////////////////////////////////////////////////////////////////////////
53//
54G4VEmAdjointModel::~G4VEmAdjointModel()
55{;}
56////////////////////////////////////////////////////////////////////////////////
57//
58void G4VEmAdjointModel::SampleSecondaries(const G4Track& aTrack,
59                       G4bool IsScatProjToProjCase,
60                       G4ParticleChange* fParticleChange)
61{ 
62
63  const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
64  //DefineCurrentMaterial(aTrack->GetMaterialCutsCouple());
65  size_t ind=0;
66  if (!UseMatrixPerElement) ind = currentMaterialIndex;
67  //G4cout<<theAdjointPrimary<<std::endl;
68  else if (!UseOnlyOneMatrixForAllElements) { //Select Material
69        std::vector<double>* CS_Vs_Element = &CS_Vs_ElementForScatProjToProjCase;
70        if ( !IsScatProjToProjCase) CS_Vs_Element = &CS_Vs_ElementForProdToProjCase;
71        G4double rand_var= G4UniformRand();
72        G4double SumCS=0.;
73        for (size_t i=0;i<CS_Vs_Element->size();i++){
74                SumCS+=(*CS_Vs_Element)[i];
75                if (rand_var<=SumCS/lastCS){
76                        ind=i;
77                        break;
78                }
79        }
80        ind = currentMaterial->GetElement(ind)->GetIndex();
81  }
82 
83 
84 
85 //Elastic inverse scattering //not correct in all the cases
86 //---------------------------------------------------------
87 G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
88 G4double adjointPrimTotalEnergy = theAdjointPrimary->GetTotalEnergy();
89 G4double adjointPrimP =theAdjointPrimary->GetTotalMomentum();
90 //G4cout<<adjointPrimKinEnergy<<std::endl;
91 if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
92        return;
93 }
94 
95 //Sample secondary energy
96 //-----------------------
97 G4double projectileKinEnergy;
98// if (!IsIonisation  ) {
99        projectileKinEnergy = SampleAdjSecEnergyFromCSMatrix(ind,
100                                                  adjointPrimKinEnergy,
101                                                  IsScatProjToProjCase);
102 //}
103 /*else {
104        projectileKinEnergy =   SampleAdjSecEnergyFromDiffCrossSectionPerAtom(adjointPrimKinEnergy,IsScatProjToProjCase);
105        //G4cout<<projectileKinEnergy<<std::endl;
106 }*/   
107 //Weight correction
108 //-----------------------
109 
110 CorrectPostStepWeight(fParticleChange, aTrack.GetWeight(), adjointPrimKinEnergy,projectileKinEnergy);                                     
111 
112                 
113 
114 
115 
116 //Kinematic
117 //---------
118 
119 G4double projectileM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
120 G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
121 G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0;       
122 
123 
124 
125 //Companion
126 //-----------
127 G4double companionM0;
128 companionM0=(adjointPrimTotalEnergy-adjointPrimKinEnergy);
129 if (IsScatProjToProjCase) {
130        companionM0=theAdjEquivOfDirectSecondPartDef->GetPDGMass();
131 } 
132 G4double companionTotalEnergy =companionM0+ projectileKinEnergy-adjointPrimKinEnergy;
133 G4double companionP2 = companionTotalEnergy*companionTotalEnergy - companionM0*companionM0;   
134 
135 
136 //Projectile momentum
137 //--------------------
138 G4double  P_parallel = (adjointPrimP*adjointPrimP +  projectileP2 - companionP2)/(2.*adjointPrimP);
139 G4double  P_perp = std::sqrt( projectileP2 -  P_parallel*P_parallel);
140 G4ThreeVector dir_parallel=theAdjointPrimary->GetMomentumDirection();
141 G4double phi =G4UniformRand()*2.*3.1415926;
142 G4ThreeVector projectileMomentum = G4ThreeVector(P_perp*std::cos(phi),P_perp*std::sin(phi),P_parallel);
143 projectileMomentum.rotateUz(dir_parallel);
144 
145 
146 
147 if (!IsScatProjToProjCase && CorrectWeightMode){ //kill the primary and add a secondary
148        fParticleChange->ProposeTrackStatus(fStopAndKill);
149        fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum));
150        //G4cout<<"projectileMomentum "<<projectileMomentum<<std::endl;
151 }
152 else {
153        fParticleChange->ProposeEnergy(projectileKinEnergy);
154        fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
155 }     
156}
157////////////////////////////////////////////////////////////////////////////////
158//
159void G4VEmAdjointModel::CorrectPostStepWeight(G4ParticleChange* fParticleChange, G4double old_weight,  G4double , G4double ) 
160{
161 G4double new_weight=old_weight;
162 if (CorrectWeightMode) {
163        G4double w_corr =1./CS_biasing_factor;
164        //G4cout<<w_corr<<std::endl;
165       
166        /*G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection(theAdjEquivOfDirectPrimPartDef,
167                                                                               theAdjEquivOfDirectSecondPartDef,
168                                                                               adjointPrimKinEnergy,projectileKinEnergy,
169                                                                               aTrack.GetMaterialCutsCouple());
170        w_corr = projectileKinEnergy;
171        G4double Emin,Emax;
172        if (IsScatProjToProjCase) {
173                Emax = GetSecondAdjEnergyMaxForScatProjToProjCase(adjointPrimKinEnergy);
174                Emin = GetSecondAdjEnergyMinForScatProjToProjCase(adjointPrimKinEnergy, currentTcutForDirectSecond);
175               
176        }
177        else {
178                Emax = GetSecondAdjEnergyMaxForProdToProjCase(adjointPrimKinEnergy);
179                Emin = GetSecondAdjEnergyMinForProdToProjCase(adjointPrimKinEnergy);   
180        }
181        w_corr *=std::log(Emax/Emin)/(Emax-Emin); */   
182               
183        new_weight*=w_corr;
184 }
185 G4cout<< "new weight"<<new_weight<<std::endl;
186 fParticleChange->SetParentWeightByProcess(false);
187 fParticleChange->SetSecondaryWeightByProcess(false);
188 fParticleChange->ProposeParentWeight(new_weight);     
189}
190////////////////////////////////////////////////////////////////////////////////
191//                             
192G4double G4VEmAdjointModel::AdjointCrossSection(const G4MaterialCutsCouple* aCouple,
193                                G4double primEnergy,
194                                G4bool IsScatProjToProjCase)
195{ 
196  DefineCurrentMaterial(aCouple);
197  //G4double fwdCS = G4AdjointCSManager::GetAdjointCSManager()->GetTotalForwardCS(G4AdjointElectron::AdjointElectron(),primEnergy,aCouple);
198  //G4double adjCS = G4AdjointCSManager::GetAdjointCSManager()->GetTotalAdjointCS(G4AdjointElectron::AdjointElectron(), primEnergy,aCouple);
199  if (IsScatProjToProjCase){
200        lastCS = G4AdjointCSManager::GetAdjointCSManager()->ComputeAdjointCS(currentMaterial,
201                                                                        this, 
202                                                                        primEnergy,
203                                                                        currentTcutForDirectSecond,
204                                                                        true,
205                                                                        CS_Vs_ElementForScatProjToProjCase);
206        /*G4double fwdCS = G4AdjointCSManager::GetAdjointCSManager()->GetTotalForwardCS(theAdjEquivOfDirectPrimPartDef,primEnergy,aCouple);
207        G4double adjCS = G4AdjointCSManager::GetAdjointCSManager()->GetTotalAdjointCS(theAdjEquivOfDirectPrimPartDef, primEnergy,aCouple);
208        */
209        //if (adjCS >0 )lastCS *=fwdCS/adjCS;
210 
211  }
212  else {
213        lastCS = G4AdjointCSManager::GetAdjointCSManager()->ComputeAdjointCS(currentMaterial,
214                                                                        this, 
215                                                                        primEnergy,
216                                                                        currentTcutForDirectSecond,
217                                                                        false,
218                                                                        CS_Vs_ElementForProdToProjCase);
219        /*G4double fwdCS = G4AdjointCSManager::GetAdjointCSManager()->GetTotalForwardCS(theAdjEquivOfDirectSecondPartDef,primEnergy,aCouple);
220        G4double adjCS = G4AdjointCSManager::GetAdjointCSManager()->GetTotalAdjointCS(theAdjEquivOfDirectSecondPartDef, primEnergy,aCouple);
221        */
222        //if (adjCS >0 )lastCS *=fwdCS/adjCS;                                                           
223        //lastCS=0.;                                                           
224  }
225 
226 
227  return lastCS;
228                                                                       
229}                               
230////////////////////////////////////////////////////////////////////////////////
231//
232//The implementation here is correct for energy loss process, for the photoelectric and compton scattering the method should be redefine 
233G4double G4VEmAdjointModel::DiffCrossSectionPerAtomPrimToSecond(
234                                      G4double kinEnergyProj, 
235                                      G4double kinEnergyProd,
236                                      G4double Z, 
237                                      G4double A)
238{
239 G4double dSigmadEprod=0;
240 G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
241 G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
242 
243 
244 if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){ //the produced particle should have a kinetic energy smaller than the projectile
245        G4double Tmax=kinEnergyProj;
246        if (second_part_of_same_type) Tmax = kinEnergyProj/2.;
247        return Z*DiffCrossSectionMoller(kinEnergyProj,kinEnergyProd);
248        //it could be thta Tmax here should be DBLMAX
249        //Tmax=DBLMAX;
250       
251        G4double E1=kinEnergyProd;
252        G4double E2=kinEnergyProd*1.000001;
253        G4double dE=(E2-E1);
254        G4double sigma1=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E1,1.e20);
255        G4double sigma2=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E2,1.e20);
256       
257        dSigmadEprod=(sigma1-sigma2)/dE;
258        if (dSigmadEprod>1.) {
259                G4cout<<"sigma1 "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<sigma1<<std::endl;
260                G4cout<<"sigma2 "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<sigma2<<std::endl;
261                G4cout<<"dsigma "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<dSigmadEprod<<std::endl;
262               
263        }
264       
265       
266       
267 }
268 
269 
270
271       
272 return dSigmadEprod;   
273 
274 
275 
276}
277//The implementation here is correct for energy loss process, for the photoelectric and compton scattering the method should be redefine
278////////////////////////////////////////////////////////////////////////////////
279//
280G4double G4VEmAdjointModel::DiffCrossSectionPerAtomPrimToScatPrim(
281                                      G4double kinEnergyProj, 
282                                      G4double kinEnergyScatProj,
283                                      G4double Z, 
284                                      G4double A)
285{ G4double kinEnergyProd = kinEnergyProj - kinEnergyScatProj;
286  G4double dSigmadEprod;
287  if (kinEnergyProd <=0) dSigmadEprod=0;
288  else dSigmadEprod=DiffCrossSectionPerAtomPrimToSecond(kinEnergyProj,kinEnergyProd,Z,A);
289  return dSigmadEprod; 
290
291}
292
293////////////////////////////////////////////////////////////////////////////////
294//
295//The implementation here is correct for energy loss process, for the photoelectric and compton scattering the method should be redefine 
296G4double G4VEmAdjointModel::DiffCrossSectionPerVolumePrimToSecond(
297                                      const G4Material* aMaterial,
298                                      G4double kinEnergyProj, 
299                                      G4double kinEnergyProd)
300{
301 G4double dSigmadEprod=0;
302 G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
303 G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
304 
305 
306 if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){ 
307        G4double Tmax=kinEnergyProj;
308        if (second_part_of_same_type) Tmax = kinEnergyProj/2.;
309        //it could be thta Tmax here should be DBLMAX
310        //Tmax=DBLMAX;
311       
312        G4double E1=kinEnergyProd;
313       
314        G4double E2=kinEnergyProd*1.0001;
315        G4double dE=(E2-E1);
316        G4double sigma1=theDirectEMModel->CrossSectionPerVolume(aMaterial,theDirectPrimaryPartDef,kinEnergyProj,E1,E2);
317       
318        //G4double sigma2=theDirectEMModel->CrossSectionPerVolume(aMaterial,theDirectPrimaryPartDef,kinEnergyProj,E2,1.e50);
319        dSigmadEprod=sigma1/dE;
320        if (dSigmadEprod <0) { //could happen with bremstrahlung dur to suppression effect
321                G4cout<<"Halllllllllllllllllllllllllllllllllllllllllllllllo "<<kinEnergyProj<<'\t'<<E1<<'\t'<<dSigmadEprod<<std::endl;
322                E1=kinEnergyProd;
323                E2=E1*1.1;
324                dE=E2-E1;
325                sigma1=theDirectEMModel->CrossSectionPerVolume(aMaterial,theDirectPrimaryPartDef,kinEnergyProj,E1,1.e50);
326                G4double sigma2=theDirectEMModel->CrossSectionPerVolume(aMaterial,theDirectPrimaryPartDef,kinEnergyProj,E2,1.e50);
327                dSigmadEprod=(sigma1-sigma2)/dE;
328                G4cout<<dSigmadEprod<<std::endl;
329        }       
330       
331       
332 }
333
334
335       
336 return dSigmadEprod;   
337 
338 
339 
340}
341//The implementation here is correct for energy loss process, for the photoelectric and compton scattering the method should be redefine
342////////////////////////////////////////////////////////////////////////////////
343//
344G4double G4VEmAdjointModel::DiffCrossSectionPerVolumePrimToScatPrim(
345                                      const G4Material* aMaterial,
346                                      G4double kinEnergyProj, 
347                                      G4double kinEnergyScatProj)
348{ G4double kinEnergyProd = kinEnergyProj - kinEnergyScatProj;
349  G4double dSigmadEprod;
350  if (kinEnergyProd <=0) dSigmadEprod=0;
351  else dSigmadEprod=DiffCrossSectionPerVolumePrimToSecond(aMaterial,kinEnergyProj,kinEnergyProd);
352  return dSigmadEprod; 
353
354}
355///////////////////////////////////////////////////////////////////////////////////////////////////////////
356//
357G4double G4VEmAdjointModel::DiffCrossSectionFunction1(G4double kinEnergyProj){
358  //return kinEnergyProj*kinEnergyProj;
359  //ApplyBiasing=false;
360  G4double bias_factor = CS_biasing_factor*kinEnergyProdForIntegration/kinEnergyProj;   
361  if (!ApplyBiasing) bias_factor =CS_biasing_factor; 
362  //G4cout<<bias_factor<<std::endl;
363  if (UseMatrixPerElement ) {
364        return DiffCrossSectionPerAtomPrimToSecond(kinEnergyProj,kinEnergyProdForIntegration,ZSelectedNucleus,ASelectedNucleus)*bias_factor;
365  }
366  else {
367        return DiffCrossSectionPerVolumePrimToSecond(SelectedMaterial,kinEnergyProj,kinEnergyProdForIntegration)*bias_factor;
368  }     
369}
370//////////////////////////////////////////////////////////////////////////////
371//
372G4double G4VEmAdjointModel::DiffCrossSectionMoller(G4double kinEnergyProj,G4double kinEnergyProd){
373  G4double electron_mass_c2=0.51099906*MeV;
374  G4double energy = kinEnergyProj + electron_mass_c2;
375  G4double x   = kinEnergyProd/kinEnergyProj;
376  G4double gam    = energy/electron_mass_c2;
377  G4double gamma2 = gam*gam;
378  G4double beta2  = 1.0 - 1.0/gamma2;
379 
380  G4double g = (2.0*gam - 1.0)/gamma2;
381  G4double y = 1.0 - x;
382  G4double fac=twopi_mc2_rcl2/electron_mass_c2;
383  G4double dCS = fac*( 1.-g + ((1.0 - g*x)/(x*x)) + ((1.0 - g*y)/(y*y)))/(beta2*(gam-1));
384  return dCS/kinEnergyProj;
385 
386 
387
388} 
389////////////////////////////////////////////////////////////////////////////////
390//
391G4double G4VEmAdjointModel::DiffCrossSectionFunction2(G4double kinEnergyProj){
392  //return kinEnergyProj*kinEnergyProj;
393  G4double bias_factor =  CS_biasing_factor*kinEnergyScatProjForIntegration/kinEnergyProj;
394  //ApplyBiasing=false;
395 if (!ApplyBiasing) bias_factor = CS_biasing_factor;
396 //G4cout<<bias_factor<<std::endl;
397 if (UseMatrixPerElement ) {
398        return DiffCrossSectionPerAtomPrimToScatPrim(kinEnergyProj,kinEnergyScatProjForIntegration,ZSelectedNucleus,ASelectedNucleus)*bias_factor;
399 }     
400 else {
401        return DiffCrossSectionPerVolumePrimToScatPrim(SelectedMaterial,kinEnergyProj,kinEnergyScatProjForIntegration)*bias_factor;
402 
403 }     
404               
405}
406////////////////////////////////////////////////////////////////////////////////
407//
408std::vector< std::vector<G4double>* > G4VEmAdjointModel::ComputeAdjointCrossSectionVectorPerAtomForSecond(     
409                                G4double kinEnergyProd,
410                                G4double Z, 
411                                G4double A ,
412                                G4int nbin_pro_decade) //nb bins pro order of magnitude of energy
413{ G4Integrator<G4VEmAdjointModel, G4double(G4VEmAdjointModel::*)(G4double)> integral;
414  ASelectedNucleus= G4int(A);
415  ZSelectedNucleus=G4int(Z);
416  kinEnergyProdForIntegration = kinEnergyProd;
417 
418  //compute the vector of integrated cross sections
419  //-------------------
420 
421  G4double minEProj= GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
422  G4double maxEProj= GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
423  G4double E1=minEProj;
424  std::vector< G4double >*  log_ESec_vector = new  std::vector< G4double >();
425  std::vector< G4double >*  log_Prob_vector = new  std::vector< G4double >();
426  log_ESec_vector->clear();
427  log_Prob_vector->clear();
428  log_ESec_vector->push_back(std::log(E1));
429  log_Prob_vector->push_back(-50.);
430 
431  G4double E2=std::pow(10.,G4double( G4int(std::log10(minEProj)*nbin_pro_decade)+1)/nbin_pro_decade);
432  G4double fE=std::pow(10.,1./nbin_pro_decade);
433  G4double int_cross_section=0.;
434 
435  if (std::pow(fE,5.)>(maxEProj/minEProj)) fE = std::pow(maxEProj/minEProj,0.2);
436 
437  while (E1 <maxEProj*0.9999999){
438        //G4cout<<E1<<'\t'<<E2<<std::endl;
439       
440        int_cross_section +=integral.Simpson(this, &G4VEmAdjointModel::DiffCrossSectionFunction1,E1,std::min(E2,maxEProj*0.99999999), 10);
441        //G4cout<<"int_cross_section 1 "<<'\t'<<int_cross_section<<std::endl;
442        log_ESec_vector->push_back(std::log(std::min(E2,maxEProj)));
443        log_Prob_vector->push_back(std::log(int_cross_section));       
444        E1=E2;
445        E2*=fE;
446 
447  }
448  std::vector< std::vector<G4double>* > res_mat;
449  res_mat.clear();
450  if (int_cross_section >0.) {
451        res_mat.push_back(log_ESec_vector);
452        res_mat.push_back(log_Prob_vector);
453  }     
454 
455  return res_mat;
456}
457
458/////////////////////////////////////////////////////////////////////////////////////
459//
460std::vector< std::vector<G4double>* > G4VEmAdjointModel::ComputeAdjointCrossSectionVectorPerAtomForScatProj(
461                                G4double kinEnergyScatProj,
462                                G4double Z, 
463                                G4double A ,
464                                G4int nbin_pro_decade) //nb bins pro order of magnitude of energy
465{ G4Integrator<G4VEmAdjointModel, G4double(G4VEmAdjointModel::*)(G4double)> integral;
466  ASelectedNucleus=G4int(A);
467  ZSelectedNucleus=G4int(Z);
468  kinEnergyScatProjForIntegration = kinEnergyScatProj;
469 
470  //compute the vector of integrated cross sections
471  //-------------------
472 
473  G4double minEProj= GetSecondAdjEnergyMinForScatProjToProjCase(kinEnergyScatProj); 
474  G4double maxEProj= GetSecondAdjEnergyMaxForScatProjToProjCase(kinEnergyScatProj);
475  G4double dEmax=maxEProj-kinEnergyScatProj;
476  G4double dEmin=GetLowEnergyLimit();
477  G4double dE1=dEmin;
478  G4double dE2=dEmin;
479 
480 
481  std::vector< G4double >*  log_ESec_vector = new std::vector< G4double >();
482  std::vector< G4double >*  log_Prob_vector = new std::vector< G4double >();
483  log_ESec_vector->push_back(std::log(dEmin));
484  log_Prob_vector->push_back(-50.);
485  G4int nbins=std::max( G4int(std::log10(dEmax/dEmin))*nbin_pro_decade,5);
486  G4double fE=std::pow(dEmax/dEmin,1./nbins);
487 
488  G4double int_cross_section=0.;
489 
490  while (dE1 <dEmax*0.9999999999999){
491        dE2=dE1*fE;
492        int_cross_section +=integral.Simpson(this,
493        &G4VEmAdjointModel::DiffCrossSectionFunction2,minEProj+dE1,std::min(minEProj+dE2,maxEProj), 20);
494        //G4cout<<"int_cross_section "<<minEProj+dE1<<'\t'<<int_cross_section<<std::endl;
495        log_ESec_vector->push_back(std::log(std::min(dE2,maxEProj)));
496        log_Prob_vector->push_back(std::log(int_cross_section));       
497        dE1=dE2;
498
499  }
500  /*G4cout<<"total int_cross_section"<<'\t'<<int_cross_section<<std::endl;
501  G4cout<<"energy "<<kinEnergyScatProj<<std::endl;*/
502 
503 
504 
505 
506  std::vector< std::vector<G4double> *> res_mat; 
507  res_mat.clear();
508  if (int_cross_section >0.) {
509        res_mat.push_back(log_ESec_vector);
510        res_mat.push_back(log_Prob_vector);
511  }     
512 
513  return res_mat;
514}
515////////////////////////////////////////////////////////////////////////////////
516//
517std::vector< std::vector<G4double>* > G4VEmAdjointModel::ComputeAdjointCrossSectionVectorPerVolumeForSecond(     
518                                G4Material* aMaterial,
519                                G4double kinEnergyProd,
520                                G4int nbin_pro_decade) //nb bins pro order of magnitude of energy
521{ G4Integrator<G4VEmAdjointModel, G4double(G4VEmAdjointModel::*)(G4double)> integral;
522  SelectedMaterial= aMaterial;
523  kinEnergyProdForIntegration = kinEnergyProd;
524  //G4cout<<aMaterial->GetName()<<std::endl;
525  //G4cout<<kinEnergyProd/MeV<<std::endl;
526  //compute the vector of integrated cross sections
527  //-------------------
528 
529  G4double minEProj= GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
530  G4double maxEProj= GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
531  G4double E1=minEProj;
532  std::vector< G4double >*  log_ESec_vector = new  std::vector< G4double >();
533  std::vector< G4double >*  log_Prob_vector = new  std::vector< G4double >();
534  log_ESec_vector->clear();
535  log_Prob_vector->clear();
536  log_ESec_vector->push_back(std::log(E1));
537  log_Prob_vector->push_back(-50.);
538 
539  G4double E2=std::pow(10.,G4double( G4int(std::log10(minEProj)*nbin_pro_decade)+1)/nbin_pro_decade);
540  G4double fE=std::pow(10.,1./nbin_pro_decade);
541  G4double int_cross_section=0.;
542 
543  if (std::pow(fE,5.)>(maxEProj/minEProj)) fE = std::pow(maxEProj/minEProj,0.2);
544 
545  while (E1 <maxEProj*0.9999999){
546        //G4cout<<E1<<'\t'<<E2<<std::endl;
547       
548        int_cross_section +=integral.Simpson(this, &G4VEmAdjointModel::DiffCrossSectionFunction1,E1,std::min(E2,maxEProj*0.99999999), 10);
549        //G4cout<<"int_cross_section 1 "<<E1<<'\t'<<int_cross_section<<std::endl;
550        log_ESec_vector->push_back(std::log(std::min(E2,maxEProj)));
551        log_Prob_vector->push_back(std::log(int_cross_section));
552        E1=E2;
553        E2*=fE;
554 
555  }
556  std::vector< std::vector<G4double>* > res_mat;
557  res_mat.clear();
558 
559  //if (int_cross_section >0.) {
560        res_mat.push_back(log_ESec_vector);
561        res_mat.push_back(log_Prob_vector);
562  //}   
563 
564  return res_mat;
565}
566
567/////////////////////////////////////////////////////////////////////////////////////
568//
569std::vector< std::vector<G4double>* > G4VEmAdjointModel::ComputeAdjointCrossSectionVectorPerVolumeForScatProj(
570                                G4Material* aMaterial,
571                                G4double kinEnergyScatProj,
572                                G4int nbin_pro_decade) //nb bins pro order of magnitude of energy
573{ G4Integrator<G4VEmAdjointModel, G4double(G4VEmAdjointModel::*)(G4double)> integral;
574  SelectedMaterial= aMaterial;
575  kinEnergyScatProjForIntegration = kinEnergyScatProj;
576  /*G4cout<<name<<std::endl;
577  G4cout<<aMaterial->GetName()<<std::endl;
578  G4cout<<kinEnergyScatProj/MeV<<std::endl;*/
579  //compute the vector of integrated cross sections
580  //-------------------
581 
582  G4double minEProj= GetSecondAdjEnergyMinForScatProjToProjCase(kinEnergyScatProj);
583  G4double maxEProj= GetSecondAdjEnergyMaxForScatProjToProjCase(kinEnergyScatProj);
584 
585 
586  G4double dEmax=maxEProj-kinEnergyScatProj;
587  G4double dEmin=GetLowEnergyLimit();
588  G4double dE1=dEmin;
589  G4double dE2=dEmin;
590 
591 
592  std::vector< G4double >*  log_ESec_vector = new std::vector< G4double >();
593  std::vector< G4double >*  log_Prob_vector = new std::vector< G4double >();
594  log_ESec_vector->push_back(std::log(dEmin));
595  log_Prob_vector->push_back(-50.);
596  G4int nbins=std::max( G4int(std::log10(dEmax/dEmin))*nbin_pro_decade,5);
597  G4double fE=std::pow(dEmax/dEmin,1./nbins);
598 
599  G4double int_cross_section=0.;
600 
601  while (dE1 <dEmax*0.9999999999999){
602        dE2=dE1*fE;
603        int_cross_section +=integral.Simpson(this,
604        &G4VEmAdjointModel::DiffCrossSectionFunction2,minEProj+dE1,std::min(minEProj+dE2,maxEProj), 20);
605        //G4cout<<"int_cross_section "<<minEProj+dE1<<'\t'<<int_cross_section<<std::endl;
606        log_ESec_vector->push_back(std::log(std::min(dE2,maxEProj)));
607        log_Prob_vector->push_back(std::log(int_cross_section));
608        dE1=dE2;
609
610  }
611 
612 
613 
614 
615 
616  std::vector< std::vector<G4double> *> res_mat;
617  res_mat.clear();
618  if (int_cross_section >0.) {
619        res_mat.push_back(log_ESec_vector);
620        res_mat.push_back(log_Prob_vector);
621  }     
622 
623  return res_mat;
624}
625//////////////////////////////////////////////////////////////////////////////
626//                             
627G4double G4VEmAdjointModel::SampleAdjSecEnergyFromCSMatrix(size_t MatrixIndex,G4double aPrimEnergy,G4bool IsScatProjToProjCase)
628{ 
629 
630 
631  G4AdjointCSMatrix* theMatrix= (*pOnCSMatrixForProdToProjBackwardScattering)[MatrixIndex];
632  if (IsScatProjToProjCase) theMatrix= (*pOnCSMatrixForScatProjToProjBackwardScattering)[MatrixIndex];
633  std::vector< G4double >* theLogPrimEnergyVector = theMatrix->GetLogPrimEnergyVector();
634  //G4double dLog = theMatrix->GetDlog();
635 
636 
637 
638  if (theLogPrimEnergyVector->size() ==0){
639        G4cout<<"No data are contained in the given AdjointCSMatrix!"<<std::endl;
640        G4cout<<"The sampling procedure will be stopped."<<std::endl;
641        return 0.;
642       
643  }
644 
645  G4AdjointInterpolator* theInterpolator=G4AdjointInterpolator::GetInstance();
646  G4double aLogPrimEnergy = std::log(aPrimEnergy);
647  size_t ind =theInterpolator->FindPositionForLogVector(aLogPrimEnergy,*theLogPrimEnergyVector);
648 
649 
650  G4double aLogPrimEnergy1,aLogPrimEnergy2;
651  G4double aLogCS1,aLogCS2;
652   G4double log01,log02;
653  std::vector< G4double>* aLogSecondEnergyVector1 =0;
654  std::vector< G4double>* aLogSecondEnergyVector2  =0;
655  std::vector< G4double>* aLogProbVector1=0;
656  std::vector< G4double>* aLogProbVector2=0; 
657  std::vector< size_t>* aLogProbVectorIndex1=0;
658  std::vector< size_t>* aLogProbVectorIndex2=0;
659                                                                     
660  theMatrix->GetData(ind, aLogPrimEnergy1,aLogCS1,log01, aLogSecondEnergyVector1,aLogProbVector1,aLogProbVectorIndex1);
661  theMatrix->GetData(ind+1, aLogPrimEnergy2,aLogCS2,log02, aLogSecondEnergyVector2,aLogProbVector2,aLogProbVectorIndex2);
662 
663  G4double rand_var = G4UniformRand();
664  G4double log_rand_var= std::log(rand_var);
665  G4double log_Tcut =std::log(currentTcutForDirectSecond);
666  G4double Esec=0;
667  G4double log_dE1,log_dE2;
668  G4double log_rand_var1,log_rand_var2;
669  G4double log_E1,log_E2;
670  log_rand_var1=log_rand_var;
671  log_rand_var2=log_rand_var;
672 
673  G4double Emin=0.;
674  G4double Emax=0.;
675  if (theMatrix->IsScatProjToProjCase()){ //case where Tcut plays a role
676        //G4cout<<"Here "<<std::endl;
677        if (ApplyCutInRange) {
678                if (second_part_of_same_type && currentTcutForDirectSecond>aPrimEnergy) return aPrimEnergy;
679                /*if (IsIonisation){
680                        G4double inv_Tcut= 1./currentTcutForDirectSecond;
681                        G4double inv_dE=inv_Tcut-rand_var*(inv_Tcut-1./aPrimEnergy);
682                        Esec= aPrimEnergy+1./inv_dE;
683                        //return Esec;
684                        G4double dE1=currentTcutForDirectSecond;
685                        G4double dE2=currentTcutForDirectSecond*1.00001;
686                        G4double dCS1=DiffCrossSectionMoller(aPrimEnergy+dE1,dE1);
687                        G4double dCS2=DiffCrossSectionMoller(aPrimEnergy+dE2,dE2);
688                        G4double alpha1=std::log(dCS1/dCS2)/std::log(dE1/dE2);
689                        G4double a1=dCS1/std::pow(dE1,alpha1);
690                        dCS1=DiffCrossSectionMoller(aPrimEnergy+dE1,dE1);
691                        dCS2=DiffCrossSectionMoller(aPrimEnergy+dE2,dE2);
692                       
693                        return Esec;
694                       
695                       
696                       
697                        dE1=aPrimEnergy/1.00001;
698                        dE2=aPrimEnergy;
699                        dCS1=DiffCrossSectionMoller(aPrimEnergy+dE1,dE1);
700                        dCS2=DiffCrossSectionMoller(aPrimEnergy+dE2,dE2);
701                        G4double alpha2=std::log(dCS1/dCS2)/std::log(dE1/dE2);
702                        G4double a2=dCS1/std::pow(dE1,alpha1);
703                        return Esec;
704                       
705                       
706                       
707                       
708                }*/     
709                log_rand_var1=log_rand_var+theInterpolator->InterpolateForLogVector(log_Tcut,*aLogSecondEnergyVector1,*aLogProbVector1);
710                log_rand_var2=log_rand_var+theInterpolator->InterpolateForLogVector(log_Tcut,*aLogSecondEnergyVector2,*aLogProbVector2);
711               
712        }       
713        log_dE1 = theInterpolator->Interpolate(log_rand_var1,*aLogProbVector1,*aLogSecondEnergyVector1,"Lin");
714        log_dE2 = theInterpolator->Interpolate(log_rand_var2,*aLogProbVector2,*aLogSecondEnergyVector2,"Lin");
715       
716        /*log_dE1 = theInterpolator->InterpolateWithIndexVector(log_rand_var1,*aLogProbVector1,*aLogSecondEnergyVector1,*aLogProbVectorIndex1,log01,dLog);
717        log_dE2 = theInterpolator->InterpolateWithIndexVector(log_rand_var1,*aLogProbVector1,*aLogSecondEnergyVector1,*aLogProbVectorIndex1,log02,dLog);
718        */                                 
719                                           
720       
721       
722       
723       
724        Esec = aPrimEnergy +
725        std::exp(theInterpolator->LinearInterpolation(aLogPrimEnergy,aLogPrimEnergy1,aLogPrimEnergy2,log_dE1,log_dE2));
726       
727        Emin=GetSecondAdjEnergyMinForScatProjToProjCase(aPrimEnergy);
728        Emax=GetSecondAdjEnergyMaxForScatProjToProjCase(aPrimEnergy);
729        Esec=std::max(Esec,Emin);
730        Esec=std::min(Esec,Emax);
731       
732       
733        //G4cout<<"Esec "<<Esec<<std::endl;
734        //if (Esec > 2.*aPrimEnergy && second_part_of_same_type) Esec = 2.*aPrimEnergy;
735  }
736  else { //Tcut condition is already full-filled
737        /*G4cout<<"Start "<<std::endl;
738        G4cout<<std::exp((*aLogProbVector1)[0])<<std::endl;
739        G4cout<<std::exp((*aLogProbVector2)[0])<<std::endl;*/
740        /*G4double inv_E1= .5/aPrimEnergy;
741               
742        G4double inv_E=inv_E1-rand_var*(inv_E1-0.00001);
743        Esec= 1./inv_E;
744        return Esec;*/
745        log_E1 = theInterpolator->Interpolate(log_rand_var,*aLogProbVector1,*aLogSecondEnergyVector1,"Lin");
746        log_E2 = theInterpolator->Interpolate(log_rand_var,*aLogProbVector2,*aLogSecondEnergyVector2,"Lin");
747        /*log_E1 = theInterpolator->InterpolateWithIndexVector(log_rand_var1,*aLogProbVector1,*aLogSecondEnergyVector1,*aLogProbVectorIndex1,log01,dLog);
748        log_E2 = theInterpolator->InterpolateWithIndexVector(log_rand_var1,*aLogProbVector1,*aLogSecondEnergyVector1,*aLogProbVectorIndex1,log02,dLog);
749        */
750       
751       
752        /*G4cout<<std::exp(log_E1)<<std::endl;
753        G4cout<<std::exp(log_E2)<<std::endl;*/
754       
755        Esec = std::exp(theInterpolator->LinearInterpolation(aLogPrimEnergy,aLogPrimEnergy1,aLogPrimEnergy2,log_E1,log_E2));
756        Emin=GetSecondAdjEnergyMinForProdToProjCase(aPrimEnergy);
757        Emax=GetSecondAdjEnergyMaxForProdToProjCase(aPrimEnergy);
758        Esec=std::max(Esec,Emin);
759        Esec=std::min(Esec,Emax);
760       
761  }
762       
763  return Esec;
764 
765 
766 
767 
768                                                                                   
769}
770//////////////////////////////////////////////////////////////////////////////
771//                             
772G4double G4VEmAdjointModel::SampleAdjSecEnergyFromDiffCrossSectionPerAtom(G4double prim_energy,G4bool IsScatProjToProjCase)
773{ 
774  // here we try to use the rejection method
775  //-----------------------------------------
776 
777  G4double E=0;
778  G4double x,xmin,greject,q;
779  if ( IsScatProjToProjCase){
780        G4double Emax = GetSecondAdjEnergyMaxForScatProjToProjCase(prim_energy);
781        G4double Emin= prim_energy+currentTcutForDirectSecond;
782        xmin=Emin/Emax;
783        G4double grejmax = DiffCrossSectionPerAtomPrimToScatPrim(Emin,prim_energy,1)*prim_energy; 
784
785        do {
786                q = G4UniformRand();
787                x = 1./(q*(1./xmin -1.) +1.);
788                E=x*Emax;
789                greject = DiffCrossSectionPerAtomPrimToScatPrim( E,prim_energy ,1)*prim_energy; 
790               
791        }       
792       
793        while( greject < G4UniformRand()*grejmax );     
794       
795  }
796  else {
797        G4double Emax = GetSecondAdjEnergyMaxForProdToProjCase(prim_energy);
798        G4double Emin=  GetSecondAdjEnergyMinForProdToProjCase(prim_energy);;
799        xmin=Emin/Emax;
800        G4double grejmax = DiffCrossSectionPerAtomPrimToSecond(Emin,prim_energy,1); 
801        do {
802                q = G4UniformRand();
803                x = std::pow(xmin, q);
804                E=x*Emax;
805                greject = DiffCrossSectionPerAtomPrimToSecond( E,prim_energy ,1); 
806               
807        }       
808       
809        while( greject < G4UniformRand()*grejmax );
810 
811 
812 
813  }
814 
815  return E;
816 
817 
818 
819 
820                                                                                   
821}
822//////////////////////////////////////////////////////////////////////////////
823//                             
824G4double G4VEmAdjointModel::GetSecondAdjEnergyMaxForScatProjToProjCase(G4double kinEnergyScatProj)
825{ G4double maxEProj= HighEnergyLimit;
826  if (second_part_of_same_type)  maxEProj=std::min(kinEnergyScatProj*2.,HighEnergyLimit);
827  return maxEProj;
828}
829//////////////////////////////////////////////////////////////////////////////
830//
831G4double G4VEmAdjointModel::GetSecondAdjEnergyMinForScatProjToProjCase(G4double PrimAdjEnergy,G4double Tcut)
832{ return PrimAdjEnergy+Tcut;
833}
834//////////////////////////////////////////////////////////////////////////////
835//                             
836G4double G4VEmAdjointModel::GetSecondAdjEnergyMaxForProdToProjCase(G4double )
837{ return HighEnergyLimit;
838}
839//////////////////////////////////////////////////////////////////////////////
840//
841G4double G4VEmAdjointModel::GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy)
842{ G4double minEProj=PrimAdjEnergy;
843  if (second_part_of_same_type)  minEProj=PrimAdjEnergy*2.;
844  return minEProj;
845}
846////////////////////////////////////////////////////////////////////////////////////////////
847//
848void  G4VEmAdjointModel::DefineCurrentMaterial(const G4MaterialCutsCouple* couple)
849{ if(couple != currentCouple) {
850        currentCouple   = const_cast<G4MaterialCutsCouple*> (couple);
851        currentMaterial = const_cast<G4Material*> (couple->GetMaterial());
852        currentCoupleIndex = couple->GetIndex();
853        currentMaterialIndex = currentMaterial->GetIndex();
854        size_t idx=56;
855   
856        if (theAdjEquivOfDirectPrimPartDef) {
857                if (theAdjEquivOfDirectPrimPartDef->GetParticleName() == "adj_gamma") idx = 0;
858                else if (theAdjEquivOfDirectPrimPartDef->GetParticleName() == "adj_e-") idx = 1;
859                else if (theAdjEquivOfDirectPrimPartDef->GetParticleName() == "adj_e+") idx = 2;
860                const std::vector<G4double>* aVec = G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(idx);
861                currentTcutForDirectPrim=(*aVec)[currentCoupleIndex];
862        }       
863   
864   
865        if (theAdjEquivOfDirectPrimPartDef == theAdjEquivOfDirectSecondPartDef) {
866                currentTcutForDirectSecond = currentTcutForDirectPrim;
867        }
868        else { 
869                if (theAdjEquivOfDirectSecondPartDef){
870                        if (theAdjEquivOfDirectSecondPartDef->GetParticleName() == "adj_gamma") idx = 0;
871                        else if (theAdjEquivOfDirectSecondPartDef->GetParticleName() == "adj_e-") idx = 1;
872                        else if (theAdjEquivOfDirectSecondPartDef->GetParticleName() == "adj_e+") idx = 2;
873                        const std::vector<G4double>* aVec = G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(idx);
874                        currentTcutForDirectSecond=(*aVec)[currentCoupleIndex];
875                }
876        }
877  }
878}
Note: See TracBrowser for help on using the repository browser.