G4AngularDistribution.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // hpw: done, but low quality at present.
00027 
00028 #include "globals.hh"
00029 #include "G4SystemOfUnits.hh"
00030 #include "G4AngularDistribution.hh"
00031 #include "Randomize.hh"
00032 
00033 G4AngularDistribution::G4AngularDistribution(G4bool symmetrize)
00034   : sym(symmetrize)
00035 {
00036   // The following are parameters of the model - not to be confused with the PDG values!
00037 
00038       mSigma = 0.55; 
00039       cmSigma = 1.20; 
00040       gSigma  = 9.4;
00041 
00042       mOmega = 0.783; 
00043       cmOmega = 0.808; 
00044       gOmega = 10.95;
00045 
00046       mPion = 0.138; 
00047       cmPion = 0.51; 
00048       gPion  = 7.27;
00049 
00050       mNucleon = 0.938;   
00051 
00052       // Definition of constants for pion-Term (no s-dependence)
00053 
00054       m42   = 4. * mNucleon * mNucleon;
00055       mPion2  = mPion * mPion;           
00056       cmPion2 = cmPion * cmPion;
00057       dPion1 = cmPion2-mPion2;
00058       dPion2 = dPion1 * dPion1;
00059       cm6gp = 1.5 * (cmPion2*cmPion2*cmPion2) * (gPion*gPion*gPion*gPion) * m42 * m42 / dPion2;
00060 
00061       cPion_3 = -(cm6gp/3.);
00062       cPion_2 = -(cm6gp * mPion2/dPion1);
00063       cPion_1 = -(cm6gp * mPion2 * (2. * cmPion2 + mPion2) / dPion2);
00064       cPion_m = -(cm6gp * cmPion2 * mPion2 / dPion2);
00065       cPion_L = -(cm6gp * 2. * cmPion2 * mPion2 * (cmPion2 + mPion2) / dPion2 / dPion1);
00066       cPion_0 = -(cPion_3 + cPion_2 + cPion_1 + cPion_m);
00067     
00068      // Definition of constants for sigma-Term (no s-dependence)
00069 
00070       G4double gSigmaSq = gSigma * gSigma; 
00071 
00072       mSigma2  = mSigma * mSigma;
00073       cmSigma2 = cmSigma * cmSigma;
00074       cmSigma4 = cmSigma2 * cmSigma2;
00075       cmSigma6 = cmSigma2 * cmSigma4;
00076       dSigma1 = m42 - cmSigma2;
00077       dSigma2 = m42 - mSigma2;
00078       dSigma3 = cmSigma2 - mSigma2;
00079 
00080       G4double dSigma1Sq = dSigma1 * dSigma1;
00081       G4double dSigma2Sq = dSigma2 * dSigma2;
00082       G4double dSigma3Sq = dSigma3 * dSigma3;
00083 
00084       cm2gs = 0.5 * cmSigma2 * gSigmaSq*gSigmaSq / dSigma3Sq;     
00085 
00086 
00087       cSigma_3 = -(cm2gs * dSigma1Sq / 3.);
00088       cSigma_2 = -(cm2gs * cmSigma2 * dSigma1 * dSigma2 / dSigma3);
00089       cSigma_1 = -(cm2gs * cmSigma4 * (2. * dSigma1 + dSigma2) * dSigma2 / dSigma3Sq);
00090       cSigma_m = -(cm2gs * cmSigma6 * dSigma2Sq / mSigma2 / dSigma3Sq);
00091       cSigma_L = -(cm2gs * cmSigma6 * dSigma2 * (dSigma1 + dSigma2) * 2. / (dSigma3 * dSigma3Sq));
00092       cSigma_0 = -(cSigma_3 + cSigma_2 + cSigma_1 + cSigma_m);
00093 
00094       // Definition of constants for omega-Term
00095 
00096       G4double gOmegaSq = gOmega * gOmega;
00097 
00098       mOmega2  = mOmega * mOmega;
00099       cmOmega2 = cmOmega * cmOmega;
00100       cmOmega4 = cmOmega2 * cmOmega2;
00101       cmOmega6 = cmOmega2 * cmOmega4;
00102       dOmega1 = m42 - cmOmega2;
00103       dOmega2 = m42 - mOmega2;
00104       dOmega3 = cmOmega2 - mOmega2;
00105       sOmega1 = cmOmega2 + mOmega2;
00106 
00107       G4double dOmega3Sq = dOmega3 * dOmega3;
00108 
00109       cm2go = 0.5 * cmOmega2 * gOmegaSq * gOmegaSq / dOmega3Sq;
00110 
00111       cOmega_3 =  cm2go / 3.;
00112       cOmega_2 = -(cm2go * cmOmega2 / dOmega3);
00113       cOmega_1 =  cm2go * cmOmega4 / dOmega3Sq;
00114       cOmega_m =  cm2go * cmOmega6 / (dOmega3Sq * mOmega2);
00115       cOmega_L = -(cm2go * cmOmega6 * 4. / (dOmega3 * dOmega3Sq));
00116 
00117       // Definition of constants for mix-Term
00118 
00119       G4double fac1Tmp = (gSigma * gOmega * cmSigma2 * cmOmega2);
00120       fac1 = -(fac1Tmp *  fac1Tmp * m42);  
00121       dMix1 = cmOmega2 - cmSigma2;
00122       dMix2 = cmOmega2 - mSigma2;
00123       dMix3 = cmSigma2 - mOmega2;
00124 
00125       G4double dMix1Sq = dMix1 * dMix1;
00126       G4double dMix2Sq = dMix2 * dMix2;
00127       G4double dMix3Sq = dMix3 * dMix3;
00128 
00129       cMix_o1 =    fac1 / (cmOmega2  * dMix1Sq * dMix2 * dOmega3);
00130       cMix_s1 =    fac1 / (cmSigma2  * dMix1Sq * dMix3 * dSigma3);
00131       cMix_Omega = fac1 / (dOmega3Sq * dMix3Sq * (mOmega2 - mSigma2));
00132       cMix_sm =    fac1 / (dSigma3Sq * dMix2Sq * (mSigma2 - mOmega2)); 
00133       fac2 = (-fac1) / (dMix1*dMix1Sq * dOmega3Sq * dMix2Sq);
00134       fac3 = (-fac1) / (dMix1*dMix1Sq * dSigma3Sq * dMix3Sq); 
00135       
00136       cMix_oLc = fac2 * (3. * cmOmega2*cmOmega4            - cmOmega4 * cmSigma2 
00137                        - 2. * cmOmega4 * mOmega2           - 2. * cmOmega4 * mSigma2 
00138                        + cmOmega2 * mOmega2 * mSigma2      + cmSigma2 * mOmega2 * mSigma2 
00139                        - 4. * cmOmega4 * m42               + 2. * cmOmega2 * cmSigma2 * m42 
00140                        + 3. * cmOmega2 * mOmega2 * m42     - cmSigma2 * mOmega2 * m42 
00141                        + 3. * cmOmega2 * mSigma2 * m42     - cmSigma2 * mSigma2 * m42 
00142                        - 2. * mOmega2 * mSigma2 * m42);
00143       
00144       cMix_oLs = fac2 * (8. * cmOmega4                     - 4. * cmOmega2 * cmSigma2 
00145                          - 6. * cmOmega2 * mOmega2         + 2. * cmSigma2 * mOmega2 
00146                          - 6. * cmOmega2 * mSigma2         + 2. * cmSigma2 * mSigma2 
00147                          + 4. * mOmega2 * mSigma2);
00148       
00149       cMix_sLc = fac3 * (cmOmega2 * cmSigma4               - 3. * cmSigma6    
00150                          + 2. * cmSigma4 * mOmega2         + 2. * cmSigma4 * mSigma2 
00151                          - cmOmega2 * mOmega2 * mSigma2    - cmSigma2 * mOmega2 * mSigma2 
00152                          - 2. * cmOmega2 * cmSigma2 * m42  + 4. * cmSigma4 * m42 
00153                          + cmOmega2 * mOmega2 * m42        - 3. * cmSigma2 * mOmega2 * m42 
00154                          + cmOmega2 * mSigma2 * m42        - 3. * cmSigma2 * mSigma2 * m42 
00155                          + 2. * mOmega2 * mSigma2 * m42);
00156 
00157       cMix_sLs = fac3 * (4. * cmOmega2 * cmSigma2          - 8. * cmSigma4
00158                        - 2. * cmOmega2 * mOmega2           + 6. * cmSigma2 * mOmega2 
00159                        - 2. * cmOmega2 * mSigma2           + 6. * cmSigma2 * mSigma2 
00160                        - 4. * mOmega2 * mSigma2);
00161 }
00162 
00163 
00164 G4AngularDistribution::~
00165 G4AngularDistribution()
00166 { }
00167 
00168 
00169 G4double G4AngularDistribution::CosTheta(G4double S, G4double m_1, G4double m_2) const
00170 {
00171    G4double random = G4UniformRand();
00172    G4double dCosTheta = 2.;
00173    G4double cosTheta = -1.;
00174 
00175    // For jmax=12 the accuracy is better than 0.1 degree 
00176    G4int jMax = 12;
00177 
00178    for (G4int j = 1; j <= jMax; ++j)
00179    {
00180       // Accuracy is 2^-jmax 
00181       dCosTheta *= 0.5;
00182       G4double cosTh = cosTheta + dCosTheta;
00183       if(DifferentialCrossSection(S, m_1, m_2, cosTh) <= random) cosTheta = cosTh;
00184     }
00185 
00186    // Randomize in final interval in order to avoid discrete angles 
00187    cosTheta += G4UniformRand() * dCosTheta;
00188 
00189 
00190    if (cosTheta > 1. || cosTheta < -1.)
00191      throw G4HadronicException(__FILE__, __LINE__, "G4AngularDistribution::CosTheta - std::cos(theta) outside allowed range");
00192 
00193    return cosTheta;
00194 }
00195 
00196 
00197 G4double G4AngularDistribution::DifferentialCrossSection(G4double sIn, G4double m_1, G4double m_2, 
00198                                                          G4double cosTheta) const
00199 {
00200 // local calculus is in GeV, ie. normalize input
00201   sIn = sIn/sqr(GeV)+m42/2.;
00202   m_1  = m_1/GeV;
00203   m_2  = m_2/GeV;
00204 //  G4cout << "Here we go"<<sIn << " "<<m1 << " " << m2 <<" " m42<< G4endl;
00205 // scaling from masses other than p,p.
00206   G4double S = sIn - (m_1+m_2) * (m_1+m_2) + m42;
00207   G4double tMax = S - m42;
00208   G4double tp = 0.5 * (cosTheta + 1.) * tMax; 
00209   G4double twoS = 2. * S;
00210   
00211   // Define s-dependent stuff for omega-Term
00212   G4double brak1 = (twoS-m42) * (twoS-m42);
00213   G4double bOmega_3 = cOmega_3 * (-2. * cmOmega4 - 2. * cmOmega2 * twoS - brak1);
00214   G4double bOmega_2 = cOmega_2 * ( 2. * cmOmega2 * mOmega2 + sOmega1 * twoS + brak1);
00215   G4double bOmega_1 = cOmega_1 * (-4. * cmOmega2 * mOmega2 
00216                          - 2. * mOmega2*mOmega2 
00217                          - 2. * (cmOmega2 + 2 * mOmega2) * twoS 
00218                          - 3. * brak1);
00219   G4double bOmega_m = cOmega_m * (-2. * mOmega2*mOmega2 - 2. * mOmega2 * twoS - brak1);
00220   G4double bOmega_L = cOmega_L * (sOmega1 * mOmega2 + (cmOmega2 + 3. * mOmega2) * S + brak1);
00221   G4double bOmega_0 = -(bOmega_3 + bOmega_2 + bOmega_1 + bOmega_m);
00222   
00223   // Define s-dependent stuff for mix-Term            
00224   G4double bMix_o1 = cMix_o1 * (dOmega1 - twoS);
00225   G4double bMix_s1 = cMix_s1 * (dSigma1 - twoS);
00226   G4double bMix_Omega = cMix_Omega * (dOmega2 - twoS);
00227   G4double bMix_sm = cMix_sm * (dSigma2 - twoS);
00228   G4double bMix_oL = cMix_oLc + cMix_oLs * S;
00229   G4double bMix_sL = cMix_sLc + cMix_sLs * S;
00230   
00231   G4double t1_Pion = 1. / (1. + tMax / cmPion2);
00232   G4double t2_Pion = 1. + tMax / mPion2;
00233   G4double t1_Sigma = 1. / (1. + tMax / cmSigma2);
00234   G4double t2_Sigma = 1. + tMax / mSigma2;
00235   G4double t1_Omega = 1. / (1. + tMax / cmOmega2);
00236   G4double t2_Omega = 1. + tMax / mOmega2;
00237   
00238   G4double norm = Cross(t1_Pion, t1_Sigma, t1_Omega,
00239                         t2_Pion, t2_Sigma, t2_Omega,
00240                         bMix_o1, bMix_s1, bMix_Omega,
00241                         bMix_sm, bMix_oL, bMix_sL,
00242                         bOmega_0, bOmega_1, bOmega_2,
00243                         bOmega_3, bOmega_m, bOmega_L);
00244   
00245   t1_Pion = 1. / (1. + tp / cmPion2);
00246   t2_Pion = 1. + tp / mPion2;
00247   t1_Sigma = 1. / (1. + tp / cmSigma2);
00248   t2_Sigma = 1. + tp / mSigma2;
00249   t1_Omega = 1. / (1. + tp / cmOmega2);
00250   t2_Omega = 1. + tp / mOmega2;
00251   
00252   G4double dSigma;
00253   if (sym) 
00254     { 
00255       G4double to;
00256       norm = 2. * norm;
00257       to = tMax - tp;
00258       G4double t3_Pion = 1. / (1. + to / cmPion2);
00259       G4double t4_Pion = 1. + to / mPion2;
00260       G4double t3_Sigma = 1. / (1. + to / cmSigma2);
00261       G4double t4_Sigma = 1. + to / mSigma2;
00262       G4double t3_Omega = 1. / (1. + to / cmOmega2);
00263       G4double t4_Omega = 1. + to / mOmega2;
00264       
00265       dSigma = ( Cross(t1_Pion, t1_Sigma, t1_Omega,
00266                        t2_Pion,t2_Sigma, t2_Omega,
00267                        bMix_o1, bMix_s1, bMix_Omega,
00268                        bMix_sm, bMix_oL, bMix_sL,
00269                        bOmega_0, bOmega_1, bOmega_2,
00270                        bOmega_3, bOmega_m, bOmega_L) - 
00271                  Cross(t3_Pion,t3_Sigma, t3_Omega,
00272                        t4_Pion, t4_Sigma, t4_Omega,
00273                        bMix_o1, bMix_s1, bMix_Omega,
00274                        bMix_sm, bMix_oL, bMix_sL,
00275                        bOmega_0, bOmega_1, bOmega_2,
00276                        bOmega_3, bOmega_m, bOmega_L) ) 
00277         / norm + 0.5;
00278     }
00279   else
00280     {
00281       dSigma = Cross(t1_Pion, t1_Sigma, t1_Omega, 
00282                      t2_Pion, t2_Sigma, t2_Omega,
00283                      bMix_o1, bMix_s1, bMix_Omega,
00284                      bMix_sm, bMix_oL, bMix_sL,
00285                      bOmega_0, bOmega_1, bOmega_2,
00286                      bOmega_3, bOmega_m, bOmega_L) 
00287         / norm;
00288     }
00289   
00290   return dSigma;
00291 }
00292 
00293 
00294 G4double G4AngularDistribution::Cross(G4double tpPion,
00295                                       G4double tpSigma,
00296                                       G4double tpOmega,
00297                                       G4double tmPion,
00298                                       G4double tmSigma,
00299                                       G4double tmOmega,
00300                                       G4double bMix_o1,
00301                                       G4double bMix_s1,
00302                                       G4double bMix_Omega,
00303                                       G4double bMix_sm,
00304                                       G4double bMix_oL,
00305                                       G4double bMix_sL,
00306                                       G4double bOmega_0,
00307                                       G4double bOmega_1,
00308                                       G4double bOmega_2,
00309                                       G4double bOmega_3,
00310                                       G4double bOmega_m,
00311                                       G4double bOmega_L) const
00312 {
00313   G4double cross = 0;
00314      //  Pion
00315     cross += ((cPion_3 * tpPion  + cPion_2)  * tpPion  + cPion_1)  * tpPion  + cPion_m/tmPion   + cPion_0  +  cPion_L * std::log(tpPion*tmPion);
00316 //    G4cout << "cross1 "<< cross<<G4endl;
00317     //  Sigma
00318     cross += ((cSigma_3 * tpSigma + cSigma_2) * tpSigma + cSigma_1) * tpSigma + cSigma_m/tmSigma + cSigma_0 + cSigma_L * std::log(tpSigma*tmSigma);
00319 //    G4cout << "cross2 "<< cross<<G4endl;
00320     // Omega
00321     cross += ((bOmega_3 * tpOmega + bOmega_2) * tpOmega + bOmega_1) * tpOmega + bOmega_m/tmOmega + bOmega_0 + bOmega_L * std::log(tpOmega*tmOmega)
00322     // Mix
00323     +  bMix_o1 * (tpOmega - 1.)
00324     +  bMix_s1 * (tpSigma - 1.)
00325     +  bMix_Omega * std::log(tmOmega)
00326     +  bMix_sm * std::log(tmSigma)
00327     +  bMix_oL * std::log(tpOmega)
00328     +  bMix_sL * std::log(tpSigma);
00329 /*      G4cout << "cross3 "<< cross<<" "
00330              <<bMix_o1<<" "
00331              <<bMix_s1<<" "
00332              <<bMix_Omega<<" "
00333              <<bMix_sm<<" "
00334              <<bMix_oL<<" "
00335              <<bMix_sL<<" "
00336              <<tpOmega<<" "
00337              <<tpSigma<<" "
00338              <<tmOmega<<" "
00339              <<tmSigma<<" "
00340              <<tpOmega<<" "
00341              <<tpSigma
00342              <<G4endl;
00343 */
00344   return cross;
00345 
00346 }

Generated on Mon May 27 17:47:39 2013 for Geant4 by  doxygen 1.4.7