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 #include <fstream> >> 29 #include <iomanip> >> 30 28 #include "G4AdjointCSManager.hh" 31 #include "G4AdjointCSManager.hh" 29 32 >> 33 #include "G4PhysicalConstants.hh" >> 34 #include "G4SystemOfUnits.hh" 30 #include "G4AdjointCSMatrix.hh" 35 #include "G4AdjointCSMatrix.hh" 31 #include "G4AdjointElectron.hh" << 32 #include "G4AdjointGamma.hh" << 33 #include "G4AdjointInterpolator.hh" 36 #include "G4AdjointInterpolator.hh" 34 #include "G4AdjointProton.hh" << 37 #include "G4AdjointCSMatrix.hh" 35 #include "G4Electron.hh" << 38 #include "G4VEmAdjointModel.hh" 36 #include "G4Element.hh" << 37 #include "G4ElementTable.hh" 39 #include "G4ElementTable.hh" 38 #include "G4Gamma.hh" << 40 #include "G4Element.hh" 39 #include "G4ParticleDefinition.hh" 41 #include "G4ParticleDefinition.hh" 40 #include "G4PhysicalConstants.hh" << 42 #include "G4Element.hh" 41 #include "G4PhysicsLogVector.hh" << 42 #include "G4PhysicsTable.hh" << 43 #include "G4ProductionCutsTable.hh" << 44 #include "G4Proton.hh" << 45 #include "G4SystemOfUnits.hh" << 46 #include "G4VEmAdjointModel.hh" << 47 #include "G4VEmProcess.hh" 43 #include "G4VEmProcess.hh" 48 #include "G4VEnergyLossProcess.hh" 44 #include "G4VEnergyLossProcess.hh" >> 45 #include "G4PhysicsTable.hh" >> 46 #include "G4PhysicsLogVector.hh" >> 47 #include "G4PhysicsTableHelper.hh" >> 48 #include "G4Electron.hh" >> 49 #include "G4Gamma.hh" >> 50 #include "G4Proton.hh" >> 51 #include "G4AdjointElectron.hh" >> 52 #include "G4AdjointGamma.hh" >> 53 #include "G4AdjointProton.hh" >> 54 #include "G4ProductionCutsTable.hh" >> 55 #include "G4ProductionCutsTable.hh" 49 56 50 G4ThreadLocal G4AdjointCSManager* G4AdjointCSM << 57 G4ThreadLocal G4AdjointCSManager* G4AdjointCSManager::theInstance = nullptr; 51 << 52 constexpr G4double G4AdjointCSManager::fTmin; << 53 constexpr G4double G4AdjointCSManager::fTmax; << 54 constexpr G4int G4AdjointCSManager::fNbins; << 55 << 56 ////////////////////////////////////////////// 58 /////////////////////////////////////////////////////// >> 59 // 57 G4AdjointCSManager* G4AdjointCSManager::GetAdj 60 G4AdjointCSManager* G4AdjointCSManager::GetAdjointCSManager() 58 { << 61 { 59 if(fInstance == nullptr) << 62 if(theInstance == nullptr) { 60 { << 61 static G4ThreadLocalSingleton<G4AdjointCSM 63 static G4ThreadLocalSingleton<G4AdjointCSManager> inst; 62 fInstance = inst.Instance(); << 64 theInstance = inst.Instance(); 63 } 65 } 64 return fInstance; << 66 return theInstance; 65 } 67 } 66 68 67 ////////////////////////////////////////////// 69 /////////////////////////////////////////////////////// >> 70 // 68 G4AdjointCSManager::G4AdjointCSManager() 71 G4AdjointCSManager::G4AdjointCSManager() 69 { << 72 { CrossSectionMatrixesAreBuilt=false; >> 73 TotalSigmaTableAreBuilt=false; >> 74 theTotalForwardSigmaTableVector.clear(); >> 75 theTotalAdjointSigmaTableVector.clear(); >> 76 listOfForwardEmProcess.clear(); >> 77 listOfForwardEnergyLossProcess.clear(); >> 78 theListOfAdjointParticlesInAction.clear(); >> 79 EminForFwdSigmaTables.clear(); >> 80 EminForAdjSigmaTables.clear(); >> 81 EkinofFwdSigmaMax.clear(); >> 82 EkinofAdjSigmaMax.clear(); >> 83 listSigmaTableForAdjointModelScatProjToProj.clear(); >> 84 listSigmaTableForAdjointModelProdToProj.clear(); >> 85 Tmin=0.1*keV; >> 86 Tmax=100.*TeV; >> 87 nbins=320; //probably this should be decrease, that was choosen to avoid error in the CS value closed to CS jump.(For example at Tcut) >> 88 70 RegisterAdjointParticle(G4AdjointElectron::A 89 RegisterAdjointParticle(G4AdjointElectron::AdjointElectron()); 71 RegisterAdjointParticle(G4AdjointGamma::Adjo 90 RegisterAdjointParticle(G4AdjointGamma::AdjointGamma()); 72 RegisterAdjointParticle(G4AdjointProton::Adj 91 RegisterAdjointParticle(G4AdjointProton::AdjointProton()); 73 } << 92 >> 93 verbose = 1; >> 94 currentParticleIndex = 0; >> 95 currentMatIndex = 0; >> 96 eindex = 0; >> 97 >> 98 lastPartDefForCS = nullptr; >> 99 LastEkinForCS = lastPrimaryEnergy = lastTcut = 0.; >> 100 LastCSCorrectionFactor = massRatio = 1.; >> 101 >> 102 forward_CS_is_used = true; >> 103 forward_CS_mode = true; >> 104 >> 105 currentParticleDef = nullptr; >> 106 currentCouple =nullptr; >> 107 currentMaterial=nullptr; >> 108 lastMaterial=nullptr; 74 109 >> 110 theAdjIon = nullptr; >> 111 theFwdIon = nullptr; >> 112 >> 113 PreadjCS = PostadjCS = PrefwdCS = PostfwdCS = 0.0; >> 114 } 75 ////////////////////////////////////////////// 115 /////////////////////////////////////////////////////// >> 116 // 76 G4AdjointCSManager::~G4AdjointCSManager() 117 G4AdjointCSManager::~G4AdjointCSManager() 77 { << 118 {; 78 for (auto& p : fAdjointCSMatricesForProdToPr << 79 for (auto p1 : p) { << 80 if (p1) { << 81 delete p1; << 82 p1 = nullptr; << 83 } << 84 } << 85 p.clear(); << 86 } << 87 fAdjointCSMatricesForProdToProj.clear(); << 88 << 89 for (auto& p : fAdjointCSMatricesForScatProj << 90 for (auto p1 : p) { << 91 if (p1) { << 92 delete p1; << 93 p1 = nullptr; << 94 } << 95 } << 96 p.clear(); << 97 } << 98 fAdjointCSMatricesForScatProjToProj.clear(); << 99 << 100 for (auto p : fAdjointModels) { << 101 if (p) { << 102 delete p; << 103 p = nullptr; << 104 } << 105 } << 106 fAdjointModels.clear(); << 107 << 108 for (auto p : fTotalAdjSigmaTable) { << 109 p->clearAndDestroy(); << 110 delete p; << 111 p = nullptr; << 112 } << 113 fTotalAdjSigmaTable.clear(); << 114 << 115 for (auto p : fSigmaTableForAdjointModelScat << 116 p->clearAndDestroy(); << 117 delete p; << 118 p = nullptr; << 119 } << 120 fSigmaTableForAdjointModelScatProjToProj.cle << 121 << 122 for (auto p : fSigmaTableForAdjointModelProd << 123 p->clearAndDestroy(); << 124 delete p; << 125 p = nullptr; << 126 } << 127 fSigmaTableForAdjointModelProdToProj.clear() << 128 << 129 for (auto p : fTotalFwdSigmaTable) { << 130 p->clearAndDestroy(); << 131 delete p; << 132 p = nullptr; << 133 } << 134 fTotalFwdSigmaTable.clear(); << 135 << 136 for (auto p : fForwardProcesses) { << 137 delete p; << 138 p = nullptr; << 139 } << 140 fForwardProcesses.clear(); << 141 << 142 for (auto p : fForwardLossProcesses) { << 143 delete p; << 144 p = nullptr; << 145 } << 146 fForwardLossProcesses.clear(); << 147 } 119 } 148 << 149 ////////////////////////////////////////////// 120 /////////////////////////////////////////////////////// 150 std::size_t G4AdjointCSManager::RegisterEmAdjo << 121 // 151 { << 122 size_t G4AdjointCSManager::RegisterEmAdjointModel(G4VEmAdjointModel* aModel) 152 fAdjointModels.push_back(aModel); << 123 {listOfAdjointEMModel.push_back(aModel); 153 fSigmaTableForAdjointModelScatProjToProj.pus << 124 listSigmaTableForAdjointModelScatProjToProj.push_back(new G4PhysicsTable); 154 fSigmaTableForAdjointModelProdToProj.push_ba << 125 listSigmaTableForAdjointModelProdToProj.push_back(new G4PhysicsTable); 155 return fAdjointModels.size() - 1; << 126 return listOfAdjointEMModel.size() -1; >> 127 156 } 128 } 157 << 158 ////////////////////////////////////////////// 129 /////////////////////////////////////////////////////// 159 void G4AdjointCSManager::RegisterEmProcess(G4V << 130 // 160 G4P << 131 void G4AdjointCSManager::RegisterEmProcess(G4VEmProcess* aProcess, G4ParticleDefinition* aFwdPartDef) 161 { << 132 { 162 G4ParticleDefinition* anAdjPartDef = << 133 G4ParticleDefinition* anAdjPartDef = GetAdjointParticleEquivalent(aFwdPartDef); 163 GetAdjointParticleEquivalent(aFwdPartDef); << 134 if (anAdjPartDef && aProcess){ 164 if(anAdjPartDef && aProcess) << 135 RegisterAdjointParticle(anAdjPartDef); 165 { << 136 G4int index=-1; 166 RegisterAdjointParticle(anAdjPartDef); << 137 167 << 138 for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){ 168 for(std::size_t i = 0; i < fAdjointParticl << 139 if (anAdjPartDef->GetParticleName() == theListOfAdjointParticlesInAction[i]->GetParticleName()) index=i; 169 { << 140 } 170 if(anAdjPartDef->GetParticleName() == << 141 listOfForwardEmProcess[index]->push_back(aProcess); 171 fAdjointParticlesInAction[i]->GetPart << 172 fForwardProcesses[i]->push_back(aProce << 173 } << 174 } 142 } 175 } 143 } 176 << 177 ////////////////////////////////////////////// 144 /////////////////////////////////////////////////////// 178 void G4AdjointCSManager::RegisterEnergyLossPro << 145 // 179 G4VEnergyLossProcess* aProcess, G4ParticleDe << 146 void G4AdjointCSManager::RegisterEnergyLossProcess(G4VEnergyLossProcess* aProcess, G4ParticleDefinition* aFwdPartDef) 180 { << 147 { 181 G4ParticleDefinition* anAdjPartDef = << 148 G4ParticleDefinition* anAdjPartDef = GetAdjointParticleEquivalent(aFwdPartDef); 182 GetAdjointParticleEquivalent(aFwdPartDef); << 149 if (anAdjPartDef && aProcess){ 183 if(anAdjPartDef && aProcess) << 150 RegisterAdjointParticle(anAdjPartDef); 184 { << 151 G4int index=-1; 185 RegisterAdjointParticle(anAdjPartDef); << 152 for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){ 186 for(std::size_t i = 0; i < fAdjointParticl << 153 if (anAdjPartDef->GetParticleName() == theListOfAdjointParticlesInAction[i]->GetParticleName()) index=i; 187 { << 154 } 188 if(anAdjPartDef->GetParticleName() == << 155 listOfForwardEnergyLossProcess[index]->push_back(aProcess); 189 fAdjointParticlesInAction[i]->GetPart << 156 } 190 fForwardLossProc << 191 << 192 } << 193 } << 194 } 157 } 195 << 196 ////////////////////////////////////////////// 158 /////////////////////////////////////////////////////// >> 159 // 197 void G4AdjointCSManager::RegisterAdjointPartic 160 void G4AdjointCSManager::RegisterAdjointParticle(G4ParticleDefinition* aPartDef) 198 { << 161 { G4int index=-1; 199 G4bool found = false; << 162 for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){ 200 for(auto p : fAdjointParticlesInAction) << 163 if (aPartDef->GetParticleName() == theListOfAdjointParticlesInAction[i]->GetParticleName()) index=i; 201 { << 164 } 202 if(p->GetParticleName() == aPartDef->GetPa << 165 203 { << 166 if (index ==-1){ 204 found = true; << 167 listOfForwardEnergyLossProcess.push_back(new std::vector<G4VEnergyLossProcess*>()); 205 } << 168 theTotalForwardSigmaTableVector.push_back(new G4PhysicsTable); 206 } << 169 theTotalAdjointSigmaTableVector.push_back(new G4PhysicsTable); 207 if(!found) << 170 listOfForwardEmProcess.push_back(new std::vector<G4VEmProcess*>()); 208 { << 171 theListOfAdjointParticlesInAction.push_back(aPartDef); 209 fForwardLossProcesses.push_back(new std::v << 172 EminForFwdSigmaTables.push_back(std::vector<G4double> ()); 210 fTotalFwdSigmaTable.push_back(new G4Physic << 173 EminForAdjSigmaTables.push_back(std::vector<G4double> ()); 211 fTotalAdjSigmaTable.push_back(new G4Physic << 174 EkinofFwdSigmaMax.push_back(std::vector<G4double> ()); 212 fForwardProcesses.push_back(new std::vecto << 175 EkinofAdjSigmaMax.push_back(std::vector<G4double> ()); 213 fAdjointParticlesInAction.push_back(aPartD << 176 214 fEminForFwdSigmaTables.push_back(std::vect << 177 } 215 fEminForAdjSigmaTables.push_back(std::vect << 216 fEkinofFwdSigmaMax.push_back(std::vector<G << 217 fEkinofAdjSigmaMax.push_back(std::vector<G << 218 } << 219 } 178 } 220 << 221 ////////////////////////////////////////////// 179 /////////////////////////////////////////////////////// >> 180 // 222 void G4AdjointCSManager::BuildCrossSectionMatr 181 void G4AdjointCSManager::BuildCrossSectionMatrices() 223 { << 182 { 224 if(fCSMatricesBuilt) << 183 if (CrossSectionMatrixesAreBuilt) return; 225 return; << 184 //Tcut, Tmax 226 // The Tcut and Tmax matrices will be comput << 185 //The matrices will be computed probably just once 227 // When Tcut changes, some PhysicsTable will << 186 //When Tcut will change some PhysicsTable will be recomputed 228 // for each MaterialCutCouple but not all th << 187 // for each MaterialCutCouple but not all the matrices 229 // The Tcut defines a lower limit in the ene << 188 //The Tcut defines a lower limit in the energy of the Projectile before the scattering 230 // scattering. In the Projectile to Scattere << 189 //In the Projectile to Scattered Projectile case we have 231 // E_ScatProj<E_Proj-Tcut << 190 // E_ScatProj<E_Proj-Tcut 232 // Therefore in the adjoint case we have << 191 //Therefore in the adjoint case we have 233 // Eproj> E_ScatProj+Tcut << 192 // Eproj> E_ScatProj+Tcut 234 // This implies that when computing the adjo << 193 //This implies that when computing the adjoint CS we should integrate over Epro 235 // Epro from E_ScatProj+Tcut to Emax << 194 // from E_ScatProj+Tcut to Emax 236 // In the Projectile to Secondary case Tcut << 195 //In the Projectile to Secondary case Tcut plays a role only in the fact that 237 // Esecond should be greater than Tcut to ha << 196 // Esecond should be greater than Tcut to have the possibility to have any adjoint 238 // adjoint process. << 197 //process 239 // To avoid recomputing matrices for all cha << 198 //To avoid to recompute the matrices for all changes of MaterialCutCouple 240 // we propose to compute the matrices only o << 199 //We propose to compute the matrices only once for the minimum possible Tcut and then 241 // and then to interpolate the probability f << 200 //to interpolate the probability for a new Tcut (implemented in G4VAdjointEmModel) 242 // G4VAdjointEmModel) << 201 243 << 202 244 fAdjointCSMatricesForScatProjToProj.clear(); << 203 theAdjointCSMatricesForScatProjToProj.clear(); 245 fAdjointCSMatricesForProdToProj.clear(); << 204 theAdjointCSMatricesForProdToProj.clear(); 246 const G4ElementTable* theElementTable = G4 << 205 const G4ElementTable* theElementTable = G4Element::GetElementTable(); 247 const G4MaterialTable* theMaterialTable = G4 << 206 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable(); 248 << 207 249 G4cout << "========== Computation of cross s << 208 G4cout<<"========== Computation of cross section matrices for adjoint models =========="<<G4endl; 250 "models ==========" << 209 for (size_t i=0; i<listOfAdjointEMModel.size();i++){ 251 << G4endl; << 210 G4VEmAdjointModel* aModel =listOfAdjointEMModel[i]; 252 for(const auto& aModel : fAdjointModels) << 211 G4cout<<"Build adjoint cross section matrices for "<<aModel->GetName()<<G4endl; 253 { << 212 if (aModel->GetUseMatrix()){ 254 G4cout << "Build adjoint cross section mat << 213 std::vector<G4AdjointCSMatrix*>* aListOfMat1 = new std::vector<G4AdjointCSMatrix*>(); 255 << G4endl; << 214 std::vector<G4AdjointCSMatrix*>* aListOfMat2 = new std::vector<G4AdjointCSMatrix*>(); 256 if(aModel->GetUseMatrix()) << 215 aListOfMat1->clear(); 257 { << 216 aListOfMat2->clear(); 258 std::vector<G4AdjointCSMatrix*>* aListOf << 217 if (aModel->GetUseMatrixPerElement()){ 259 new std::vector<G4AdjointCSMatrix*>(); << 218 if (aModel->GetUseOnlyOneMatrixForAllElements()){ 260 std::vector<G4AdjointCSMatrix*>* aListOf << 219 std::vector<G4AdjointCSMatrix*> 261 new std::vector<G4AdjointCSMatrix*>(); << 220 two_matrices=BuildCrossSectionsMatricesForAGivenModelAndElement(aModel,1, 1, 80); 262 if(aModel->GetUseMatrixPerElement()) << 221 aListOfMat1->push_back(two_matrices[0]); 263 { << 222 aListOfMat2->push_back(two_matrices[1]); 264 if(aModel->GetUseOnlyOneMatrixForAllEl << 223 } 265 { << 224 else { 266 std::vector<G4AdjointCSMatrix*> two_ << 225 for (size_t j=0; j<theElementTable->size();j++){ 267 BuildCrossSectionsModelAndElement( << 226 G4Element* anElement=(*theElementTable)[j]; 268 aListOfMat1->push_back(two_matrices[ << 227 G4int Z = G4lrint(anElement->GetZ()); 269 aListOfMat2->push_back(two_matrices[ << 228 G4int A = G4lrint(anElement->GetN()); 270 } << 229 std::vector<G4AdjointCSMatrix*> 271 else << 230 two_matrices=BuildCrossSectionsMatricesForAGivenModelAndElement(aModel,Z, A, 40); 272 { << 231 aListOfMat1->push_back(two_matrices[0]); 273 for(const auto& anElement : *theElem << 232 aListOfMat2->push_back(two_matrices[1]); 274 { << 233 } 275 G4int Z = G4lrint(anElement->GetZ( << 234 } 276 G4int A = G4lrint(anElement->GetN( << 235 } 277 std::vector<G4AdjointCSMatrix*> tw << 236 else { //Per material case 278 BuildCrossSectionsModelAndElemen << 237 for (size_t j=0; j<theMaterialTable->size();j++){ 279 aListOfMat1->push_back(two_matrice << 238 G4Material* aMaterial=(*theMaterialTable)[j]; 280 aListOfMat2->push_back(two_matrice << 239 std::vector<G4AdjointCSMatrix*> 281 } << 240 two_matrices=BuildCrossSectionsMatricesForAGivenModelAndMaterial(aModel,aMaterial, 40); 282 } << 241 aListOfMat1->push_back(two_matrices[0]); 283 } << 242 aListOfMat2->push_back(two_matrices[1]); 284 else << 243 } 285 { // Per material case << 244 286 for(const auto& aMaterial : *theMateri << 245 } 287 { << 246 theAdjointCSMatricesForProdToProj.push_back(*aListOfMat1); 288 std::vector<G4AdjointCSMatrix*> two_ << 247 theAdjointCSMatricesForScatProjToProj.push_back(*aListOfMat2); 289 BuildCrossSectionsModelAndMaterial << 248 aModel->SetCSMatrices(aListOfMat1, aListOfMat2); 290 aListOfMat1->push_back(two_matrices[ << 249 } 291 aListOfMat2->push_back(two_matrices[ << 250 else { G4cout<<"The model "<<aModel->GetName()<<" does not use cross section matrices"<<G4endl; 292 } << 251 std::vector<G4AdjointCSMatrix*> two_empty_matrices; 293 } << 252 theAdjointCSMatricesForProdToProj.push_back(two_empty_matrices); 294 fAdjointCSMatricesForProdToProj.push_bac << 253 theAdjointCSMatricesForScatProjToProj.push_back(two_empty_matrices); 295 fAdjointCSMatricesForScatProjToProj.push << 254 296 aModel->SetCSMatrices(aListOfMat1, aList << 255 } 297 } << 256 } 298 else << 257 G4cout<<" All adjoint cross section matrices are computed!"<<G4endl; 299 { << 258 G4cout<<"======================================================================"<<G4endl; 300 G4cout << "The model " << aModel->GetNam << 259 301 << " does not use cross section m << 260 CrossSectionMatrixesAreBuilt = true; 302 std::vector<G4AdjointCSMatrix*> two_empt << 303 fAdjointCSMatricesForProdToProj.push_bac << 304 fAdjointCSMatricesForScatProjToProj.push << 305 } << 306 } << 307 G4cout << " All adjoint cross s << 308 << G4endl; << 309 G4cout << 310 << "====================================== << 311 << G4endl; << 312 261 313 fCSMatricesBuilt = true; << 314 } << 315 << 316 ////////////////////////////////////////////// << 317 void G4AdjointCSManager::BuildTotalSigmaTables << 318 { << 319 if(fSigmaTableBuilt) << 320 return; << 321 << 322 const G4ProductionCutsTable* theCoupleTable << 323 G4ProductionCutsTable::GetProductionCutsTa << 324 262 325 // Prepare the Sigma table for all AdjointEM << 326 for(std::size_t i = 0; i < fAdjointModels.si << 327 { << 328 fSigmaTableForAdjointModelScatProjToProj[i << 329 fSigmaTableForAdjointModelProdToProj[i]->c << 330 for(std::size_t j = 0; j < theCoupleTable- << 331 { << 332 fSigmaTableForAdjointModelScatProjToProj << 333 new G4PhysicsLogVector(fTmin, fTmax, f << 334 fSigmaTableForAdjointModelProdToProj[i]- << 335 new G4PhysicsLogVector(fTmin, fTmax, f << 336 } << 337 } << 338 << 339 for(std::size_t i = 0; i < fAdjointParticles << 340 { << 341 G4ParticleDefinition* thePartDef = fAdjoin << 342 DefineCurrentParticle(thePartDef); << 343 fTotalFwdSigmaTable[i]->clearAndDestroy(); << 344 fTotalAdjSigmaTable[i]->clearAndDestroy(); << 345 fEminForFwdSigmaTables[i].clear(); << 346 fEminForAdjSigmaTables[i].clear(); << 347 fEkinofFwdSigmaMax[i].clear(); << 348 fEkinofAdjSigmaMax[i].clear(); << 349 << 350 for(std::size_t j = 0; j < theCoupleTable- << 351 { << 352 const G4MaterialCutsCouple* couple = << 353 theCoupleTable->GetMaterialCutsCouple( << 354 << 355 // make first the total fwd CS table for << 356 G4PhysicsVector* aVector = new G4Physics << 357 G4bool Emin_found = false; << 358 G4double sigma_max = 0.; << 359 G4double e_sigma_max = 0.; << 360 for(std::size_t l = 0; l < fNbins; ++l) << 361 { << 362 G4double totCS = 0.; << 363 G4double e = aVector->Energy(l); << 364 for(std::size_t k = 0; k < fForwardPro << 365 { << 366 totCS += (*fForwardProcesses[i])[k]- << 367 } << 368 for(std::size_t k = 0; k < fForwardLos << 369 { << 370 if(thePartDef == fAdjIon) << 371 { // e is considered already as the << 372 std::size_t mat_index = couple->Ge << 373 G4VEmModel* currentModel = << 374 (*fForwardLossProcesses[i])[k]-> << 375 << 376 G4double chargeSqRatio = currentMo << 377 fFwdIon, couple->GetMaterial(), << 378 (*fForwardLossProcesses[i])[k]->Se << 379 << 380 } << 381 G4double e1 = e / fMassRatio; << 382 totCS += (*fForwardLossProcesses[i]) << 383 } << 384 aVector->PutValue(l, totCS); << 385 if(totCS > sigma_max) << 386 { << 387 sigma_max = totCS; << 388 e_sigma_max = e; << 389 } << 390 if(totCS > 0 && !Emin_found) << 391 { << 392 fEminForFwdSigmaTables[i].push_back( << 393 Emin_found = true; << 394 } << 395 } << 396 << 397 fEkinofFwdSigmaMax[i].push_back(e_sigma_ << 398 << 399 if(!Emin_found) << 400 fEminForFwdSigmaTables[i].push_back(fT << 401 << 402 fTotalFwdSigmaTable[i]->push_back(aVecto << 403 << 404 Emin_found = false; << 405 sigma_max = 0; << 406 e_sigma_max = 0.; << 407 G4PhysicsVector* aVector1 = new G4Physic << 408 for(std::size_t eindex = 0; eindex < fNb << 409 { << 410 G4double e = aVector1->Energy(eind << 411 G4double totCS = ComputeTotalAdjointCS << 412 couple, thePartDef, << 413 e * 0.9999999 / fMassRatio); // fMa << 414 aVector1->PutValue(eindex, totCS); << 415 if(totCS > sigma_max) << 416 { << 417 sigma_max = totCS; << 418 e_sigma_max = e; << 419 } << 420 if(totCS > 0 && !Emin_found) << 421 { << 422 fEminForAdjSigmaTables[i].push_back( << 423 Emin_found = true; << 424 } << 425 } << 426 fEkinofAdjSigmaMax[i].push_back(e_sigma_ << 427 if(!Emin_found) << 428 fEminForAdjSigmaTables[i].push_back(fT << 429 << 430 fTotalAdjSigmaTable[i]->push_back(aVecto << 431 } << 432 } << 433 fSigmaTableBuilt = true; << 434 } 263 } 435 264 436 ////////////////////////////////////////////// << 437 G4double G4AdjointCSManager::GetTotalAdjointCS << 438 G4ParticleDefinition* aPartDef, G4double Eki << 439 const G4MaterialCutsCouple* aCouple) << 440 { << 441 DefineCurrentMaterial(aCouple); << 442 DefineCurrentParticle(aPartDef); << 443 return (((*fTotalAdjSigmaTable[fCurrentParti << 444 ->Value(Ekin * fMassRatio)); << 445 } << 446 265 447 ////////////////////////////////////////////// 266 /////////////////////////////////////////////////////// 448 G4double G4AdjointCSManager::GetTotalForwardCS << 267 // 449 G4ParticleDefinition* aPartDef, G4double Eki << 268 void G4AdjointCSManager::BuildTotalSigmaTables() 450 const G4MaterialCutsCouple* aCouple) << 269 { if (TotalSigmaTableAreBuilt) return; 451 { << 452 DefineCurrentMaterial(aCouple); << 453 DefineCurrentParticle(aPartDef); << 454 return (((*fTotalFwdSigmaTable[fCurrentParti << 455 ->Value(Ekin * fMassRatio)); << 456 } << 457 << 458 ////////////////////////////////////////////// << 459 G4double G4AdjointCSManager::GetAdjointSigma( << 460 G4double Ekin_nuc, std::size_t index_model, << 461 const G4MaterialCutsCouple* aCouple) << 462 { << 463 DefineCurrentMaterial(aCouple); << 464 if(is_scat_proj_to_proj) << 465 return (((*fSigmaTableForAdjointModelScatP << 466 [fCurrentMatIndex])->Value(Ekin << 467 else << 468 return ( << 469 ((*fSigmaTableForAdjointModelProdToProj[ << 470 ->Value(Ekin_nuc)); << 471 } << 472 270 >> 271 >> 272 const G4ProductionCutsTable* theCoupleTable= G4ProductionCutsTable::GetProductionCutsTable(); >> 273 >> 274 >> 275 //Prepare the Sigma table for all AdjointEMModel, will be filled later on >> 276 for (size_t i=0; i<listOfAdjointEMModel.size();i++){ >> 277 listSigmaTableForAdjointModelScatProjToProj[i]->clearAndDestroy(); >> 278 listSigmaTableForAdjointModelProdToProj[i]->clearAndDestroy(); >> 279 for (size_t j=0;j<theCoupleTable->GetTableSize();j++){ >> 280 listSigmaTableForAdjointModelScatProjToProj[i]->push_back(new G4PhysicsLogVector(Tmin, Tmax, nbins)); >> 281 listSigmaTableForAdjointModelProdToProj[i]->push_back(new G4PhysicsLogVector(Tmin, Tmax, nbins)); >> 282 } >> 283 } >> 284 >> 285 >> 286 >> 287 for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){ >> 288 G4ParticleDefinition* thePartDef = theListOfAdjointParticlesInAction[i]; >> 289 DefineCurrentParticle(thePartDef); >> 290 theTotalForwardSigmaTableVector[i]->clearAndDestroy(); >> 291 theTotalAdjointSigmaTableVector[i]->clearAndDestroy(); >> 292 EminForFwdSigmaTables[i].clear(); >> 293 EminForAdjSigmaTables[i].clear(); >> 294 EkinofFwdSigmaMax[i].clear(); >> 295 EkinofAdjSigmaMax[i].clear(); >> 296 //G4cout<<thePartDef->GetParticleName(); >> 297 >> 298 for (size_t j=0;j<theCoupleTable->GetTableSize();j++){ >> 299 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(j); >> 300 >> 301 /* >> 302 G4String file_name1=couple->GetMaterial()->GetName()+"_"+thePartDef->GetParticleName()+"_adj_totCS.txt"; >> 303 G4String file_name2=couple->GetMaterial()->GetName()+"_"+thePartDef->GetParticleName()+"_fwd_totCS.txt"; >> 304 >> 305 std::fstream FileOutputAdjCS(file_name1, std::ios::out); >> 306 std::fstream FileOutputFwdCS(file_name2, std::ios::out); >> 307 >> 308 >> 309 >> 310 FileOutputAdjCS<<std::setiosflags(std::ios::scientific); >> 311 FileOutputAdjCS<<std::setprecision(6); >> 312 FileOutputFwdCS<<std::setiosflags(std::ios::scientific); >> 313 FileOutputFwdCS<<std::setprecision(6); >> 314 */ >> 315 >> 316 >> 317 //make first the total fwd CS table for FwdProcess >> 318 G4PhysicsVector* aVector = new G4PhysicsLogVector(Tmin, Tmax, nbins); >> 319 G4bool Emin_found=false; >> 320 G4double sigma_max =0.; >> 321 G4double e_sigma_max =0.; >> 322 for(size_t l=0; l<aVector->GetVectorLength(); l++) { >> 323 G4double totCS=0.; >> 324 G4double e=aVector->GetLowEdgeEnergy(l); >> 325 for (size_t k=0; k<listOfForwardEmProcess[i]->size(); k++){ >> 326 totCS+=(*listOfForwardEmProcess[i])[k]->GetLambda(e, couple); >> 327 } >> 328 for (size_t k=0; k<listOfForwardEnergyLossProcess[i]->size(); k++){ >> 329 if (thePartDef == theAdjIon) { // e is considered already as the scaled energy >> 330 size_t mat_index = couple->GetIndex(); >> 331 G4VEmModel* currentModel = (*listOfForwardEnergyLossProcess[i])[k]->SelectModelForMaterial(e,mat_index); >> 332 G4double chargeSqRatio = currentModel->GetChargeSquareRatio(theFwdIon,couple->GetMaterial(),e/massRatio); >> 333 (*listOfForwardEnergyLossProcess[i])[k]->SetDynamicMassCharge(massRatio,chargeSqRatio); >> 334 } >> 335 G4double e1=e/massRatio; >> 336 totCS+=(*listOfForwardEnergyLossProcess[i])[k]->GetLambda(e1, couple); >> 337 } >> 338 aVector->PutValue(l,totCS); >> 339 if (totCS>sigma_max){ >> 340 sigma_max=totCS; >> 341 e_sigma_max = e; >> 342 >> 343 } >> 344 //FileOutputFwdCS<<e<<'\t'<<totCS<<G4endl; >> 345 >> 346 if (totCS>0 && !Emin_found) { >> 347 EminForFwdSigmaTables[i].push_back(e); >> 348 Emin_found=true; >> 349 } >> 350 >> 351 >> 352 } >> 353 //FileOutputFwdCS.close(); >> 354 >> 355 EkinofFwdSigmaMax[i].push_back(e_sigma_max); >> 356 >> 357 >> 358 if(!Emin_found) EminForFwdSigmaTables[i].push_back(Tmax); >> 359 >> 360 theTotalForwardSigmaTableVector[i]->push_back(aVector); >> 361 >> 362 >> 363 Emin_found=false; >> 364 sigma_max=0; >> 365 e_sigma_max =0.; >> 366 G4PhysicsVector* aVector1 = new G4PhysicsLogVector(Tmin, Tmax, nbins); >> 367 for(eindex=0; eindex<aVector->GetVectorLength(); eindex++) { >> 368 G4double e=aVector->GetLowEdgeEnergy(eindex); >> 369 G4double totCS =ComputeTotalAdjointCS(couple,thePartDef,e*0.9999999/massRatio); //massRatio needed for ions >> 370 aVector1->PutValue(eindex,totCS); >> 371 if (totCS>sigma_max){ >> 372 sigma_max=totCS; >> 373 e_sigma_max = e; >> 374 >> 375 } >> 376 //FileOutputAdjCS<<e<<'\t'<<totCS<<G4endl; >> 377 if (totCS>0 && !Emin_found) { >> 378 EminForAdjSigmaTables[i].push_back(e); >> 379 Emin_found=true; >> 380 } >> 381 >> 382 } >> 383 //FileOutputAdjCS.close(); >> 384 EkinofAdjSigmaMax[i].push_back(e_sigma_max); >> 385 if(!Emin_found) EminForAdjSigmaTables[i].push_back(Tmax); >> 386 >> 387 theTotalAdjointSigmaTableVector[i]->push_back(aVector1); >> 388 >> 389 } >> 390 } >> 391 TotalSigmaTableAreBuilt =true; >> 392 >> 393 } >> 394 /////////////////////////////////////////////////////// >> 395 // >> 396 G4double G4AdjointCSManager::GetTotalAdjointCS(G4ParticleDefinition* aPartDef, G4double Ekin, >> 397 const G4MaterialCutsCouple* aCouple) >> 398 { DefineCurrentMaterial(aCouple); >> 399 DefineCurrentParticle(aPartDef); >> 400 G4bool b; >> 401 return (((*theTotalAdjointSigmaTableVector[currentParticleIndex])[currentMatIndex])->GetValue(Ekin*massRatio, b)); >> 402 >> 403 >> 404 >> 405 } >> 406 /////////////////////////////////////////////////////// >> 407 // >> 408 G4double G4AdjointCSManager::GetTotalForwardCS(G4ParticleDefinition* aPartDef, G4double Ekin, >> 409 const G4MaterialCutsCouple* aCouple) >> 410 { DefineCurrentMaterial(aCouple); >> 411 DefineCurrentParticle(aPartDef); >> 412 G4bool b; >> 413 return (((*theTotalForwardSigmaTableVector[currentParticleIndex])[currentMatIndex])->GetValue(Ekin*massRatio, b)); >> 414 >> 415 >> 416 } >> 417 /////////////////////////////////////////////////////// >> 418 // >> 419 G4double G4AdjointCSManager::GetAdjointSigma(G4double Ekin_nuc, size_t index_model,G4bool is_scat_proj_to_proj, >> 420 const G4MaterialCutsCouple* aCouple) >> 421 { DefineCurrentMaterial(aCouple); >> 422 G4bool b; >> 423 if (is_scat_proj_to_proj) return (((*listSigmaTableForAdjointModelScatProjToProj[index_model])[currentMatIndex])->GetValue(Ekin_nuc, b)); >> 424 else return (((*listSigmaTableForAdjointModelProdToProj[index_model])[currentMatIndex])->GetValue(Ekin_nuc, b)); >> 425 } 473 ////////////////////////////////////////////// 426 /////////////////////////////////////////////////////// >> 427 // 474 void G4AdjointCSManager::GetEminForTotalCS(G4P 428 void G4AdjointCSManager::GetEminForTotalCS(G4ParticleDefinition* aPartDef, 475 con << 429 const G4MaterialCutsCouple* aCouple, G4double& emin_adj, G4double& emin_fwd) 476 G4d << 430 { DefineCurrentMaterial(aCouple); 477 G4d << 478 { << 479 DefineCurrentMaterial(aCouple); << 480 DefineCurrentParticle(aPartDef); 431 DefineCurrentParticle(aPartDef); 481 emin_adj = fEminForAdjSigmaTables[fCurrentPa << 432 emin_adj = EminForAdjSigmaTables[currentParticleIndex][currentMatIndex]/massRatio; 482 fMassRatio; << 433 emin_fwd = EminForFwdSigmaTables[currentParticleIndex][currentMatIndex]/massRatio; 483 emin_fwd = fEminForFwdSigmaTables[fCurrentPa << 434 484 fMassRatio; << 435 >> 436 485 } 437 } 486 << 487 ////////////////////////////////////////////// 438 /////////////////////////////////////////////////////// >> 439 // 488 void G4AdjointCSManager::GetMaxFwdTotalCS(G4Pa 440 void G4AdjointCSManager::GetMaxFwdTotalCS(G4ParticleDefinition* aPartDef, 489 cons << 441 const G4MaterialCutsCouple* aCouple, G4double& e_sigma_max, G4double& sigma_max) 490 G4do << 442 { DefineCurrentMaterial(aCouple); 491 G4do << 492 { << 493 DefineCurrentMaterial(aCouple); << 494 DefineCurrentParticle(aPartDef); 443 DefineCurrentParticle(aPartDef); 495 e_sigma_max = fEkinofFwdSigmaMax[fCurrentPar << 444 e_sigma_max = EkinofFwdSigmaMax[currentParticleIndex][currentMatIndex]; 496 sigma_max = ((*fTotalFwdSigmaTable[fCurrentP << 445 G4bool b; 497 ->Value(e_sigma_max); << 446 sigma_max =((*theTotalForwardSigmaTableVector[currentParticleIndex])[currentMatIndex])->GetValue(e_sigma_max, b); 498 e_sigma_max /= fMassRatio; << 447 e_sigma_max/=massRatio; >> 448 >> 449 499 } 450 } 500 << 501 ////////////////////////////////////////////// 451 /////////////////////////////////////////////////////// >> 452 // 502 void G4AdjointCSManager::GetMaxAdjTotalCS(G4Pa 453 void G4AdjointCSManager::GetMaxAdjTotalCS(G4ParticleDefinition* aPartDef, 503 cons << 454 const G4MaterialCutsCouple* aCouple, G4double& e_sigma_max, G4double& sigma_max) 504 G4do << 455 { DefineCurrentMaterial(aCouple); 505 G4do << 506 { << 507 DefineCurrentMaterial(aCouple); << 508 DefineCurrentParticle(aPartDef); 456 DefineCurrentParticle(aPartDef); 509 e_sigma_max = fEkinofAdjSigmaMax[fCurrentPar << 457 e_sigma_max = EkinofAdjSigmaMax[currentParticleIndex][currentMatIndex]; 510 sigma_max = ((*fTotalAdjSigmaTable[fCurrentP << 458 G4bool b; 511 ->Value(e_sigma_max); << 459 sigma_max =((*theTotalAdjointSigmaTableVector[currentParticleIndex])[currentMatIndex])->GetValue(e_sigma_max, b); 512 e_sigma_max /= fMassRatio; << 460 e_sigma_max/=massRatio; 513 } << 461 514 << 462 515 ////////////////////////////////////////////// << 463 } 516 G4double G4AdjointCSManager::GetCrossSectionCo << 464 /////////////////////////////////////////////////////// 517 G4ParticleDefinition* aPartDef, G4double Pre << 465 // 518 const G4MaterialCutsCouple* aCouple, G4bool& << 466 G4double G4AdjointCSManager::GetCrossSectionCorrection(G4ParticleDefinition* aPartDef,G4double PreStepEkin,const G4MaterialCutsCouple* aCouple, G4bool& fwd_is_used, 519 { << 467 G4double& fwd_TotCS) 520 static G4double lastEkin = 0.; << 468 { G4double corr_fac = 1.; 521 static G4ParticleDefinition* lastPartDef; << 469 if (forward_CS_mode && aPartDef ) { 522 << 470 fwd_TotCS=PrefwdCS; 523 G4double corr_fac = 1.; << 471 if (LastEkinForCS != PreStepEkin || aPartDef != lastPartDefForCS || aCouple!=currentCouple) { 524 if(fForwardCSMode && aPartDef) << 472 DefineCurrentMaterial(aCouple); 525 { << 473 PreadjCS = GetTotalAdjointCS(aPartDef, PreStepEkin,aCouple); 526 if(lastEkin != PreStepEkin || aPartDef != << 474 PrefwdCS = GetTotalForwardCS(aPartDef, PreStepEkin,aCouple); 527 aCouple != fCurrentCouple) << 475 LastEkinForCS = PreStepEkin; 528 { << 476 lastPartDefForCS = aPartDef; 529 DefineCurrentMaterial(aCouple); << 477 if (PrefwdCS >0. && PreadjCS >0.) { 530 G4double preadjCS = GetTotalAdjointCS(aP << 478 forward_CS_is_used = true; 531 G4double prefwdCS = GetTotalForwardCS(aP << 479 LastCSCorrectionFactor = PrefwdCS/PreadjCS; 532 lastEkin = PreStepEkin; << 480 } 533 lastPartDef = aPartDef; << 481 else { 534 if(prefwdCS > 0. && preadjCS > 0.) << 482 forward_CS_is_used = false; 535 { << 483 LastCSCorrectionFactor = 1.; 536 fForwardCSUsed = true; << 484 537 fLastCSCorrectionFactor = prefwdCS / p << 485 } 538 } << 486 539 else << 487 } 540 { << 488 corr_fac =LastCSCorrectionFactor; 541 fForwardCSUsed = false; << 489 542 fLastCSCorrectionFactor = 1.; << 490 543 } << 491 544 } << 492 } 545 corr_fac = fLastCSCorrectionFactor; << 493 else { 546 } << 494 forward_CS_is_used = false; 547 else << 495 LastCSCorrectionFactor = 1.; 548 { << 496 } 549 fForwardCSUsed = false; << 497 fwd_TotCS=PrefwdCS; 550 fLastCSCorrectionFactor = 1.; << 498 fwd_is_used = forward_CS_is_used; 551 } << 499 return corr_fac; 552 fwd_is_used = fForwardCSUsed; << 500 } 553 return corr_fac; << 501 /////////////////////////////////////////////////////// 554 } << 502 // 555 << 503 G4double G4AdjointCSManager::GetContinuousWeightCorrection(G4ParticleDefinition* aPartDef, G4double PreStepEkin,G4double AfterStepEkin, 556 ////////////////////////////////////////////// << 504 const G4MaterialCutsCouple* aCouple, G4double step_length) 557 G4double G4AdjointCSManager::GetContinuousWeig << 505 { G4double corr_fac = 1.; 558 G4ParticleDefinition* aPartDef, G4double Pre << 506 //return corr_fac; 559 const G4MaterialCutsCouple* aCouple, G4doubl << 507 //G4double after_adjCS = GetTotalAdjointCS(aPartDef, AfterStepEkin,aCouple); 560 { << 508 G4double after_fwdCS = GetTotalForwardCS(aPartDef, AfterStepEkin,aCouple); 561 G4double corr_fac = 1.; << 509 G4double pre_adjCS = GetTotalAdjointCS(aPartDef, PreStepEkin,aCouple); 562 G4double after_fwdCS = GetTotalForwardCS(aPa << 510 if (!forward_CS_is_used || pre_adjCS ==0. || after_fwdCS==0.) { 563 G4double pre_adjCS = GetTotalAdjointCS(aPa << 511 forward_CS_is_used=false; 564 if(!fForwardCSUsed || pre_adjCS == 0. || aft << 512 G4double pre_fwdCS = GetTotalForwardCS(aPartDef, PreStepEkin,aCouple); 565 { << 513 corr_fac *=std::exp((pre_adjCS-pre_fwdCS)*step_length); 566 G4double pre_fwdCS = GetTotalForwardCS(aPa << 514 LastCSCorrectionFactor = 1.; 567 corr_fac *= std::exp((pre_adjCS - pre_fwdC << 515 } 568 fLastCSCorrectionFactor = 1.; << 516 else { 569 } << 517 LastCSCorrectionFactor = after_fwdCS/pre_adjCS; 570 else << 518 } 571 { << 519 572 fLastCSCorrectionFactor = after_fwdCS / pr << 573 } << 574 return corr_fac; << 575 } << 576 << 577 ////////////////////////////////////////////// << 578 G4double G4AdjointCSManager::GetPostStepWeight << 579 { << 580 return 1. / fLastCSCorrectionFactor; << 581 } << 582 520 >> 521 >> 522 return corr_fac; >> 523 } 583 ////////////////////////////////////////////// 524 /////////////////////////////////////////////////////// 584 G4double G4AdjointCSManager::ComputeAdjointCS( << 525 // 585 G4Material* aMaterial, G4VEmAdjointModel* aM << 526 G4double G4AdjointCSManager::GetPostStepWeightCorrection( ) 586 G4double Tcut, G4bool isScatProjToProj, std: << 527 {//return 1.; 587 { << 528 return 1./LastCSCorrectionFactor; 588 G4double EminSec = 0.; << 529 589 G4double EmaxSec = 0.; << 530 } 590 << 531 /////////////////////////////////////////////////////// 591 static G4double lastPrimaryEnergy = 0.; << 532 // 592 static G4double lastTcut = 0.; << 533 G4double G4AdjointCSManager::ComputeAdjointCS(G4Material* aMaterial, 593 static G4Material* lastMaterial = nullptr; << 534 G4VEmAdjointModel* aModel, 594 << 535 G4double PrimEnergy, 595 if(isScatProjToProj) << 536 G4double Tcut, 596 { << 537 G4bool IsScatProjToProjCase, 597 EminSec = aModel->GetSecondAdjEnergyMinFor << 538 std::vector<G4double>& CS_Vs_Element) 598 EmaxSec = aModel->GetSecondAdjEnergyMaxFor << 539 { 599 } << 540 600 else if(PrimEnergy > Tcut || !aModel->GetApp << 541 G4double EminSec=0; 601 { << 542 G4double EmaxSec=0; 602 EminSec = aModel->GetSecondAdjEnergyMinFor << 543 603 EmaxSec = aModel->GetSecondAdjEnergyMaxFor << 544 if (IsScatProjToProjCase){ 604 } << 545 EminSec= aModel->GetSecondAdjEnergyMinForScatProjToProjCase(PrimEnergy,Tcut); 605 if(EminSec >= EmaxSec) << 546 EmaxSec= aModel->GetSecondAdjEnergyMaxForScatProjToProjCase(PrimEnergy); 606 return 0.; << 547 } 607 << 548 else if (PrimEnergy > Tcut || !aModel->GetApplyCutInRange()) { 608 G4bool need_to_compute = false; << 549 EminSec= aModel->GetSecondAdjEnergyMinForProdToProjCase(PrimEnergy); 609 if(aMaterial != lastMaterial || PrimEnergy ! << 550 EmaxSec= aModel->GetSecondAdjEnergyMaxForProdToProjCase(PrimEnergy); 610 Tcut != lastTcut) << 551 } 611 { << 552 if (EminSec >= EmaxSec) return 0.; 612 lastMaterial = aMaterial; << 553 613 lastPrimaryEnergy = PrimEnergy; << 554 614 lastTcut = Tcut; << 555 G4bool need_to_compute=false; 615 fIndexOfAdjointEMModelInAction.clear(); << 556 if ( aMaterial!= lastMaterial || PrimEnergy != lastPrimaryEnergy || Tcut != lastTcut){ 616 fIsScatProjToProj.clear(); << 557 lastMaterial =aMaterial; 617 fLastAdjointCSVsModelsAndElements.clear(); << 558 lastPrimaryEnergy = PrimEnergy; 618 need_to_compute = true; << 559 lastTcut=Tcut; 619 } << 560 listOfIndexOfAdjointEMModelInAction.clear(); 620 << 561 listOfIsScatProjToProjCase.clear(); 621 std::size_t ind = 0; << 562 lastAdjointCSVsModelsAndElements.clear(); 622 if(!need_to_compute) << 563 need_to_compute=true; 623 { << 564 624 need_to_compute = true; << 565 } 625 for(std::size_t i = 0; i < fIndexOfAdjoint << 566 size_t ind=0; 626 { << 567 if (!need_to_compute){ 627 std::size_t ind1 = fIndexOfAdjointEMMode << 568 need_to_compute=true; 628 if(aModel == fAdjointModels[ind1] && << 569 for (size_t i=0;i<listOfIndexOfAdjointEMModelInAction.size();i++){ 629 isScatProjToProj == fIsScatProjToProj << 570 size_t ind1=listOfIndexOfAdjointEMModelInAction[i]; 630 { << 571 if (aModel == listOfAdjointEMModel[ind1] && IsScatProjToProjCase == listOfIsScatProjToProjCase[i]){ 631 need_to_compute = false; << 572 need_to_compute=false; 632 CS_Vs_Element = fLastAdjointCSVsMode << 573 CS_Vs_Element = lastAdjointCSVsModelsAndElements[ind]; 633 } << 574 } 634 ++ind; << 575 ind++; 635 } << 576 } 636 } << 577 } 637 << 578 638 if(need_to_compute) << 579 if (need_to_compute){ 639 { << 580 size_t ind_model=0; 640 std::size_t ind_model = 0; << 581 for (size_t i=0;i<listOfAdjointEMModel.size();i++){ 641 for(std::size_t i = 0; i < fAdjointModels. << 582 if (aModel == listOfAdjointEMModel[i]){ 642 { << 583 ind_model=i; 643 if(aModel == fAdjointModels[i]) << 584 break; 644 { << 585 } 645 ind_model = i; << 586 } 646 break; << 587 G4double Tlow=Tcut; 647 } << 588 if (!listOfAdjointEMModel[ind_model]->GetApplyCutInRange()) Tlow =listOfAdjointEMModel[ind_model]->GetLowEnergyLimit(); 648 } << 589 listOfIndexOfAdjointEMModelInAction.push_back(ind_model); 649 G4double Tlow = Tcut; << 590 listOfIsScatProjToProjCase.push_back(IsScatProjToProjCase); 650 if(!fAdjointModels[ind_model]->GetApplyCut << 591 CS_Vs_Element.clear(); 651 Tlow = fAdjointModels[ind_model]->GetLow << 592 if (!aModel->GetUseMatrix()){ 652 fIndexOfAdjointEMModelInAction.push_back(i << 593 CS_Vs_Element.push_back(aModel->AdjointCrossSection(currentCouple,PrimEnergy,IsScatProjToProjCase)); 653 fIsScatProjToProj.push_back(isScatProjToPr << 594 654 CS_Vs_Element.clear(); << 595 655 if(!aModel->GetUseMatrix()) << 596 } 656 { << 597 else if (aModel->GetUseMatrixPerElement()){ 657 CS_Vs_Element.push_back(aModel->AdjointC << 598 size_t n_el = aMaterial->GetNumberOfElements(); 658 fCurrentCouple, PrimEnergy, isScatProj << 599 if (aModel->GetUseOnlyOneMatrixForAllElements()){ 659 } << 600 G4AdjointCSMatrix* theCSMatrix; 660 else if(aModel->GetUseMatrixPerElement()) << 601 if (IsScatProjToProjCase){ 661 { << 602 theCSMatrix=theAdjointCSMatricesForScatProjToProj[ind_model][0]; 662 std::size_t n_el = aMaterial->GetNumberO << 603 } 663 if(aModel->GetUseOnlyOneMatrixForAllElem << 604 else theCSMatrix=theAdjointCSMatricesForProdToProj[ind_model][0]; 664 { << 605 G4double CS =0.; 665 G4AdjointCSMatrix* theCSMatrix; << 606 if (PrimEnergy > Tlow) 666 if(isScatProjToProj) << 607 CS = ComputeAdjointCS(PrimEnergy,theCSMatrix,Tlow); 667 { << 608 G4double factor=0.; 668 theCSMatrix = fAdjointCSMatricesForS << 609 for (size_t i=0;i<n_el;i++){ //this could be computed only once 669 } << 610 //size_t ind_el = aMaterial->GetElement(i)->GetIndex(); 670 else << 611 factor+=aMaterial->GetElement(i)->GetZ()*aMaterial->GetVecNbOfAtomsPerVolume()[i]; 671 theCSMatrix = fAdjointCSMatricesForP << 612 } 672 G4double CS = 0.; << 613 CS *=factor; 673 if(PrimEnergy > Tlow) << 614 CS_Vs_Element.push_back(CS); 674 CS = ComputeAdjointCS(PrimEnergy, th << 615 675 G4double factor = 0.; << 616 } 676 for(G4int i = 0; i < (G4int)n_el; ++i) << 617 else { 677 { // this could be computed only once << 618 for (size_t i=0;i<n_el;i++){ 678 factor += aMaterial->GetElement(i)-> << 619 size_t ind_el = aMaterial->GetElement(i)->GetIndex(); 679 aMaterial->GetVecNbOfAtoms << 620 //G4cout<<aMaterial->GetName()<<G4endl; 680 } << 621 G4AdjointCSMatrix* theCSMatrix; 681 CS *= factor; << 622 if (IsScatProjToProjCase){ 682 CS_Vs_Element.push_back(CS); << 623 theCSMatrix=theAdjointCSMatricesForScatProjToProj[ind_model][ind_el]; 683 } << 624 } 684 else << 625 else theCSMatrix=theAdjointCSMatricesForProdToProj[ind_model][ind_el]; 685 { << 626 G4double CS =0.; 686 for(G4int i = 0; i < (G4int)n_el; ++i) << 627 if (PrimEnergy > Tlow) 687 { << 628 CS = ComputeAdjointCS(PrimEnergy,theCSMatrix,Tlow); 688 std::size_t ind_el = aMaterial->GetE << 629 //G4cout<<CS<<G4endl; 689 G4AdjointCSMatrix* theCSMatrix; << 630 CS_Vs_Element.push_back(CS*(aMaterial->GetVecNbOfAtomsPerVolume()[i])); 690 if(isScatProjToProj) << 631 } 691 { << 632 } 692 theCSMatrix = << 633 693 fAdjointCSMatricesForScatProjToP << 634 } 694 } << 635 else { 695 else << 636 size_t ind_mat = aMaterial->GetIndex(); 696 theCSMatrix = fAdjointCSMatricesFo << 637 G4AdjointCSMatrix* theCSMatrix; 697 G4double CS = 0.; << 638 if (IsScatProjToProjCase){ 698 if(PrimEnergy > Tlow) << 639 theCSMatrix=theAdjointCSMatricesForScatProjToProj[ind_model][ind_mat]; 699 CS = ComputeAdjointCS(PrimEnergy, << 640 } 700 CS_Vs_Element.push_back(CS * << 641 else theCSMatrix=theAdjointCSMatricesForProdToProj[ind_model][ind_mat]; 701 (aMaterial-> << 642 G4double CS =0.; 702 } << 643 if (PrimEnergy > Tlow) 703 } << 644 CS = ComputeAdjointCS(PrimEnergy,theCSMatrix,Tlow); 704 } << 645 CS_Vs_Element.push_back(CS); 705 else << 646 706 { << 647 707 std::size_t ind_mat = aMaterial->GetInde << 648 } 708 G4AdjointCSMatrix* theCSMatrix; << 649 lastAdjointCSVsModelsAndElements.push_back(CS_Vs_Element); 709 if(isScatProjToProj) << 650 710 { << 651 } 711 theCSMatrix = fAdjointCSMatricesForSca << 652 712 } << 653 713 else << 654 G4double CS=0; 714 theCSMatrix = fAdjointCSMatricesForPro << 655 for (size_t i=0;i<CS_Vs_Element.size();i++){ 715 G4double CS = 0.; << 656 CS+=CS_Vs_Element[i]; //We could put the progressive sum of the CS instead of the CS of an element itself 716 if(PrimEnergy > Tlow) << 657 717 CS = ComputeAdjointCS(PrimEnergy, theC << 718 CS_Vs_Element.push_back(CS); << 719 } << 720 fLastAdjointCSVsModelsAndElements.push_bac << 721 } << 722 << 723 G4double CS = 0.; << 724 for(const auto& cs_vs_el : CS_Vs_Element) << 725 { << 726 // We could put the progressive sum of the << 727 // element itself << 728 CS += cs_vs_el; << 729 } 658 } 730 return CS; 659 return CS; 731 } << 660 } 732 << 733 ////////////////////////////////////////////// 661 /////////////////////////////////////////////////////// 734 G4Element* G4AdjointCSManager::SampleElementFr << 662 // 735 G4Material* aMaterial, G4VEmAdjointModel* aM << 663 G4Element* G4AdjointCSManager::SampleElementFromCSMatrices(G4Material* aMaterial, 736 G4double Tcut, G4bool isScatProjToProj) << 664 G4VEmAdjointModel* aModel, 737 { << 665 G4double PrimEnergy, 738 std::vector<G4double> CS_Vs_Element; << 666 G4double Tcut, 739 G4double CS = ComputeAdjointCS(aMaterial, << 667 G4bool IsScatProjToProjCase) 740 isScatProjToP << 668 { std::vector<G4double> CS_Vs_Element; 741 G4double SumCS = 0.; << 669 G4double CS = ComputeAdjointCS(aMaterial,aModel,PrimEnergy,Tcut,IsScatProjToProjCase,CS_Vs_Element); 742 std::size_t ind = 0; << 670 G4double rand_var= G4UniformRand(); 743 for(std::size_t i = 0; i < CS_Vs_Element.siz << 671 G4double SumCS=0.; 744 { << 672 size_t ind=0; 745 SumCS += CS_Vs_Element[i]; << 673 for (size_t i=0;i<CS_Vs_Element.size();i++){ 746 if(G4UniformRand() <= SumCS / CS) << 674 SumCS+=CS_Vs_Element[i]; 747 { << 675 if (rand_var<=SumCS/CS){ 748 ind = i; << 676 ind=i; 749 break; << 677 break; 750 } << 678 } 751 } 679 } 752 << 680 753 return const_cast<G4Element*>(aMaterial->Get << 681 return const_cast<G4Element*>(aMaterial->GetElement(ind)); 754 } << 682 755 << 683 >> 684 >> 685 } 756 ////////////////////////////////////////////// 686 /////////////////////////////////////////////////////// 757 G4double G4AdjointCSManager::ComputeTotalAdjoi << 687 // 758 const G4MaterialCutsCouple* aCouple, G4Parti << 688 G4double G4AdjointCSManager::ComputeTotalAdjointCS(const G4MaterialCutsCouple* aCouple, 759 G4double Ekin) << 689 G4ParticleDefinition* aPartDef, 760 { << 690 G4double Ekin) 761 G4double TotalCS = 0.; << 691 { 762 << 692 G4double TotalCS=0.; 763 DefineCurrentMaterial(aCouple); << 693 764 << 694 DefineCurrentMaterial(aCouple); 765 std::vector<G4double> CS_Vs_Element; << 695 766 G4double CS; << 696 767 G4VEmAdjointModel* adjModel = nullptr; << 697 std::vector<G4double> CS_Vs_Element; 768 for(std::size_t i = 0; i < fAdjointModels.si << 698 G4double CS; 769 { << 699 for (size_t i=0; i<listOfAdjointEMModel.size();i++){ 770 G4double Tlow = 0.; << 700 771 adjModel = fAdjointModels[i]; << 701 G4double Tlow=0; 772 if(!adjModel->GetApplyCutInRange()) << 702 if (!listOfAdjointEMModel[i]->GetApplyCutInRange()) Tlow =listOfAdjointEMModel[i]->GetLowEnergyLimit(); 773 Tlow = adjModel->GetLowEnergyLimit(); << 703 else { 774 else << 704 G4ParticleDefinition* theDirSecondPartDef = 775 { << 705 GetForwardParticleEquivalent(listOfAdjointEMModel[i]->GetAdjointEquivalentOfDirectSecondaryParticleDefinition()); 776 G4ParticleDefinition* theDirSecondPartDe << 706 size_t idx=56; 777 adjModel->GetAdjointEquivalentOfDirect << 707 if (theDirSecondPartDef->GetParticleName() == "gamma") idx = 0; 778 std::size_t idx = 56; << 708 else if (theDirSecondPartDef->GetParticleName() == "e-") idx = 1; 779 if(theDirSecondPartDef->GetParticleName( << 709 else if (theDirSecondPartDef->GetParticleName() == "e+") idx = 2; 780 idx = 0; << 710 if (idx <56) { 781 else if(theDirSecondPartDef->GetParticle << 711 const std::vector<G4double>* aVec = G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(idx); 782 idx = 1; << 712 Tlow =(*aVec)[aCouple->GetIndex()]; 783 else if(theDirSecondPartDef->GetParticle << 713 } 784 idx = 2; << 714 785 if(idx < 56) << 715 786 { << 716 } 787 const std::vector<G4double>* aVec = << 717 if ( Ekin<=listOfAdjointEMModel[i]->GetHighEnergyLimit() && Ekin>=listOfAdjointEMModel[i]->GetLowEnergyLimit()){ 788 G4ProductionCutsTable::GetProduction << 718 if (aPartDef == listOfAdjointEMModel[i]->GetAdjointEquivalentOfDirectPrimaryParticleDefinition()){ 789 idx); << 719 CS=ComputeAdjointCS(currentMaterial, 790 Tlow = (*aVec)[aCouple->GetIndex()]; << 720 listOfAdjointEMModel[i], 791 } << 721 Ekin, Tlow,true,CS_Vs_Element); 792 } << 722 TotalCS += CS; 793 if(Ekin <= adjModel->GetHighEnergyLimit() << 723 (*listSigmaTableForAdjointModelScatProjToProj[i])[currentMatIndex]->PutValue(eindex,CS); 794 Ekin >= adjModel->GetLowEnergyLimit()) << 724 } 795 { << 725 if (aPartDef == listOfAdjointEMModel[i]->GetAdjointEquivalentOfDirectSecondaryParticleDefinition()){ 796 if(aPartDef == << 726 CS = ComputeAdjointCS(currentMaterial, 797 adjModel->GetAdjointEquivalentOfDirec << 727 listOfAdjointEMModel[i], 798 { << 728 Ekin, Tlow,false, CS_Vs_Element); 799 CS = ComputeAdjointCS(fCurrentMaterial << 729 TotalCS += CS; 800 CS_Vs_Element); << 730 (*listSigmaTableForAdjointModelProdToProj[i])[currentMatIndex]->PutValue(eindex,CS); 801 TotalCS += CS; << 731 } 802 (*fSigmaTableForAdjointModelScatProjTo << 732 803 ->PutValue(fNbins, CS); << 733 } 804 } << 734 else { 805 if(aPartDef == << 735 (*listSigmaTableForAdjointModelScatProjToProj[i])[currentMatIndex]->PutValue(eindex,0.); 806 adjModel->GetAdjointEquivalentOfDirec << 736 (*listSigmaTableForAdjointModelProdToProj[i])[currentMatIndex]->PutValue(eindex,0.); 807 { << 737 808 CS = ComputeAdjointCS(fCurrentMaterial << 738 } 809 CS_Vs_Element); << 739 } 810 TotalCS += CS; << 740 return TotalCS; 811 (*fSigmaTableForAdjointModelProdToProj << 741 812 fNbins, CS); << 742 813 } << 743 } 814 } << 815 else << 816 { << 817 (*fSigmaTableForAdjointModelScatProjToPr << 818 ->PutValue(fNbins, 0.); << 819 (*fSigmaTableForAdjointModelProdToProj[i << 820 fNbins, 0.); << 821 } << 822 } << 823 return TotalCS; << 824 } << 825 << 826 ////////////////////////////////////////////// 744 /////////////////////////////////////////////////////// >> 745 // 827 std::vector<G4AdjointCSMatrix*> 746 std::vector<G4AdjointCSMatrix*> 828 G4AdjointCSManager::BuildCrossSectionsModelAnd << 747 G4AdjointCSManager::BuildCrossSectionsMatricesForAGivenModelAndElement(G4VEmAdjointModel* aModel,G4int Z,G4int A, 829 << 748 G4int nbin_pro_decade) 830 << 749 { 831 { << 750 G4AdjointCSMatrix* theCSMatForProdToProjBackwardScattering = new G4AdjointCSMatrix(false); 832 G4AdjointCSMatrix* theCSMatForProdToProjBack << 751 G4AdjointCSMatrix* theCSMatForScatProjToProjBackwardScattering = new G4AdjointCSMatrix(true); 833 new G4AdjointCSMatrix(false); << 752 834 G4AdjointCSMatrix* theCSMatForScatProjToProj << 753 835 new G4AdjointCSMatrix(true); << 754 //make the vector of primary energy of the adjoint particle, could try to make this just once ? 836 << 755 837 // make the vector of primary energy of the << 756 G4double EkinMin =aModel->GetLowEnergyLimit(); 838 G4double EkinMin = aModel->GetLowEner << 757 G4double EkinMaxForScat =aModel->GetHighEnergyLimit()*0.999; 839 G4double EkinMaxForScat = aModel->GetHighEne << 758 G4double EkinMaxForProd =aModel->GetHighEnergyLimit()*0.999; 840 G4double EkinMaxForProd = aModel->GetHighEne << 759 if (aModel->GetSecondPartOfSameType() )EkinMaxForProd =EkinMaxForProd/2.; 841 if(aModel->GetSecondPartOfSameType()) << 760 842 EkinMaxForProd = EkinMaxForProd / 2.; << 761 843 << 762 //Product to projectile backward scattering 844 // Product to projectile backward scattering << 763 //----------------------------------------- 845 G4double dE = std::pow(10., 1. / nbin_pro_de << 764 G4double fE=std::pow(10.,1./nbin_pro_decade); 846 G4double E2 = << 765 G4double E2=std::pow(10.,double( int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE; 847 std::pow(10., double(int(std::log10(EkinMi << 766 G4double E1=EkinMin; 848 nbin_pro_decade) / dE; << 767 // Loop checking, 07-Aug-2015, Vladimir Ivanchenko 849 G4double E1 = EkinMin; << 768 while (E1 <EkinMaxForProd){ 850 // Loop checking, 07-Aug-2015, Vladimir Ivan << 769 E1=std::max(EkinMin,E2); 851 while(E1 < EkinMaxForProd) << 770 E1=std::min(EkinMaxForProd,E1); 852 { << 771 std::vector< std::vector< double>* > aMat= aModel->ComputeAdjointCrossSectionVectorPerAtomForSecond(E1,Z,A,nbin_pro_decade); 853 E1 = std::max(EkinMin, E2); << 772 if (aMat.size()>=2) { 854 E1 = std::min(EkinMaxForProd, E1); << 773 std::vector< double>* log_ESecVec=aMat[0]; 855 std::vector<std::vector<double>*> aMat = << 774 std::vector< double>* log_CSVec=aMat[1]; 856 aModel->ComputeAdjointCrossSectionVector << 775 G4double log_adjointCS=log_CSVec->back(); 857 << 776 //normalise CSVec such that it becomes a probability vector 858 if(aMat.size() >= 2) << 777 for (size_t j=0;j<log_CSVec->size();j++) { 859 { << 778 if (j==0) (*log_CSVec)[j] = 0.; 860 std::vector<double>* log_ESecVec = aMat[ << 779 else (*log_CSVec)[j]=std::log(1.-std::exp((*log_CSVec)[j]-log_adjointCS) +1e-50); 861 std::vector<double>* log_CSVec = aMat[ << 780 } 862 G4double log_adjointCS = log_C << 781 (*log_CSVec)[log_CSVec->size()-1]=(*log_CSVec)[log_CSVec->size()-2]-std::log(1000.); 863 // normalise CSVec such that it becomes << 782 theCSMatForProdToProjBackwardScattering->AddData(std::log(E1),log_adjointCS,log_ESecVec,log_CSVec,0); 864 for(std::size_t j = 0; j < log_CSVec->si << 783 } 865 { << 784 E1=E2; 866 if(j == 0) << 785 E2*=fE; 867 (*log_CSVec)[j] = 0.; << 786 } 868 else << 787 869 (*log_CSVec)[j] = << 788 //Scattered projectile to projectile backward scattering 870 std::log(1. - std::exp((*log_CSVec << 789 //----------------------------------------- 871 } << 790 872 (*log_CSVec)[log_CSVec->size() - 1] = << 791 E2=std::pow(10.,double( int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE; 873 (*log_CSVec)[log_CSVec->size() - 2] - << 792 E1=EkinMin; 874 theCSMatForProdToProjBackwardScattering- << 793 // Loop checking, 07-Aug-2015, Vladimir Ivanchenko 875 std::log(E1), log_adjointCS, log_ESecV << 794 while (E1 <EkinMaxForScat){ 876 } << 795 E1=std::max(EkinMin,E2); 877 E1 = E2; << 796 E1=std::min(EkinMaxForScat,E1); 878 E2 *= dE; << 797 std::vector< std::vector< double>* > aMat= aModel->ComputeAdjointCrossSectionVectorPerAtomForScatProj(E1,Z,A,nbin_pro_decade); 879 } << 798 if (aMat.size()>=2) { 880 << 799 std::vector< double>* log_ESecVec=aMat[0]; 881 // Scattered projectile to projectile backwa << 800 std::vector< double>* log_CSVec=aMat[1]; 882 E2 = std::pow(10., G4double(int(std::log10(E << 801 G4double log_adjointCS=log_CSVec->back(); 883 nbin_pro_decade) / dE; << 802 //normalise CSVec such that it becomes a probability vector 884 E1 = EkinMin; << 803 for (size_t j=0;j<log_CSVec->size();j++) { 885 // Loop checking, 07-Aug-2015, Vladimir Ivan << 804 if (j==0) (*log_CSVec)[j] = 0.; 886 while(E1 < EkinMaxForScat) << 805 else (*log_CSVec)[j]=std::log(1.-std::exp((*log_CSVec)[j]-log_adjointCS)+1e-50); 887 { << 806 } 888 E1 = std::max(EkinMin, E2); << 807 (*log_CSVec)[log_CSVec->size()-1]=(*log_CSVec)[log_CSVec->size()-2]-std::log(1000.); 889 E1 = std::min(EkinMaxForScat, E1); << 808 theCSMatForScatProjToProjBackwardScattering->AddData(std::log(E1),log_adjointCS,log_ESecVec,log_CSVec,0); 890 std::vector<std::vector<G4double>*> aMat = << 809 } 891 aModel->ComputeAdjointCrossSectionVector << 810 E1=E2; 892 E1, Z, A, nbin_pro_decade); << 811 E2*=fE; 893 if(aMat.size() >= 2) << 812 } 894 { << 813 895 std::vector<G4double>* log_ESecVec = aMa << 814 896 std::vector<G4double>* log_CSVec = aMa << 897 G4double log_adjointCS = log << 898 // normalise CSVec such that it becomes << 899 for(std::size_t j = 0; j < log_CSVec->si << 900 { << 901 if(j == 0) << 902 (*log_CSVec)[j] = 0.; << 903 else << 904 (*log_CSVec)[j] = << 905 std::log(1. - std::exp((*log_CSVec << 906 } << 907 (*log_CSVec)[log_CSVec->size() - 1] = << 908 (*log_CSVec)[log_CSVec->size() - 2] - << 909 theCSMatForScatProjToProjBackwardScatter << 910 std::log(E1), log_adjointCS, log_ESecV << 911 } << 912 E1 = E2; << 913 E2 *= dE; << 914 } << 915 << 916 std::vector<G4AdjointCSMatrix*> res; 815 std::vector<G4AdjointCSMatrix*> res; >> 816 res.clear(); 917 res.push_back(theCSMatForProdToProjBackwardS 817 res.push_back(theCSMatForProdToProjBackwardScattering); 918 res.push_back(theCSMatForScatProjToProjBackw 818 res.push_back(theCSMatForScatProjToProjBackwardScattering); >> 819 >> 820 >> 821 /* >> 822 G4String file_name; >> 823 std::stringstream astream; >> 824 G4String str_Z; >> 825 astream<<Z; >> 826 astream>>str_Z; >> 827 theCSMatForProdToProjBackwardScattering->Write(aModel->GetName()+G4String("_CSMat_Z")+str_Z+"_ProdToProj.txt"); >> 828 theCSMatForScatProjToProjBackwardScattering->Write(aModel->GetName()+G4String("_CSMat_Z")+str_Z+"_ScatProjToProj.txt"); >> 829 >> 830 */ 919 831 >> 832 920 return res; 833 return res; >> 834 >> 835 921 } 836 } 922 << 923 ////////////////////////////////////////////// 837 /////////////////////////////////////////////////////// >> 838 // 924 std::vector<G4AdjointCSMatrix*> 839 std::vector<G4AdjointCSMatrix*> 925 G4AdjointCSManager::BuildCrossSectionsModelAnd << 840 G4AdjointCSManager::BuildCrossSectionsMatricesForAGivenModelAndMaterial(G4VEmAdjointModel* aModel, 926 G4VEmAdjointModel* aModel, G4Material* aMate << 841 G4Material* aMaterial, 927 { << 842 G4int nbin_pro_decade) 928 G4AdjointCSMatrix* theCSMatForProdToProjBack << 843 { 929 new G4AdjointCSMatrix(false); << 844 G4AdjointCSMatrix* theCSMatForProdToProjBackwardScattering = new G4AdjointCSMatrix(false); 930 G4AdjointCSMatrix* theCSMatForScatProjToProj << 845 G4AdjointCSMatrix* theCSMatForScatProjToProjBackwardScattering = new G4AdjointCSMatrix(true); 931 new G4AdjointCSMatrix(true); << 846 932 << 847 933 G4double EkinMin = aModel->GetLowEner << 848 //make the vector of primary energy of the adjoint particle, could try to make this just once ? 934 G4double EkinMaxForScat = aModel->GetHighEne << 849 935 G4double EkinMaxForProd = aModel->GetHighEne << 850 G4double EkinMin =aModel->GetLowEnergyLimit(); 936 if(aModel->GetSecondPartOfSameType()) << 851 G4double EkinMaxForScat =aModel->GetHighEnergyLimit()*0.999; 937 EkinMaxForProd /= 2.; << 852 G4double EkinMaxForProd =aModel->GetHighEnergyLimit()*0.999; 938 << 853 if (aModel->GetSecondPartOfSameType() )EkinMaxForProd =EkinMaxForProd/2.; 939 // Product to projectile backward scattering << 854 940 G4double dE = std::pow(10., 1. / nbin_pro_de << 855 941 G4double E2 = << 856 942 std::pow(10., double(int(std::log10(EkinMi << 857 943 nbin_pro_decade) / dE; << 858 944 G4double E1 = EkinMin; << 859 945 // Loop checking, 07-Aug-2015, Vladimir Ivan << 860 946 while(E1 < EkinMaxForProd) << 861 //Product to projectile backward scattering 947 { << 862 //----------------------------------------- 948 E1 = std::max(EkinMin, E2); << 863 G4double fE=std::pow(10.,1./nbin_pro_decade); 949 E1 = std::min(EkinMaxForProd, E1); << 864 G4double E2=std::pow(10.,double( int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE; 950 std::vector<std::vector<G4double>*> aMat = << 865 G4double E1=EkinMin; 951 aModel->ComputeAdjointCrossSectionVector << 866 // Loop checking, 07-Aug-2015, Vladimir Ivanchenko 952 aMaterial, E1, nbin_pro_decade); << 867 while (E1 <EkinMaxForProd){ 953 if(aMat.size() >= 2) << 868 E1=std::max(EkinMin,E2); 954 { << 869 E1=std::min(EkinMaxForProd,E1); 955 std::vector<G4double>* log_ESecVec = aMa << 870 std::vector< std::vector< double>* > aMat= aModel->ComputeAdjointCrossSectionVectorPerVolumeForSecond(aMaterial,E1,nbin_pro_decade); 956 std::vector<G4double>* log_CSVec = aMa << 871 if (aMat.size()>=2) { 957 G4double log_adjointCS = log << 872 std::vector< double>* log_ESecVec=aMat[0]; 958 << 873 std::vector< double>* log_CSVec=aMat[1]; 959 // normalise CSVec such that it becomes << 874 G4double log_adjointCS=log_CSVec->back(); 960 for(std::size_t j = 0; j < log_CSVec->si << 875 961 { << 876 //normalise CSVec such that it becomes a probability vector 962 if(j == 0) << 877 for (size_t j=0;j<log_CSVec->size();j++) { 963 (*log_CSVec)[j] = 0.; << 878 //G4cout<<"CSMan1 "<<(*log_CSVec)[j]<<G4endl; 964 else << 879 if (j==0) (*log_CSVec)[j] = 0.; 965 (*log_CSVec)[j] = << 880 else (*log_CSVec)[j]=std::log(1.-std::exp((*log_CSVec)[j]-log_adjointCS)); 966 std::log(1. - std::exp((*log_CSVec << 881 //G4cout<<"CSMan2 "<<(*log_CSVec)[j]<<G4endl; 967 } << 882 } 968 (*log_CSVec)[log_CSVec->size() - 1] = << 883 (*log_CSVec)[log_CSVec->size()-1]=(*log_CSVec)[log_CSVec->size()-2]-std::log(1000.); 969 (*log_CSVec)[log_CSVec->size() - 2] - << 884 theCSMatForProdToProjBackwardScattering->AddData(std::log(E1),log_adjointCS,log_ESecVec,log_CSVec,0); 970 theCSMatForProdToProjBackwardScattering- << 885 } 971 std::log(E1), log_adjointCS, log_ESecV << 886 972 } << 887 973 << 888 974 E1 = E2; << 889 E1=E2; 975 E2 *= dE; << 890 E2*=fE; 976 } << 891 } 977 << 892 978 // Scattered projectile to projectile backwa << 893 //Scattered projectile to projectile backward scattering 979 E2 = << 894 //----------------------------------------- 980 std::pow(10., G4double(G4int(std::log10(Ek << 895 981 nbin_pro_decade) / << 896 E2=std::pow(10.,double( int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE; 982 dE; << 897 E1=EkinMin; 983 E1 = EkinMin; << 898 // Loop checking, 07-Aug-2015, Vladimir Ivanchenko 984 while(E1 < EkinMaxForScat) << 899 while (E1 <EkinMaxForScat){ 985 { << 900 E1=std::max(EkinMin,E2); 986 E1 = std::max(EkinMin, E2); << 901 E1=std::min(EkinMaxForScat,E1); 987 E1 = std::min(EkinMaxForScat, E1); << 902 std::vector< std::vector< double>* > aMat= aModel->ComputeAdjointCrossSectionVectorPerVolumeForScatProj(aMaterial,E1,nbin_pro_decade); 988 std::vector<std::vector<G4double>*> aMat = << 903 if (aMat.size()>=2) { 989 aModel->ComputeAdjointCrossSectionVector << 904 std::vector< double>* log_ESecVec=aMat[0]; 990 aMaterial, E1, nbin_pro_decade); << 905 std::vector< double>* log_CSVec=aMat[1]; 991 if(aMat.size() >= 2) << 906 G4double log_adjointCS=log_CSVec->back(); 992 { << 907 993 std::vector<G4double>* log_ESecVec = aMa << 908 for (size_t j=0;j<log_CSVec->size();j++) { 994 std::vector<G4double>* log_CSVec = aMa << 909 //G4cout<<"CSMan1 "<<(*log_CSVec)[j]<<G4endl; 995 G4double log_adjointCS = log << 910 if (j==0) (*log_CSVec)[j] = 0.; 996 << 911 else (*log_CSVec)[j]=std::log(1.-std::exp((*log_CSVec)[j]-log_adjointCS)); 997 for(std::size_t j = 0; j < log_CSVec->si << 912 //G4cout<<"CSMan2 "<<(*log_CSVec)[j]<<G4endl;if (theAdjPartDef->GetParticleName() == "adj_gamma") return G4Gamma::Gamma(); 998 { << 913 999 if(j == 0) << 914 } 1000 (*log_CSVec)[j] = 0.; << 915 (*log_CSVec)[log_CSVec->size()-1]=(*log_CSVec)[log_CSVec->size()-2]-std::log(1000.); 1001 else << 916 1002 (*log_CSVec)[j] = << 917 theCSMatForScatProjToProjBackwardScattering->AddData(std::log(E1),log_adjointCS,log_ESecVec,log_CSVec,0); 1003 std::log(1. - std::exp((*log_CSVe << 918 } 1004 } << 919 E1=E2; 1005 (*log_CSVec)[log_CSVec->size() - 1] = << 920 E2*=fE; 1006 (*log_CSVec)[log_CSVec->size() - 2] - << 921 } 1007 << 922 1008 theCSMatForScatProjToProjBackwardScatte << 923 1009 std::log(E1), log_adjointCS, log_ESec << 924 1010 } << 925 1011 E1 = E2; << 926 1012 E2 *= dE; << 927 1013 } << 928 1014 << 1015 std::vector<G4AdjointCSMatrix*> res; 929 std::vector<G4AdjointCSMatrix*> res; >> 930 res.clear(); >> 931 1016 res.push_back(theCSMatForProdToProjBackward 932 res.push_back(theCSMatForProdToProjBackwardScattering); 1017 res.push_back(theCSMatForScatProjToProjBack << 933 res.push_back(theCSMatForScatProjToProjBackwardScattering); >> 934 >> 935 /* >> 936 theCSMatForProdToProjBackwardScattering->Write(aModel->GetName()+"_CSMat_"+aMaterial->GetName()+"_ProdToProj.txt"); >> 937 theCSMatForScatProjToProjBackwardScattering->Write(aModel->GetName()+"_CSMat_"+aMaterial->GetName()+"_ScatProjToProj.txt"); >> 938 */ >> 939 1018 940 1019 return res; 941 return res; >> 942 >> 943 1020 } 944 } 1021 945 1022 ///////////////////////////////////////////// 946 /////////////////////////////////////////////////////// 1023 G4ParticleDefinition* G4AdjointCSManager::Get << 947 // 1024 G4ParticleDefinition* theFwdPartDef) << 948 G4ParticleDefinition* G4AdjointCSManager::GetAdjointParticleEquivalent(G4ParticleDefinition* theFwdPartDef) 1025 { 949 { 1026 if(theFwdPartDef->GetParticleName() == "e-" << 950 if (theFwdPartDef->GetParticleName() == "e-") return G4AdjointElectron::AdjointElectron(); 1027 return G4AdjointElectron::AdjointElectron << 951 else if (theFwdPartDef->GetParticleName() == "gamma") return G4AdjointGamma::AdjointGamma(); 1028 else if(theFwdPartDef->GetParticleName() == << 952 else if (theFwdPartDef->GetParticleName() == "proton") return G4AdjointProton::AdjointProton(); 1029 return G4AdjointGamma::AdjointGamma(); << 953 else if (theFwdPartDef ==theFwdIon) return theAdjIon; 1030 else if(theFwdPartDef->GetParticleName() == << 954 1031 return G4AdjointProton::AdjointProton(); << 955 return 0; 1032 else if(theFwdPartDef == fFwdIon) << 1033 return fAdjIon; << 1034 return nullptr; << 1035 } 956 } 1036 << 1037 ///////////////////////////////////////////// 957 /////////////////////////////////////////////////////// 1038 G4ParticleDefinition* G4AdjointCSManager::Get << 958 // 1039 G4ParticleDefinition* theAdjPartDef) << 959 G4ParticleDefinition* G4AdjointCSManager::GetForwardParticleEquivalent(G4ParticleDefinition* theAdjPartDef) 1040 { 960 { 1041 if(theAdjPartDef->GetParticleName() == "adj << 961 if (theAdjPartDef->GetParticleName() == "adj_e-") return G4Electron::Electron(); 1042 return G4Electron::Electron(); << 962 else if (theAdjPartDef->GetParticleName() == "adj_gamma") return G4Gamma::Gamma(); 1043 else if(theAdjPartDef->GetParticleName() == << 963 else if (theAdjPartDef->GetParticleName() == "adj_proton") return G4Proton::Proton(); 1044 return G4Gamma::Gamma(); << 964 else if (theAdjPartDef == theAdjIon) return theFwdIon; 1045 else if(theAdjPartDef->GetParticleName() == << 965 return 0; 1046 return G4Proton::Proton(); << 966 } 1047 else if(theAdjPartDef == fAdjIon) << 967 /////////////////////////////////////////////////////// 1048 return fFwdIon; << 968 // 1049 return nullptr; << 969 void G4AdjointCSManager::DefineCurrentMaterial(const G4MaterialCutsCouple* couple) >> 970 { >> 971 if(couple != currentCouple) { >> 972 currentCouple = const_cast<G4MaterialCutsCouple*> (couple); >> 973 currentMaterial = const_cast<G4Material*> (couple->GetMaterial()); >> 974 currentMatIndex = couple->GetIndex(); >> 975 lastPartDefForCS =0; >> 976 LastEkinForCS =0; >> 977 LastCSCorrectionFactor =1.; >> 978 } 1050 } 979 } 1051 980 1052 ///////////////////////////////////////////// 981 /////////////////////////////////////////////////////// 1053 void G4AdjointCSManager::DefineCurrentMateria << 982 // 1054 const G4MaterialCutsCouple* couple) << 983 void G4AdjointCSManager::DefineCurrentParticle(const G4ParticleDefinition* aPartDef) 1055 { << 984 { 1056 if(couple != fCurrentCouple) << 985 if(aPartDef != currentParticleDef) { 1057 { << 986 1058 fCurrentCouple = const_cast<G4Ma << 987 currentParticleDef= const_cast< G4ParticleDefinition* > (aPartDef); 1059 fCurrentMaterial = const_cast<G4Ma << 988 massRatio=1; 1060 fCurrentMatIndex = couple->GetInde << 989 if (aPartDef == theAdjIon) massRatio = proton_mass_c2/aPartDef->GetPDGMass(); 1061 fLastCSCorrectionFactor = 1.; << 990 currentParticleIndex=1000000; 1062 } << 991 for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){ >> 992 if (aPartDef == theListOfAdjointParticlesInAction[i]) currentParticleIndex=i; >> 993 } >> 994 >> 995 } 1063 } 996 } 1064 997 1065 ///////////////////////////////////////////// << 1066 void G4AdjointCSManager::DefineCurrentParticl << 1067 const G4ParticleDefinition* aPartDef) << 1068 { << 1069 << 1070 if(aPartDef != fCurrentParticleDef) << 1071 { << 1072 fCurrentParticleDef = const_cast<G4Partic << 1073 fMassRatio = 1.; << 1074 if(aPartDef == fAdjIon) << 1075 fMassRatio = proton_mass_c2 / aPartDef- << 1076 fCurrentParticleIndex = 1000000; << 1077 for(std::size_t i = 0; i < fAdjointPartic << 1078 { << 1079 if(aPartDef == fAdjointParticlesInActio << 1080 fCurrentParticleIndex = i; << 1081 } << 1082 } << 1083 } << 1084 << 1085 ///////////////////////////////////////////// << 1086 G4double G4AdjointCSManager::ComputeAdjointCS << 1087 G4double aPrimEnergy, G4AdjointCSMatrix* an << 1088 { << 1089 std::vector<double>* theLogPrimEnergyVector << 1090 anAdjointCSMatrix->GetLogPrimEnergyVector << 1091 if(theLogPrimEnergyVector->empty()) << 1092 { << 1093 G4cout << "No data are contained in the g << 1094 return 0.; << 1095 } << 1096 G4double log_Tcut = std::log(Tcut); << 1097 G4double log_E = std::log(aPrimEnergy); << 1098 << 1099 if(aPrimEnergy <= Tcut || log_E > theLogPri << 1100 return 0.; << 1101 998 1102 G4AdjointInterpolator* theInterpolator = G4 << 1103 999 1104 std::size_t ind = << 1000 ///////////////////////////////////////////////////////////////////////////////////////////////// 1105 theInterpolator->FindPositionForLogVector << 1001 // 1106 G4double aLogPrimEnergy1, aLogPrimEnergy2; << 1002 G4double G4AdjointCSManager::ComputeAdjointCS(G4double aPrimEnergy,G4AdjointCSMatrix* 1107 G4double aLogCS1, aLogCS2; << 1003 anAdjointCSMatrix,G4double Tcut) 1108 G4double log01, log02; << 1004 { 1109 std::vector<G4double>* aLogSecondEnergyVect << 1005 std::vector< double> *theLogPrimEnergyVector = anAdjointCSMatrix->GetLogPrimEnergyVector(); 1110 std::vector<G4double>* aLogSecondEnergyVect << 1006 if (theLogPrimEnergyVector->size() ==0){ 1111 std::vector<G4double>* aLogProbVector1 << 1007 G4cout<<"No data are contained in the given AdjointCSMatrix!"<<G4endl; 1112 std::vector<G4double>* aLogProbVector2 << 1008 G4cout<<"The s"<<G4endl; 1113 std::vector<std::size_t>* aLogProbVectorInd << 1009 return 0.; 1114 std::vector<std::size_t>* aLogProbVectorInd << 1010 1115 << 1116 anAdjointCSMatrix->GetData((G4int)ind, aLog << 1117 aLogSecondEnergy << 1118 aLogProbVectorIn << 1119 anAdjointCSMatrix->GetData(G4int(ind + 1), << 1120 aLogSecondEnergy << 1121 aLogProbVectorIn << 1122 if (! (aLogProbVector1 && aLogProbVector2 & << 1123 aLogSecondEnergyVector1 && aLogS << 1124 return 0.; << 1125 } 1011 } >> 1012 G4double log_Tcut = std::log(Tcut); >> 1013 G4double log_E =std::log(aPrimEnergy); >> 1014 >> 1015 if (aPrimEnergy <= Tcut || log_E > theLogPrimEnergyVector->back()) return 0.; >> 1016 >> 1017 1126 1018 1127 if(anAdjointCSMatrix->IsScatProjToProj()) << 1019 G4AdjointInterpolator* theInterpolator=G4AdjointInterpolator::GetInstance(); 1128 { // case where the Tcut plays a role << 1020 1129 G4double log_minimum_prob1, log_minimum_p << 1021 size_t ind =theInterpolator->FindPositionForLogVector(log_E,*theLogPrimEnergyVector); 1130 log_minimum_prob1 = theInterpolator->Inte << 1022 G4double aLogPrimEnergy1,aLogPrimEnergy2; 1131 log_Tcut, *aLogSecondEnergyVector1, *aL << 1023 G4double aLogCS1,aLogCS2; 1132 log_minimum_prob2 = theInterpolator->Inte << 1024 G4double log01,log02; 1133 log_Tcut, *aLogSecondEnergyVector2, *aL << 1025 std::vector< double>* aLogSecondEnergyVector1 =0; 1134 aLogCS1 += log_minimum_prob1; << 1026 std::vector< double>* aLogSecondEnergyVector2 =0; 1135 aLogCS2 += log_minimum_prob2; << 1027 std::vector< double>* aLogProbVector1=0; >> 1028 std::vector< double>* aLogProbVector2=0; >> 1029 std::vector< size_t>* aLogProbVectorIndex1=0; >> 1030 std::vector< size_t>* aLogProbVectorIndex2=0; >> 1031 >> 1032 >> 1033 anAdjointCSMatrix->GetData(ind, aLogPrimEnergy1,aLogCS1,log01, aLogSecondEnergyVector1,aLogProbVector1,aLogProbVectorIndex1); >> 1034 anAdjointCSMatrix->GetData(ind+1, aLogPrimEnergy2,aLogCS2,log02, aLogSecondEnergyVector2,aLogProbVector2,aLogProbVectorIndex2); >> 1035 if (anAdjointCSMatrix->IsScatProjToProjCase()){ //case where the Tcut plays a role >> 1036 G4double log_minimum_prob1, log_minimum_prob2; >> 1037 log_minimum_prob1=theInterpolator->InterpolateForLogVector(log_Tcut,*aLogSecondEnergyVector1,*aLogProbVector1); >> 1038 log_minimum_prob2=theInterpolator->InterpolateForLogVector(log_Tcut,*aLogSecondEnergyVector2,*aLogProbVector2); >> 1039 aLogCS1+= log_minimum_prob1; >> 1040 aLogCS2+= log_minimum_prob2; 1136 } 1041 } 1137 << 1042 1138 G4double log_adjointCS = theInterpolator->L << 1043 G4double log_adjointCS = theInterpolator->LinearInterpolation(log_E,aLogPrimEnergy1,aLogPrimEnergy2,aLogCS1,aLogCS2); 1139 log_E, aLogPrimEnergy1, aLogPrimEnergy2, << 1044 return std::exp(log_adjointCS); 1140 return std::exp(log_adjointCS); << 1045 1141 } << 1046 >> 1047 } 1142 1048