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$ >> 27 // 26 // ------------------------------------------- 28 // ------------------------------------------------------------------- 27 // 29 // 28 // GEANT4 Class file 30 // GEANT4 Class file 29 // 31 // 30 // 32 // 31 // File name: G4LossTableBuilder 33 // File name: G4LossTableBuilder 32 // 34 // 33 // Author: Vladimir Ivanchenko 35 // Author: Vladimir Ivanchenko 34 // 36 // 35 // Creation date: 03.01.2002 37 // Creation date: 03.01.2002 36 // 38 // 37 // Modifications: 39 // Modifications: 38 // 40 // 39 // 23-01-03 V.Ivanchenko Cut per region 41 // 23-01-03 V.Ivanchenko Cut per region 40 // 21-07-04 V.Ivanchenko Fix problem of range 42 // 21-07-04 V.Ivanchenko Fix problem of range for dedx=0 41 // 08-11-04 Migration to new interface of Stor 43 // 08-11-04 Migration to new interface of Store/Retrieve tables (V.Ivanchenko) 42 // 07-12-04 Fix of BuildDEDX table (V.Ivanchen 44 // 07-12-04 Fix of BuildDEDX table (V.Ivanchenko) 43 // 27-03-06 Add bool options isIonisation (V.I 45 // 27-03-06 Add bool options isIonisation (V.Ivanchenko) 44 // 16-01-07 Fill new (not old) DEDX table (V.I 46 // 16-01-07 Fill new (not old) DEDX table (V.Ivanchenko) 45 // 12-02-07 Use G4LPhysicsFreeVector for the i 47 // 12-02-07 Use G4LPhysicsFreeVector for the inverse range table (V.Ivanchenko) 46 // 24-06-09 Removed hidden bin in G4PhysicsVec 48 // 24-06-09 Removed hidden bin in G4PhysicsVector (V.Ivanchenko) 47 // 49 // 48 // Class Description: 50 // Class Description: 49 // 51 // 50 // ------------------------------------------- 52 // ------------------------------------------------------------------- 51 // 53 // 52 //....oooOO0OOooo........oooOO0OOooo........oo 54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 53 //....oooOO0OOooo........oooOO0OOooo........oo 55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 54 56 55 #include "G4LossTableBuilder.hh" 57 #include "G4LossTableBuilder.hh" 56 #include "G4SystemOfUnits.hh" 58 #include "G4SystemOfUnits.hh" 57 #include "G4PhysicsTable.hh" 59 #include "G4PhysicsTable.hh" 58 #include "G4PhysicsLogVector.hh" 60 #include "G4PhysicsLogVector.hh" 59 #include "G4PhysicsTableHelper.hh" 61 #include "G4PhysicsTableHelper.hh" 60 #include "G4PhysicsFreeVector.hh" << 62 #include "G4LPhysicsFreeVector.hh" 61 #include "G4ProductionCutsTable.hh" 63 #include "G4ProductionCutsTable.hh" 62 #include "G4MaterialCutsCouple.hh" 64 #include "G4MaterialCutsCouple.hh" 63 #include "G4Material.hh" 65 #include "G4Material.hh" 64 #include "G4VEmModel.hh" 66 #include "G4VEmModel.hh" 65 #include "G4ParticleDefinition.hh" 67 #include "G4ParticleDefinition.hh" 66 #include "G4LossTableManager.hh" 68 #include "G4LossTableManager.hh" 67 #include "G4EmParameters.hh" << 68 << 69 G4bool G4LossTableBuilder::baseMatFlag = false << 70 std::vector<G4double>* G4LossTableBuilder::the << 71 std::vector<G4int>* G4LossTableBuilder::the << 72 std::vector<G4bool>* G4LossTableBuilder::the << 73 69 74 //....oooOO0OOooo........oooOO0OOooo........oo 70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 75 71 76 G4LossTableBuilder::G4LossTableBuilder(G4bool << 72 G4LossTableBuilder::G4LossTableBuilder() 77 : isInitializer(master) << 78 { 73 { 79 theParameters = G4EmParameters::Instance(); << 74 splineFlag = true; 80 if (nullptr == theFlag) { << 75 isInitialized = false; 81 theDensityFactor = new std::vector<G4doubl << 82 theDensityIdx = new std::vector<G4int>; << 83 theFlag = new std::vector<G4bool>; << 84 } << 85 } << 86 << 87 //....oooOO0OOooo........oooOO0OOooo........oo << 88 76 89 G4LossTableBuilder::~G4LossTableBuilder() << 77 theDensityFactor = new std::vector<G4double>; 90 { << 78 theDensityIdx = new std::vector<G4int>; 91 if (isInitializer) { << 79 theFlag = new std::vector<G4bool>; 92 delete theDensityFactor; << 93 delete theDensityIdx; << 94 delete theFlag; << 95 theDensityFactor = nullptr; << 96 theDensityIdx = nullptr; << 97 theFlag = nullptr; << 98 } << 99 } << 100 << 101 //....oooOO0OOooo........oooOO0OOooo........oo << 102 << 103 const std::vector<G4int>* G4LossTableBuilder:: << 104 { << 105 return theDensityIdx; << 106 } << 107 << 108 //....oooOO0OOooo........oooOO0OOooo........oo << 109 << 110 const std::vector<G4double>* G4LossTableBuilde << 111 { << 112 return theDensityFactor; << 113 } 80 } 114 81 115 //....oooOO0OOooo........oooOO0OOooo........oo 82 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 116 83 117 G4bool G4LossTableBuilder::GetFlag(std::size_t << 84 G4LossTableBuilder::~G4LossTableBuilder() 118 { << 119 return (idx < theFlag->size()) ? (*theFlag)[ << 120 } << 121 << 122 //....oooOO0OOooo........oooOO0OOooo........oo << 123 << 124 G4bool G4LossTableBuilder::GetBaseMaterialFlag << 125 { 85 { 126 return baseMatFlag; << 86 delete theDensityFactor; >> 87 delete theDensityIdx; >> 88 delete theFlag; 127 } 89 } 128 90 129 //....oooOO0OOooo........oooOO0OOooo........oo 91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 130 92 131 void 93 void 132 G4LossTableBuilder::BuildDEDXTable(G4PhysicsTa 94 G4LossTableBuilder::BuildDEDXTable(G4PhysicsTable* dedxTable, 133 const std:: << 95 const std::vector<G4PhysicsTable*>& list) 134 { 96 { 135 InitialiseBaseMaterials(dedxTable); << 97 size_t n_processes = list.size(); 136 std::size_t n_processes = list.size(); << 98 //G4cout << "Nproc= " << n_processes << " Ncoup= " << dedxTable->size() << G4endl; 137 if(1 >= n_processes) { return; } 99 if(1 >= n_processes) { return; } 138 100 139 std::size_t nCouples = dedxTable->size(); << 101 size_t nCouples = dedxTable->size(); 140 //G4cout << "Nproc= " << n_processes << " nC << 141 // << dedxTable->size() << G4endl; << 142 if(0 >= nCouples) { return; } 102 if(0 >= nCouples) { return; } 143 103 144 for (std::size_t i=0; i<nCouples; ++i) { << 104 for (size_t i=0; i<nCouples; ++i) { 145 auto pv0 = static_cast<G4PhysicsLogVector* << 105 // if ((*theFlag)[i]) { 146 //if (0 == i) G4cout << i << ". pv0=" << p << 106 G4PhysicsLogVector* pv0 = static_cast<G4PhysicsLogVector*>((*(list[0]))[i]); 147 if(pv0 == nullptr) { continue; } << 107 if(pv0) { 148 std::size_t npoints = pv0->GetVectorLength << 108 size_t npoints = pv0->GetVectorLength(); 149 auto pv = new G4PhysicsLogVector(*pv0); << 109 G4PhysicsLogVector* pv = new G4PhysicsLogVector(*pv0); 150 for (std::size_t j=0; j<npoints; ++j) { << 110 // pv = new G4PhysicsLogVector(elow, ehigh, npoints-1); 151 G4double dedx = 0.0; << 111 pv->SetSpline(splineFlag); 152 for (std::size_t k=0; k<n_processes; ++k << 112 for (size_t j=0; j<npoints; ++j) { 153 const G4PhysicsVector* pv1 = (*(list[k]))[i] << 113 G4double dedx = 0.0; 154 //if (0 == i) G4cout << " " << k << ". p << 114 for (size_t k=0; k<n_processes; ++k) { 155 dedx += (*pv1)[j]; << 115 G4PhysicsVector* pv1 = (*(list[k]))[i]; >> 116 dedx += (*pv1)[j]; >> 117 } >> 118 pv->PutValue(j, dedx); 156 } 119 } 157 pv->PutValue(j, dedx); << 120 if(splineFlag) { pv->FillSecondDerivatives(); } >> 121 G4PhysicsTableHelper::SetPhysicsVector(dedxTable, i, pv); 158 } 122 } 159 if(splineFlag) { pv->FillSecondDerivatives << 160 G4PhysicsTableHelper::SetPhysicsVector(ded << 161 } 123 } 162 //G4cout << "### G4LossTableBuilder::BuildDE << 163 //G4cout << *dedxTable << G4endl; << 164 } 124 } 165 125 166 //....oooOO0OOooo........oooOO0OOooo........oo 126 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 167 127 168 void G4LossTableBuilder::BuildRangeTable(const 128 void G4LossTableBuilder::BuildRangeTable(const G4PhysicsTable* dedxTable, 169 G4Phy << 129 G4PhysicsTable* rangeTable, >> 130 G4bool isIonisation) 170 // Build range table from the energy loss tabl 131 // Build range table from the energy loss table 171 { 132 { 172 //G4cout << "### G4LossTableBuilder::BuildRa << 133 size_t nCouples = dedxTable->size(); 173 //G4cout << *const_cast<G4PhysicsTable*>(ded << 174 const std::size_t nCouples = dedxTable->size << 175 if(0 >= nCouples) { return; } 134 if(0 >= nCouples) { return; } 176 135 177 const std::size_t n = 100; << 136 size_t n = 100; 178 const G4double del = 1.0/(G4double)n; << 137 G4double del = 1.0/(G4double)n; 179 138 180 for (std::size_t i=0; i<nCouples; ++i) { << 139 for (size_t i=0; i<nCouples; ++i) { 181 auto pv = static_cast<G4PhysicsLogVector*> << 140 if(isIonisation) { 182 if((pv == nullptr) || (isBaseMatActive && << 141 if( !(*theFlag)[i] ) { continue; } 183 std::size_t npoints = pv->GetVectorLength( << 142 } 184 std::size_t bin0 = 0; << 143 G4PhysicsLogVector* pv = static_cast<G4PhysicsLogVector*>((*dedxTable)[i]); >> 144 size_t npoints = pv->GetVectorLength(); >> 145 size_t bin0 = 0; 185 G4double elow = pv->Energy(0); 146 G4double elow = pv->Energy(0); 186 G4double ehigh = pv->Energy(npoints-1); 147 G4double ehigh = pv->Energy(npoints-1); 187 G4double dedx1 = (*pv)[0]; 148 G4double dedx1 = (*pv)[0]; 188 149 >> 150 //G4cout << "i= " << i << "npoints= " << npoints << " dedx1= " << dedx1 << G4endl; >> 151 189 // protection for specific cases dedx=0 152 // protection for specific cases dedx=0 190 if(dedx1 == 0.0) { 153 if(dedx1 == 0.0) { 191 for (std::size_t k=1; k<npoints; ++k) { << 154 for (size_t k=1; k<npoints; ++k) { 192 ++bin0; << 155 bin0++; 193 elow = pv->Energy(k); 156 elow = pv->Energy(k); 194 dedx1 = (*pv)[k]; 157 dedx1 = (*pv)[k]; 195 if(dedx1 > 0.0) { break; } 158 if(dedx1 > 0.0) { break; } 196 } 159 } 197 npoints -= bin0; 160 npoints -= bin0; 198 } 161 } >> 162 //G4cout<<"New Range vector" << G4endl; >> 163 //G4cout<<"nbins= "<<npoints-1<<" elow= "<<elow<<" ehigh= "<<ehigh >> 164 // <<" bin0= " << bin0 <<G4endl; 199 165 200 // initialisation of a new vector 166 // initialisation of a new vector 201 if(npoints < 3) { npoints = 3; } << 167 if(npoints < 2) { npoints = 2; } 202 168 203 delete (*rangeTable)[i]; 169 delete (*rangeTable)[i]; 204 G4PhysicsLogVector* v; 170 G4PhysicsLogVector* v; 205 if(0 == bin0) { v = new G4PhysicsLogVector 171 if(0 == bin0) { v = new G4PhysicsLogVector(*pv); } 206 else { v = new G4PhysicsLogVector(elow, eh << 172 else { v = new G4PhysicsLogVector(elow, ehigh, npoints-1); } >> 173 >> 174 // dedx is exact zero cannot build range table >> 175 if(2 == npoints) { >> 176 v->PutValue(0,1000.); >> 177 v->PutValue(1,2000.); >> 178 G4PhysicsTableHelper::SetPhysicsVector(rangeTable, i, v); >> 179 return; >> 180 } >> 181 v->SetSpline(splineFlag); 207 182 208 // assumed dedx proportional to beta 183 // assumed dedx proportional to beta 209 G4double energy1 = v->Energy(0); 184 G4double energy1 = v->Energy(0); 210 G4double range = 2.*energy1/dedx1; 185 G4double range = 2.*energy1/dedx1; 211 /* << 186 //G4cout << "range0= " << range << G4endl; 212 G4cout << "New Range vector Npoints=" << v << 213 << " coupleIdx=" << i << " spline=" << v- << 214 << " Elow=" << v->GetMinEnergy() <<" Ehig << 215 << " DEDX(Elow)=" << dedx1 << " R(Elow)=" << 216 */ << 217 v->PutValue(0,range); 187 v->PutValue(0,range); 218 188 219 for (std::size_t j=1; j<npoints; ++j) { << 189 for (size_t j=1; j<npoints; ++j) { 220 190 221 G4double energy2 = v->Energy(j); 191 G4double energy2 = v->Energy(j); 222 G4double de = (energy2 - energy1) * 192 G4double de = (energy2 - energy1) * del; 223 G4double energy = energy2 + de*0.5; 193 G4double energy = energy2 + de*0.5; 224 G4double sum = 0.0; 194 G4double sum = 0.0; 225 std::size_t idx = j - 1; << 195 //G4cout << "j= " << j << " e1= " << energy1 << " e2= " << energy2 226 for (std::size_t k=0; k<n; ++k) { << 196 // << " n= " << n << G4endl; >> 197 for (size_t k=0; k<n; ++k) { 227 energy -= de; 198 energy -= de; 228 dedx1 = pv->Value(energy, idx); << 199 dedx1 = pv->Value(energy); 229 if(dedx1 > 0.0) { sum += de/dedx1; } 200 if(dedx1 > 0.0) { sum += de/dedx1; } 230 } 201 } 231 range += sum; 202 range += sum; 232 /* << 233 if(energy < 10.) << 234 G4cout << "j= " << j << " e1= " << energy1 < << 235 << " n= " << n << " range=" << range< << 236 */ << 237 v->PutValue(j,range); 203 v->PutValue(j,range); 238 energy1 = energy2; 204 energy1 = energy2; 239 } 205 } 240 if(splineFlag) { v->FillSecondDerivatives( 206 if(splineFlag) { v->FillSecondDerivatives(); } 241 G4PhysicsTableHelper::SetPhysicsVector(ran 207 G4PhysicsTableHelper::SetPhysicsVector(rangeTable, i, v); 242 } 208 } 243 //G4cout << "### Range table" << G4endl; << 244 //G4cout << *rangeTable << G4endl; << 245 } 209 } 246 210 247 //....oooOO0OOooo........oooOO0OOooo........oo 211 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 248 212 249 void << 213 void G4LossTableBuilder::BuildInverseRangeTable(const G4PhysicsTable* rangeTable, 250 G4LossTableBuilder::BuildInverseRangeTable(con << 214 G4PhysicsTable* invRangeTable, 251 G4P << 215 G4bool isIonisation) 252 // Build inverse range table from the energy l 216 // Build inverse range table from the energy loss table 253 { 217 { 254 std::size_t nCouples = rangeTable->size(); << 218 size_t nCouples = rangeTable->size(); 255 if(0 >= nCouples) { return; } 219 if(0 >= nCouples) { return; } 256 220 257 for (std::size_t i=0; i<nCouples; ++i) { << 221 for (size_t i=0; i<nCouples; ++i) { >> 222 >> 223 if(isIonisation) { >> 224 if( !(*theFlag)[i] ) { continue; } >> 225 } 258 G4PhysicsVector* pv = (*rangeTable)[i]; 226 G4PhysicsVector* pv = (*rangeTable)[i]; 259 if((pv == nullptr) || (isBaseMatActive && << 227 size_t npoints = pv->GetVectorLength(); 260 std::size_t npoints = pv->GetVectorLength( << 228 G4double rlow = (*pv)[0]; >> 229 G4double rhigh = (*pv)[npoints-1]; 261 230 262 delete (*invRangeTable)[i]; 231 delete (*invRangeTable)[i]; 263 auto v = new G4PhysicsFreeVector(npoints, << 232 G4LPhysicsFreeVector* v = new G4LPhysicsFreeVector(npoints,rlow,rhigh); >> 233 v->SetSpline(splineFlag); 264 234 265 for (std::size_t j=0; j<npoints; ++j) { << 235 for (size_t j=0; j<npoints; ++j) { 266 G4double e = pv->Energy(j); 236 G4double e = pv->Energy(j); 267 G4double r = (*pv)[j]; 237 G4double r = (*pv)[j]; 268 v->PutValues(j,r,e); 238 v->PutValues(j,r,e); 269 } 239 } 270 if (splineFlag) { v->FillSecondDerivatives << 240 if(splineFlag) { v->FillSecondDerivatives(); } 271 v->EnableLogBinSearch(theParameters->Numbe << 272 241 273 G4PhysicsTableHelper::SetPhysicsVector(inv 242 G4PhysicsTableHelper::SetPhysicsVector(invRangeTable, i, v); 274 } 243 } 275 //G4cout << "### Inverse range table" << G4e << 276 //G4cout << *invRangeTable << G4endl; << 277 } 244 } 278 245 279 //....oooOO0OOooo........oooOO0OOooo........oo 246 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 280 247 281 void G4LossTableBuilder::InitialiseBaseMateria << 248 void >> 249 G4LossTableBuilder::InitialiseBaseMaterials(G4PhysicsTable* table) 282 { 250 { 283 if(!isInitializer) { return; } << 251 size_t nCouples = table->size(); >> 252 size_t nFlags = theFlag->size(); >> 253 >> 254 if(nCouples == nFlags && isInitialized) { return; } >> 255 >> 256 isInitialized = true; >> 257 >> 258 //G4cout << "%%%%%% G4LossTableBuilder::InitialiseBaseMaterials Ncouples= " >> 259 // << nCouples << " FlagSize= " << nFlags << G4endl; >> 260 >> 261 // variable density check 284 const G4ProductionCutsTable* theCoupleTable= 262 const G4ProductionCutsTable* theCoupleTable= 285 G4ProductionCutsTable::GetProductionCutsTa 263 G4ProductionCutsTable::GetProductionCutsTable(); 286 std::size_t nCouples = theCoupleTable->GetTa << 264 287 std::size_t nFlags = theFlag->size(); << 288 /* 265 /* 289 G4cout << "### InitialiseBaseMaterials: nCou << 266 for(size_t i=0; i<nFlags; ++i) { 290 << " nFlags=" << nFlags << " isInit:" << is << 267 G4cout << "CoupleIdx= " << i << " Flag= " << (*theFlag)[i] 291 << " baseMat:" << baseMatFlag << G4endl; << 268 << " tableFlag= " << table->GetFlag(i) << " " >> 269 << theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial()->GetName() >> 270 << G4endl; >> 271 } 292 */ 272 */ 293 // define base material flag << 273 294 if(isBaseMatActive && !baseMatFlag) { << 274 // expand vectors 295 for(G4int i=0; i<(G4int)nCouples; ++i) { << 275 if(nFlags < nCouples) { 296 if(nullptr != theCoupleTable->GetMateria << 276 for(size_t i=nFlags; i<nCouples; ++i) { theDensityFactor->push_back(1.0); } 297 baseMatFlag = true; << 277 for(size_t i=nFlags; i<nCouples; ++i) { theDensityIdx->push_back(-1); } 298 isInitialized = false; << 278 for(size_t i=nFlags; i<nCouples; ++i) { theFlag->push_back(true); } 299 break; << 279 } >> 280 for(size_t i=0; i<nCouples; ++i) { >> 281 >> 282 // base material is needed only for a couple which is not >> 283 // initialised and for which tables will be computed >> 284 (*theFlag)[i] = table->GetFlag(i); >> 285 if ((*theDensityIdx)[i] < 0) { >> 286 (*theDensityIdx)[i] = i; >> 287 const G4MaterialCutsCouple* couple = >> 288 theCoupleTable->GetMaterialCutsCouple(i); >> 289 const G4ProductionCuts* pcuts = couple->GetProductionCuts(); >> 290 const G4Material* mat = couple->GetMaterial(); >> 291 const G4Material* bmat = mat->GetBaseMaterial(); >> 292 >> 293 // base material exists - find it and check if it can be reused >> 294 if(bmat) { >> 295 for(size_t j=0; j<nCouples; ++j) { >> 296 >> 297 if(j == i) { continue; } >> 298 const G4MaterialCutsCouple* bcouple = >> 299 theCoupleTable->GetMaterialCutsCouple(j); >> 300 >> 301 if(bcouple->GetMaterial() == bmat && >> 302 bcouple->GetProductionCuts() == pcuts) { >> 303 >> 304 // based couple exist in the same region >> 305 (*theDensityIdx)[i] = j; >> 306 (*theDensityFactor)[i] = mat->GetDensity()/bmat->GetDensity(); >> 307 (*theFlag)[i] = false; >> 308 >> 309 // ensure that there will no double initialisation >> 310 (*theDensityIdx)[j] = j; >> 311 (*theDensityFactor)[j] = 1.0; >> 312 (*theFlag)[j] = true; >> 313 break; >> 314 } >> 315 } 300 } 316 } 301 } 317 } 302 } 318 } >> 319 /* >> 320 for(size_t i=0; i<nCouples; ++i) { >> 321 G4cout << "CoupleIdx= " << i << " Flag= " << (*theFlag)[i] >> 322 << " TableFlag= " << table->GetFlag(i) << " " >> 323 << theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial()->GetName() >> 324 << G4endl; >> 325 } >> 326 G4cout << "%%%%%% G4LossTableBuilder::InitialiseBaseMaterials end" >> 327 << G4endl; >> 328 */ >> 329 } 303 330 304 if(nFlags != nCouples) { isInitialized = fal << 331 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 305 if(isInitialized) { return; } << 332 >> 333 void G4LossTableBuilder::InitialiseCouples() >> 334 { >> 335 isInitialized = true; 306 336 307 // reserve memory << 337 // variable density initialisation for the cas without tables >> 338 const G4ProductionCutsTable* theCoupleTable= >> 339 G4ProductionCutsTable::GetProductionCutsTable(); >> 340 size_t nCouples = theCoupleTable->GetTableSize(); >> 341 //G4cout << "%%%%%% G4LossTableBuilder::InitialiseCouples() nCouples= " >> 342 // << nCouples << G4endl; >> 343 >> 344 theDensityFactor->resize(nCouples, 1.0); >> 345 theDensityIdx->resize(nCouples, -1); 308 theFlag->resize(nCouples, true); 346 theFlag->resize(nCouples, true); 309 theDensityFactor->resize(nCouples,1.0); << 310 theDensityIdx->resize(nCouples, 0); << 311 347 312 // define default flag and index of used mat << 348 for(size_t i=0; i<nCouples; ++i) { 313 for (G4int i=0; i<(G4int)nCouples; ++i) { << 314 (*theFlag)[i] = (nullptr == table) ? true << 315 (*theDensityIdx)[i] = i; << 316 } << 317 isInitialized = true; << 318 if (!baseMatFlag) { return; } << 319 349 320 // use base materials << 321 for (G4int i=0; i<(G4int)nCouples; ++i) { << 322 // base material is needed only for a coup 350 // base material is needed only for a couple which is not 323 // initialised and for which tables will b 351 // initialised and for which tables will be computed 324 auto couple = theCoupleTable->GetMaterialC << 352 if ((*theDensityIdx)[i] < 0) { 325 auto pcuts = couple->GetProductionCuts(); << 353 (*theDensityIdx)[i] = i; 326 auto mat = couple->GetMaterial(); << 354 const G4MaterialCutsCouple* couple = 327 auto bmat = mat->GetBaseMaterial(); << 355 theCoupleTable->GetMaterialCutsCouple(i); 328 << 356 const G4ProductionCuts* pcuts = couple->GetProductionCuts(); 329 // base material exists - find it and chec << 357 const G4Material* mat = couple->GetMaterial(); 330 if(nullptr != bmat) { << 358 const G4Material* bmat = mat->GetBaseMaterial(); 331 for(G4int j=0; j<(G4int)nCouples; ++j) { << 359 332 if(j == i) { continue; } << 360 // base material exists - find it and check if it can be reused 333 auto bcouple = theCoupleTable->GetMaterialCu << 361 if(bmat) { 334 << 362 for(size_t j=0; j<nCouples; ++j) { 335 if(bcouple->GetMaterial() == bmat && << 363 336 bcouple->GetProductionCuts() == pcuts) { << 364 if(j == i) { continue; } 337 << 365 const G4MaterialCutsCouple* bcouple = 338 // based couple exist in the same region << 366 theCoupleTable->GetMaterialCutsCouple(j); 339 (*theDensityFactor)[i] = mat->GetDensity() << 367 340 (*theDensityIdx)[i] = j; << 368 if(bcouple->GetMaterial() == bmat && 341 (*theFlag)[i] = false; << 369 bcouple->GetProductionCuts() == pcuts) { 342 << 370 343 // ensure that there will no double initia << 371 // based couple exist in the same region 344 (*theDensityFactor)[j] = 1.0; << 372 (*theDensityIdx)[i] = j; 345 (*theDensityIdx)[j] = j; << 373 (*theDensityFactor)[i] = mat->GetDensity()/bmat->GetDensity(); 346 (*theFlag)[j] = true; << 374 (*theFlag)[i] = false; 347 break; << 375 >> 376 // ensure that there will no double initialisation >> 377 (*theDensityIdx)[j] = j; >> 378 (*theDensityFactor)[j] = 1.0; >> 379 (*theFlag)[j] = true; >> 380 break; >> 381 } 348 } 382 } 349 } 383 } 350 } 384 } 351 } 385 } >> 386 /* >> 387 for(size_t i=0; i<nCouples; ++i) { >> 388 G4cout << "CoupleIdx= " << i << " Flag= " << (*theFlag)[i] << " " >> 389 << theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial()->GetName() >> 390 << " DensityFactor= " << (*theDensityFactor)[i] >> 391 << G4endl; >> 392 } >> 393 G4cout << "%%%%%% G4LossTableBuilder::InitialiseCouples() end" << G4endl; >> 394 */ 352 } 395 } 353 396 354 //....oooOO0OOooo........oooOO0OOooo........oo 397 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 355 398 356 G4PhysicsTable* 399 G4PhysicsTable* 357 G4LossTableBuilder::BuildTableForModel(G4Physi 400 G4LossTableBuilder::BuildTableForModel(G4PhysicsTable* aTable, 358 G4VEmMo << 401 G4VEmModel* model, 359 const G << 402 const G4ParticleDefinition* part, 360 G4doubl << 403 G4double emin, G4double emax, 361 G4bool << 404 G4bool spline) 362 { 405 { 363 // check input 406 // check input 364 G4PhysicsTable* table = G4PhysicsTableHelper 407 G4PhysicsTable* table = G4PhysicsTableHelper::PreparePhysicsTable(aTable); 365 if (nullptr == table) { return table; } << 408 if(!table) { return table; } 366 if (aTable != nullptr && aTable != table) { << 409 if(emin >= emax) { 367 aTable->clearAndDestroy(); << 410 table->clearAndDestroy(); 368 delete aTable; << 411 delete table; >> 412 table = 0; >> 413 return table; 369 } 414 } 370 << 371 InitialiseBaseMaterials(table); 415 InitialiseBaseMaterials(table); 372 G4int nbins = theParameters->NumberOfBinsPer << 416 >> 417 G4int nbins = G4int(std::log10(emax/emin) + 0.5) >> 418 *G4LossTableManager::Instance()->GetNumberOfBinsPerDecade(); >> 419 if(nbins < 3) { nbins = 3; } 373 420 374 // Access to materials 421 // Access to materials 375 const G4ProductionCutsTable* theCoupleTable= 422 const G4ProductionCutsTable* theCoupleTable= 376 G4ProductionCutsTable::GetProductionCu 423 G4ProductionCutsTable::GetProductionCutsTable(); 377 std::size_t numOfCouples = theCoupleTable->G << 424 size_t numOfCouples = theCoupleTable->GetTableSize(); 378 /* << 425 379 G4cout << " G4LossTableBuilder::BuildTable << 426 G4PhysicsLogVector* aVector = 0; 380 << " isMaster=" << isInitializer << " model << 427 G4PhysicsLogVector* bVector = 0; 381 << " " << part->GetParticleName() << G4end << 428 382 */ << 429 for(size_t i=0; i<numOfCouples; ++i) { 383 G4PhysicsLogVector* aVector = nullptr; << 430 >> 431 //G4cout<< "i= " << i << " Flag= " << GetFlag(i) << G4endl; 384 432 385 for(G4int i=0; i<(G4int)numOfCouples; ++i) { << 433 if (GetFlag(i)) { 386 //G4cout << i << ". " << (*theFlag)[i] << << 387 if (table->GetFlag(i)) { << 388 434 389 // create physics vector and fill it 435 // create physics vector and fill it 390 auto couple = theCoupleTable->GetMateria << 436 const G4MaterialCutsCouple* couple = >> 437 theCoupleTable->GetMaterialCutsCouple(i); 391 delete (*table)[i]; 438 delete (*table)[i]; 392 439 393 // if start from zero then change the sc 440 // if start from zero then change the scale >> 441 394 const G4Material* mat = couple->GetMater 442 const G4Material* mat = couple->GetMaterial(); 395 443 396 G4double tmin = std::max(emin, model->Mi << 444 G4double tmin = std::max(emin,model->MinPrimaryEnergy(mat,part)); 397 if(0.0 >= tmin) { tmin = CLHEP::eV; } << 445 if(0.0 >= tmin) { tmin = eV; } 398 G4int n = nbins; << 446 G4int n = nbins + 1; 399 447 400 if(tmin >= emax) { 448 if(tmin >= emax) { 401 aVector = nullptr; << 449 aVector = 0; >> 450 } else if(tmin > emin) { >> 451 G4int bin = nbins*G4int(std::log10(emax/tmin) + 0.5); >> 452 if(bin < 3) { bin = 3; } >> 453 n = bin + 1; >> 454 aVector = new G4PhysicsLogVector(tmin, emax, bin); >> 455 >> 456 } else if(!bVector) { >> 457 aVector = new G4PhysicsLogVector(emin, emax, nbins); >> 458 bVector = aVector; >> 459 402 } else { 460 } else { 403 n *= G4lrint(std::log10(emax/tmin)); << 461 aVector = new G4PhysicsLogVector(*bVector); 404 n = std::max(n, 3); << 405 aVector = new G4PhysicsLogVector(tmin, << 406 } 462 } 407 463 408 if(nullptr != aVector) { << 464 if(aVector) { 409 //G4cout << part->GetParticleName() << << 465 aVector->SetSpline(spline); 410 // << " emin= " << tmin << " emax=" << e << 466 for(G4int j=0; j<n; ++j) { 411 for(G4int j=0; j<=n; ++j) { << 467 aVector->PutValue(j, model->Value(couple, part, aVector->Energy(j))); 412 G4double e = aVector->Energy(j); << 468 } 413 G4double y = model->Value(couple, part, e) << 469 if(spline) { aVector->FillSecondDerivatives(); } 414 //G4cout << " " << j << ") E=" << e < << 415 aVector->PutValue(j, y); << 416 } << 417 if(spline) { aVector->FillSecondDeriva << 418 } 470 } 419 G4PhysicsTableHelper::SetPhysicsVector(t 471 G4PhysicsTableHelper::SetPhysicsVector(table, i, aVector); 420 } 472 } 421 } 473 } 422 /* 474 /* 423 G4cout << "G4LossTableBuilder::BuildTableFor 475 G4cout << "G4LossTableBuilder::BuildTableForModel done for " 424 << part->GetParticleName() << " and " << 476 << part->GetParticleName() << " and "<< model->GetName() 425 << " " << table << G4endl; << 477 << " " << table << G4endl; 426 */ 478 */ 427 //G4cout << *table << G4endl; << 428 return table; 479 return table; 429 } 480 } 430 481 431 //....oooOO0OOooo........oooOO0OOooo........oo 482 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... >> 483 432 484