Changeset 991 for trunk/source/processes/electromagnetic/utils/src
- Timestamp:
- Apr 17, 2009, 12:17:14 PM (15 years ago)
- Location:
- trunk/source/processes/electromagnetic/utils/src
- Files:
-
- 10 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/source/processes/electromagnetic/utils/src/G4EmCalculator.cc
r961 r991 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4EmCalculator.cc,v 1.4 6 2009/02/24 09:56:03vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-02 -ref-02$26 // $Id: G4EmCalculator.cc,v 1.44 2008/08/03 18:47:15 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-02 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 755 755 if(currentProcess) currentProcessName = currentProcess->GetProcessName(); 756 756 757 if(p->GetParticleType() == "nucleus" 758 && currentParticleName != "deuteron" 759 && currentParticleName != "triton" 760 && currentParticleName != "alpha+" 761 && currentParticleName != "helium" 762 && currentParticleName != "hydrogen" 763 ) { 757 if(p->GetParticleType() == "nucleus" && 758 currentParticleName != "deuteron" && currentParticleName != "triton") { 764 759 baseParticle = theGenericIon; 765 760 massRatio = baseParticle->GetPDGMass()/p->GetPDGMass(); -
trunk/source/processes/electromagnetic/utils/src/G4EmConfigurator.cc
r968 r991 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4EmConfigurator.cc,v 1. 4 2009/02/26 11:33:33vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-02 -ref-02$26 // $Id: G4EmConfigurator.cc,v 1.3 2008/11/21 12:30:29 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-02 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 117 117 G4String fname = ""; 118 118 if(fm) fname = fm->GetName(); 119 G4String mname = ""; 120 if(mod) mname = mod->GetName(); 121 AddModelForRegion(particleName, processName, mname, regionName, 119 AddModelForRegion(particleName, processName, mod->GetName(), regionName, 122 120 emin, emax, fname); 123 121 } … … 205 203 206 204 for(G4int i=0; i<nm; i++) { 207 G4String mname = ""; 208 if(modelList[i]) mname = modelList[i]->GetName(); 209 G4String fname = ""; 210 if(flucModelList[i]) fname = flucModelList[i]->GetName(); 211 if(modelName == mname && flucModelName == fname && 205 if(modelName == modelList[i]->GetName() && 212 206 (particleList[i] == "" || particleList[i] == particleName) ) { 213 207 mod = modelList[i]; … … 220 214 221 215 if(!mod) { 222 223 // set fluctuation model for ionisation processes 224 if(fluc && ptype == eloss) { 225 G4VEnergyLossProcess* p = reinterpret_cast<G4VEnergyLossProcess*>(proc); 226 p->SetFluctModel(fluc); 227 228 } else { 229 G4cout << "### G4EmConfigurator WARNING: fails to find a model <" 230 << modelName << "> for process <" 231 << processName << "> and " << particleName 232 << G4endl; 233 if(flucModelName != "") { 234 G4cout << " fluctuation model <" 235 << flucModelName << G4endl; 236 } 237 } 216 G4cout << "### G4EmConfigurator WARNING: fails to find a model <" 217 << modelName << "> for process <" 218 << processName << "> and " << particleName 219 << G4endl; 220 if(flucModelName != "") 221 G4cout << " fluctuation model <" 222 << flucModelName << G4endl; 238 223 } else { 239 224 -
trunk/source/processes/electromagnetic/utils/src/G4EmProcessOptions.cc
r961 r991 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4EmProcessOptions.cc,v 1.2 6 2009/02/18 14:43:27 vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-02 -ref-02$26 // $Id: G4EmProcessOptions.cc,v 1.24 2008/04/17 10:33:27 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-02 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 59 59 #include "G4VMultipleScattering.hh" 60 60 #include "G4Region.hh" 61 #include "G4RegionStore.hh"62 61 63 62 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... … … 393 392 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 394 393 395 void G4EmProcessOptions::ActivateDeexcitation(const G4String& pname, 396 G4bool val, 397 const G4String& reg) 398 { 399 G4RegionStore* regionStore = G4RegionStore::GetInstance(); 400 const G4Region* r = 0; 401 if(reg == "" || reg == "World") { 402 r = regionStore->GetRegion("DefaultRegionForTheWorld", false); 403 } else { 404 r = regionStore->GetRegion(reg, false); 405 } 406 if(!r) { 407 G4cout << "G4EmProcessOptions::ActivateDeexcitation ERROR: G4Region <" 408 << reg << "> not found, the command ignored" << G4endl; 409 return; 410 } 411 412 const std::vector<G4VEnergyLossProcess*>& v = 413 theManager->GetEnergyLossProcessVector(); 414 std::vector<G4VEnergyLossProcess*>::const_iterator itr; 415 for(itr = v.begin(); itr != v.end(); itr++) { 416 G4VEnergyLossProcess* p = *itr; 417 if(p) { 418 if(pname == p->GetProcessName()) p->ActivateDeexcitation(val,r); 419 } 420 } 421 const std::vector<G4VEmProcess*>& w = 422 theManager->GetEmProcessVector(); 423 std::vector<G4VEmProcess*>::const_iterator itp; 424 for(itp = w.begin(); itp != w.end(); itp++) { 425 G4VEmProcess* q = *itp; 426 if(q) { 427 if(pname == q->GetProcessName()) q->ActivateDeexcitation(val,r); 428 } 394 void G4EmProcessOptions::ActivateDeexcitation(G4bool val, const G4Region* r) 395 { 396 const std::vector<G4VEnergyLossProcess*>& v = 397 theManager->GetEnergyLossProcessVector(); 398 std::vector<G4VEnergyLossProcess*>::const_iterator itr; 399 for(itr = v.begin(); itr != v.end(); itr++) { 400 G4VEnergyLossProcess* p = *itr; 401 if(p) p->ActivateDeexcitation(val,r); 402 } 403 const std::vector<G4VEmProcess*>& w = 404 theManager->GetEmProcessVector(); 405 std::vector<G4VEmProcess*>::const_iterator itp; 406 for(itp = w.begin(); itp != w.end(); itp++) { 407 G4VEmProcess* q = *itp; 408 if(q) q->ActivateDeexcitation(val,r); 429 409 } 430 410 } -
trunk/source/processes/electromagnetic/utils/src/G4EnergyLossMessenger.cc
r961 r991 25 25 // 26 26 // 27 // $Id: G4EnergyLossMessenger.cc,v 1.3 7 2009/02/18 14:43:27vnivanch Exp $28 // GEANT4 tag $Name: geant4-09-02 -ref-02$27 // $Id: G4EnergyLossMessenger.cc,v 1.35 2008/10/20 13:27:45 vnivanch Exp $ 28 // GEANT4 tag $Name: geant4-09-02 $ 29 29 // 30 30 // ------------------------------------------------------------------- … … 178 178 aplCmd->AvailableForStates(G4State_PreInit,G4State_Idle); 179 179 180 deexCmd = new G4UIcommand("/process/em/deexcitation",this);181 deexCmd->SetGuidance("Set deexcitation flag per process and G4Region.");182 deexCmd->SetGuidance(" procName : process name");183 deexCmd->SetGuidance(" flag : flag");184 deexCmd->SetGuidance(" regName : G4Region name");185 186 G4UIparameter* pName = new G4UIparameter("pName",'s',false);187 deexCmd->SetParameter(pName);188 189 G4UIparameter* flag = new G4UIparameter("flag",'s',false);190 deexCmd->SetParameter(flag);191 192 G4UIparameter* regName = new G4UIparameter("regName",'s',false);193 deexCmd->SetParameter(regName);194 195 180 dedxCmd = new G4UIcmdWithAnInteger("/process/eLoss/binsDEDX",this); 196 181 dedxCmd->SetGuidance("Set number of bins for DEDX tables"); … … 274 259 delete MinSubSecCmd; 275 260 delete StepFuncCmd; 276 delete deexCmd;277 261 delete eLossDirectory; 278 262 delete mscDirectory; … … 336 320 } 337 321 338 if (command == deexCmd) {339 G4String s1 (""), s2(""), s3("");340 G4bool b = false;341 std::istringstream is(newValue);342 is >> s1 >> s2 >> s3;343 if(s2 == "true") b = true;344 opt->ActivateDeexcitation(s1,b,s3);345 }346 347 322 if (command == mscCmd) { 348 323 if(newValue == "Minimal") -
trunk/source/processes/electromagnetic/utils/src/G4LossTableBuilder.cc
r961 r991 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4LossTableBuilder.cc,v 1.2 8 2009/02/18 16:24:47vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-02 -ref-02$26 // $Id: G4LossTableBuilder.cc,v 1.27 2008/07/22 15:55:15 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-02 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 132 132 G4double dedx1 = pv->GetValue(elow, b); 133 133 134 //G4cout << "nbins= " << nbins << " dedx1= " << dedx1 << G4endl;135 136 134 // protection for specific cases dedx=0 137 135 if(dedx1 == 0.0) { 138 for (size_t k=1; k<nbins -1; k++) {136 for (size_t k=1; k<nbins; k++) { 139 137 bin0++; 140 138 elow = pv->GetLowEdgeEnergy(k); … … 145 143 } 146 144 147 //G4cout << "nbins= " << nbins << " elow= " << elow << " ehigh= " << ehigh << G4endl;148 145 // initialisation of a new vector 149 146 G4PhysicsLogVector* v = new G4PhysicsLogVector(elow, ehigh, nbins); 150 // dedx is exect zero151 if(nbins <= 2) {152 v->PutValue(0,1000.);153 v->PutValue(1,2000.);154 G4PhysicsTableHelper::SetPhysicsVector(rangeTable, i, v);155 return;156 }157 147 v->SetSpline(splineFlag); 158 148 -
trunk/source/processes/electromagnetic/utils/src/G4VEmFluctuationModel.cc
r961 r991 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4VEmFluctuationModel.cc,v 1. 4 2009/02/19 11:25:50vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-02 -ref-02$26 // $Id: G4VEmFluctuationModel.cc,v 1.3 2008/07/15 16:56:39 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-02 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 65 65 } 66 66 67 void G4VEmFluctuationModel::InitialiseMe(const G4ParticleDefinition*)68 {}69 67 70 void G4VEmFluctuationModel::SetParticleAndCharge(const G4ParticleDefinition*,71 G4double)72 {}73 74 -
trunk/source/processes/electromagnetic/utils/src/G4VEmModel.cc
r961 r991 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4VEmModel.cc,v 1.2 5 2009/02/19 09:57:36 vnivanchExp $27 // GEANT4 tag $Name: geant4-09-02 -ref-02$26 // $Id: G4VEmModel.cc,v 1.20 2008/11/13 23:13:18 schaelic Exp $ 27 // GEANT4 tag $Name: geant4-09-02 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 41 41 // 25.10.2005 Set default highLimit=100.TeV (V.Ivanchenko) 42 42 // 06.02.2006 add method ComputeMeanFreePath() (mma) 43 // 16.02.2009 Move implementations of virtual methods to source (VI)44 43 // 45 44 // … … 61 60 fluc(0), name(nam), lowLimit(0.1*keV), highLimit(100.0*TeV), 62 61 polarAngleLimit(0.0),secondaryThreshold(DBL_MAX),theLPMflag(false), 63 pParticleChange(0),nuclearStopping(false), 64 currentCouple(0),currentElement(0), 65 nsec(5),flagDeexcitation(false) 62 pParticleChange(0),nuclearStopping(false),nsec(5) 66 63 { 67 64 xsec.resize(nsec); … … 85 82 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 86 83 84 G4double G4VEmModel::CrossSectionPerVolume(const G4Material* material, 85 const G4ParticleDefinition* p, 86 G4double ekin, 87 G4double emin, 88 G4double emax) 89 { 90 SetupForMaterial(p, material, ekin); 91 G4double cross = 0.0; 92 const G4ElementVector* theElementVector = material->GetElementVector(); 93 const G4double* theAtomNumDensityVector = material->GetVecNbOfAtomsPerVolume(); 94 G4int nelm = material->GetNumberOfElements(); 95 if(nelm > nsec) { 96 xsec.resize(nelm); 97 nsec = nelm; 98 } 99 for (G4int i=0; i<nelm; i++) { 100 cross += theAtomNumDensityVector[i]* 101 ComputeCrossSectionPerAtom(p,(*theElementVector)[i],ekin,emin,emax); 102 xsec[i] = cross; 103 } 104 return cross; 105 } 106 107 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 108 109 G4double G4VEmModel::ComputeMeanFreePath(const G4ParticleDefinition* p, 110 G4double ekin, 111 const G4Material* material, 112 G4double emin, 113 G4double emax) 114 { 115 G4double mfp = DBL_MAX; 116 G4double cross = CrossSectionPerVolume(material,p,ekin,emin,emax); 117 if (cross > DBL_MIN) mfp = 1./cross; 118 return mfp; 119 } 120 121 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 122 87 123 void G4VEmModel::InitialiseElementSelectors(const G4ParticleDefinition* p, 88 124 const G4DataVector& cuts) 89 125 { 90 // initialise before run91 flagDeexcitation = false;92 93 126 G4int nbins = G4int(std::log10(highLimit/lowLimit) + 0.5); 94 127 if(nbins < 3) nbins = 3; … … 129 162 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 130 163 131 G4double G4VEmModel::ComputeDEDXPerVolume(const G4Material*,132 const G4ParticleDefinition*,133 G4double,G4double)134 {135 return 0.0;136 }137 164 138 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......139 140 G4double G4VEmModel::CrossSectionPerVolume(const G4Material* material,141 const G4ParticleDefinition* p,142 G4double ekin,143 G4double emin,144 G4double emax)145 {146 SetupForMaterial(p, material, ekin);147 G4double cross = 0.0;148 const G4ElementVector* theElementVector = material->GetElementVector();149 const G4double* theAtomNumDensityVector = material->GetVecNbOfAtomsPerVolume();150 G4int nelm = material->GetNumberOfElements();151 if(nelm > nsec) {152 xsec.resize(nelm);153 nsec = nelm;154 }155 for (G4int i=0; i<nelm; i++) {156 cross += theAtomNumDensityVector[i]*157 ComputeCrossSectionPerAtom(p,(*theElementVector)[i],ekin,emin,emax);158 xsec[i] = cross;159 }160 return cross;161 }162 163 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......164 165 G4double G4VEmModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,166 G4double, G4double, G4double,167 G4double, G4double)168 {169 return 0.0;170 }171 172 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......173 174 G4double G4VEmModel::MinEnergyCut(const G4ParticleDefinition*,175 const G4MaterialCutsCouple*)176 {177 return 0.0;178 }179 180 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......181 182 G4double G4VEmModel::GetChargeSquareRatio(const G4ParticleDefinition* p,183 const G4Material*, G4double)184 {185 G4double q = p->GetPDGCharge()/CLHEP::eplus;186 return q*q;187 }188 189 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......190 191 G4double G4VEmModel::GetParticleCharge(const G4ParticleDefinition* p,192 const G4Material*, G4double)193 {194 return p->GetPDGCharge();195 }196 197 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......198 199 void G4VEmModel::CorrectionsAlongStep(const G4MaterialCutsCouple*,200 const G4DynamicParticle*,201 G4double&,G4double&,G4double)202 {}203 204 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......205 206 void G4VEmModel::SampleDeexcitationAlongStep(const G4Material*,207 const G4Track&,208 G4double& )209 {}210 211 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......212 213 G4double G4VEmModel::MaxSecondaryEnergy(const G4ParticleDefinition*,214 G4double kineticEnergy)215 {216 return kineticEnergy;217 }218 219 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......220 221 void G4VEmModel::SampleScattering(const G4DynamicParticle*, G4double)222 {}223 224 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......225 226 G4double G4VEmModel::ComputeTruePathLengthLimit(const G4Track&,227 G4PhysicsTable*,228 G4double)229 {230 return DBL_MAX;231 }232 233 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......234 235 G4double G4VEmModel::ComputeGeomPathLength(G4double truePathLength)236 {237 return truePathLength;238 }239 240 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......241 242 G4double G4VEmModel::ComputeTrueStepLength(G4double geomPathLength)243 {244 return geomPathLength;245 }246 247 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......248 249 void G4VEmModel::DefineForRegion(const G4Region*)250 {}251 252 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......253 254 void G4VEmModel::SetupForMaterial(const G4ParticleDefinition*,255 const G4Material*, G4double)256 {}257 258 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... -
trunk/source/processes/electromagnetic/utils/src/G4VEmProcess.cc
r961 r991 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4VEmProcess.cc,v 1.6 2 2009/02/19 09:57:36 vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-02 -ref-02$26 // $Id: G4VEmProcess.cc,v 1.60 2008/10/17 14:46:16 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-02 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 91 91 applyCuts(false), 92 92 startFromNull(true), 93 useDeexcitation(false), 94 nDERegions(0), 95 idxDERegions(0), 96 currentModel(0), 93 nRegions(0), 94 selectedModel(0), 97 95 particle(0), 98 96 currentCouple(0) … … 137 135 delete modelManager; 138 136 (G4LossTableManager::Instance())->DeRegister(this); 139 }140 141 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....142 143 void G4VEmProcess::Clear()144 {145 delete [] theEnergyOfCrossSectionMax;146 delete [] theCrossSectionMax;147 delete [] idxDERegions;148 theEnergyOfCrossSectionMax = 0;149 theCrossSectionMax = 0;150 idxDERegions = 0;151 currentCouple = 0;152 preStepLambda = 0.0;153 mfpKinEnergy = DBL_MAX;154 deRegions.clear();155 nDERegions = 0;156 137 } 157 138 … … 181 162 theLambdaTable = G4PhysicsTableHelper::PreparePhysicsTable(theLambdaTable); 182 163 } 183 // Sub Cutoff and Deexcitation 184 if (nDERegions>0) { 185 186 const G4ProductionCutsTable* theCoupleTable= 187 G4ProductionCutsTable::GetProductionCutsTable(); 188 size_t numOfCouples = theCoupleTable->GetTableSize(); 189 190 idxDERegions = new G4bool[numOfCouples]; 191 192 for (size_t j=0; j<numOfCouples; j++) { 193 194 const G4MaterialCutsCouple* couple = 195 theCoupleTable->GetMaterialCutsCouple(j); 196 const G4ProductionCuts* pcuts = couple->GetProductionCuts(); 197 G4bool reg = false; 198 for(G4int i=0; i<nDERegions; i++) { 199 if(deRegions[i]) { 200 if(pcuts == deRegions[i]->GetProductionCuts()) reg = true; 201 } 202 } 203 idxDERegions[j] = reg; 204 } 205 } 206 if (1 < verboseLevel && nDERegions>0) { 207 G4cout << " Deexcitation is activated for regions: " << G4endl; 208 for (G4int i=0; i<nDERegions; i++) { 209 const G4Region* r = deRegions[i]; 210 G4cout << " " << r->GetName() << G4endl; 211 } 212 } 164 } 165 166 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 167 168 void G4VEmProcess::Clear() 169 { 170 if(theEnergyOfCrossSectionMax) delete [] theEnergyOfCrossSectionMax; 171 if(theCrossSectionMax) delete [] theCrossSectionMax; 172 theEnergyOfCrossSectionMax = 0; 173 theCrossSectionMax = 0; 174 currentCouple = 0; 175 preStepLambda = 0.0; 176 mfpKinEnergy = DBL_MAX; 213 177 } 214 178 … … 273 237 G4cout << *theLambdaTable << G4endl; 274 238 } 275 }276 }277 278 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....279 280 void G4VEmProcess::PrintInfoDefinition()281 {282 if(verboseLevel > 0) {283 G4cout << G4endl << GetProcessName() << ": for "284 << particle->GetParticleName();285 if(integral) G4cout << ", integral: 1 ";286 if(applyCuts) G4cout << ", applyCuts: 1 ";287 G4cout << " SubType= " << GetProcessSubType() << G4endl;288 if(buildLambdaTable) {289 G4cout << " Lambda tables from "290 << G4BestUnit(minKinEnergy,"Energy")291 << " to "292 << G4BestUnit(maxKinEnergy,"Energy")293 << " in " << nLambdaBins << " bins, spline: "294 << (G4LossTableManager::Instance())->SplineFlag()295 << G4endl;296 }297 PrintInfo();298 modelManager->DumpModelList(verboseLevel);299 }300 301 if(verboseLevel > 2 && buildLambdaTable) {302 G4cout << " LambdaTable address= " << theLambdaTable << G4endl;303 if(theLambdaTable) G4cout << (*theLambdaTable) << G4endl;304 239 } 305 240 } … … 369 304 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 370 305 306 G4double G4VEmProcess::GetMeanFreePath(const G4Track& track, 307 G4double, 308 G4ForceCondition* condition) 309 { 310 *condition = NotForced; 311 return G4VEmProcess::MeanFreePath(track); 312 } 313 314 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 315 371 316 G4VParticleChange* G4VEmProcess::PostStepDoIt(const G4Track& track, 372 317 const G4Step&) … … 397 342 } 398 343 399 SelectModel(finalT); 400 if(useDeexcitation) { 401 currentModel->SetDeexcitationFlag(idxDERegions[currentMaterialIndex]); 402 } 344 G4VEmModel* currentModel = SelectModel(finalT); 345 403 346 /* 404 347 if(0 < verboseLevel) { … … 461 404 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 462 405 406 void G4VEmProcess::PrintInfoDefinition() 407 { 408 if(verboseLevel > 0) { 409 G4cout << G4endl << GetProcessName() << ": for " 410 << particle->GetParticleName(); 411 if(integral) G4cout << ", integral: 1 "; 412 if(applyCuts) G4cout << ", applyCuts: 1 "; 413 G4cout << " SubType= " << GetProcessSubType() << G4endl; 414 if(buildLambdaTable) { 415 G4cout << " Lambda tables from " 416 << G4BestUnit(minKinEnergy,"Energy") 417 << " to " 418 << G4BestUnit(maxKinEnergy,"Energy") 419 << " in " << nLambdaBins << " bins, spline: " 420 << (G4LossTableManager::Instance())->SplineFlag() 421 << G4endl; 422 } 423 PrintInfo(); 424 modelManager->DumpModelList(verboseLevel); 425 } 426 427 if(verboseLevel > 2 && buildLambdaTable) { 428 G4cout << " LambdaTable address= " << theLambdaTable << G4endl; 429 if(theLambdaTable) G4cout << (*theLambdaTable) << G4endl; 430 } 431 } 432 433 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 434 435 G4double G4VEmProcess::CrossSectionPerVolume(G4double kineticEnergy, 436 const G4MaterialCutsCouple* couple) 437 { 438 // Cross section per atom is calculated 439 DefineMaterial(couple); 440 G4double cross = 0.0; 441 G4bool b; 442 if(theLambdaTable) { 443 cross = (((*theLambdaTable)[currentMaterialIndex])-> 444 GetValue(kineticEnergy, b)); 445 } else { 446 G4VEmModel* model = SelectModel(kineticEnergy); 447 cross = 448 model->CrossSectionPerVolume(currentMaterial,particle,kineticEnergy); 449 } 450 451 return cross; 452 } 453 454 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 455 463 456 G4bool G4VEmProcess::StorePhysicsTable(const G4ParticleDefinition* part, 464 457 const G4String& directory, … … 528 521 529 522 return yes; 530 }531 532 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....533 534 void G4VEmProcess::ActivateDeexcitation(G4bool val, const G4Region* r)535 {536 G4RegionStore* regionStore = G4RegionStore::GetInstance();537 const G4Region* reg = r;538 if (!reg) {reg = regionStore->GetRegion("DefaultRegionForTheWorld", false);}539 540 // the region is in the list541 if (nDERegions) {542 for (G4int i=0; i<nDERegions; i++) {543 if (reg == deRegions[i]) {544 if(!val) deRegions[i] = 0;545 return;546 }547 }548 }549 550 // new region551 if(val) {552 useDeexcitation = true;553 deRegions.push_back(reg);554 nDERegions++;555 } else {556 useDeexcitation = false;557 }558 }559 560 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....561 562 G4double G4VEmProcess::CrossSectionPerVolume(G4double kineticEnergy,563 const G4MaterialCutsCouple* couple)564 {565 // Cross section per atom is calculated566 DefineMaterial(couple);567 G4double cross = 0.0;568 G4bool b;569 if(theLambdaTable) {570 cross = (((*theLambdaTable)[currentMaterialIndex])->571 GetValue(kineticEnergy, b));572 } else {573 SelectModel(kineticEnergy);574 cross = currentModel->CrossSectionPerVolume(currentMaterial,575 particle,kineticEnergy);576 }577 578 return cross;579 }580 581 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....582 583 G4double G4VEmProcess::GetMeanFreePath(const G4Track& track,584 G4double,585 G4ForceCondition* condition)586 {587 *condition = NotForced;588 return G4VEmProcess::MeanFreePath(track);589 523 } 590 524 -
trunk/source/processes/electromagnetic/utils/src/G4VEnergyLossProcess.cc
r961 r991 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4VEnergyLossProcess.cc,v 1.14 6 2009/02/19 11:25:50vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-02 -ref-02$26 // $Id: G4VEnergyLossProcess.cc,v 1.143 2008/10/17 14:46:16 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-02 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 150 150 secondaryParticle(0), 151 151 nSCoffRegions(0), 152 nDERegions(0),153 152 idxSCoffRegions(0), 154 idxDERegions(0),155 153 nProcesses(0), 156 154 theDEDXTable(0), … … 178 176 isIonisation(true), 179 177 useSubCutoff(false), 180 useDeexcitation(false),181 178 particle(0), 182 179 currentCouple(0), … … 240 237 << G4endl; 241 238 delete vstrag; 242 Clea n();239 Clear(); 243 240 244 241 if ( !baseParticle ) { … … 294 291 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 295 292 296 void G4VEnergyLossProcess::Clea n()293 void G4VEnergyLossProcess::Clear() 297 294 { 298 295 if(1 < verboseLevel) { … … 305 302 delete [] theCrossSectionMax; 306 303 delete [] idxSCoffRegions; 307 delete [] idxDERegions;308 304 309 305 theDEDXAtMaxEnergy = 0; … … 316 312 scProcesses.clear(); 317 313 nProcesses = 0; 318 }319 320 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....321 322 G4double G4VEnergyLossProcess::MinPrimaryEnergy(const G4ParticleDefinition*,323 const G4Material*,324 G4double cut)325 {326 return cut;327 314 } 328 315 … … 377 364 } 378 365 379 Clea n();366 Clear(); 380 367 381 368 // Base particle and set of models can be defined here … … 420 407 minSubRange, verboseLevel); 421 408 422 // Sub Cutoff and Deexcitation423 if (nSCoffRegions>0 || nDERegions>0) {409 // Sub Cutoff Regime 410 if (nSCoffRegions>0) { 424 411 theSubCuts = modelManager->SubCutoff(); 425 412 … … 427 414 G4ProductionCutsTable::GetProductionCutsTable(); 428 415 size_t numOfCouples = theCoupleTable->GetTableSize(); 429 430 if(nSCoffRegions>0) idxSCoffRegions = new G4bool[numOfCouples]; 431 if(nDERegions>0) idxDERegions = new G4bool[numOfCouples]; 416 idxSCoffRegions = new G4int[numOfCouples]; 432 417 433 418 for (size_t j=0; j<numOfCouples; j++) { … … 436 421 theCoupleTable->GetMaterialCutsCouple(j); 437 422 const G4ProductionCuts* pcuts = couple->GetProductionCuts(); 438 439 if(nSCoffRegions>0) { 440 G4bool reg = false; 441 for(G4int i=0; i<nSCoffRegions; i++) { 442 if( pcuts == scoffRegions[i]->GetProductionCuts()) reg = true; 443 } 444 idxSCoffRegions[j] = reg; 445 } 446 if(nDERegions>0) { 447 G4bool reg = false; 448 for(G4int i=0; i<nDERegions; i++) { 449 if( pcuts == deRegions[i]->GetProductionCuts()) reg = true; 450 } 451 idxDERegions[j] = reg; 452 } 423 G4int reg = 0; 424 for(G4int i=0; i<nSCoffRegions; i++) { 425 if( pcuts == scoffRegions[i]->GetProductionCuts()) reg = 1; 426 } 427 idxSCoffRegions[j] = reg; 453 428 } 454 429 } … … 470 445 } 471 446 } 472 if (nDERegions) {473 G4cout << " Deexcitation is ON for regions: " << G4endl;474 for (G4int i=0; i<nDERegions; i++) {475 const G4Region* r = deRegions[i];476 G4cout << " " << r->GetName() << G4endl;477 }478 }479 447 } 480 448 } … … 511 479 if(isIonisation) G4cout << " isIonisation flag = 1"; 512 480 G4cout << G4endl; 481 } 482 } 483 484 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 485 486 void G4VEnergyLossProcess::ActivateSubCutoff(G4bool val, const G4Region* r) 487 { 488 G4RegionStore* regionStore = G4RegionStore::GetInstance(); 489 if(val) { 490 useSubCutoff = true; 491 if (!r) {r = regionStore->GetRegion("DefaultRegionForTheWorld", false);} 492 if (nSCoffRegions) { 493 for (G4int i=0; i<nSCoffRegions; i++) { 494 if (r == scoffRegions[i]) return; 495 } 496 } 497 scoffRegions.push_back(r); 498 nSCoffRegions++; 499 } else { 500 useSubCutoff = false; 513 501 } 514 502 } … … 652 640 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 653 641 654 void G4VEnergyLossProcess::PrintInfoDefinition() 655 { 656 if(0 < verboseLevel) { 657 G4cout << G4endl << GetProcessName() << ": for " 658 << particle->GetParticleName() 659 << " SubType= " << GetProcessSubType() 660 << G4endl 661 << " dE/dx and range tables from " 662 << G4BestUnit(minKinEnergy,"Energy") 663 << " to " << G4BestUnit(maxKinEnergy,"Energy") 664 << " in " << nBins << " bins" << G4endl 665 << " Lambda tables from threshold to " 666 << G4BestUnit(maxKinEnergy,"Energy") 667 << " in " << nBins << " bins, spline: " 668 << (G4LossTableManager::Instance())->SplineFlag() 669 << G4endl; 670 if(theRangeTableForLoss && isIonisation) { 671 G4cout << " finalRange(mm)= " << finalRange/mm 672 << ", dRoverRange= " << dRoverRange 673 << ", integral: " << integral 674 << ", fluct: " << lossFluctuationFlag 675 << ", linLossLimit= " << linLossLimit 676 << G4endl; 677 } 678 PrintInfo(); 679 modelManager->DumpModelList(verboseLevel); 680 if(theCSDARangeTable && isIonisation) { 681 G4cout << " CSDA range table up" 682 << " to " << G4BestUnit(maxKinEnergyCSDA,"Energy") 683 << " in " << nBinsCSDA << " bins" << G4endl; 684 } 685 if(nSCoffRegions>0 && isIonisation) { 686 G4cout << " Subcutoff sampling in " << nSCoffRegions 687 << " regions" << G4endl; 688 } 689 if(2 < verboseLevel) { 690 G4cout << " DEDXTable address= " << theDEDXTable << G4endl; 691 if(theDEDXTable && isIonisation) G4cout << (*theDEDXTable) << G4endl; 692 G4cout << "non restricted DEDXTable address= " 693 << theDEDXunRestrictedTable << G4endl; 694 if(theDEDXunRestrictedTable && isIonisation) { 695 G4cout << (*theDEDXunRestrictedTable) << G4endl; 696 } 697 if(theDEDXSubTable && isIonisation) { 698 G4cout << (*theDEDXSubTable) << G4endl; 699 } 700 G4cout << " CSDARangeTable address= " << theCSDARangeTable 701 << G4endl; 702 if(theCSDARangeTable && isIonisation) { 703 G4cout << (*theCSDARangeTable) << G4endl; 704 } 705 G4cout << " RangeTableForLoss address= " << theRangeTableForLoss 706 << G4endl; 707 if(theRangeTableForLoss && isIonisation) { 708 G4cout << (*theRangeTableForLoss) << G4endl; 709 } 710 G4cout << " InverseRangeTable address= " << theInverseRangeTable 711 << G4endl; 712 if(theInverseRangeTable && isIonisation) { 713 G4cout << (*theInverseRangeTable) << G4endl; 714 } 715 G4cout << " LambdaTable address= " << theLambdaTable << G4endl; 716 if(theLambdaTable && isIonisation) { 717 G4cout << (*theLambdaTable) << G4endl; 718 } 719 G4cout << " SubLambdaTable address= " << theSubLambdaTable << G4endl; 720 if(theSubLambdaTable && isIonisation) { 721 G4cout << (*theSubLambdaTable) << G4endl; 722 } 723 } 724 } 725 } 726 727 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 728 729 void G4VEnergyLossProcess::ActivateSubCutoff(G4bool val, const G4Region* r) 730 { 731 G4RegionStore* regionStore = G4RegionStore::GetInstance(); 732 const G4Region* reg = r; 733 if (!reg) {reg = regionStore->GetRegion("DefaultRegionForTheWorld", false);} 734 735 // the region is in the list 736 if (nSCoffRegions) { 737 for (G4int i=0; i<nSCoffRegions; i++) { 738 if (reg == scoffRegions[i]) { 739 if(!val) deRegions[i] = 0; 740 return; 741 } 742 } 743 } 744 745 // new region 746 if(val) { 747 useSubCutoff = true; 748 scoffRegions.push_back(reg); 749 nSCoffRegions++; 750 } else { 751 useSubCutoff = false; 752 } 753 } 754 755 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 756 757 void G4VEnergyLossProcess::ActivateDeexcitation(G4bool val, const G4Region* r) 758 { 759 G4RegionStore* regionStore = G4RegionStore::GetInstance(); 760 const G4Region* reg = r; 761 if (!reg) {reg = regionStore->GetRegion("DefaultRegionForTheWorld", false);} 762 763 // the region is in the list 764 if (nDERegions) { 765 for (G4int i=0; i<nDERegions; i++) { 766 if (reg == deRegions[i]) { 767 if(!val) deRegions[i] = 0; 768 return; 769 } 770 } 771 } 772 773 // new region 774 if(val) { 775 useDeexcitation = true; 776 deRegions.push_back(reg); 777 nDERegions++; 778 } else { 779 useDeexcitation = false; 780 } 642 G4double G4VEnergyLossProcess::GetContinuousStepLimit( 643 const G4Track&, 644 G4double, G4double, G4double&) 645 { 646 return DBL_MAX; 781 647 } 782 648 … … 808 674 // <<" stepLimit= "<<x<<G4endl; 809 675 return x; 676 } 677 678 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 679 680 G4double G4VEnergyLossProcess::GetMeanFreePath( 681 const G4Track& track, 682 G4double, 683 G4ForceCondition* condition) 684 685 { 686 *condition = NotForced; 687 return MeanFreePath(track); 810 688 } 811 689 … … 928 806 if (length >= fRange) { 929 807 eloss = preStepKinEnergy; 930 if (useDeexcitation) { 931 if(idxDERegions[currentMaterialIndex]) { 932 currentModel->SampleDeexcitationAlongStep(currentMaterial, track, eloss); 933 if(eloss < 0.0) eloss = 0.0; 934 } 935 } 808 currentModel->CorrectionsAlongStep(currentCouple, dynParticle, 809 eloss, esecdep, length); 936 810 fParticleChange.SetProposedKineticEnergy(0.0); 937 fParticleChange.ProposeLocalEnergyDeposit( eloss);811 fParticleChange.ProposeLocalEnergyDeposit(preStepKinEnergy); 938 812 return &fParticleChange; 939 813 } … … 1060 934 G4double tmax = 1061 935 std::min(currentModel->MaxSecondaryKinEnergy(dynParticle),cut); 1062 G4double emean = eloss;1063 936 eloss = fluc->SampleFluctuations(currentMaterial,dynParticle, 1064 tmax,length,e mean);937 tmax,length,eloss); 1065 938 /* 1066 939 if(-1 < verboseLevel) … … 1077 950 eloss += esecdep; 1078 951 if(eloss < 0.0) eloss = 0.0; 1079 1080 // deexcitation1081 else if (useDeexcitation) {1082 if(idxDERegions[currentMaterialIndex]) {1083 currentModel->SampleDeexcitationAlongStep(currentMaterial, track, eloss);1084 if(eloss < 0.0) eloss = 0.0;1085 }1086 }1087 952 1088 953 // Energy balanse … … 1241 1106 1242 1107 SelectModel(postStepScaledEnergy); 1243 if(useDeexcitation) {1244 currentModel->SetDeexcitationFlag(idxDERegions[currentMaterialIndex]);1245 }1246 1247 1108 const G4DynamicParticle* dynParticle = track.GetDynamicParticle(); 1248 1109 G4double tcut = (*theCuts)[currentMaterialIndex]; … … 1279 1140 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1280 1141 1281 G4bool G4VEnergyLossProcess::StorePhysicsTable( 1282 const G4ParticleDefinition* part, const G4String& directory, 1283 G4bool ascii) 1284 { 1285 G4bool res = true; 1286 if ( baseParticle || part != particle ) return res; 1287 1288 if(!StoreTable(part,theDEDXTable,ascii,directory,"DEDX")) 1289 {res = false;} 1290 1291 if(!StoreTable(part,theDEDXunRestrictedTable,ascii,directory,"DEDXnr")) 1292 {res = false;} 1293 1294 if(!StoreTable(part,theDEDXSubTable,ascii,directory,"SubDEDX")) 1295 {res = false;} 1296 1297 if(!StoreTable(part,theIonisationTable,ascii,directory,"Ionisation")) 1298 {res = false;} 1299 1300 if(!StoreTable(part,theIonisationSubTable,ascii,directory,"SubIonisation")) 1301 {res = false;} 1302 1303 if(isIonisation && 1304 !StoreTable(part,theCSDARangeTable,ascii,directory,"CSDARange")) 1305 {res = false;} 1306 1307 if(isIonisation && 1308 !StoreTable(part,theRangeTableForLoss,ascii,directory,"Range")) 1309 {res = false;} 1310 1311 if(isIonisation && 1312 !StoreTable(part,theInverseRangeTable,ascii,directory,"InverseRange")) 1313 {res = false;} 1314 1315 if(!StoreTable(part,theLambdaTable,ascii,directory,"Lambda")) 1316 {res = false;} 1317 1318 if(!StoreTable(part,theSubLambdaTable,ascii,directory,"SubLambda")) 1319 {res = false;} 1320 1321 if ( res ) { 1322 if(0 < verboseLevel) { 1323 G4cout << "Physics tables are stored for " << particle->GetParticleName() 1324 << " and process " << GetProcessName() 1325 << " in the directory <" << directory 1326 << "> " << G4endl; 1327 } 1328 } else { 1329 G4cout << "Fail to store Physics Tables for " 1330 << particle->GetParticleName() 1331 << " and process " << GetProcessName() 1332 << " in the directory <" << directory 1333 << "> " << G4endl; 1334 } 1335 return res; 1336 } 1337 1338 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... 1339 1340 G4bool G4VEnergyLossProcess::RetrievePhysicsTable( 1341 const G4ParticleDefinition* part, const G4String& directory, 1342 G4bool ascii) 1343 { 1344 G4bool res = true; 1345 const G4String particleName = part->GetParticleName(); 1346 1347 if(1 < verboseLevel) { 1348 G4cout << "G4VEnergyLossProcess::RetrievePhysicsTable() for " 1349 << particleName << " and process " << GetProcessName() 1350 << "; tables_are_built= " << tablesAreBuilt 1142 void G4VEnergyLossProcess::PrintInfoDefinition() 1143 { 1144 if(0 < verboseLevel) { 1145 G4cout << G4endl << GetProcessName() << ": for " 1146 << particle->GetParticleName() 1147 << " SubType= " << GetProcessSubType() 1148 << G4endl 1149 << " dE/dx and range tables from " 1150 << G4BestUnit(minKinEnergy,"Energy") 1151 << " to " << G4BestUnit(maxKinEnergy,"Energy") 1152 << " in " << nBins << " bins" << G4endl 1153 << " Lambda tables from threshold to " 1154 << G4BestUnit(maxKinEnergy,"Energy") 1155 << " in " << nBins << " bins, spline: " 1156 << (G4LossTableManager::Instance())->SplineFlag() 1351 1157 << G4endl; 1352 } 1353 if(particle == part) { 1354 1355 // G4bool yes = true; 1356 if ( !baseParticle ) { 1357 1358 G4bool fpi = true; 1359 if(!RetrieveTable(part,theDEDXTable,ascii,directory,"DEDX",fpi)) 1360 {fpi = false;} 1361 1362 if(!RetrieveTable(part,theIonisationTable,ascii,directory,"Ionisation",false)) 1363 {fpi = false;} 1364 1365 if(!RetrieveTable(part,theRangeTableForLoss,ascii,directory,"Range",fpi)) 1366 {res = false;} 1367 1368 if(!RetrieveTable(part,theDEDXunRestrictedTable,ascii,directory,"DEDXnr",false)) 1369 {res = false;} 1370 1371 if(!RetrieveTable(part,theCSDARangeTable,ascii,directory,"CSDARange",false)) 1372 {res = false;} 1373 1374 if(!RetrieveTable(part,theInverseRangeTable,ascii,directory,"InverseRange",fpi)) 1375 {res = false;} 1376 1377 if(!RetrieveTable(part,theLambdaTable,ascii,directory,"Lambda",true)) 1378 {res = false;} 1379 1380 G4bool yes = false; 1381 if(nSCoffRegions > 0) {yes = true;} 1382 1383 if(!RetrieveTable(part,theDEDXSubTable,ascii,directory,"SubDEDX",yes)) 1384 {res = false;} 1385 1386 if(!RetrieveTable(part,theSubLambdaTable,ascii,directory,"SubLambda",yes)) 1387 {res = false;} 1388 1389 if(!fpi) yes = false; 1390 if(!RetrieveTable(part,theIonisationSubTable,ascii,directory,"SubIonisation",yes)) 1391 {res = false;} 1392 } 1393 } 1394 1395 return res; 1396 } 1397 1398 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... 1399 1400 G4bool G4VEnergyLossProcess::StoreTable(const G4ParticleDefinition* part, 1401 G4PhysicsTable* aTable, G4bool ascii, 1402 const G4String& directory, 1403 const G4String& tname) 1404 { 1405 G4bool res = true; 1406 if ( aTable ) { 1407 const G4String name = GetPhysicsTableFileName(part,directory,tname,ascii); 1408 if( !aTable->StorePhysicsTable(name,ascii)) res = false; 1409 } 1410 return res; 1411 } 1412 1413 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... 1414 1415 G4bool G4VEnergyLossProcess::RetrieveTable(const G4ParticleDefinition* part, 1416 G4PhysicsTable* aTable, G4bool ascii, 1417 const G4String& directory, 1418 const G4String& tname, 1419 G4bool mandatory) 1420 { 1421 G4bool res = true; 1422 G4String filename = GetPhysicsTableFileName(part,directory,tname,ascii); 1423 G4bool yes = aTable->ExistPhysicsTable(filename); 1424 if(yes) { 1425 yes = G4PhysicsTableHelper::RetrievePhysicsTable(aTable,filename,ascii); 1426 if((G4LossTableManager::Instance())->SplineFlag()) { 1427 size_t n = aTable->length(); 1428 for(size_t i=0; i<n; i++) {(*aTable)[i]->SetSpline(true);} 1429 } 1430 } 1431 if(yes) { 1432 if (0 < verboseLevel) { 1433 G4cout << tname << " table for " << part->GetParticleName() 1434 << " is Retrieved from <" << filename << ">" 1158 if(theRangeTableForLoss && isIonisation) { 1159 G4cout << " finalRange(mm)= " << finalRange/mm 1160 << ", dRoverRange= " << dRoverRange 1161 << ", integral: " << integral 1162 << ", fluct: " << lossFluctuationFlag 1163 << ", linLossLimit= " << linLossLimit 1164 << G4endl; 1165 } 1166 PrintInfo(); 1167 modelManager->DumpModelList(verboseLevel); 1168 if(theCSDARangeTable && isIonisation) { 1169 G4cout << " CSDA range table up" 1170 << " to " << G4BestUnit(maxKinEnergyCSDA,"Energy") 1171 << " in " << nBinsCSDA << " bins" << G4endl; 1172 } 1173 if(nSCoffRegions>0 && isIonisation) { 1174 G4cout << " Subcutoff sampling in " << nSCoffRegions 1175 << " regions" << G4endl; 1176 } 1177 if(2 < verboseLevel) { 1178 G4cout << " DEDXTable address= " << theDEDXTable << G4endl; 1179 if(theDEDXTable && isIonisation) G4cout << (*theDEDXTable) << G4endl; 1180 G4cout << "non restricted DEDXTable address= " 1181 << theDEDXunRestrictedTable << G4endl; 1182 if(theDEDXunRestrictedTable && isIonisation) { 1183 G4cout << (*theDEDXunRestrictedTable) << G4endl; 1184 } 1185 if(theDEDXSubTable && isIonisation) { 1186 G4cout << (*theDEDXSubTable) << G4endl; 1187 } 1188 G4cout << " CSDARangeTable address= " << theCSDARangeTable 1435 1189 << G4endl; 1436 } 1437 } else { 1438 if(mandatory) res = false; 1439 if(mandatory || 1 < verboseLevel) { 1440 G4cout << tname << " table for " << part->GetParticleName() 1441 << " from file <" 1442 << filename << "> is not Retrieved" 1190 if(theCSDARangeTable && isIonisation) { 1191 G4cout << (*theCSDARangeTable) << G4endl; 1192 } 1193 G4cout << " RangeTableForLoss address= " << theRangeTableForLoss 1443 1194 << G4endl; 1444 } 1445 } 1446 return res; 1447 } 1448 1449 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1450 1451 G4double G4VEnergyLossProcess::GetDEDXDispersion( 1452 const G4MaterialCutsCouple *couple, 1453 const G4DynamicParticle* dp, 1454 G4double length) 1455 { 1456 DefineMaterial(couple); 1457 G4double ekin = dp->GetKineticEnergy(); 1458 SelectModel(ekin*massRatio); 1459 G4double tmax = currentModel->MaxSecondaryKinEnergy(dp); 1460 tmax = std::min(tmax,(*theCuts)[currentMaterialIndex]); 1461 G4double d = 0.0; 1462 G4VEmFluctuationModel* fm = currentModel->GetModelOfFluctuations(); 1463 if(fm) d = fm->Dispersion(currentMaterial,dp,tmax,length); 1464 return d; 1465 } 1466 1467 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1468 1469 G4double G4VEnergyLossProcess::CrossSectionPerVolume( 1470 G4double kineticEnergy, const G4MaterialCutsCouple* couple) 1471 { 1472 // Cross section per volume is calculated 1473 DefineMaterial(couple); 1474 G4double cross = 0.0; 1475 G4bool b; 1476 if(theLambdaTable) { 1477 cross = 1478 ((*theLambdaTable)[currentMaterialIndex])->GetValue(kineticEnergy, b); 1479 } else { 1480 SelectModel(kineticEnergy); 1481 cross = 1482 currentModel->CrossSectionPerVolume(currentMaterial, 1483 particle, kineticEnergy, 1484 (*theCuts)[currentMaterialIndex]); 1485 } 1486 1487 return cross; 1488 } 1489 1490 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1491 1492 G4double G4VEnergyLossProcess::MeanFreePath(const G4Track& track) 1493 { 1494 DefineMaterial(track.GetMaterialCutsCouple()); 1495 preStepLambda = GetLambdaForScaledEnergy(track.GetKineticEnergy()*massRatio); 1496 G4double x = DBL_MAX; 1497 if(DBL_MIN < preStepLambda) x = 1.0/preStepLambda; 1498 return x; 1499 } 1500 1501 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1502 1503 G4double G4VEnergyLossProcess::ContinuousStepLimit(const G4Track& track, 1504 G4double x, G4double y, 1505 G4double& z) 1506 { 1507 G4GPILSelection sel; 1508 return AlongStepGetPhysicalInteractionLength(track, x, y, z, &sel); 1509 } 1510 1511 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1512 1513 G4double G4VEnergyLossProcess::GetMeanFreePath( 1514 const G4Track& track, 1515 G4double, 1516 G4ForceCondition* condition) 1517 1518 { 1519 *condition = NotForced; 1520 return MeanFreePath(track); 1521 } 1522 1523 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1524 1525 G4double G4VEnergyLossProcess::GetContinuousStepLimit( 1526 const G4Track&, 1527 G4double, G4double, G4double&) 1528 { 1529 return DBL_MAX; 1530 } 1531 1532 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1533 1534 G4PhysicsVector* G4VEnergyLossProcess::LambdaPhysicsVector( 1535 const G4MaterialCutsCouple* couple, G4double cut) 1536 { 1537 // G4double cut = (*theCuts)[couple->GetIndex()]; 1538 // G4int nbins = nLambdaBins; 1539 G4double tmin = 1540 std::max(MinPrimaryEnergy(particle, couple->GetMaterial(), cut), 1541 minKinEnergy); 1542 if(tmin >= maxKinEnergy) tmin = 0.5*maxKinEnergy; 1543 G4PhysicsVector* v = new G4PhysicsLogVector(tmin, maxKinEnergy, nBins); 1544 v->SetSpline((G4LossTableManager::Instance())->SplineFlag()); 1545 1546 return v; 1547 } 1548 1549 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1550 1551 void G4VEnergyLossProcess::AddCollaborativeProcess( 1552 G4VEnergyLossProcess* p) 1553 { 1554 G4bool add = true; 1555 if(p->GetProcessName() != "eBrem") add = false; 1556 if(add && nProcesses > 0) { 1557 for(G4int i=0; i<nProcesses; i++) { 1558 if(p == scProcesses[i]) { 1559 add = false; 1560 break; 1561 } 1562 } 1563 } 1564 if(add) { 1565 scProcesses.push_back(p); 1566 nProcesses++; 1567 if (1 < verboseLevel) { 1568 G4cout << "### The process " << p->GetProcessName() 1569 << " is added to the list of collaborative processes of " 1570 << GetProcessName() << G4endl; 1195 if(theRangeTableForLoss && isIonisation) { 1196 G4cout << (*theRangeTableForLoss) << G4endl; 1197 } 1198 G4cout << " InverseRangeTable address= " << theInverseRangeTable 1199 << G4endl; 1200 if(theInverseRangeTable && isIonisation) { 1201 G4cout << (*theInverseRangeTable) << G4endl; 1202 } 1203 G4cout << " LambdaTable address= " << theLambdaTable << G4endl; 1204 if(theLambdaTable && isIonisation) { 1205 G4cout << (*theLambdaTable) << G4endl; 1206 } 1207 G4cout << " SubLambdaTable address= " << theSubLambdaTable << G4endl; 1208 if(theSubLambdaTable && isIonisation) { 1209 G4cout << (*theSubLambdaTable) << G4endl; 1210 } 1571 1211 } 1572 1212 } … … 1737 1377 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1738 1378 1379 G4PhysicsVector* G4VEnergyLossProcess::LambdaPhysicsVector( 1380 const G4MaterialCutsCouple* couple, G4double cut) 1381 { 1382 // G4double cut = (*theCuts)[couple->GetIndex()]; 1383 // G4int nbins = nLambdaBins; 1384 G4double tmin = 1385 std::max(MinPrimaryEnergy(particle, couple->GetMaterial(), cut), 1386 minKinEnergy); 1387 if(tmin >= maxKinEnergy) tmin = 0.5*maxKinEnergy; 1388 G4PhysicsVector* v = new G4PhysicsLogVector(tmin, maxKinEnergy, nBins); 1389 v->SetSpline((G4LossTableManager::Instance())->SplineFlag()); 1390 1391 return v; 1392 } 1393 1394 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1395 1396 G4double G4VEnergyLossProcess::CrossSectionPerVolume( 1397 G4double kineticEnergy, const G4MaterialCutsCouple* couple) 1398 { 1399 // Cross section per volume is calculated 1400 DefineMaterial(couple); 1401 G4double cross = 0.0; 1402 G4bool b; 1403 if(theLambdaTable) { 1404 cross = 1405 ((*theLambdaTable)[currentMaterialIndex])->GetValue(kineticEnergy, b); 1406 } else { 1407 SelectModel(kineticEnergy); 1408 cross = 1409 currentModel->CrossSectionPerVolume(currentMaterial, 1410 particle, kineticEnergy, 1411 (*theCuts)[currentMaterialIndex]); 1412 } 1413 1414 return cross; 1415 } 1416 1417 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1418 1419 G4bool G4VEnergyLossProcess::StorePhysicsTable( 1420 const G4ParticleDefinition* part, const G4String& directory, 1421 G4bool ascii) 1422 { 1423 G4bool res = true; 1424 if ( baseParticle || part != particle ) return res; 1425 1426 if(!StoreTable(part,theDEDXTable,ascii,directory,"DEDX")) 1427 {res = false;} 1428 1429 if(!StoreTable(part,theDEDXunRestrictedTable,ascii,directory,"DEDXnr")) 1430 {res = false;} 1431 1432 if(!StoreTable(part,theDEDXSubTable,ascii,directory,"SubDEDX")) 1433 {res = false;} 1434 1435 if(!StoreTable(part,theIonisationTable,ascii,directory,"Ionisation")) 1436 {res = false;} 1437 1438 if(!StoreTable(part,theIonisationSubTable,ascii,directory,"SubIonisation")) 1439 {res = false;} 1440 1441 if(isIonisation && 1442 !StoreTable(part,theCSDARangeTable,ascii,directory,"CSDARange")) 1443 {res = false;} 1444 1445 if(isIonisation && 1446 !StoreTable(part,theRangeTableForLoss,ascii,directory,"Range")) 1447 {res = false;} 1448 1449 if(isIonisation && 1450 !StoreTable(part,theInverseRangeTable,ascii,directory,"InverseRange")) 1451 {res = false;} 1452 1453 if(!StoreTable(part,theLambdaTable,ascii,directory,"Lambda")) 1454 {res = false;} 1455 1456 if(!StoreTable(part,theSubLambdaTable,ascii,directory,"SubLambda")) 1457 {res = false;} 1458 1459 if ( res ) { 1460 if(0 < verboseLevel) { 1461 G4cout << "Physics tables are stored for " << particle->GetParticleName() 1462 << " and process " << GetProcessName() 1463 << " in the directory <" << directory 1464 << "> " << G4endl; 1465 } 1466 } else { 1467 G4cout << "Fail to store Physics Tables for " 1468 << particle->GetParticleName() 1469 << " and process " << GetProcessName() 1470 << " in the directory <" << directory 1471 << "> " << G4endl; 1472 } 1473 return res; 1474 } 1475 1476 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... 1477 1478 G4bool G4VEnergyLossProcess::StoreTable(const G4ParticleDefinition* part, 1479 G4PhysicsTable* aTable, G4bool ascii, 1480 const G4String& directory, 1481 const G4String& tname) 1482 { 1483 G4bool res = true; 1484 if ( aTable ) { 1485 const G4String name = GetPhysicsTableFileName(part,directory,tname,ascii); 1486 if( !aTable->StorePhysicsTable(name,ascii)) res = false; 1487 } 1488 return res; 1489 } 1490 1491 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... 1492 1493 G4bool G4VEnergyLossProcess::RetrievePhysicsTable( 1494 const G4ParticleDefinition* part, const G4String& directory, 1495 G4bool ascii) 1496 { 1497 G4bool res = true; 1498 const G4String particleName = part->GetParticleName(); 1499 1500 if(1 < verboseLevel) { 1501 G4cout << "G4VEnergyLossProcess::RetrievePhysicsTable() for " 1502 << particleName << " and process " << GetProcessName() 1503 << "; tables_are_built= " << tablesAreBuilt 1504 << G4endl; 1505 } 1506 if(particle == part) { 1507 1508 // G4bool yes = true; 1509 if ( !baseParticle ) { 1510 1511 G4bool fpi = true; 1512 if(!RetrieveTable(part,theDEDXTable,ascii,directory,"DEDX",fpi)) 1513 {fpi = false;} 1514 1515 if(!RetrieveTable(part,theIonisationTable,ascii,directory,"Ionisation",false)) 1516 {fpi = false;} 1517 1518 if(!RetrieveTable(part,theRangeTableForLoss,ascii,directory,"Range",fpi)) 1519 {res = false;} 1520 1521 if(!RetrieveTable(part,theDEDXunRestrictedTable,ascii,directory,"DEDXnr",false)) 1522 {res = false;} 1523 1524 if(!RetrieveTable(part,theCSDARangeTable,ascii,directory,"CSDARange",false)) 1525 {res = false;} 1526 1527 if(!RetrieveTable(part,theInverseRangeTable,ascii,directory,"InverseRange",fpi)) 1528 {res = false;} 1529 1530 if(!RetrieveTable(part,theLambdaTable,ascii,directory,"Lambda",true)) 1531 {res = false;} 1532 1533 G4bool yes = false; 1534 if(nSCoffRegions > 0) {yes = true;} 1535 1536 if(!RetrieveTable(part,theDEDXSubTable,ascii,directory,"SubDEDX",yes)) 1537 {res = false;} 1538 1539 if(!RetrieveTable(part,theSubLambdaTable,ascii,directory,"SubLambda",yes)) 1540 {res = false;} 1541 1542 if(!fpi) yes = false; 1543 if(!RetrieveTable(part,theIonisationSubTable,ascii,directory,"SubIonisation",yes)) 1544 {res = false;} 1545 } 1546 } 1547 1548 return res; 1549 } 1550 1551 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... 1552 1553 G4bool G4VEnergyLossProcess::RetrieveTable(const G4ParticleDefinition* part, 1554 G4PhysicsTable* aTable, G4bool ascii, 1555 const G4String& directory, 1556 const G4String& tname, 1557 G4bool mandatory) 1558 { 1559 G4bool res = true; 1560 G4String filename = GetPhysicsTableFileName(part,directory,tname,ascii); 1561 G4bool yes = aTable->ExistPhysicsTable(filename); 1562 if(yes) { 1563 yes = G4PhysicsTableHelper::RetrievePhysicsTable(aTable,filename,ascii); 1564 if((G4LossTableManager::Instance())->SplineFlag()) { 1565 size_t n = aTable->length(); 1566 for(size_t i=0; i<n; i++) {(*aTable)[i]->SetSpline(true);} 1567 } 1568 } 1569 if(yes) { 1570 if (0 < verboseLevel) { 1571 G4cout << tname << " table for " << part->GetParticleName() 1572 << " is Retrieved from <" << filename << ">" 1573 << G4endl; 1574 } 1575 } else { 1576 if(mandatory) res = false; 1577 if(mandatory || 1 < verboseLevel) { 1578 G4cout << tname << " table for " << part->GetParticleName() 1579 << " from file <" 1580 << filename << "> is not Retrieved" 1581 << G4endl; 1582 } 1583 } 1584 return res; 1585 } 1586 1587 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1588 1589 void G4VEnergyLossProcess::AddCollaborativeProcess( 1590 G4VEnergyLossProcess* p) 1591 { 1592 G4bool add = true; 1593 if(p->GetProcessName() != "eBrem") add = false; 1594 if(add && nProcesses > 0) { 1595 for(G4int i=0; i<nProcesses; i++) { 1596 if(p == scProcesses[i]) { 1597 add = false; 1598 break; 1599 } 1600 } 1601 } 1602 if(add) { 1603 scProcesses.push_back(p); 1604 nProcesses++; 1605 if (1 < verboseLevel) { 1606 G4cout << "### The process " << p->GetProcessName() 1607 << " is added to the list of collaborative processes of " 1608 << GetProcessName() << G4endl; 1609 } 1610 } 1611 } 1612 1613 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1614 1615 G4double G4VEnergyLossProcess::GetDEDXDispersion( 1616 const G4MaterialCutsCouple *couple, 1617 const G4DynamicParticle* dp, 1618 G4double length) 1619 { 1620 DefineMaterial(couple); 1621 G4double ekin = dp->GetKineticEnergy(); 1622 SelectModel(ekin*massRatio); 1623 G4double tmax = currentModel->MaxSecondaryKinEnergy(dp); 1624 tmax = std::min(tmax,(*theCuts)[currentMaterialIndex]); 1625 G4double d = 0.0; 1626 G4VEmFluctuationModel* fm = currentModel->GetModelOfFluctuations(); 1627 if(fm) d = fm->Dispersion(currentMaterial,dp,tmax,length); 1628 return d; 1629 } 1630 1631 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1632 1633 void G4VEnergyLossProcess::ActivateDeexcitation(G4bool, const G4Region*) 1634 {} 1635 1636 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1637 -
trunk/source/processes/electromagnetic/utils/src/G4VMscModel.cc
r968 r991 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4VMscModel.cc,v 1. 10 2009/02/24 09:56:03vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-02 -ref-02$26 // $Id: G4VMscModel.cc,v 1.4 2008/03/10 18:39:45 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-02 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 49 49 50 50 #include "G4VMscModel.hh" 51 #include "G4ParticleChangeForMSC.hh"52 #include "G4TransportationManager.hh"53 51 54 52 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... … … 57 55 G4VMscModel::G4VMscModel(const G4String& nam): 58 56 G4VEmModel(nam), 59 safetyHelper(0),60 57 facrange(0.02), 61 58 facgeom(2.5), … … 64 61 dtrl(0.05), 65 62 lambdalimit(mm), 66 geommax(1.e50*mm),67 63 steppingAlgorithm(fUseSafety), 68 64 samplez(false), … … 76 72 77 73 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 78 79 void G4VMscModel::InitialiseSafetyHelper()80 {81 if(!safetyHelper) {82 safetyHelper = G4TransportationManager::GetTransportationManager()83 ->GetSafetyHelper();84 safetyHelper->InitialiseHelper();85 }86 }87 88 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......89 90 void G4VMscModel::ComputeDisplacement(G4ParticleChangeForMSC* fParticleChange,91 const G4ThreeVector& dir,92 G4double displacement,93 G4double postsafety)94 {95 const G4ThreeVector* pos = fParticleChange->GetProposedPosition();96 G4double r = displacement;97 if(r > postsafety) {98 G4double newsafety = safetyHelper->ComputeSafety(*pos);99 if(r > newsafety) r = newsafety;100 }101 if(r > 0.) {102 103 // compute new endpoint of the Step104 G4ThreeVector newPosition = *pos + r*dir;105 106 // definitely not on boundary107 if(displacement == r) {108 safetyHelper->ReLocateWithinVolume(newPosition);109 110 } else {111 // check safety after displacement112 G4double postsafety = safetyHelper->ComputeSafety(newPosition);113 114 // displacement to boundary115 if(postsafety <= 0.0) {116 safetyHelper->Locate(newPosition,117 *fParticleChange->GetProposedMomentumDirection());118 119 // not on the boundary120 } else {121 safetyHelper->ReLocateWithinVolume(newPosition);122 }123 }124 fParticleChange->ProposePosition(newPosition);125 }126 }127 128 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......129 130 void G4VMscModel::SampleSecondaries(std::vector<G4DynamicParticle*>*,131 const G4MaterialCutsCouple*,132 const G4DynamicParticle*,133 G4double, G4double)134 {}135 136 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Note: See TracChangeset
for help on using the changeset viewer.