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 // 26 // 27 // ------------------------------------------- 27 // -------------------------------------------------------------- 28 // 28 // 29 // File name: G4PenelopeBremsstrahlungAngu 29 // File name: G4PenelopeBremsstrahlungAngular 30 // 30 // 31 // Author: Luciano Pandola 31 // Author: Luciano Pandola 32 // 32 // 33 // Creation date: November 2010 33 // Creation date: November 2010 34 // 34 // 35 // History: 35 // History: 36 // ----------- 36 // ----------- 37 // 23 Nov 2010 L. Pandola 1st implement 37 // 23 Nov 2010 L. Pandola 1st implementation 38 // 24 May 2011 L. Pandola Renamed (make 38 // 24 May 2011 L. Pandola Renamed (make v2008 as default Penelope) 39 // 13 Mar 2012 L. Pandola Made a derive 39 // 13 Mar 2012 L. Pandola Made a derived class of G4VEmAngularDistribution 40 // and update th 40 // and update the interface accordingly 41 // 18 Jul 2012 L. Pandola Migrated to t 41 // 18 Jul 2012 L. Pandola Migrated to the new basic interface of G4VEmAngularDistribution 42 // Now returns a 42 // Now returns a G4ThreeVector and takes care of the rotation 43 // 03 Oct 2013 L. Pandola Migrated to M 43 // 03 Oct 2013 L. Pandola Migrated to MT: only the master model handles tables 44 // 17 Oct 2013 L. Pandola Partially rev 44 // 17 Oct 2013 L. Pandola Partially revert MT migration. The angular generator is kept as 45 // thread-local 45 // thread-local, and each worker has full access to it. 46 // 46 // 47 //-------------------------------------------- 47 //---------------------------------------------------------------- 48 48 49 #include "G4PenelopeBremsstrahlungAngular.hh" 49 #include "G4PenelopeBremsstrahlungAngular.hh" 50 50 51 #include "globals.hh" 51 #include "globals.hh" 52 #include "G4PhysicalConstants.hh" 52 #include "G4PhysicalConstants.hh" 53 #include "G4SystemOfUnits.hh" 53 #include "G4SystemOfUnits.hh" 54 #include "G4PhysicsFreeVector.hh" 54 #include "G4PhysicsFreeVector.hh" 55 #include "G4PhysicsTable.hh" 55 #include "G4PhysicsTable.hh" 56 #include "G4Material.hh" 56 #include "G4Material.hh" 57 #include "Randomize.hh" 57 #include "Randomize.hh" 58 #include "G4Exp.hh" 58 #include "G4Exp.hh" 59 59 60 G4PenelopeBremsstrahlungAngular::G4PenelopeBre 60 G4PenelopeBremsstrahlungAngular::G4PenelopeBremsstrahlungAngular() : 61 G4VEmAngularDistribution("Penelope"), fEffec << 61 G4VEmAngularDistribution("Penelope"), theEffectiveZSq(0), 62 fLorentzTables1(nullptr),fLorentzTables2(nul << 62 theLorentzTables1(0),theLorentzTables2(0) 63 63 64 { 64 { 65 fDataRead = false; << 65 dataRead = false; 66 fVerbosityLevel = 0; << 66 verbosityLevel = 0; 67 } 67 } 68 68 69 //....oooOO0OOooo........oooOO0OOooo........oo 69 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 70 70 71 G4PenelopeBremsstrahlungAngular::~G4PenelopeBr 71 G4PenelopeBremsstrahlungAngular::~G4PenelopeBremsstrahlungAngular() 72 { 72 { 73 ClearTables(); 73 ClearTables(); 74 } 74 } 75 75 76 //....oooOO0OOooo........oooOO0OOooo........oo 76 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 77 77 78 void G4PenelopeBremsstrahlungAngular::Initiali 78 void G4PenelopeBremsstrahlungAngular::Initialize() 79 { 79 { 80 ClearTables(); 80 ClearTables(); 81 } 81 } 82 82 83 //....oooOO0OOooo........oooOO0OOooo........oo 83 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 84 84 85 void G4PenelopeBremsstrahlungAngular::ClearTab 85 void G4PenelopeBremsstrahlungAngular::ClearTables() 86 { 86 { 87 if (fLorentzTables1) << 87 if (theLorentzTables1) 88 { 88 { 89 for (auto j=fLorentzTables1->begin(); j << 89 for (auto j = theLorentzTables1->begin(); j != theLorentzTables1->end(); j++) 90 { 90 { 91 G4PhysicsTable* tab = j->second; 91 G4PhysicsTable* tab = j->second; 92 tab->clearAndDestroy(); << 92 //tab->clearAndDestroy(); 93 delete tab; 93 delete tab; 94 } 94 } 95 fLorentzTables1->clear(); << 95 delete theLorentzTables1; 96 delete fLorentzTables1; << 96 theLorentzTables1 = nullptr; 97 fLorentzTables1 = nullptr; << 98 } 97 } 99 98 100 if (fLorentzTables2) << 99 if (theLorentzTables2) 101 { 100 { 102 for (auto j=fLorentzTables2->begin(); j << 101 for (auto j=theLorentzTables2->begin(); j != theLorentzTables2->end(); j++) 103 { 102 { 104 G4PhysicsTable* tab = j->second; 103 G4PhysicsTable* tab = j->second; 105 tab->clearAndDestroy(); << 104 //tab->clearAndDestroy(); 106 delete tab; 105 delete tab; 107 } 106 } 108 fLorentzTables2->clear(); << 107 delete theLorentzTables2; 109 delete fLorentzTables2; << 108 theLorentzTables2 = nullptr; 110 fLorentzTables2 = nullptr; << 111 } 109 } 112 if (fEffectiveZSq) << 110 if (theEffectiveZSq) 113 { 111 { 114 delete fEffectiveZSq; << 112 delete theEffectiveZSq; 115 fEffectiveZSq = nullptr; << 113 theEffectiveZSq = nullptr; 116 } 114 } 117 } 115 } 118 116 119 //....oooOO0OOooo........oooOO0OOooo........oo 117 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 120 118 121 void G4PenelopeBremsstrahlungAngular::ReadData 119 void G4PenelopeBremsstrahlungAngular::ReadDataFile() 122 { 120 { 123 //Read information from DataBase file 121 //Read information from DataBase file 124 const char* path = G4FindDataDir("G4LEDATA") << 122 char* path = std::getenv("G4LEDATA"); 125 if (!path) 123 if (!path) 126 { 124 { 127 G4String excep = 125 G4String excep = 128 "G4PenelopeBremsstrahlungAngular - G4LEDATA 126 "G4PenelopeBremsstrahlungAngular - G4LEDATA environment variable not set!"; 129 G4Exception("G4PenelopeBremsstrahlungAng 127 G4Exception("G4PenelopeBremsstrahlungAngular::ReadDataFile()", 130 "em0006",FatalException,excep); 128 "em0006",FatalException,excep); 131 return; 129 return; 132 } 130 } 133 G4String pathString(path); 131 G4String pathString(path); 134 G4String pathFile = pathString + "/penelope/ 132 G4String pathFile = pathString + "/penelope/bremsstrahlung/pdbrang.p08"; 135 std::ifstream file(pathFile); 133 std::ifstream file(pathFile); 136 134 137 if (!file.is_open()) 135 if (!file.is_open()) 138 { 136 { 139 G4String excep = "G4PenelopeBremsstrahlu 137 G4String excep = "G4PenelopeBremsstrahlungAngular - data file " + pathFile + " not found!"; 140 G4Exception("G4PenelopeBremsstrahlungAng 138 G4Exception("G4PenelopeBremsstrahlungAngular::ReadDataFile()", 141 "em0003",FatalException,excep); 139 "em0003",FatalException,excep); 142 return; 140 return; 143 } 141 } 144 G4int i=0,j=0,k=0; // i=index for Z, j=index 142 G4int i=0,j=0,k=0; // i=index for Z, j=index for E, k=index for K 145 143 146 for (k=0;k<fNumberofKPoints;k++) << 144 for (k=0;k<NumberofKPoints;k++) 147 for (i=0;i<fNumberofZPoints;i++) << 145 for (i=0;i<NumberofZPoints;i++) 148 for (j=0;j<fNumberofEPoints;j++) << 146 for (j=0;j<NumberofEPoints;j++) 149 { 147 { 150 G4double a1,a2; 148 G4double a1,a2; 151 G4int ik1,iz1,ie1; 149 G4int ik1,iz1,ie1; 152 G4double zr,er,kr; 150 G4double zr,er,kr; 153 file >> iz1 >> ie1 >> ik1 >> zr >> er >> k 151 file >> iz1 >> ie1 >> ik1 >> zr >> er >> kr >> a1 >> a2; 154 //check the data are correct 152 //check the data are correct 155 if ((iz1-1 == i) && (ik1-1 == k) && (ie1-1 153 if ((iz1-1 == i) && (ik1-1 == k) && (ie1-1 == j)) 156 { 154 { 157 fQQ1[i][j][k]=a1; << 155 QQ1[i][j][k]=a1; 158 fQQ2[i][j][k]=a2; << 156 QQ2[i][j][k]=a2; 159 } 157 } 160 else 158 else 161 { 159 { 162 G4ExceptionDescription ed; 160 G4ExceptionDescription ed; 163 ed << "Corrupted data file " << pathFi 161 ed << "Corrupted data file " << pathFile << "?" << G4endl; 164 G4Exception("G4PenelopeBremsstrahlungA 162 G4Exception("G4PenelopeBremsstrahlungAngular::ReadDataFile()", 165 "em0005",FatalException,ed); 163 "em0005",FatalException,ed); 166 } 164 } 167 } 165 } 168 file.close(); 166 file.close(); 169 fDataRead = true; << 167 dataRead = true; 170 } 168 } 171 169 172 //....oooOO0OOooo........oooOO0OOooo........oo 170 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 173 171 174 void G4PenelopeBremsstrahlungAngular::PrepareT 172 void G4PenelopeBremsstrahlungAngular::PrepareTables(const G4Material* material,G4bool /*isMaster*/ ) 175 { 173 { 176 //Unused at the moment: the G4PenelopeBremss 174 //Unused at the moment: the G4PenelopeBremsstrahlungAngular is thread-local, so each worker 177 //builds its own version of the tables. 175 //builds its own version of the tables. 178 << 176 /* >> 177 if (!isMaster) >> 178 //Should not be here! >> 179 G4Exception("G4PenelopeBremsstrahlungAngular::PrepareTables()", >> 180 "em0100",FatalException,"Worker thread in this method"); >> 181 */ >> 182 179 //Check if data file has already been read 183 //Check if data file has already been read 180 if (!fDataRead) << 184 if (!dataRead) 181 { 185 { 182 ReadDataFile(); 186 ReadDataFile(); 183 if (!fDataRead) << 187 if (!dataRead) 184 G4Exception("G4PenelopeBremsstrahlungAngular 188 G4Exception("G4PenelopeBremsstrahlungAngular::PrepareInterpolationTables()", 185 "em2001",FatalException,"Unable to bui 189 "em2001",FatalException,"Unable to build interpolation table"); 186 } 190 } 187 191 188 if (!fLorentzTables1) << 192 if (!theLorentzTables1) 189 fLorentzTables1 = new std::map<G4double, << 193 theLorentzTables1 = new std::map<G4double,G4PhysicsTable*>; 190 if (!fLorentzTables2) << 194 if (!theLorentzTables2) 191 fLorentzTables2 = new std::map<G4double,G4 << 195 theLorentzTables2 = new std::map<G4double,G4PhysicsTable*>; 192 196 193 G4double Zmat = CalculateEffectiveZ(material 197 G4double Zmat = CalculateEffectiveZ(material); 194 198 195 const G4int reducedEnergyGrid=21; 199 const G4int reducedEnergyGrid=21; 196 //Support arrays. 200 //Support arrays. 197 G4double betas[fNumberofEPoints]; //betas fo << 201 G4double betas[NumberofEPoints]; //betas for interpolation 198 //tables for interpolation 202 //tables for interpolation 199 G4double Q1[fNumberofEPoints][fNumberofKPoin << 203 G4double Q1[NumberofEPoints][NumberofKPoints]; 200 G4double Q2[fNumberofEPoints][fNumberofKPoin << 204 G4double Q2[NumberofEPoints][NumberofKPoints]; 201 //expanded tables for interpolation 205 //expanded tables for interpolation 202 G4double Q1E[fNumberofEPoints][reducedEnergy << 206 G4double Q1E[NumberofEPoints][reducedEnergyGrid]; 203 G4double Q2E[fNumberofEPoints][reducedEnergy << 207 G4double Q2E[NumberofEPoints][reducedEnergyGrid]; 204 G4double pZ[fNumberofZPoints] = {2.0,8.0,13. << 208 G4double pZ[NumberofZPoints] = {2.0,8.0,13.0,47.0,79.0,92.0}; 205 209 206 G4int i=0,j=0,k=0; // i=index for Z, j=index 210 G4int i=0,j=0,k=0; // i=index for Z, j=index for E, k=index for K 207 //Interpolation in Z 211 //Interpolation in Z 208 for (i=0;i<fNumberofEPoints;i++) << 212 for (i=0;i<NumberofEPoints;i++) 209 { 213 { 210 for (j=0;j<fNumberofKPoints;j++) << 214 for (j=0;j<NumberofKPoints;j++) 211 { 215 { 212 G4PhysicsFreeVector* QQ1vector = << 216 G4PhysicsFreeVector* QQ1vector = new G4PhysicsFreeVector(NumberofZPoints); 213 new G4PhysicsFreeVector(fNumberofZPoints << 217 G4PhysicsFreeVector* QQ2vector = new G4PhysicsFreeVector(NumberofZPoints); 214 G4PhysicsFreeVector* QQ2vector = << 215 new G4PhysicsFreeVector(fNumberofZPoints << 216 218 217 //fill vectors 219 //fill vectors 218 for (k=0;k<fNumberofZPoints;k++) << 220 for (k=0;k<NumberofZPoints;k++) 219 { 221 { 220 QQ1vector->PutValues(k,pZ[k],G4Log(fQQ << 222 QQ1vector->PutValue(k,pZ[k],std::log(QQ1[k][i][j])); 221 QQ2vector->PutValues(k,pZ[k],fQQ2[k][i << 223 QQ2vector->PutValue(k,pZ[k],QQ2[k][i][j]); 222 } 224 } 223 //Filled: now calculate derivatives << 225 224 QQ1vector->FillSecondDerivatives(); << 226 QQ1vector->SetSpline(true); 225 QQ2vector->FillSecondDerivatives(); << 227 QQ2vector->SetSpline(true); 226 228 227 Q1[i][j]= G4Exp(QQ1vector->Value(Zmat)); 229 Q1[i][j]= G4Exp(QQ1vector->Value(Zmat)); 228 Q2[i][j]=QQ2vector->Value(Zmat); 230 Q2[i][j]=QQ2vector->Value(Zmat); 229 delete QQ1vector; 231 delete QQ1vector; 230 delete QQ2vector; 232 delete QQ2vector; 231 } 233 } 232 } 234 } 233 G4double pE[fNumberofEPoints] = {1.0e-03*MeV << 235 G4double pE[NumberofEPoints] = {1.0e-03*MeV,5.0e-03*MeV,1.0e-02*MeV,5.0e-02*MeV, 234 1.0e-01*MeV,5.0e-01*MeV}; 236 1.0e-01*MeV,5.0e-01*MeV}; 235 G4double pK[fNumberofKPoints] = {0.0,0.6,0.8 << 237 G4double pK[NumberofKPoints] = {0.0,0.6,0.8,0.95}; 236 G4double ppK[reducedEnergyGrid]; 238 G4double ppK[reducedEnergyGrid]; 237 239 238 for(i=0;i<reducedEnergyGrid;i++) 240 for(i=0;i<reducedEnergyGrid;i++) 239 ppK[i]=((G4double) i) * 0.05; 241 ppK[i]=((G4double) i) * 0.05; 240 242 241 243 242 for(i=0;i<fNumberofEPoints;i++) << 244 for(i=0;i<NumberofEPoints;i++) 243 betas[i]=std::sqrt(pE[i]*(pE[i]+2*electron 245 betas[i]=std::sqrt(pE[i]*(pE[i]+2*electron_mass_c2))/(pE[i]+electron_mass_c2); 244 246 245 247 246 for (i=0;i<fNumberofEPoints;i++) << 248 for (i=0;i<NumberofEPoints;i++) 247 { 249 { 248 for (j=0;j<fNumberofKPoints;j++) << 250 for (j=0;j<NumberofKPoints;j++) 249 Q1[i][j]=Q1[i][j]/Zmat; 251 Q1[i][j]=Q1[i][j]/Zmat; 250 } 252 } 251 253 252 //Expanded table of distribution parameters 254 //Expanded table of distribution parameters 253 for (i=0;i<fNumberofEPoints;i++) << 255 for (i=0;i<NumberofEPoints;i++) 254 { 256 { 255 G4PhysicsFreeVector* Q1vector = new G4Ph << 257 G4PhysicsFreeVector* Q1vector = new G4PhysicsFreeVector(NumberofKPoints); 256 G4PhysicsFreeVector* Q2vector = new G4Ph << 258 G4PhysicsFreeVector* Q2vector = new G4PhysicsFreeVector(NumberofKPoints); 257 259 258 for (j=0;j<fNumberofKPoints;j++) << 260 for (j=0;j<NumberofKPoints;j++) 259 { 261 { 260 Q1vector->PutValues(j,pK[j],G4Log(Q1[i][j] << 262 Q1vector->PutValue(j,pK[j],std::log(Q1[i][j])); //logarithmic 261 Q2vector->PutValues(j,pK[j],Q2[i][j]); << 263 Q2vector->PutValue(j,pK[j],Q2[i][j]); 262 } 264 } 263 265 264 for (j=0;j<reducedEnergyGrid;j++) 266 for (j=0;j<reducedEnergyGrid;j++) 265 { 267 { 266 Q1E[i][j]=Q1vector->Value(ppK[j]); 268 Q1E[i][j]=Q1vector->Value(ppK[j]); 267 Q2E[i][j]=Q2vector->Value(ppK[j]); 269 Q2E[i][j]=Q2vector->Value(ppK[j]); 268 } 270 } 269 delete Q1vector; 271 delete Q1vector; 270 delete Q2vector; 272 delete Q2vector; 271 } 273 } 272 // 274 // 273 //TABLES to be stored 275 //TABLES to be stored 274 // 276 // 275 G4PhysicsTable* theTable1 = new G4PhysicsTab 277 G4PhysicsTable* theTable1 = new G4PhysicsTable(); 276 G4PhysicsTable* theTable2 = new G4PhysicsTab 278 G4PhysicsTable* theTable2 = new G4PhysicsTable(); 277 // the table will contain reducedEnergyGrid 279 // the table will contain reducedEnergyGrid G4PhysicsFreeVectors with different 278 // values of k, 280 // values of k, 279 // Each of the G4PhysicsFreeVectors has a pr 281 // Each of the G4PhysicsFreeVectors has a profile of 280 // y vs. E 282 // y vs. E 281 // 283 // 282 //reserve space of the vectors. 284 //reserve space of the vectors. 283 for (j=0;j<reducedEnergyGrid;j++) 285 for (j=0;j<reducedEnergyGrid;j++) 284 { 286 { 285 G4PhysicsFreeVector* thevec = new G4Phys << 287 G4PhysicsFreeVector* thevec = new G4PhysicsFreeVector(NumberofEPoints); 286 theTable1->push_back(thevec); 288 theTable1->push_back(thevec); 287 G4PhysicsFreeVector* thevec2 = new G4Phy << 289 G4PhysicsFreeVector* thevec2 = new G4PhysicsFreeVector(NumberofEPoints); 288 theTable2->push_back(thevec2); 290 theTable2->push_back(thevec2); 289 } 291 } 290 292 291 for (j=0;j<reducedEnergyGrid;j++) 293 for (j=0;j<reducedEnergyGrid;j++) 292 { 294 { 293 G4PhysicsFreeVector* thevec = (G4Physics 295 G4PhysicsFreeVector* thevec = (G4PhysicsFreeVector*) (*theTable1)[j]; 294 G4PhysicsFreeVector* thevec2 = (G4Physic 296 G4PhysicsFreeVector* thevec2 = (G4PhysicsFreeVector*) (*theTable2)[j]; 295 for (i=0;i<fNumberofEPoints;i++) << 297 for (i=0;i<NumberofEPoints;i++) 296 { 298 { 297 thevec->PutValues(i,betas[i],Q1E[i][j]); << 299 thevec->PutValue(i,betas[i],Q1E[i][j]); 298 thevec2->PutValues(i,betas[i],Q2E[i][j]); << 300 thevec2->PutValue(i,betas[i],Q2E[i][j]); 299 } 301 } 300 //Vectors are filled: calculate derivati << 302 thevec->SetSpline(true); 301 thevec->FillSecondDerivatives(); << 303 thevec2->SetSpline(true); 302 thevec2->FillSecondDerivatives(); << 303 } 304 } 304 305 305 if (fLorentzTables1 && fLorentzTables2) << 306 if (theLorentzTables1 && theLorentzTables2) 306 { 307 { 307 fLorentzTables1->insert(std::make_pair(Z << 308 theLorentzTables1->insert(std::make_pair(Zmat,theTable1)); 308 fLorentzTables2->insert(std::make_pair(Z << 309 theLorentzTables2->insert(std::make_pair(Zmat,theTable2)); 309 } 310 } 310 else 311 else 311 { 312 { 312 G4ExceptionDescription ed; 313 G4ExceptionDescription ed; 313 ed << "Unable to create tables of Lorent 314 ed << "Unable to create tables of Lorentz coefficients for " << G4endl; 314 ed << "<Z>= " << Zmat << " in G4Penelop 315 ed << "<Z>= " << Zmat << " in G4PenelopeBremsstrahlungAngular" << G4endl; 315 delete theTable1; 316 delete theTable1; 316 delete theTable2; 317 delete theTable2; 317 G4Exception("G4PenelopeBremsstrahlungAng 318 G4Exception("G4PenelopeBremsstrahlungAngular::PrepareInterpolationTables()", 318 "em2005",FatalException,ed); 319 "em2005",FatalException,ed); 319 } 320 } 320 321 321 return; 322 return; 322 } 323 } 323 324 324 //....oooOO0OOooo........oooOO0OOooo........oo 325 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 325 326 326 G4ThreeVector& G4PenelopeBremsstrahlungAngular 327 G4ThreeVector& G4PenelopeBremsstrahlungAngular::SampleDirection(const G4DynamicParticle* dp, 327 G4double eGamma, 328 G4double eGamma, 328 G4int, 329 G4int, 329 const G4Material* material) 330 const G4Material* material) 330 { 331 { 331 if (!material) 332 if (!material) 332 { 333 { 333 G4Exception("G4PenelopeBremsstrahlungAng 334 G4Exception("G4PenelopeBremsstrahlungAngular::SampleDirection()", 334 "em2040",FatalException,"The pointer to 335 "em2040",FatalException,"The pointer to G4Material* is nullptr"); 335 return fLocalDirection; 336 return fLocalDirection; 336 } 337 } 337 338 338 //Retrieve the effective Z 339 //Retrieve the effective Z 339 G4double Zmat = 0; 340 G4double Zmat = 0; 340 341 341 //The model might be initialized incorrectly << 342 if (!theEffectiveZSq) 342 //G4PenelopeBremsstrahungModel: make sure it << 343 if (!fEffectiveZSq) << 344 { 343 { 345 G4Exception("G4PenelopeBremsstrahlungAng 344 G4Exception("G4PenelopeBremsstrahlungAngular::SampleDirection()", 346 "em2040",JustWarning,"EffectiveZSq table << 345 "em2040",FatalException,"EffectiveZ table not available"); 347 PrepareTables(material,false); << 346 return fLocalDirection; 348 //return fLocalDirection; << 347 } 349 } << 350 348 351 //found in the table: return it 349 //found in the table: return it 352 if (fEffectiveZSq->count(material)) << 350 if (theEffectiveZSq->count(material)) 353 Zmat = fEffectiveZSq->find(material)->seco << 351 Zmat = theEffectiveZSq->find(material)->second; 354 else //this can happen in unit tests or when << 352 else 355 //models other than G4PenelopeBremsstrahun << 353 { 356 { << 357 G4Exception("G4PenelopeBremsstrahlungAng 354 G4Exception("G4PenelopeBremsstrahlungAngular::SampleDirection()", 358 "em2040",JustWarning,"Material not found << 355 "em2040",FatalException,"Material not found in the effectiveZ table"); 359 PrepareTables(material,false); << 356 return fLocalDirection; 360 Zmat = fEffectiveZSq->find(material)->se << 361 // return fLocalDirection; << 362 } 357 } 363 358 364 if (fVerbosityLevel > 0) << 359 if (verbosityLevel > 0) 365 { 360 { 366 G4cout << "Effective <Z> for material : 361 G4cout << "Effective <Z> for material : " << material->GetName() << 367 " = " << Zmat << G4endl; 362 " = " << Zmat << G4endl; 368 } 363 } 369 364 370 G4double ePrimary = dp->GetKineticEnergy(); 365 G4double ePrimary = dp->GetKineticEnergy(); 371 366 372 G4double beta = std::sqrt(ePrimary*(ePrimary 367 G4double beta = std::sqrt(ePrimary*(ePrimary+2*electron_mass_c2))/ 373 (ePrimary+electron_mass_c2); 368 (ePrimary+electron_mass_c2); 374 G4double cdt = 0; 369 G4double cdt = 0; 375 G4double sinTheta = 0; 370 G4double sinTheta = 0; 376 G4double phi = 0; 371 G4double phi = 0; 377 372 378 //Use a pure dipole distribution for energy 373 //Use a pure dipole distribution for energy above 500 keV 379 if (ePrimary > 500*keV) 374 if (ePrimary > 500*keV) 380 { 375 { 381 cdt = 2.0*G4UniformRand() - 1.0; 376 cdt = 2.0*G4UniformRand() - 1.0; 382 if (G4UniformRand() > 0.75) 377 if (G4UniformRand() > 0.75) 383 { 378 { 384 if (cdt<0) 379 if (cdt<0) 385 cdt = -1.0*std::pow(-cdt,1./3.); 380 cdt = -1.0*std::pow(-cdt,1./3.); 386 else 381 else 387 cdt = std::pow(cdt,1./3.); 382 cdt = std::pow(cdt,1./3.); 388 } 383 } 389 cdt = (cdt+beta)/(1.0+beta*cdt); 384 cdt = (cdt+beta)/(1.0+beta*cdt); 390 //Get primary kinematics 385 //Get primary kinematics 391 sinTheta = std::sqrt(1. - cdt*cdt); 386 sinTheta = std::sqrt(1. - cdt*cdt); 392 phi = twopi * G4UniformRand(); 387 phi = twopi * G4UniformRand(); 393 fLocalDirection.set(sinTheta* std::cos(p 388 fLocalDirection.set(sinTheta* std::cos(phi), 394 sinTheta* std::sin(phi), 389 sinTheta* std::sin(phi), 395 cdt); 390 cdt); 396 //rotate 391 //rotate 397 fLocalDirection.rotateUz(dp->GetMomentum 392 fLocalDirection.rotateUz(dp->GetMomentumDirection()); 398 //return 393 //return 399 return fLocalDirection; 394 return fLocalDirection; 400 } 395 } 401 396 402 if (!(fLorentzTables1->count(Zmat)) || !(fLo << 397 if (!(theLorentzTables1->count(Zmat)) || !(theLorentzTables2->count(Zmat))) 403 { 398 { 404 G4ExceptionDescription ed; 399 G4ExceptionDescription ed; 405 ed << "Unable to retrieve Lorentz tables 400 ed << "Unable to retrieve Lorentz tables for Z= " << Zmat << G4endl; 406 G4Exception("G4PenelopeBremsstrahlungAng 401 G4Exception("G4PenelopeBremsstrahlungAngular::SampleDirection()", 407 "em2006",FatalException,ed); 402 "em2006",FatalException,ed); 408 } 403 } 409 404 410 //retrieve actual tables 405 //retrieve actual tables 411 const G4PhysicsTable* theTable1 = fLorentzTa << 406 const G4PhysicsTable* theTable1 = theLorentzTables1->find(Zmat)->second; 412 const G4PhysicsTable* theTable2 = fLorentzTa << 407 const G4PhysicsTable* theTable2 = theLorentzTables2->find(Zmat)->second; 413 408 414 G4double RK=20.0*eGamma/ePrimary; 409 G4double RK=20.0*eGamma/ePrimary; 415 G4int ik=std::min((G4int) RK,19); 410 G4int ik=std::min((G4int) RK,19); 416 411 417 G4double P10=0,P11=0,P1=0; 412 G4double P10=0,P11=0,P1=0; 418 G4double P20=0,P21=0,P2=0; 413 G4double P20=0,P21=0,P2=0; 419 414 420 //First coefficient 415 //First coefficient 421 const G4PhysicsFreeVector* v1 = (G4PhysicsFr 416 const G4PhysicsFreeVector* v1 = (G4PhysicsFreeVector*) (*theTable1)[ik]; 422 const G4PhysicsFreeVector* v2 = (G4PhysicsFr 417 const G4PhysicsFreeVector* v2 = (G4PhysicsFreeVector*) (*theTable1)[ik+1]; 423 P10 = v1->Value(beta); 418 P10 = v1->Value(beta); 424 P11 = v2->Value(beta); 419 P11 = v2->Value(beta); 425 P1=P10+(RK-(G4double) ik)*(P11-P10); 420 P1=P10+(RK-(G4double) ik)*(P11-P10); 426 421 427 //Second coefficient 422 //Second coefficient 428 const G4PhysicsFreeVector* v3 = (G4PhysicsFr 423 const G4PhysicsFreeVector* v3 = (G4PhysicsFreeVector*) (*theTable2)[ik]; 429 const G4PhysicsFreeVector* v4 = (G4PhysicsFr 424 const G4PhysicsFreeVector* v4 = (G4PhysicsFreeVector*) (*theTable2)[ik+1]; 430 P20=v3->Value(beta); 425 P20=v3->Value(beta); 431 P21=v4->Value(beta); 426 P21=v4->Value(beta); 432 P2=P20+(RK-(G4double) ik)*(P21-P20); 427 P2=P20+(RK-(G4double) ik)*(P21-P20); 433 428 434 //Sampling from the Lorenz-trasformed dipole 429 //Sampling from the Lorenz-trasformed dipole distributions 435 P1=std::min(G4Exp(P1)/beta,1.0); 430 P1=std::min(G4Exp(P1)/beta,1.0); 436 G4double betap = std::min(std::max(beta*(1.0 431 G4double betap = std::min(std::max(beta*(1.0+P2/beta),0.0),0.9999); 437 432 438 G4double testf=0; 433 G4double testf=0; 439 434 440 if (G4UniformRand() < P1) 435 if (G4UniformRand() < P1) 441 { 436 { 442 do{ 437 do{ 443 cdt = 2.0*G4UniformRand()-1.0; 438 cdt = 2.0*G4UniformRand()-1.0; 444 testf=2.0*G4UniformRand()-(1.0+cdt*cdt); 439 testf=2.0*G4UniformRand()-(1.0+cdt*cdt); 445 }while(testf>0); 440 }while(testf>0); 446 } 441 } 447 else 442 else 448 { 443 { 449 do{ 444 do{ 450 cdt = 2.0*G4UniformRand()-1.0; 445 cdt = 2.0*G4UniformRand()-1.0; 451 testf=G4UniformRand()-(1.0-cdt*cdt); 446 testf=G4UniformRand()-(1.0-cdt*cdt); 452 }while(testf>0); 447 }while(testf>0); 453 } 448 } 454 cdt = (cdt+betap)/(1.0+betap*cdt); 449 cdt = (cdt+betap)/(1.0+betap*cdt); 455 450 456 //Get primary kinematics 451 //Get primary kinematics 457 sinTheta = std::sqrt(1. - cdt*cdt); 452 sinTheta = std::sqrt(1. - cdt*cdt); 458 phi = twopi * G4UniformRand(); 453 phi = twopi * G4UniformRand(); 459 fLocalDirection.set(sinTheta* std::cos(phi), 454 fLocalDirection.set(sinTheta* std::cos(phi), 460 sinTheta* std::sin(phi), 455 sinTheta* std::sin(phi), 461 cdt); 456 cdt); 462 //rotate 457 //rotate 463 fLocalDirection.rotateUz(dp->GetMomentumDire 458 fLocalDirection.rotateUz(dp->GetMomentumDirection()); 464 //return 459 //return 465 return fLocalDirection; 460 return fLocalDirection; 466 461 467 } 462 } 468 463 469 //....oooOO0OOooo........oooOO0OOooo........oo 464 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 470 465 >> 466 G4double G4PenelopeBremsstrahlungAngular::PolarAngle(const G4double , >> 467 const G4double , >> 468 const G4int ) >> 469 { >> 470 G4cout << "WARNING: G4PenelopeBremsstrahlungAngular() does NOT support PolarAngle()" << G4endl; >> 471 G4cout << "Please use the alternative interface SampleDirection()" << G4endl; >> 472 G4Exception("G4PenelopeBremsstrahlungAngular::PolarAngle()", >> 473 "em0005",FatalException,"Unsupported interface"); >> 474 return 0; >> 475 } >> 476 >> 477 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... >> 478 471 G4double G4PenelopeBremsstrahlungAngular::Calc 479 G4double G4PenelopeBremsstrahlungAngular::CalculateEffectiveZ(const G4Material* material) 472 { 480 { 473 if (!fEffectiveZSq) << 481 if (!theEffectiveZSq) 474 fEffectiveZSq = new std::map<const G4Mater << 482 theEffectiveZSq = new std::map<const G4Material*,G4double>; 475 483 476 //Already exists: return it 484 //Already exists: return it 477 if (fEffectiveZSq->count(material)) << 485 if (theEffectiveZSq->count(material)) 478 return fEffectiveZSq->find(material)->seco << 486 return theEffectiveZSq->find(material)->second; 479 487 480 //Helper for the calculation 488 //Helper for the calculation 481 std::vector<G4double> *StechiometricFactors 489 std::vector<G4double> *StechiometricFactors = new std::vector<G4double>; 482 G4int nElements = (G4int)material->GetNumber << 490 G4int nElements = material->GetNumberOfElements(); 483 const G4ElementVector* elementVector = mater 491 const G4ElementVector* elementVector = material->GetElementVector(); 484 const G4double* fractionVector = material->G 492 const G4double* fractionVector = material->GetFractionVector(); 485 for (G4int i=0;i<nElements;++i) << 493 for (G4int i=0;i<nElements;i++) 486 { 494 { 487 G4double fraction = fractionVector[i]; 495 G4double fraction = fractionVector[i]; 488 G4double atomicWeigth = (*elementVector) 496 G4double atomicWeigth = (*elementVector)[i]->GetA()/(g/mole); 489 StechiometricFactors->push_back(fraction 497 StechiometricFactors->push_back(fraction/atomicWeigth); 490 } 498 } 491 //Find max 499 //Find max 492 G4double MaxStechiometricFactor = 0.; 500 G4double MaxStechiometricFactor = 0.; 493 for (G4int i=0;i<nElements;++i) << 501 for (G4int i=0;i<nElements;i++) 494 { 502 { 495 if ((*StechiometricFactors)[i] > MaxStec 503 if ((*StechiometricFactors)[i] > MaxStechiometricFactor) 496 MaxStechiometricFactor = (*Stechiometr 504 MaxStechiometricFactor = (*StechiometricFactors)[i]; 497 } 505 } 498 //Normalize 506 //Normalize 499 for (G4int i=0;i<nElements;++i) << 507 for (G4int i=0;i<nElements;i++) 500 (*StechiometricFactors)[i] /= MaxStechiom 508 (*StechiometricFactors)[i] /= MaxStechiometricFactor; 501 509 502 G4double sumz2 = 0; 510 G4double sumz2 = 0; 503 G4double sums = 0; 511 G4double sums = 0; 504 for (G4int i=0;i<nElements;++i) << 512 for (G4int i=0;i<nElements;i++) 505 { 513 { 506 G4double Z = (*elementVector)[i]->GetZ() 514 G4double Z = (*elementVector)[i]->GetZ(); 507 sumz2 += (*StechiometricFactors)[i]*Z*Z; 515 sumz2 += (*StechiometricFactors)[i]*Z*Z; 508 sums += (*StechiometricFactors)[i]; 516 sums += (*StechiometricFactors)[i]; 509 } 517 } 510 delete StechiometricFactors; 518 delete StechiometricFactors; 511 519 512 G4double ZBR = std::sqrt(sumz2/sums); 520 G4double ZBR = std::sqrt(sumz2/sums); 513 fEffectiveZSq->insert(std::make_pair(materia << 521 theEffectiveZSq->insert(std::make_pair(material,ZBR)); 514 522 515 return ZBR; 523 return ZBR; 516 } 524 } 517 525