ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4NRESP71M03.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4NRESP71M03.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 #include "G4NRESP71M03.hh"
27 
28 #include "G4Geantino.hh"
29 #include "G4Neutron.hh"
30 #include "G4Alpha.hh"
31 #include "G4IonTable.hh"
32 #include "G4RotationMatrix.hh"
33 
34 #include "G4PhysicalConstants.hh"
35 #include "G4SystemOfUnits.hh"
36 #include "Randomize.hh"
37 
38 // Energies for which angular distributions are given.
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. };
40 // Angular distribution of alpha particles for each energy value in BEN2.
41 const G4double G4NRESP71M03::B2[ndist][nrhos] = {
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 */
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 */
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 */
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 */
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 */
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 */
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 */
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 */
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 */
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 */
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 */
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 */
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 */
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 */
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 */
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 */
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 */
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 */
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 */
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 */
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 */
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 */
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 */
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 */
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 */
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 */
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 */
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 */
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 */
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 */
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 */
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*/
74  };
75 
76 void G4NRESP71M03::DKINMA(G4ReactionProduct *p1, G4ReactionProduct *p2, G4ReactionProduct *p3, G4ReactionProduct *p4, const G4double Q, const G4double costhcm3)
77  {
78  G4ReactionProduct CM;
79  G4double TotalEnergyCM;
80 
81  if ( p2 ) // If it is not a decay reaction...
82  {
83  // Calculating (total momentum, energy and mass) of the center of mass.
84  const G4ThreeVector TotalMomentumLAB = p1->GetMomentum()+p2->GetMomentum();
85  CM.SetMomentum(TotalMomentumLAB);
86 
87  const G4double TotalEnergyLAB = p1->GetTotalEnergy()+p2->GetTotalEnergy();
88  CM.SetTotalEnergy(TotalEnergyLAB);
89 
90  CM.SetMass(std::sqrt(TotalEnergyLAB*TotalEnergyLAB-TotalMomentumLAB*TotalMomentumLAB));
91 
92  // Transforming primary particles from laboratory to center of mass system.
93  p1->Lorentz(*p1, CM);
94  p2->Lorentz(*p2, CM);
95 
96  TotalEnergyCM = p1->GetTotalEnergy()+p2->GetTotalEnergy();
97 
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).
99  p4->SetMass(m4);
100  }
101  else // If it is a decay reaction...
102  {
103  const G4ThreeVector TotalMomentumLAB = p1->GetMomentum();
104  CM.SetMomentum(TotalMomentumLAB);
105 
106  const G4double TotalEnergyLAB = p1->GetTotalEnergy();
107  CM.SetTotalEnergy(TotalEnergyLAB);
108 
109  CM.SetMass(std::sqrt(TotalEnergyLAB*TotalEnergyLAB-TotalMomentumLAB*TotalMomentumLAB));
110 
111  // Transforming primary particles from laboratory to center of mass system (not really necessary in this case).
112  p1->Lorentz(*p1, CM);
113 
114  const G4double m4 = p1->GetMass()-(p3->GetMass()+Q); // Mass of the residual nucleus in the excited state (not in the ground state).
115  p4->SetMass(m4);
116 
117  TotalEnergyCM = p1->GetTotalEnergy();
118  }
119 
120  // Calculating momentum and total energy of the reaction products.
121 
122  const G4ThreeVector p1unit = p1->GetMomentum().unit();
123 
124  G4RotationMatrix rot(std::acos(p1unit*G4ThreeVector(0., 1., 0.)), std::acos(p1unit*G4ThreeVector(0., 0., 1.)), 0.);
125  rot.inverse();
126 
127  const G4double theta = std::acos(costhcm3);
128  const G4double phi = twopi*G4UniformRand();
129 
130  const G4double Energy3CM = (std::pow(TotalEnergyCM, 2.)+std::pow(p3->GetMass(), 2.)-std::pow(p4->GetMass(), 2.))/(2.*TotalEnergyCM);
131  p3->SetTotalEnergy(Energy3CM);
132 
133  const G4double Momentum3CM = std::sqrt(std::pow(Energy3CM, 2.)-std::pow(p3->GetMass(), 2.));
134  p3->SetMomentum(rot*G4ThreeVector(Momentum3CM*std::sin(theta)*std::cos(phi), Momentum3CM*std::sin(theta)*std::sin(phi), Momentum3CM*costhcm3));
135 
136  const G4double Energy4CM = TotalEnergyCM-Energy3CM;
137  p4->SetTotalEnergy(Energy4CM);
138 
139  const G4double Momentum4CM = std::sqrt(std::pow(Energy4CM, 2.)-std::pow(p4->GetMass(), 2.));
140  p4->SetMomentum(-Momentum4CM*p3->GetMomentum().unit());
141 
142  // Transforming reaction products from center of mass to laboratory system.
143  p3->Lorentz(*p3, -1.*CM);
144  p4->Lorentz(*p4, -1.*CM);
145  }
146 
147 G4int G4NRESP71M03::ApplyMechanismI_NBeA2A(G4ReactionProduct &neut, G4ReactionProduct &carb, G4ReactionProduct *theProds, const G4double QI)
148  {
149  // N+12C --> A+9BE*
150  G4ReactionProduct p4;
151 
152  theProds[0].SetDefinition(G4Alpha::Alpha());
153 
154  DKINMA(&neut, &carb, &(theProds[0]), &p4, QI, 2.*G4UniformRand()-1.);
155 
156  // 9BE* --> N+8BE
157  G4ReactionProduct p1(p4);
158 
159  theProds[1].SetDefinition(G4Neutron::Neutron());
160 
161  DKINMA(&p1, NULL, &(theProds[1]), &p4, -QI-7.369, 2.*G4UniformRand()-1.);
162 
163  // 8BE --> 2*A
164  p1 = p4;
165 
166  theProds[2].SetDefinition(G4Alpha::Alpha());
167  theProds[3].SetDefinition(G4Alpha::Alpha());
168 
169  DKINMA(&p1, NULL, &(theProds[2]), &(theProds[3]), 0.09538798439007223351, 2.*G4UniformRand()-1.);
170 
171  return 0;
172  }
173 
174 // Apply kinematics for excited level selected by GEANT4 (it).
175 G4int G4NRESP71M03::ApplyMechanismII_ACN2A(G4ReactionProduct &neut, G4ReactionProduct &carb, G4ReactionProduct *theProds, const G4double QI)
176  {
177  // 12C(N,N')12C'
178  G4ReactionProduct p4;
179 
180  theProds[0].SetDefinition(G4Neutron::Neutron());
181 
182  DKINMA(&neut, &carb, &(theProds[0]), &p4, QI, 2.*G4UniformRand()-1.);
183 
184  // 12C' --> ALPHA+8BE'
185  G4ReactionProduct p1(p4);
186 
187  theProds[1].SetDefinition(G4Alpha::Alpha());
188 
189  DKINMA(&p1, NULL, &(theProds[1]), &p4, -QI-7.369, 2.*G4UniformRand()-1.);
190 
191  // 8BE --> 2*ALPHA
192  p1 = p4;
193 
194  theProds[2].SetDefinition(G4Alpha::Alpha());
195  theProds[3].SetDefinition(G4Alpha::Alpha());
196 
197  DKINMA(&p1, NULL, &(theProds[2]), &(theProds[3]), 0.09538798439007223351, 2.*G4UniformRand()-1.);
198 
199  return 0;
200  }
201 
202 G4int G4NRESP71M03::ApplyMechanismABE(G4ReactionProduct &neut, G4ReactionProduct &carb, G4ReactionProduct *theProds)
203  {
204  G4double CosThetaCMAlpha = 0.; // Cosine of the angle of emission of the alpha particle (theta).
205 
206  G4double Kn = neut.GetKineticEnergy(); // Neutron energy in the center of mass system.
207  if ( Kn > 5.7*MeV )
208  {
209  // Sorting.
210  for ( G4int i=1; i<ndist; i++ )
211  {
212  if ( BEN2[i] >= Kn/keV )
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];
217 
218  // Performing energy and angle interpolation.
219 
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).
221  G4double DJTETA = FRA-G4int(FRA); // Distance in bin units (DJTETA) from the low edge of the bin with Rho.
222  G4int JTETA = G4int(FRA)+1; // Getting the bin (JTETA) next to the bin with the value of Rho.
223 
224  // Calculating the angle from the cumulative distribution at energy i.
225 
226  G4double B11 = B2[i-1][JTETA-1];
227  G4double B12 = B2[i-1][JTETA];
228 
229  G4double TCM1 = B11+(B12-B11)*DJTETA; // Angle at energy i.
230 
231  // Calculating the angle from the cumulative distribution at energy i-1.
232 
233  G4double B21 = B2[i][JTETA-1];
234  G4double B22 = B2[i][JTETA];
235 
236  G4double TCM2 = B21+(B22-B21)*DJTETA; // Angle at energy i-1.
237 
238  // Interpolating the angle between energy values i and i-1.
239  G4double TCM = (TCM1+(TCM2-TCM1)*(Kn/keV-BE1)/(BE2-BE1))*1.E-4;
240  CosThetaCMAlpha = std::cos(TCM);
241 
242  break;
243  }
244  }
245  }
246  else
247  {
248  // Isotropic distribution in CM.
249  CosThetaCMAlpha = 1.-2.*G4UniformRand();
250  }
251 
252  //N+12C --> A+9BE
253  theProds[0].SetDefinition(G4Alpha::Alpha());
254  theProds[1].SetDefinition(G4IonTable::GetIonTable()->GetIon(4, 9, 0. ));
255 
256  DKINMA(&neut, &carb, &(theProds[0]), &(theProds[1]), -5.71*MeV, CosThetaCMAlpha);
257 
258  return 0;
259  }