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 // particle_hp -- source file 26 // particle_hp -- source file 27 // J.P. Wellisch, Nov-1996 27 // J.P. Wellisch, Nov-1996 28 // A prototype of the low energy neutron trans 28 // A prototype of the low energy neutron transport model. 29 // 29 // 30 // T.Koi 02-08-06 Modified Harmonise to resolv << 30 // 02-08-06 Modified Harmonise to reslove cross section trouble at high-end. T. KOI 31 // 31 // 32 // P. Arce, June-2014 Conversion neutron_hp to 32 // P. Arce, June-2014 Conversion neutron_hp to particle_hp 33 // V. Ivanchenko July-2023 converted back capt << 34 // 33 // 35 #include "G4ParticleHPElementData.hh" 34 #include "G4ParticleHPElementData.hh" 36 35 37 G4ParticleHPElementData::G4ParticleHPElementDa << 36 G4ParticleHPElementData::G4ParticleHPElementData() 38 { << 37 { 39 precision = 0.02; << 38 precision = 0.02; 40 theFissionData = new G4ParticleHPVector; << 39 theFissionData = new G4ParticleHPVector; 41 theCaptureData = new G4ParticleHPVector; << 40 theCaptureData = new G4ParticleHPVector; 42 theElasticData = new G4ParticleHPVector; << 41 theElasticData = new G4ParticleHPVector; 43 theInelasticData = new G4ParticleHPVector; << 42 theInelasticData = new G4ParticleHPVector; 44 theIsotopeWiseData = nullptr; << 43 theIsotopeWiseData = 0; 45 theBuffer = nullptr; << 44 theBuffer = NULL; 46 } << 45 } 47 << 46 48 G4ParticleHPElementData::~G4ParticleHPElementD << 47 G4ParticleHPElementData::~G4ParticleHPElementData() 49 { << 48 { 50 delete theFissionData; << 49 delete theFissionData; 51 delete theCaptureData; << 50 delete theCaptureData; 52 delete theElasticData; << 51 delete theElasticData; 53 delete theInelasticData; << 52 delete theInelasticData; 54 delete[] theIsotopeWiseData; << 53 delete [] theIsotopeWiseData; 55 } << 56 << 57 void G4ParticleHPElementData::Init(G4Element* << 58 G4ParticleD << 59 const char* << 60 { << 61 auto nIso = (G4int)theElement->GetNumberOfIs << 62 auto Z = theElement->GetZasInt(); << 63 const std::size_t dsize = nIso > 0 ? nIso : << 64 theIsotopeWiseData = new G4ParticleHPIsoData << 65 << 66 for (G4int i1 = 0; i1 < nIso; ++i1) { << 67 G4int A = theElement->GetIsotope(i1)->GetN << 68 G4int M = theElement->GetIsotope(i1)->Getm << 69 G4double frac = theElement->GetRelativeAbu << 70 UpdateData(A, Z, M, i1, frac, projectile, << 71 } 54 } 72 //theElasticData->ThinOut(precision); << 55 73 //if (projectile == G4Neutron::Neutron()) th << 56 void G4ParticleHPElementData::Init(G4Element * theElement, G4ParticleDefinition* projectile, const char* dataDirVariable ) 74 << 57 { 75 //theCaptureData->ThinOut(precision); << 58 G4int count = theElement->GetNumberOfIsotopes(); 76 //theFissionData->ThinOut(precision); << 59 if(count == 0) count += 77 } << 60 theStableOnes.GetNumberOfIsotopes(static_cast<G4int>(theElement->GetZ())); 78 << 61 theIsotopeWiseData = new G4ParticleHPIsoData[count]; 79 void G4ParticleHPElementData::UpdateData(G4int << 62 // filename = ein data-set je isotope. 80 G4dou << 63 count = 0; 81 G4Par << 64 G4int nIso = theElement->GetNumberOfIsotopes(); 82 const << 65 G4int Z = static_cast<G4int> (theElement->GetZ()); 83 { << 66 // G4int i1; 84 // Reads in the Data, using G4ParticleHPIsoD << 67 if(nIso!=0) 85 theIsotopeWiseData[index].Init(A, Z, M, abun << 68 { 86 << 69 for (G4int i1=0; i1<nIso; i1++) 87 theBuffer = theIsotopeWiseData[index].MakeEl << 70 { 88 Harmonise(theElasticData, theBuffer); << 71 // G4cout <<" Init: normal case"<<G4endl; 89 delete theBuffer; << 72 G4int A = theElement->GetIsotope(i1)->GetN(); 90 << 73 G4int M = theElement->GetIsotope(i1)->Getm(); 91 theBuffer = theIsotopeWiseData[index].MakeIn << 74 G4double frac = theElement->GetRelativeAbundanceVector()[i1]/CLHEP::perCent; 92 Harmonise(theInelasticData, theBuffer); << 75 //UpdateData(A, Z, count++, frac); 93 delete theBuffer; << 76 UpdateData(A, Z, M, count++, frac, projectile, dataDirVariable); 94 << 77 } 95 theBuffer = theIsotopeWiseData[index].MakeCa << 78 }else{ 96 Harmonise(theCaptureData, theBuffer); << 79 // G4cout <<" Init: theStableOnes case: Z="<<Z<<G4endl; 97 delete theBuffer; << 80 G4int first = theStableOnes.GetFirstIsotope(Z); 98 << 81 // G4cout <<"first="<<first<<" "<<theStableOnes.GetNumberOfIsotopes(theElement->GetZ())<<G4endl; 99 theBuffer = theIsotopeWiseData[index].MakeFi << 82 for(G4int i1=0; 100 Harmonise(theFissionData, theBuffer); << 83 i1<theStableOnes.GetNumberOfIsotopes(static_cast<G4int>(theElement->GetZ()) ); 101 delete theBuffer; << 84 i1++) 102 } << 85 { 103 << 86 // G4cout <<" Init: theStableOnes in the loop"<<G4endl; 104 void G4ParticleHPElementData::Harmonise(G4Part << 87 G4int A = theStableOnes.GetIsotopeNucleonCount(first+i1); 105 { << 88 G4double frac = theStableOnes.GetAbundance(first+i1); 106 if (theNew == nullptr) { << 89 // G4cout <<" Init: theStableOnes in the loop: "<<A<<G4endl; 107 return; << 90 UpdateData(A, Z, count++, frac, projectile, dataDirVariable); 108 } << 91 } 109 G4int s_tmp = 0, n = 0, m_tmp = 0; << 110 auto theMerge = new G4ParticleHPVector(theSt << 111 while (theStore->GetEnergy(s_tmp) < theNew-> << 112 && s_tmp < theStore->GetVectorLength( << 113 { << 114 theMerge->SetData(m_tmp++, theStore->GetEn << 115 ++s_tmp; << 116 } << 117 G4ParticleHPVector* active = theStore; << 118 G4ParticleHPVector* passive = theNew; << 119 G4ParticleHPVector* tmp; << 120 G4int a = s_tmp, p = n, t; << 121 while (a < active->GetVectorLength() << 122 && p < passive->GetVectorLength()) / << 123 { << 124 if (active->GetEnergy(a) <= passive->GetEn << 125 theMerge->SetData(m_tmp, active->GetEner << 126 G4double x = theMerge->GetEnergy(m_tmp); << 127 G4double y = std::max(0., passive->GetXs << 128 theMerge->SetData(m_tmp, x, theMerge->Ge << 129 ++m_tmp; << 130 ++a; << 131 } << 132 else { << 133 tmp = active; << 134 t = a; << 135 active = passive; << 136 a = p; << 137 passive = tmp; << 138 p = t; << 139 } 92 } >> 93 theElasticData->ThinOut(precision); >> 94 if( projectile == G4Neutron::Neutron() ) theInelasticData->ThinOut(precision); >> 95 >> 96 theCaptureData->ThinOut(precision); >> 97 theFissionData->ThinOut(precision); 140 } 98 } 141 while (a != active->GetVectorLength()) // L << 99 >> 100 //void G4ParticleHPElementData::UpdateData(G4int A, G4int Z, G4int index, G4double abundance) >> 101 void G4ParticleHPElementData::UpdateData(G4int A, G4int Z, G4int M, G4int index, G4double abundance, G4ParticleDefinition* projectile, const char* dataDirVariable ) 142 { 102 { 143 theMerge->SetData(m_tmp++, active->GetEner << 103 //Reads in the Data, using G4ParticleHPIsoData[], and its Init 144 ++a; << 104 // G4cout << "entered: ElementWiseData::UpdateData"<<G4endl; >> 105 //theIsotopeWiseData[index].Init(A, Z, abundance); >> 106 theIsotopeWiseData[index].Init(A, Z, M, abundance,projectile, dataDirVariable); >> 107 // G4cout << "ElementWiseData::UpdateData Init finished"<<G4endl; >> 108 >> 109 theBuffer = theIsotopeWiseData[index].MakeElasticData(); >> 110 // G4cout << "ElementWiseData::UpdateData MakeElasticData finished: " >> 111 // <<theBuffer->GetVectorLength()<<G4endl; >> 112 Harmonise(theElasticData, theBuffer); >> 113 // G4cout << "ElementWiseData::UpdateData Harmonise finished: " >> 114 // <<theElasticData->GetVectorLength()<<G4endl; >> 115 delete theBuffer; >> 116 >> 117 theBuffer = theIsotopeWiseData[index].MakeInelasticData(); >> 118 // G4cout << "ElementWiseData::UpdateData MakeInelasticData finished: " >> 119 // <<theBuffer->GetVectorLength()<<G4endl; >> 120 Harmonise(theInelasticData, theBuffer); >> 121 // G4cout << "ElementWiseData::UpdateData Harmonise finished: " >> 122 // <<theInelasticData->GetVectorLength()<<G4endl; >> 123 delete theBuffer; >> 124 >> 125 theBuffer = theIsotopeWiseData[index].MakeCaptureData(); >> 126 // G4cout << "ElementWiseData::UpdateData MakeCaptureData finished: " >> 127 // <<theBuffer->GetVectorLength()<<G4endl; >> 128 Harmonise(theCaptureData, theBuffer); >> 129 // G4cout << "ElementWiseData::UpdateData Harmonise finished: " >> 130 // <<theCaptureData->GetVectorLength()<<G4endl; >> 131 delete theBuffer; >> 132 >> 133 theBuffer = theIsotopeWiseData[index].MakeFissionData(); >> 134 // G4cout << "ElementWiseData::UpdateData MakeFissionData finished: " >> 135 // <<theBuffer->GetVectorLength()<<G4endl; >> 136 Harmonise(theFissionData, theBuffer); >> 137 // G4cout << "ElementWiseData::UpdateData Harmonise finished: " >> 138 // <<theFissionData->GetVectorLength()<<G4endl; >> 139 delete theBuffer; >> 140 >> 141 // G4cout << "ElementWiseData::UpdateData finished"<endl; 145 } 142 } 146 while (p != passive->GetVectorLength()) // << 143 >> 144 void G4ParticleHPElementData::Harmonise(G4ParticleHPVector *& theStore, G4ParticleHPVector * theNew) 147 { 145 { 148 G4double x = passive->GetEnergy(p); << 146 if(theNew == 0) { return; } 149 G4double y = std::max(0., active->GetXsec( << 147 G4int s_tmp = 0, n=0, m_tmp=0; 150 theMerge->SetData(m_tmp++, x, passive->Get << 148 G4ParticleHPVector * theMerge = new G4ParticleHPVector(theStore->GetVectorLength()); 151 ++p; << 149 // G4cout << "Harmonise 1: "<<theStore->GetEnergy(s_tmp)<<" "<<theNew->GetEnergy(0)<<G4endl; >> 150 while ( theStore->GetEnergy(s_tmp)<theNew->GetEnergy(0)&&s_tmp<theStore->GetVectorLength() ) // Loop checking, 11.05.2015, T. Koi >> 151 { >> 152 theMerge->SetData(m_tmp++, theStore->GetEnergy(s_tmp), theStore->GetXsec(s_tmp)); >> 153 s_tmp++; >> 154 } >> 155 G4ParticleHPVector *active = theStore; >> 156 G4ParticleHPVector * passive = theNew; >> 157 G4ParticleHPVector * tmp; >> 158 G4int a = s_tmp, p = n, t; >> 159 // G4cout << "Harmonise 2: "<<active->GetVectorLength()<<" "<<passive->GetVectorLength()<<G4endl; >> 160 while (a<active->GetVectorLength()&&p<passive->GetVectorLength()) // Loop checking, 11.05.2015, T. Koi >> 161 { >> 162 if(active->GetEnergy(a) <= passive->GetEnergy(p)) >> 163 { >> 164 theMerge->SetData(m_tmp, active->GetEnergy(a), active->GetXsec(a)); >> 165 G4double x = theMerge->GetEnergy(m_tmp); >> 166 G4double y = std::max(0., passive->GetXsec(x)); >> 167 theMerge->SetData(m_tmp, x, theMerge->GetXsec(m_tmp)+y); >> 168 m_tmp++; >> 169 a++; >> 170 } else { >> 171 // G4cout << "swapping in Harmonise"<<G4endl; >> 172 tmp = active; t=a; >> 173 active = passive; a=p; >> 174 passive = tmp; p=t; >> 175 } >> 176 } >> 177 // G4cout << "Harmonise 3: "<< a <<" "<<active->GetVectorLength()<<" "<<m<<G4endl; >> 178 while (a!=active->GetVectorLength()) // Loop checking, 11.05.2015, T. Koi >> 179 { >> 180 theMerge->SetData(m_tmp++, active->GetEnergy(a), active->GetXsec(a)); >> 181 a++; >> 182 } >> 183 // G4cout << "Harmonise 4: "<< p <<" "<<passive->GetVectorLength()<<" "<<m<<G4endl; >> 184 while (p!=passive->GetVectorLength()) // Loop checking, 11.05.2015, T. Koi >> 185 { >> 186 // Modified by T. KOI >> 187 //theMerge->SetData(m++, passive->GetEnergy(p), passive->GetXsec(p)); >> 188 G4double x = passive->GetEnergy(p); >> 189 G4double y = std::max(0., active->GetXsec(x)); >> 190 theMerge->SetData(m_tmp++, x, passive->GetXsec(p)+y); >> 191 p++; >> 192 } >> 193 // G4cout <<"Harmonise 5: "<< theMerge->GetVectorLength() << " " << m << G4endl; >> 194 delete theStore; >> 195 theStore = theMerge; >> 196 // G4cout <<"Harmonise 6: "<< theStore->GetVectorLength() << " " << m << G4endl; 152 } 197 } 153 delete theStore; << 154 theStore = theMerge; << 155 } << 156 << 157 G4ParticleHPVector* << 158 G4ParticleHPElementData::MakePhysicsVector(G4E << 159 G4P << 160 G4P << 161 cha << 162 { << 163 Init(theElement, projectile, dataDirVariable << 164 return GetData(theSet); << 165 } << 166 198 167 G4ParticleHPVector* << 199 G4ParticleHPVector * G4ParticleHPElementData::MakePhysicsVector(G4Element * theElement, 168 G4ParticleHPElementData::MakePhysicsVector(G4E << 200 G4ParticleDefinition * projectile, 169 G4P << 201 G4ParticleHPFissionData* theSet, 170 G4N << 202 char* dataDirVariable) 171 cha << 203 { 172 { << 204 if(projectile != G4Neutron::Neutron()) throw G4HadronicException(__FILE__, __LINE__, "not a neutron"); 173 Init(theElement, projectile, dataDirVariable << 205 Init ( theElement, projectile, dataDirVariable ); 174 return GetData(theSet); << 206 return GetData(theSet); 175 } << 207 } 176 << 208 G4ParticleHPVector * G4ParticleHPElementData::MakePhysicsVector(G4Element * theElement, 177 G4ParticleHPVector* << 209 G4ParticleDefinition * projectile, 178 G4ParticleHPElementData::MakePhysicsVector(G4E << 210 G4ParticleHPCaptureData * theSet, 179 G4P << 211 char* dataDirVariable) 180 G4P << 212 { 181 cha << 213 if(projectile != G4Neutron::Neutron()) throw G4HadronicException(__FILE__, __LINE__, "not a neutron"); 182 { << 214 Init ( theElement, projectile, dataDirVariable ); 183 Init(theElement, projectile, dataDirVariable << 215 return GetData(theSet); 184 return GetData(theSet); << 216 } 185 } << 217 G4ParticleHPVector * G4ParticleHPElementData::MakePhysicsVector(G4Element * theElement, 186 G4ParticleHPVector* << 218 G4ParticleDefinition * projectile, 187 G4ParticleHPElementData::MakePhysicsVector(G4E << 219 G4ParticleHPElasticData * theSet, 188 G4P << 220 char* dataDirVariable) 189 G4P << 221 { 190 cha << 222 if(projectile != G4Neutron::Neutron()) throw G4HadronicException(__FILE__, __LINE__, "not a neutron"); 191 { << 223 Init ( theElement, projectile, dataDirVariable ); 192 Init(theElement, projectile, dataDirVariable << 224 return GetData(theSet); 193 return GetData(theSet); << 225 } 194 } << 226 G4ParticleHPVector * G4ParticleHPElementData::MakePhysicsVector(G4Element * theElement, >> 227 G4ParticleDefinition * projectile, >> 228 G4ParticleHPInelasticData * theSet, >> 229 char* dataDirVariable) >> 230 { >> 231 if(projectile != G4Neutron::Neutron()) throw G4HadronicException(__FILE__, __LINE__, "not a neutron"); >> 232 Init ( theElement, projectile, dataDirVariable ); >> 233 return GetData(theSet); >> 234 } 195 235