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 "G4AdjointBremsstrahlungModel.hh"
00029 #include "G4AdjointCSManager.hh"
00030
00031 #include "G4PhysicalConstants.hh"
00032 #include "G4SystemOfUnits.hh"
00033
00034 #include "G4Integrator.hh"
00035 #include "G4TrackStatus.hh"
00036 #include "G4ParticleChange.hh"
00037 #include "G4AdjointElectron.hh"
00038 #include "G4AdjointGamma.hh"
00039 #include "G4Electron.hh"
00040 #include "G4Timer.hh"
00041 #include "G4SeltzerBergerModel.hh"
00042
00043
00045
00046 G4AdjointBremsstrahlungModel::G4AdjointBremsstrahlungModel(G4VEmModel* aModel):
00047 G4VEmAdjointModel("AdjointeBremModel")
00048 {
00049 SetUseMatrix(false);
00050 SetUseMatrixPerElement(false);
00051
00052 theDirectStdBremModel = aModel;
00053 theDirectEMModel=theDirectStdBremModel;
00054 theEmModelManagerForFwdModels = new G4EmModelManager();
00055 isDirectModelInitialised = false;
00056 G4VEmFluctuationModel* f=0;
00057 G4Region* r=0;
00058 theEmModelManagerForFwdModels->AddEmModel(1, theDirectStdBremModel, f, r);
00059
00060 SetApplyCutInRange(true);
00061 highKinEnergy= 100.*TeV;
00062 lowKinEnergy = 1.0*keV;
00063
00064 lastCZ =0.;
00065
00066
00067 theAdjEquivOfDirectPrimPartDef =G4AdjointElectron::AdjointElectron();
00068 theAdjEquivOfDirectSecondPartDef=G4AdjointGamma::AdjointGamma();
00069 theDirectPrimaryPartDef=G4Electron::Electron();
00070 second_part_of_same_type=false;
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086 }
00088
00089 G4AdjointBremsstrahlungModel::G4AdjointBremsstrahlungModel():
00090 G4VEmAdjointModel("AdjointeBremModel")
00091 {
00092 SetUseMatrix(false);
00093 SetUseMatrixPerElement(false);
00094
00095 theDirectStdBremModel = new G4SeltzerBergerModel();
00096 theDirectEMModel=theDirectStdBremModel;
00097 theEmModelManagerForFwdModels = new G4EmModelManager();
00098 isDirectModelInitialised = false;
00099 G4VEmFluctuationModel* f=0;
00100 G4Region* r=0;
00101 theEmModelManagerForFwdModels->AddEmModel(1, theDirectStdBremModel, f, r);
00102
00103 SetApplyCutInRange(true);
00104 highKinEnergy= 1.*GeV;
00105 lowKinEnergy = 1.0*keV;
00106 lastCZ =0.;
00107 theAdjEquivOfDirectPrimPartDef =G4AdjointElectron::AdjointElectron();
00108 theAdjEquivOfDirectSecondPartDef=G4AdjointGamma::AdjointGamma();
00109 theDirectPrimaryPartDef=G4Electron::Electron();
00110 second_part_of_same_type=false;
00111
00112 }
00114
00115 G4AdjointBremsstrahlungModel::~G4AdjointBremsstrahlungModel()
00116 {if (theDirectStdBremModel) delete theDirectStdBremModel;
00117 if (theEmModelManagerForFwdModels) delete theEmModelManagerForFwdModels;
00118 }
00119
00121
00122 void G4AdjointBremsstrahlungModel::SampleSecondaries(const G4Track& aTrack,
00123 G4bool IsScatProjToProjCase,
00124 G4ParticleChange* fParticleChange)
00125 {
00126 if (!UseMatrix) return RapidSampleSecondaries(aTrack,IsScatProjToProjCase,fParticleChange);
00127
00128 const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
00129 DefineCurrentMaterial(aTrack.GetMaterialCutsCouple());
00130
00131
00132 G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
00133 G4double adjointPrimTotalEnergy = theAdjointPrimary->GetTotalEnergy();
00134
00135 if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
00136 return;
00137 }
00138
00139 G4double projectileKinEnergy = SampleAdjSecEnergyFromCSMatrix(adjointPrimKinEnergy,
00140 IsScatProjToProjCase);
00141
00142
00143 CorrectPostStepWeight(fParticleChange,
00144 aTrack.GetWeight(),
00145 adjointPrimKinEnergy,
00146 projectileKinEnergy,
00147 IsScatProjToProjCase);
00148
00149
00150
00151
00152 G4double projectileM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
00153 G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
00154 G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0;
00155 G4double projectileP = std::sqrt(projectileP2);
00156
00157
00158
00159
00160 G4double u;
00161 const G4double a1 = 0.625 , a2 = 3.*a1 , d = 27. ;
00162
00163 if (9./(9.+d) > G4UniformRand()) u = - std::log(G4UniformRand()*G4UniformRand())/a1;
00164 else u = - std::log(G4UniformRand()*G4UniformRand())/a2;
00165
00166 G4double theta = u*electron_mass_c2/projectileTotalEnergy;
00167
00168 G4double sint = std::sin(theta);
00169 G4double cost = std::cos(theta);
00170
00171 G4double phi = twopi * G4UniformRand() ;
00172
00173 G4ThreeVector projectileMomentum;
00174 projectileMomentum=G4ThreeVector(std::cos(phi)*sint,std::sin(phi)*sint,cost)*projectileP;
00175 if (IsScatProjToProjCase) {
00176 G4ThreeVector gammaMomentum = (projectileTotalEnergy-adjointPrimTotalEnergy)*G4ThreeVector(0.,0.,1.);
00177 G4ThreeVector dirProd=projectileMomentum-gammaMomentum;
00178 G4double cost1 = std::cos(dirProd.angle(projectileMomentum));
00179 G4double sint1 = std::sqrt(1.-cost1*cost1);
00180 projectileMomentum=G4ThreeVector(std::cos(phi)*sint1,std::sin(phi)*sint1,cost1)*projectileP;
00181
00182 }
00183
00184 projectileMomentum.rotateUz(theAdjointPrimary->GetMomentumDirection());
00185
00186
00187
00188 if (!IsScatProjToProjCase ){
00189 fParticleChange->ProposeTrackStatus(fStopAndKill);
00190 fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum));
00191 }
00192 else {
00193 fParticleChange->ProposeEnergy(projectileKinEnergy);
00194 fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
00195
00196 }
00197 }
00199
00200 void G4AdjointBremsstrahlungModel::RapidSampleSecondaries(const G4Track& aTrack,
00201 G4bool IsScatProjToProjCase,
00202 G4ParticleChange* fParticleChange)
00203 {
00204
00205 const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
00206 DefineCurrentMaterial(aTrack.GetMaterialCutsCouple());
00207
00208
00209 G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
00210 G4double adjointPrimTotalEnergy = theAdjointPrimary->GetTotalEnergy();
00211
00212 if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
00213 return;
00214 }
00215
00216 G4double projectileKinEnergy =0.;
00217 G4double gammaEnergy=0.;
00218 G4double diffCSUsed=0.;
00219 if (!IsScatProjToProjCase){
00220 gammaEnergy=adjointPrimKinEnergy;
00221 G4double Emax = GetSecondAdjEnergyMaxForProdToProjCase(adjointPrimKinEnergy);
00222 G4double Emin= GetSecondAdjEnergyMinForProdToProjCase(adjointPrimKinEnergy);;
00223 if (Emin>=Emax) return;
00224 projectileKinEnergy=Emin*std::pow(Emax/Emin,G4UniformRand());
00225 diffCSUsed=lastCZ/projectileKinEnergy;
00226
00227 }
00228 else { G4double Emax = GetSecondAdjEnergyMaxForScatProjToProjCase(adjointPrimKinEnergy);
00229 G4double Emin = GetSecondAdjEnergyMinForScatProjToProjCase(adjointPrimKinEnergy,currentTcutForDirectSecond);
00230 if (Emin>=Emax) return;
00231 G4double f1=(Emin-adjointPrimKinEnergy)/Emin;
00232 G4double f2=(Emax-adjointPrimKinEnergy)/Emax/f1;
00233
00234 projectileKinEnergy=adjointPrimKinEnergy/(1.-f1*std::pow(f2,G4UniformRand()));
00235 gammaEnergy=projectileKinEnergy-adjointPrimKinEnergy;
00236 diffCSUsed=lastCZ*adjointPrimKinEnergy/projectileKinEnergy/gammaEnergy;
00237
00238 }
00239
00240
00241
00242
00243
00244
00245
00246 G4double w_corr=G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection();
00247
00248
00249
00250
00251
00252 G4double diffCS = DiffCrossSectionPerVolumePrimToSecond(currentMaterial, projectileKinEnergy, gammaEnergy);
00253 w_corr*=diffCS/diffCSUsed;
00254
00255 G4double new_weight = aTrack.GetWeight()*w_corr;
00256 fParticleChange->SetParentWeightByProcess(false);
00257 fParticleChange->SetSecondaryWeightByProcess(false);
00258 fParticleChange->ProposeParentWeight(new_weight);
00259
00260
00261
00262 G4double projectileM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
00263 G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
00264 G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0;
00265 G4double projectileP = std::sqrt(projectileP2);
00266
00267
00268
00269
00270 G4double u;
00271 const G4double a1 = 0.625 , a2 = 3.*a1 , d = 27. ;
00272
00273 if (9./(9.+d) > G4UniformRand()) u = - std::log(G4UniformRand()*G4UniformRand())/a1;
00274 else u = - std::log(G4UniformRand()*G4UniformRand())/a2;
00275
00276 G4double theta = u*electron_mass_c2/projectileTotalEnergy;
00277
00278 G4double sint = std::sin(theta);
00279 G4double cost = std::cos(theta);
00280
00281 G4double phi = twopi * G4UniformRand() ;
00282
00283 G4ThreeVector projectileMomentum;
00284 projectileMomentum=G4ThreeVector(std::cos(phi)*sint,std::sin(phi)*sint,cost)*projectileP;
00285 if (IsScatProjToProjCase) {
00286 G4ThreeVector gammaMomentum = (projectileTotalEnergy-adjointPrimTotalEnergy)*G4ThreeVector(0.,0.,1.);
00287 G4ThreeVector dirProd=projectileMomentum-gammaMomentum;
00288 G4double cost1 = std::cos(dirProd.angle(projectileMomentum));
00289 G4double sint1 = std::sqrt(1.-cost1*cost1);
00290 projectileMomentum=G4ThreeVector(std::cos(phi)*sint1,std::sin(phi)*sint1,cost1)*projectileP;
00291
00292 }
00293
00294 projectileMomentum.rotateUz(theAdjointPrimary->GetMomentumDirection());
00295
00296
00297
00298 if (!IsScatProjToProjCase ){
00299 fParticleChange->ProposeTrackStatus(fStopAndKill);
00300 fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum));
00301 }
00302 else {
00303 fParticleChange->ProposeEnergy(projectileKinEnergy);
00304 fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
00305
00306 }
00307 }
00309
00310 G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecond(const G4Material* aMaterial,
00311 G4double kinEnergyProj,
00312 G4double kinEnergyProd
00313 )
00314 {if (!isDirectModelInitialised) {
00315 theEmModelManagerForFwdModels->Initialise(G4Electron::Electron(),G4Gamma::Gamma(),1.,0);
00316 isDirectModelInitialised =true;
00317 }
00318
00319 return DiffCrossSectionPerVolumePrimToSecondApproximated2(aMaterial,
00320 kinEnergyProj,
00321 kinEnergyProd);
00322
00323
00324
00325 }
00326
00328
00329 G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecondApproximated1(
00330 const G4Material* aMaterial,
00331 G4double kinEnergyProj,
00332 G4double kinEnergyProd
00333 )
00334 {
00335 G4double dCrossEprod=0.;
00336 G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
00337 G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348 if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){
00349 G4double sigma=theDirectEMModel->CrossSectionPerVolume(aMaterial,theDirectPrimaryPartDef,kinEnergyProj,1.*keV);
00350 dCrossEprod=sigma/kinEnergyProd/std::log(kinEnergyProj/keV);
00351 }
00352 return dCrossEprod;
00353
00354 }
00355
00357
00358 G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecondApproximated2(
00359 const G4Material* material,
00360 G4double kinEnergyProj,
00361 G4double kinEnergyProd
00362 )
00363 {
00364
00365
00366
00367 G4double dCrossEprod=0.;
00368
00369 const G4ElementVector* theElementVector = material->GetElementVector();
00370 const double* theAtomNumDensityVector = material->GetAtomicNumDensityVector();
00371 G4double dum=0.;
00372 G4double E1=kinEnergyProd,E2=kinEnergyProd*1.001;
00373 G4double dE=E2-E1;
00374 for (size_t i=0; i<material->GetNumberOfElements(); i++) {
00375 G4double C1=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,(*theElementVector)[i]->GetZ(),dum ,E1);
00376 G4double C2=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,(*theElementVector)[i]->GetZ(),dum,E2);
00377 dCrossEprod += theAtomNumDensityVector[i] * (C1-C2)/dE;
00378
00379 }
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394 return dCrossEprod;
00395
00396 }
00398
00399 G4double G4AdjointBremsstrahlungModel::AdjointCrossSection(const G4MaterialCutsCouple* aCouple,
00400 G4double primEnergy,
00401 G4bool IsScatProjToProjCase)
00402 { if (!isDirectModelInitialised) {
00403 theEmModelManagerForFwdModels->Initialise(G4Electron::Electron(),G4Gamma::Gamma(),1.,0);
00404 isDirectModelInitialised =true;
00405 }
00406 if (UseMatrix) return G4VEmAdjointModel::AdjointCrossSection(aCouple,primEnergy,IsScatProjToProjCase);
00407 DefineCurrentMaterial(aCouple);
00408 G4double Cross=0.;
00409 lastCZ=theDirectEMModel->CrossSectionPerVolume(aCouple->GetMaterial(),theDirectPrimaryPartDef,100.*MeV,100.*MeV/std::exp(1.));
00410
00411 if (!IsScatProjToProjCase ){
00412 G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(primEnergy);
00413 G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(primEnergy);
00414 if (Emax_proj>Emin_proj && primEnergy > currentTcutForDirectSecond) Cross= lastCZ*std::log(Emax_proj/Emin_proj);
00415 }
00416 else {
00417 G4double Emax_proj = GetSecondAdjEnergyMaxForScatProjToProjCase(primEnergy);
00418 G4double Emin_proj = GetSecondAdjEnergyMinForScatProjToProjCase(primEnergy,currentTcutForDirectSecond);
00419 if (Emax_proj>Emin_proj) Cross= lastCZ*std::log((Emax_proj-primEnergy)*Emin_proj/Emax_proj/(Emin_proj-primEnergy));
00420
00421 }
00422 return Cross;
00423 }
00424
00425 G4double G4AdjointBremsstrahlungModel::GetAdjointCrossSection(const G4MaterialCutsCouple* aCouple,
00426 G4double primEnergy,
00427 G4bool IsScatProjToProjCase)
00428 {
00429 return AdjointCrossSection(aCouple, primEnergy,IsScatProjToProjCase);
00430 lastCZ=theDirectEMModel->CrossSectionPerVolume(aCouple->GetMaterial(),theDirectPrimaryPartDef,100.*MeV,100.*MeV/std::exp(1.));
00431 return G4VEmAdjointModel::GetAdjointCrossSection(aCouple, primEnergy,IsScatProjToProjCase);
00432
00433 }
00434
00435
00436