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
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050 #include "G4eeToHadronsMultiModel.hh"
00051 #include "G4PhysicalConstants.hh"
00052 #include "G4SystemOfUnits.hh"
00053 #include "G4eeToTwoPiModel.hh"
00054 #include "G4eeTo3PiModel.hh"
00055 #include "G4eeToPGammaModel.hh"
00056 #include "G4ee2KNeutralModel.hh"
00057 #include "G4ee2KChargedModel.hh"
00058 #include "G4eeCrossSections.hh"
00059 #include "G4Vee2hadrons.hh"
00060
00061
00062
00063 using namespace std;
00064
00065 G4eeToHadronsMultiModel::G4eeToHadronsMultiModel(G4int ver, const G4String& mname)
00066 : G4VEmModel(mname),
00067 csFactor(1.0),
00068 nModels(0),
00069 verbose(ver),
00070 isInitialised(false)
00071 {
00072 thKineticEnergy = DBL_MAX;
00073 maxKineticEnergy = 1.2*GeV;
00074 fParticleChange = 0;
00075 cross = 0;
00076 }
00077
00078
00079
00080 G4eeToHadronsMultiModel::~G4eeToHadronsMultiModel()
00081 {
00082 G4int n = models.size();
00083 if(n>0) {
00084 for(G4int i=0; i<n; i++) {
00085 delete models[i];
00086 }
00087 }
00088 delete cross;
00089 }
00090
00091
00092
00093 void G4eeToHadronsMultiModel::Initialise(const G4ParticleDefinition*,
00094 const G4DataVector&)
00095 {
00096 if(!isInitialised) {
00097 isInitialised = true;
00098
00099 cross = new G4eeCrossSections();
00100
00101 G4eeToTwoPiModel* m2pi = new G4eeToTwoPiModel(cross);
00102 m2pi->SetHighEnergy(maxKineticEnergy);
00103 AddEEModel(m2pi);
00104
00105 G4eeTo3PiModel* m3pi1 = new G4eeTo3PiModel(cross);
00106 m3pi1->SetHighEnergy(0.95*GeV);
00107 AddEEModel(m3pi1);
00108
00109 G4eeTo3PiModel* m3pi2 = new G4eeTo3PiModel(cross);
00110 m3pi2->SetLowEnergy(0.95*GeV);
00111 m3pi2->SetHighEnergy(maxKineticEnergy);
00112 AddEEModel(m3pi2);
00113
00114 G4ee2KChargedModel* m2kc = new G4ee2KChargedModel(cross);
00115 m2kc->SetHighEnergy(maxKineticEnergy);
00116 AddEEModel(m2kc);
00117
00118 G4ee2KNeutralModel* m2kn = new G4ee2KNeutralModel(cross);
00119 m2kn->SetHighEnergy(maxKineticEnergy);
00120 AddEEModel(m2kn);
00121
00122 G4eeToPGammaModel* mpg1 = new G4eeToPGammaModel(cross,"pi0");
00123 mpg1->SetLowEnergy(0.7*GeV);
00124 mpg1->SetHighEnergy(maxKineticEnergy);
00125 AddEEModel(mpg1);
00126
00127 G4eeToPGammaModel* mpg2 = new G4eeToPGammaModel(cross,"eta");
00128 mpg2->SetLowEnergy(0.7*GeV);
00129 mpg2->SetHighEnergy(maxKineticEnergy);
00130 AddEEModel(mpg2);
00131
00132 nModels = models.size();
00133
00134 fParticleChange = GetParticleChangeForGamma();
00135 }
00136 }
00137
00138
00139
00140 void G4eeToHadronsMultiModel::AddEEModel(G4Vee2hadrons* mod)
00141 {
00142 G4eeToHadronsModel* model = new G4eeToHadronsModel(mod, verbose);
00143 model->SetLowEnergyLimit(LowEnergyLimit());
00144 model->SetHighEnergyLimit(HighEnergyLimit());
00145 models.push_back(model);
00146 G4double elow = mod->ThresholdEnergy();
00147 ekinMin.push_back(elow);
00148 if(thKineticEnergy > elow) thKineticEnergy = elow;
00149 ekinMax.push_back(mod->HighEnergy());
00150 ekinPeak.push_back(mod->PeakEnergy());
00151 cumSum.push_back(0.0);
00152 }
00153
00154
00155
00156 G4double G4eeToHadronsMultiModel::CrossSectionPerVolume(
00157 const G4Material* mat,
00158 const G4ParticleDefinition* p,
00159 G4double kineticEnergy,
00160 G4double, G4double)
00161 {
00162 return mat->GetElectronDensity()*
00163 ComputeCrossSectionPerElectron(p, kineticEnergy);
00164 }
00165
00166
00167
00168 G4double G4eeToHadronsMultiModel::ComputeCrossSectionPerAtom(
00169 const G4ParticleDefinition* p,
00170 G4double kineticEnergy,
00171 G4double Z, G4double,
00172 G4double, G4double)
00173 {
00174 return Z*ComputeCrossSectionPerElectron(p, kineticEnergy);
00175 }
00176
00177
00178
00179
00180 void G4eeToHadronsMultiModel::SampleSecondaries(std::vector<G4DynamicParticle*>* newp,
00181 const G4MaterialCutsCouple* couple,
00182 const G4DynamicParticle* dp,
00183 G4double, G4double)
00184 {
00185 G4double kinEnergy = dp->GetKineticEnergy();
00186 if (kinEnergy > thKineticEnergy) {
00187 G4double q = cumSum[nModels-1]*G4UniformRand();
00188 for(G4int i=0; i<nModels; i++) {
00189 if(q <= cumSum[i]) {
00190 (models[i])->SampleSecondaries(newp, couple,dp);
00191 if(newp->size() > 0) fParticleChange->ProposeTrackStatus(fStopAndKill);
00192 break;
00193 }
00194 }
00195 }
00196 }
00197
00198
00199
00200 void G4eeToHadronsMultiModel::PrintInfo()
00201 {
00202 if(verbose > 0) {
00203 G4double e1 = 0.5*thKineticEnergy*thKineticEnergy/electron_mass_c2
00204 - 2.0*electron_mass_c2;
00205 G4double e2 = 0.5*maxKineticEnergy*maxKineticEnergy/electron_mass_c2
00206 - 2.0*electron_mass_c2;
00207 G4cout << " e+ annihilation into hadrons active from "
00208 << e1/GeV << " GeV to " << e2/GeV << " GeV"
00209 << G4endl;
00210 }
00211 }
00212
00213
00214
00215 void G4eeToHadronsMultiModel::SetCrossSecFactor(G4double fac)
00216 {
00217 if(fac > 1.0) {
00218 csFactor = fac;
00219 if(verbose > 0)
00220 G4cout << "### G4eeToHadronsMultiModel: The cross section for G4eeToHadronsMultiModel "
00221 << " is increased by the Factor= " << csFactor << G4endl;
00222 }
00223 }
00224
00225