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