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 #include "G4NRESP71M03.hh" 26 #include "G4NRESP71M03.hh" 27 27 28 #include "G4Alpha.hh" << 29 #include "G4Geantino.hh" 28 #include "G4Geantino.hh" 30 #include "G4IonTable.hh" << 31 #include "G4Neutron.hh" 29 #include "G4Neutron.hh" 32 #include "G4PhysicalConstants.hh" << 30 #include "G4Alpha.hh" >> 31 #include "G4IonTable.hh" 33 #include "G4RotationMatrix.hh" 32 #include "G4RotationMatrix.hh" >> 33 >> 34 #include "G4PhysicalConstants.hh" 34 #include "G4SystemOfUnits.hh" 35 #include "G4SystemOfUnits.hh" 35 #include "Randomize.hh" 36 #include "Randomize.hh" 36 37 37 // Energies for which angular distributions ar 38 // Energies for which angular distributions are given. 38 const G4double G4NRESP71M03::BEN2[ndist] = { << 39 const G4double G4NRESP71M03::BEN2[ndist] = { 5700., 8000., 8640., 8990., 9220., 9410., 9830., 10400., 10800., 11250., 11460., 11870., 12140., 12320., 12570., 12940., 13420., 13760., 14020., 14200., 14440., 14620., 14820., 15050., 15660., 15980., 16470., 16940., 17970., 18000., 19000., 20000. }; 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 40 // Angular distribution of alpha particles for each energy value in BEN2. 43 const G4double G4NRESP71M03::B2[ndist][nrhos] 41 const G4double G4NRESP71M03::B2[ndist][nrhos] = { 44 {0.000, 2839., 4027., 4951., 5736., 643 << 42 { 0.000, 2839., 4027., 4951., 5736., 6437., 7074., 7671., 8230., 8765., 9273., 9767., 10241., 10703., 11152., 11595., 12026., 12453., 12871., 13285., 13697., 14102., 14507., 14909., 15308., 15711., 16110., 16509., 16911., 17316., 17721., 18133., 18547., 18965., 19393., 19823., 20266., 20715., 21177., 21651., 22145., 22654., 23188., 23747., 24344., 24981., 25682., 26467., 27391., 28579., 31416. }, /* En=5700 keV */ 45 9767., 10241., 10703., 11152., 11595., 120 << 43 { 0.000, 1771., 2528., 3129., 3653., 4128., 4574., 4998., 5403., 5802., 6192., 6578., 6961., 7345., 7725., 8108., 8494., 8887., 9280., 9676., 10078., 10483., 10891., 11306., 11721., 12142., 12563., 12987., 13411., 13841., 14272., 14705., 15142., 15585., 16034., 16490., 16958., 17438., 17935., 18453., 18994., 19572., 20187., 20857., 21595., 22424., 23373., 24482., 25820., 27539., 31416. }, /* En=8000 keV */ 46 14507., 14909., 15308., 15711., 16110., 165 << 44 { 0.000, 3518., 4863., 5805., 6540., 7143., 7665., 8121., 8535., 8912., 9261., 9588., 9896., 10191., 10470., 10744., 11004., 11259., 11507., 11749., 11988., 12223., 12456., 12685., 12915., 13144., 13373., 13603., 13832., 14064., 14300., 14539., 14784., 15032., 15286., 15547., 15821., 16100., 16399., 16710., 17043., 17401., 17790., 18221., 18711., 19283., 19980., 20897., 22251., 24554., 31416. }, /* En=8640 keV */ 47 18965., 19393., 19823., 20266., 20715., 211 << 45 { 0.000, 2287., 3364., 4319., 5277., 6327., 7536., 8868., 10043., 10964., 11696., 12305., 12839., 13317., 13757., 14168., 14558., 14931., 15293., 15645., 15993., 16339., 16681., 17024., 17372., 17724., 18086., 18453., 18833., 19226., 19638., 20065., 20511., 20976., 21463., 21965., 22481., 23005., 23527., 24048., 24560., 25069., 25569., 26068., 26568., 27080., 27614., 28183., 28820., 29612., 31416. }, /* En=8990 keV */ 48 24344., 24981., 25682., 26467., 27391., 285 << 46 { 0.000, 1690., 2456., 3097., 3697., 4294., 4919., 5614., 6433., 7517., 9076., 10722., 11906., 12795., 13515., 14134., 14686., 15189., 15654., 16094., 16515., 16920., 17313., 17699., 18082., 18463., 18849., 19239., 19638., 20049., 20480., 20932., 21406., 21912., 22446., 23002., 23571., 24133., 24677., 25195., 25682., 26147., 26590., 27023., 27451., 27878., 28321., 28789., 29314., 29955., 31416. }, /* En=9220 keV */ 49 {0.000, 1771., 2528., 3129., 3653., 412 << 47 { 0.000, 1674., 2434., 3078., 3685., 4297., 4951., 5695., 6619., 7913., 9597., 11014., 12069., 12918., 13640., 14278., 14853., 15378., 15865., 16323., 16757., 17171., 17574., 17966., 18353., 18733., 19116., 19502., 19892., 20294., 20712., 21146., 21601., 22085., 22591., 23122., 23665., 24212., 24743., 25248., 25735., 26194., 26637., 27067., 27492., 27919., 28355., 28820., 29336., 29973., 31416. }, /* En=9410 keV */ 50 6578., 6961., 7345., 7725., 8108., 849 << 48 { 0.000, 2302., 3245., 3967., 4577., 5117., 5617., 6085., 6534., 6971., 7401., 7838., 8278., 8739., 9229., 9767., 10382., 11140., 12186., 13637., 14818., 15604., 16198., 16688., 17121., 17514., 17881., 18230., 18569., 18902., 19232., 19568., 19911., 20269., 20646., 21051., 21497., 22009., 22619., 23364., 24253., 25148., 25921., 26571., 27146., 27677., 28189., 28710., 29267., 29936., 31416. }, /* En=9830 keV */ 51 10891., 11306., 11721., 12142., 12563., 129 << 49 { 0.000, 2195., 2934., 3458., 3879., 4244., 4567., 4866., 5142., 5406., 5658., 5899., 6138., 6371., 6600., 6832., 7062., 7291., 7527., 7766., 8011., 8265., 8529., 8809., 9104., 9424., 9773., 10159., 10602., 11124., 11762., 12572., 13643., 15060., 16801., 18312., 19352., 20140., 20794., 21378., 21925., 22449., 22971., 23505., 24067., 24674., 25355., 26153., 27137., 28443., 31416. }, /* En=10400 keV */ 52 15585., 16034., 16490., 16958., 17438., 179 << 50 { 0.000, 1712., 2406., 2937., 3380., 3773., 4134., 4470., 4787., 5095., 5394., 5689., 5978., 6270., 6565., 6864., 7175., 7495., 7831., 8190., 8579., 9003., 9484., 10034., 10675., 11422., 12258., 13147., 14055., 14994., 15996., 17077., 18164., 19138., 19971., 20696., 21353., 21965., 22553., 23128., 23703., 24281., 24865., 25462., 26068., 26684., 27316., 27972., 28689., 29543., 31416. }, /* En=10800 keV */ 53 21595., 22424., 23373., 24482., 25820., 275 << 51 { 0.000, 1853., 2651., 3289., 3851., 4366., 4853., 5324., 5786., 6245., 6707., 7172., 7646., 8139., 8645., 9179., 9735., 10326., 10945., 11592., 12261., 12943., 13621., 14287., 14928., 15544., 16128., 16691., 17228., 17746., 18249., 18736., 19213., 19682., 20143., 20602., 21061., 21523., 21984., 22452., 22933., 23423., 23929., 24457., 25013., 25603., 26247., 26964., 27803., 28874., 31416. }, /* En=11250 keV */ 54 {0.000, 3518., 4863., 5805., 6540., 714 << 52 { 0.000, 1959., 2770., 3397., 3934., 4417., 4866., 5294., 5709., 6118., 6529., 6948., 7382., 7840., 8333., 8876., 9484., 10167., 10902., 11629., 12292., 12877., 13398., 13871., 14309., 14725., 15125., 15517., 15906., 16298., 16697., 17109., 17540., 17992., 18473., 18982., 19519., 20077., 20646., 21217., 21783., 22344., 22906., 23477., 24069., 24699., 25391., 26185., 27146., 28422., 31416. }, /* En=11460 keV */ 55 9588., 9896., 10191., 10470., 10744., 110 << 53 { 0.000, 2290., 3391., 4447., 5824., 8351., 9448., 10077., 10553., 10952., 11306., 11630., 11935., 12227., 12510., 12789., 13066., 13344., 13628., 13918., 14220., 14538., 14875., 15239., 15637., 16077., 16565., 17101., 17671., 18251., 18820., 19368., 19897., 20413., 20927., 21449., 21989., 22558., 23162., 23802., 24459., 25110., 25731., 26317., 26873., 27410., 27944., 28493., 29091., 29810., 31416. }, /* En=11870 keV */ 56 12456., 12685., 12915., 13144., 13373., 136 << 54 { 0.000, 2119., 3072., 3874., 4627., 5378., 6155., 6966., 7786., 8568., 9281., 9921., 10501., 11034., 11531., 12002., 12452., 12886., 13306., 13715., 14114., 14504., 14888., 15266., 15641., 16013., 16386., 16761., 17142., 17534., 17940., 18367., 18821., 19312., 19849., 20438., 21072., 21724., 22360., 22957., 23511., 24030., 24525., 25008., 25491., 25987., 26514., 27098., 27791., 28730., 31416. }, /* En=12140 keV */ 57 15032., 15286., 15547., 15821., 16100., 163 << 55 { 0.000, 2425., 3423., 4211., 4909., 5563., 6200., 6829., 7455., 8070., 8664., 9226., 9751., 10239., 10694., 11120., 11522., 11904., 12271., 12625., 12969., 13305., 13636., 13962., 14287., 14610., 14934., 15261., 15591., 15928., 16273., 16630., 17002., 17393., 17810., 18263., 18764., 19330., 19985., 20742., 21575., 22399., 23153., 23836., 24471., 25085., 25704., 26365., 27125., 28138., 31416. }, /* En=12320 keV */ 58 18711., 19283., 19980., 20897., 22251., 245 << 56 { 0.000, 2661., 3692., 4487., 5182., 5829., 6457., 7081., 7706., 8329., 8933., 9504., 10030., 10509., 10946., 11346., 11717., 12063., 12390., 12701., 13000., 13289., 13570., 13846., 14118., 14388., 14658., 14929., 15203., 15482., 15767., 16062., 16370., 16694., 17041., 17417., 17834., 18311., 18879., 19596., 20545., 21666., 22662., 23468., 24158., 24791., 25408., 26049., 26769., 27705., 31416. }, /* En=12570 keV */ 59 {0.000, 2287., 3364., 4319., 5277., 632 << 57 { 0.000, 1941., 2896., 3780., 4713., 5731., 6720., 7532., 8172., 8693., 9137., 9528., 9880., 10204., 10506., 10792., 11064., 11325., 11578., 11824., 12064., 12301., 12534., 12766., 12997., 13227., 13458., 13691., 13926., 14165., 14409., 14659., 14916., 15182., 15461., 15753., 16064., 16399., 16765., 17175., 17647., 18217., 18962., 20057., 21618., 23016., 24127., 25137., 26190., 27504., 31416. }, /* En=12940 keV */ 60 12305., 12839., 13317., 13757., 14168., 145 << 58 { 0.000, 2028., 2996., 3875., 4801., 5899., 7158., 8182., 8905., 9458., 9914., 10310., 10665., 10992., 11298., 11588., 11868., 12138., 12403., 12664., 12922., 13179., 13436., 13695., 13957., 14223., 14494., 14772., 15057., 15352., 15658., 15977., 16310., 16660., 17031., 17425., 17848., 18305., 18805., 19358., 19977., 20675., 21451., 22281., 23121., 23943., 24754., 25586., 26507., 27691., 31416. }, /* En=13420 keV */ 61 16681., 17024., 17372., 17724., 18086., 184 << 59 { 0.000, 2099., 3016., 3762., 4437., 5088., 5747., 6450., 7233., 8116., 9014., 9779., 10388., 10886., 11312., 11690., 12034., 12354., 12657., 12948., 13231., 13509., 13784., 14060., 14337., 14619., 14908., 15206., 15516., 15842., 16187., 16554., 16950., 17378., 17846., 18361., 18932., 19572., 20296., 21115., 22008., 22896., 23707., 24434., 25101., 25740., 26380., 27056., 27825., 28822., 31416. }, /* En=13760 keV */ 62 20976., 21463., 21965., 22481., 23005., 235 << 60 { 0.000, 2216., 3042., 3663., 4192., 4675., 5137., 5595., 6065., 6565., 7121., 7768., 8547., 9421., 10221., 10870., 11403., 11859., 12264., 12634., 12980., 13309., 13626., 13935., 14240., 14543., 14846., 15152., 15463., 15783., 16113., 16456., 16817., 17199., 17610., 18056., 18549., 19104., 19743., 20491., 21350., 22257., 23110., 23873., 24563., 25214., 25858., 26532., 27296., 28301., 31416. }, /* En=14020 keV */ 63 26568., 27080., 27614., 28183., 28820., 296 << 61 { 0.000, 2249., 3077., 3698., 4228., 4711., 5172., 5628., 6094., 6586., 7123., 7731., 8436., 9232., 10033., 10745., 11353., 11879., 12348., 12776., 13175., 13553., 13916., 14268., 14613., 14955., 15296., 15639., 15986., 16342., 16709., 17090., 17492., 17919., 18379., 18881., 19436., 20053., 20732., 21456., 22190., 22904., 23585., 24235., 24868., 25501., 26155., 26860., 27672., 28727., 31416. }, /* En=14200 keV */ 64 {0.000, 1690., 2456., 3097., 3697., 429 << 62 { 0.000, 2206., 3020., 3628., 4143., 4607., 5043., 5466., 5884., 6307., 6745., 7205., 7698., 8237., 8833., 9491., 10201., 10929., 11634., 12293., 12903., 13472., 14009., 14525., 15026., 15519., 16009., 16500., 16995., 17496., 18005., 18524., 19053., 19590., 20134., 20681., 21225., 21762., 22289., 22805., 23314., 23821., 24333., 24858., 25408., 25997., 26643., 27368., 28202., 29223., 31416. }, /* En=14440 keV */ 65 10722., 11906., 12795., 13515., 14134., 146 << 63 { 0.000, 2489., 3349., 3983., 4517., 4998., 5447., 5879., 6304., 6729., 7160., 7600., 8054., 8523., 9008., 9508., 10019., 10539., 11065., 11593., 12124., 12657., 13192., 13731., 14275., 14826., 15388., 15964., 16559., 17177., 17821., 18488., 19164., 19832., 20472., 21076., 21644., 22180., 22692., 23187., 23673., 24156., 24643., 25141., 25657., 26203., 26790., 27441., 28196., 29158., 31416. }, /* En=14620 keV */ 66 17313., 17699., 18082., 18463., 18849., 192 << 64 { 0.000, 2945., 3913., 4642., 5270., 5848., 6400., 6937., 7468., 7992., 8510., 9019., 9517., 10001., 10472., 10930., 11375., 11811., 12238., 12658., 13073., 13486., 13898., 14312., 14731., 15159., 15600., 16059., 16544., 17063., 17627., 18245., 18921., 19641., 20371., 21073., 21726., 22329., 22888., 23413., 23913., 24397., 24871., 25344., 25825., 26323., 26854., 27442., 28132., 29045., 31416. }, /* En=14820 keV */ 67 21912., 22446., 23002., 23571., 24133., 246 << 65 { 0.000, 3308., 4355., 5183., 5930., 6638., 7311., 7938., 8510., 9030., 9504., 9941., 10350., 10738., 11111., 11474., 11831., 12184., 12538., 12894., 13255., 13622., 13996., 14380., 14773., 15177., 15594., 16024., 16472., 16941., 17437., 17970., 18549., 19186., 19878., 20604., 21317., 21982., 22593., 23158., 23691., 24203., 24706., 25211., 25728., 26270., 26850., 27493., 28240., 29191., 31416. }, /* En=15050 keV */ 68 27451., 27878., 28321., 28789., 29314., 299 << 66 { 0.000, 2774., 3718., 4452., 5121., 5796., 6548., 7503., 8803., 9927., 10670., 11225., 11684., 12085., 12450., 12790., 13113., 13425., 13729., 14030., 14329., 14629., 14933., 15242., 15558., 15885., 16223., 16575., 16943., 17330., 17738., 18170., 18628., 19115., 19633., 20184., 20765., 21375., 22002., 22634., 23260., 23871., 24467., 25050., 25629., 26216., 26825., 27482., 28229., 29174., 31416. }, /* En=15660 keV */ 69 {0.000, 1674., 2434., 3078., 3685., 429 << 67 { 0.000, 2572., 3380., 3971., 4470., 4924., 5357., 5787., 6230., 6709., 7260., 7971., 9158., 10803., 11664., 12234., 12684., 13071., 13419., 13742., 14048., 14344., 14633., 14919., 15207., 15499., 15798., 16110., 16437., 16787., 17167., 17587., 18060., 18603., 19225., 19919., 20643., 21351., 22020., 22648., 23246., 23823., 24387., 24948., 25516., 26099., 26713., 27381., 28144., 29113., 31416. }, /* En=15980 keV */ 70 11014., 12069., 12918., 13640., 14278., 148 << 68 { 0.000, 2876., 3700., 4293., 4786., 5221., 5622., 6001., 6366., 6726., 7086., 7452., 7832., 8235., 8674., 9173., 9772., 10542., 11506., 12399., 13077., 13612., 14066., 14471., 14845., 15199., 15542., 15880., 16219., 16564., 16919., 17292., 17689., 18123., 18606., 19163., 19825., 20623., 21523., 22383., 23124., 23768., 24352., 24906., 25454., 26017., 26620., 27294., 28088., 29106., 31416. }, /* En=16470 keV */ 71 17574., 17966., 18353., 18733., 19116., 195 << 69 { 0.000, 3534., 4408., 5069., 5638., 6157., 6649., 7126., 7595., 8063., 8537., 9020., 9518., 10033., 10565., 11110., 11658., 12195., 12707., 13189., 13639., 14059., 14453., 14827., 15182., 15525., 15856., 16181., 16500., 16817., 17134., 17453., 17779., 18115., 18465., 18836., 19237., 19682., 20194., 20809., 21581., 22515., 23441., 24240., 24941., 25597., 26249., 26940., 27735., 28775., 31416. }, /* En=16940 keV */ 72 22085., 22591., 23122., 23665., 24212., 247 << 70 { 0.000, 3063., 3895., 4531., 5092., 5621., 6144., 6674., 7219., 7776., 8334., 8879., 9398., 9888., 10352., 10795., 11221., 11638., 12048., 12455., 12863., 13272., 13685., 14100., 14517., 14935., 15353., 15768., 16181., 16590., 16995., 17397., 17797., 18195., 18594., 18994., 19398., 19810., 20232., 20671., 21133., 21628., 22168., 22773., 23466., 24265., 25162., 26124., 27154., 28382., 31416. }, /* En=17970 keV */ 73 27492., 27919., 28355., 28820., 29336., 299 << 71 { 0.000, 2839., 4027., 4951., 5736., 6437., 7074., 7671., 8230., 8765., 9273., 9767., 10241., 10703., 11152., 11595., 12026., 12453., 12871., 13285., 13697., 14102., 14507., 14909., 15308., 15711., 16110., 16509., 16911., 17316., 17721., 18133., 18547., 18965., 19393., 19823., 20266., 20715., 21177., 21651., 22145., 22654., 23188., 23747., 24344., 24981., 25682., 26467., 27391., 28579., 31416. }, /* En=18000 keV */ 74 {0.000, 2302., 3245., 3967., 4577., 511 << 72 { 0.000, 2839., 4027., 4951., 5736., 6437., 7074., 7671., 8230., 8765., 9273., 9767., 10241., 10703., 11152., 11595., 12026., 12453., 12871., 13285., 13697., 14102., 14507., 14909., 15308., 15711., 16110., 16509., 16911., 17316., 17721., 18133., 18547., 18965., 19393., 19823., 20266., 20715., 21177., 21651., 22145., 22654., 23188., 23747., 24344., 24981., 25682., 26467., 27391., 28579., 31416. }, /* En=19000 keV */ 75 7838., 8278., 8739., 9229., 9767., 103 << 73 { 0.000, 2839., 4027., 4951., 5736., 6437., 7074., 7671., 8230., 8765., 9273., 9767., 10241., 10703., 11152., 11595., 12026., 12453., 12871., 13285., 13697., 14102., 14507., 14909., 15308., 15711., 16110., 16509., 16911., 17316., 17721., 18133., 18547., 18965., 19393., 19823., 20266., 20715., 21177., 21651., 22145., 22654., 23188., 23747., 24344., 24981., 25682., 26467., 27391., 28579., 31416. } /* En=20000 keV*/ 76 16198., 16688., 17121., 17514., 17881., 182 << 74 }; 77 20269., 20646., 21051., 21497., 22009., 226 << 75 78 27146., 27677., 28189., 28710., 29267., 299 << 76 void G4NRESP71M03::DKINMA(G4ReactionProduct *p1, G4ReactionProduct *p2, G4ReactionProduct *p3, G4ReactionProduct *p4, const G4double Q, const G4double costhcm3) 79 {0.000, 2195., 2934., 3458., 3879., 424 << 77 { 80 5899., 6138., 6371., 6600., 6832., 706 << 78 G4ReactionProduct CM; 81 8529., 8809., 9104., 9424., 9773., 101 << 79 G4double TotalEnergyCM; 82 15060., 16801., 18312., 19352., 20140., 207 << 80 83 24067., 24674., 25355., 26153., 27137., 284 << 81 if ( p2 ) // If it is not a decay reaction... 84 {0.000, 1712., 2406., 2937., 3380., 377 << 82 { 85 5689., 5978., 6270., 6565., 6864., 717 << 83 // Calculating (total momentum, energy and mass) of the center of mass. 86 9484., 10034., 10675., 11422., 12258., 131 << 84 const G4ThreeVector TotalMomentumLAB = p1->GetMomentum()+p2->GetMomentum(); 87 19138., 19971., 20696., 21353., 21965., 225 << 85 CM.SetMomentum(TotalMomentumLAB); 88 26068., 26684., 27316., 27972., 28689., 295 << 86 89 {0.000, 1853., 2651., 3289., 3851., 436 << 87 const G4double TotalEnergyLAB = p1->GetTotalEnergy()+p2->GetTotalEnergy(); 90 7172., 7646., 8139., 8645., 9179., 973 << 88 CM.SetTotalEnergy(TotalEnergyLAB); 91 13621., 14287., 14928., 15544., 16128., 166 << 89 92 19682., 20143., 20602., 21061., 21523., 219 << 90 CM.SetMass(std::sqrt(TotalEnergyLAB*TotalEnergyLAB-TotalMomentumLAB*TotalMomentumLAB)); 93 25013., 25603., 26247., 26964., 27803., 288 << 91 94 {0.000, 1959., 2770., 3397., 3934., 441 << 92 // Transforming primary particles from laboratory to center of mass system. 95 6948., 7382., 7840., 8333., 8876., 948 << 93 p1->Lorentz(*p1, CM); 96 13398., 13871., 14309., 14725., 15125., 155 << 94 p2->Lorentz(*p2, CM); 97 17992., 18473., 18982., 19519., 20077., 206 << 95 98 24069., 24699., 25391., 26185., 27146., 284 << 96 TotalEnergyCM = p1->GetTotalEnergy()+p2->GetTotalEnergy(); 99 {0.000, 2290., 3391., 4447., 5824., 835 << 97 100 11630., 11935., 12227., 12510., 12789., 130 << 98 const G4double m4 = (p1->GetMass()+p2->GetMass())-(p3->GetMass()+Q); // Mass of the residual nucleus in the excited state (not in the ground state). 101 14875., 15239., 15637., 16077., 16565., 171 << 99 p4->SetMass(m4); 102 20413., 20927., 21449., 21989., 22558., 231 << 100 } 103 26873., 27410., 27944., 28493., 29091., 298 << 101 else // If it is a decay reaction... 104 {0.000, 2119., 3072., 3874., 4627., 537 << 102 { 105 9921., 10501., 11034., 11531., 12002., 124 << 103 const G4ThreeVector TotalMomentumLAB = p1->GetMomentum(); 106 14888., 15266., 15641., 16013., 16386., 167 << 104 CM.SetMomentum(TotalMomentumLAB); 107 19312., 19849., 20438., 21072., 21724., 223 << 105 108 25491., 25987., 26514., 27098., 27791., 287 << 106 const G4double TotalEnergyLAB = p1->GetTotalEnergy(); 109 {0.000, 2425., 3423., 4211., 4909., 556 << 107 CM.SetTotalEnergy(TotalEnergyLAB); 110 9226., 9751., 10239., 10694., 11120., 115 << 108 111 13636., 13962., 14287., 14610., 14934., 152 << 109 CM.SetMass(std::sqrt(TotalEnergyLAB*TotalEnergyLAB-TotalMomentumLAB*TotalMomentumLAB)); 112 17393., 17810., 18263., 18764., 19330., 199 << 110 113 24471., 25085., 25704., 26365., 27125., 281 << 111 // Transforming primary particles from laboratory to center of mass system (not really necessary in this case). 114 {0.000, 2661., 3692., 4487., 5182., 582 << 112 p1->Lorentz(*p1, CM); 115 9504., 10030., 10509., 10946., 11346., 117 << 113 116 13570., 13846., 14118., 14388., 14658., 149 << 114 const G4double m4 = p1->GetMass()-(p3->GetMass()+Q); // Mass of the residual nucleus in the excited state (not in the ground state). 117 16694., 17041., 17417., 17834., 18311., 188 << 115 p4->SetMass(m4); 118 24158., 24791., 25408., 26049., 26769., 277 << 116 119 {0.000, 1941., 2896., 3780., 4713., 573 << 117 TotalEnergyCM = p1->GetTotalEnergy(); 120 9528., 9880., 10204., 10506., 10792., 110 << 118 } 121 12534., 12766., 12997., 13227., 13458., 136 << 119 122 15182., 15461., 15753., 16064., 16399., 167 << 120 // Calculating momentum and total energy of the reaction products. 123 21618., 23016., 24127., 25137., 26190., 275 << 121 124 {0.000, 2028., 2996., 3875., 4801., 589 << 122 const G4ThreeVector p1unit = p1->GetMomentum().unit(); 125 10310., 10665., 10992., 11298., 11588., 118 << 123 126 13436., 13695., 13957., 14223., 14494., 147 << 124 G4RotationMatrix rot(std::acos(p1unit*G4ThreeVector(0., 1., 0.)), std::acos(p1unit*G4ThreeVector(0., 0., 1.)), 0.); 127 16660., 17031., 17425., 17848., 18305., 188 << 125 rot.inverse(); 128 23121., 23943., 24754., 25586., 26507., 276 << 126 129 {0.000, 2099., 3016., 3762., 4437., 508 << 127 const G4double theta = std::acos(costhcm3); 130 9779., 10388., 10886., 11312., 11690., 120 << 128 const G4double phi = twopi*G4UniformRand(); 131 13784., 14060., 14337., 14619., 14908., 152 << 129 132 17378., 17846., 18361., 18932., 19572., 202 << 130 const G4double Energy3CM = (std::pow(TotalEnergyCM, 2.)+std::pow(p3->GetMass(), 2.)-std::pow(p4->GetMass(), 2.))/(2.*TotalEnergyCM); 133 25101., 25740., 26380., 27056., 27825., 288 << 131 p3->SetTotalEnergy(Energy3CM); 134 {0.000, 2216., 3042., 3663., 4192., 467 << 132 135 7768., 8547., 9421., 10221., 10870., 114 << 133 const G4double Momentum3CM = std::sqrt(std::pow(Energy3CM, 2.)-std::pow(p3->GetMass(), 2.)); 136 13626., 13935., 14240., 14543., 14846., 151 << 134 p3->SetMomentum(rot*G4ThreeVector(Momentum3CM*std::sin(theta)*std::cos(phi), Momentum3CM*std::sin(theta)*std::sin(phi), Momentum3CM*costhcm3)); 137 17199., 17610., 18056., 18549., 19104., 197 << 135 138 24563., 25214., 25858., 26532., 27296., 283 << 136 const G4double Energy4CM = TotalEnergyCM-Energy3CM; 139 {0.000, 2249., 3077., 3698., 4228., 471 << 137 p4->SetTotalEnergy(Energy4CM); 140 7731., 8436., 9232., 10033., 10745., 113 << 138 141 13916., 14268., 14613., 14955., 15296., 156 << 139 const G4double Momentum4CM = std::sqrt(std::pow(Energy4CM, 2.)-std::pow(p4->GetMass(), 2.)); 142 17919., 18379., 18881., 19436., 20053., 207 << 140 p4->SetMomentum(-Momentum4CM*p3->GetMomentum().unit()); 143 24868., 25501., 26155., 26860., 27672., 287 << 141 144 {0.000, 2206., 3020., 3628., 4143., 460 << 142 // Transforming reaction products from center of mass to laboratory system. 145 7205., 7698., 8237., 8833., 9491., 102 << 143 p3->Lorentz(*p3, -1.*CM); 146 14009., 14525., 15026., 15519., 16009., 165 << 144 p4->Lorentz(*p4, -1.*CM); 147 19590., 20134., 20681., 21225., 21762., 222 << 145 } 148 25408., 25997., 26643., 27368., 28202., 292 << 146 149 {0.000, 2489., 3349., 3983., 4517., 499 << 147 G4int G4NRESP71M03::ApplyMechanismI_NBeA2A(G4ReactionProduct &neut, G4ReactionProduct &carb, G4ReactionProduct *theProds, const G4double QI) 150 7600., 8054., 8523., 9008., 9508., 100 << 148 { 151 13192., 13731., 14275., 14826., 15388., 159 << 149 // N+12C --> A+9BE* 152 19832., 20472., 21076., 21644., 22180., 226 << 150 G4ReactionProduct p4; 153 25657., 26203., 26790., 27441., 28196., 291 << 151 154 {0.000, 2945., 3913., 4642., 5270., 584 << 152 theProds[0].SetDefinition(G4Alpha::Alpha()); 155 9019., 9517., 10001., 10472., 10930., 113 << 153 156 13898., 14312., 14731., 15159., 15600., 160 << 154 DKINMA(&neut, &carb, &(theProds[0]), &p4, QI, 2.*G4UniformRand()-1.); 157 19641., 20371., 21073., 21726., 22329., 228 << 155 158 25825., 26323., 26854., 27442., 28132., 290 << 156 // 9BE* --> N+8BE 159 {0.000, 3308., 4355., 5183., 5930., 663 << 157 G4ReactionProduct p1(p4); 160 9941., 10350., 10738., 11111., 11474., 118 << 158 161 13996., 14380., 14773., 15177., 15594., 160 << 159 theProds[1].SetDefinition(G4Neutron::Neutron()); 162 19186., 19878., 20604., 21317., 21982., 225 << 160 163 25728., 26270., 26850., 27493., 28240., 291 << 161 DKINMA(&p1, NULL, &(theProds[1]), &p4, -QI-7.369, 2.*G4UniformRand()-1.); 164 {0.000, 2774., 3718., 4452., 5121., 579 << 162 165 11225., 11684., 12085., 12450., 12790., 131 << 163 // 8BE --> 2*A 166 14933., 15242., 15558., 15885., 16223., 165 << 164 p1 = p4; 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 165 311 theProds[2].SetDefinition(G4Alpha::Alpha()); << 166 theProds[2].SetDefinition(G4Alpha::Alpha()); 312 theProds[3].SetDefinition(G4Alpha::Alpha()); << 167 theProds[3].SetDefinition(G4Alpha::Alpha()); 313 168 314 DKINMA(&p1, nullptr, &(theProds[2]), &(thePr << 169 DKINMA(&p1, NULL, &(theProds[2]), &(theProds[3]), 0.09538798439007223351, 2.*G4UniformRand()-1.); 315 2. * G4UniformRand() - 1.); << 316 170 317 return 0; << 171 return 0; 318 } << 172 } 319 173 320 // Apply kinematics for excited level selected 174 // Apply kinematics for excited level selected by GEANT4 (it). 321 G4int G4NRESP71M03::ApplyMechanismII_ACN2A(G4R << 175 G4int G4NRESP71M03::ApplyMechanismII_ACN2A(G4ReactionProduct &neut, G4ReactionProduct &carb, G4ReactionProduct *theProds, const G4double QI) 322 G4R << 176 { 323 { << 177 // 12C(N,N')12C' 324 // 12C(N,N')12C' << 178 G4ReactionProduct p4; 325 G4ReactionProduct p4; << 326 179 327 theProds[0].SetDefinition(G4Neutron::Neutron << 180 theProds[0].SetDefinition(G4Neutron::Neutron()); 328 181 329 DKINMA(&neut, &carb, &(theProds[0]), &p4, QI << 182 DKINMA(&neut, &carb, &(theProds[0]), &p4, QI, 2.*G4UniformRand()-1.); 330 183 331 // 12C' --> ALPHA+8BE' << 184 // 12C' --> ALPHA+8BE' 332 G4ReactionProduct p1(p4); << 185 G4ReactionProduct p1(p4); 333 186 334 theProds[1].SetDefinition(G4Alpha::Alpha()); << 187 theProds[1].SetDefinition(G4Alpha::Alpha()); 335 188 336 DKINMA(&p1, nullptr, &(theProds[1]), &p4, -Q << 189 DKINMA(&p1, NULL, &(theProds[1]), &p4, -QI-7.369, 2.*G4UniformRand()-1.); 337 190 338 // 8BE --> 2*ALPHA << 191 // 8BE --> 2*ALPHA 339 p1 = p4; << 192 p1 = p4; 340 193 341 theProds[2].SetDefinition(G4Alpha::Alpha()); << 194 theProds[2].SetDefinition(G4Alpha::Alpha()); 342 theProds[3].SetDefinition(G4Alpha::Alpha()); << 195 theProds[3].SetDefinition(G4Alpha::Alpha()); 343 196 344 DKINMA(&p1, nullptr, &(theProds[2]), &(thePr << 197 DKINMA(&p1, NULL, &(theProds[2]), &(theProds[3]), 0.09538798439007223351, 2.*G4UniformRand()-1.); 345 2. * G4UniformRand() - 1.); << 346 198 347 return 0; << 199 return 0; 348 } << 200 } 349 201 350 G4int G4NRESP71M03::ApplyMechanismABE(G4Reacti << 202 G4int G4NRESP71M03::ApplyMechanismABE(G4ReactionProduct &neut, G4ReactionProduct &carb, G4ReactionProduct *theProds) 351 G4Reacti << 203 { 352 { << 204 G4double CosThetaCMAlpha = 0.; // Cosine of the angle of emission of the alpha particle (theta). 353 G4double CosThetaCMAlpha = 0.; // Cosine of << 354 205 355 G4double Kn = neut.GetKineticEnergy(); // N << 206 G4double Kn = neut.GetKineticEnergy(); // Neutron energy in the center of mass system. 356 if (Kn > 5.7 * MeV) { << 207 if ( Kn > 5.7*MeV ) 357 // Sorting. << 208 { 358 for (G4int i = 1; i < ndist; i++) { << 209 // Sorting. 359 if (BEN2[i] >= Kn / keV) { << 210 for ( G4int i=1; i<ndist; i++ ) 360 // Ok. The neutron energy is between v << 211 { 361 G4double BE1 = BEN2[i - 1]; << 212 if ( BEN2[i] >= Kn/keV ) 362 G4double BE2 = BEN2[i]; << 213 { >> 214 // Ok. The neutron energy is between values i-1 and i of BEN2. Taking them. >> 215 G4double BE1 = BEN2[i-1]; >> 216 G4double BE2 = BEN2[i]; 363 217 364 // Performing energy and angle interpo << 218 // Performing energy and angle interpolation. 365 219 366 G4double FRA = << 220 G4double FRA = G4UniformRand()*49.99999999; // Sorting the bin of the cumulative probability FRA (Rho). There are 51 values of Rho from 0 to 1 (0; 0.02; 0.04 ... 1). 367 G4UniformRand() << 221 G4double DJTETA = FRA-G4int(FRA); // Distance in bin units (DJTETA) from the low edge of the bin with Rho. 368 * 49.99999999; // Sorting the bin o << 222 G4int JTETA = G4int(FRA)+1; // Getting the bin (JTETA) next to the bin with the value of Rho. 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 223 376 // Calculating the angle from the cumu << 224 // Calculating the angle from the cumulative distribution at energy i. 377 225 378 G4double B11 = B2[i - 1][JTETA - 1]; << 226 G4double B11 = B2[i-1][JTETA-1]; 379 G4double B12 = B2[i - 1][JTETA]; << 227 G4double B12 = B2[i-1][JTETA]; 380 228 381 G4double TCM1 = B11 + (B12 - B11) * DJ << 229 G4double TCM1 = B11+(B12-B11)*DJTETA; // Angle at energy i. 382 230 383 // Calculating the angle from the cumu << 231 // Calculating the angle from the cumulative distribution at energy i-1. 384 232 385 G4double B21 = B2[i][JTETA - 1]; << 233 G4double B21 = B2[i][JTETA-1]; 386 G4double B22 = B2[i][JTETA]; << 234 G4double B22 = B2[i][JTETA]; 387 235 388 G4double TCM2 = B21 + (B22 - B21) * DJ << 236 G4double TCM2 = B21+(B22-B21)*DJTETA; // Angle at energy i-1. 389 237 390 // Interpolating the angle between ene << 238 // Interpolating the angle between energy values i and i-1. 391 G4double TCM = (TCM1 + (TCM2 - TCM1) * << 239 G4double TCM = (TCM1+(TCM2-TCM1)*(Kn/keV-BE1)/(BE2-BE1))*1.E-4; 392 CosThetaCMAlpha = std::cos(TCM); << 240 CosThetaCMAlpha = std::cos(TCM); 393 241 394 break; << 242 break; 395 } << 243 } 396 } << 244 } 397 } << 245 } 398 else { << 246 else 399 // Isotropic distribution in CM. << 247 { 400 CosThetaCMAlpha = 1. - 2. * G4UniformRand( << 248 // Isotropic distribution in CM. 401 } << 249 CosThetaCMAlpha = 1.-2.*G4UniformRand(); >> 250 } 402 251 403 // N+12C --> A+9BE << 252 //N+12C --> A+9BE 404 theProds[0].SetDefinition(G4Alpha::Alpha()); << 253 theProds[0].SetDefinition(G4Alpha::Alpha()); 405 theProds[1].SetDefinition(G4IonTable::GetIon << 254 theProds[1].SetDefinition(G4IonTable::GetIonTable()->GetIon(4, 9, 0. )); 406 255 407 DKINMA(&neut, &carb, &(theProds[0]), &(thePr << 256 DKINMA(&neut, &carb, &(theProds[0]), &(theProds[1]), -5.71*MeV, CosThetaCMAlpha); 408 257 409 return 0; << 258 return 0; 410 } << 259 } 411 260