source: trunk/source/processes/hadronic/models/im_r_matrix/src/G4AngularDistribution.cc @ 836

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

import all except CVS
File size: 12.6 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// hpw: done, but low quality at present.
27
28#include "globals.hh"
29#include "G4AngularDistribution.hh"
30#include "Randomize.hh"
31
32G4AngularDistribution::G4AngularDistribution(G4bool symmetrize)
33  : sym(symmetrize)
34{
35  // The following are parameters of the model - not to be confused with the PDG values!
36
37      mSigma = 0.55; 
38      cmSigma = 1.20; 
39      gSigma  = 9.4;
40
41      mOmega = 0.783; 
42      cmOmega = 0.808; 
43      gOmega = 10.95;
44
45      mPion = 0.138; 
46      cmPion = 0.51; 
47      gPion  = 7.27;
48
49      mNucleon = 0.938;   
50
51      // Definition of constants for pion-Term (no s-dependence)
52
53      m42   = 4. * mNucleon * mNucleon;
54      mPion2  = mPion * mPion;           
55      cmPion2 = cmPion * cmPion;
56      dPion1 = cmPion2-mPion2;
57      dPion2 = dPion1 * dPion1;
58      cm6gp = 1.5 * (cmPion2*cmPion2*cmPion2) * (gPion*gPion*gPion*gPion) * m42 * m42 / dPion2;
59
60      cPion_3 = -(cm6gp/3.);
61      cPion_2 = -(cm6gp * mPion2/dPion1);
62      cPion_1 = -(cm6gp * mPion2 * (2. * cmPion2 + mPion2) / dPion2);
63      cPion_m = -(cm6gp * cmPion2 * mPion2 / dPion2);
64      cPion_L = -(cm6gp * 2. * cmPion2 * mPion2 * (cmPion2 + mPion2) / dPion2 / dPion1);
65      cPion_0 = -(cPion_3 + cPion_2 + cPion_1 + cPion_m);
66   
67     // Definition of constants for sigma-Term (no s-dependence)
68
69      G4double gSigmaSq = gSigma * gSigma; 
70
71      mSigma2  = mSigma * mSigma;
72      cmSigma2 = cmSigma * cmSigma;
73      cmSigma4 = cmSigma2 * cmSigma2;
74      cmSigma6 = cmSigma2 * cmSigma4;
75      dSigma1 = m42 - cmSigma2;
76      dSigma2 = m42 - mSigma2;
77      dSigma3 = cmSigma2 - mSigma2;
78
79      G4double dSigma1Sq = dSigma1 * dSigma1;
80      G4double dSigma2Sq = dSigma2 * dSigma2;
81      G4double dSigma3Sq = dSigma3 * dSigma3;
82
83      cm2gs = 0.5 * cmSigma2 * gSigmaSq*gSigmaSq / dSigma3Sq;     
84
85
86      cSigma_3 = -(cm2gs * dSigma1Sq / 3.);
87      cSigma_2 = -(cm2gs * cmSigma2 * dSigma1 * dSigma2 / dSigma3);
88      cSigma_1 = -(cm2gs * cmSigma4 * (2. * dSigma1 + dSigma2) * dSigma2 / dSigma3Sq);
89      cSigma_m = -(cm2gs * cmSigma6 * dSigma2Sq / mSigma2 / dSigma3Sq);
90      cSigma_L = -(cm2gs * cmSigma6 * dSigma2 * (dSigma1 + dSigma2) * 2. / (dSigma3 * dSigma3Sq));
91      cSigma_0 = -(cSigma_3 + cSigma_2 + cSigma_1 + cSigma_m);
92
93      // Definition of constants for omega-Term
94
95      G4double gOmegaSq = gOmega * gOmega;
96
97      mOmega2  = mOmega * mOmega;
98      cmOmega2 = cmOmega * cmOmega;
99      cmOmega4 = cmOmega2 * cmOmega2;
100      cmOmega6 = cmOmega2 * cmOmega4;
101      dOmega1 = m42 - cmOmega2;
102      dOmega2 = m42 - mOmega2;
103      dOmega3 = cmOmega2 - mOmega2;
104      sOmega1 = cmOmega2 + mOmega2;
105
106      G4double dOmega3Sq = dOmega3 * dOmega3;
107
108      cm2go = 0.5 * cmOmega2 * gOmegaSq * gOmegaSq / dOmega3Sq;
109
110      cOmega_3 =  cm2go / 3.;
111      cOmega_2 = -(cm2go * cmOmega2 / dOmega3);
112      cOmega_1 =  cm2go * cmOmega4 / dOmega3Sq;
113      cOmega_m =  cm2go * cmOmega6 / (dOmega3Sq * mOmega2);
114      cOmega_L = -(cm2go * cmOmega6 * 4. / (dOmega3 * dOmega3Sq));
115
116      // Definition of constants for mix-Term
117
118      G4double fac1Tmp = (gSigma * gOmega * cmSigma2 * cmOmega2);
119      fac1 = -(fac1Tmp *  fac1Tmp * m42); 
120      dMix1 = cmOmega2 - cmSigma2;
121      dMix2 = cmOmega2 - mSigma2;
122      dMix3 = cmSigma2 - mOmega2;
123
124      G4double dMix1Sq = dMix1 * dMix1;
125      G4double dMix2Sq = dMix2 * dMix2;
126      G4double dMix3Sq = dMix3 * dMix3;
127
128      cMix_o1 =    fac1 / (cmOmega2  * dMix1Sq * dMix2 * dOmega3);
129      cMix_s1 =    fac1 / (cmSigma2  * dMix1Sq * dMix3 * dSigma3);
130      cMix_Omega = fac1 / (dOmega3Sq * dMix3Sq * (mOmega2 - mSigma2));
131      cMix_sm =    fac1 / (dSigma3Sq * dMix2Sq * (mSigma2 - mOmega2)); 
132      fac2 = (-fac1) / (dMix1*dMix1Sq * dOmega3Sq * dMix2Sq);
133      fac3 = (-fac1) / (dMix1*dMix1Sq * dSigma3Sq * dMix3Sq); 
134     
135      cMix_oLc = fac2 * (3. * cmOmega2*cmOmega4            - cmOmega4 * cmSigma2
136                       - 2. * cmOmega4 * mOmega2           - 2. * cmOmega4 * mSigma2
137                       + cmOmega2 * mOmega2 * mSigma2      + cmSigma2 * mOmega2 * mSigma2
138                       - 4. * cmOmega4 * m42               + 2. * cmOmega2 * cmSigma2 * m42
139                       + 3. * cmOmega2 * mOmega2 * m42     - cmSigma2 * mOmega2 * m42
140                       + 3. * cmOmega2 * mSigma2 * m42     - cmSigma2 * mSigma2 * m42
141                       - 2. * mOmega2 * mSigma2 * m42);
142     
143      cMix_oLs = fac2 * (8. * cmOmega4                     - 4. * cmOmega2 * cmSigma2
144                         - 6. * cmOmega2 * mOmega2         + 2. * cmSigma2 * mOmega2
145                         - 6. * cmOmega2 * mSigma2         + 2. * cmSigma2 * mSigma2
146                         + 4. * mOmega2 * mSigma2);
147     
148      cMix_sLc = fac3 * (cmOmega2 * cmSigma4               - 3. * cmSigma6   
149                         + 2. * cmSigma4 * mOmega2         + 2. * cmSigma4 * mSigma2
150                         - cmOmega2 * mOmega2 * mSigma2    - cmSigma2 * mOmega2 * mSigma2
151                         - 2. * cmOmega2 * cmSigma2 * m42  + 4. * cmSigma4 * m42
152                         + cmOmega2 * mOmega2 * m42        - 3. * cmSigma2 * mOmega2 * m42
153                         + cmOmega2 * mSigma2 * m42        - 3. * cmSigma2 * mSigma2 * m42
154                         + 2. * mOmega2 * mSigma2 * m42);
155
156      cMix_sLs = fac3 * (4. * cmOmega2 * cmSigma2          - 8. * cmSigma4
157                       - 2. * cmOmega2 * mOmega2           + 6. * cmSigma2 * mOmega2
158                       - 2. * cmOmega2 * mSigma2           + 6. * cmSigma2 * mSigma2
159                       - 4. * mOmega2 * mSigma2);
160}
161
162
163G4AngularDistribution::~
164G4AngularDistribution()
165{ }
166
167
168G4double G4AngularDistribution::CosTheta(G4double s, G4double m1, G4double m2) const
169{
170   G4double random = G4UniformRand();
171   G4double dCosTheta = 2.;
172   G4double cosTheta = -1.;
173
174   // For jmax=12 the accuracy is better than 0.1 degree
175   G4int jMax = 12;
176
177   for (G4int j = 1; j <= jMax; ++j)
178   {
179      // Accuracy is 2^-jmax
180      dCosTheta *= 0.5;
181      G4double cosTh = cosTheta + dCosTheta;
182      if(DifferentialCrossSection(s, m1, m2, cosTh) <= random) cosTheta = cosTh;
183    }
184
185   // Randomize in final interval in order to avoid discrete angles
186   cosTheta += G4UniformRand() * dCosTheta;
187
188
189   if (cosTheta > 1. || cosTheta < -1.)
190     throw G4HadronicException(__FILE__, __LINE__, "G4AngularDistribution::CosTheta - std::cos(theta) outside allowed range");
191
192   return cosTheta;
193}
194
195
196G4double G4AngularDistribution::DifferentialCrossSection(G4double sIn, G4double m1, G4double m2, 
197                                                         G4double cosTheta) const
198{
199// local calculus is in GeV, ie. normalize input
200  sIn = sIn/sqr(GeV)+m42/2.;
201  m1  = m1/GeV;
202  m2  = m2/GeV;
203//  G4cout << "Here we go"<<sIn << " "<<m1 << " " << m2 <<" " m42<< G4endl;
204// scaling from masses other than p,p.
205  G4double s = sIn - (m1+m2) * (m1+m2) + m42;
206  G4double tMax = s - m42;
207  G4double tp = 0.5 * (cosTheta + 1.) * tMax; 
208  G4double twoS = 2. * s;
209 
210  // Define s-dependent stuff for omega-Term
211  G4double brak1 = (twoS-m42) * (twoS-m42);
212  G4double bOmega_3 = cOmega_3 * (-2. * cmOmega4 - 2. * cmOmega2 * twoS - brak1);
213  G4double bOmega_2 = cOmega_2 * ( 2. * cmOmega2 * mOmega2 + sOmega1 * twoS + brak1);
214  G4double bOmega_1 = cOmega_1 * (-4. * cmOmega2 * mOmega2
215                         - 2. * mOmega2*mOmega2
216                         - 2. * (cmOmega2 + 2 * mOmega2) * twoS
217                         - 3. * brak1);
218  G4double bOmega_m = cOmega_m * (-2. * mOmega2*mOmega2 - 2. * mOmega2 * twoS - brak1);
219  G4double bOmega_L = cOmega_L * (sOmega1 * mOmega2 + (cmOmega2 + 3. * mOmega2) * s + brak1);
220  G4double bOmega_0 = -(bOmega_3 + bOmega_2 + bOmega_1 + bOmega_m);
221 
222  // Define s-dependent stuff for mix-Term           
223  G4double bMix_o1 = cMix_o1 * (dOmega1 - twoS);
224  G4double bMix_s1 = cMix_s1 * (dSigma1 - twoS);
225  G4double bMix_Omega = cMix_Omega * (dOmega2 - twoS);
226  G4double bMix_sm = cMix_sm * (dSigma2 - twoS);
227  G4double bMix_oL = cMix_oLc + cMix_oLs * s;
228  G4double bMix_sL = cMix_sLc + cMix_sLs * s;
229 
230  G4double t1_Pion = 1. / (1. + tMax / cmPion2);
231  G4double t2_Pion = 1. + tMax / mPion2;
232  G4double t1_Sigma = 1. / (1. + tMax / cmSigma2);
233  G4double t2_Sigma = 1. + tMax / mSigma2;
234  G4double t1_Omega = 1. / (1. + tMax / cmOmega2);
235  G4double t2_Omega = 1. + tMax / mOmega2;
236 
237  G4double norm = Cross(t1_Pion, t1_Sigma, t1_Omega,
238                        t2_Pion, t2_Sigma, t2_Omega,
239                        bMix_o1, bMix_s1, bMix_Omega,
240                        bMix_sm, bMix_oL, bMix_sL,
241                        bOmega_0, bOmega_1, bOmega_2,
242                        bOmega_3, bOmega_m, bOmega_L);
243 
244  t1_Pion = 1. / (1. + tp / cmPion2);
245  t2_Pion = 1. + tp / mPion2;
246  t1_Sigma = 1. / (1. + tp / cmSigma2);
247  t2_Sigma = 1. + tp / mSigma2;
248  t1_Omega = 1. / (1. + tp / cmOmega2);
249  t2_Omega = 1. + tp / mOmega2;
250 
251  G4double dSigma;
252  if (sym) 
253    { 
254      G4double to;
255      norm = 2. * norm;
256      to = tMax - tp;
257      G4double t3_Pion = 1. / (1. + to / cmPion2);
258      G4double t4_Pion = 1. + to / mPion2;
259      G4double t3_Sigma = 1. / (1. + to / cmSigma2);
260      G4double t4_Sigma = 1. + to / mSigma2;
261      G4double t3_Omega = 1. / (1. + to / cmOmega2);
262      G4double t4_Omega = 1. + to / mOmega2;
263     
264      dSigma = ( Cross(t1_Pion, t1_Sigma, t1_Omega,
265                       t2_Pion,t2_Sigma, t2_Omega,
266                       bMix_o1, bMix_s1, bMix_Omega,
267                       bMix_sm, bMix_oL, bMix_sL,
268                       bOmega_0, bOmega_1, bOmega_2,
269                       bOmega_3, bOmega_m, bOmega_L) - 
270                 Cross(t3_Pion,t3_Sigma, t3_Omega,
271                       t4_Pion, t4_Sigma, t4_Omega,
272                       bMix_o1, bMix_s1, bMix_Omega,
273                       bMix_sm, bMix_oL, bMix_sL,
274                       bOmega_0, bOmega_1, bOmega_2,
275                       bOmega_3, bOmega_m, bOmega_L) ) 
276        / norm + 0.5;
277    }
278  else
279    {
280      dSigma = Cross(t1_Pion, t1_Sigma, t1_Omega, 
281                     t2_Pion, t2_Sigma, t2_Omega,
282                     bMix_o1, bMix_s1, bMix_Omega,
283                     bMix_sm, bMix_oL, bMix_sL,
284                     bOmega_0, bOmega_1, bOmega_2,
285                     bOmega_3, bOmega_m, bOmega_L) 
286        / norm;
287    }
288 
289  return dSigma;
290}
291
292
293G4double G4AngularDistribution::Cross(G4double tpPion,
294                                      G4double tpSigma,
295                                      G4double tpOmega,
296                                      G4double tmPion,
297                                      G4double tmSigma,
298                                      G4double tmOmega,
299                                      G4double bMix_o1,
300                                      G4double bMix_s1,
301                                      G4double bMix_Omega,
302                                      G4double bMix_sm,
303                                      G4double bMix_oL,
304                                      G4double bMix_sL,
305                                      G4double bOmega_0,
306                                      G4double bOmega_1,
307                                      G4double bOmega_2,
308                                      G4double bOmega_3,
309                                      G4double bOmega_m,
310                                      G4double bOmega_L) const
311{
312  G4double cross = 0;
313     //  Pion
314    cross += ((cPion_3 * tpPion  + cPion_2)  * tpPion  + cPion_1)  * tpPion  + cPion_m/tmPion   + cPion_0  +  cPion_L * std::log(tpPion*tmPion);
315//    G4cout << "cross1 "<< cross<<G4endl;
316    //  Sigma
317    cross += ((cSigma_3 * tpSigma + cSigma_2) * tpSigma + cSigma_1) * tpSigma + cSigma_m/tmSigma + cSigma_0 + cSigma_L * std::log(tpSigma*tmSigma);
318//    G4cout << "cross2 "<< cross<<G4endl;
319    // Omega
320    cross += ((bOmega_3 * tpOmega + bOmega_2) * tpOmega + bOmega_1) * tpOmega + bOmega_m/tmOmega + bOmega_0 + bOmega_L * std::log(tpOmega*tmOmega)
321    // Mix
322    +  bMix_o1 * (tpOmega - 1.)
323    +  bMix_s1 * (tpSigma - 1.)
324    +  bMix_Omega * std::log(tmOmega)
325    +  bMix_sm * std::log(tmSigma)
326    +  bMix_oL * std::log(tpOmega)
327    +  bMix_sL * std::log(tpSigma);
328/*      G4cout << "cross3 "<< cross<<" "
329             <<bMix_o1<<" "
330             <<bMix_s1<<" "
331             <<bMix_Omega<<" "
332             <<bMix_sm<<" "
333             <<bMix_oL<<" "
334             <<bMix_sL<<" "
335             <<tpOmega<<" "
336             <<tpSigma<<" "
337             <<tmOmega<<" "
338             <<tmSigma<<" "
339             <<tpOmega<<" "
340             <<tpSigma
341             <<G4endl;
342*/
343  return cross;
344
345}
Note: See TracBrowser for help on using the repository browser.