00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
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
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
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108 }
00110
00111
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){
00124
00125
00126
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
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
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
00172
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
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)
00243 {
00244 G4Integrator<G4VEmAdjointModel, double(G4VEmAdjointModel::*)(double)> integral;
00245 ASelectedNucleus= int(A);
00246 ZSelectedNucleus=int(Z);
00247 kinEnergyProdForIntegration = kinEnergyProd;
00248
00249
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
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)
00296 { G4Integrator<G4VEmAdjointModel, double(G4VEmAdjointModel::*)(double)> integral;
00297 ASelectedNucleus=int(A);
00298 ZSelectedNucleus=int(Z);
00299 kinEnergyScatProjForIntegration = kinEnergyScatProj;
00300
00301
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
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)
00352 { G4Integrator<G4VEmAdjointModel, double(G4VEmAdjointModel::*)(double)> integral;
00353 SelectedMaterial= aMaterial;
00354 kinEnergyProdForIntegration = kinEnergyProd;
00355
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)
00403 { G4Integrator<G4VEmAdjointModel, double(G4VEmAdjointModel::*)(double)> integral;
00404 SelectedMaterial= aMaterial;
00405 kinEnergyScatProjForIntegration = kinEnergyScatProj;
00406
00407
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()){
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 {
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) {
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
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){
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
00646 new_weight*=projectileKinEnergy/adjointPrimKinEnergy;
00647
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 }