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 #include "G4NRESP71M03.hh" 27 28 #include "G4Alpha.hh" 29 #include "G4Geantino.hh" 30 #include "G4IonTable.hh" 31 #include "G4Neutron.hh" 32 #include "G4PhysicalConstants.hh" 33 #include "G4RotationMatrix.hh" 34 #include "G4SystemOfUnits.hh" 35 #include "Randomize.hh" 36 37 // Energies for which angular distributions ar 38 const G4double G4NRESP71M03::BEN2[ndist] = { 39 5700., 8000., 8640., 8990., 9220., 9410 40 11870., 12140., 12320., 12570., 12940., 1342 41 14820., 15050., 15660., 15980., 16470., 1694 42 // Angular distribution of alpha particles for 43 const G4double G4NRESP71M03::B2[ndist][nrhos] 44 {0.000, 2839., 4027., 4951., 5736., 643 45 9767., 10241., 10703., 11152., 11595., 120 46 14507., 14909., 15308., 15711., 16110., 165 47 18965., 19393., 19823., 20266., 20715., 211 48 24344., 24981., 25682., 26467., 27391., 285 49 {0.000, 1771., 2528., 3129., 3653., 412 50 6578., 6961., 7345., 7725., 8108., 849 51 10891., 11306., 11721., 12142., 12563., 129 52 15585., 16034., 16490., 16958., 17438., 179 53 21595., 22424., 23373., 24482., 25820., 275 54 {0.000, 3518., 4863., 5805., 6540., 714 55 9588., 9896., 10191., 10470., 10744., 110 56 12456., 12685., 12915., 13144., 13373., 136 57 15032., 15286., 15547., 15821., 16100., 163 58 18711., 19283., 19980., 20897., 22251., 245 59 {0.000, 2287., 3364., 4319., 5277., 632 60 12305., 12839., 13317., 13757., 14168., 145 61 16681., 17024., 17372., 17724., 18086., 184 62 20976., 21463., 21965., 22481., 23005., 235 63 26568., 27080., 27614., 28183., 28820., 296 64 {0.000, 1690., 2456., 3097., 3697., 429 65 10722., 11906., 12795., 13515., 14134., 146 66 17313., 17699., 18082., 18463., 18849., 192 67 21912., 22446., 23002., 23571., 24133., 246 68 27451., 27878., 28321., 28789., 29314., 299 69 {0.000, 1674., 2434., 3078., 3685., 429 70 11014., 12069., 12918., 13640., 14278., 148 71 17574., 17966., 18353., 18733., 19116., 195 72 22085., 22591., 23122., 23665., 24212., 247 73 27492., 27919., 28355., 28820., 29336., 299 74 {0.000, 2302., 3245., 3967., 4577., 511 75 7838., 8278., 8739., 9229., 9767., 103 76 16198., 16688., 17121., 17514., 17881., 182 77 20269., 20646., 21051., 21497., 22009., 226 78 27146., 27677., 28189., 28710., 29267., 299 79 {0.000, 2195., 2934., 3458., 3879., 424 80 5899., 6138., 6371., 6600., 6832., 706 81 8529., 8809., 9104., 9424., 9773., 101 82 15060., 16801., 18312., 19352., 20140., 207 83 24067., 24674., 25355., 26153., 27137., 284 84 {0.000, 1712., 2406., 2937., 3380., 377 85 5689., 5978., 6270., 6565., 6864., 717 86 9484., 10034., 10675., 11422., 12258., 131 87 19138., 19971., 20696., 21353., 21965., 225 88 26068., 26684., 27316., 27972., 28689., 295 89 {0.000, 1853., 2651., 3289., 3851., 436 90 7172., 7646., 8139., 8645., 9179., 973 91 13621., 14287., 14928., 15544., 16128., 166 92 19682., 20143., 20602., 21061., 21523., 219 93 25013., 25603., 26247., 26964., 27803., 288 94 {0.000, 1959., 2770., 3397., 3934., 441 95 6948., 7382., 7840., 8333., 8876., 948 96 13398., 13871., 14309., 14725., 15125., 155 97 17992., 18473., 18982., 19519., 20077., 206 98 24069., 24699., 25391., 26185., 27146., 284 99 {0.000, 2290., 3391., 4447., 5824., 835 100 11630., 11935., 12227., 12510., 12789., 130 101 14875., 15239., 15637., 16077., 16565., 171 102 20413., 20927., 21449., 21989., 22558., 231 103 26873., 27410., 27944., 28493., 29091., 298 104 {0.000, 2119., 3072., 3874., 4627., 537 105 9921., 10501., 11034., 11531., 12002., 124 106 14888., 15266., 15641., 16013., 16386., 167 107 19312., 19849., 20438., 21072., 21724., 223 108 25491., 25987., 26514., 27098., 27791., 287 109 {0.000, 2425., 3423., 4211., 4909., 556 110 9226., 9751., 10239., 10694., 11120., 115 111 13636., 13962., 14287., 14610., 14934., 152 112 17393., 17810., 18263., 18764., 19330., 199 113 24471., 25085., 25704., 26365., 27125., 281 114 {0.000, 2661., 3692., 4487., 5182., 582 115 9504., 10030., 10509., 10946., 11346., 117 116 13570., 13846., 14118., 14388., 14658., 149 117 16694., 17041., 17417., 17834., 18311., 188 118 24158., 24791., 25408., 26049., 26769., 277 119 {0.000, 1941., 2896., 3780., 4713., 573 120 9528., 9880., 10204., 10506., 10792., 110 121 12534., 12766., 12997., 13227., 13458., 136 122 15182., 15461., 15753., 16064., 16399., 167 123 21618., 23016., 24127., 25137., 26190., 275 124 {0.000, 2028., 2996., 3875., 4801., 589 125 10310., 10665., 10992., 11298., 11588., 118 126 13436., 13695., 13957., 14223., 14494., 147 127 16660., 17031., 17425., 17848., 18305., 188 128 23121., 23943., 24754., 25586., 26507., 276 129 {0.000, 2099., 3016., 3762., 4437., 508 130 9779., 10388., 10886., 11312., 11690., 120 131 13784., 14060., 14337., 14619., 14908., 152 132 17378., 17846., 18361., 18932., 19572., 202 133 25101., 25740., 26380., 27056., 27825., 288 134 {0.000, 2216., 3042., 3663., 4192., 467 135 7768., 8547., 9421., 10221., 10870., 114 136 13626., 13935., 14240., 14543., 14846., 151 137 17199., 17610., 18056., 18549., 19104., 197 138 24563., 25214., 25858., 26532., 27296., 283 139 {0.000, 2249., 3077., 3698., 4228., 471 140 7731., 8436., 9232., 10033., 10745., 113 141 13916., 14268., 14613., 14955., 15296., 156 142 17919., 18379., 18881., 19436., 20053., 207 143 24868., 25501., 26155., 26860., 27672., 287 144 {0.000, 2206., 3020., 3628., 4143., 460 145 7205., 7698., 8237., 8833., 9491., 102 146 14009., 14525., 15026., 15519., 16009., 165 147 19590., 20134., 20681., 21225., 21762., 222 148 25408., 25997., 26643., 27368., 28202., 292 149 {0.000, 2489., 3349., 3983., 4517., 499 150 7600., 8054., 8523., 9008., 9508., 100 151 13192., 13731., 14275., 14826., 15388., 159 152 19832., 20472., 21076., 21644., 22180., 226 153 25657., 26203., 26790., 27441., 28196., 291 154 {0.000, 2945., 3913., 4642., 5270., 584 155 9019., 9517., 10001., 10472., 10930., 113 156 13898., 14312., 14731., 15159., 15600., 160 157 19641., 20371., 21073., 21726., 22329., 228 158 25825., 26323., 26854., 27442., 28132., 290 159 {0.000, 3308., 4355., 5183., 5930., 663 160 9941., 10350., 10738., 11111., 11474., 118 161 13996., 14380., 14773., 15177., 15594., 160 162 19186., 19878., 20604., 21317., 21982., 225 163 25728., 26270., 26850., 27493., 28240., 291 164 {0.000, 2774., 3718., 4452., 5121., 579 165 11225., 11684., 12085., 12450., 12790., 131 166 14933., 15242., 15558., 15885., 16223., 165 167 19115., 19633., 20184., 20765., 21375., 220 168 25629., 26216., 26825., 27482., 28229., 291 169 {0.000, 2572., 3380., 3971., 4470., 492 170 7971., 9158., 10803., 11664., 12234., 126 171 14633., 14919., 15207., 15499., 15798., 161 172 18603., 19225., 19919., 20643., 21351., 220 173 25516., 26099., 26713., 27381., 28144., 291 174 {0.000, 2876., 3700., 4293., 4786., 522 175 7452., 7832., 8235., 8674., 9173., 977 176 14066., 14471., 14845., 15199., 15542., 158 177 18123., 18606., 19163., 19825., 20623., 215 178 25454., 26017., 26620., 27294., 28088., 291 179 {0.000, 3534., 4408., 5069., 5638., 615 180 9020., 9518., 10033., 10565., 11110., 116 181 14453., 14827., 15182., 15525., 15856., 161 182 18115., 18465., 18836., 19237., 19682., 201 183 24941., 25597., 26249., 26940., 27735., 287 184 {0.000, 3063., 3895., 4531., 5092., 562 185 8879., 9398., 9888., 10352., 10795., 112 186 13685., 14100., 14517., 14935., 15353., 157 187 18195., 18594., 18994., 19398., 19810., 202 188 23466., 24265., 25162., 26124., 27154., 283 189 {0.000, 2839., 4027., 4951., 5736., 643 190 9767., 10241., 10703., 11152., 11595., 120 191 14507., 14909., 15308., 15711., 16110., 165 192 18965., 19393., 19823., 20266., 20715., 211 193 24344., 24981., 25682., 26467., 27391., 285 194 {0.000, 2839., 4027., 4951., 5736., 643 195 9767., 10241., 10703., 11152., 11595., 120 196 14507., 14909., 15308., 15711., 16110., 165 197 18965., 19393., 19823., 20266., 20715., 211 198 24344., 24981., 25682., 26467., 27391., 285 199 {0.000, 2839., 4027., 4951., 5736., 643 200 9767., 10241., 10703., 11152., 11595., 120 201 14507., 14909., 15308., 15711., 16110., 165 202 18965., 19393., 19823., 20266., 20715., 211 203 24344., 24981., 25682., 26467., 27391., 285 204 }; 205 206 void G4NRESP71M03::DKINMA(G4ReactionProduct* p 207 G4ReactionProduct* p 208 { 209 G4ReactionProduct CM; 210 G4double TotalEnergyCM; 211 212 if (p2 != nullptr) // If it is not a decay 213 { 214 // Calculating (total momentum, energy and 215 const G4ThreeVector TotalMomentumLAB = p1- 216 CM.SetMomentum(TotalMomentumLAB); 217 218 const G4double TotalEnergyLAB = p1->GetTot 219 CM.SetTotalEnergy(TotalEnergyLAB); 220 221 CM.SetMass(std::sqrt(TotalEnergyLAB * Tota 222 223 // Transforming primary particles from lab 224 p1->Lorentz(*p1, CM); 225 p2->Lorentz(*p2, CM); 226 227 TotalEnergyCM = p1->GetTotalEnergy() + p2- 228 229 const G4double m4 = 230 (p1->GetMass() + p2->GetMass()) 231 - (p3->GetMass() 232 + Q); // Mass of the residual nucleu 233 p4->SetMass(m4); 234 } 235 else // If it is a decay reaction... 236 { 237 const G4ThreeVector TotalMomentumLAB = p1- 238 CM.SetMomentum(TotalMomentumLAB); 239 240 const G4double TotalEnergyLAB = p1->GetTot 241 CM.SetTotalEnergy(TotalEnergyLAB); 242 243 CM.SetMass(std::sqrt(TotalEnergyLAB * Tota 244 245 // Transforming primary particles from lab 246 // in this case). 247 p1->Lorentz(*p1, CM); 248 249 const G4double m4 = 250 p1->GetMass() 251 - (p3->GetMass() 252 + Q); // Mass of the residual nucleu 253 p4->SetMass(m4); 254 255 TotalEnergyCM = p1->GetTotalEnergy(); 256 } 257 258 // Calculating momentum and total energy of 259 260 const G4ThreeVector p1unit = p1->GetMomentum 261 262 G4RotationMatrix rot(std::acos(p1unit * G4Th 263 std::acos(p1unit * G4Th 264 rot.inverse(); 265 266 const G4double theta = std::acos(costhcm3); 267 const G4double phi = twopi * G4UniformRand() 268 269 const G4double Energy3CM = 270 (std::pow(TotalEnergyCM, 2.) + std::pow(p3 271 / (2. * TotalEnergyCM); 272 p3->SetTotalEnergy(Energy3CM); 273 274 const G4double Momentum3CM = std::sqrt(std:: 275 p3->SetMomentum(rot 276 * G4ThreeVector(Momentum3CM 277 Momentum3CM 278 Momentum3CM 279 280 const G4double Energy4CM = TotalEnergyCM - E 281 p4->SetTotalEnergy(Energy4CM); 282 283 const G4double Momentum4CM = std::sqrt(std:: 284 p4->SetMomentum(-Momentum4CM * p3->GetMoment 285 286 // Transforming reaction products from cente 287 p3->Lorentz(*p3, -1. * CM); 288 p4->Lorentz(*p4, -1. * CM); 289 } 290 291 G4int G4NRESP71M03::ApplyMechanismI_NBeA2A(G4R 292 G4R 293 { 294 // N+12C --> A+9BE* 295 G4ReactionProduct p4; 296 297 theProds[0].SetDefinition(G4Alpha::Alpha()); 298 299 DKINMA(&neut, &carb, &(theProds[0]), &p4, QI 300 301 // 9BE* --> N+8BE 302 G4ReactionProduct p1(p4); 303 304 theProds[1].SetDefinition(G4Neutron::Neutron 305 306 DKINMA(&p1, nullptr, &(theProds[1]), &p4, -Q 307 308 // 8BE --> 2*A 309 p1 = p4; 310 311 theProds[2].SetDefinition(G4Alpha::Alpha()); 312 theProds[3].SetDefinition(G4Alpha::Alpha()); 313 314 DKINMA(&p1, nullptr, &(theProds[2]), &(thePr 315 2. * G4UniformRand() - 1.); 316 317 return 0; 318 } 319 320 // Apply kinematics for excited level selected 321 G4int G4NRESP71M03::ApplyMechanismII_ACN2A(G4R 322 G4R 323 { 324 // 12C(N,N')12C' 325 G4ReactionProduct p4; 326 327 theProds[0].SetDefinition(G4Neutron::Neutron 328 329 DKINMA(&neut, &carb, &(theProds[0]), &p4, QI 330 331 // 12C' --> ALPHA+8BE' 332 G4ReactionProduct p1(p4); 333 334 theProds[1].SetDefinition(G4Alpha::Alpha()); 335 336 DKINMA(&p1, nullptr, &(theProds[1]), &p4, -Q 337 338 // 8BE --> 2*ALPHA 339 p1 = p4; 340 341 theProds[2].SetDefinition(G4Alpha::Alpha()); 342 theProds[3].SetDefinition(G4Alpha::Alpha()); 343 344 DKINMA(&p1, nullptr, &(theProds[2]), &(thePr 345 2. * G4UniformRand() - 1.); 346 347 return 0; 348 } 349 350 G4int G4NRESP71M03::ApplyMechanismABE(G4Reacti 351 G4Reacti 352 { 353 G4double CosThetaCMAlpha = 0.; // Cosine of 354 355 G4double Kn = neut.GetKineticEnergy(); // N 356 if (Kn > 5.7 * MeV) { 357 // Sorting. 358 for (G4int i = 1; i < ndist; i++) { 359 if (BEN2[i] >= Kn / keV) { 360 // Ok. The neutron energy is between v 361 G4double BE1 = BEN2[i - 1]; 362 G4double BE2 = BEN2[i]; 363 364 // Performing energy and angle interpo 365 366 G4double FRA = 367 G4UniformRand() 368 * 49.99999999; // Sorting the bin o 369 // values of Rho fro 370 G4double DJTETA = 371 FRA 372 - G4int(FRA); // Distance in bin un 373 G4int JTETA = 374 G4int(FRA) + 1; // Getting the bin 375 376 // Calculating the angle from the cumu 377 378 G4double B11 = B2[i - 1][JTETA - 1]; 379 G4double B12 = B2[i - 1][JTETA]; 380 381 G4double TCM1 = B11 + (B12 - B11) * DJ 382 383 // Calculating the angle from the cumu 384 385 G4double B21 = B2[i][JTETA - 1]; 386 G4double B22 = B2[i][JTETA]; 387 388 G4double TCM2 = B21 + (B22 - B21) * DJ 389 390 // Interpolating the angle between ene 391 G4double TCM = (TCM1 + (TCM2 - TCM1) * 392 CosThetaCMAlpha = std::cos(TCM); 393 394 break; 395 } 396 } 397 } 398 else { 399 // Isotropic distribution in CM. 400 CosThetaCMAlpha = 1. - 2. * G4UniformRand( 401 } 402 403 // N+12C --> A+9BE 404 theProds[0].SetDefinition(G4Alpha::Alpha()); 405 theProds[1].SetDefinition(G4IonTable::GetIon 406 407 DKINMA(&neut, &carb, &(theProds[0]), &(thePr 408 409 return 0; 410 } 411