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