Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/externals/clhep/src/RandBreitWigner.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 /externals/clhep/src/RandBreitWigner.cc (Version 11.3.0) and /externals/clhep/src/RandBreitWigner.cc (Version 9.1)


  1 // -*- C++ -*-                                      1 
  2 //                                                
  3 // -------------------------------------------    
  4 //                             HEP Random         
  5 //                       --- RandBreitWigner -    
  6 //                      class implementation f    
  7 // -------------------------------------------    
  8 // This file is part of Geant4 (simulation too    
  9                                                   
 10 // ===========================================    
 11 // Gabriele Cosmo - Created: 5th September 199    
 12 //                - Added methods to shoot arr    
 13 // J.Marraffino   - Added default arguments as    
 14 //                  operator() with arguments:    
 15 // M Fischler     - put and get to/from stream    
 16 // M Fischler       - put/get to/from streams     
 17 //      + storing doubles avoid problems with     
 18 //      4/14/05                                   
 19 // ===========================================    
 20                                                   
 21 #include "CLHEP/Random/RandBreitWigner.h"         
 22 #include "CLHEP/Units/PhysicalConstants.h"        
 23 #include "CLHEP/Random/DoubConv.h"                
 24 #include <algorithm>  // for min() and max()      
 25 #include <cmath>                                  
 26 #include <iostream>                               
 27 #include <string>                                 
 28 #include <vector>                                 
 29                                                   
 30 namespace CLHEP {                                 
 31                                                   
 32 std::string RandBreitWigner::name() const {ret    
 33 HepRandomEngine & RandBreitWigner::engine() {r    
 34                                                   
 35 RandBreitWigner::~RandBreitWigner() {             
 36 }                                                 
 37                                                   
 38 double RandBreitWigner::operator()() {            
 39    return fire( defaultA, defaultB );             
 40 }                                                 
 41                                                   
 42 double RandBreitWigner::operator()( double a,     
 43    return fire( a, b );                           
 44 }                                                 
 45                                                   
 46 double RandBreitWigner::operator()( double a,     
 47    return fire( a, b, c );                        
 48 }                                                 
 49                                                   
 50 double RandBreitWigner::shoot(double mean, dou    
 51 {                                                 
 52    double rval, displ;                            
 53                                                   
 54    rval = 2.0*HepRandom::getTheEngine()->flat(    
 55    displ = 0.5*gamma*std::tan(rval*CLHEP::half    
 56                                                   
 57    return mean + displ;                           
 58 }                                                 
 59                                                   
 60 double RandBreitWigner::shoot(double mean, dou    
 61 {                                                 
 62    double val, rval, displ;                       
 63                                                   
 64    if ( gamma == 0.0 ) return mean;               
 65    val = std::atan(2.0*cut/gamma);                
 66    rval = 2.0*HepRandom::getTheEngine()->flat(    
 67    displ = 0.5*gamma*std::tan(rval*val);          
 68                                                   
 69    return mean + displ;                           
 70 }                                                 
 71                                                   
 72 double RandBreitWigner::shootM2(double mean, d    
 73 {                                                 
 74    double val, rval, displ;                       
 75                                                   
 76    if ( gamma == 0.0 ) return mean;               
 77    val = std::atan(-mean/gamma);                  
 78    rval = RandFlat::shoot(val, CLHEP::halfpi);    
 79    displ = gamma*std::tan(rval);                  
 80                                                   
 81    return std::sqrt(mean*mean + mean*displ);      
 82 }                                                 
 83                                                   
 84 double RandBreitWigner::shootM2(double mean, d    
 85 {                                                 
 86    double rval, displ;                            
 87    double lower, upper, tmp;                      
 88                                                   
 89    if ( gamma == 0.0 ) return mean;               
 90    tmp = std::max(0.0,(mean-cut));                
 91    lower = std::atan( (tmp*tmp-mean*mean)/(mea    
 92    upper = std::atan( ((mean+cut)*(mean+cut)-m    
 93    rval = RandFlat::shoot(lower, upper);          
 94    displ = gamma*std::tan(rval);                  
 95                                                   
 96    return std::sqrt(std::max(0.0, mean*mean +     
 97 }                                                 
 98                                                   
 99 void RandBreitWigner::shootArray ( const int s    
100 {                                                 
101   for( double* v = vect; v != vect + size; ++v    
102     *v = shoot( 1.0, 0.2 );                       
103 }                                                 
104                                                   
105 void RandBreitWigner::shootArray ( const int s    
106                                    double a, d    
107 {                                                 
108   for( double* v = vect; v != vect + size; ++v    
109     *v = shoot( a, b );                           
110 }                                                 
111                                                   
112 void RandBreitWigner::shootArray ( const int s    
113                                    double a, d    
114                                    double c )     
115 {                                                 
116   for( double* v = vect; v != vect + size; ++v    
117     *v = shoot( a, b, c );                        
118 }                                                 
119                                                   
120 //----------------                                
121                                                   
122 double RandBreitWigner::shoot(HepRandomEngine*    
123                                  double mean,     
124 {                                                 
125    double rval, displ;                            
126                                                   
127    rval = 2.0*anEngine->flat()-1.0;               
128    displ = 0.5*gamma*std::tan(rval*CLHEP::half    
129                                                   
130    return mean + displ;                           
131 }                                                 
132                                                   
133 double RandBreitWigner::shoot(HepRandomEngine*    
134                                  double mean,     
135 {                                                 
136    double val, rval, displ;                       
137                                                   
138    if ( gamma == 0.0 ) return mean;               
139    val = std::atan(2.0*cut/gamma);                
140    rval = 2.0*anEngine->flat()-1.0;               
141    displ = 0.5*gamma*std::tan(rval*val);          
142                                                   
143    return mean + displ;                           
144 }                                                 
145                                                   
146 double RandBreitWigner::shootM2(HepRandomEngin    
147                                    double mean    
148 {                                                 
149    double val, rval, displ;                       
150                                                   
151    if ( gamma == 0.0 ) return mean;               
152    val = std::atan(-mean/gamma);                  
153    rval = RandFlat::shoot(anEngine,val, CLHEP:    
154    displ = gamma*std::tan(rval);                  
155                                                   
156    return std::sqrt(mean*mean + mean*displ);      
157 }                                                 
158                                                   
159 double RandBreitWigner::shootM2(HepRandomEngin    
160                                    double mean    
161 {                                                 
162    double rval, displ;                            
163    double lower, upper, tmp;                      
164                                                   
165    if ( gamma == 0.0 ) return mean;               
166    tmp = std::max(0.0,(mean-cut));                
167    lower = std::atan( (tmp*tmp-mean*mean)/(mea    
168    upper = std::atan( ((mean+cut)*(mean+cut)-m    
169    rval = RandFlat::shoot(anEngine, lower, upp    
170    displ = gamma*std::tan(rval);                  
171                                                   
172    return std::sqrt( std::max(0.0, mean*mean+m    
173 }                                                 
174                                                   
175 void RandBreitWigner::shootArray ( HepRandomEn    
176                                    const int s    
177 {                                                 
178   for( double* v = vect; v != vect + size; ++v    
179     *v = shoot( anEngine, 1.0, 0.2 );             
180 }                                                 
181                                                   
182 void RandBreitWigner::shootArray ( HepRandomEn    
183                                    const int s    
184                                    double a, d    
185 {                                                 
186   for( double* v = vect; v != vect + size; ++v    
187     *v = shoot( anEngine, a, b );                 
188 }                                                 
189                                                   
190 void RandBreitWigner::shootArray ( HepRandomEn    
191                                    const int s    
192                                    double a, d    
193 {                                                 
194   for( double* v = vect; v != vect + size; ++v    
195     *v = shoot( anEngine, a, b, c );              
196 }                                                 
197                                                   
198 //----------------                                
199                                                   
200 double RandBreitWigner::fire()                    
201 {                                                 
202   return fire( defaultA, defaultB );              
203 }                                                 
204                                                   
205 double RandBreitWigner::fire(double mean, doub    
206 {                                                 
207    double rval, displ;                            
208                                                   
209    rval = 2.0*localEngine->flat()-1.0;            
210    displ = 0.5*gamma*std::tan(rval*CLHEP::half    
211                                                   
212    return mean + displ;                           
213 }                                                 
214                                                   
215 double RandBreitWigner::fire(double mean, doub    
216 {                                                 
217    double val, rval, displ;                       
218                                                   
219    if ( gamma == 0.0 ) return mean;               
220    val = std::atan(2.0*cut/gamma);                
221    rval = 2.0*localEngine->flat()-1.0;            
222    displ = 0.5*gamma*std::tan(rval*val);          
223                                                   
224    return mean + displ;                           
225 }                                                 
226                                                   
227 double RandBreitWigner::fireM2()                  
228 {                                                 
229   return fireM2( defaultA, defaultB );            
230 }                                                 
231                                                   
232 double RandBreitWigner::fireM2(double mean, do    
233 {                                                 
234    double val, rval, displ;                       
235                                                   
236    if ( gamma == 0.0 ) return mean;               
237    val = std::atan(-mean/gamma);                  
238    rval = RandFlat::shoot(localEngine.get(),va    
239    displ = gamma*std::tan(rval);                  
240                                                   
241    return std::sqrt(mean*mean + mean*displ);      
242 }                                                 
243                                                   
244 double RandBreitWigner::fireM2(double mean, do    
245 {                                                 
246    double rval, displ;                            
247    double lower, upper, tmp;                      
248                                                   
249    if ( gamma == 0.0 ) return mean;               
250    tmp = std::max(0.0,(mean-cut));                
251    lower = std::atan( (tmp*tmp-mean*mean)/(mea    
252    upper = std::atan( ((mean+cut)*(mean+cut)-m    
253    rval = RandFlat::shoot(localEngine.get(),lo    
254    displ = gamma*std::tan(rval);                  
255                                                   
256    return std::sqrt(std::max(0.0, mean*mean +     
257 }                                                 
258                                                   
259 void RandBreitWigner::fireArray ( const int si    
260 {                                                 
261   for( double* v = vect; v != vect + size; ++v    
262     *v = fire(defaultA, defaultB );               
263 }                                                 
264                                                   
265 void RandBreitWigner::fireArray ( const int si    
266                                   double a, do    
267 {                                                 
268   for( double* v = vect; v != vect + size; ++v    
269     *v = fire( a, b );                            
270 }                                                 
271                                                   
272 void RandBreitWigner::fireArray ( const int si    
273                                   double a, do    
274 {                                                 
275   for( double* v = vect; v != vect + size; ++v    
276     *v = fire( a, b, c );                         
277 }                                                 
278                                                   
279                                                   
280 std::ostream & RandBreitWigner::put ( std::ost    
281   long pr=os.precision(20);                       
282   std::vector<unsigned long> t(2);                
283   os << " " << name() << "\n";                    
284   os << "Uvec" << "\n";                           
285   t = DoubConv::dto2longs(defaultA);              
286   os << defaultA << " " << t[0] << " " << t[1]    
287   t = DoubConv::dto2longs(defaultB);              
288   os << defaultB << " " << t[0] << " " << t[1]    
289   os.precision(pr);                               
290   return os;                                      
291 }                                                 
292                                                   
293 std::istream & RandBreitWigner::get ( std::ist    
294   std::string inName;                             
295   is >> inName;                                   
296   if (inName != name()) {                         
297     is.clear(std::ios::badbit | is.rdstate());    
298     std::cerr << "Mismatch when expecting to r    
299             << name() << " distribution\n"        
300         << "Name found was " << inName            
301         << "\nistream is left in the badbit st    
302     return is;                                    
303   }                                               
304   if (possibleKeywordInput(is, "Uvec", default    
305     std::vector<unsigned long> t(2);              
306     is >> defaultA >> t[0] >> t[1]; defaultA =    
307     is >> defaultB >> t[0] >> t[1]; defaultB =    
308     return is;                                    
309   }                                               
310   // is >> defaultA encompassed by possibleKey    
311   is >> defaultB;                                 
312   return is;                                      
313 }                                                 
314                                                   
315                                                   
316 }  // namespace CLHEP                             
317                                                   
318