// // ******************************************************************** // * License and Disclaimer * // * * // * The Geant4 software is copyright of the Copyright Holders of * // * the Geant4 Collaboration. It is provided under the terms and * // * conditions of the Geant4 Software License, included in the file * // * LICENSE and available at http://cern.ch/geant4/license . These * // * include a list of copyright holders. * // * * // * Neither the authors of this software system, nor their employing * // * institutes,nor the agencies providing financial support for this * // * work make any representation or warranty, express or implied, * // * regarding this software system or assume any liability for its * // * use. Please see the license in the file LICENSE and URL above * // * for the full disclaimer and the limitation of liability. * // * * // * This code implementation is the result of the scientific and * // * technical work of the GEANT4 collaboration. * // * By using, copying, modifying or distributing the software (or * // * any work based on the software) you agree to acknowledge its * // * use in resulting scientific publications, and indicate your * // * acceptance of all terms of the Geant4 Software license. * // ******************************************************************** // // $Id: G4AdjointBremsstrahlungModel.cc,v 1.5 2009/12/16 17:50:01 gunter Exp $ // GEANT4 tag $Name: geant4-09-04-beta-01 $ // #include "G4AdjointBremsstrahlungModel.hh" #include "G4AdjointCSManager.hh" #include "G4Integrator.hh" #include "G4TrackStatus.hh" #include "G4ParticleChange.hh" #include "G4AdjointElectron.hh" #include "G4AdjointGamma.hh" #include "G4Electron.hh" #include "G4Timer.hh" //#include "G4PenelopeBremsstrahlungModel.hh" //////////////////////////////////////////////////////////////////////////////// // G4AdjointBremsstrahlungModel::G4AdjointBremsstrahlungModel(): G4VEmAdjointModel("AdjointeBremModel"), MigdalConstant(classic_electr_radius*electron_Compton_length*electron_Compton_length*4.0*pi) { SetUseMatrix(false); SetUseMatrixPerElement(false); theDirectStdBremModel = new G4eBremsstrahlungModel(G4Electron::Electron(),"TheDirecteBremModel"); theDirectEMModel=theDirectStdBremModel; // theDirectPenelopeBremModel =0; SetApplyCutInRange(true); highKinEnergy= 100.*TeV; lowKinEnergy = 1.0*keV; theTimer =new G4Timer(); theAdjEquivOfDirectPrimPartDef =G4AdjointElectron::AdjointElectron(); theAdjEquivOfDirectSecondPartDef=G4AdjointGamma::AdjointGamma(); theDirectPrimaryPartDef=G4Electron::Electron(); second_part_of_same_type=false; /*UsePenelopeModel=false; if (UsePenelopeModel) { G4PenelopeBremsstrahlungModel* thePenelopeModel = new G4PenelopeBremsstrahlungModel(G4Electron::Electron(),"PenelopeBrem"); theEmModelManagerForFwdModels = new G4EmModelManager(); isPenelopeModelInitialised = false; G4VEmFluctuationModel* f=0; G4Region* r=0; theDirectEMModel=thePenelopeModel; theEmModelManagerForFwdModels->AddEmModel(1, thePenelopeModel, f, r); } */ } //////////////////////////////////////////////////////////////////////////////// // void G4AdjointBremsstrahlungModel::SampleSecondaries(const G4Track& aTrack, G4bool IsScatProjToProjCase, G4ParticleChange* fParticleChange) { if (!UseMatrix) return RapidSampleSecondaries(aTrack,IsScatProjToProjCase,fParticleChange); const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle(); DefineCurrentMaterial(aTrack.GetMaterialCutsCouple()); G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy(); G4double adjointPrimTotalEnergy = theAdjointPrimary->GetTotalEnergy(); if (adjointPrimKinEnergy>HighEnergyLimit*0.999){ return; } G4double projectileKinEnergy = SampleAdjSecEnergyFromCSMatrix(adjointPrimKinEnergy, IsScatProjToProjCase); //Weight correction //----------------------- CorrectPostStepWeight(fParticleChange, aTrack.GetWeight(), adjointPrimKinEnergy, projectileKinEnergy, IsScatProjToProjCase); //Kinematic //--------- G4double projectileM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass(); G4double projectileTotalEnergy = projectileM0+projectileKinEnergy; G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0; G4double projectileP = std::sqrt(projectileP2); //Angle of the gamma direction with the projectile taken from G4eBremsstrahlungModel //------------------------------------------------ G4double u; const G4double a1 = 0.625 , a2 = 3.*a1 , d = 27. ; if (9./(9.+d) > G4UniformRand()) u = - std::log(G4UniformRand()*G4UniformRand())/a1; else u = - std::log(G4UniformRand()*G4UniformRand())/a2; G4double theta = u*electron_mass_c2/projectileTotalEnergy; G4double sint = std::sin(theta); G4double cost = std::cos(theta); G4double phi = twopi * G4UniformRand() ; G4ThreeVector projectileMomentum; projectileMomentum=G4ThreeVector(std::cos(phi)*sint,std::sin(phi)*sint,cost)*projectileP; //gamma frame if (IsScatProjToProjCase) {//the adjoint primary is the scattered e- G4ThreeVector gammaMomentum = (projectileTotalEnergy-adjointPrimTotalEnergy)*G4ThreeVector(0.,0.,1.); G4ThreeVector dirProd=projectileMomentum-gammaMomentum; G4double cost1 = std::cos(dirProd.angle(projectileMomentum)); G4double sint1 = std::sqrt(1.-cost1*cost1); projectileMomentum=G4ThreeVector(std::cos(phi)*sint1,std::sin(phi)*sint1,cost1)*projectileP; } projectileMomentum.rotateUz(theAdjointPrimary->GetMomentumDirection()); if (!IsScatProjToProjCase ){ //kill the primary and add a secondary fParticleChange->ProposeTrackStatus(fStopAndKill); fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum)); } else { fParticleChange->ProposeEnergy(projectileKinEnergy); fParticleChange->ProposeMomentumDirection(projectileMomentum.unit()); } } //////////////////////////////////////////////////////////////////////////////// // void G4AdjointBremsstrahlungModel::RapidSampleSecondaries(const G4Track& aTrack, G4bool IsScatProjToProjCase, G4ParticleChange* fParticleChange) { const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle(); DefineCurrentMaterial(aTrack.GetMaterialCutsCouple()); G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy(); G4double adjointPrimTotalEnergy = theAdjointPrimary->GetTotalEnergy(); if (adjointPrimKinEnergy>HighEnergyLimit*0.999){ return; } G4double projectileKinEnergy =0.; G4double gammaEnergy=0.; G4double diffCSUsed=0.; if (!IsScatProjToProjCase){ gammaEnergy=adjointPrimKinEnergy; G4double Emax = GetSecondAdjEnergyMaxForProdToProjCase(adjointPrimKinEnergy); G4double Emin= GetSecondAdjEnergyMinForProdToProjCase(adjointPrimKinEnergy);; if (Emin>=Emax) return; projectileKinEnergy=Emin*std::pow(Emax/Emin,G4UniformRand()); diffCSUsed=lastCZ/projectileKinEnergy; } else { G4double Emax = GetSecondAdjEnergyMaxForScatProjToProjCase(adjointPrimKinEnergy); G4double Emin = GetSecondAdjEnergyMinForScatProjToProjCase(adjointPrimKinEnergy,currentTcutForDirectSecond); if (Emin>=Emax) return; G4double f1=(Emin-adjointPrimKinEnergy)/Emin; G4double f2=(Emax-adjointPrimKinEnergy)/Emax/f1; //G4cout<<"f1 and f2 "<GetPostStepWeightCorrection(); //Then another correction is needed due to the fact that a biaised differential CS has been used rather than the one consistent with the direct model //Here we consider the true diffCS as the one obtained by the numericla differentiation over Tcut of the direct CS, corrected by the Migdal term. //Basically any other differential CS diffCS could be used here (example Penelope). G4double diffCS = DiffCrossSectionPerVolumePrimToSecond(currentMaterial, projectileKinEnergy, gammaEnergy); w_corr*=diffCS/diffCSUsed; G4double new_weight = aTrack.GetWeight()*w_corr; fParticleChange->SetParentWeightByProcess(false); fParticleChange->SetSecondaryWeightByProcess(false); fParticleChange->ProposeParentWeight(new_weight); //Kinematic //--------- G4double projectileM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass(); G4double projectileTotalEnergy = projectileM0+projectileKinEnergy; G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0; G4double projectileP = std::sqrt(projectileP2); //Angle of the gamma direction with the projectile taken from G4eBremsstrahlungModel //------------------------------------------------ G4double u; const G4double a1 = 0.625 , a2 = 3.*a1 , d = 27. ; if (9./(9.+d) > G4UniformRand()) u = - std::log(G4UniformRand()*G4UniformRand())/a1; else u = - std::log(G4UniformRand()*G4UniformRand())/a2; G4double theta = u*electron_mass_c2/projectileTotalEnergy; G4double sint = std::sin(theta); G4double cost = std::cos(theta); G4double phi = twopi * G4UniformRand() ; G4ThreeVector projectileMomentum; projectileMomentum=G4ThreeVector(std::cos(phi)*sint,std::sin(phi)*sint,cost)*projectileP; //gamma frame if (IsScatProjToProjCase) {//the adjoint primary is the scattered e- G4ThreeVector gammaMomentum = (projectileTotalEnergy-adjointPrimTotalEnergy)*G4ThreeVector(0.,0.,1.); G4ThreeVector dirProd=projectileMomentum-gammaMomentum; G4double cost1 = std::cos(dirProd.angle(projectileMomentum)); G4double sint1 = std::sqrt(1.-cost1*cost1); projectileMomentum=G4ThreeVector(std::cos(phi)*sint1,std::sin(phi)*sint1,cost1)*projectileP; } projectileMomentum.rotateUz(theAdjointPrimary->GetMomentumDirection()); if (!IsScatProjToProjCase ){ //kill the primary and add a secondary fParticleChange->ProposeTrackStatus(fStopAndKill); fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum)); } else { fParticleChange->ProposeEnergy(projectileKinEnergy); fParticleChange->ProposeMomentumDirection(projectileMomentum.unit()); } } //////////////////////////////////////////////////////////////////////////////// // G4AdjointBremsstrahlungModel::~G4AdjointBremsstrahlungModel() {;} //////////////////////////////////////////////////////////////////////////////// // G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecond(const G4Material* aMaterial, G4double kinEnergyProj, // kinetic energy of the primary particle before the interaction G4double kinEnergyProd // kinetic energy of the secondary particle ) {/*if (UsePenelopeModel && !isPenelopeModelInitialised) { theEmModelManagerForFwdModels->Initialise(G4Electron::Electron(),G4Gamma::Gamma(),1.,0); isPenelopeModelInitialised =true; } */ return DiffCrossSectionPerVolumePrimToSecondApproximated2(aMaterial, kinEnergyProj, kinEnergyProd); /*return G4VEmAdjointModel::DiffCrossSectionPerVolumePrimToSecond(aMaterial, kinEnergyProj, kinEnergyProd);*/ } //////////////////////////////////////////////////////////////////////////////// // G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecondApproximated1( const G4Material* aMaterial, G4double kinEnergyProj, // kinetic energy of the primary particle before the interaction G4double kinEnergyProd // kinetic energy of the secondary particle ) { G4double dCrossEprod=0.; G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd); G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd); //In this approximation we consider that the secondary gammas are sampled with 1/Egamma energy distribution //This is what is applied in the discrete standard model before the rejection test that make a cooerction //The application of the same rejection function is not possble here. //The differentiation of the CS over Ecut does not produce neither a good differential CS. That is due to the // fact that in the discrete model the differential CS and the integrated CS are both fitted but separatly and // therefore do not allow a correct numerical differentiation of the integrated CS to get the differential one. // In the future we plan to use the brem secondary spectra from the G4Penelope implementation if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){ G4double sigma=theDirectEMModel->CrossSectionPerVolume(aMaterial,theDirectPrimaryPartDef,kinEnergyProj,1.*keV); dCrossEprod=sigma/kinEnergyProd/std::log(kinEnergyProj/keV); } return dCrossEprod; } //////////////////////////////////////////////////////////////////////////////// // G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecondApproximated2( const G4Material* material, G4double kinEnergyProj, // kinetic energy of the primary particle before the interaction G4double kinEnergyProd // kinetic energy of the secondary particle ) { //In this approximation we derive the direct cross section over Tcut=gamma energy, en after apply the Migdla correction factor //used in the direct model G4double dCrossEprod=0.; const G4ElementVector* theElementVector = material->GetElementVector(); const double* theAtomNumDensityVector = material->GetAtomicNumDensityVector(); G4double dum=0.; G4double E1=kinEnergyProd,E2=kinEnergyProd*1.001; G4double dE=E2-E1; for (size_t i=0; iGetNumberOfElements(); i++) { G4double C1=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,(*theElementVector)[i]->GetZ(),dum ,E1); G4double C2=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,(*theElementVector)[i]->GetZ(),dum,E2); dCrossEprod += theAtomNumDensityVector[i] * (C1-C2)/dE; } //Now the Migdal correction G4double totalEnergy = kinEnergyProj+electron_mass_c2 ; G4double kp2 = MigdalConstant*totalEnergy*totalEnergy *(material->GetElectronDensity()); G4double MigdalFactor = 1./(1.+kp2/(kinEnergyProd*kinEnergyProd)); // its seems that the factor used in the CS compuation i the direct //model is different than the one used in the secondary sampling by a //factor (1.+kp2) To be checked! dCrossEprod*=MigdalFactor; return dCrossEprod; } //////////////////////////////////////////////////////////////////////////////// // G4double G4AdjointBremsstrahlungModel::AdjointCrossSection(const G4MaterialCutsCouple* aCouple, G4double primEnergy, G4bool IsScatProjToProjCase) {/* if (UsePenelopeModel && !isPenelopeModelInitialised) { theEmModelManagerForFwdModels->Initialise(G4Electron::Electron(),G4Gamma::Gamma(),1.,0); isPenelopeModelInitialised =true; } */ if (UseMatrix) return G4VEmAdjointModel::AdjointCrossSection(aCouple,primEnergy,IsScatProjToProjCase); DefineCurrentMaterial(aCouple); G4double Cross=0.; lastCZ=theDirectEMModel->CrossSectionPerVolume(aCouple->GetMaterial(),theDirectPrimaryPartDef,100.*MeV,100.*MeV/std::exp(1.));//this give the constant above if (!IsScatProjToProjCase ){ G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(primEnergy); G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(primEnergy); if (Emax_proj>Emin_proj && primEnergy > currentTcutForDirectSecond) Cross= lastCZ*std::log(Emax_proj/Emin_proj); } else { G4double Emax_proj = GetSecondAdjEnergyMaxForScatProjToProjCase(primEnergy); G4double Emin_proj = GetSecondAdjEnergyMinForScatProjToProjCase(primEnergy,currentTcutForDirectSecond); if (Emax_proj>Emin_proj) Cross= lastCZ*std::log((Emax_proj-primEnergy)*Emin_proj/Emax_proj/(Emin_proj-primEnergy)); } return Cross; }