Jpp  18.4.0
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JScatteringModel.hh
Go to the documentation of this file.
1 #ifndef H_J_SCATTERING_MODEL
2 #define H_J_SCATTERING_MODEL
3 
4 /**
5  * \author mjongen
6  */
7 
8 // c++ standard library includes
9 #include<iostream>
10 #include<iomanip>
11 #include<cmath>
12 #include<cstdlib>
13 
14 // root includes
15 #include "TObject.h"
16 #include "TH1F.h"
17 #include "TH2F.h"
18 #include "TRandom3.h"
19 
20 // JPP includes
21 // these include files are hidden from rootcint
22 #include "JROOT/JRoot.hh"
23 #include "JMarkov/JPhotonPath.hh"
24 #include "JGeometry3D/JAxis3D.hh"
26 
27 #ifndef __CINT__
28  #include "JPhysics/KM3NeT.hh"
29 #endif
30 
31 namespace JMARKOV {
32 
33  // for rootcint, all we need is this class declaration
34  #ifdef __CINT__
35  class JPhotonPath ;
36  #endif
37 }
38 
39 namespace JMARKOV {
40  using namespace JGEOMETRY3D ;
41 
42  /**
43  Virtual base class for a scattering model.
44  **/
46  public :
47 
48  // constructor
50 
51  // destructor
52  virtual ~JScatteringModel() {}
53 
54  /**
55  Return the probability density as a function of cos(theta)
56 
57  dP / dOmega = dP / dcosTheta dPhi
58 
59  to scatter in a given direction. Theta=0 is forward scattering.
60  **/
61  virtual double getScatteringProbability( double ct ) = 0 ;
62 
63  /**
64  Return a randomly generated direction according to the scattering
65  probability distribution.
66 
67  This uses gRandom.
68  **/
69  virtual JVersor3D generateDirection() = 0 ;
70 
71  // return the scattering length in [m]
72  virtual double getScatteringLength() { return lambda_scat ; }
73 
74  protected :
75 
76  double lambda_scat ; // scattering length in [m]
77  } ;
78 
79  /**
80  Virtual base class for a light source.
81  **/
82  class JSourceModel {
83  public:
84 
86 
87  virtual ~JSourceModel() {}
88 
89  /**
90  Return the probability density
91 
92  dP / dOmega = dP / dCosTheta dPhi
93 
94  that a photon from this source is emitted in a given direction,
95  given that a photon is emitted.
96  **/
97  virtual double getEmissionProbability( JVersor3D dir ) = 0 ;
98 
99  /**
100  Return a randomly generated direction according to the emission
101  distribution.
102 
103  This uses gRandom.
104  **/
105  virtual JVersor3D generateDirection() = 0 ;
106 
107  void setPosition( JPosition3D& _pos ) { pos = _pos ; }
108 
109  const JPosition3D& getPosition() const { return pos ; }
110 
111  protected:
112 
114  } ;
115 
116  /**
117  Implementation of the JSourceModel class that represents an
118  isotropic source.
119  **/
121  public :
122 
124 
126  return 1.0/(4.0*M_PI) ;
127  }
128 
130  double dx, dy, dz ;
131  gRandom->Sphere(dx,dy,dz,1) ;
132  return JVersor3D(dx,dy,dz) ;
133  }
134  } ;
135 
136  /**
137  Implementation of the JSourceModel class that represents a
138  directed source with a flat intensity distribution.
139 
140  When _cutSurface = true, photons that would be emitted into
141  the surface (defined by the surface normal) will not be emitted.
142  **/
143  class JDirectedSource : public JSourceModel {
144  public :
145 
146  /**
147  Constructor.
148  _source_dir is the source direction
149  _alpha is the opening angle in radians
150 
151  note that _alpha should be > -1
152  **/
153  JDirectedSource( JVersor3D _source_dir, double _alpha, bool _cutSurface=false, const JVersor3D _surface_normal=JVersor3D() ) :
154  source_dir(_source_dir), ctmin(cos(0.5*_alpha)), alpha(_alpha), cutSurface(_cutSurface), surface_normal(_surface_normal) {
155  // we have to find out which fraction of photons would be emitted into the surface, given
156  // the direction and opening angle
157  // as far as I can tell, this is a non-trivial integral, so I solve it numerically
158 
159  // delta is the angle of the beam direction over the surface
160  // (negative means the beam points into the surface)
161  delta = 0.5*M_PI - acos(surface_normal.getDot(source_dir)) ;
162  omega = 2.0*M_PI*(1.0-ctmin) ; // space angle of beam ("spherical cap")
163 
164  if( cutSurface ) {
165  if( delta>0.0 && delta>0.5*_alpha ) {
166  // no part of the beam points into the surface
167  //fraction_into_surface = 0.0 ;
168  } else if( delta<=0 && fabs(delta)>=0.5*alpha ) {
169  // the beam points completely into the surface, throw an error
170  //fraction_into_surface = 1.0 ;
171  omega = 0.0 ;
172  cerr << "FATAL ERROR in JDirectedSource: beam points completely into the surface." << endl ;
173  exit(1) ;
174  } else {
175  // part of the beam points into the surface, which effectively reduces
176  // the space angle omega
177  // here we numerically compute omega
178 
179  // only at theta=delta to alpha/2 is part of the beam absorbed
180  // we calculate the space angle of the absorbed beam part
181  double omega_absorbed = 0.0 ;
182  const int __nct = 100000 ;
183  const double __ctmin = cos(0.5*alpha) ;
184  const double __ctmax = cos(delta) ;
185  const double __dct = (__ctmax - __ctmin)/__nct ;
186  for( int i=0; i<__nct; ++i ) {
187  double __ct = __ctmin + (i+0.5) * __dct ;
188  double __theta = acos(__ct) ;
189  double __phirange = 2.0*(M_PI-getPhiMax(__theta)) ;
190  omega_absorbed += __phirange ;
191  }
192  omega_absorbed *= __dct ;
193  omega -= omega_absorbed ;
194 
195  // if the beam direction is into the surface, this part gets absorbed completely
196  if( delta<0 ) {
197  omega -= 2.0 * M_PI * (1.0-cos(delta)) ;
198  }
199  }
200  }
201  }
202 
203  /// return the opening angle in radians
204  double getOpeningAngle() const { return alpha ; }
205 
206  /// return the source direction
207  JVersor3D getDirection() const { return source_dir ; }
208 
210  double ct = dir.getDot(source_dir) ;
211  if( ct>ctmin ) {
212  if( cutSurface && dir.getDot(surface_normal)<0 ) return 0.0 ;
213  return 1.0/omega ;
214  }
215  return 0.0 ;
216  }
217 
219  while(true) {
220  double phi = gRandom->Uniform(0,2*M_PI) ;
221  double ct = gRandom->Uniform(ctmin,1) ;
222  JVersor3D _dir( JAngle3D(acos(ct),phi) ) ;
223  JRotation3D R( source_dir ) ;
224  JVersor3D dir(_dir.rotate_back(R)) ;
225  if( !cutSurface ) return dir ;
226  if( dir.getDot(surface_normal)>0 ) return dir ;
227  }
228  }
229 
230  protected :
231 
232  /**
233  For a given angle theta with respect to the beam direction, and a given surface
234  normal, we define phi as the rotation angle around the beam, where phi=0
235  corresponds to the direction sticking out over the surface a much as possible.
236 
237  This function returns 'phi max', which is the largest value of phi for which the
238  direction does not point into the surface.
239 
240  **/
241  double getPhiMax( double theta ) {
242  if( delta>0 ) {
243  if( theta<delta ) return M_PI ;
244  if( theta+delta > M_PI ) return 0.0 ;
245  }
246  if( delta<0 ) {
247  if( theta<fabs(delta) ) return 0.0 ;
248  if( theta+fabs(delta) > M_PI ) return M_PI ;
249  }
250  double phimax = acos(-tan(delta)/tan(theta)) ;
251  return phimax ;
252  }
253 
254 
256  double ctmin ;
257  double alpha ;
258  bool cutSurface ;
260  double delta ;
261  //double fraction_into_surface ; // the fraction of photons that would be emitted into the surface
262  double omega ; // integrated phase angle of photon emission directions
263  } ;
264 
265 
266  /**
267  Virtual base class for a light detector ("photon target").
268  **/
269  class JTargetModel {
270  public:
271 
273  r = 0 ; // infinitely small target by default
274  }
275 
276  virtual ~JTargetModel() {}
277 
278  /**
279  Return the efficiency, which is defined as the probability that
280  a photon with the given direction hitting the target will be
281  registered.
282 
283  Note that we assume by convention that the direction is the
284  PHOTON direction, NOT the direction that you would see the
285  photon coming from!
286 
287  By convention the highest efficiency is 1.
288  **/
289  virtual double getEfficiency( JVersor3D dir ) const = 0 ;
290 
291  /**
292  Return the effective area, i.e. the integral over the whole
293  surface of the target, weighted by the efficiency.
294  **/
295  virtual double getEffectiveArea() { return 0.0 ; }
296 
297  // set the target radius in [m]
298  void setRadius( double _r ) { r = _r ; }
299 
300  // get the target radius in [m]
301  double getRadius() const { return r ; }
302 
303  void setPosition( JPosition3D& _pos ) { axis = JAxis3D(_pos,axis.getDirection()) ; }
304  const JPosition3D& getPosition() const { return axis.getPosition() ; }
305 
306  void setDirection( JVersor3D& _dir ) { axis = JAxis3D(axis.getPosition(),_dir) ; }
307  const JVersor3D& getDirection() const { return axis.getDirection() ; }
308 
309  protected:
310 
311  double r ;
313  } ;
314 
315  /**
316  Implementation of the JTargetModel class that represents a
317  single PMT on a DOM
318  **/
319  class JPMTTarget : public JTargetModel {
320  public :
321 
322  /**
323  Constructor
324 
325  _r is DOM radius, default = 0.2159 [m] (=17 inch diameter glass sphere)
326 
327  _alpha is the opening angle in radians
328  **/
329  JPMTTarget( JAxis3D _axis, double _alpha, double _r=0.2159 ) {
330  setRadius(_r) ;
331  axis = _axis ;
332  alpha = _alpha ;
333 
334  // minimal dot product of PMT direction vs hit direction on DOM
335  ctmin = cos(0.5*_alpha) ;
336  }
337 
338  /// return the PMT opening angle in radians
339  double getOpeningAngle() const { return alpha ; }
340 
341  double getEfficiency( JVersor3D dir ) const {
342  const double ct = axis.getDirection().getDot(dir) ;
343  if( ct > ctmin ) {
344  return 1.0 ;
345  } else {
346  return 0.0 ;
347  }
348  }
349 
350  double getEffectiveArea() const {
351  return 2.0*r*r*M_PI*(1-ctmin) ;
352  }
353 
354  double ctmin ;
355 
356  protected :
357 
358  double alpha ;
359  } ;
360 
361  /**
362  Implementation of the JTargetModel class that represents a
363  spherically symmetric target
364  **/
366  public :
367 
369 
370  double getEfficiency( JVersor3D dir ) const { return 1.0 ; }
371  } ;
372 
373  /**
374  Implementation of the JTargetModel class that represents a
375  directed target with a cos(theta) acceptance.
376  **/
377  class JCosineTarget : public JTargetModel {
378  public :
379 
380  /**
381  Constructor.
382  _target_dir is the direction where the photon has
383  the highest chance to be detected.
384  **/
385  JCosineTarget( JVersor3D _target_dir ) :
386  target_dir(_target_dir) {}
387 
388  double getEfficiency( JVersor3D dir ) const {
389  return max( dir.getDot(target_dir), 0.0 ) ;
390  }
391 
392  protected :
393 
395  } ;
396 
397 
398  /**
399  Implementation of the JScatteringModel interface with scattering
400  according to the Henyey-Greenstein function.
401  **/
403 
404  public:
405 
406  /**
407  Constructor.
408 
409  g is a parameter that determines how forward the scattering is.
410  g = 0 corresponds to isotropic scattering
411  g = 1 corresponds to fully forward scattering
412  **/
413  JHenyeyGreensteinScattering( double _lambda_scat, double _g ) {
414  lambda_scat = _lambda_scat ;
415  g = _g ;
416  }
417 
418  double getScatteringProbability( double ct ) {
419  //double ct = cos(dir.getTheta()) ;
420  double den = 1 + g*g - 2*g*ct ;
421  return (1-g*g)/(4*M_PI*den*sqrt(den)) ;
422  }
423 
425  double phi = gRandom->Uniform(0,2*M_PI) ;
426  double ct ;
427  if( g>0 ) {
428  double xi = gRandom->Uniform(0,1) ;
429  ct = 1.0/(2*g) * ( 1 + g*g - (1-g*g)*(1-g*g)/( (1-g+2*g*xi)*(1-g+2*g*xi) ) ) ;
430  } else {
431  ct = gRandom->Uniform(-1,1) ;
432  }
433  return JVersor3D( JAngle3D( acos(ct),phi) ) ;
434  }
435 
436  protected :
437 
438  double g ;
439 
440  } ;
441 
442  /**
443  Implementation of the JScatteringModel interface with
444  isotropic scattering
445  **/
447 
448  public:
449 
450  /**
451  Constructor.
452  **/
453  JIsotropicScattering( double _lambda_scat ) {
454  lambda_scat = _lambda_scat ;
455  }
456 
457  double getScatteringProbability( double ct ) {
458  return 1.0/(4.0*M_PI) ;
459  }
460 
462  double phi = gRandom->Uniform(0,2*M_PI) ;
463  double ct = gRandom->Uniform(-1,1) ;
464  return JVersor3D( JAngle3D( acos(ct),phi) ) ;
465  }
466  } ;
467 
468  /**
469  Implementation of the JScatteringModel interface with
470  Rayleigh scattering
471  **/
473 
474  public:
475 
476  /**
477  Constructor.
478 
479  note that one should never set a<0
480  **/
481  JRayleighScattering( double _lambda_scat, double _a ) {
482  lambda_scat = _lambda_scat ;
483  a = _a ;
484  if( fabs(a)<0 ) {
485  cerr << "Fatal error in initialization of JRayleighScattering: a<0 !" << endl ;
486  exit(1) ;
487  }
488  }
489 
490  double getScatteringProbability( double ct ) {
491  //double ct = cos(dir.getTheta()) ;
492  return (1+a*ct*ct)/(4.0*M_PI*(1+a/3.0)) ;
493  }
494 
496  double phi = gRandom->Uniform(0,2*M_PI) ;
497  double ct ;
498  if( a>0 ) {
499  double xi = gRandom->Uniform(0,1) ;
500  double d0 = -9/a ;
501  double d1 = 27.0*(1-2*xi)*(3+a)/a ;
502  double C = pow( 0.5*(d1 + sqrt(d1*d1-4*d0*d0*d0) ), 1.0/3.0 ) ;
503  ct = -1.0/(3.0)*(C+d0/C) ;
504  } else {
505  ct = gRandom->Uniform(-1,1) ;
506  }
507  return JVersor3D( JAngle3D( acos(ct),phi) ) ;
508  }
509 
510  protected:
511 
512  double a ;
513  } ;
514 
515  /**
516  Implementation of the JScatteringModel interface with
517  that combines two scattering models into one effective model.
518  **/
520 
521  public:
522 
523  /**
524  Constructor.
525 
526  Note that this only copies only the pointers to the JScatteringModel
527  instances. So if _sm1 or _sm2 are deleted or changed, it will also
528  affect the behaviour of this instance.
529  **/
530  JCombinedScattering( JScatteringModel* _sm1, JScatteringModel* _sm2 ) : sm1(_sm1), sm2(_sm2) { }
531 
532  virtual double getScatteringLength() {
533  // calculate effective scattering length
534  double l1 = sm1->getScatteringLength() ;
535  double l2 = sm2->getScatteringLength() ;
536  return 1.0 / ( 1.0/l1 + 1.0/l2 ) ;
537  }
538 
539  double getScatteringProbability( double ct ) {
540  double val1 = sm1->getScatteringProbability(ct) ;
541  double val2 = sm2->getScatteringProbability(ct) ;
542  double l1 = sm1->getScatteringLength() ;
543  double l2 = sm2->getScatteringLength() ;
544  double lambda = getScatteringLength() ;
545  // effective scattering probability
546  return lambda * ( val1/l1 + val2/l2 ) ;
547  }
548 
550  double l1 = sm1->getScatteringLength() ;
551  double l2 = sm2->getScatteringLength() ;
552  // probability to scatter with process 1
553  double P = l2/(l1+l2) ;
554  if( gRandom->Uniform(0,1)<P ) {
555  return sm1->generateDirection() ;
556  } else {
557  return sm2->generateDirection() ;
558  }
559  }
560 
561  protected:
562 
563  double a ;
566  } ;
567 
568  /**
569  Return the probability density for a photon path with the given
570  ingredients. These are
571  - a model of the source direction distribution
572  - a scattering model
573  - a target model
574  - the absorption length (note that this is allowed to be +INFTY)
575 
576  Multiply by dV1, dV2, etc. and the target cross-section sigma
577  to get a probability.
578 
579  The first (last) vertex of the photon path is used as the source
580  (target) position.
581  **/
582  double getPhotonPathProbabilityDensity( JPhotonPath& p, JSourceModel* src, JScatteringModel* sm, JTargetModel* trg, double lambda_abs ) {
583  //const bool verbose = true ;
584 
585  // calculate inverse scattering and absorption lengths
586  const double lambda_scat_inv = 1.0/sm->getScatteringLength() ;
587  double lambda_abs_inv = 0 ;
588  // if the absorption length is not infinite
589  if( lambda_abs == lambda_abs ) lambda_abs_inv = 1.0/lambda_abs ;
590 
591  double rho = 1 ;
592  const int nvert = p.n ; // number of vertices in the path
593  const int nscat = p.n-2 ; // number of scattering vertices in the path
594 
595  //if( verbose ) cout << "$ rho = " << rho << " (before calculating path segment lengths)" << endl ;
596 
597  // calculate path segment lengths
598  // index i corresponds to the path segment connecting
599  // vertex i and vertex i+1
600  double Ltotal = 0 ;
601  vector<double> lengths(nvert-1) ;
602  for( int i=0 ; i<nvert-1 ; ++i ) {
603  lengths[i] = (p[i+1]-p[i]).getLength() ;
604  Ltotal += lengths[i] ;
605  }
606 
607  //if( verbose ) cout << "$ Total length = " << Ltotal << " m" << endl ;
608 
609  // cosine of scattering angles (angle between previous segment and new segment)
610  // index i is the cosine of the angle of segment i+1 w.r.t segment i
611  vector<double> cts(nscat) ;
612 
613  // scattering angles
614  for( int i=0 ; i<nscat ; ++i ) {
615  JVersor3D seg1(p[i+1]-p[i]) ;
616  JVersor3D seg2(p[i+2]-p[i+1]) ;
617  cts[i] = seg1.getDot(seg2) ;
618  }
619 
620  // term for absorption length
621  // (=probability to not be absorbed within Ltotal)
622  rho *= exp(-Ltotal*lambda_abs_inv ) ;
623 
624  //if( verbose ) cout << "$ rho = " << rho << " (after absorption length)" << endl ;
625 
626  // term for emission profile
627  rho *= src->getEmissionProbability( JVersor3D(p[1]-p[0]) ) ;
628 
629  //if( verbose ) cout << "$ rho = " << rho << " (after emission probability)" << endl ;
630 
631  // term for target efficiency
632  if( trg->getRadius()>0 ) {
633  // finite target (spherical), the efficiency is a measure of how good the
634  // target is at a given spot
635  rho *= trg->getEfficiency( JVersor3D(p[nscat+1]-trg->getPosition()) ) ;
636  // this factor accounts for the angle of incidence
637  double ct = JVersor3D(p[nscat]-p[nscat+1]).getDot( JVersor3D(p[nscat+1]-trg->getPosition()) ) ;
638  rho *= max(0.0,ct) ;
639 
640  //if( verbose ) cout << "$ rho = " << rho << " (after target efficiency and angle of incidence)" << endl ;
641 
642  // check whether the photon path hits the sphere
643  // we do not need to consider the last vertex, because if that intersects the sphere
644  // the angle of incidence factor above would be 0.
645  if( p.n>2 ) {
646  JPhotonPath pshort(p) ;
647  pshort.resize( pshort.size()-1 ) ;
648  --pshort.n ;
649  if( pshort.hitsSphere(trg->getPosition(),trg->getRadius()) ) {
650  /*cout << "Path hits sphere at "
651  << JPosition3D(pshort.getSphereHitPosition(trg->getPosition(),trg->getRadius())-trg->getPosition())
652  << endl ;*/
653  rho *= 0 ;
654  }
655  }
656  } else {
657  rho *= trg->getEfficiency( p[nscat+1]-p[nscat] ) ;
658  }
659 
660  //if( verbose ) cout << "$ rho = " << rho << " (after target efficiency)" << endl ;
661 
662  // terms for scattering directions
663  for( int i=0 ; i<nscat ; ++i ) {
664  rho *= sm->getScatteringProbability(cts[i]) ;
665  }
666 
667  //if( verbose ) cout << "$ rho = " << rho << " (after scattering directions)" << endl ;
668 
669  // phase space terms for scattering to each of the scattering
670  // vertices
671  // this is the probability/dV to scatter into a volume element
672  // dV at a given distance
673  // (although it is missing the direction-dependent term, which is
674  // the emission probability or scattering probability)
675  for( int i=0; i<nscat; ++i ) {
676  // reminder: lengths[i] is the length from vertex i to vertex i+1
677  rho *= lambda_scat_inv / (lengths[i]*lengths[i]) * exp(-lengths[i]*lambda_scat_inv) ;
678  }
679 
680  //if( verbose ) cout << "$ rho = " << rho << " (after scattering phase space term)" << endl ;
681 
682  // phase space term for scattering to the target
683  // this is the probability/dsigma to hit a target at a given
684  // distance, without scattering along the way
685  rho *= exp(-lengths[nscat]*lambda_scat_inv) / ( lengths[nscat]*lengths[nscat] ) ;
686 
687  //if( verbose ) cout << "$ rho = " << rho << " (after target phase space term)" << endl ;
688 
689  // check for NaN
690  if( rho!=rho ) {
691  cerr << "NaN for probability density" << endl ;
692  cerr << "rho = " << rho << endl ;
693  exit(1) ;
694  }
695 
696  //if( verbose ) cout << "$ rho = " << rho << endl ;
697 
698  return rho ;
699  }
700 
701 }
702 
703 #endif
Data structure for angles in three dimensions.
Definition: JAngle3D.hh:33
Implementation of the JSourceModel class that represents an isotropic source.
virtual double getScatteringLength()
JVersor3D generateDirection()
Return a randomly generated direction according to the emission distribution.
JVersor3D generateDirection()
Return a randomly generated direction according to the scattering probability distribution.
void setRadius(double _r)
JDirectedSource(JVersor3D _source_dir, double _alpha, bool _cutSurface=false, const JVersor3D _surface_normal=JVersor3D())
Constructor.
JCosineTarget(JVersor3D _target_dir)
Constructor.
Implementation of the JScatteringModel interface with isotropic scattering.
double getRadius() const
const JPosition3D & getPosition() const
JVersor3D generateDirection()
Return a randomly generated direction according to the scattering probability distribution.
A photon path.
Definition: JPhotonPath.hh:38
Rotation matrix.
Definition: JRotation3D.hh:111
Implementation of the JTargetModel class that represents a spherically symmetric target.
double getScatteringProbability(double ct)
Return the probability density as a function of cos(theta)
double getEfficiency(JVersor3D dir) const
Return the efficiency, which is defined as the probability that a photon with the given direction hit...
double getScatteringProbability(const double x)
Function to describe light scattering in water.
Definition: Antares.hh:210
exit
Definition: JPizza.sh:36
void setDirection(JVersor3D &_dir)
This include file serves the purpose of hiding ROOT dependencies and circumphere namespace problems w...
data_type r[M+1]
Definition: JPolint.hh:868
JPMTTarget(JAxis3D _axis, double _alpha, double _r=0.2159)
Constructor.
void setPosition(JPosition3D &_pos)
Implementation of the JScatteringModel interface with scattering according to the Henyey-Greenstein f...
double getEfficiency(JVersor3D dir) const
Return the efficiency, which is defined as the probability that a photon with the given direction hit...
Axis object.
Definition: JAxis3D.hh:38
bool hitsSphere(const JPosition3D &pos, double r)
Returns whether the photon path intersects a sphere of radius r at position pos.
Definition: JPhotonPath.hh:176
Implementation of the JScatteringModel interface with Rayleigh scattering.
static const double C
Physics constants.
double getEmissionProbability(JVersor3D dir)
Return the probability density.
const JPosition3D & getPosition() const
Properties of KM3NeT PMT and deep-sea water.
virtual double getEmissionProbability(JVersor3D dir)=0
Return the probability density.
JVersor3D getDirection() const
return the source direction
Implementation of the JSourceModel class that represents a directed source with a flat intensity dist...
double getOpeningAngle() const
return the opening angle in radians
double getScatteringProbability(double ct)
Return the probability density as a function of cos(theta)
then JCalibrateToT a
Definition: JTuneHV.sh:113
const JVersor3D & getDirection() const
Implementation of the JTargetModel class that represents a directed target with a cos(theta) acceptan...
double getOpeningAngle() const
return the PMT opening angle in radians
JVersor3D generateDirection()
Return a randomly generated direction according to the scattering probability distribution.
virtual double getEfficiency(JVersor3D dir) const =0
Return the efficiency, which is defined as the probability that a photon with the given direction hit...
void setPosition(JPosition3D &_pos)
Virtual base class for a light source.
T pow(const T &x, const double y)
Power .
Definition: JMath.hh:97
JHenyeyGreensteinScattering(double _lambda_scat, double _g)
Constructor.
const JPosition3D & getPosition() const
Get position.
Definition: JPosition3D.hh:130
double getEffectiveArea() const
virtual double getEffectiveArea()
Return the effective area, i.e.
Implementation of the JTargetModel class that represents a single PMT on a DOM.
double getScatteringProbability(double ct)
Return the probability density as a function of cos(theta)
then JCookie sh JDataQuality D $DETECTOR_ID R
Definition: JDataQuality.sh:41
Implementation of the JScatteringModel interface with that combines two scattering models into one ef...
double getPhiMax(double theta)
For a given angle theta with respect to the beam direction, and a given surface normal, we define phi as the rotation angle around the beam, where phi=0 corresponds to the direction sticking out over the surface a much as possible.
double getScatteringProbability(double ct)
Return the probability density as a function of cos(theta)
JVersor3D generateDirection()
Return a randomly generated direction according to the scattering probability distribution.
double getScatteringLength(const double lambda)
Get scattering length.
Definition: Antares.hh:148
virtual double getScatteringProbability(double ct)=0
Return the probability density as a function of cos(theta)
JCombinedScattering(JScatteringModel *_sm1, JScatteringModel *_sm2)
Constructor.
double getEfficiency(JVersor3D dir) const
Return the efficiency, which is defined as the probability that a photon with the given direction hit...
Virtual base class for a scattering model.
double getEmissionProbability(JVersor3D dir)
Return the probability density.
JRayleighScattering(double _lambda_scat, double _a)
Constructor.
Virtual base class for a light detector (&quot;photon target&quot;).
double getDot(const JVersor3D &versor) const
Get dot product.
Definition: JVersor3D.hh:156
Data structure for position in three dimensions.
Definition: JPosition3D.hh:36
virtual double getScatteringLength()
Data structure for normalised vector in three dimensions.
Definition: JVersor3D.hh:26
JIsotropicScattering(double _lambda_scat)
Constructor.
double getPhotonPathProbabilityDensity(JPhotonPath &p, JSourceModel *src, JScatteringModel *sm, JTargetModel *trg, double lambda_abs)
Return the probability density for a photon path with the given ingredients.
then fatal Wrong number of arguments fi set_variable DETECTOR $argv[1] set_variable STRING $argv[2] set_array QUANTILES set_variable FORMULA *[0] exp(-0.5 *(x-[1])*(x-[1])/([2]*[2]))" set_variable MODULE `getModule -a $DETECTOR -L "$STRING 0"` source JAcousticsToolkit.sh typeset -A TRIPODS get_tripods $WORKDIR/tripod.txt TRIPODS XMEAN
JVersor3D generateDirection()
Return a randomly generated direction according to the emission distribution.
then $DIR JPlotNPE PDG P
Definition: JPlotNPE-PDG.sh:62