source: trunk/source/processes/hadronic/models/rpg/src/G4RPGTwoBody.cc @ 910

Last change on this file since 910 was 819, checked in by garnier, 16 years ago

import all except CVS
File size: 10.4 KB
Line 
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// $Id: G4RPGTwoBody.cc,v 1.2 2007/08/15 20:38:37 dennis Exp $
27// GEANT4 tag $Name: geant4-09-01-patch-02 $
28//
29
30#include "G4RPGTwoBody.hh"
31#include "Randomize.hh"
32#include "G4Poisson.hh"
33#include <iostream>
34#include "G4HadReentrentException.hh"
35#include <signal.h>
36
37
38G4RPGTwoBody::G4RPGTwoBody()
39  : G4RPGReaction() {}
40
41
42G4bool G4RPGTwoBody::
43ReactionStage(const G4HadProjectile* /*originalIncident*/,
44              G4ReactionProduct& modifiedOriginal,
45              G4bool& /*incidentHasChanged*/,
46              const G4DynamicParticle* originalTarget,
47              G4ReactionProduct& targetParticle,
48              G4bool& /*targetHasChanged*/,
49              const G4Nucleus& targetNucleus,
50              G4ReactionProduct& currentParticle,
51              G4FastVector<G4ReactionProduct,256>& vec,
52              G4int& vecLen,
53              G4bool /*leadFlag*/,
54              G4ReactionProduct& /*leadingStrangeParticle*/)
55{
56  //
57  // Derived from H. Fesefeldt's original FORTRAN code TWOB
58  //
59  // Generation of momenta for elastic and quasi-elastic 2 body reactions
60  //
61  // The simple formula ds/d|t| = s0* std::exp(-b*|t|) is used.
62  // The b values are parametrizations from experimental data.
63  // Unavailable values are taken from those of similar reactions.
64  //
65   
66  const G4double ekOriginal = modifiedOriginal.GetKineticEnergy()/GeV;
67  const G4double etOriginal = modifiedOriginal.GetTotalEnergy()/GeV;
68  const G4double mOriginal = modifiedOriginal.GetMass()/GeV;
69  const G4double pOriginal = modifiedOriginal.GetMomentum().mag()/GeV;
70  G4double currentMass = currentParticle.GetMass()/GeV;
71  G4double targetMass = targetParticle.GetDefinition()->GetPDGMass()/GeV;
72
73  targetMass = targetParticle.GetMass()/GeV;
74  const G4double atomicWeight = targetNucleus.GetN();
75   
76  G4double etCurrent = currentParticle.GetTotalEnergy()/GeV;
77  G4double pCurrent = currentParticle.GetTotalMomentum()/GeV;
78
79  G4double cmEnergy = std::sqrt( currentMass*currentMass +
80                            targetMass*targetMass +
81                            2.0*targetMass*etCurrent );  // in GeV
82
83  if( (pCurrent < 0.1) || (cmEnergy < 0.01) ) // 2-body scattering not possible
84  {
85    targetParticle.SetMass( 0.0 );  // flag that the target particle doesn't exist
86  }
87  else
88  {
89    // Projectile momentum in cm
90
91    G4double pf = targetMass*pCurrent/cmEnergy;
92
93    //
94    // Set beam and target in centre of mass system
95    //
96    G4ReactionProduct pseudoParticle[3];
97     
98    if (targetParticle.GetDefinition()->GetParticleSubType() == "kaon" ||
99        targetParticle.GetDefinition()->GetParticleSubType() == "pi") {
100
101      //      G4double pf1 = pOriginal*mOriginal/std::sqrt(2.*mOriginal*(mOriginal+etOriginal));
102
103      pseudoParticle[0].SetMass( targetMass*GeV );
104      pseudoParticle[0].SetTotalEnergy( etOriginal*GeV );
105      pseudoParticle[0].SetMomentum( 0.0, 0.0, pOriginal*GeV );
106     
107      pseudoParticle[1].SetMomentum( 0.0, 0.0, 0.0 );
108      pseudoParticle[1].SetMass( mOriginal*GeV );
109      pseudoParticle[1].SetKineticEnergy( 0.0 );
110
111    } else {
112      pseudoParticle[0].SetMass( currentMass*GeV );
113      pseudoParticle[0].SetTotalEnergy( etCurrent*GeV );
114      pseudoParticle[0].SetMomentum( 0.0, 0.0, pCurrent*GeV );
115     
116      pseudoParticle[1].SetMomentum( 0.0, 0.0, 0.0 );
117      pseudoParticle[1].SetMass( targetMass*GeV );
118      pseudoParticle[1].SetKineticEnergy( 0.0 );
119    }
120    //
121    // Transform into center of mass system
122    //
123    pseudoParticle[2] = pseudoParticle[0] + pseudoParticle[1];
124    pseudoParticle[0].Lorentz( pseudoParticle[0], pseudoParticle[2] );
125    pseudoParticle[1].Lorentz( pseudoParticle[1], pseudoParticle[2] );
126    //
127    // Set final state masses and energies in centre of mass system
128    //
129    currentParticle.SetTotalEnergy( std::sqrt(pf*pf+currentMass*currentMass)*GeV );
130    targetParticle.SetTotalEnergy( std::sqrt(pf*pf+targetMass*targetMass)*GeV );
131
132    //
133    // Calculate slope b for elastic scattering on proton/neutron
134    //
135    const G4double cb = 0.01;
136    const G4double b1 = 4.225;
137    const G4double b2 = 1.795;
138    G4double b = std::max( cb, b1+b2*std::log(pOriginal) );
139     
140    //
141    // Get cm scattering angle by sampling t from tmin to tmax
142    //
143    G4double btrang = b * 4.0 * pf * pseudoParticle[0].GetMomentum().mag()/GeV;
144    G4double exindt = std::exp(-btrang) - 1.0;
145    G4double costheta = 1.0 + 2*std::log( 1.0+G4UniformRand()*exindt ) /btrang;
146    costheta = std::max(-1., std::min(1., costheta) );
147    G4double sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
148    G4double phi = twopi * G4UniformRand();
149
150    //
151    // Calculate final state momenta in centre of mass system
152    //
153    if (targetParticle.GetDefinition()->GetParticleSubType() == "kaon" ||
154        targetParticle.GetDefinition()->GetParticleSubType() == "pi") {
155
156      currentParticle.SetMomentum( -pf*sintheta*std::cos(phi)*GeV,
157                                   -pf*sintheta*std::sin(phi)*GeV,
158                                   -pf*costheta*GeV );
159    } else {
160
161      currentParticle.SetMomentum( pf*sintheta*std::cos(phi)*GeV,
162                                   pf*sintheta*std::sin(phi)*GeV,
163                                   pf*costheta*GeV );
164    }
165    targetParticle.SetMomentum( -currentParticle.GetMomentum() );
166
167    //
168    // Transform into lab system
169    //
170    currentParticle.Lorentz( currentParticle, pseudoParticle[1] );
171    targetParticle.Lorentz( targetParticle, pseudoParticle[1] );
172
173    // Rotate final state particle vectors wrt incident momentum
174     
175    Defs1( modifiedOriginal, currentParticle, targetParticle, vec, vecLen );
176
177    // Subtract binding energy
178
179    G4double pp, pp1, ekin;
180    if( atomicWeight >= 1.5 )
181    {
182      const G4double cfa = 0.025*((atomicWeight-1.)/120.)*std::exp(-(atomicWeight-1.)/120.);
183      pp1 = currentParticle.GetMomentum().mag()/MeV;
184      if( pp1 >= 1.0 )
185      {
186        ekin = currentParticle.GetKineticEnergy()/MeV - cfa*(1.0+0.5*normal())*GeV;
187        ekin = std::max( 0.0001*GeV, ekin );
188        currentParticle.SetKineticEnergy( ekin*MeV );
189        pp = currentParticle.GetTotalMomentum()/MeV;
190        currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
191      }
192      pp1 = targetParticle.GetMomentum().mag()/MeV;
193      if( pp1 >= 1.0 )
194      {
195        ekin = targetParticle.GetKineticEnergy()/MeV - cfa*(1.0+normal()/2.)*GeV;
196        ekin = std::max( 0.0001*GeV, ekin );
197        targetParticle.SetKineticEnergy( ekin*MeV );
198        pp = targetParticle.GetTotalMomentum()/MeV;
199        targetParticle.SetMomentum( targetParticle.GetMomentum() * (pp/pp1) );
200      }
201    }
202  }
203
204  // Get number of final state nucleons and nucleons remaining in
205  // target nucleus
206
207  std::pair<G4int, G4int> finalStateNucleons = 
208    GetFinalStateNucleons(originalTarget, vec, vecLen);
209  G4int protonsInFinalState = finalStateNucleons.first;
210  G4int neutronsInFinalState = finalStateNucleons.second;
211
212  G4int PinNucleus = std::max(0, 
213    G4int(targetNucleus.GetZ()) - protonsInFinalState);
214  G4int NinNucleus = std::max(0,
215    G4int(targetNucleus.GetN()-targetNucleus.GetZ()) - neutronsInFinalState);
216
217  if( atomicWeight >= 1.5 )
218  {
219    // Add black track particles
220    // npnb: number of proton/neutron black track particles
221    // ndta: number of deuterons, tritons, and alphas produced
222    // epnb: kinetic energy available for proton/neutron black track
223    //   particles
224    // edta: kinetic energy available for deuteron/triton/alpha particles
225
226    G4double epnb, edta;
227    G4int npnb=0, ndta=0;
228     
229    epnb = targetNucleus.GetPNBlackTrackEnergy();   // was enp1 in fortran code
230    edta = targetNucleus.GetDTABlackTrackEnergy();  // was enp3 in fortran code
231    const G4double pnCutOff = 0.0001;       // GeV
232    const G4double dtaCutOff = 0.0001;      // GeV
233    const G4double kineticMinimum = 0.0001;
234    const G4double kineticFactor = -0.010;
235    G4double sprob = 0.0; // sprob = probability of self-absorption in heavy molecules
236    if( epnb >= pnCutOff )
237    {
238      npnb = G4Poisson( epnb/0.02 );
239      if( npnb > atomicWeight )npnb = G4int(atomicWeight);
240      if( (epnb > pnCutOff) && (npnb <= 0) )npnb = 1;
241      npnb = std::min( npnb, 127-vecLen );
242    }
243    if( edta >= dtaCutOff )
244    {
245      ndta = G4int(2.0 * std::log(atomicWeight));
246      ndta = std::min( ndta, 127-vecLen );
247    }
248
249    if (npnb == 0 && ndta == 0) npnb = 1;
250
251    AddBlackTrackParticles(epnb, npnb, edta, ndta, sprob, kineticMinimum, 
252                           kineticFactor, modifiedOriginal, 
253                           PinNucleus, NinNucleus, targetNucleus,
254                           vec, vecLen);
255  }
256
257  //
258  //  calculate time delay for nuclear reactions
259  //
260  if( (atomicWeight >= 1.5) && (atomicWeight <= 230.0) && (ekOriginal <= 0.2) )
261    currentParticle.SetTOF( 1.0-500.0*std::exp(-ekOriginal/0.04)*std::log(G4UniformRand()) );
262  else
263    currentParticle.SetTOF( 1.0 );
264
265  return true;
266}
267 
268/* end of file */
Note: See TracBrowser for help on using the repository browser.