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 #include "G4VLEPTSModel.hh" 26 #include "G4VLEPTSModel.hh" 27 27 28 #include "CLHEP/Units/SystemOfUnits.h" 28 #include "CLHEP/Units/SystemOfUnits.h" 29 29 30 //....oooOO0OOooo........oooOO0OOooo........oo 30 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 31 G4VLEPTSModel::G4VLEPTSModel(const G4String& m << 31 G4VLEPTSModel::G4VLEPTSModel(const G4String& modelName) : G4VEmModel(modelName),isInitialised(false) 32 { 32 { 33 theMeanFreePathTable=nullptr; << 33 theMeanFreePathTable=NULL; 34 34 35 theNumbBinTable=100; 35 theNumbBinTable=100; 36 36 37 verboseLevel = 0; 37 verboseLevel = 0; 38 << 39 theLowestEnergyLimit = 0.5*CLHEP::eV; << 40 << 41 theHighestEnergyLimit = 1.0*CLHEP::MeV; << 42 << 43 theXSType = XSEnergy; << 44 } 38 } 45 39 46 40 47 //....oooOO0OOooo........oooOO0OOooo........oo 41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 48 G4VLEPTSModel::~G4VLEPTSModel() 42 G4VLEPTSModel::~G4VLEPTSModel() 49 { 43 { 50 44 51 if(theMeanFreePathTable != nullptr) { << 45 if(theMeanFreePathTable) { 52 theMeanFreePathTable->clearAndDestroy(); 46 theMeanFreePathTable->clearAndDestroy(); 53 delete theMeanFreePathTable; 47 delete theMeanFreePathTable; 54 } 48 } 55 } 49 } 56 50 57 51 58 //....oooOO0OOooo........oooOO0OOooo........oo 52 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 59 void G4VLEPTSModel::Init() 53 void G4VLEPTSModel::Init() 60 { 54 { 61 theLowestEnergyLimit = 0.5*CLHEP::eV; 55 theLowestEnergyLimit = 0.5*CLHEP::eV; 62 theHighestEnergyLimit = 1.0*CLHEP::MeV; 56 theHighestEnergyLimit = 1.0*CLHEP::MeV; 63 //t theHighestEnergyLimit = 15.0*CLHEP::M 57 //t theHighestEnergyLimit = 15.0*CLHEP::MeV; 64 SetLowEnergyLimit(theLowestEnergyLimit); 58 SetLowEnergyLimit(theLowestEnergyLimit); 65 SetHighEnergyLimit(theHighestEnergyLimit); 59 SetHighEnergyLimit(theHighestEnergyLimit); 66 theNumbBinTable = 100; 60 theNumbBinTable = 100; 67 61 68 } 62 } 69 63 70 64 71 65 72 //....oooOO0OOooo........oooOO0OOooo........oo 66 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 73 G4double G4VLEPTSModel::GetMeanFreePath(const 67 G4double G4VLEPTSModel::GetMeanFreePath(const G4Material* aMaterial, 74 const G4ParticleDefinition* , 68 const G4ParticleDefinition* , 75 G4double kineticEnergy ) 69 G4double kineticEnergy ) 76 { 70 { >> 71 77 G4double MeanFreePath; 72 G4double MeanFreePath; >> 73 G4bool isOutRange ; 78 74 79 if( verboseLevel >= 3 ) G4cout << aMaterial- 75 if( verboseLevel >= 3 ) G4cout << aMaterial->GetIndex() << " G4VLEPTSModel::GetMeanFreePath " << kineticEnergy << " > " << theHighestEnergyLimit << " < " << theLowestEnergyLimit << G4endl; 80 if (kineticEnergy > theHighestEnergyLimit || 76 if (kineticEnergy > theHighestEnergyLimit || kineticEnergy < theLowestEnergyLimit) 81 MeanFreePath = DBL_MAX; 77 MeanFreePath = DBL_MAX; 82 else 78 else 83 MeanFreePath = (*theMeanFreePathTable)(aMa << 79 MeanFreePath = (*theMeanFreePathTable)(aMaterial->GetIndex())-> >> 80 GetValue(kineticEnergy, isOutRange); 84 81 85 return MeanFreePath; 82 return MeanFreePath; 86 } 83 } 87 84 88 85 89 //....oooOO0OOooo........oooOO0OOooo........oo 86 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 90 void G4VLEPTSModel::BuildPhysicsTable(const G4 87 void G4VLEPTSModel::BuildPhysicsTable(const G4ParticleDefinition& aParticleType) 91 { 88 { 92 //CHECK IF PATH VARIABLE IS DEFINED 89 //CHECK IF PATH VARIABLE IS DEFINED 93 const char* path = G4FindDataDir("G4LEDATA") << 90 char* path = getenv("G4LEDATA"); 94 if( path == nullptr ) { << 91 if( !path ) { 95 G4Exception("G4VLEPTSModel", 92 G4Exception("G4VLEPTSModel", 96 "", 93 "", 97 FatalException, 94 FatalException, 98 "variable G4LEDATA not defined"); 95 "variable G4LEDATA not defined"); 99 } 96 } 100 97 101 // Build microscopic cross section table and 98 // Build microscopic cross section table and mean free path table 102 99 103 G4String aParticleName = aParticleType.GetPa 100 G4String aParticleName = aParticleType.GetParticleName(); 104 101 105 if (theMeanFreePathTable != nullptr) { << 102 if (theMeanFreePathTable) { 106 theMeanFreePathTable->clearAndDestroy(); 103 theMeanFreePathTable->clearAndDestroy(); 107 delete theMeanFreePathTable; 104 delete theMeanFreePathTable; 108 } 105 } 109 106 110 theMeanFreePathTable = new G4PhysicsTable(G4 107 theMeanFreePathTable = new G4PhysicsTable(G4Material::GetNumberOfMaterials()); 111 108 112 //LOOP TO MATERIALS IN GEOMETRY 109 //LOOP TO MATERIALS IN GEOMETRY 113 const G4MaterialTable * materialTable = G4Ma 110 const G4MaterialTable * materialTable = G4Material::GetMaterialTable() ; 114 std::vector<G4Material*>::const_iterator mat 111 std::vector<G4Material*>::const_iterator matite; 115 for( matite = materialTable->cbegin(); matit << 112 for( matite = materialTable->begin(); matite != materialTable->end(); matite++ ) { 116 const G4Material * aMaterial = (*matite); 113 const G4Material * aMaterial = (*matite); 117 G4String mateName = aMaterial->GetName(); 114 G4String mateName = aMaterial->GetName(); 118 115 119 //READ PARAMETERS FOR THIS MATERIAL 116 //READ PARAMETERS FOR THIS MATERIAL 120 std::string dirName = std::string(path) + << 117 string dirName = string(path) + "/lepts/"; 121 std::string fnParam = dirName + mateName << 118 string fnParam = dirName + mateName + "." + aParticleName + ".param.dat"; 122 std::string baseName = std::string(path) + << 119 string baseName = string(path) + "/lepts/" + mateName + "." + aParticleName; 123 G4bool bData = ReadParam( fnParam, aMateri 120 G4bool bData = ReadParam( fnParam, aMaterial ); 124 if( !bData ) continue; // MATERIAL NOT EX 121 if( !bData ) continue; // MATERIAL NOT EXISTING, DO NOT READ OTHER FILES 125 122 126 //READ INTEGRAL CROSS SECTION FOR THIS MAT 123 //READ INTEGRAL CROSS SECTION FOR THIS MATERIAL 127 std::string fnIXS = baseName + ".IXS.dat" << 124 string fnIXS = baseName + ".IXS.dat"; 128 125 129 std::map< G4int, std::vector<G4double> > i 126 std::map< G4int, std::vector<G4double> > integralXS = ReadIXS(fnIXS, aMaterial); 130 if( verboseLevel >= 2 ) G4cout << GetName( 127 if( verboseLevel >= 2 ) G4cout << GetName() << " : " << theXSType << " " << mateName << " INTEGRALXS " << integralXS.size() << G4endl; 131 128 132 if( integralXS.empty() ) { << 129 if( integralXS.size() == 0 ) { 133 G4cerr << " Integral cross sections will 130 G4cerr << " Integral cross sections will be set to 0. for material " << mateName << G4endl; 134 auto ptrVector = new G4PhysicsLogVector << 131 G4PhysicsLogVector* ptrVector = new G4PhysicsLogVector(theLowestEnergyLimit, theHighestEnergyLimit, 2); 135 ptrVector->PutValue(0, DBL_MAX); 132 ptrVector->PutValue(0, DBL_MAX); 136 ptrVector->PutValue(1, DBL_MAX); 133 ptrVector->PutValue(1, DBL_MAX); 137 134 138 std::size_t matIdx = aMaterial->GetIndex << 135 unsigned int matIdx = aMaterial->GetIndex(); 139 theMeanFreePathTable->insertAt( matIdx , 136 theMeanFreePathTable->insertAt( matIdx , ptrVector ) ; 140 137 141 } else { 138 } else { 142 139 143 if( verboseLevel >= 2 ) { 140 if( verboseLevel >= 2 ) { 144 std::map< G4int, std::vector<G4double> >::co 141 std::map< G4int, std::vector<G4double> >::const_iterator itei; 145 for( itei = integralXS.begin(); itei != inte 142 for( itei = integralXS.begin(); itei != integralXS.end(); itei++ ){ 146 G4cout << GetName() << " : " << (*itei).fi 143 G4cout << GetName() << " : " << (*itei).first << " INTEGRALXS NDATA " << (*itei).second.size() << G4endl; 147 } 144 } 148 } 145 } 149 146 150 BuildMeanFreePathTable( aMaterial, integ 147 BuildMeanFreePathTable( aMaterial, integralXS ); 151 148 152 std::string fnDXS = baseName + ".DXS.dat << 149 string fnDXS = baseName + ".DXS.dat"; 153 std::string fnRMT = baseName + ".RMT.dat << 150 string fnRMT = baseName + ".RMT.dat"; 154 std::string fnEloss = baseName + ".Eloss << 151 string fnEloss = baseName + ".Eloss.dat"; 155 std::string fnEloss2 = baseName + ".Elos << 152 string fnEloss2 = baseName + ".Eloss2.dat"; 156 153 157 theDiffXS[aMaterial] = new G4LEPTSDiffXS 154 theDiffXS[aMaterial] = new G4LEPTSDiffXS(fnDXS); 158 if( !theDiffXS[aMaterial]->IsFileFound() 155 if( !theDiffXS[aMaterial]->IsFileFound() ) { 159 G4Exception("G4VLEPTSModel::BuildPhysicsTabl 156 G4Exception("G4VLEPTSModel::BuildPhysicsTable", 160 "", 157 "", 161 FatalException, 158 FatalException, 162 (G4String("File not found :" + fnDXS). 159 (G4String("File not found :" + fnDXS).c_str())); 163 } 160 } 164 161 165 theRMTDistr[aMaterial] = new G4LEPTSDist 162 theRMTDistr[aMaterial] = new G4LEPTSDistribution(); 166 theRMTDistr[aMaterial]->ReadFile(fnRMT); 163 theRMTDistr[aMaterial]->ReadFile(fnRMT); 167 164 168 theElostDistr[aMaterial] = new G4LEPTSEl 165 theElostDistr[aMaterial] = new G4LEPTSElossDistr(fnEloss); 169 if( !theElostDistr[aMaterial]->IsFileFou 166 if( !theElostDistr[aMaterial]->IsFileFound() ) { 170 G4Exception("G4VLEPTSModel::BuildPhysicsTabl 167 G4Exception("G4VLEPTSModel::BuildPhysicsTable", 171 "", 168 "", 172 FatalException, 169 FatalException, 173 (G4String("File not found :" + fnEloss 170 (G4String("File not found :" + fnEloss).c_str())); 174 } 171 } 175 } 172 } 176 173 177 } 174 } 178 175 179 } 176 } 180 177 181 void G4VLEPTSModel::BuildMeanFreePathTable( co 178 void G4VLEPTSModel::BuildMeanFreePathTable( const G4Material* aMaterial, std::map< G4int, std::vector<G4double> >& integralXS ) 182 { 179 { 183 G4double LowEdgeEnergy, fValue; 180 G4double LowEdgeEnergy, fValue; 184 181 185 //BUILD MEAN FREE PATH TABLE FROM INTEGRAL C 182 //BUILD MEAN FREE PATH TABLE FROM INTEGRAL CROSS SECTION 186 std::size_t matIdx = aMaterial->GetIndex(); << 183 unsigned int matIdx = aMaterial->GetIndex(); 187 auto ptrVector = new G4PhysicsLogVector(the << 184 G4PhysicsLogVector* ptrVector = new G4PhysicsLogVector(theLowestEnergyLimit, theHighestEnergyLimit, theNumbBinTable); 188 185 189 for (G4int ii=0; ii < theNumbBinTable; ++ii) << 186 for (G4int ii=0; ii < theNumbBinTable; ii++) { 190 LowEdgeEnergy = ptrVector->Energy(ii); << 187 LowEdgeEnergy = ptrVector->GetLowEdgeEnergy(ii); 191 if( verboseLevel >= 2 ) G4cout << GetName( << 188 if( verboseLevel >= 2 ) G4cout << GetName() << " " << ii << " LowEdgeEnergy " << LowEdgeEnergy << " > " << theLowestEnergyLimit << " < " << theHighestEnergyLimit << G4endl; 192 //- fValue = ComputeMFP(LowEdgeEnergy 189 //- fValue = ComputeMFP(LowEdgeEnergy, material, aParticleName); 193 fValue = 0.; 190 fValue = 0.; 194 if( LowEdgeEnergy >= theLowestEnergyLimit 191 if( LowEdgeEnergy >= theLowestEnergyLimit && 195 LowEdgeEnergy <= theHighestEnergyLimit) { 192 LowEdgeEnergy <= theHighestEnergyLimit) { 196 G4double NbOfMoleculesPerVolume = aMater 193 G4double NbOfMoleculesPerVolume = aMaterial->GetDensity()/theMolecularMass[aMaterial]*CLHEP::Avogadro; 197 194 198 G4double SIGMA = 0. ; 195 G4double SIGMA = 0. ; 199 //- for ( std::size_t elm=0 ; elm < << 196 //- for ( size_t elm=0 ; elm < aMaterial->GetNumberOfElements() ; elm++ ) { 200 G4double crossSection = 0.; 197 G4double crossSection = 0.; 201 198 202 G4double eVEnergy = LowEdgeEnergy/CLHEP::eV; 199 G4double eVEnergy = LowEdgeEnergy/CLHEP::eV; 203 200 204 //- if( verboseLevel >= 2 ) G4cout << " 201 //- if( verboseLevel >= 2 ) G4cout << " eVEnergy " << eVEnergy << " LowEdgeE " << LowEdgeEnergy << " " << integralXS[theXSType][1] << G4endl; 205 202 206 if(eVEnergy < integralXS[0][1] ) { 203 if(eVEnergy < integralXS[0][1] ) { 207 crossSection = 0.; 204 crossSection = 0.; 208 } else { 205 } else { 209 G4int Bin = 0; // locate bin 206 G4int Bin = 0; // locate bin 210 G4double aa, bb; 207 G4double aa, bb; 211 for( G4int jj=1; jj<theNXSdat[aMaterial]; 208 for( G4int jj=1; jj<theNXSdat[aMaterial]; jj++) { // Extrapolate for E > Emax !!! 212 if( verboseLevel >= 3 ) G4cout << " GET 209 if( verboseLevel >= 3 ) G4cout << " GET BIN " << jj << " "<< eVEnergy << " > " << integralXS[0][jj] << G4endl; 213 if( eVEnergy > integralXS[0][jj]) { 210 if( eVEnergy > integralXS[0][jj]) { 214 Bin = jj; 211 Bin = jj; 215 } else { 212 } else { 216 break; 213 break; 217 } 214 } 218 } 215 } 219 aa = integralXS[0][Bin]; 216 aa = integralXS[0][Bin]; 220 bb = integralXS[0][Bin+1]; 217 bb = integralXS[0][Bin+1]; 221 crossSection = (integralXS[theXSType][Bin] 218 crossSection = (integralXS[theXSType][Bin] + (integralXS[theXSType][Bin+1]-integralXS[theXSType][Bin])/(bb-aa)*(eVEnergy-aa) ) * 1.e-16*CLHEP::cm2; 222 219 223 if( verboseLevel >= 3 ) G4cout << " crossS 220 if( verboseLevel >= 3 ) G4cout << " crossSection " << crossSection << " " <<integralXS[theXSType][Bin] << " + " << (integralXS[theXSType][Bin+1]-integralXS[theXSType][Bin]) << " / " << (bb-aa) << " *" << (eVEnergy-aa) << " * " << 1.e-16*CLHEP::cm2 << G4endl;; 224 221 225 // SIGMA += NbOfAtomsPerVolume[elm] * c 222 // SIGMA += NbOfAtomsPerVolume[elm] * crossSection; 226 SIGMA = NbOfMoleculesPerVolume * crossSect 223 SIGMA = NbOfMoleculesPerVolume * crossSection; 227 if( verboseLevel >= 2 ) G4cout << GetName( 224 if( verboseLevel >= 2 ) G4cout << GetName() << " ADDING SIGMA " << SIGMA << " NAtoms " << NbOfMoleculesPerVolume 228 << " Bin " << Bin << " TOTAL " << a 225 << " Bin " << Bin << " TOTAL " << aa << " " << bb 229 << " XS " << integralXS[theXSType][ 226 << " XS " << integralXS[theXSType][Bin] << " " << integralXS[theXSType][Bin+1] << G4endl; 230 } 227 } 231 228 232 //-} 229 //-} 233 230 234 fValue = SIGMA > DBL_MIN ? 1./SIGMA : DB 231 fValue = SIGMA > DBL_MIN ? 1./SIGMA : DBL_MAX; 235 } 232 } 236 233 237 ptrVector->PutValue(ii, fValue); 234 ptrVector->PutValue(ii, fValue); 238 if( verboseLevel >= 2 ) G4cout << GetName( 235 if( verboseLevel >= 2 ) G4cout << GetName() << " BUILDXS " << ii << " : " << LowEdgeEnergy << " = " << fValue << G4endl; 239 } 236 } 240 237 241 theMeanFreePathTable->insertAt( matIdx , ptr 238 theMeanFreePathTable->insertAt( matIdx , ptrVector ) ; 242 } 239 } 243 240 244 241 245 //....oooOO0OOooo........oooOO0OOooo........oo 242 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 246 G4double G4VLEPTSModel::SampleAngle(const G4Ma 243 G4double G4VLEPTSModel::SampleAngle(const G4Material* aMaterial, G4double e, G4double el) 247 { 244 { 248 G4double x; 245 G4double x; 249 246 250 if( e < 10001) { 247 if( e < 10001) { 251 x = theDiffXS[aMaterial]->SampleAngleMT(e, 248 x = theDiffXS[aMaterial]->SampleAngleMT(e, el); 252 } 249 } 253 else { 250 else { 254 G4double Ei = e; 251 G4double Ei = e; //incidente 255 G4double Ed = e -el; 252 G4double Ed = e -el; //dispersado 256 253 257 G4double Pi = std::sqrt( std::pow( (Ei/27. 254 G4double Pi = std::sqrt( std::pow( (Ei/27.2/137),2) +2*Ei/27.2); //incidente 258 G4double Pd = std::sqrt( std::pow( (Ed/27. 255 G4double Pd = std::sqrt( std::pow( (Ed/27.2/137),2) +2*Ed/27.2); //dispersado 259 256 260 G4double Kmin = Pi - Pd; 257 G4double Kmin = Pi - Pd; 261 G4double Kmax = Pi + Pd; 258 G4double Kmax = Pi + Pd; 262 259 263 G4double KR = theRMTDistr[aMaterial]->Samp 260 G4double KR = theRMTDistr[aMaterial]->Sample(Kmin, Kmax); //sorteo mom. transf. 264 261 265 G4double co = (Pi*Pi + Pd*Pd - KR*KR) / (2 262 G4double co = (Pi*Pi + Pd*Pd - KR*KR) / (2*Pi*Pd); //cos ang. disp. 266 if( co > 1. ) co = 1.; 263 if( co > 1. ) co = 1.; 267 x = std::acos(co); //*360/twopi; 264 x = std::acos(co); //*360/twopi; //ang. dispers. 268 } 265 } 269 return(x); 266 return(x); 270 } 267 } 271 268 272 //....oooOO0OOooo........oooOO0OOooo........oo 269 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 273 G4ThreeVector G4VLEPTSModel::SampleNewDirectio 270 G4ThreeVector G4VLEPTSModel::SampleNewDirection(const G4Material* aMaterial, G4ThreeVector P0Dir, G4double e, G4double el) { 274 271 275 G4double x = SampleAngle(aMaterial, e, el); 272 G4double x = SampleAngle(aMaterial, e, el); 276 273 277 G4double cosTeta = std::cos(x); //*twopi/360 274 G4double cosTeta = std::cos(x); //*twopi/360.0); 278 G4double sinTeta = std::sqrt(1.0-cosTeta*cos 275 G4double sinTeta = std::sqrt(1.0-cosTeta*cosTeta); 279 G4double Phi = CLHEP::twopi * G4UniformR 276 G4double Phi = CLHEP::twopi * G4UniformRand() ; 280 G4double dirx = sinTeta*std::cos(Phi) , d 277 G4double dirx = sinTeta*std::cos(Phi) , diry = sinTeta*std::sin(Phi) , dirz = cosTeta ; 281 278 282 G4ThreeVector P1Dir(dirx, diry, dirz); 279 G4ThreeVector P1Dir(dirx, diry, dirz); 283 #ifdef DEBUG_LEPTS 280 #ifdef DEBUG_LEPTS 284 if( verboseLevel >= 2 ) G4cout << " G4VLEPTS 281 if( verboseLevel >= 2 ) G4cout << " G4VLEPTSModel::SampleNewDirection " <<P1Dir << G4endl; 285 #endif 282 #endif 286 P1Dir.rotateUz(P0Dir); 283 P1Dir.rotateUz(P0Dir); 287 #ifdef DEBUG_LEPTS 284 #ifdef DEBUG_LEPTS 288 if( verboseLevel >= 2 ) G4cout << " G4VLEPTS 285 if( verboseLevel >= 2 ) G4cout << " G4VLEPTSModel::SampleNewDirection rotated " <<P1Dir << " " << P0Dir << G4endl; 289 #endif 286 #endif 290 287 291 return(P1Dir); 288 return(P1Dir); 292 } 289 } 293 290 294 291 295 //....oooOO0OOooo........oooOO0OOooo........oo 292 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 296 G4ThreeVector G4VLEPTSModel::SampleNewDirectio 293 G4ThreeVector G4VLEPTSModel::SampleNewDirection(G4ThreeVector P0Dir, G4double x) 297 { 294 { 298 G4double cosTeta = std::cos(x); //*twopi/360 295 G4double cosTeta = std::cos(x); //*twopi/360.0); 299 G4double sinTeta = std::sqrt(1.0-cosTeta*cos 296 G4double sinTeta = std::sqrt(1.0-cosTeta*cosTeta); 300 G4double Phi = CLHEP::twopi * G4UniformR 297 G4double Phi = CLHEP::twopi * G4UniformRand() ; 301 G4double dirx = sinTeta*std::cos(Phi) , d 298 G4double dirx = sinTeta*std::cos(Phi) , diry = sinTeta*std::sin(Phi) , dirz = cosTeta ; 302 299 303 G4ThreeVector P1Dir( dirx,diry,dirz ); 300 G4ThreeVector P1Dir( dirx,diry,dirz ); 304 P1Dir.rotateUz(P0Dir); 301 P1Dir.rotateUz(P0Dir); 305 302 306 return(P1Dir); 303 return(P1Dir); 307 } 304 } 308 305 309 306 310 //....oooOO0OOooo........oooOO0OOooo........oo 307 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 311 G4double G4VLEPTSModel::SampleEnergyLoss(const 308 G4double G4VLEPTSModel::SampleEnergyLoss(const G4Material* aMaterial, G4double eMin, G4double eMax) 312 { 309 { 313 G4double el; 310 G4double el; 314 el = theElostDistr[aMaterial]->Sample(eMin/C 311 el = theElostDistr[aMaterial]->Sample(eMin/CLHEP::eV, eMax/CLHEP::eV)*CLHEP::eV; 315 312 316 #ifdef DEBUG_LEPTS 313 #ifdef DEBUG_LEPTS 317 if( verboseLevel >= 2 ) G4cout << aMaterial- 314 if( verboseLevel >= 2 ) G4cout << aMaterial->GetName() <<"SampleEnergyLoss/eV " << el/CLHEP::eV << " eMax/eV " << eMax/CLHEP::eV << " " 318 << " " << GetName() << G4endl; 315 << " " << GetName() << G4endl; 319 #endif 316 #endif 320 return el; 317 return el; 321 } 318 } 322 319 323 320 324 //....oooOO0OOooo........oooOO0OOooo........oo 321 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 325 G4bool G4VLEPTSModel::ReadParam(const G4String << 322 G4bool G4VLEPTSModel::ReadParam(G4String fnParam, const G4Material* aMaterial ) 326 { 323 { 327 std::ifstream fin(fnParam); 324 std::ifstream fin(fnParam); 328 if (!fin.is_open()) { 325 if (!fin.is_open()) { 329 G4Exception("G4VLEPTSModel::ReadParam", 326 G4Exception("G4VLEPTSModel::ReadParam", 330 "", 327 "", 331 JustWarning, 328 JustWarning, 332 (G4String("File not found: ")+ fnParam).c_ 329 (G4String("File not found: ")+ fnParam).c_str()); 333 return false; 330 return false; 334 } 331 } 335 332 336 G4double IonisPot, IonisPotInt; 333 G4double IonisPot, IonisPotInt; 337 334 338 fin >> IonisPot >> IonisPotInt; 335 fin >> IonisPot >> IonisPotInt; 339 if( verboseLevel >= 1 ) G4cout << "Read para 336 if( verboseLevel >= 1 ) G4cout << "Read param (" << fnParam << ")\t IonisPot: " << IonisPot 340 << " IonisPotInt: " << IonisPotInt << G4en 337 << " IonisPotInt: " << IonisPotInt << G4endl; 341 338 342 theIonisPot[aMaterial] = IonisPot * CLHEP::e 339 theIonisPot[aMaterial] = IonisPot * CLHEP::eV; 343 theIonisPotInt[aMaterial] = IonisPotInt * CL 340 theIonisPotInt[aMaterial] = IonisPotInt * CLHEP::eV; 344 341 345 G4double MolecularMass = 0; 342 G4double MolecularMass = 0; 346 auto nelem = (G4int)aMaterial->GetNumberOfE << 343 size_t nelem = aMaterial->GetNumberOfElements(); 347 const G4int* atomsV = aMaterial->GetAtomsVe 344 const G4int* atomsV = aMaterial->GetAtomsVector(); 348 for( G4int ii = 0; ii < nelem; ++ii ) { << 345 for( size_t ii = 0; ii < nelem; ii++ ) { 349 MolecularMass += aMaterial->GetElement(ii) 346 MolecularMass += aMaterial->GetElement(ii)->GetA()*atomsV[ii]/CLHEP::g; 350 // G4cout << " MMASS1 " << mmass/CLHEP: 347 // G4cout << " MMASS1 " << mmass/CLHEP::g << " " << aMaterial->GetElement(ii)->GetName() << " " << aMaterial->GetElement(ii)->GetA()/CLHEP::g << G4endl; 351 } 348 } 352 // G4cout << " MMASS " << MolecularMass << 349 // G4cout << " MMASS " << MolecularMass << " " << MolecularMass*CLHEP::g << " ME " << mmass << " " << mmass/CLHEP::g << G4endl; 353 theMolecularMass[aMaterial] = MolecularMass* 350 theMolecularMass[aMaterial] = MolecularMass* CLHEP::g/CLHEP::mole; 354 // theMolecularMass[aMaterial] = aMaterial- 351 // theMolecularMass[aMaterial] = aMaterial->GetMassOfMolecule()*CLHEP::Avogadro; // Material mixtures do not calculate molecular mass 355 352 356 if( verboseLevel >= 1) G4cout << " IonisPot: 353 if( verboseLevel >= 1) G4cout << " IonisPot: " << IonisPot/CLHEP::eV << " eV " 357 << " IonisPotInt: " << IonisPotInt/CLH 354 << " IonisPotInt: " << IonisPotInt/CLHEP::eV << " eV" 358 << " MolecularMass " << MolecularMass/ 355 << " MolecularMass " << MolecularMass/(CLHEP::g/CLHEP::mole) << " g/mole" << G4endl; 359 356 360 return true; 357 return true; 361 } 358 } 362 359 363 //....oooOO0OOooo........oooOO0OOooo........oo 360 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 364 std::map< G4int, std::vector<G4double> > G4VLE << 361 std::map< G4int, std::vector<G4double> > G4VLEPTSModel::ReadIXS(G4String fnIXS, const G4Material* aMaterial ) 365 { 362 { 366 std::map< G4int, std::vector<G4double> > int 363 std::map< G4int, std::vector<G4double> > integralXS; // process type - energy 367 //G4cout << "fnIXS (" << fnIXS << ")" << G4e 364 //G4cout << "fnIXS (" << fnIXS << ")" << G4endl; 368 365 369 std::ifstream fin(fnIXS); 366 std::ifstream fin(fnIXS); 370 if (!fin.is_open()) { 367 if (!fin.is_open()) { 371 G4Exception("G4VLEPTSModel::ReadIXS", 368 G4Exception("G4VLEPTSModel::ReadIXS", 372 "", 369 "", 373 JustWarning, 370 JustWarning, 374 (G4String("File not found: ")+ fnIXS).c_st 371 (G4String("File not found: ")+ fnIXS).c_str()); 375 return integralXS; 372 return integralXS; 376 } 373 } 377 374 378 G4int nXSdat, nXSsub; 375 G4int nXSdat, nXSsub; 379 fin >> nXSdat >> nXSsub; 376 fin >> nXSdat >> nXSsub; 380 if( verboseLevel >= 1 ) G4cout << "Read IXS 377 if( verboseLevel >= 1 ) G4cout << "Read IXS (" << fnIXS << ")\t nXSdat: " << nXSdat 381 << " nXSsub: " << nXSsub << G4endl; 378 << " nXSsub: " << nXSsub << G4endl; 382 theNXSdat[aMaterial] = nXSdat; 379 theNXSdat[aMaterial] = nXSdat; 383 theNXSsub[aMaterial] = nXSsub; 380 theNXSsub[aMaterial] = nXSsub; 384 381 385 G4double xsdat; 382 G4double xsdat; 386 for (G4int ip=0; ip<=nXSsub; ip++) { 383 for (G4int ip=0; ip<=nXSsub; ip++) { 387 integralXS[ip].push_back(0.); 384 integralXS[ip].push_back(0.); 388 } 385 } 389 for (G4int ie=1; ie<=nXSdat; ie++) { 386 for (G4int ie=1; ie<=nXSdat; ie++) { 390 for (G4int ip=0; ip<=nXSsub; ip++) { 387 for (G4int ip=0; ip<=nXSsub; ip++) { 391 fin >> xsdat; 388 fin >> xsdat; 392 integralXS[ip].push_back(xsdat); 389 integralXS[ip].push_back(xsdat); 393 if( verboseLevel >= 3 ) G4cout << GetNa 390 if( verboseLevel >= 3 ) G4cout << GetName() << " FILL IXS " << ip << " " << ie << " = " << integralXS[ip][ie] << " " << xsdat << G4endl; 394 // xsdat 1e-16*cm2 391 // xsdat 1e-16*cm2 395 } 392 } 396 } 393 } 397 fin.close(); 394 fin.close(); 398 395 399 return integralXS; 396 return integralXS; 400 } 397 } 401 398 402 399