Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 // particle_hp -- source file 27 // J.P. Wellisch, Nov-1996 28 // A prototype of the low energy neutron trans 29 // 30 // 080612 Bug fix contribution from Benoit Pir 31 // 080709 Bug fix Sampling Legendre expansion 32 // 101110 Bug fix in MF=6, LAW=2 case; contrib 33 // 34 // P. Arce, June-2014 Conversion neutron_hp to 35 // 36 #include "G4ParticleHPDiscreteTwoBody.hh" 37 38 #include "G4Alpha.hh" 39 #include "G4Deuteron.hh" 40 #include "G4Electron.hh" 41 #include "G4Gamma.hh" 42 #include "G4He3.hh" 43 #include "G4Neutron.hh" 44 #include "G4ParticleHPLegendreStore.hh" 45 #include "G4ParticleHPVector.hh" 46 #include "G4ParticleHPManager.hh" 47 #include "G4Positron.hh" 48 #include "G4Proton.hh" 49 #include "G4Triton.hh" 50 51 G4ParticleHPDiscreteTwoBody::G4ParticleHPDiscr 52 { 53 if (G4ParticleHPManager::GetInstance()->GetP 54 bCheckDiffCoeffRepr = false; 55 } 56 57 G4ParticleHPDiscreteTwoBody::~G4ParticleHPDisc 58 { 59 delete[] theCoeff; 60 } 61 62 void G4ParticleHPDiscreteTwoBody::Init(std::is 63 { 64 aDataFile >> nEnergy; 65 theManager.Init(aDataFile); 66 const std::size_t tsize = nEnergy > 0 ? nEne 67 theCoeff = new G4ParticleHPLegendreTable[tsi 68 for (std::size_t i = 0; i < tsize; ++i) { 69 G4double energy; 70 G4int aRep, nCoeff; 71 aDataFile >> energy >> aRep >> nCoeff; 72 // G4cout << this << " " << i << " G4Parti 73 //<< " " << aRep << " " << nCoeff << G4end 74 energy *= CLHEP::eV; 75 G4int nPoints = nCoeff; 76 if (aRep > 0) nPoints *= 2; 77 theCoeff[i].Init(energy, nPoints - 1); 78 theCoeff[i].SetRepresentation(aRep); 79 for (G4int ii = 0; ii < nPoints; ++ii) { 80 G4double y; 81 aDataFile >> y; 82 theCoeff[i].SetCoeff(ii, y); 83 } 84 } 85 } 86 87 G4ReactionProduct* G4ParticleHPDiscreteTwoBody 88 89 { // Interpolation still only for the most us 90 auto result = new G4ReactionProduct; 91 auto Z = static_cast<G4int>(massCode / 1000) 92 auto A = static_cast<G4int>(massCode - 1000 93 94 if (massCode == 0) { 95 result->SetDefinition(G4Gamma::Gamma()); 96 } 97 else if (A == 0) { 98 result->SetDefinition(G4Electron::Electron 99 if (Z == 1) result->SetDefinition(G4Positr 100 } 101 else if (A == 1) { 102 result->SetDefinition(G4Neutron::Neutron() 103 if (Z == 1) result->SetDefinition(G4Proton 104 } 105 else if (A == 2) { 106 result->SetDefinition(G4Deuteron::Deuteron 107 } 108 else if (A == 3) { 109 result->SetDefinition(G4Triton::Triton()); 110 if (Z == 2) result->SetDefinition(G4He3::H 111 } 112 else if (A == 4) { 113 result->SetDefinition(G4Alpha::Alpha()); 114 if (Z != 2) throw G4HadronicException(__FI 115 } 116 else { 117 throw G4HadronicException(__FILE__, __LINE 118 "G4ParticleHPDis 119 } 120 121 // get cosine(theta) 122 G4int i(0), it(0); 123 G4double cosTh(0); 124 for (i = 0; i < nEnergy; i++) { 125 it = i; 126 if (theCoeff[i].GetEnergy() > anEnergy) br 127 } 128 if (it == 0 || it == nEnergy - 1) { 129 if (theCoeff[it].GetRepresentation() == 0) 130 // TK Legendre expansion 131 G4ParticleHPLegendreStore theStore(1); 132 theStore.SetCoeff(0, theCoeff); 133 theStore.SetManager(theManager); 134 // cosTh = theStore.SampleMax(anEnergy); 135 // 080612TK contribution from Benoit Pir 136 cosTh = theStore.SampleDiscreteTwoBody(a 137 } 138 else if (theCoeff[it].GetRepresentation() 139 { 140 G4ParticleHPVector theStore; 141 G4InterpolationManager aManager; 142 aManager.Init(LINLIN, theCoeff[it].GetNu 143 theStore.SetInterpolationManager(aManage 144 for (i = 0; i < theCoeff[it].GetNumberOf 145 // 101110 146 // theStore.SetX(i, theCoeff[it].GetCo 147 // theStore.SetY(i, theCoeff[it].GetCo 148 theStore.SetX(i / 2, theCoeff[it].GetC 149 theStore.SetY(i / 2, theCoeff[it].GetC 150 } 151 cosTh = theStore.Sample(); 152 } 153 else if (theCoeff[it].GetRepresentation() 154 { 155 G4ParticleHPVector theStore; 156 G4InterpolationManager aManager; 157 aManager.Init(LOGLIN, theCoeff[it].GetNu 158 theStore.SetInterpolationManager(aManage 159 for (i = 0; i < theCoeff[it].GetNumberOf 160 // 101110 161 // theStore.SetX(i, theCoeff[it].GetCo 162 // theStore.SetY(i, theCoeff[it].GetCo 163 theStore.SetX(i / 2, theCoeff[it].GetC 164 theStore.SetY(i / 2, theCoeff[it].GetC 165 } 166 cosTh = theStore.Sample(); 167 } 168 else { 169 throw G4HadronicException(__FILE__, __LI 170 "unknown repre 171 } 172 } 173 else { 174 if (!bCheckDiffCoeffRepr 175 || theCoeff[it].GetRepresentation() == 176 { 177 if (theCoeff[it].GetRepresentation() == 178 // TK Legendre expansion 179 G4ParticleHPLegendreStore theStore(2); 180 theStore.SetCoeff(0, &(theCoeff[it - 1 181 theStore.SetCoeff(1, &(theCoeff[it])); 182 G4InterpolationManager aManager; 183 aManager.Init(theManager.GetScheme(it) 184 theStore.SetManager(aManager); 185 // 080709 TKDB 186 cosTh = theStore.SampleDiscreteTwoBody 187 } 188 else if (theCoeff[it].GetRepresentation( 189 { 190 G4ParticleHPVector theBuff1; 191 G4InterpolationManager aManager1; 192 aManager1.Init(LINLIN, theCoeff[it - 1 193 theBuff1.SetInterpolationManager(aMana 194 for (i = 0; i < theCoeff[it - 1].GetNu 195 // 101110 196 theBuff1.SetX(i / 2, theCoeff[it - 1 197 theBuff1.SetY(i / 2, theCoeff[it - 1 198 } 199 G4ParticleHPVector theBuff2; 200 G4InterpolationManager aManager2; 201 aManager2.Init(LINLIN, theCoeff[it].Ge 202 theBuff2.SetInterpolationManager(aMana 203 for (i = 0; i < theCoeff[it].GetNumber 204 theBuff2.SetX(i / 2, theCoeff[it].Ge 205 theBuff2.SetY(i / 2, theCoeff[it].Ge 206 } 207 208 G4double x1 = theCoeff[it - 1].GetEner 209 G4double x2 = theCoeff[it].GetEnergy() 210 G4double x = anEnergy; 211 G4double y1, y2, y, mu; 212 213 G4ParticleHPVector theStore1; 214 theStore1.SetInterpolationManager(aMan 215 G4ParticleHPVector theStore2; 216 theStore2.SetInterpolationManager(aMan 217 G4ParticleHPVector theStore; 218 219 // for fixed mu get p1, p2 and interpo 220 for (i = 0; i < theBuff1.GetVectorLeng 221 mu = theBuff1.GetX(i); 222 y1 = theBuff1.GetY(i); 223 y2 = theBuff2.GetY(mu); 224 y = theInt.Interpolate(theManager.Ge 225 theStore1.SetData(i, mu, y); 226 } 227 for (i = 0; i < theBuff2.GetVectorLeng 228 mu = theBuff2.GetX(i); 229 y1 = theBuff2.GetY(i); 230 y2 = theBuff1.GetY(mu); 231 y = theInt.Interpolate(theManager.Ge 232 theStore2.SetData(i, mu, y); 233 } 234 theStore.Merge(&theStore1, &theStore2) 235 cosTh = theStore.Sample(); 236 } 237 else if (theCoeff[it].GetRepresentation( 238 { 239 G4ParticleHPVector theBuff1; 240 G4InterpolationManager aManager1; 241 aManager1.Init(LOGLIN, theCoeff[it - 1 242 theBuff1.SetInterpolationManager(aMana 243 for (i = 0; i < theCoeff[it - 1].GetNu 244 // 101110 245 theBuff1.SetX(i / 2, theCoeff[it - 1 246 theBuff1.SetY(i / 2, theCoeff[it - 1 247 } 248 249 G4ParticleHPVector theBuff2; 250 G4InterpolationManager aManager2; 251 aManager2.Init(LOGLIN, theCoeff[it].Ge 252 theBuff2.SetInterpolationManager(aMana 253 for (i = 0; i < theCoeff[it].GetNumber 254 // 101110 255 theBuff2.SetX(i / 2, theCoeff[it].Ge 256 theBuff2.SetY(i / 2, theCoeff[it].Ge 257 } 258 259 G4double x1 = theCoeff[it - 1].GetEner 260 G4double x2 = theCoeff[it].GetEnergy() 261 G4double x = anEnergy; 262 G4double y1, y2, y, mu; 263 264 G4ParticleHPVector theStore1; 265 theStore1.SetInterpolationManager(aMan 266 G4ParticleHPVector theStore2; 267 theStore2.SetInterpolationManager(aMan 268 G4ParticleHPVector theStore; 269 270 // for fixed mu get p1, p2 and interpo 271 for (i = 0; i < theBuff1.GetVectorLeng 272 mu = theBuff1.GetX(i); 273 y1 = theBuff1.GetY(i); 274 y2 = theBuff2.GetY(mu); 275 y = theInt.Interpolate(theManager.Ge 276 theStore1.SetData(i, mu, y); 277 } 278 for (i = 0; i < theBuff2.GetVectorLeng 279 mu = theBuff2.GetX(i); 280 y1 = theBuff2.GetY(i); 281 y2 = theBuff1.GetY(mu); 282 y = theInt.Interpolate(theManager.Ge 283 theStore2.SetData(i, mu, y); 284 } 285 theStore.Merge(&theStore1, &theStore2) 286 cosTh = theStore.Sample(); 287 } 288 else { 289 throw G4HadronicException(__FILE__, __ 290 "Two neighbo 291 } 292 } 293 else { 294 G4cout << " theCoeff[it].GetRepresent ME 295 << &theCoeff[it - 1] << G4endl; 296 G4cout << " theCoeff[it].GetRepresent " 297 << " != " << theCoeff[it - 1].Get 298 299 throw G4HadronicException(__FILE__, __LI 300 "unknown repre 301 } 302 } 303 304 // now get the energy from kinematics and Q- 305 306 // G4double restEnergy = anEnergy+GetQValue( 307 308 // assumed to be in CMS @@@@@@@@@@@@@@@@@ 309 310 // 080612TK contribution from Benoit Pirard 311 // G4double residualMass = GetTarget()->Ge 312 // - result->GetMass 313 // G4double kinE = restEnergy/(1+result->Get 314 G4double A1 = GetTarget()->GetMass() / GetPr 315 G4double A1prim = result->GetMass() / GetPro 316 // G4double E1 = (A1+1)*(A1+1)/A1/A1*an 317 // Bug fix Bugzilla #1815 318 G4double E1 = anEnergy; 319 G4double kinE = (A1 + 1 - A1prim) / (A1 + 1) 320 321 result->SetKineticEnergy(kinE); // non rela 322 if (cosTh > 1.0) { cosTh = 1.0; } 323 else if(cosTh < -1.0) { cosTh = -1.0; } 324 G4double phi = CLHEP::twopi * G4UniformRand( 325 G4double sinth = std::sqrt((1.0 + cosTh)*(1. 326 G4double mtot = result->GetTotalMomentum(); 327 G4ThreeVector tempVector(mtot * sinth * std: 328 mtot * cosTh); 329 result->SetMomentum(tempVector); 330 return result; 331 } 332