Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/im_r_matrix/src/G4XNNElasticLowE.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 ]

Diff markup

Differences between /processes/hadronic/models/im_r_matrix/src/G4XNNElasticLowE.cc (Version 11.3.0) and /processes/hadronic/models/im_r_matrix/src/G4XNNElasticLowE.cc (Version 7.0.p1)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                    <<   3 // * DISCLAIMER                                                       *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th <<   5 // * The following disclaimer summarizes all the specific disclaimers *
  6 // * the Geant4 Collaboration.  It is provided <<   6 // * of contributors to this software. The specific disclaimers,which *
  7 // * conditions of the Geant4 Software License <<   7 // * govern, are listed with their locations in:                      *
  8 // * LICENSE and available at  http://cern.ch/ <<   8 // *   http://cern.ch/geant4/license                                  *
  9 // * include a list of copyright holders.      << 
 10 // *                                                9 // *                                                                  *
 11 // * Neither the authors of this software syst     10 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     11 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     12 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     13 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  <<  14 // * use.                                                             *
 16 // * for the full disclaimer and the limitatio << 
 17 // *                                               15 // *                                                                  *
 18 // * This  code  implementation is the result  <<  16 // * This  code  implementation is the  intellectual property  of the *
 19 // * technical work of the GEANT4 collaboratio <<  17 // * GEANT4 collaboration.                                            *
 20 // * By using,  copying,  modifying or  distri <<  18 // * By copying,  distributing  or modifying the Program (or any work *
 21 // * any work based  on the software)  you  ag <<  19 // * based  on  the Program)  you indicate  your  acceptance of  this *
 22 // * use  in  resulting  scientific  publicati <<  20 // * statement, and all its terms.                                    *
 23 // * acceptance of all terms of the Geant4 Sof << 
 24 // *******************************************     21 // ********************************************************************
 25 //                                                 22 //
 26 #include "globals.hh"                              23 #include "globals.hh"
 27 #include "G4ios.hh"                                24 #include "G4ios.hh"
 28 #include "G4Exp.hh"                            << 
 29 #include "G4Log.hh"                            << 
 30 #include "G4SystemOfUnits.hh"                  << 
 31 #include "G4XNNElasticLowE.hh"                     25 #include "G4XNNElasticLowE.hh"
 32 #include "G4KineticTrack.hh"                       26 #include "G4KineticTrack.hh"
 33 #include "G4ParticleDefinition.hh"                 27 #include "G4ParticleDefinition.hh"
 34 #include "G4PhysicsVector.hh"                      28 #include "G4PhysicsVector.hh"
 35 #include "G4PhysicsLogVector.hh"               <<  29 #include "G4PhysicsLnVector.hh"
                                                   >>  30 #include "G4ParticleDefinition.hh"
 36 #include "G4Proton.hh"                             31 #include "G4Proton.hh"
 37 #include "G4Neutron.hh"                            32 #include "G4Neutron.hh"
 38                                                    33 
                                                   >>  34 
 39 const G4double G4XNNElasticLowE::_lowLimit = 0     35 const G4double G4XNNElasticLowE::_lowLimit = 0.;
 40 const G4double G4XNNElasticLowE::_highLimit =      36 const G4double G4XNNElasticLowE::_highLimit = 3.*GeV;
 41                                                    37 
 42 // Low energy limit of the cross-section table     38 // Low energy limit of the cross-section table (in GeV)
 43 // Units are assigned when filling the Physics <<  39 // Units are assigned while filling the PhysicsVector
 44 const G4double G4XNNElasticLowE::_eMinTable =      40 const G4double G4XNNElasticLowE::_eMinTable = 1.8964808;
 45 const G4double G4XNNElasticLowE::_eStepLog = 0     41 const G4double G4XNNElasticLowE::_eStepLog = 0.01;
 46                                                    42 
 47 // Cross-sections in mb                            43 // Cross-sections in mb
 48 // Units are assigned when filling the Physics <<  44 // Units are assigned while filling the PhysicsVector
 49                                                    45 
 50 const G4int G4XNNElasticLowE::tableSize = 101;     46 const G4int G4XNNElasticLowE::tableSize = 101;
 51                                                    47 
 52 const G4double G4XNNElasticLowE::ppTable[101]      48 const G4double G4XNNElasticLowE::ppTable[101] = 
 53 {                                                  49 { 
 54   60.00, //was 0.                                  50   60.00, //was 0.
 55   33.48, 26.76, 25.26, 24.55, 23.94, 23.77, 23     51   33.48, 26.76, 25.26, 24.55, 23.94, 23.77, 23.72, 23.98,
 56   25.48, 27.52, 27.72, 27.21, 25.80, 26.00, 24     52   25.48, 27.52, 27.72, 27.21, 25.80, 26.00, 24.32, 23.81,
 57   24.37, 24.36, 23.13, 22.43, 21.71, 21.01, 20     53   24.37, 24.36, 23.13, 22.43, 21.71, 21.01, 20.83, 20.74,
 58   20.25, 20.10, 20.59, 20.04, 20.83, 20.84, 21     54   20.25, 20.10, 20.59, 20.04, 20.83, 20.84, 21.07, 20.83,
 59   20.79, 21.88, 21.15, 20.92, 19.00, 18.60, 17     55   20.79, 21.88, 21.15, 20.92, 19.00, 18.60, 17.30, 17.00,
 60   16.70, 16.50, 16.20, 15.80, 15.57, 15.20, 15     56   16.70, 16.50, 16.20, 15.80, 15.57, 15.20, 15.00, 14.60,
 61   14.20, 14.00, 13.80, 13.60, 13.40, 13.20, 13     57   14.20, 14.00, 13.80, 13.60, 13.40, 13.20, 13.00, 12.85,
 62   12.70, 12.60, 12.50, 12.40, 12.30, 12.20, 12     58   12.70, 12.60, 12.50, 12.40, 12.30, 12.20, 12.10, 12.00,
 63   11.90, 11.80, 11.75, 11.70, 11.64, 11.53, 11     59   11.90, 11.80, 11.75, 11.70, 11.64, 11.53, 11.41, 11.31,
 64   11.22, 11.13, 11.05, 10.97, 10.89, 10.82, 10     60   11.22, 11.13, 11.05, 10.97, 10.89, 10.82, 10.75, 10.68,
 65   10.61, 10.54, 10.48, 10.41, 10.35, 10.28, 10     61   10.61, 10.54, 10.48, 10.41, 10.35, 10.28, 10.22, 10.16,
 66   10.13, 10.10, 10.08, 10.05, 10.02,  9.99,  9     62   10.13, 10.10, 10.08, 10.05, 10.02,  9.99,  9.96,  9.93,
 67   9.90,  9.87,  9.84,  9.80                        63   9.90,  9.87,  9.84,  9.80       
 68 };                                                 64 };
 69                                                    65 
 70 const G4double G4XNNElasticLowE::npTable[101]      66 const G4double G4XNNElasticLowE::npTable[101] = 
 71 {                                                  67 { 
 72   1500.00, // was 0.                               68   1500.00, // was 0.
 73   248.20, 93.38, 55.26, 44.50, 41.33, 38.48, 3     69   248.20, 93.38, 55.26, 44.50, 41.33, 38.48, 37.20, 35.98,
 74   35.02, 34.47, 32.48, 30.76, 29.46, 28.53, 27     70   35.02, 34.47, 32.48, 30.76, 29.46, 28.53, 27.84, 27.20,
 75   26.53, 25.95, 25.59, 25.46, 25.00, 24.49, 24     71   26.53, 25.95, 25.59, 25.46, 25.00, 24.49, 24.08, 23.86,
 76   23.17, 22.70, 21.88, 21.48, 20.22, 19.75, 18     72   23.17, 22.70, 21.88, 21.48, 20.22, 19.75, 18.97, 18.39,
 77   17.98, 17.63, 17.21, 16.72, 16.68, 16.58, 16     73   17.98, 17.63, 17.21, 16.72, 16.68, 16.58, 16.42, 16.22,
 78   15.98, 15.71, 15.42, 15.14, 14.87, 14.65, 14     74   15.98, 15.71, 15.42, 15.14, 14.87, 14.65, 14.44, 14.26,
 79   14.10, 13.95, 13.80, 13.64, 13.47, 13.29, 13     75   14.10, 13.95, 13.80, 13.64, 13.47, 13.29, 13.09, 12.89,
 80   12.68, 12.47, 12.27, 12.06, 11.84, 11.76, 11     76   12.68, 12.47, 12.27, 12.06, 11.84, 11.76, 11.69, 11.60,
 81   11.50, 11.41, 11.29, 11.17, 11.06, 10.93, 10     77   11.50, 11.41, 11.29, 11.17, 11.06, 10.93, 10.81, 10.68,
 82   10.56, 10.44, 10.33, 10.21, 10.12, 10.03,  9     78   10.56, 10.44, 10.33, 10.21, 10.12, 10.03,  9.96,  9.89,
 83   9.83,  9.80,  9.77,  9.75,  9.74,  9.74,  9.     79   9.83,  9.80,  9.77,  9.75,  9.74,  9.74,  9.74,  9.76,
 84   9.73,  9.70,  9.68,  9.65,  9.63,  9.60,  9.     80   9.73,  9.70,  9.68,  9.65,  9.63,  9.60,  9.57,  9.55,
 85   9.52,  9.49,  9.46,  9.43                        81   9.52,  9.49,  9.46,  9.43
 86 };                                                 82 };
 87                                                    83 
 88                                                    84 
 89 G4XNNElasticLowE::G4XNNElasticLowE()               85 G4XNNElasticLowE::G4XNNElasticLowE() 
 90 {                                                  86 { 
 91   // Cross-sections are available in the range     87   // Cross-sections are available in the range (_eMin,_eMax)
 92                                                    88 
 93   _eMin = _eMinTable * GeV;                        89   _eMin = _eMinTable * GeV;
 94   _eMax = G4Exp(G4Log(_eMinTable) + tableSize  <<  90   _eMax = std::exp(std::log(_eMinTable) + tableSize * _eStepLog) * GeV;
 95   if (_eMin < _lowLimit)                           91   if (_eMin < _lowLimit)
 96     throw G4HadronicException(__FILE__, __LINE     92     throw G4HadronicException(__FILE__, __LINE__, "G4XNNElasticLowE::G4XNNElasticLowE - Low energy limit not valid");    
 97   if (_highLimit > _eMax)                          93   if (_highLimit > _eMax)
 98     throw G4HadronicException(__FILE__, __LINE     94     throw G4HadronicException(__FILE__, __LINE__, "G4XNNElasticLowE::G4XNNElasticLowE - High energy limit not valid");    
 99   G4PhysicsVector* pp = new G4PhysicsLogVector <<  95   G4PhysicsVector* pp = new G4PhysicsLnVector(_eMin,_eMax,tableSize);
100                                                    96 
101   _eMin = G4Exp(G4Log(_eMinTable)-_eStepLog)*G <<  97   _eMin = std::exp(std::log(_eMinTable)-_eStepLog)*GeV;
102   if (_eMin < _lowLimit)                           98   if (_eMin < _lowLimit)
103     throw G4HadronicException(__FILE__, __LINE     99     throw G4HadronicException(__FILE__, __LINE__, "G4XNNElasticLowE::G4XNNElasticLowE - Low energy limit not valid");
104   G4PhysicsVector* np = new G4PhysicsLogVector << 100   G4PhysicsVector* np = new G4PhysicsLnVector(_eMin,_eMax,tableSize);
105                                                   101 
106   G4int i;                                        102   G4int i;
107   for (i=0; i<tableSize; i++)                     103   for (i=0; i<tableSize; i++)
108     {                                             104     {
109       G4double value = ppTable[i] * millibarn;    105       G4double value = ppTable[i] * millibarn;
110       pp->PutValue(i,value);                      106       pp->PutValue(i,value);
111       value = npTable[i] * millibarn;             107       value = npTable[i] * millibarn;
112       np->PutValue(i,value);                      108       np->PutValue(i,value);
113     }                                             109     }
114   xMap[G4Proton::ProtonDefinition()] = pp;     << 110   xMap[G4Proton::ProtonDefinition()->GetParticleName()] = pp;
115   xMap[G4Neutron::NeutronDefinition()] = np;   << 111   xMap[G4Neutron::NeutronDefinition()->GetParticleName()] = np;
116 }                                                 112 }
117                                                   113 
118                                                   114 
119 G4XNNElasticLowE::~G4XNNElasticLowE()             115 G4XNNElasticLowE::~G4XNNElasticLowE()
120 {                                                 116 {
121   delete xMap[G4Proton::ProtonDefinition()];   << 117   delete xMap[G4Proton::ProtonDefinition()->GetParticleName()];
122   delete xMap[G4Neutron::NeutronDefinition()]; << 118   delete xMap[G4Neutron::NeutronDefinition()->GetParticleName()];
123 }                                                 119 }
124                                                   120 
125                                                   121 
126 G4bool G4XNNElasticLowE::operator==(const G4XN    122 G4bool G4XNNElasticLowE::operator==(const G4XNNElasticLowE &right) const
127 {                                                 123 {
128   return (this == (G4XNNElasticLowE *) &right)    124   return (this == (G4XNNElasticLowE *) &right);
129 }                                                 125 }
130                                                   126 
131                                                   127 
132 G4bool G4XNNElasticLowE::operator!=(const G4XN    128 G4bool G4XNNElasticLowE::operator!=(const G4XNNElasticLowE &right) const
133 {                                                 129 {
134                                                   130 
135   return (this != (G4XNNElasticLowE *) &right)    131   return (this != (G4XNNElasticLowE *) &right);
136 }                                                 132 }
137                                                   133 
138                                                   134 
139 G4double G4XNNElasticLowE::CrossSection(const     135 G4double G4XNNElasticLowE::CrossSection(const G4KineticTrack& trk1, const G4KineticTrack& trk2) const  
140 {                                                 136 {
141   G4double sigma = 0.;                            137   G4double sigma = 0.;
142   G4double sqrtS = (trk1.Get4Momentum() + trk2    138   G4double sqrtS = (trk1.Get4Momentum() + trk2.Get4Momentum()).mag();
143   G4bool dummy = false;                           139   G4bool dummy = false;
144                                                   140 
145   const G4ParticleDefinition * key = FindKeyPa << 141   G4String key = FindKeyParticle(trk1,trk2);
146                                                   142 
147   typedef std::map <const G4ParticleDefinition << 143   typedef std::map <G4String, G4PhysicsVector*, std::less<G4String> > StringPhysMap;
148                                                   144 
149   if (xMap.find(key)!= xMap.end())                145   if (xMap.find(key)!= xMap.end())
150     {                                             146     {
151                                                   147 
152       StringPhysMap::const_iterator iter;         148       StringPhysMap::const_iterator iter;
153       for (iter = xMap.begin(); iter != xMap.e    149       for (iter = xMap.begin(); iter != xMap.end(); ++iter)
154   {                                               150   {
155     const G4ParticleDefinition * str = (*iter) << 151     G4String str = (*iter).first;
156           if (str == key)                         152           if (str == key)
157       {                                           153       {
158         G4PhysicsVector* physVector = (*iter).    154         G4PhysicsVector* physVector = (*iter).second; 
159         //     G4PhysicsVector* physVector = x    155         //     G4PhysicsVector* physVector = xMap[key];
160         if (sqrtS >= _eMin && sqrtS <= _eMax)     156         if (sqrtS >= _eMin && sqrtS <= _eMax)
161     {                                             157     {
162       sigma = physVector->GetValue(sqrtS,dummy    158       sigma = physVector->GetValue(sqrtS,dummy);
163     } else if ( sqrtS < _eMin )                   159     } else if ( sqrtS < _eMin )
164                 {                                 160                 {
165                   sigma = physVector->GetValue    161                   sigma = physVector->GetValue(_eMin,dummy);
166                 }                                 162                 }
167     //G4cout << " sqrtS / sigma " << sqrtS/GeV << 
168     //          sigma/millibarn << G4endl;     << 
169       }                                           163       }
170   }                                               164   }
171     }                                             165     }
172   return sigma;                                   166   return sigma;
173 }                                                 167 }
174                                                   168 
175                                                   169 
176 void G4XNNElasticLowE::Print() const              170 void G4XNNElasticLowE::Print() const
177 {                                                 171 {
178   // Dump the pp cross-section table              172   // Dump the pp cross-section table
179                                                   173 
180   G4cout << Name() << ", pp cross-section: " <    174   G4cout << Name() << ", pp cross-section: " << G4endl;
181                                                   175 
182   G4bool dummy = false;                           176   G4bool dummy = false;
183   G4int i;                                        177   G4int i;
184   const G4ParticleDefinition * key = G4Proton: << 178   G4String key = G4Proton::ProtonDefinition()->GetParticleName();
185   G4PhysicsVector* pp = 0;                        179   G4PhysicsVector* pp = 0;
186                                                   180 
187   typedef std::map <const G4ParticleDefinition << 181   typedef std::map <G4String, G4PhysicsVector*, std::less<G4String> > StringPhysMap;
188   StringPhysMap::const_iterator iter;             182   StringPhysMap::const_iterator iter;
189                                                   183 
190   for (iter = xMap.begin(); iter != xMap.end()    184   for (iter = xMap.begin(); iter != xMap.end(); ++iter)
191     {                                             185     {
192       const G4ParticleDefinition * str = (*ite << 186       G4String str = (*iter).first;
193       if (str == key)                             187       if (str == key)
194   {                                               188   {
195     pp = (*iter).second;                          189     pp = (*iter).second; 
196   }                                               190   }
197     }                                             191     }
198                                                   192   
199   if (pp != 0)                                    193   if (pp != 0)
200     {                                             194     { 
201       for (i=0; i<tableSize; i++)                 195       for (i=0; i<tableSize; i++)
202   {                                               196   {
203     G4double e = pp->GetLowEdgeEnergy(i);         197     G4double e = pp->GetLowEdgeEnergy(i);
204     G4double sigma = pp->GetValue(e,dummy) / m    198     G4double sigma = pp->GetValue(e,dummy) / millibarn;
205     G4cout << i << ") e = " << e / GeV << " Ge    199     G4cout << i << ") e = " << e / GeV << " GeV ---- Cross section = " << sigma << " mb " << G4endl;
206   }                                               200   }
207     }                                             201     }
208                                                   202   
209   // Dump the np cross-section table              203   // Dump the np cross-section table
210                                                   204 
211   G4cout << Name() << ", np cross-section: " <    205   G4cout << Name() << ", np cross-section: " << G4endl;
212                                                   206 
213   key = G4Neutron::NeutronDefinition();        << 207   key = G4Neutron::NeutronDefinition()->GetParticleName();
214   G4PhysicsVector* np = 0;                        208   G4PhysicsVector* np = 0;
215   for (iter = xMap.begin(); iter != xMap.end()    209   for (iter = xMap.begin(); iter != xMap.end(); ++iter)
216     {                                             210     {
217       const G4ParticleDefinition * str = (*ite << 211       G4String str = (*iter).first;
218       if (str == key)                             212       if (str == key)
219   {                                               213   {
220     np = (*iter).second;                          214     np = (*iter).second; 
221   }                                               215   }
222     }                                             216     }
223                                                   217   
224   //  G4PhysicsVector* np = xMap[G4Neutron::Ne    218   //  G4PhysicsVector* np = xMap[G4Neutron::NeutronDefinition()->GetParticleName()];
225                                                   219   
226   if (np != 0)                                    220   if (np != 0)
227     {                                             221     { 
228       for (i=0; i<tableSize; i++)                 222       for (i=0; i<tableSize; i++)
229   {                                               223   {
230     G4double e = np->GetLowEdgeEnergy(i);         224     G4double e = np->GetLowEdgeEnergy(i);
231     G4double sigma = np->GetValue(e,dummy) / m    225     G4double sigma = np->GetValue(e,dummy) / millibarn;
232     G4cout << i << ") e = " << e / GeV << " Ge    226     G4cout << i << ") e = " << e / GeV << " GeV ---- Cross section = " << sigma << " mb " << G4endl;
233   }                                               227   }
234     }                                             228     }
235   G4VCrossSectionSource::Print();                 229   G4VCrossSectionSource::Print();
236 }                                                 230 }
237                                                   231 
238                                                   232 
239 G4String G4XNNElasticLowE::Name() const           233 G4String G4XNNElasticLowE::Name() const
240 {                                                 234 {
241   G4String name("NNElasticLowE");                 235   G4String name("NNElasticLowE");
242   return name;                                    236   return name;
243 }                                                 237 }
244                                                   238 
245                                                   239 
246                                                   240 
247 G4bool G4XNNElasticLowE::IsValid(G4double e) c    241 G4bool G4XNNElasticLowE::IsValid(G4double e) const
248 {                                                 242 {
249   G4bool answer = InLimits(e,_lowLimit,_highLi    243   G4bool answer = InLimits(e,_lowLimit,_highLimit);
250                                                   244 
251   return answer;                                  245   return answer;
252 }                                                 246 }
253                                                   247 
254                                                   248 
255                                                   249