Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // >> 26 // $Id: G4eeToHadronsMultiModel.cc 66996 2013-01-29 14:50:52Z gcosmo $ 26 // 27 // 27 // ------------------------------------------- 28 // ------------------------------------------------------------------- 28 // 29 // 29 // GEANT4 Class file 30 // GEANT4 Class file 30 // 31 // 31 // 32 // 32 // File name: G4eeToHadronsMultiModel 33 // File name: G4eeToHadronsMultiModel 33 // 34 // 34 // Author: Vladimir Ivanchenko 35 // Author: Vladimir Ivanchenko 35 // 36 // 36 // Creation date: 02.08.2004 37 // Creation date: 02.08.2004 37 // 38 // 38 // Modifications: 39 // Modifications: 39 // 08-11-04 Migration to new interface of Stor 40 // 08-11-04 Migration to new interface of Store/Retrieve tables (V.Ivantchenko) 40 // 08-04-05 Major optimisation of internal int 41 // 08-04-05 Major optimisation of internal interfaces (V.Ivantchenko) 41 // 42 // 42 43 43 // 44 // 44 // ------------------------------------------- 45 // ------------------------------------------------------------------- 45 // 46 // 46 //....oooOO0OOooo........oooOO0OOooo........oo 47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 47 //....oooOO0OOooo........oooOO0OOooo........oo 48 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 48 49 49 #include "G4eeToHadronsMultiModel.hh" 50 #include "G4eeToHadronsMultiModel.hh" 50 #include "G4PhysicalConstants.hh" 51 #include "G4PhysicalConstants.hh" 51 #include "G4SystemOfUnits.hh" 52 #include "G4SystemOfUnits.hh" 52 #include "G4eeToTwoPiModel.hh" 53 #include "G4eeToTwoPiModel.hh" 53 #include "G4eeTo3PiModel.hh" 54 #include "G4eeTo3PiModel.hh" 54 #include "G4eeToPGammaModel.hh" 55 #include "G4eeToPGammaModel.hh" 55 #include "G4ee2KNeutralModel.hh" 56 #include "G4ee2KNeutralModel.hh" 56 #include "G4ee2KChargedModel.hh" 57 #include "G4ee2KChargedModel.hh" 57 #include "G4eeCrossSections.hh" 58 #include "G4eeCrossSections.hh" 58 #include "G4Vee2hadrons.hh" 59 #include "G4Vee2hadrons.hh" 59 #include "G4Positron.hh" << 60 60 61 //....oooOO0OOooo........oooOO0OOooo........oo 61 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 62 62 63 using namespace std; 63 using namespace std; 64 64 65 G4eeToHadronsMultiModel::G4eeToHadronsMultiMod << 65 G4eeToHadronsMultiModel::G4eeToHadronsMultiModel(G4int ver, const G4String& mname) 66 const G4String& mname) : G4VEmModel(mname), << 66 : G4VEmModel(mname), 67 { << 67 csFactor(1.0), 68 maxKineticEnergy = 4.521*CLHEP::GeV; //crre << 68 nModels(0), 69 delta = 1.0*CLHEP::MeV; //for bin width << 69 verbose(ver), >> 70 isInitialised(false) >> 71 { >> 72 thKineticEnergy = DBL_MAX; >> 73 maxKineticEnergy = 1.2*GeV; >> 74 fParticleChange = 0; >> 75 cross = 0; 70 } 76 } 71 77 72 //....oooOO0OOooo........oooOO0OOooo........oo 78 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 73 79 74 G4eeToHadronsMultiModel::~G4eeToHadronsMultiMo 80 G4eeToHadronsMultiModel::~G4eeToHadronsMultiModel() 75 { 81 { >> 82 G4int n = models.size(); >> 83 if(n>0) { >> 84 for(G4int i=0; i<n; i++) { >> 85 delete models[i]; >> 86 } >> 87 } 76 delete cross; 88 delete cross; 77 } 89 } 78 90 79 //....oooOO0OOooo........oooOO0OOooo........oo 91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 80 92 81 void G4eeToHadronsMultiModel::Initialise(const 93 void G4eeToHadronsMultiModel::Initialise(const G4ParticleDefinition*, 82 const G4DataVector& cuts) << 94 const G4DataVector&) 83 { 95 { 84 if(!isInitialised) { 96 if(!isInitialised) { 85 isInitialised = true; 97 isInitialised = true; 86 98 87 //G4cout<<"###Initialise in HadronMultiMod << 88 << 89 cross = new G4eeCrossSections(); 99 cross = new G4eeCrossSections(); 90 100 91 G4eeToTwoPiModel* m2pi = << 101 G4eeToTwoPiModel* m2pi = new G4eeToTwoPiModel(cross); 92 new G4eeToTwoPiModel(cross,maxKineticEne << 102 m2pi->SetHighEnergy(maxKineticEnergy); 93 AddEEModel(m2pi,cuts); << 103 AddEEModel(m2pi); 94 << 104 95 G4eeTo3PiModel* m3pi = << 105 G4eeTo3PiModel* m3pi1 = new G4eeTo3PiModel(cross); 96 new G4eeTo3PiModel(cross,maxKineticEnerg << 106 m3pi1->SetHighEnergy(0.95*GeV); 97 AddEEModel(m3pi,cuts); << 107 AddEEModel(m3pi1); 98 << 108 99 G4ee2KChargedModel* m2kc = << 109 G4eeTo3PiModel* m3pi2 = new G4eeTo3PiModel(cross); 100 new G4ee2KChargedModel(cross,maxKineticE << 110 m3pi2->SetLowEnergy(0.95*GeV); 101 AddEEModel(m2kc,cuts); << 111 m3pi2->SetHighEnergy(maxKineticEnergy); 102 << 112 AddEEModel(m3pi2); 103 G4ee2KNeutralModel* m2kn = << 113 104 new G4ee2KNeutralModel(cross,maxKineticE << 114 G4ee2KChargedModel* m2kc = new G4ee2KChargedModel(cross); 105 AddEEModel(m2kn,cuts); << 115 m2kc->SetHighEnergy(maxKineticEnergy); 106 << 116 AddEEModel(m2kc); 107 G4eeToPGammaModel* mpg1 = << 117 108 new G4eeToPGammaModel(cross,"pi0",maxKin << 118 G4ee2KNeutralModel* m2kn = new G4ee2KNeutralModel(cross); 109 AddEEModel(mpg1,cuts); << 119 m2kn->SetHighEnergy(maxKineticEnergy); 110 << 120 AddEEModel(m2kn); 111 G4eeToPGammaModel* mpg2 = << 121 112 new G4eeToPGammaModel(cross,"eta",maxKin << 122 G4eeToPGammaModel* mpg1 = new G4eeToPGammaModel(cross,"pi0"); 113 AddEEModel(mpg2,cuts); << 123 mpg1->SetLowEnergy(0.7*GeV); >> 124 mpg1->SetHighEnergy(maxKineticEnergy); >> 125 AddEEModel(mpg1); >> 126 >> 127 G4eeToPGammaModel* mpg2 = new G4eeToPGammaModel(cross,"eta"); >> 128 mpg2->SetLowEnergy(0.7*GeV); >> 129 mpg2->SetHighEnergy(maxKineticEnergy); >> 130 AddEEModel(mpg2); 114 131 115 nModels = (G4int)models.size(); << 132 nModels = models.size(); 116 133 117 fParticleChange = GetParticleChangeForGamm 134 fParticleChange = GetParticleChangeForGamma(); 118 } 135 } 119 } 136 } 120 137 121 //....oooOO0OOooo........oooOO0OOooo........oo 138 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 122 139 123 void G4eeToHadronsMultiModel::AddEEModel(G4Vee << 140 void G4eeToHadronsMultiModel::AddEEModel(G4Vee2hadrons* mod) 124 const G4DataVector& cuts) << 125 { 141 { 126 G4eeToHadronsModel* model = new G4eeToHadron 142 G4eeToHadronsModel* model = new G4eeToHadronsModel(mod, verbose); >> 143 model->SetLowEnergyLimit(LowEnergyLimit()); >> 144 model->SetHighEnergyLimit(HighEnergyLimit()); 127 models.push_back(model); 145 models.push_back(model); 128 G4double elow = mod->LowEnergy(); << 146 G4double elow = mod->ThresholdEnergy(); 129 ekinMin.push_back(elow); 147 ekinMin.push_back(elow); 130 if(thKineticEnergy > elow) { thKineticEnergy << 148 if(thKineticEnergy > elow) thKineticEnergy = elow; 131 ekinMax.push_back(mod->HighEnergy()); << 149 ekinMax.push_back(mod->HighEnergy()); 132 ekinPeak.push_back(mod->PeakEnergy()); 150 ekinPeak.push_back(mod->PeakEnergy()); 133 cumSum.push_back(0.0); 151 cumSum.push_back(0.0); 134 << 135 const G4ParticleDefinition* positron = G4Pos << 136 model->Initialise(positron,cuts); << 137 } 152 } 138 153 139 //....oooOO0OOooo........oooOO0OOooo........oo 154 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 140 155 141 G4double G4eeToHadronsMultiModel::CrossSection 156 G4double G4eeToHadronsMultiModel::CrossSectionPerVolume( 142 const G4Material* mat, 157 const G4Material* mat, 143 const G4ParticleDefinition* p, 158 const G4ParticleDefinition* p, 144 G4double kineticEnergy, 159 G4double kineticEnergy, 145 G4double, G4double) 160 G4double, G4double) 146 { 161 { 147 return mat->GetElectronDensity()* 162 return mat->GetElectronDensity()* 148 ComputeCrossSectionPerElectron(p, kineticE 163 ComputeCrossSectionPerElectron(p, kineticEnergy); 149 } 164 } 150 165 151 //....oooOO0OOooo........oooOO0OOooo........oo 166 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 152 167 153 G4double G4eeToHadronsMultiModel::ComputeCross 168 G4double G4eeToHadronsMultiModel::ComputeCrossSectionPerAtom( 154 const G4 169 const G4ParticleDefinition* p, 155 G4double kineticEnergy, 170 G4double kineticEnergy, 156 G4double Z, G4double, 171 G4double Z, G4double, 157 G4double, G4double) 172 G4double, G4double) 158 { 173 { 159 return Z*ComputeCrossSectionPerElectron(p, k 174 return Z*ComputeCrossSectionPerElectron(p, kineticEnergy); 160 } 175 } 161 176 162 //....oooOO0OOooo........oooOO0OOooo........oo << 163 << 164 G4double G4eeToHadronsMultiModel::ComputeCross << 165 G4double kineticEnergy, << 166 G4double, G4double) << 167 { << 168 G4double res = 0.0; << 169 << 170 G4double energy = LabToCM(kineticEnergy); << 171 << 172 if (energy > thKineticEnergy) { << 173 for(G4int i=0; i<nModels; i++) { << 174 if(energy >= ekinMin[i] && energy <= eki << 175 res += (models[i])->ComputeCrossSectio << 176 } << 177 cumSum[i] = res; << 178 } << 179 } << 180 return res*csFactor; << 181 } << 182 177 183 //....oooOO0OOooo........oooOO0OOooo........oo 178 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 184 179 185 void G4eeToHadronsMultiModel::SampleSecondarie << 180 void G4eeToHadronsMultiModel::SampleSecondaries(std::vector<G4DynamicParticle*>* newp, 186 std::vector<G4DynamicParticle*>* ne << 181 const G4MaterialCutsCouple* couple, 187 const G4MaterialCutsCouple* couple, << 182 const G4DynamicParticle* dp, 188 const G4DynamicParticle* dp, << 183 G4double, G4double) 189 G4double, G4double) << 190 { 184 { 191 G4double kinEnergy = dp->GetKineticEnergy(); 185 G4double kinEnergy = dp->GetKineticEnergy(); 192 G4double energy = LabToCM(kinEnergy); << 186 if (kinEnergy > thKineticEnergy) { 193 if (energy > thKineticEnergy) { << 194 G4double q = cumSum[nModels-1]*G4UniformRa 187 G4double q = cumSum[nModels-1]*G4UniformRand(); 195 for(G4int i=0; i<nModels; ++i) { << 188 for(G4int i=0; i<nModels; i++) { 196 if(q <= cumSum[i]) { 189 if(q <= cumSum[i]) { 197 (models[i])->SampleSecondaries(newp, c 190 (models[i])->SampleSecondaries(newp, couple,dp); 198 if(newp->size() > 0) { << 191 if(newp->size() > 0) fParticleChange->ProposeTrackStatus(fStopAndKill); 199 fParticleChange->ProposeTrackStatus(fStopA << 200 } << 201 break; 192 break; 202 } 193 } 203 } 194 } 204 } 195 } 205 } 196 } 206 197 207 //....oooOO0OOooo........oooOO0OOooo........oo 198 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 208 199 209 void G4eeToHadronsMultiModel::ModelDescription << 200 void G4eeToHadronsMultiModel::PrintInfo() 210 { 201 { 211 if(verbose > 0) { 202 if(verbose > 0) { 212 G4double e1 = 0.5*thKineticEnergy*thKineti 203 G4double e1 = 0.5*thKineticEnergy*thKineticEnergy/electron_mass_c2 213 - 2.0*electron_mass_c2; 204 - 2.0*electron_mass_c2; 214 G4double e2 = 0.5*maxKineticEnergy*maxKine 205 G4double e2 = 0.5*maxKineticEnergy*maxKineticEnergy/electron_mass_c2 215 - 2.0*electron_mass_c2; 206 - 2.0*electron_mass_c2; 216 outFile << " e+ annihilation into had << 207 G4cout << " e+ annihilation into hadrons active from " 217 << e1/GeV << " GeV to " << e2/GeV << " G << 208 << e1/GeV << " GeV to " << e2/GeV << " GeV" >> 209 << G4endl; 218 } 210 } 219 } 211 } 220 212 221 //....oooOO0OOooo........oooOO0OOooo........oo 213 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 222 214 223 void G4eeToHadronsMultiModel::SetCrossSecFacto 215 void G4eeToHadronsMultiModel::SetCrossSecFactor(G4double fac) 224 { 216 { 225 if(fac > 1.0) { 217 if(fac > 1.0) { 226 csFactor = fac; 218 csFactor = fac; 227 if(verbose > 0) { << 219 if(verbose > 0) 228 G4cout << "### G4eeToHadronsMultiModel: << 220 G4cout << "### G4eeToHadronsMultiModel: The cross section for G4eeToHadronsMultiModel " 229 << "G4eeToHadronsMultiModel is increase << 221 << " is increased by the Factor= " << csFactor << G4endl; 230 << csFactor << " times" << G4endl; << 231 } << 232 } 222 } 233 } 223 } 234 224 235 //....oooOO0OOooo........oooOO0OOooo........oo 225 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 236 226