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 // neutron_hp -- source file 26 // neutron_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 // 25-08-06 New Final State type (refFlag==3 , << 30 // 25-08-06 New Final State type (refFlag==3 , Legendre (Low Energy) + Probability (High Energy) ) 31 // is added by T. KOI 31 // is added by T. KOI 32 // 080904 Add Protection for negative energy r << 32 // 080904 Add Protection for negative energy results in very low energy ( 1E-6 eV ) scattering by T. Koi 33 // Koi << 34 // 33 // 35 // P. Arce, June-2014 Conversion neutron_hp to 34 // P. Arce, June-2014 Conversion neutron_hp to particle_hp 36 // 35 // 37 #include "G4ParticleHPElasticFS.hh" 36 #include "G4ParticleHPElasticFS.hh" >> 37 #include "G4ParticleHPManager.hh" 38 38 39 #include "G4Alpha.hh" << 39 #include "G4PhysicalConstants.hh" >> 40 #include "G4SystemOfUnits.hh" >> 41 #include "G4ReactionProduct.hh" >> 42 #include "G4Nucleus.hh" >> 43 #include "G4Proton.hh" 40 #include "G4Deuteron.hh" 44 #include "G4Deuteron.hh" 41 #include "G4HadronicParameters.hh" << 45 #include "G4Triton.hh" 42 #include "G4IonTable.hh" << 46 #include "G4Alpha.hh" >> 47 #include "G4ThreeVector.hh" 43 #include "G4LorentzVector.hh" 48 #include "G4LorentzVector.hh" 44 #include "G4Nucleus.hh" << 49 #include "G4IonTable.hh" 45 #include "G4ParticleHPDataUsed.hh" 50 #include "G4ParticleHPDataUsed.hh" 46 #include "G4ParticleHPManager.hh" << 47 #include "G4PhysicalConstants.hh" << 48 #include "G4PhysicsModelCatalog.hh" << 49 #include "G4Pow.hh" 51 #include "G4Pow.hh" 50 #include "G4Proton.hh" << 51 #include "G4ReactionProduct.hh" << 52 #include "G4SystemOfUnits.hh" << 53 #include "G4ThreeVector.hh" << 54 #include "G4Triton.hh" << 55 << 56 #include "zlib.h" 52 #include "zlib.h" 57 53 58 G4ParticleHPElasticFS::G4ParticleHPElasticFS() << 54 void G4ParticleHPElasticFS::Init (G4double A, G4double Z, G4int M, G4String & dirName, G4String &, G4ParticleDefinition* ) 59 { << 55 { 60 svtEmax = 0.0; << 56 G4String tString = "/FS"; 61 dbrcEmax = 0.0; << 57 G4bool dbool; 62 dbrcEmin = 0.0; << 58 G4ParticleHPDataUsed aFile = theNames.GetName(static_cast<G4int>(A), static_cast<G4int>(Z), M, dirName, tString, dbool); 63 dbrcAmin = 0.0; << 59 G4String filename = aFile.GetName(); 64 dbrcUse = false; << 60 SetAZMs( A, Z, M, aFile ); 65 xsForDBRC = nullptr; << 61 //theBaseA = aFile.GetA(); 66 << 62 //theBaseZ = aFile.GetZ(); 67 secID = G4PhysicsModelCatalog::GetModelID("m << 63 if(!dbool) 68 << 64 { 69 hasXsec = false; << 65 hasAnyData = false; 70 theCoefficients = nullptr; << 66 hasFSData = false; 71 theProbArray = nullptr; << 67 hasXsec = false; 72 << 68 return; 73 repFlag = 0; << 74 tE_of_repFlag3 = 0.0; << 75 targetMass = 0.0; << 76 frameFlag = 0; << 77 } << 78 << 79 void G4ParticleHPElasticFS::Init(G4double A, G << 80 const G4Strin << 81 G4ParticleDef << 82 { << 83 G4String tString = "/FS"; << 84 G4bool dbool = true; << 85 SetA_Z(A, Z, M); << 86 const G4ParticleHPDataUsed& aFile = << 87 theNames.GetName(theBaseA, theBaseZ, M, di << 88 const G4String& filename = aFile.GetName(); << 89 SetAZMs(aFile); << 90 if (!dbool) { << 91 hasAnyData = false; << 92 hasFSData = false; << 93 hasXsec = false; << 94 return; << 95 } << 96 << 97 // 130205 For compressed data files << 98 std::istringstream theData(std::ios::in); << 99 G4ParticleHPManager::GetInstance()->GetDataS << 100 // 130205 END << 101 theData >> repFlag >> targetMass >> frameFla << 102 << 103 if (repFlag == 1) { << 104 G4int nEnergy; << 105 theData >> nEnergy; << 106 theCoefficients = new G4ParticleHPLegendre << 107 theCoefficients->InitInterpolation(theData << 108 G4double temp, energy; << 109 G4int tempdep, nLegendre; << 110 G4int i, ii; << 111 for (i = 0; i < nEnergy; i++) { << 112 theData >> temp >> energy >> tempdep >> << 113 energy *= eV; << 114 theCoefficients->Init(i, energy, nLegend << 115 theCoefficients->SetTemperature(i, temp) << 116 G4double coeff = 0; << 117 for (ii = 0; ii < nLegendre; ii++) { << 118 // load legendre coefficients. << 119 theData >> coeff; << 120 theCoefficients->SetCoeff(i, ii + 1, c << 121 } << 122 } 69 } 123 } << 70 //130205 For compressed data files 124 else if (repFlag == 2) { << 71 std::istringstream theData(std::ios::in); 125 G4int nEnergy; << 72 G4ParticleHPManager::GetInstance()->GetDataStream(filename,theData); 126 theData >> nEnergy; << 73 //130205 END 127 theProbArray = new G4ParticleHPPartial(nEn << 74 theData >> repFlag >> targetMass >> frameFlag; 128 theProbArray->InitInterpolation(theData); << 75 if(repFlag==1) 129 G4double temp, energy; << 76 { 130 G4int tempdep, nPoints; << 77 G4int nEnergy; 131 for (G4int i = 0; i < nEnergy; i++) { << 78 theData >> nEnergy; 132 theData >> temp >> energy >> tempdep >> << 79 theCoefficients = new G4ParticleHPLegendreStore(nEnergy); 133 energy *= eV; << 80 theCoefficients->InitInterpolation(theData); 134 theProbArray->InitInterpolation(i, theDa << 81 G4double temp, energy; 135 theProbArray->SetT(i, temp); << 82 G4int tempdep, nLegendre; 136 theProbArray->SetX(i, energy); << 83 G4int i, ii; 137 G4double prob, costh; << 84 for (i=0; i<nEnergy; i++) 138 for (G4int ii = 0; ii < nPoints; ii++) { << 85 { 139 // fill probability arrays. << 86 theData >> temp >> energy >> tempdep >> nLegendre; 140 theData >> costh >> prob; << 87 energy *=eV; 141 theProbArray->SetX(i, ii, costh); << 88 theCoefficients->Init(i, energy, nLegendre); 142 theProbArray->SetY(i, ii, prob); << 89 theCoefficients->SetTemperature(i, temp); >> 90 G4double coeff=0; >> 91 for(ii=0; ii<nLegendre; ii++) >> 92 { >> 93 // load legendre coefficients. >> 94 theData >> coeff; >> 95 theCoefficients->SetCoeff(i, ii+1, coeff); // @@@HPW@@@ >> 96 } 143 } 97 } 144 theProbArray->DoneSetXY(i); << 145 } 98 } 146 } << 99 else if (repFlag==2) 147 else if (repFlag == 3) { << 100 { 148 G4int nEnergy_Legendre; << 101 G4int nEnergy; 149 theData >> nEnergy_Legendre; << 102 theData >> nEnergy; 150 if (nEnergy_Legendre <= 0) { << 103 theProbArray = new G4ParticleHPPartial(nEnergy, nEnergy); 151 std::stringstream iss; << 104 theProbArray->InitInterpolation(theData); 152 iss << "G4ParticleHPElasticFS::Init Data << 105 G4double temp, energy; 153 iss << "Z, A and M of problematic file i << 106 G4int tempdep, nPoints; 154 << theNDLDataM << " respectively."; << 107 for(G4int i=0; i<nEnergy; i++) 155 throw G4HadronicException(__FILE__, __LI << 108 { 156 } << 109 theData >> temp >> energy >> tempdep >> nPoints; 157 theCoefficients = new G4ParticleHPLegendre << 110 energy *= eV; 158 theCoefficients->InitInterpolation(theData << 111 theProbArray->InitInterpolation(i, theData); 159 G4double temp, energy; << 112 theProbArray->SetT(i, temp); 160 G4int tempdep, nLegendre; << 113 theProbArray->SetX(i, energy); 161 << 114 G4double prob, costh; 162 for (G4int i = 0; i < nEnergy_Legendre; i+ << 115 for(G4int ii=0; ii<nPoints; ii++) 163 theData >> temp >> energy >> tempdep >> << 116 { 164 energy *= eV; << 117 // fill probability arrays. 165 theCoefficients->Init(i, energy, nLegend << 118 theData >> costh >> prob; 166 theCoefficients->SetTemperature(i, temp) << 119 theProbArray->SetX(i, ii, costh); 167 G4double coeff = 0; << 120 theProbArray->SetY(i, ii, prob); 168 for (G4int ii = 0; ii < nLegendre; ii++) << 121 } 169 // load legendre coefficients. << 122 theProbArray->DoneSetXY( i ); 170 theData >> coeff; << 171 theCoefficients->SetCoeff(i, ii + 1, c << 172 } 123 } 173 } 124 } 174 << 125 else if ( repFlag==3 ) 175 tE_of_repFlag3 = energy; << 126 { 176 << 127 G4int nEnergy_Legendre; 177 G4int nEnergy_Prob; << 128 theData >> nEnergy_Legendre; 178 theData >> nEnergy_Prob; << 129 if ( nEnergy_Legendre == 0 ) { 179 theProbArray = new G4ParticleHPPartial(nEn << 130 std::stringstream iss; 180 theProbArray->InitInterpolation(theData); << 131 iss << "G4ParticleHPElasticFS::Init Data Error repFlag is 3 but nEnergy_Legendre is 0."; 181 G4int nPoints; << 132 iss << "Z, A and M of problematic file is " << theNDLDataZ << ", " << theNDLDataA << " and " << theNDLDataM << " respectively."; 182 for (G4int i = 0; i < nEnergy_Prob; i++) { << 133 throw G4HadronicException(__FILE__, __LINE__, iss.str() ); 183 theData >> temp >> energy >> tempdep >> << 134 } 184 energy *= eV; << 135 theCoefficients = new G4ParticleHPLegendreStore( nEnergy_Legendre ); 185 << 136 theCoefficients->InitInterpolation( theData ); 186 // consistency check << 137 G4double temp, energy; 187 if (i == 0) << 138 G4int tempdep, nLegendre; 188 // if ( energy != tE_of_repFlag3 ) //1 << 139 //G4int i, ii; 189 if (std::abs(energy - tE_of_repFlag3) << 140 for ( G4int i = 0 ; i < nEnergy_Legendre ; i++ ) 190 G4cout << "Warning Transition Energy << 141 { 191 << 142 theData >> temp >> energy >> tempdep >> nLegendre; 192 theProbArray->InitInterpolation(i, theDa << 143 energy *=eV; 193 theProbArray->SetT(i, temp); << 144 theCoefficients->Init( i , energy , nLegendre ); 194 theProbArray->SetX(i, energy); << 145 theCoefficients->SetTemperature( i , temp ); 195 G4double prob, costh; << 146 G4double coeff = 0; 196 for (G4int ii = 0; ii < nPoints; ii++) { << 147 for (G4int ii = 0 ; ii < nLegendre ; ii++ ) 197 // fill probability arrays. << 148 { 198 theData >> costh >> prob; << 149 // load legendre coefficients. 199 theProbArray->SetX(i, ii, costh); << 150 theData >> coeff; 200 theProbArray->SetY(i, ii, prob); << 151 theCoefficients->SetCoeff(i, ii+1, coeff); // @@@HPW@@@ 201 } << 152 } 202 theProbArray->DoneSetXY(i); << 153 } >> 154 >> 155 tE_of_repFlag3 = energy; >> 156 >> 157 G4int nEnergy_Prob; >> 158 theData >> nEnergy_Prob; >> 159 theProbArray = new G4ParticleHPPartial( nEnergy_Prob , nEnergy_Prob ); >> 160 theProbArray->InitInterpolation( theData ); >> 161 G4int nPoints; >> 162 for ( G4int i=0 ; i < nEnergy_Prob ; i++ ) >> 163 { >> 164 theData >> temp >> energy >> tempdep >> nPoints; >> 165 >> 166 energy *= eV; >> 167 >> 168 // consistency check >> 169 if ( i == 0 ) >> 170 //if ( energy != tE_of_repFlag3 ) //110620TK This is too tight for 32bit machines >> 171 if ( std::abs( energy - tE_of_repFlag3 ) / tE_of_repFlag3 > 1.0e-15 ) >> 172 G4cout << "Warning Transition Energy of repFlag3 is not consistent." << G4endl; >> 173 >> 174 theProbArray->InitInterpolation( i , theData ); >> 175 theProbArray->SetT( i , temp ); >> 176 theProbArray->SetX( i , energy ); >> 177 G4double prob, costh; >> 178 for( G4int ii = 0 ; ii < nPoints ; ii++ ) >> 179 { >> 180 // fill probability arrays. >> 181 theData >> costh >> prob; >> 182 theProbArray->SetX( i , ii , costh ); >> 183 theProbArray->SetY( i , ii , prob ); >> 184 } >> 185 theProbArray->DoneSetXY( i ); >> 186 } 203 } 187 } >> 188 else if (repFlag==0) >> 189 { >> 190 theData >> frameFlag; >> 191 } >> 192 else >> 193 { >> 194 G4cout << "unusable number for repFlag: repFlag="<<repFlag<<G4endl; >> 195 throw G4HadronicException(__FILE__, __LINE__, "G4ParticleHPElasticFS::Init -- unusable number for repFlag"); >> 196 } >> 197 //130205 For compressed data files(theData changed from ifstream to istringstream) >> 198 //theData.close(); 204 } 199 } 205 else if (repFlag == 0) { << 200 G4HadFinalState * G4ParticleHPElasticFS::ApplyYourself(const G4HadProjectile & theTrack) 206 theData >> frameFlag; << 201 { 207 } << 202 // G4cout << "G4ParticleHPElasticFS::ApplyYourself+"<<G4endl; 208 else { << 203 if ( theResult.Get() == NULL ) theResult.Put( new G4HadFinalState ); 209 G4cout << "unusable number for repFlag: re << 204 theResult.Get()->Clear(); 210 throw G4HadronicException(__FILE__, __LINE << 205 G4double eKinetic = theTrack.GetKineticEnergy(); 211 "G4ParticleHPEla << 206 const G4HadProjectile *incidentParticle = &theTrack; 212 } << 207 G4ReactionProduct theNeutron( const_cast<G4ParticleDefinition *>(incidentParticle->GetDefinition() )); 213 // 130205 For compressed data files(theData << 208 theNeutron.SetMomentum( incidentParticle->Get4Momentum().vect() ); 214 // theData.close(); << 209 theNeutron.SetKineticEnergy( eKinetic ); 215 } << 210 // G4cout << "G4ParticleHPElasticFS::ApplyYourself++"<<eKinetic<<" "<<G4endl; 216 << 211 // G4cout << "CMSVALUES 0 "<<theNeutron.GetTotalMomentum()<<G4endl; 217 G4HadFinalState* G4ParticleHPElasticFS::ApplyY << 212 218 { << 213 G4ReactionProduct theTarget; 219 if (theResult.Get() == nullptr) theResult.Pu << 214 G4Nucleus aNucleus; 220 theResult.Get()->Clear(); << 215 G4ThreeVector neuVelo = (1./incidentParticle->GetDefinition()->GetPDGMass())*theNeutron.GetMomentum(); 221 G4double eKinetic = theTrack.GetKineticEnerg << 216 theTarget = aNucleus.GetBiasedThermalNucleus( targetMass, neuVelo, theTrack.GetMaterial()->GetTemperature()); 222 const G4HadProjectile* incidentParticle = &t << 217 //t theTarget.SetDefinition( G4IonTable::GetIonTable()->GetIon( G4int(theBaseZ), G4int(theBaseA) , 0.0 ) ); //TESTPHP 223 G4ReactionProduct theNeutron( << 218 // G4cout << "Nucleus-test"<<" "<<targetMass<<" "; 224 const_cast<G4ParticleDefinition*>(incident << 219 // G4cout << theTarget.GetMomentum().x()<<" "; 225 theNeutron.SetMomentum(incidentParticle->Get << 220 // G4cout << theTarget.GetMomentum().y()<<" "; 226 theNeutron.SetKineticEnergy(eKinetic); << 221 // G4cout << theTarget.GetMomentum().z()<<G4endl; 227 << 222 228 G4ThreeVector neuVelo = << 223 // neutron and target defined as reaction products. 229 (1. / incidentParticle->GetDefinition()->G << 224 230 G4ReactionProduct theTarget = << 225 // prepare lorentz-transformation to Lab. 231 GetBiasedThermalNucleus(targetMass, neuVel << 226 232 << 227 G4ThreeVector the3Neutron = theNeutron.GetMomentum(); 233 // Neutron and target defined as G4ReactionP << 228 G4double nEnergy = theNeutron.GetTotalEnergy(); 234 // Prepare Lorentz transformation to lab << 229 G4ThreeVector the3Target = theTarget.GetMomentum(); 235 << 230 // cout << "@@@" << the3Target<<G4endl; 236 G4ThreeVector the3Neutron = theNeutron.GetMo << 231 G4double tEnergy = theTarget.GetTotalEnergy(); 237 G4double nEnergy = theNeutron.GetTotalEnergy << 232 G4ReactionProduct theCMS; 238 G4ThreeVector the3Target = theTarget.GetMome << 233 G4double totE = nEnergy+tEnergy; 239 G4double tEnergy = theTarget.GetTotalEnergy( << 234 G4ThreeVector the3CMS = the3Target+the3Neutron; 240 G4ReactionProduct theCMS; << 235 theCMS.SetMomentum(the3CMS); 241 G4double totE = nEnergy + tEnergy; << 236 G4double cmsMom = std::sqrt(the3CMS*the3CMS); 242 G4ThreeVector the3CMS = the3Target + the3Neu << 237 G4double sqrts = std::sqrt((totE-cmsMom)*(totE+cmsMom)); 243 theCMS.SetMomentum(the3CMS); << 238 theCMS.SetMass(sqrts); 244 G4double cmsMom = std::sqrt(the3CMS * the3CM << 239 theCMS.SetTotalEnergy(totE); 245 G4double sqrts = std::sqrt((totE - cmsMom) * << 240 246 theCMS.SetMass(sqrts); << 241 // data come as fcn of n-energy in nuclear rest frame 247 theCMS.SetTotalEnergy(totE); << 242 G4ReactionProduct boosted; 248 << 243 boosted.Lorentz(theNeutron, theTarget); 249 // Data come as function of n-energy in nucl << 244 eKinetic = boosted.GetKineticEnergy(); // get kinetic energy for scattering 250 G4ReactionProduct boosted; << 245 G4double cosTh = -2; 251 boosted.Lorentz(theNeutron, theTarget); << 246 if(repFlag == 1) 252 eKinetic = boosted.GetKineticEnergy(); // g << 247 { 253 G4double cosTh = -2; << 254 << 255 if (repFlag == 1) { << 256 cosTh = theCoefficients->SampleElastic(eKi << 257 } << 258 else if (repFlag == 2) { << 259 cosTh = theProbArray->Sample(eKinetic); << 260 } << 261 else if (repFlag == 3) { << 262 if (eKinetic <= tE_of_repFlag3) { << 263 cosTh = theCoefficients->SampleElastic(e 248 cosTh = theCoefficients->SampleElastic(eKinetic); 264 } 249 } 265 else { << 250 >> 251 else if (repFlag==2) >> 252 { 266 cosTh = theProbArray->Sample(eKinetic); 253 cosTh = theProbArray->Sample(eKinetic); 267 } 254 } 268 } << 255 else if (repFlag==3) 269 else if (repFlag == 0) { << 256 { 270 cosTh = 2. * G4UniformRand() - 1.; << 257 if ( eKinetic <= tE_of_repFlag3 ) 271 } << 258 { 272 else { << 259 cosTh = theCoefficients->SampleElastic(eKinetic); 273 G4cout << "Unusable number for repFlag: re << 260 } 274 throw G4HadronicException(__FILE__, __LINE << 261 else 275 "G4ParticleHPEla << 262 { 276 } << 263 cosTh = theProbArray->Sample(eKinetic); 277 << 264 } 278 if (cosTh < -1.1) { << 279 return nullptr; << 280 } << 281 << 282 G4double phi = twopi * G4UniformRand(); << 283 G4double cosPhi = std::cos(phi); << 284 G4double sinPhi = std::sin(phi); << 285 G4double theta = std::acos(cosTh); << 286 G4double sinth = std::sin(theta); << 287 << 288 if (frameFlag == 1) { << 289 // Projectile scattering values cosTh are << 290 // In this frame, do relativistic calculat << 291 // target 4-momenta << 292 << 293 theNeutron.Lorentz(theNeutron, theTarget); << 294 G4double mN = theNeutron.GetMass(); << 295 G4double Pinit = theNeutron.GetTotalMoment << 296 G4double Einit = theNeutron.GetTotalEnergy << 297 G4double mT = theTarget.GetMass(); << 298 << 299 G4double ratio = mT / mN; << 300 G4double sqt = std::sqrt(ratio * ratio - 1 << 301 G4double beta = Pinit / (mT + Einit); // << 302 G4double denom = 1. - beta * beta * cosTh << 303 G4double term1 = cosTh * (Einit * ratio + << 304 G4double pN = beta * mN * (term1 + sqt) / << 305 << 306 // Get the scattered momentum and rotate i << 307 G4ThreeVector pDir = theNeutron.GetMomentu << 308 G4double px = pN * pDir.x(); << 309 G4double py = pN * pDir.y(); << 310 G4double pz = pN * pDir.z(); << 311 << 312 G4ThreeVector pcmRot; << 313 pcmRot.setX(px * cosTh * cosPhi - py * sin << 314 pcmRot.setY(px * cosTh * sinPhi + py * cos << 315 pcmRot.setZ(-px * sinth + pz * cosTh); << 316 theNeutron.SetMomentum(pcmRot); << 317 G4double eN = std::sqrt(pN * pN + mN * mN) << 318 theNeutron.SetTotalEnergy(eN); << 319 << 320 // Get the scattered target momentum << 321 G4ReactionProduct toLab(-1. * theTarget); << 322 theTarget.SetMomentum(pDir * Pinit - pcmRo << 323 G4double eT = Einit - eN + mT; << 324 theTarget.SetTotalEnergy(eT); << 325 << 326 // Now back to lab frame << 327 theNeutron.Lorentz(theNeutron, toLab); << 328 theTarget.Lorentz(theTarget, toLab); << 329 << 330 // 111005 Protection for not producing 0 k << 331 if (theNeutron.GetKineticEnergy() <= 0) << 332 theNeutron.SetTotalEnergy(theNeutron.Get << 333 * (1. + G4Pow: << 334 if (theTarget.GetKineticEnergy() <= 0) << 335 theTarget.SetTotalEnergy(theTarget.GetMa << 336 } << 337 else if (frameFlag == 2) { << 338 // Projectile scattering values cosTh take << 339 << 340 G4LorentzVector proj(nEnergy, the3Neutron) << 341 G4LorentzVector targ(tEnergy, the3Target); << 342 G4ThreeVector boostToCM = proj.findBoostTo << 343 proj.boost(boostToCM); << 344 targ.boost(boostToCM); << 345 << 346 // Rotate projectile and target momenta by << 347 // Note: at this point collision axis is n << 348 // momentum given target nucleus by << 349 G4double px = proj.px(); << 350 G4double py = proj.py(); << 351 G4double pz = proj.pz(); << 352 << 353 G4ThreeVector pcmRot; << 354 pcmRot.setX(px * cosTh * cosPhi - py * sin << 355 pcmRot.setY(px * cosTh * sinPhi + py * cos << 356 pcmRot.setZ(-px * sinth + pz * cosTh); << 357 proj.setVect(pcmRot); << 358 targ.setVect(-pcmRot); << 359 << 360 // Back to lab frame << 361 proj.boost(-boostToCM); << 362 targ.boost(-boostToCM); << 363 << 364 theNeutron.SetMomentum(proj.vect()); << 365 theNeutron.SetTotalEnergy(proj.e()); << 366 << 367 theTarget.SetMomentum(targ.vect()); << 368 theTarget.SetTotalEnergy(targ.e()); << 369 << 370 // 080904 Add Protection for very low ener << 371 if (theNeutron.GetKineticEnergy() <= 0) { << 372 theNeutron.SetTotalEnergy(theNeutron.Get << 373 * (1. + G4Pow: << 374 } << 375 << 376 // 080904 Add Protection for very low ener << 377 if (theTarget.GetKineticEnergy() <= 0) { << 378 theTarget.SetTotalEnergy(theTarget.GetMa << 379 } 265 } 380 } << 266 else if (repFlag==0) 381 else { << 267 { 382 G4cout << "Value of frameFlag (1=LAB, 2=CM << 268 cosTh = 2.*G4UniformRand()-1.; 383 throw G4HadronicException(__FILE__, __LINE << 269 } 384 "G4ParticleHPEla << 270 else 385 } << 271 { 386 << 272 G4cout << "unusable number for repFlag: repFlag="<<repFlag<<G4endl; 387 // Everything is now in the lab frame << 273 throw G4HadronicException(__FILE__, __LINE__, "G4ParticleHPElasticFS::Init -- unusable number for repFlag"); 388 // Set energy change and momentum change << 274 } 389 theResult.Get()->SetEnergyChange(theNeutron. << 275 if(cosTh<-1.1) { return 0; } 390 theResult.Get()->SetMomentumChange(theNeutro << 276 G4double phi = twopi*G4UniformRand(); 391 << 277 G4double theta = std::acos(cosTh); 392 // Make recoil a G4DynamicParticle << 278 G4double sinth = std::sin(theta); 393 auto theRecoil = new G4DynamicParticle; << 279 if (frameFlag == 1) // final state data given in target rest frame. 394 theRecoil->SetDefinition(G4IonTable::GetIonT << 280 { 395 << 281 // we have the scattering angle, now we need the energy, then do the 396 theRecoil->SetMomentum(theTarget.GetMomentum << 282 // boosting. 397 theResult.Get()->AddSecondary(theRecoil, sec << 283 // relativistic elastic scattering energy angular correlation: 398 << 284 theNeutron.Lorentz(theNeutron, theTarget); 399 // Postpone the tracking of the primary neut << 285 G4double e0 = theNeutron.GetTotalEnergy(); 400 theResult.Get()->SetStatusChange(suspend); << 286 G4double p0 = theNeutron.GetTotalMomentum(); 401 return theResult.Get(); << 287 G4double mN = theNeutron.GetMass(); 402 } << 288 G4double mT = theTarget.GetMass(); 403 << 289 G4double eE = e0+mT; 404 void G4ParticleHPElasticFS::InitializeScatteri << 290 G4double ap = (mT+eE)*(mT-eE) + (p0+mN)*(p0-mN); >> 291 G4double a = 4*(eE+p0*cosTh)*(eE-p0*cosTh); >> 292 G4double b = 4*ap*p0*cosTh; >> 293 G4double c = (2.*eE*mN-ap)*(2.*eE*mN+ap); >> 294 G4double en = (-b+std::sqrt(b*b - 4*a*c) )/(2*a); >> 295 G4ThreeVector tempVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) ); >> 296 theNeutron.SetMomentum(tempVector); >> 297 theNeutron.SetTotalEnergy(std::sqrt(en*en+theNeutron.GetMass()*theNeutron.GetMass())); >> 298 // first to lab >> 299 theNeutron.Lorentz(theNeutron, -1.*theTarget); >> 300 // now to CMS >> 301 theNeutron.Lorentz(theNeutron, theCMS); >> 302 theTarget.SetMomentum(-theNeutron.GetMomentum()); >> 303 theTarget.SetTotalEnergy(theNeutron.GetTotalEnergy()); >> 304 // and back to lab >> 305 theNeutron.Lorentz(theNeutron, -1.*theCMS); >> 306 theTarget.Lorentz(theTarget, -1.*theCMS); >> 307 //111005 Protection for not producing 0 kinetic energy target >> 308 if ( theNeutron.GetKineticEnergy() <= 0 ) theNeutron.SetTotalEnergy ( theNeutron.GetMass() * ( 1 + G4Pow::GetInstance()->powA( 10 , -15.65 ) ) ); >> 309 if ( theTarget.GetKineticEnergy() <= 0 ) theTarget.SetTotalEnergy ( theTarget.GetMass() * ( 1 + G4Pow::GetInstance()->powA( 10 , -15.65 ) ) ); >> 310 } >> 311 else if (frameFlag == 2) // CMS >> 312 { >> 313 theNeutron.Lorentz(theNeutron, theCMS); >> 314 theTarget.Lorentz(theTarget, theCMS); >> 315 G4double en = theNeutron.GetTotalMomentum(); // already in CMS. >> 316 G4ThreeVector cmsMom_tmp=theNeutron.GetMomentum(); // for neutron direction in CMS >> 317 G4double cms_theta=cmsMom_tmp.theta(); >> 318 G4double cms_phi=cmsMom_tmp.phi(); >> 319 G4ThreeVector tempVector; >> 320 tempVector.setX(std::cos(theta)*std::sin(cms_theta)*std::cos(cms_phi) >> 321 +std::sin(theta)*std::cos(phi)*std::cos(cms_theta)*std::cos(cms_phi) >> 322 -std::sin(theta)*std::sin(phi)*std::sin(cms_phi) ); >> 323 tempVector.setY(std::cos(theta)*std::sin(cms_theta)*std::sin(cms_phi) >> 324 +std::sin(theta)*std::cos(phi)*std::cos(cms_theta)*std::sin(cms_phi) >> 325 +std::sin(theta)*std::sin(phi)*std::cos(cms_phi) ); >> 326 tempVector.setZ(std::cos(theta)*std::cos(cms_theta) >> 327 -std::sin(theta)*std::cos(phi)*std::sin(cms_theta) ); >> 328 tempVector *= en; >> 329 theNeutron.SetMomentum(tempVector); >> 330 theTarget.SetMomentum(-tempVector); >> 331 G4double tP = theTarget.GetTotalMomentum(); >> 332 G4double tM = theTarget.GetMass(); >> 333 theTarget.SetTotalEnergy(std::sqrt((tP+tM)*(tP+tM)-2.*tP*tM)); >> 334 >> 335 /* >> 336 For debug purpose. >> 337 Same transformation G4ReactionProduct.Lorentz() by 4vectors 405 { 338 { 406 // Initialize DBRC variables << 339 G4LorentzVector n4p = G4LorentzVector ( theNeutron.GetMomentum() , theNeutron.GetKineticEnergy() + theNeutron.GetMass() ); 407 svtEmax = G4HadronicParameters::Instance()-> << 340 G4cout << "before " << ( n4p.e() - n4p.m() ) / eV<< G4endl; 408 G4ParticleHPManager* manager = G4ParticleHPM << 341 G4LorentzVector cm4p = G4LorentzVector ( theCMS.GetMomentum() , theCMS.GetKineticEnergy() + theCMS.GetMass() ); 409 dbrcUse = manager->GetUseDBRC(); << 342 n4p.boost( cm4p.boostVector() ); 410 dbrcEmax = manager->GetMaxEnergyDBRC(); << 343 G4cout << cm4p/eV << G4endl; 411 dbrcEmin = manager->GetMinEnergyDBRC(); << 344 G4cout << "after " << ( n4p.e() - n4p.m() ) / eV<< G4endl; 412 dbrcAmin = manager->GetMinADBRC(); << 413 } 345 } >> 346 */ 414 347 415 G4ReactionProduct G4ParticleHPElasticFS::GetBi << 348 theNeutron.Lorentz(theNeutron, -1.*theCMS); 416 << 349 //080904 Add Protection for very low energy (1e-6eV) scattering 417 << 350 if ( theNeutron.GetKineticEnergy() <= 0 ) 418 { << 351 { 419 // This new method implements the DBRC (Dopp << 352 //theNeutron.SetMomentum( G4ThreeVector(0) ); 420 // on top of the SVT (Sampling of the Veloci << 353 //theNeutron.SetTotalEnergy ( theNeutron.GetMass() ); 421 // The SVT algorithm was written by Loic Thu << 354 //110822 Protection for not producing 0 kinetic energy neutron 422 // the method G4Nucleus::GetBiasedThermalNuc << 355 theNeutron.SetTotalEnergy ( theNeutron.GetMass() * ( 1 + G4Pow::GetInstance()->powA( 10 , -15.65 ) ) ); 423 // implemented the DBRC algorithm on top of << 424 // While the SVT algorithm is still present << 425 // the DBRC algorithm on top of the SVT one << 426 // order to avoid a cycle dependency between << 427 << 428 InitializeScatteringKernelParameters(); << 429 << 430 // Set threshold for SVT algorithm << 431 G4double E_threshold = svtEmax; << 432 if (svtEmax == -1.) { << 433 // If E_neutron <= 400*kB*T (400 is a comm << 434 // then apply the Sampling ot the Velocity << 435 // else consider the target nucleus being << 436 E_threshold = 400.0 * 8.617333262E-11 * te << 437 } << 438 << 439 // If DBRC is enabled and the nucleus is hea << 440 if (dbrcUse && aMass >= dbrcAmin) { << 441 E_threshold = std::max(svtEmax, dbrcEmax); << 442 } << 443 << 444 G4double E_neutron = 0.5 * aVelocity.mag2() << 445 << 446 G4bool dbrcIsOn = dbrcUse && E_neutron >= db << 447 << 448 G4Nucleus aNucleus; << 449 if (E_neutron > E_threshold || !dbrcIsOn) { << 450 // Apply only the SVT algorithm, not the D << 451 return aNucleus.GetBiasedThermalNucleus(ta << 452 } << 453 << 454 G4ReactionProduct result; << 455 result.SetMass(aMass * G4Neutron::Neutron()- << 456 << 457 // Beta = sqrt(m/2kT) << 458 G4double beta = << 459 std::sqrt(result.GetMass() << 460 / (2. * 8.617333262E-11 * temp)) << 461 << 462 // Neutron speed vn << 463 G4double vN_norm = aVelocity.mag(); << 464 G4double vN_norm2 = vN_norm * vN_norm; << 465 G4double y = beta * vN_norm; << 466 << 467 // Normalize neutron velocity << 468 aVelocity = (1. / vN_norm) * aVelocity; << 469 << 470 // Variables for sampling of target speed an << 471 G4double x2; << 472 G4double randThresholdSVT; << 473 G4double vT_norm, vT_norm2, mu; << 474 G4double acceptThresholdSVT; << 475 G4double vRelativeSpeed; << 476 G4double cdf0 = 2. / (2. + std::sqrt(CLHEP:: << 477 << 478 // DBRC variables << 479 G4double xsRelative = -99.; << 480 G4double randThresholdDBRC = 0.; << 481 // Calculate max cross-section in interval f << 482 G4double eMin = << 483 0.5 * G4Neutron::Neutron()->GetPDGMass() * << 484 G4double eMax = << 485 0.5 * G4Neutron::Neutron()->GetPDGMass() * << 486 G4double xsMax = xsForDBRC->GetMaxY(eMin, eM << 487 << 488 do { << 489 do { << 490 // Sample the target velocity vT in the << 491 if (G4UniformRand() < cdf0) { << 492 // Sample in C45 from https://laws.lan << 493 x2 = -std::log(G4UniformRand() * G4Uni << 494 } 356 } 495 else { << 496 // Sample in C61 from https://laws.lan << 497 G4double ampl = std::cos(CLHEP::pi / 2 << 498 x2 = -std::log(G4UniformRand()) - std: << 499 } << 500 << 501 vT_norm = std::sqrt(x2) / beta; << 502 vT_norm2 = vT_norm * vT_norm; << 503 357 504 // Sample cosine between the incident ne << 358 theTarget.Lorentz(theTarget, -1.*theCMS); 505 mu = 2. * G4UniformRand() - 1.; << 359 //080904 Add Protection for very low energy (1e-6eV) scattering 506 << 360 if ( theTarget.GetKineticEnergy() < 0 ) 507 // Define acceptance threshold for SVT << 361 { 508 vRelativeSpeed = std::sqrt(vN_norm2 + vT << 362 //theTarget.SetMomentum( G4ThreeVector(0) ); 509 acceptThresholdSVT = vRelativeSpeed / (v << 363 //theTarget.SetTotalEnergy ( theTarget.GetMass() ); 510 randThresholdSVT = G4UniformRand(); << 364 //110822 Protection for not producing 0 kinetic energy target 511 } while (randThresholdSVT >= acceptThresho << 365 theTarget.SetTotalEnergy ( theTarget.GetMass() * ( 1 + G4Pow::GetInstance()->powA( 10 , -15.65 ) ) ); 512 << 366 } 513 // Apply DBRC rejection << 367 } 514 xsRelative = xsForDBRC->GetXsec(0.5 * G4Ne << 368 else 515 * vRelativ << 369 { 516 randThresholdDBRC = G4UniformRand(); << 370 G4cout <<"Value of frameFlag (1=LAB, 2=CMS): "<<frameFlag; 517 << 371 throw G4HadronicException(__FILE__, __LINE__, "G4ParticleHPElasticFS::ApplyYourSelf frameflag incorrect"); 518 } while (randThresholdDBRC >= xsRelative / x << 372 } 519 << 373 // now all in Lab 520 aNucleus.DoKinematicsOfThermalNucleus(mu, vT << 374 // nun den recoil generieren...und energy change, momentum change angeben. 521 << 375 theResult.Get()->SetEnergyChange(theNeutron.GetKineticEnergy()); 522 return result; << 376 theResult.Get()->SetMomentumChange(theNeutron.GetMomentum().unit()); 523 } << 377 G4DynamicParticle* theRecoil = new G4DynamicParticle; >> 378 theRecoil->SetDefinition( G4IonTable::GetIonTable()->GetIon(static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA), 0 ) ); >> 379 theRecoil->SetMomentum(theTarget.GetMomentum()); >> 380 theResult.Get()->AddSecondary(theRecoil); >> 381 // G4cout << "G4ParticleHPElasticFS::ApplyYourself 10+"<<G4endl; >> 382 // postpone the tracking of the primary neutron >> 383 theResult.Get()->SetStatusChange(suspend); >> 384 return theResult.Get(); >> 385 } 524 386