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 // >> 27 // >> 28 // G4ForwardXrayTR class -- implementation file >> 29 // 26 // History: 30 // History: 27 // 1st version 11.09.97 V. Grichine (Vladimir. 31 // 1st version 11.09.97 V. Grichine (Vladimir.Grichine@cern.ch ) 28 // 2nd version 17.12.97 V. Grichine 32 // 2nd version 17.12.97 V. Grichine 29 // 17-09-01, migration of Materials to pure ST 33 // 17-09-01, migration of Materials to pure STL (mma) 30 // 10-03-03, migration to "cut per region" (V. 34 // 10-03-03, migration to "cut per region" (V.Ivanchenko) 31 // 03.06.03, V.Ivanchenko fix compilation warn 35 // 03.06.03, V.Ivanchenko fix compilation warnings 32 36 33 #include "G4ForwardXrayTR.hh" 37 #include "G4ForwardXrayTR.hh" 34 38 35 #include "globals.hh" 39 #include "globals.hh" 36 #include "G4Gamma.hh" << 37 #include "G4GeometryTolerance.hh" << 38 #include "G4Material.hh" << 39 #include "G4PhysicalConstants.hh" 40 #include "G4PhysicalConstants.hh" 40 #include "G4PhysicsLinearVector.hh" << 41 #include "G4SystemOfUnits.hh" 41 #include "G4PhysicsLogVector.hh" << 42 #include "G4Poisson.hh" >> 43 #include "G4Material.hh" 42 #include "G4PhysicsTable.hh" 44 #include "G4PhysicsTable.hh" 43 #include "G4PhysicsVector.hh" 45 #include "G4PhysicsVector.hh" 44 #include "G4Poisson.hh" << 46 #include "G4PhysicsLinearVector.hh" >> 47 #include "G4PhysicsLogVector.hh" 45 #include "G4ProductionCutsTable.hh" 48 #include "G4ProductionCutsTable.hh" 46 #include "G4SystemOfUnits.hh" << 49 #include "G4GeometryTolerance.hh" 47 #include "G4PhysicsModelCatalog.hh" << 50 >> 51 // Initialization of local constants >> 52 >> 53 G4int G4ForwardXrayTR::fSympsonNumber = 100; >> 54 >> 55 G4double G4ForwardXrayTR::fTheMinEnergyTR = 1.0*keV; >> 56 G4double G4ForwardXrayTR::fTheMaxEnergyTR = 100.0*keV; >> 57 G4double G4ForwardXrayTR::fTheMaxAngle = 1.0e-3; >> 58 G4double G4ForwardXrayTR::fTheMinAngle = 5.0e-6; >> 59 G4int G4ForwardXrayTR::fBinTR = 50; >> 60 >> 61 G4double G4ForwardXrayTR::fMinProtonTkin = 100.0*GeV; >> 62 G4double G4ForwardXrayTR::fMaxProtonTkin = 100.0*TeV; >> 63 G4int G4ForwardXrayTR::fTotBin = 50; >> 64 >> 65 G4double G4ForwardXrayTR::fPlasmaCof = 4.0*pi*fine_structure_const* >> 66 hbarc*hbarc*hbarc/electron_mass_c2; >> 67 >> 68 G4double G4ForwardXrayTR::fCofTR = fine_structure_const/pi; >> 69 48 70 49 ////////////////////////////////////////////// 71 ////////////////////////////////////////////////////////////////////// 50 // 72 // 51 // Constructor for creation of physics tables 73 // Constructor for creation of physics tables (angle and energy TR 52 // distributions) for a couple of selected mat 74 // distributions) for a couple of selected materials. 53 // 75 // 54 // Recommended for use in applications with ma 76 // Recommended for use in applications with many materials involved, 55 // when only few (usually couple) materials ar 77 // when only few (usually couple) materials are interested for generation 56 // of TR on the interface between them 78 // of TR on the interface between them 57 G4ForwardXrayTR::G4ForwardXrayTR(const G4Strin << 79 58 const G4Strin << 80 59 const G4Strin << 81 G4ForwardXrayTR:: 60 : G4TransitionRadiation(processName) << 82 G4ForwardXrayTR( const G4String& matName1, // G4Material* pMat1, >> 83 const G4String& matName2, // G4Material* pMat2, >> 84 const G4String& processName ) >> 85 : G4TransitionRadiation(processName) 61 { 86 { 62 secID = G4PhysicsModelCatalog::GetModelID("m << 87 fPtrGamma = nullptr; 63 fPtrGamma = nullptr; << 64 fGammaCutInKineticEnergy = nullptr; 88 fGammaCutInKineticEnergy = nullptr; 65 fGammaTkinCut = fMinEnergyTR = fMaxEnergyTR << 89 fGammaTkinCut = fMinEnergyTR = fMaxEnergyTR = fMaxThetaTR = fGamma = 66 fGamma = fSigma1 = fSigma2 = 0.0; << 90 fSigma1 = fSigma2 = 0.0; 67 fAngleDistrTable = nullptr; << 91 fAngleDistrTable = nullptr; 68 fEnergyDistrTable = nullptr; << 92 fEnergyDistrTable = nullptr; 69 fMatIndex1 = fMatIndex2 = 0; 93 fMatIndex1 = fMatIndex2 = 0; 70 94 71 // Proton energy vector initialization 95 // Proton energy vector initialization 72 fProtonEnergyVector = << 96 // 73 new G4PhysicsLogVector(fMinProtonTkin, fMa << 97 fProtonEnergyVector = new G4PhysicsLogVector(fMinProtonTkin, >> 98 fMaxProtonTkin, fTotBin ); 74 G4int iMat; 99 G4int iMat; 75 const G4ProductionCutsTable* theCoupleTable << 100 const G4ProductionCutsTable* theCoupleTable= 76 G4ProductionCutsTable::GetProductionCutsTa << 101 G4ProductionCutsTable::GetProductionCutsTable(); 77 G4int numOfCouples = (G4int)theCoupleTable-> << 102 G4int numOfCouples = theCoupleTable->GetTableSize(); 78 103 79 G4bool build = true; 104 G4bool build = true; 80 105 81 for(iMat = 0; iMat < numOfCouples; ++iMat) << 106 for(iMat=0;iMat<numOfCouples;iMat++) // check first material name 82 { 107 { 83 const G4MaterialCutsCouple* couple = << 108 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(iMat); 84 theCoupleTable->GetMaterialCutsCouple(iM << 109 if( matName1 == couple->GetMaterial()->GetName() ) 85 if(matName1 == couple->GetMaterial()->GetN << 86 { 110 { 87 fMatIndex1 = couple->GetIndex(); 111 fMatIndex1 = couple->GetIndex(); 88 break; 112 break; 89 } 113 } 90 } 114 } 91 if(iMat == numOfCouples) 115 if(iMat == numOfCouples) 92 { 116 { 93 G4Exception("G4ForwardXrayTR::G4ForwardXra << 117 G4Exception("G4ForwardXrayTR::G4ForwardXrayTR", "ForwardXrayTR01", JustWarning, "Invalid first material name in G4ForwardXrayTR constructor!"); 94 JustWarning, << 95 "Invalid first material name i << 96 build = false; 118 build = false; 97 } 119 } 98 120 99 if(build) << 121 if(build) { 100 { << 122 for(iMat=0;iMat<numOfCouples;iMat++) // check second material name 101 for(iMat = 0; iMat < numOfCouples; ++iMat) << 102 { << 103 const G4MaterialCutsCouple* couple = << 104 theCoupleTable->GetMaterialCutsCouple( << 105 if(matName2 == couple->GetMaterial()->Ge << 106 { 123 { 107 fMatIndex2 = couple->GetIndex(); << 124 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(iMat); 108 break; << 125 if( matName2 == couple->GetMaterial()->GetName() ) >> 126 { >> 127 fMatIndex2 = couple->GetIndex(); >> 128 break; >> 129 } 109 } 130 } 110 } << 111 if(iMat == numOfCouples) 131 if(iMat == numOfCouples) 112 { << 132 { 113 G4Exception( << 133 G4Exception("G4ForwardXrayTR::G4ForwardXrayTR", "ForwardXrayTR02", JustWarning, "Invalid second material name in G4ForwardXrayTR constructor!"); 114 "G4ForwardXrayTR::G4ForwardXrayTR", "F << 134 build = false; 115 "Invalid second material name in G4For << 135 } 116 build = false; << 117 } << 118 } << 119 if(build) << 120 { << 121 BuildXrayTRtables(); << 122 } 136 } >> 137 // G4cout<<"G4ForwardXray constructor is called"<<G4endl; >> 138 if(build) { BuildXrayTRtables(); } 123 } 139 } 124 140 125 ////////////////////////////////////////////// 141 ///////////////////////////////////////////////////////////////////////// >> 142 // 126 // Constructor used by X-ray transition radiat 143 // Constructor used by X-ray transition radiation parametrisation models 127 G4ForwardXrayTR::G4ForwardXrayTR(const G4Strin << 144 128 : G4TransitionRadiation(processName) << 145 G4ForwardXrayTR:: >> 146 G4ForwardXrayTR( const G4String& processName ) >> 147 : G4TransitionRadiation(processName) 129 { 148 { 130 fPtrGamma = nullptr; << 149 fPtrGamma = nullptr; 131 fGammaCutInKineticEnergy = nullptr; 150 fGammaCutInKineticEnergy = nullptr; 132 fGammaTkinCut = fMinEnergyTR = fMaxEnergyTR << 151 fGammaTkinCut = fMinEnergyTR = fMaxEnergyTR = fMaxThetaTR = fGamma = 133 fGamma = fSigma1 = fSigma2 = 0.0; << 152 fSigma1 = fSigma2 = 0.0; 134 fAngleDistrTable = nullptr; << 153 fAngleDistrTable = nullptr; 135 fEnergyDistrTable = nullptr; << 154 fEnergyDistrTable = nullptr; 136 fMatIndex1 = fMatIndex2 = 0; 155 fMatIndex1 = fMatIndex2 = 0; 137 156 138 // Proton energy vector initialization 157 // Proton energy vector initialization 139 fProtonEnergyVector = << 158 // 140 new G4PhysicsLogVector(fMinProtonTkin, fMa << 159 fProtonEnergyVector = new G4PhysicsLogVector(fMinProtonTkin, >> 160 fMaxProtonTkin, fTotBin ); 141 } 161 } 142 162 >> 163 143 ////////////////////////////////////////////// 164 ////////////////////////////////////////////////////////////////////// >> 165 // 144 // Destructor 166 // Destructor >> 167 // >> 168 145 G4ForwardXrayTR::~G4ForwardXrayTR() 169 G4ForwardXrayTR::~G4ForwardXrayTR() 146 { 170 { 147 delete fAngleDistrTable; 171 delete fAngleDistrTable; 148 delete fEnergyDistrTable; 172 delete fEnergyDistrTable; 149 delete fProtonEnergyVector; 173 delete fProtonEnergyVector; 150 } 174 } 151 175 152 void G4ForwardXrayTR::ProcessDescription(std:: << 176 G4double G4ForwardXrayTR::GetMeanFreePath(const G4Track&, 153 { << 177 G4double, 154 out << "Simulation of forward X-ray transiti << 178 G4ForceCondition* condition) 155 "relativistic charged particles cross << 156 "two materials.\n"; << 157 } << 158 << 159 G4double G4ForwardXrayTR::GetMeanFreePath(cons << 160 G4Fo << 161 { 179 { 162 *condition = Forced; 180 *condition = Forced; 163 return DBL_MAX; // so TR doesn't limit mean << 181 return DBL_MAX; // so TR doesn't limit mean free path 164 } 182 } 165 183 166 ////////////////////////////////////////////// 184 ////////////////////////////////////////////////////////////////////////////// >> 185 // 167 // Build physics tables for energy and angular 186 // Build physics tables for energy and angular distributions of X-ray TR photon >> 187 168 void G4ForwardXrayTR::BuildXrayTRtables() 188 void G4ForwardXrayTR::BuildXrayTRtables() 169 { 189 { 170 G4int iMat, jMat, iTkin, iTR, iPlace; 190 G4int iMat, jMat, iTkin, iTR, iPlace; 171 const G4ProductionCutsTable* theCoupleTable << 191 const G4ProductionCutsTable* theCoupleTable= 172 G4ProductionCutsTable::GetProductionCutsTa << 192 G4ProductionCutsTable::GetProductionCutsTable(); 173 G4int numOfCouples = (G4int)theCoupleTable-> << 193 G4int numOfCouples = theCoupleTable->GetTableSize(); 174 194 175 fGammaCutInKineticEnergy = theCoupleTable->G 195 fGammaCutInKineticEnergy = theCoupleTable->GetEnergyCutsVector(idxG4GammaCut); 176 196 177 fAngleDistrTable = new G4PhysicsTable(2 * f << 197 fAngleDistrTable = new G4PhysicsTable(2*fTotBin); 178 fEnergyDistrTable = new G4PhysicsTable(2 * f << 198 fEnergyDistrTable = new G4PhysicsTable(2*fTotBin); >> 199 179 200 180 for(iMat = 0; iMat < numOfCouples; << 201 for(iMat=0;iMat<numOfCouples;iMat++) // loop over pairs of different materials 181 ++iMat) // loop over pairs of different << 182 { 202 { 183 if(iMat != fMatIndex1 && iMat != fMatIndex << 203 if( iMat != fMatIndex1 && iMat != fMatIndex2 ) continue; 184 continue; << 185 204 186 for(jMat = 0; jMat < numOfCouples; ++jMat) << 205 for(jMat=0;jMat<numOfCouples;jMat++) // transition iMat -> jMat !!! 187 { 206 { 188 if(iMat == jMat || (jMat != fMatIndex1 & << 207 if( iMat == jMat || ( jMat != fMatIndex1 && jMat != fMatIndex2 ) ) 189 { 208 { 190 continue; 209 continue; 191 } 210 } 192 else 211 else 193 { 212 { 194 const G4MaterialCutsCouple* iCouple = << 213 const G4MaterialCutsCouple* iCouple = theCoupleTable->GetMaterialCutsCouple(iMat); 195 theCoupleTable->GetMaterialCutsCoupl << 214 const G4MaterialCutsCouple* jCouple = theCoupleTable->GetMaterialCutsCouple(jMat); 196 const G4MaterialCutsCouple* jCouple = << 197 theCoupleTable->GetMaterialCutsCoupl << 198 const G4Material* mat1 = iCouple->GetM 215 const G4Material* mat1 = iCouple->GetMaterial(); 199 const G4Material* mat2 = jCouple->GetM 216 const G4Material* mat2 = jCouple->GetMaterial(); 200 217 201 fSigma1 = fPlasmaCof * (mat1->GetElect << 218 fSigma1 = fPlasmaCof*(mat1->GetElectronDensity()); 202 fSigma2 = fPlasmaCof * (mat2->GetElect << 219 fSigma2 = fPlasmaCof*(mat2->GetElectronDensity()); >> 220 >> 221 // fGammaTkinCut = fGammaCutInKineticEnergy[jMat]; // TR photon in jMat ! 203 222 204 fGammaTkinCut = 0.0; 223 fGammaTkinCut = 0.0; 205 224 206 if(fGammaTkinCut > fTheMinEnergyTR) / << 225 if(fGammaTkinCut > fTheMinEnergyTR) // setting of min/max TR energies 207 { << 226 { 208 fMinEnergyTR = fGammaTkinCut; 227 fMinEnergyTR = fGammaTkinCut; 209 } << 228 } 210 else 229 else 211 { << 230 { 212 fMinEnergyTR = fTheMinEnergyTR; 231 fMinEnergyTR = fTheMinEnergyTR; 213 } << 232 } 214 if(fGammaTkinCut > fTheMaxEnergyTR) 233 if(fGammaTkinCut > fTheMaxEnergyTR) 215 { << 234 { 216 fMaxEnergyTR = 2.0 * fGammaTkinCut; << 235 fMaxEnergyTR = 2.0*fGammaTkinCut; // usually very low TR rate 217 } << 236 } 218 else 237 else 219 { << 238 { 220 fMaxEnergyTR = fTheMaxEnergyTR; 239 fMaxEnergyTR = fTheMaxEnergyTR; 221 } << 240 } 222 for(iTkin = 0; iTkin < fTotBin; ++iTki << 241 for(iTkin=0;iTkin<fTotBin;iTkin++) // Lorentz factor loop 223 { << 242 { 224 auto energyVector = << 243 G4PhysicsLogVector* 225 new G4PhysicsLogVector(fMinEnergyT << 244 energyVector = new G4PhysicsLogVector( fMinEnergyTR, >> 245 fMaxEnergyTR, >> 246 fBinTR ); 226 247 227 fGamma = 1.0 + (fProtonEnergyVector- << 248 fGamma = 1.0 + (fProtonEnergyVector-> 228 proton_mass_c2); << 249 GetLowEdgeEnergy(iTkin)/proton_mass_c2); 229 250 230 fMaxThetaTR = 10000.0 / (fGamma * fG << 251 fMaxThetaTR = 10000.0/(fGamma*fGamma); 231 252 232 if(fMaxThetaTR > fTheMaxAngle) 253 if(fMaxThetaTR > fTheMaxAngle) 233 { 254 { 234 fMaxThetaTR = fTheMaxAngle; 255 fMaxThetaTR = fTheMaxAngle; 235 } << 256 } 236 else 257 else 237 { << 258 { 238 if(fMaxThetaTR < fTheMinAngle) 259 if(fMaxThetaTR < fTheMinAngle) 239 { << 260 { 240 fMaxThetaTR = fTheMinAngle; 261 fMaxThetaTR = fTheMinAngle; 241 } << 262 } 242 } << 263 } 243 auto angleVector = << 264 // G4cout<<G4endl<<"fGamma = "<<fGamma<<" fMaxThetaTR = "<<fMaxThetaTR<<G4endl; 244 new G4PhysicsLinearVector(0.0, fMa << 265 G4PhysicsLinearVector* >> 266 angleVector = new G4PhysicsLinearVector( 0.0, >> 267 fMaxThetaTR, >> 268 fBinTR ); 245 G4double energySum = 0.0; 269 G4double energySum = 0.0; 246 G4double angleSum = 0.0; 270 G4double angleSum = 0.0; 247 271 248 energyVector->PutValue(fBinTR - 1, e << 272 energyVector->PutValue(fBinTR-1,energySum); 249 angleVector->PutValue(fBinTR - 1, an << 273 angleVector->PutValue(fBinTR-1,angleSum); 250 274 251 for(iTR = fBinTR - 2; iTR >= 0; --iT << 275 for(iTR=fBinTR-2;iTR>=0;iTR--) 252 { << 276 { 253 energySum += << 277 energySum += fCofTR*EnergySum(energyVector->GetLowEdgeEnergy(iTR), 254 fCofTR * EnergySum(energyVector- << 278 energyVector->GetLowEdgeEnergy(iTR+1)); 255 energyVector- << 279 256 << 280 angleSum += fCofTR*AngleSum(angleVector->GetLowEdgeEnergy(iTR), 257 angleSum += << 281 angleVector->GetLowEdgeEnergy(iTR+1)); 258 fCofTR * AngleSum(angleVector->G << 282 259 angleVector->G << 283 energyVector->PutValue(iTR,energySum); 260 << 284 angleVector ->PutValue(iTR,angleSum); 261 energyVector->PutValue(iTR, energy << 285 } 262 angleVector->PutValue(iTR, angleSu << 286 // G4cout<<"sumE = "<<energySum<<"; sumA = "<<angleSum<<G4endl; 263 } << 264 287 265 if(jMat < iMat) 288 if(jMat < iMat) 266 { << 289 { 267 iPlace = fTotBin + iTkin; << 290 iPlace = fTotBin+iTkin; // (iMat*(numOfMat-1)+jMat)* 268 } << 291 } 269 else // jMat > iMat right part of m << 292 else // jMat > iMat right part of matrices (jMat-1) ! 270 { << 293 { 271 iPlace = iTkin; << 294 iPlace = iTkin; // (iMat*(numOfMat-1)+jMat-1)*fTotBin+ 272 } << 295 } 273 fEnergyDistrTable->insertAt(iPlace, << 296 fEnergyDistrTable->insertAt(iPlace,energyVector); 274 fAngleDistrTable->insertAt(iPlace, a << 297 fAngleDistrTable->insertAt(iPlace,angleVector); 275 } // iTkin << 298 } // iTkin 276 } // jMat != iMat << 299 } // jMat != iMat 277 } // jMat << 300 } // jMat 278 } // iMat << 301 } // iMat >> 302 // G4cout<<"G4ForwardXrayTR::BuildXrayTRtables have been called"<<G4endl; 279 } 303 } 280 304 281 ////////////////////////////////////////////// 305 /////////////////////////////////////////////////////////////////////// 282 // 306 // 283 // This function returns the spectral and angl 307 // This function returns the spectral and angle density of TR quanta 284 // in X-ray energy region generated forward wh 308 // in X-ray energy region generated forward when a relativistic 285 // charged particle crosses interface between 309 // charged particle crosses interface between two materials. 286 // The high energy small theta approximation i 310 // The high energy small theta approximation is applied. 287 // (matter1 -> matter2) 311 // (matter1 -> matter2) 288 // varAngle =2* (1 - std::cos(Theta)) or appro 312 // varAngle =2* (1 - std::cos(Theta)) or approximately = Theta*Theta 289 // 313 // 290 G4double G4ForwardXrayTR::SpectralAngleTRdensi << 314 291 << 315 G4double 292 { << 316 G4ForwardXrayTR::SpectralAngleTRdensity( G4double energy, 293 G4double formationLength1, formationLength2; << 317 G4double varAngle ) const 294 formationLength1 = << 318 { 295 1.0 / (1.0 / (fGamma * fGamma) + fSigma1 / << 319 G4double formationLength1, formationLength2; 296 formationLength2 = << 320 formationLength1 = 1.0/ 297 1.0 / (1.0 / (fGamma * fGamma) + fSigma2 / << 321 (1.0/(fGamma*fGamma) 298 return (varAngle / energy) * (formationLengt << 322 + fSigma1/(energy*energy) 299 (formationLength1 - formationLength2) << 323 + varAngle); >> 324 formationLength2 = 1.0/ >> 325 (1.0/(fGamma*fGamma) >> 326 + fSigma2/(energy*energy) >> 327 + varAngle); >> 328 return (varAngle/energy)*(formationLength1 - formationLength2) >> 329 *(formationLength1 - formationLength2); >> 330 300 } 331 } 301 332 >> 333 302 ////////////////////////////////////////////// 334 ////////////////////////////////////////////////////////////////// >> 335 // 303 // Analytical formula for angular density of X 336 // Analytical formula for angular density of X-ray TR photons 304 G4double G4ForwardXrayTR::AngleDensity(G4doubl << 337 // >> 338 >> 339 G4double G4ForwardXrayTR::AngleDensity( G4double energy, >> 340 G4double varAngle ) const 305 { 341 { 306 G4double x, x2, c, d, f, a2, b2, a4, b4; << 342 G4double x, x2, /*a, b,*/ c, d, f, a2, b2, a4, b4; 307 G4double cof1, cof2, cof3; 343 G4double cof1, cof2, cof3; 308 x = 1.0 / energy; << 344 x = 1.0/energy; 309 x2 = x * x; << 345 x2 = x*x; 310 c = 1.0 / fSigma1; << 346 c = 1.0/fSigma1; 311 d = 1.0 / fSigma2; << 347 d = 1.0/fSigma2; 312 f = (varAngle + 1.0 / (fGamma * fGamma)); << 348 f = (varAngle + 1.0/(fGamma*fGamma)); 313 a2 = c * f; << 349 a2 = c*f; 314 b2 = d * f; << 350 b2 = d*f; 315 a4 = a2 * a2; << 351 a4 = a2*a2; 316 b4 = b2 * b2; << 352 b4 = b2*b2; 317 cof1 = c * c * (0.5 / (a2 * (x2 + a2)) + 0.5 << 353 //a = std::sqrt(a2); 318 cof3 = d * d * (0.5 / (b2 * (x2 + b2)) + 0.5 << 354 // b = std::sqrt(b2); 319 cof2 = -c * d * << 355 cof1 = c*c*(0.5/(a2*(x2 +a2)) +0.5*std::log(x2/(x2 +a2))/a4); 320 (std::log(x2 / (x2 + b2)) / b2 - std: << 356 cof3 = d*d*(0.5/(b2*(x2 +b2)) +0.5*std::log(x2/(x2 +b2))/b4); 321 (a2 - b2); << 357 cof2 = -c*d*(std::log(x2/(x2 +b2))/b2 - std::log(x2/(x2 +a2))/a2)/(a2 - b2) ; 322 return -varAngle * (cof1 + cof2 + cof3); << 358 return -varAngle*(cof1 + cof2 + cof3); 323 } 359 } 324 360 325 ////////////////////////////////////////////// 361 ///////////////////////////////////////////////////////////////////// >> 362 // 326 // Definite integral of X-ray TR spectral-angl 363 // Definite integral of X-ray TR spectral-angle density from energy1 327 // to energy2 364 // to energy2 328 G4double G4ForwardXrayTR::EnergyInterval(G4dou << 365 // 329 G4dou << 366 >> 367 G4double G4ForwardXrayTR::EnergyInterval( G4double energy1, >> 368 G4double energy2, >> 369 G4double varAngle ) const 330 { 370 { 331 return AngleDensity(energy2, varAngle) - Ang << 371 return AngleDensity(energy2,varAngle) >> 372 - AngleDensity(energy1,varAngle); 332 } 373 } 333 374 334 ////////////////////////////////////////////// 375 ////////////////////////////////////////////////////////////////////// >> 376 // 335 // Integral angle distribution of X-ray TR pho 377 // Integral angle distribution of X-ray TR photons based on analytical 336 // formula for angle density 378 // formula for angle density 337 G4double G4ForwardXrayTR::AngleSum(G4double va << 379 // >> 380 >> 381 G4double G4ForwardXrayTR::AngleSum( G4double varAngle1, >> 382 G4double varAngle2 ) const 338 { 383 { 339 G4int i; 384 G4int i; 340 G4double h, sumEven = 0.0, sumOdd = 0.0; << 385 G4double h , sumEven = 0.0 , sumOdd = 0.0; 341 h = 0.5 * (varAngle2 - varAngle1) / fSympson << 386 h = 0.5*(varAngle2 - varAngle1)/fSympsonNumber; 342 for(i = 1; i < fSympsonNumber; ++i) << 387 for(i=1;i<fSympsonNumber;i++) 343 { << 388 { 344 sumEven += << 389 sumEven += EnergyInterval(fMinEnergyTR,fMaxEnergyTR,varAngle1 + 2*i*h ); 345 EnergyInterval(fMinEnergyTR, fMaxEnergyT << 390 sumOdd += EnergyInterval(fMinEnergyTR,fMaxEnergyTR, 346 sumOdd += << 391 varAngle1 + (2*i - 1)*h ); 347 EnergyInterval(fMinEnergyTR, fMaxEnergyT << 392 } 348 } << 393 sumOdd += EnergyInterval(fMinEnergyTR,fMaxEnergyTR, 349 sumOdd += EnergyInterval(fMinEnergyTR, fMaxE << 394 varAngle1 + (2*fSympsonNumber - 1)*h ); 350 varAngle1 + (2 * fS << 395 351 << 396 return h*(EnergyInterval(fMinEnergyTR,fMaxEnergyTR,varAngle1) 352 return h * << 397 + EnergyInterval(fMinEnergyTR,fMaxEnergyTR,varAngle2) 353 (EnergyInterval(fMinEnergyTR, fMaxEne << 398 + 4.0*sumOdd + 2.0*sumEven )/3.0; 354 EnergyInterval(fMinEnergyTR, fMaxEne << 355 2.0 * sumEven) / << 356 3.0; << 357 } 399 } 358 400 359 ////////////////////////////////////////////// 401 ///////////////////////////////////////////////////////////////////// >> 402 // 360 // Analytical Expression for spectral densit 403 // Analytical Expression for spectral density of Xray TR photons 361 // x = 2*(1 - std::cos(Theta)) ~ Theta^2 404 // x = 2*(1 - std::cos(Theta)) ~ Theta^2 362 G4double G4ForwardXrayTR::SpectralDensity(G4do << 405 // >> 406 >> 407 G4double G4ForwardXrayTR::SpectralDensity( G4double energy, >> 408 G4double x ) const 363 { 409 { 364 G4double a, b; 410 G4double a, b; 365 a = 1.0 / (fGamma * fGamma) + fSigma1 / (ene << 411 a = 1.0/(fGamma*fGamma) 366 b = 1.0 / (fGamma * fGamma) + fSigma2 / (ene << 412 + fSigma1/(energy*energy) ; 367 return ((a + b) * std::log((x + b) / (x + a) << 413 b = 1.0/(fGamma*fGamma) 368 b / (x + b)) / << 414 + fSigma2/(energy*energy) ; 369 energy; << 415 return ( (a + b)*std::log((x + b)/(x + a))/(a - b) >> 416 + a/(x + a) + b/(x + b) )/energy; >> 417 370 } 418 } 371 419 372 ////////////////////////////////////////////// 420 //////////////////////////////////////////////////////////////////// >> 421 // 373 // The spectral density in some angle interva 422 // The spectral density in some angle interval from varAngle1 to 374 // varAngle2 423 // varAngle2 375 G4double G4ForwardXrayTR::AngleInterval(G4doub << 424 // 376 G4doub << 425 >> 426 G4double G4ForwardXrayTR::AngleInterval( G4double energy, >> 427 G4double varAngle1, >> 428 G4double varAngle2 ) const 377 { 429 { 378 return SpectralDensity(energy, varAngle2) - << 430 return SpectralDensity(energy,varAngle2) 379 SpectralDensity(energy, varAngle1); << 431 - SpectralDensity(energy,varAngle1); 380 } 432 } 381 433 382 ////////////////////////////////////////////// 434 //////////////////////////////////////////////////////////////////// >> 435 // 383 // Integral spectral distribution of X-ray TR 436 // Integral spectral distribution of X-ray TR photons based on 384 // analytical formula for spectral density 437 // analytical formula for spectral density 385 G4double G4ForwardXrayTR::EnergySum(G4double e << 438 // >> 439 >> 440 G4double G4ForwardXrayTR::EnergySum( G4double energy1, >> 441 G4double energy2 ) const 386 { 442 { 387 G4int i; 443 G4int i; 388 G4double h, sumEven = 0.0, sumOdd = 0.0; << 444 G4double h , sumEven = 0.0 , sumOdd = 0.0; 389 h = 0.5 * (energy2 - energy1) / fSympsonNumb << 445 h = 0.5*(energy2 - energy1)/fSympsonNumber; 390 for(i = 1; i < fSympsonNumber; ++i) << 446 for(i=1;i<fSympsonNumber;i++) 391 { << 447 { 392 sumEven += AngleInterval(energy1 + 2 * i * << 448 sumEven += AngleInterval(energy1 + 2*i*h,0.0,fMaxThetaTR); 393 sumOdd += AngleInterval(energy1 + (2 * i - << 449 sumOdd += AngleInterval(energy1 + (2*i - 1)*h,0.0,fMaxThetaTR); 394 } << 450 } 395 sumOdd += << 451 sumOdd += AngleInterval(energy1 + (2*fSympsonNumber - 1)*h, 396 AngleInterval(energy1 + (2 * fSympsonNumbe << 452 0.0,fMaxThetaTR); 397 << 453 398 return h * << 454 return h*( AngleInterval(energy1,0.0,fMaxThetaTR) 399 (AngleInterval(energy1, 0.0, fMaxThet << 455 + AngleInterval(energy2,0.0,fMaxThetaTR) 400 AngleInterval(energy2, 0.0, fMaxThet << 456 + 4.0*sumOdd + 2.0*sumEven )/3.0; 401 2.0 * sumEven) / << 402 3.0; << 403 } 457 } 404 458 405 ////////////////////////////////////////////// 459 ///////////////////////////////////////////////////////////////////////// >> 460 // 406 // PostStepDoIt function for creation of forwa 461 // PostStepDoIt function for creation of forward X-ray photons in TR process 407 // on boundary between two materials with real << 462 // on boubndary between two materials with really different plasma energies >> 463 // >> 464 408 G4VParticleChange* G4ForwardXrayTR::PostStepDo 465 G4VParticleChange* G4ForwardXrayTR::PostStepDoIt(const G4Track& aTrack, 409 << 466 const G4Step& aStep) 410 { 467 { 411 aParticleChange.Initialize(aTrack); 468 aParticleChange.Initialize(aTrack); >> 469 // G4cout<<"call G4ForwardXrayTR::PostStepDoIt"<<G4endl; 412 G4int iMat, jMat, iTkin, iPlace, numOfTR, iT 470 G4int iMat, jMat, iTkin, iPlace, numOfTR, iTR, iTransfer; 413 471 414 G4double energyPos, anglePos, energyTR, thet 472 G4double energyPos, anglePos, energyTR, theta, phi, dirX, dirY, dirZ; 415 G4double W, W1, W2, E1, E2; 473 G4double W, W1, W2, E1, E2; 416 474 417 G4StepPoint* pPreStepPoint = aStep.GetPreSt 475 G4StepPoint* pPreStepPoint = aStep.GetPreStepPoint(); 418 G4StepPoint* pPostStepPoint = aStep.GetPostS 476 G4StepPoint* pPostStepPoint = aStep.GetPostStepPoint(); 419 G4double tol = << 477 G4double tol=0.5*G4GeometryTolerance::GetInstance()->GetSurfaceTolerance(); 420 0.5 * G4GeometryTolerance::GetInstance()-> << 421 478 422 if(pPostStepPoint->GetStepStatus() != fGeomB << 479 if (pPostStepPoint->GetStepStatus() != fGeomBoundary) 423 { 480 { 424 return G4VDiscreteProcess::PostStepDoIt(aT 481 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep); 425 } 482 } 426 if(aTrack.GetStepLength() <= tol) << 483 if (aTrack.GetStepLength() <= tol) 427 { 484 { 428 return G4VDiscreteProcess::PostStepDoIt(aT 485 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep); 429 } 486 } 430 // Arrived at boundary, so begin to try TR << 487 // Come on boundary, so begin to try TR 431 488 432 const G4MaterialCutsCouple* iCouple = pPreSt << 489 const G4MaterialCutsCouple* iCouple = pPreStepPoint ->GetPhysicalVolume()-> 433 ->Ge << 490 GetLogicalVolume()->GetMaterialCutsCouple(); 434 ->Ge << 491 const G4MaterialCutsCouple* jCouple = pPostStepPoint ->GetPhysicalVolume()-> 435 const G4MaterialCutsCouple* jCouple = pPostS << 492 GetLogicalVolume()->GetMaterialCutsCouple(); 436 ->Ge << 437 ->Ge << 438 const G4Material* iMaterial = iCouple->GetMa 493 const G4Material* iMaterial = iCouple->GetMaterial(); 439 const G4Material* jMaterial = jCouple->GetMa 494 const G4Material* jMaterial = jCouple->GetMaterial(); 440 iMat = iCouple->GetIn << 495 iMat = iCouple->GetIndex(); 441 jMat = jCouple->GetIn << 496 jMat = jCouple->GetIndex(); 442 497 443 // The case of equal or approximate (in term 498 // The case of equal or approximate (in terms of plasma energy) materials 444 // No TR photons ?! 499 // No TR photons ?! 445 500 446 if(iMat == jMat || << 501 if ( iMat == jMat 447 ((fMatIndex1 >= 0 && fMatIndex2 >= 0) && << 502 || ( (fMatIndex1 >= 0 && fMatIndex2 >= 0) 448 (iMat != fMatIndex1 && iMat != fMatIndex << 503 && ( iMat != fMatIndex1 && iMat != fMatIndex2 ) 449 (jMat != fMatIndex1 && jMat != fMatIndex << 504 && ( jMat != fMatIndex1 && jMat != fMatIndex2 ) ) 450 505 451 || iMaterial->GetState() == jMaterial->Ge << 506 || iMaterial->GetState() == jMaterial->GetState() 452 507 453 || (iMaterial->GetState() == kStateSolid << 508 ||(iMaterial->GetState() == kStateSolid && jMaterial->GetState() == kStateLiquid ) 454 jMaterial->GetState() == kStateLiquid << 455 509 456 || (iMaterial->GetState() == kStateLiquid << 510 ||(iMaterial->GetState() == kStateLiquid && jMaterial->GetState() == kStateSolid ) ) 457 jMaterial->GetState() == kStateSolid) << 458 { 511 { 459 return G4VDiscreteProcess::PostStepDoIt(aT 512 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep); 460 } 513 } 461 514 462 const G4DynamicParticle* aParticle = aTrack. 515 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle(); 463 G4double charge = aParticle->GetDefinition() 516 G4double charge = aParticle->GetDefinition()->GetPDGCharge(); 464 517 465 if(charge == 0.0) // Uncharged particle doe << 518 if(charge == 0.0) // Uncharged particle doesn't Generate TR photons 466 { 519 { 467 return G4VDiscreteProcess::PostStepDoIt(aT 520 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep); 468 } 521 } 469 // Now we are ready to Generate TR photons 522 // Now we are ready to Generate TR photons 470 523 471 G4double chargeSq = charge * charge; << 524 G4double chargeSq = charge*charge; 472 G4double kinEnergy = aParticle->GetKineticEn << 525 G4double kinEnergy = aParticle->GetKineticEnergy(); 473 G4double massRatio = << 526 G4double massRatio = proton_mass_c2/aParticle->GetDefinition()->GetPDGMass(); 474 proton_mass_c2 / aParticle->GetDefinition( << 527 G4double TkinScaled = kinEnergy*massRatio; 475 G4double TkinScaled = kinEnergy * massRatio; << 528 for(iTkin=0;iTkin<fTotBin;iTkin++) 476 for(iTkin = 0; iTkin < fTotBin; ++iTkin) << 477 { 529 { 478 if(TkinScaled < fProtonEnergyVector->GetLo << 530 if(TkinScaled < fProtonEnergyVector->GetLowEdgeEnergy(iTkin)) // <= ? 479 { 531 { 480 break; 532 break; 481 } 533 } 482 } 534 } 483 if(jMat < iMat) 535 if(jMat < iMat) 484 { 536 { 485 iPlace = fTotBin + iTkin - 1; << 537 iPlace = fTotBin + iTkin - 1; // (iMat*(numOfMat - 1) + jMat)* 486 } 538 } 487 else 539 else 488 { 540 { 489 iPlace = iTkin - 1; << 541 iPlace = iTkin - 1; // (iMat*(numOfMat - 1) + jMat - 1)*fTotBin + 490 } 542 } >> 543 // G4PhysicsVector* energyVector1 = (*fEnergyDistrTable)(iPlace) ; >> 544 // G4PhysicsVector* energyVector2 = (*fEnergyDistrTable)(iPlace + 1); >> 545 >> 546 // G4PhysicsVector* angleVector1 = (*fAngleDistrTable)(iPlace) ; >> 547 // G4PhysicsVector* angleVector2 = (*fAngleDistrTable)(iPlace + 1) ; 491 548 492 G4ParticleMomentum particleDir = aParticle-> 549 G4ParticleMomentum particleDir = aParticle->GetMomentumDirection(); 493 550 494 if(iTkin == fTotBin) // TR plato, try from << 551 if(iTkin == fTotBin) // TR plato, try from left 495 { 552 { 496 numOfTR = (G4int)G4Poisson( << 553 // G4cout<<iTkin<<" mean TR number = "<<( (*(*fEnergyDistrTable)(iPlace))(0) + 497 ((*(*fEnergyDistrTable)(iPlace))(0) + (* << 554 // (*(*fAngleDistrTable)(iPlace))(0) ) 498 chargeSq * 0.5); << 555 // *chargeSq*0.5<<G4endl; >> 556 >> 557 numOfTR = G4Poisson( ( (*(*fEnergyDistrTable)(iPlace))(0) + >> 558 (*(*fAngleDistrTable)(iPlace))(0) ) >> 559 *chargeSq*0.5 ); 499 if(numOfTR == 0) 560 if(numOfTR == 0) 500 { 561 { 501 return G4VDiscreteProcess::PostStepDoIt( 562 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep); 502 } 563 } 503 else 564 else 504 { 565 { >> 566 // G4cout<<"Number of X-ray TR photons = "<<numOfTR<<G4endl; >> 567 505 aParticleChange.SetNumberOfSecondaries(n 568 aParticleChange.SetNumberOfSecondaries(numOfTR); 506 569 507 for(iTR = 0; iTR < numOfTR; ++iTR) << 570 for(iTR=0;iTR<numOfTR;iTR++) 508 { 571 { 509 energyPos = (*(*fEnergyDistrTable)(iPl << 572 energyPos = (*(*fEnergyDistrTable)(iPlace))(0)*G4UniformRand(); 510 for(iTransfer = 0; iTransfer < fBinTR << 573 for(iTransfer=0;iTransfer<fBinTR-1;iTransfer++) 511 { << 574 { 512 if(energyPos >= (*(*fEnergyDistrTabl << 575 if(energyPos >= (*(*fEnergyDistrTable)(iPlace))(iTransfer)) break; 513 break; << 576 } 514 } << 515 energyTR = (*fEnergyDistrTable)(iPlace 577 energyTR = (*fEnergyDistrTable)(iPlace)->GetLowEdgeEnergy(iTransfer); 516 578 >> 579 // G4cout<<"energyTR = "<<energyTR/keV<<"keV"<<G4endl; >> 580 517 kinEnergy -= energyTR; 581 kinEnergy -= energyTR; 518 aParticleChange.ProposeEnergy(kinEnerg 582 aParticleChange.ProposeEnergy(kinEnergy); 519 583 520 anglePos = (*(*fAngleDistrTable)(iPlac << 584 anglePos = (*(*fAngleDistrTable)(iPlace))(0)*G4UniformRand(); 521 for(iTransfer = 0; iTransfer < fBinTR << 585 for(iTransfer=0;iTransfer<fBinTR-1;iTransfer++) 522 { << 586 { 523 if(anglePos > (*(*fAngleDistrTable)( << 587 if(anglePos > (*(*fAngleDistrTable)(iPlace))(iTransfer)) break; 524 break; << 588 } 525 } << 589 theta = std::sqrt((*fAngleDistrTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1)); 526 theta = std::sqrt( << 590 527 (*fAngleDistrTable)(iPlace)->GetLowE << 591 // G4cout<<iTransfer<<" : theta = "<<theta<<G4endl; 528 << 592 529 phi = twopi * G4UniformRand(); << 593 phi = twopi*G4UniformRand(); 530 dirX = std::sin(theta) * std::cos(phi) << 594 dirX = std::sin(theta)*std::cos(phi) ; 531 dirY = std::sin(theta) * std::sin(phi) << 595 dirY = std::sin(theta)*std::sin(phi) ; 532 dirZ = std::cos(theta); << 596 dirZ = std::cos(theta) ; 533 G4ThreeVector directionTR(dirX, dirY, << 597 G4ThreeVector directionTR(dirX,dirY,dirZ); 534 directionTR.rotateUz(particleDir); 598 directionTR.rotateUz(particleDir); 535 auto aPhotonTR = new G4DynamicParticle << 599 G4DynamicParticle* aPhotonTR = new G4DynamicParticle(G4Gamma::Gamma(), 536 << 600 directionTR, 537 // Create the G4Track << 601 energyTR ); 538 auto aSecondaryTrack = new G4Track(aPhotonTR << 602 aParticleChange.AddSecondary(aPhotonTR); 539 aSecondaryTrack->SetTouchableHandle(aStep.Ge << 540 aSecondaryTrack->SetParentID(aTrack.GetTrack << 541 aSecondaryTrack->SetCreatorModelID(secID); << 542 aParticleChange.AddSecondary(aSecondaryTrack << 543 } 603 } 544 } 604 } 545 } 605 } 546 else 606 else 547 { 607 { 548 if(iTkin == 0) // Tkin is too small, negl << 608 if(iTkin == 0) // Tkin is too small, neglect of TR photon generation 549 { 609 { 550 return G4VDiscreteProcess::PostStepDoIt( 610 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep); 551 } 611 } 552 else // general case: Tkin between two ve << 612 else // general case: Tkin between two vectors of the material 553 { 613 { 554 E1 = fProtonEnergyVector->GetLowEdgeEner 614 E1 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin - 1); 555 E2 = fProtonEnergyVector->GetLowEdgeEner << 615 E2 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin) ; 556 W = 1.0 / (E2 - E1); << 616 W = 1.0/(E2 - E1); 557 W1 = (E2 - TkinScaled) * W; << 617 W1 = (E2 - TkinScaled)*W; 558 W2 = (TkinScaled - E1) * W; << 618 W2 = (TkinScaled - E1)*W; 559 << 619 560 numOfTR = (G4int)G4Poisson((((*(*fEnergy << 620 // G4cout<<iTkin<<" mean TR number = "<<(((*(*fEnergyDistrTable)(iPlace))(0)+ 561 (*(*fAngleD << 621 // (*(*fAngleDistrTable)(iPlace))(0))*W1 + 562 W1 + << 622 // ((*(*fEnergyDistrTable)(iPlace + 1))(0)+ 563 ((*(*fEnergy << 623 // (*(*fAngleDistrTable)(iPlace + 1))(0))*W2) 564 (*(*fAngleD << 624 // *chargeSq*0.5<<G4endl; 565 W2) * << 625 566 chargeSq * 0. << 626 numOfTR = G4Poisson((((*(*fEnergyDistrTable)(iPlace))(0)+ >> 627 (*(*fAngleDistrTable)(iPlace))(0))*W1 + >> 628 ((*(*fEnergyDistrTable)(iPlace + 1))(0)+ >> 629 (*(*fAngleDistrTable)(iPlace + 1))(0))*W2) >> 630 *chargeSq*0.5 ); 567 if(numOfTR == 0) 631 if(numOfTR == 0) 568 { 632 { 569 return G4VDiscreteProcess::PostStepDoI 633 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep); 570 } 634 } 571 else 635 else 572 { 636 { >> 637 // G4cout<<"Number of X-ray TR photons = "<<numOfTR<<G4endl; >> 638 573 aParticleChange.SetNumberOfSecondaries 639 aParticleChange.SetNumberOfSecondaries(numOfTR); 574 for(iTR = 0; iTR < numOfTR; ++iTR) << 640 for(iTR=0;iTR<numOfTR;iTR++) 575 { 641 { 576 energyPos = ((*(*fEnergyDistrTable)( << 642 energyPos = ((*(*fEnergyDistrTable)(iPlace))(0)*W1+ 577 (*(*fEnergyDistrTable)( << 643 (*(*fEnergyDistrTable)(iPlace + 1))(0)*W2)*G4UniformRand(); 578 G4UniformRand(); << 644 for(iTransfer=0;iTransfer<fBinTR-1;iTransfer++) 579 for(iTransfer = 0; iTransfer < fBinT << 645 { 580 { << 646 if(energyPos >= ((*(*fEnergyDistrTable)(iPlace))(iTransfer)*W1+ 581 if(energyPos >= << 647 (*(*fEnergyDistrTable)(iPlace + 1))(iTransfer)*W2)) break; 582 ((*(*fEnergyDistrTable)(iPlace) << 648 } 583 (*(*fEnergyDistrTable)(iPlace << 649 energyTR = ((*fEnergyDistrTable)(iPlace)->GetLowEdgeEnergy(iTransfer))*W1+ 584 break; << 650 ((*fEnergyDistrTable)(iPlace + 1)->GetLowEdgeEnergy(iTransfer))*W2; 585 } << 651 586 energyTR = << 652 // G4cout<<"energyTR = "<<energyTR/keV<<"keV"<<G4endl; 587 ((*fEnergyDistrTable)(iPlace)->Get << 588 ((*fEnergyDistrTable)(iPlace + 1)- << 589 W2; << 590 653 591 kinEnergy -= energyTR; 654 kinEnergy -= energyTR; 592 aParticleChange.ProposeEnergy(kinEne 655 aParticleChange.ProposeEnergy(kinEnergy); 593 656 594 anglePos = ((*(*fAngleDistrTable)(iP << 657 anglePos = ((*(*fAngleDistrTable)(iPlace))(0)*W1+ 595 (*(*fAngleDistrTable)(iP << 658 (*(*fAngleDistrTable)(iPlace + 1))(0)*W2)*G4UniformRand(); 596 G4UniformRand(); << 659 for(iTransfer=0;iTransfer<fBinTR-1;iTransfer++) 597 for(iTransfer = 0; iTransfer < fBinT << 660 { 598 { << 661 if(anglePos > ((*(*fAngleDistrTable)(iPlace))(iTransfer)*W1+ 599 if(anglePos > ((*(*fAngleDistrTabl << 662 (*(*fAngleDistrTable)(iPlace + 1))(iTransfer)*W2)) break; 600 (*(*fAngleDistrTabl << 663 } 601 break; << 664 theta = std::sqrt(((*fAngleDistrTable)(iPlace)-> 602 } << 665 GetLowEdgeEnergy(iTransfer-1))*W1+ 603 theta = std::sqrt( << 666 ((*fAngleDistrTable)(iPlace + 1)-> 604 ((*fAngleDistrTable)(iPlace)->GetL << 667 GetLowEdgeEnergy(iTransfer-1))*W2); 605 W1 + << 668 606 ((*fAngleDistrTable)(iPlace + 1)-> << 669 // G4cout<<iTransfer<<" : theta = "<<theta<<G4endl; 607 W2); << 670 608 << 671 phi = twopi*G4UniformRand(); 609 phi = twopi * G4UniformRand(); << 672 dirX = std::sin(theta)*std::cos(phi) ; 610 dirX = std::sin(theta) * std::cos(ph << 673 dirY = std::sin(theta)*std::sin(phi) ; 611 dirY = std::sin(theta) * std::sin(ph << 674 dirZ = std::cos(theta) ; 612 dirZ = std::cos(theta); << 675 G4ThreeVector directionTR(dirX,dirY,dirZ); 613 G4ThreeVector directionTR(dirX, dirY << 614 directionTR.rotateUz(particleDir); 676 directionTR.rotateUz(particleDir); 615 auto aPhotonTR = << 677 G4DynamicParticle* aPhotonTR = new G4DynamicParticle(G4Gamma::Gamma(), 616 new G4DynamicParticle(G4Gamma::Gam << 678 directionTR, 617 << 679 energyTR ); 618 // Create the G4Track << 680 aParticleChange.AddSecondary(aPhotonTR); 619 G4Track* aSecondaryTrack = new G4Track(aPh << 620 aSecondaryTrack->SetTouchableHandle(aStep. << 621 aSecondaryTrack->SetParentID(aTrack.GetTra << 622 aSecondaryTrack->SetCreatorModelID(secID); << 623 aParticleChange.AddSecondary(aSecondaryTra << 624 } 681 } 625 } 682 } 626 } 683 } 627 } 684 } 628 return &aParticleChange; 685 return &aParticleChange; 629 } 686 } 630 687 631 ////////////////////////////////////////////// 688 //////////////////////////////////////////////////////////////////////////// >> 689 // 632 // Test function for checking of PostStepDoIt 690 // Test function for checking of PostStepDoIt random preparation of TR photon 633 // energy 691 // energy 634 G4double G4ForwardXrayTR::GetEnergyTR(G4int iM << 692 // >> 693 >> 694 G4double >> 695 G4ForwardXrayTR::GetEnergyTR(G4int iMat, G4int jMat, G4int iTkin) const 635 { 696 { 636 G4int iPlace, numOfTR, iTR, iTransfer; << 697 G4int iPlace, numOfTR, iTR, iTransfer; 637 G4double energyTR = 0.0; // return this val << 698 G4double energyTR = 0.0; // return this value for no TR photons 638 G4double energyPos; << 699 G4double energyPos ; 639 G4double W1, W2; << 700 G4double W1, W2; 640 << 701 641 const G4ProductionCutsTable* theCoupleTable << 702 const G4ProductionCutsTable* theCoupleTable= 642 G4ProductionCutsTable::GetProductionCutsTa << 703 G4ProductionCutsTable::GetProductionCutsTable(); 643 G4int numOfCouples = (G4int)theCoupleTable-> << 704 G4int numOfCouples = theCoupleTable->GetTableSize(); 644 705 645 // The case of equal or approximate (in term 706 // The case of equal or approximate (in terms of plasma energy) materials 646 // No TR photons ?! 707 // No TR photons ?! 647 708 648 const G4MaterialCutsCouple* iCouple = << 709 const G4MaterialCutsCouple* iCouple = theCoupleTable->GetMaterialCutsCouple(iMat); 649 theCoupleTable->GetMaterialCutsCouple(iMat << 710 const G4MaterialCutsCouple* jCouple = theCoupleTable->GetMaterialCutsCouple(jMat); 650 const G4MaterialCutsCouple* jCouple = << 651 theCoupleTable->GetMaterialCutsCouple(jMat << 652 const G4Material* iMaterial = iCouple->GetMa 711 const G4Material* iMaterial = iCouple->GetMaterial(); 653 const G4Material* jMaterial = jCouple->GetMa 712 const G4Material* jMaterial = jCouple->GetMaterial(); 654 713 655 if(iMat == jMat << 714 if ( iMat == jMat 656 715 657 || iMaterial->GetState() == jMaterial->Ge << 716 || iMaterial->GetState() == jMaterial->GetState() 658 717 659 || (iMaterial->GetState() == kStateSolid << 718 ||(iMaterial->GetState() == kStateSolid && jMaterial->GetState() == kStateLiquid ) 660 jMaterial->GetState() == kStateLiquid << 661 719 662 || (iMaterial->GetState() == kStateLiquid << 720 ||(iMaterial->GetState() == kStateLiquid && jMaterial->GetState() == kStateSolid ) ) 663 jMaterial->GetState() == kStateSolid) << 664 721 665 { 722 { 666 return energyTR; 723 return energyTR; 667 } 724 } 668 725 669 if(jMat < iMat) 726 if(jMat < iMat) 670 { 727 { 671 iPlace = (iMat * (numOfCouples - 1) + jMat << 728 iPlace = (iMat*(numOfCouples - 1) + jMat)*fTotBin + iTkin - 1; 672 } 729 } 673 else 730 else 674 { 731 { 675 iPlace = (iMat * (numOfCouples - 1) + jMat << 732 iPlace = (iMat*(numOfCouples - 1) + jMat - 1)*fTotBin + iTkin - 1; 676 } 733 } 677 G4PhysicsVector* energyVector1 = (*fEnergyDi << 734 G4PhysicsVector* energyVector1 = (*fEnergyDistrTable)(iPlace) ; 678 G4PhysicsVector* energyVector2 = (*fEnergyDi << 735 G4PhysicsVector* energyVector2 = (*fEnergyDistrTable)(iPlace + 1); 679 736 680 if(iTkin == fTotBin) // TR plato, try from << 737 if(iTkin == fTotBin) // TR plato, try from left 681 { 738 { 682 numOfTR = (G4int)G4Poisson((*energyVector1 << 739 numOfTR = G4Poisson( (*energyVector1)(0) ); 683 if(numOfTR == 0) 740 if(numOfTR == 0) 684 { 741 { 685 return energyTR; 742 return energyTR; 686 } 743 } 687 else 744 else 688 { 745 { 689 for(iTR = 0; iTR < numOfTR; ++iTR) << 746 for(iTR=0;iTR<numOfTR;iTR++) 690 { 747 { 691 energyPos = (*energyVector1)(0) * G4Un << 748 energyPos = (*energyVector1)(0)*G4UniformRand(); 692 for(iTransfer = 0; iTransfer < fBinTR << 749 for(iTransfer=0;iTransfer<fBinTR-1;iTransfer++) 693 { << 750 { 694 if(energyPos >= (*energyVector1)(iTr << 751 if(energyPos >= (*energyVector1)(iTransfer)) break; 695 break; << 752 } 696 } << 697 energyTR += energyVector1->GetLowEdgeE 753 energyTR += energyVector1->GetLowEdgeEnergy(iTransfer); 698 } 754 } 699 } 755 } 700 } 756 } 701 else 757 else 702 { 758 { 703 if(iTkin == 0) // Tkin is too small, negl << 759 if(iTkin == 0) // Tkin is too small, neglect of TR photon generation 704 { 760 { 705 return energyTR; 761 return energyTR; 706 } 762 } 707 else // general case: Tkin between two ve << 763 else // general case: Tkin between two vectors of the material 708 { // use trivial mean half/half << 764 { // use trivial mean half/half 709 W1 = 0.5; << 765 W1 = 0.5; 710 W2 = 0.5; << 766 W2 = 0.5; 711 numOfTR = (G4int)G4Poisson((*energyVecto << 767 numOfTR = G4Poisson( (*energyVector1)(0)*W1 + >> 768 (*energyVector2)(0)*W2 ); 712 if(numOfTR == 0) 769 if(numOfTR == 0) 713 { 770 { 714 return energyTR; 771 return energyTR; 715 } 772 } 716 else 773 else 717 { 774 { 718 G4cout << "It is still OK in GetEnergy << 775 G4cout<<"It is still OK in GetEnergyTR(int,int,int)"<<G4endl; 719 for(iTR = 0; iTR < numOfTR; ++iTR) << 776 for(iTR=0;iTR<numOfTR;iTR++) 720 { 777 { 721 energyPos = ((*energyVector1)(0) * W << 778 energyPos = ((*energyVector1)(0)*W1+ 722 G4UniformRand(); << 779 (*energyVector2)(0)*W2)*G4UniformRand(); 723 for(iTransfer = 0; iTransfer < fBinT << 780 for(iTransfer=0;iTransfer<fBinTR-1;iTransfer++) 724 { << 781 { 725 if(energyPos >= ((*energyVector1)( << 782 if(energyPos >= ((*energyVector1)(iTransfer)*W1+ 726 (*energyVector2)( << 783 (*energyVector2)(iTransfer)*W2)) break; 727 break; << 784 } 728 } << 785 energyTR += (energyVector1->GetLowEdgeEnergy(iTransfer))*W1+ 729 energyTR += (energyVector1->GetLowEd << 786 (energyVector2->GetLowEdgeEnergy(iTransfer))*W2; 730 (energyVector2->GetLowEd << 787 731 } 788 } 732 } 789 } 733 } 790 } 734 } 791 } 735 792 736 return energyTR; 793 return energyTR; 737 } 794 } 738 795 739 ////////////////////////////////////////////// 796 //////////////////////////////////////////////////////////////////////////// >> 797 // 740 // Test function for checking of PostStepDoIt 798 // Test function for checking of PostStepDoIt random preparation of TR photon 741 // theta angle relative to particle direction 799 // theta angle relative to particle direction 742 G4double G4ForwardXrayTR::GetThetaTR(G4int, G4 << 800 // >> 801 >> 802 G4double >> 803 G4ForwardXrayTR::GetThetaTR(G4int, G4int, G4int) const >> 804 { >> 805 return 0.0; >> 806 } 743 807 744 G4int G4ForwardXrayTR::GetSympsonNumber() { re << 808 G4int G4ForwardXrayTR::GetSympsonNumber() >> 809 { >> 810 return fSympsonNumber; >> 811 } 745 812 746 G4int G4ForwardXrayTR::GetBinTR() { return fBi << 813 G4int G4ForwardXrayTR::GetBinTR() >> 814 { >> 815 return fBinTR; >> 816 } 747 817 748 G4double G4ForwardXrayTR::GetMinProtonTkin() { << 818 G4double G4ForwardXrayTR::GetMinProtonTkin() >> 819 { >> 820 return fMinProtonTkin; >> 821 } 749 822 750 G4double G4ForwardXrayTR::GetMaxProtonTkin() { << 823 G4double G4ForwardXrayTR::GetMaxProtonTkin() >> 824 { >> 825 return fMaxProtonTkin; >> 826 } 751 827 752 G4int G4ForwardXrayTR::GetTotBin() { return fT << 828 G4int G4ForwardXrayTR::GetTotBin() >> 829 { >> 830 return fTotBin; >> 831 } >> 832 >> 833 // end of G4ForwardXrayTR implementation file >> 834 // >> 835 /////////////////////////////////////////////////////////////////////////// 753 836