G4VEmAdjointModel.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // $Id: G4VEmAdjointModel.cc 69844 2013-05-16 09:19:33Z gcosmo $
00027 //
00028 #include "G4VEmAdjointModel.hh"
00029 #include "G4AdjointCSManager.hh"
00030 #include "G4Integrator.hh"
00031 #include "G4TrackStatus.hh"
00032 #include "G4ParticleChange.hh"
00033 #include "G4AdjointElectron.hh"
00034 #include "G4AdjointGamma.hh"
00035 #include "G4AdjointPositron.hh"
00036 #include "G4AdjointInterpolator.hh"
00037 #include "G4PhysicsTable.hh"
00038 
00040 //
00041 G4VEmAdjointModel::G4VEmAdjointModel(const G4String& nam):
00042 name(nam)
00043 // lowLimit(0.1*keV), highLimit(100.0*TeV), fluc(0), name(nam), pParticleChange(0)
00044 { 
00045   model_index = G4AdjointCSManager::GetAdjointCSManager()->RegisterEmAdjointModel(this);
00046   second_part_of_same_type =false;
00047   theDirectEMModel=0;
00048   mass_ratio_product=1.;
00049   mass_ratio_projectile=1.;
00050   currentCouple=0;
00051 }
00053 //
00054 G4VEmAdjointModel::~G4VEmAdjointModel()
00055 {;}
00057 //                              
00058 G4double G4VEmAdjointModel::AdjointCrossSection(const G4MaterialCutsCouple* aCouple,
00059                                 G4double primEnergy,
00060                                 G4bool IsScatProjToProjCase)
00061 { 
00062   DefineCurrentMaterial(aCouple);
00063   preStepEnergy=primEnergy;
00064   
00065   std::vector<G4double>* CS_Vs_Element = &CS_Vs_ElementForProdToProjCase;
00066   if (IsScatProjToProjCase)   CS_Vs_Element = &CS_Vs_ElementForScatProjToProjCase;
00067   lastCS = G4AdjointCSManager::GetAdjointCSManager()->ComputeAdjointCS(currentMaterial,
00068                                                                         this, 
00069                                                                         primEnergy,
00070                                                                         currentTcutForDirectSecond,
00071                                                                         IsScatProjToProjCase,
00072                                                                         *CS_Vs_Element);
00073   if (IsScatProjToProjCase) lastAdjointCSForScatProjToProjCase = lastCS;
00074   else lastAdjointCSForProdToProjCase =lastCS;                                                          
00075         
00076  
00077  
00078   return lastCS;
00079                                                                         
00080 }
00082 //
00083 G4double G4VEmAdjointModel::GetAdjointCrossSection(const G4MaterialCutsCouple* aCouple,
00084                                              G4double primEnergy,
00085                                              G4bool IsScatProjToProjCase)
00086 {
00087   return AdjointCrossSection(aCouple, primEnergy,
00088                                 IsScatProjToProjCase);
00089   
00090   /*
00091   //To continue
00092   DefineCurrentMaterial(aCouple);
00093   preStepEnergy=primEnergy;
00094   if (IsScatProjToProjCase){
00095         G4double ekin=primEnergy*mass_ratio_projectile;
00096         lastCS = G4AdjointCSManager::GetAdjointCSManager()->GetAdjointSigma(ekin, model_index,true, aCouple);
00097         lastAdjointCSForScatProjToProjCase = lastCS;
00098         //G4cout<<ekin<<std::endl;
00099   }
00100   else {
00101         G4double ekin=primEnergy*mass_ratio_product;
00102         lastCS = G4AdjointCSManager::GetAdjointCSManager()->GetAdjointSigma(ekin, model_index,false, aCouple);
00103         lastAdjointCSForProdToProjCase = lastCS;
00104         //G4cout<<ekin<<std::endl;
00105   }
00106   return lastCS;
00107   */
00108 }                                                                       
00110 //
00111 //General implementation correct for energy loss process, for the photoelectric and compton scattering the method should be redefine  
00112 G4double G4VEmAdjointModel::DiffCrossSectionPerAtomPrimToSecond(
00113                                       G4double kinEnergyProj, 
00114                                       G4double kinEnergyProd,
00115                                       G4double Z, 
00116                                       G4double A)
00117 {
00118  G4double dSigmadEprod=0;
00119  G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
00120  G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
00121  
00122  
00123  if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){ //the produced particle should have a kinetic energy smaller than the projectile 
00124         
00125         /*G4double Tmax=kinEnergyProj;
00126         if (second_part_of_same_type) Tmax = kinEnergyProj/2.;*/
00127 
00128         G4double E1=kinEnergyProd;
00129         G4double E2=kinEnergyProd*1.000001;
00130         G4double dE=(E2-E1);
00131         G4double sigma1=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E1,1.e20);
00132         G4double sigma2=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E2,1.e20);
00133         
00134         dSigmadEprod=(sigma1-sigma2)/dE;
00135  }
00136  return dSigmadEprod;   
00137  
00138  
00139  
00140 }
00141 //The implementation here is correct for energy loss process, for the photoelectric and compton scattering the method should be redefine 
00143 //
00144 G4double G4VEmAdjointModel::DiffCrossSectionPerAtomPrimToScatPrim(
00145                                       G4double kinEnergyProj, 
00146                                       G4double kinEnergyScatProj,
00147                                       G4double Z, 
00148                                       G4double A)
00149 { G4double kinEnergyProd = kinEnergyProj - kinEnergyScatProj;
00150   G4double dSigmadEprod;
00151   if (kinEnergyProd <=0) dSigmadEprod=0;
00152   else dSigmadEprod=DiffCrossSectionPerAtomPrimToSecond(kinEnergyProj,kinEnergyProd,Z,A);
00153   return dSigmadEprod;  
00154 
00155 }
00156 
00158 //
00159 //The implementation here is correct for energy loss process, for the photoelectric and compton scattering the method should be redefine  
00160 G4double G4VEmAdjointModel::DiffCrossSectionPerVolumePrimToSecond(
00161                                       const G4Material* aMaterial,
00162                                       G4double kinEnergyProj, 
00163                                       G4double kinEnergyProd)
00164 {
00165  G4double dSigmadEprod=0;
00166  G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
00167  G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
00168  
00169  
00170  if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){ 
00171         /*G4double Tmax=kinEnergyProj;
00172         if (second_part_of_same_type) Tmax = kinEnergyProj/2.;*/
00173         G4double E1=kinEnergyProd;
00174         G4double E2=kinEnergyProd*1.0001;
00175         G4double dE=(E2-E1);
00176         G4double sigma1=theDirectEMModel->CrossSectionPerVolume(aMaterial,theDirectPrimaryPartDef,kinEnergyProj,E1,1.e20);
00177     G4double sigma2=theDirectEMModel->CrossSectionPerVolume(aMaterial,theDirectPrimaryPartDef,kinEnergyProj,E2,1.e20);
00178         dSigmadEprod=(sigma1-sigma2)/dE;
00179  }
00180  return dSigmadEprod;   
00181  
00182  
00183  
00184 }
00185 //The implementation here is correct for energy loss process, for the photoelectric and compton scattering the method should be redefine 
00187 //
00188 G4double G4VEmAdjointModel::DiffCrossSectionPerVolumePrimToScatPrim(
00189                                       const G4Material* aMaterial,
00190                                       G4double kinEnergyProj, 
00191                                       G4double kinEnergyScatProj)
00192 { G4double kinEnergyProd = kinEnergyProj - kinEnergyScatProj;
00193   G4double dSigmadEprod;
00194   if (kinEnergyProd <=0) dSigmadEprod=0;
00195   else dSigmadEprod=DiffCrossSectionPerVolumePrimToSecond(aMaterial,kinEnergyProj,kinEnergyProd);
00196   return dSigmadEprod;  
00197 
00198 }
00200 //
00201 G4double G4VEmAdjointModel::DiffCrossSectionFunction1(G4double kinEnergyProj){
00202   
00203   
00204   G4double bias_factor = CS_biasing_factor*kinEnergyProdForIntegration/kinEnergyProj;   
00205 
00206 
00207   if (UseMatrixPerElement ) {
00208         return DiffCrossSectionPerAtomPrimToSecond(kinEnergyProj,kinEnergyProdForIntegration,ZSelectedNucleus,ASelectedNucleus)*bias_factor;
00209   }
00210   else  {
00211         return DiffCrossSectionPerVolumePrimToSecond(SelectedMaterial,kinEnergyProj,kinEnergyProdForIntegration)*bias_factor;
00212   }     
00213 }
00214 
00216 //
00217 G4double G4VEmAdjointModel::DiffCrossSectionFunction2(G4double kinEnergyProj){
00218   
00219  G4double bias_factor =  CS_biasing_factor*kinEnergyScatProjForIntegration/kinEnergyProj;
00220  if (UseMatrixPerElement ) {
00221         return DiffCrossSectionPerAtomPrimToScatPrim(kinEnergyProj,kinEnergyScatProjForIntegration,ZSelectedNucleus,ASelectedNucleus)*bias_factor;
00222  }      
00223  else {  
00224         return DiffCrossSectionPerVolumePrimToScatPrim(SelectedMaterial,kinEnergyProj,kinEnergyScatProjForIntegration)*bias_factor;
00225  
00226  }      
00227                 
00228 }
00230 //
00231 
00232 G4double G4VEmAdjointModel::DiffCrossSectionPerVolumeFunctionForIntegrationOverEkinProj(G4double kinEnergyProd)
00233 {
00234   return DiffCrossSectionPerVolumePrimToSecond(SelectedMaterial,kinEnergyProjForIntegration,kinEnergyProd);
00235 }
00237 //
00238 std::vector< std::vector<G4double>* > G4VEmAdjointModel::ComputeAdjointCrossSectionVectorPerAtomForSecond(      
00239                                 G4double kinEnergyProd,
00240                                 G4double Z, 
00241                                 G4double A ,
00242                                 G4int nbin_pro_decade) //nb bins pro order of magnitude of energy
00243 { 
00244   G4Integrator<G4VEmAdjointModel, double(G4VEmAdjointModel::*)(double)> integral;
00245   ASelectedNucleus= int(A);
00246   ZSelectedNucleus=int(Z);
00247   kinEnergyProdForIntegration = kinEnergyProd;
00248   
00249   //compute the vector of integrated cross sections
00250   //-------------------
00251   
00252   G4double minEProj= GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
00253   G4double maxEProj= GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
00254   G4double E1=minEProj;
00255   std::vector< double>*  log_ESec_vector = new  std::vector< double>();
00256   std::vector< double>*  log_Prob_vector = new  std::vector< double>();
00257   log_ESec_vector->clear();
00258   log_Prob_vector->clear();
00259   log_ESec_vector->push_back(std::log(E1));
00260   log_Prob_vector->push_back(-50.);
00261   
00262   G4double E2=std::pow(10.,double( int(std::log10(minEProj)*nbin_pro_decade)+1)/nbin_pro_decade);
00263   G4double fE=std::pow(10.,1./nbin_pro_decade);
00264   G4double int_cross_section=0.;
00265   
00266   if (std::pow(fE,5.)>(maxEProj/minEProj)) fE = std::pow(maxEProj/minEProj,0.2);
00267   
00268   while (E1 <maxEProj*0.9999999){
00269         //G4cout<<E1<<'\t'<<E2<<G4endl;
00270         
00271         int_cross_section +=integral.Simpson(this,
00272         &G4VEmAdjointModel::DiffCrossSectionFunction1,E1,std::min(E2,maxEProj*0.99999999), 5);
00273         log_ESec_vector->push_back(std::log(std::min(E2,maxEProj)));
00274         log_Prob_vector->push_back(std::log(int_cross_section));        
00275         E1=E2;
00276         E2*=fE;
00277   
00278   }
00279   std::vector< std::vector<G4double>* > res_mat;
00280   res_mat.clear();
00281   if (int_cross_section >0.) {
00282         res_mat.push_back(log_ESec_vector);
00283         res_mat.push_back(log_Prob_vector);
00284   }     
00285   
00286   return res_mat;
00287 }
00288 
00290 //
00291 std::vector< std::vector<G4double>* > G4VEmAdjointModel::ComputeAdjointCrossSectionVectorPerAtomForScatProj(
00292                                 G4double kinEnergyScatProj,
00293                                 G4double Z, 
00294                                 G4double A ,
00295                                 G4int nbin_pro_decade) //nb bins pro order of magnitude of energy
00296 { G4Integrator<G4VEmAdjointModel, double(G4VEmAdjointModel::*)(double)> integral;
00297   ASelectedNucleus=int(A);
00298   ZSelectedNucleus=int(Z);
00299   kinEnergyScatProjForIntegration = kinEnergyScatProj;
00300   
00301   //compute the vector of integrated cross sections
00302   //-------------------
00303   
00304   G4double minEProj= GetSecondAdjEnergyMinForScatProjToProjCase(kinEnergyScatProj); 
00305   G4double maxEProj= GetSecondAdjEnergyMaxForScatProjToProjCase(kinEnergyScatProj);
00306   G4double dEmax=maxEProj-kinEnergyScatProj;
00307   G4double dEmin=GetLowEnergyLimit();
00308   G4double dE1=dEmin;
00309   G4double dE2=dEmin;
00310   
00311   
00312   std::vector< double>*  log_ESec_vector = new std::vector< double>();
00313   std::vector< double>*  log_Prob_vector = new std::vector< double>();
00314   log_ESec_vector->push_back(std::log(dEmin));
00315   log_Prob_vector->push_back(-50.);
00316   G4int nbins=std::max( int(std::log10(dEmax/dEmin))*nbin_pro_decade,5);
00317   G4double fE=std::pow(dEmax/dEmin,1./nbins);
00318   
00319   
00320   
00321   
00322   
00323   G4double int_cross_section=0.;
00324   
00325   while (dE1 <dEmax*0.9999999999999){
00326         dE2=dE1*fE;
00327         int_cross_section +=integral.Simpson(this,
00328         &G4VEmAdjointModel::DiffCrossSectionFunction2,minEProj+dE1,std::min(minEProj+dE2,maxEProj), 5);
00329         //G4cout<<"int_cross_section "<<minEProj+dE1<<'\t'<<int_cross_section<<G4endl;
00330         log_ESec_vector->push_back(std::log(std::min(dE2,maxEProj-minEProj)));
00331         log_Prob_vector->push_back(std::log(int_cross_section));        
00332         dE1=dE2;
00333 
00334   }
00335   
00336   
00337   std::vector< std::vector<G4double> *> res_mat; 
00338   res_mat.clear();
00339   if (int_cross_section >0.) {
00340         res_mat.push_back(log_ESec_vector);
00341         res_mat.push_back(log_Prob_vector);
00342   }     
00343   
00344   return res_mat;
00345 }
00347 //
00348 std::vector< std::vector<G4double>* > G4VEmAdjointModel::ComputeAdjointCrossSectionVectorPerVolumeForSecond(      
00349                                 G4Material* aMaterial,
00350                                 G4double kinEnergyProd,
00351                                 G4int nbin_pro_decade) //nb bins pro order of magnitude of energy
00352 { G4Integrator<G4VEmAdjointModel, double(G4VEmAdjointModel::*)(double)> integral;
00353   SelectedMaterial= aMaterial;
00354   kinEnergyProdForIntegration = kinEnergyProd;
00355    //compute the vector of integrated cross sections
00356   //-------------------
00357   
00358   G4double minEProj= GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
00359   G4double maxEProj= GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
00360   G4double E1=minEProj;
00361   std::vector< double>*  log_ESec_vector = new  std::vector< double>();
00362   std::vector< double>*  log_Prob_vector = new  std::vector< double>();
00363   log_ESec_vector->clear();
00364   log_Prob_vector->clear();
00365   log_ESec_vector->push_back(std::log(E1));
00366   log_Prob_vector->push_back(-50.);
00367   
00368   G4double E2=std::pow(10.,double( int(std::log10(minEProj)*nbin_pro_decade)+1)/nbin_pro_decade);
00369   G4double fE=std::pow(10.,1./nbin_pro_decade);
00370   G4double int_cross_section=0.;
00371   
00372   if (std::pow(fE,5.)>(maxEProj/minEProj)) fE = std::pow(maxEProj/minEProj,0.2);
00373   
00374   while (E1 <maxEProj*0.9999999){
00375         
00376         int_cross_section +=integral.Simpson(this,
00377         &G4VEmAdjointModel::DiffCrossSectionFunction1,E1,std::min(E2,maxEProj*0.99999999), 5);
00378         log_ESec_vector->push_back(std::log(std::min(E2,maxEProj)));
00379         log_Prob_vector->push_back(std::log(int_cross_section));
00380         E1=E2;
00381         E2*=fE;
00382   
00383   }
00384   std::vector< std::vector<G4double>* > res_mat;
00385   res_mat.clear();
00386   
00387   if (int_cross_section >0.) {
00388         res_mat.push_back(log_ESec_vector);
00389         res_mat.push_back(log_Prob_vector);
00390   }
00391   
00392  
00393   
00394   return res_mat;
00395 }
00396 
00398 //
00399 std::vector< std::vector<G4double>* > G4VEmAdjointModel::ComputeAdjointCrossSectionVectorPerVolumeForScatProj(
00400                                 G4Material* aMaterial,
00401                                 G4double kinEnergyScatProj,
00402                                 G4int nbin_pro_decade) //nb bins pro order of magnitude of energy
00403 { G4Integrator<G4VEmAdjointModel, double(G4VEmAdjointModel::*)(double)> integral;
00404   SelectedMaterial= aMaterial;
00405   kinEnergyScatProjForIntegration = kinEnergyScatProj;
00406  
00407   //compute the vector of integrated cross sections
00408   //-------------------
00409   
00410   G4double minEProj= GetSecondAdjEnergyMinForScatProjToProjCase(kinEnergyScatProj);
00411   G4double maxEProj= GetSecondAdjEnergyMaxForScatProjToProjCase(kinEnergyScatProj);
00412  
00413   
00414   G4double dEmax=maxEProj-kinEnergyScatProj;
00415   G4double dEmin=GetLowEnergyLimit();
00416   G4double dE1=dEmin;
00417   G4double dE2=dEmin;
00418   
00419   
00420   std::vector< double>*  log_ESec_vector = new std::vector< double>();
00421   std::vector< double>*  log_Prob_vector = new std::vector< double>();
00422   log_ESec_vector->push_back(std::log(dEmin));
00423   log_Prob_vector->push_back(-50.);
00424   G4int nbins=std::max( int(std::log10(dEmax/dEmin))*nbin_pro_decade,5);
00425   G4double fE=std::pow(dEmax/dEmin,1./nbins);
00426   
00427   G4double int_cross_section=0.;
00428   
00429   while (dE1 <dEmax*0.9999999999999){
00430         dE2=dE1*fE;
00431         int_cross_section +=integral.Simpson(this,
00432         &G4VEmAdjointModel::DiffCrossSectionFunction2,minEProj+dE1,std::min(minEProj+dE2,maxEProj), 5);
00433         log_ESec_vector->push_back(std::log(std::min(dE2,maxEProj-minEProj)));
00434         log_Prob_vector->push_back(std::log(int_cross_section));
00435         dE1=dE2;
00436 
00437   }
00438   
00439   
00440   
00441   
00442   
00443   std::vector< std::vector<G4double> *> res_mat;
00444   res_mat.clear();
00445   if (int_cross_section >0.) {
00446         res_mat.push_back(log_ESec_vector);
00447         res_mat.push_back(log_Prob_vector);
00448   }     
00449   
00450   return res_mat;
00451 }
00453 //                              
00454 G4double G4VEmAdjointModel::SampleAdjSecEnergyFromCSMatrix(size_t MatrixIndex,G4double aPrimEnergy,G4bool IsScatProjToProjCase)
00455 { 
00456   
00457   
00458   G4AdjointCSMatrix* theMatrix= (*pOnCSMatrixForProdToProjBackwardScattering)[MatrixIndex];
00459   if (IsScatProjToProjCase) theMatrix= (*pOnCSMatrixForScatProjToProjBackwardScattering)[MatrixIndex];
00460   std::vector< double>* theLogPrimEnergyVector = theMatrix->GetLogPrimEnergyVector();
00461   
00462   if (theLogPrimEnergyVector->size() ==0){
00463         G4cout<<"No data are contained in the given AdjointCSMatrix!"<<G4endl;
00464         G4cout<<"The sampling procedure will be stopped."<<G4endl;
00465         return 0.;
00466         
00467   }
00468   
00469   G4AdjointInterpolator* theInterpolator=G4AdjointInterpolator::GetInstance();
00470   G4double aLogPrimEnergy = std::log(aPrimEnergy);
00471   size_t ind =theInterpolator->FindPositionForLogVector(aLogPrimEnergy,*theLogPrimEnergyVector);
00472   
00473   
00474   G4double aLogPrimEnergy1,aLogPrimEnergy2;
00475   G4double aLogCS1,aLogCS2;
00476   G4double log01,log02;
00477   std::vector< double>* aLogSecondEnergyVector1 =0;
00478   std::vector< double>* aLogSecondEnergyVector2  =0;
00479   std::vector< double>* aLogProbVector1=0;
00480   std::vector< double>* aLogProbVector2=0; 
00481   std::vector< size_t>* aLogProbVectorIndex1=0;
00482   std::vector< size_t>* aLogProbVectorIndex2=0;
00483                                                                      
00484   theMatrix->GetData(ind, aLogPrimEnergy1,aLogCS1,log01, aLogSecondEnergyVector1,aLogProbVector1,aLogProbVectorIndex1);
00485   theMatrix->GetData(ind+1, aLogPrimEnergy2,aLogCS2,log02, aLogSecondEnergyVector2,aLogProbVector2,aLogProbVectorIndex2);
00486   
00487   G4double rand_var = G4UniformRand();
00488   G4double log_rand_var= std::log(rand_var);
00489   G4double log_Tcut =std::log(currentTcutForDirectSecond);
00490   G4double Esec=0;
00491   G4double log_dE1,log_dE2;
00492   G4double log_rand_var1,log_rand_var2;
00493   G4double log_E1,log_E2;
00494   log_rand_var1=log_rand_var;
00495   log_rand_var2=log_rand_var;
00496   
00497   G4double Emin=0.;
00498   G4double Emax=0.;
00499   if (theMatrix->IsScatProjToProjCase()){ //case where Tcut plays a role
00500         Emin=GetSecondAdjEnergyMinForScatProjToProjCase(aPrimEnergy,currentTcutForDirectSecond);
00501         Emax=GetSecondAdjEnergyMaxForScatProjToProjCase(aPrimEnergy);
00502         G4double dE=0;
00503         if (Emin < Emax ){
00504             if (ApplyCutInRange) {
00505                 if (second_part_of_same_type && currentTcutForDirectSecond>aPrimEnergy) return aPrimEnergy;
00506                 
00507                 log_rand_var1=log_rand_var+theInterpolator->InterpolateForLogVector(log_Tcut,*aLogSecondEnergyVector1,*aLogProbVector1);
00508                 log_rand_var2=log_rand_var+theInterpolator->InterpolateForLogVector(log_Tcut,*aLogSecondEnergyVector2,*aLogProbVector2);
00509                 
00510             }   
00511             log_dE1 = theInterpolator->Interpolate(log_rand_var1,*aLogProbVector1,*aLogSecondEnergyVector1,"Lin");
00512             log_dE2 = theInterpolator->Interpolate(log_rand_var2,*aLogProbVector2,*aLogSecondEnergyVector2,"Lin");
00513              dE=std::exp(theInterpolator->LinearInterpolation(aLogPrimEnergy,aLogPrimEnergy1,aLogPrimEnergy2,log_dE1,log_dE2));
00514         }
00515         
00516         Esec = aPrimEnergy +dE;
00517         Esec=std::max(Esec,Emin);
00518         Esec=std::min(Esec,Emax);
00519         
00520   }
00521   else { //Tcut condition is already full-filled
00522         
00523         log_E1 = theInterpolator->Interpolate(log_rand_var,*aLogProbVector1,*aLogSecondEnergyVector1,"Lin");
00524         log_E2 = theInterpolator->Interpolate(log_rand_var,*aLogProbVector2,*aLogSecondEnergyVector2,"Lin");
00525         
00526         Esec = std::exp(theInterpolator->LinearInterpolation(aLogPrimEnergy,aLogPrimEnergy1,aLogPrimEnergy2,log_E1,log_E2));
00527         Emin=GetSecondAdjEnergyMinForProdToProjCase(aPrimEnergy);
00528         Emax=GetSecondAdjEnergyMaxForProdToProjCase(aPrimEnergy);
00529         Esec=std::max(Esec,Emin);
00530         Esec=std::min(Esec,Emax);
00531         
00532   }
00533         
00534   return Esec;
00535   
00536   
00537   
00538  
00539                                                                                    
00540 }
00541 
00543 //                              
00544 G4double G4VEmAdjointModel::SampleAdjSecEnergyFromCSMatrix(G4double aPrimEnergy,G4bool IsScatProjToProjCase)
00545 { SelectCSMatrix(IsScatProjToProjCase);
00546   return SampleAdjSecEnergyFromCSMatrix(indexOfUsedCrossSectionMatrix, aPrimEnergy, IsScatProjToProjCase);
00547 }
00549 //
00550 void G4VEmAdjointModel::SelectCSMatrix(G4bool IsScatProjToProjCase)
00551 { 
00552   indexOfUsedCrossSectionMatrix=0;
00553   if (!UseMatrixPerElement) indexOfUsedCrossSectionMatrix = currentMaterialIndex;
00554   else if (!UseOnlyOneMatrixForAllElements) { //Select Material
00555         std::vector<G4double>* CS_Vs_Element = &CS_Vs_ElementForScatProjToProjCase;
00556         lastCS=lastAdjointCSForScatProjToProjCase;
00557         if ( !IsScatProjToProjCase) {
00558                 CS_Vs_Element = &CS_Vs_ElementForProdToProjCase;
00559                 lastCS=lastAdjointCSForProdToProjCase;
00560         }       
00561         G4double rand_var= G4UniformRand();
00562         G4double SumCS=0.;
00563         size_t ind=0;
00564         for (size_t i=0;i<CS_Vs_Element->size();i++){
00565                 SumCS+=(*CS_Vs_Element)[i];
00566                 if (rand_var<=SumCS/lastCS){
00567                         ind=i;
00568                         break;
00569                 }
00570         }
00571         indexOfUsedCrossSectionMatrix = currentMaterial->GetElement(ind)->GetIndex();
00572   }
00573 }
00575 //                              
00576 G4double G4VEmAdjointModel::SampleAdjSecEnergyFromDiffCrossSectionPerAtom(G4double prim_energy,G4bool IsScatProjToProjCase)
00577 {  
00578   // here we try to use the rejection method 
00579   //-----------------------------------------
00580   
00581   G4double E=0;
00582   G4double x,xmin,greject,q;
00583   if ( IsScatProjToProjCase){
00584         G4double Emax = GetSecondAdjEnergyMaxForScatProjToProjCase(prim_energy);
00585         G4double Emin= prim_energy+currentTcutForDirectSecond;
00586         xmin=Emin/Emax;
00587         G4double grejmax = DiffCrossSectionPerAtomPrimToScatPrim(Emin,prim_energy,1)*prim_energy; 
00588 
00589         do {
00590                 q = G4UniformRand();
00591                 x = 1./(q*(1./xmin -1.) +1.);
00592                 E=x*Emax;
00593                 greject = DiffCrossSectionPerAtomPrimToScatPrim( E,prim_energy ,1)*prim_energy; 
00594                 
00595         }       
00596         
00597         while( greject < G4UniformRand()*grejmax );     
00598         
00599   }
00600   else {
00601         G4double Emax = GetSecondAdjEnergyMaxForProdToProjCase(prim_energy);
00602         G4double Emin=  GetSecondAdjEnergyMinForProdToProjCase(prim_energy);;
00603         xmin=Emin/Emax;
00604         G4double grejmax = DiffCrossSectionPerAtomPrimToSecond(Emin,prim_energy,1); 
00605         do {
00606                 q = G4UniformRand();
00607                 x = std::pow(xmin, q);
00608                 E=x*Emax;
00609                 greject = DiffCrossSectionPerAtomPrimToSecond( E,prim_energy ,1); 
00610                 
00611         }       
00612         
00613         while( greject < G4UniformRand()*grejmax );
00614   
00615   
00616   
00617   }
00618   
00619   return E;
00620 }
00621 
00623 //
00624 void G4VEmAdjointModel::CorrectPostStepWeight(G4ParticleChange* fParticleChange, 
00625                                                   G4double old_weight,  
00626                                                   G4double adjointPrimKinEnergy, 
00627                                                   G4double projectileKinEnergy,
00628                                                   G4bool IsScatProjToProjCase) 
00629 {
00630  G4double new_weight=old_weight;
00631  G4double w_corr =1./CS_biasing_factor;
00632  w_corr*=G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection();
00633  
00634  
00635  lastCS=lastAdjointCSForScatProjToProjCase;
00636  if ( !IsScatProjToProjCase) lastCS=lastAdjointCSForProdToProjCase;
00637  if ((adjointPrimKinEnergy-preStepEnergy)/preStepEnergy>0.001){ //Is that in all cases needed???
00638         G4double post_stepCS=AdjointCrossSection(currentCouple, adjointPrimKinEnergy
00639                                                  ,IsScatProjToProjCase );
00640         if (post_stepCS>0 && lastCS>0) w_corr*=post_stepCS/lastCS;
00641  }
00642         
00643  new_weight*=w_corr;
00644 
00645  //G4cout<<"Post step "<<new_weight<<'\t'<<w_corr<<'\t'<<old_weight<<G4endl;
00646  new_weight*=projectileKinEnergy/adjointPrimKinEnergy;//This is needed due to the biasing of diff CS
00647                                                         //by the factor adjointPrimKinEnergy/projectileKinEnergy
00648    
00649 
00650 
00651  fParticleChange->SetParentWeightByProcess(false);
00652  fParticleChange->SetSecondaryWeightByProcess(false);
00653  fParticleChange->ProposeParentWeight(new_weight);      
00654 }
00656 //                              
00657 G4double G4VEmAdjointModel::GetSecondAdjEnergyMaxForScatProjToProjCase(G4double kinEnergyScatProj)
00658 { G4double maxEProj= HighEnergyLimit;
00659   if (second_part_of_same_type)  maxEProj=std::min(kinEnergyScatProj*2.,HighEnergyLimit);
00660   return maxEProj;
00661 }
00663 //
00664 G4double G4VEmAdjointModel::GetSecondAdjEnergyMinForScatProjToProjCase(G4double PrimAdjEnergy,G4double Tcut)
00665 { G4double Emin=PrimAdjEnergy;
00666   if (ApplyCutInRange) Emin=PrimAdjEnergy+Tcut;
00667   return Emin;
00668 }
00670 //                              
00671 G4double G4VEmAdjointModel::GetSecondAdjEnergyMaxForProdToProjCase(G4double )
00672 { return HighEnergyLimit;
00673 }
00675 //
00676 G4double G4VEmAdjointModel::GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy)
00677 { G4double minEProj=PrimAdjEnergy;
00678   if (second_part_of_same_type)  minEProj=PrimAdjEnergy*2.;
00679   return minEProj;
00680 }
00682 //
00683 void  G4VEmAdjointModel::DefineCurrentMaterial(const G4MaterialCutsCouple* couple)
00684 { if(couple != currentCouple) {
00685         currentCouple   = const_cast<G4MaterialCutsCouple*> (couple);
00686         currentMaterial = const_cast<G4Material*> (couple->GetMaterial());
00687         currentCoupleIndex = couple->GetIndex();
00688         currentMaterialIndex = currentMaterial->GetIndex();
00689         size_t idx=56;
00690         currentTcutForDirectSecond =0.00000000001;
00691         if (theAdjEquivOfDirectSecondPartDef) {
00692                 if (theAdjEquivOfDirectSecondPartDef == G4AdjointGamma::AdjointGamma()) idx = 0;
00693                 else if (theAdjEquivOfDirectSecondPartDef == G4AdjointElectron::AdjointElectron()) idx = 1;
00694                 else if (theAdjEquivOfDirectSecondPartDef == G4AdjointPositron::AdjointPositron()) idx = 2;
00695                 if (idx <56){
00696                         const std::vector<G4double>* aVec = G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(idx);
00697                         currentTcutForDirectSecond=(*aVec)[currentCoupleIndex];
00698                 }       
00699         }       
00700     
00701         
00702   }
00703 }
00705 //
00706 void G4VEmAdjointModel::SetHighEnergyLimit(G4double aVal)
00707 { HighEnergyLimit=aVal;
00708   if (theDirectEMModel) theDirectEMModel->SetHighEnergyLimit( aVal);
00709 }
00711 //
00712 void G4VEmAdjointModel::SetLowEnergyLimit(G4double aVal)
00713 {
00714   LowEnergyLimit=aVal;
00715   if (theDirectEMModel) theDirectEMModel->SetLowEnergyLimit( aVal);
00716 }
00718 //
00719 void G4VEmAdjointModel::SetAdjointEquivalentOfDirectPrimaryParticleDefinition(G4ParticleDefinition* aPart)
00720 {
00721   theAdjEquivOfDirectPrimPartDef=aPart;
00722   if (theAdjEquivOfDirectPrimPartDef->GetParticleName() =="adj_e-")
00723                                         theDirectPrimaryPartDef=G4Electron::Electron();
00724   if (theAdjEquivOfDirectPrimPartDef->GetParticleName() =="adj_gamma")
00725                                         theDirectPrimaryPartDef=G4Gamma::Gamma();
00726 }

Generated on Mon May 27 17:50:11 2013 for Geant4 by  doxygen 1.4.7