Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/im_r_matrix/src/G4AngularDistribution.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

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