1 #ifndef H_J_SCATTERING_MODEL
2 #define H_J_SCATTERING_MODEL
40 using namespace JGEOMETRY3D ;
69 virtual JVersor3D generateDirection() = 0 ;
97 virtual double getEmissionProbability(
JVersor3D dir ) = 0 ;
105 virtual JVersor3D generateDirection() = 0 ;
126 return 1.0/(4.0*M_PI) ;
131 gRandom->Sphere(dx,dy,dz,1) ;
154 source_dir(_source_dir), ctmin(cos(0.5*_alpha)), alpha(_alpha), cutSurface(_cutSurface), surface_normal(_surface_normal) {
161 delta = 0.5*M_PI - acos(surface_normal.getDot(source_dir)) ;
162 omega = 2.0*M_PI*(1.0-ctmin) ;
165 if( delta>0.0 && delta>0.5*_alpha ) {
168 }
else if( delta<=0 && fabs(delta)>=0.5*alpha ) {
172 cerr <<
"FATAL ERROR in JDirectedSource: beam points completely into the surface." << endl ;
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 ;
192 omega_absorbed *= __dct ;
193 omega -= omega_absorbed ;
197 omega -= 2.0 * M_PI * (1.0-cos(delta)) ;
210 double ct = dir.
getDot(source_dir) ;
212 if( cutSurface && dir.
getDot(surface_normal)<0 )
return 0.0 ;
220 double phi = gRandom->Uniform(0,2*M_PI) ;
221 double ct = gRandom->Uniform(ctmin,1) ;
225 if( !cutSurface )
return dir ;
226 if( dir.getDot(surface_normal)>0 )
return dir ;
243 if( theta<delta )
return M_PI ;
244 if( theta+delta > M_PI )
return 0.0 ;
247 if( theta<fabs(delta) )
return 0.0 ;
248 if( theta+fabs(delta) > M_PI )
return M_PI ;
250 double phimax = acos(-tan(delta)/tan(theta)) ;
289 virtual double getEfficiency(
JVersor3D dir )
const = 0 ;
335 ctmin = cos(0.5*_alpha) ;
342 const double ct = axis.getDirection().getDot(dir) ;
351 return 2.0*
r*
r*M_PI*(1-ctmin) ;
386 target_dir(_target_dir) {}
389 return max( dir.
getDot(target_dir), 0.0 ) ;
414 lambda_scat = _lambda_scat ;
420 double den = 1 + g*g - 2*g*ct ;
421 return (1-g*g)/(4*M_PI*den*sqrt(den)) ;
425 double phi = gRandom->Uniform(0,2*M_PI) ;
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) ) ) ;
431 ct = gRandom->Uniform(-1,1) ;
454 lambda_scat = _lambda_scat ;
458 return 1.0/(4.0*M_PI) ;
462 double phi = gRandom->Uniform(0,2*M_PI) ;
463 double ct = gRandom->Uniform(-1,1) ;
482 lambda_scat = _lambda_scat ;
485 cerr <<
"Fatal error in initialization of JRayleighScattering: a<0 !" << endl ;
492 return (1+
a*ct*ct)/(4.0*M_PI*(1+
a/3.0)) ;
496 double phi = gRandom->Uniform(0,2*M_PI) ;
499 double xi = gRandom->Uniform(0,1) ;
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) ;
505 ct = gRandom->Uniform(-1,1) ;
534 double l1 = sm1->getScatteringLength() ;
535 double l2 = sm2->getScatteringLength() ;
536 return 1.0 / ( 1.0/l1 + 1.0/l2 ) ;
540 double val1 = sm1->getScatteringProbability(ct) ;
541 double val2 = sm2->getScatteringProbability(ct) ;
542 double l1 = sm1->getScatteringLength() ;
543 double l2 = sm2->getScatteringLength() ;
546 return lambda * ( val1/l1 + val2/l2 ) ;
550 double l1 = sm1->getScatteringLength() ;
551 double l2 = sm2->getScatteringLength() ;
553 double P = l2/(l1+l2) ;
554 if( gRandom->Uniform(0,1)<
P ) {
555 return sm1->generateDirection() ;
557 return sm2->generateDirection() ;
587 double lambda_abs_inv = 0 ;
589 if( lambda_abs == lambda_abs ) lambda_abs_inv = 1.0/lambda_abs ;
592 const int nvert = p.
n ;
593 const int nscat = p.
n-2 ;
602 for(
int i=0 ; i<nvert-1 ; ++i ) {
603 lengths[i] = (p[i+1]-p[i]).getLength() ;
604 Ltotal += lengths[i] ;
614 for(
int i=0 ; i<nscat ; ++i ) {
617 cts[i] = seg1.
getDot(seg2) ;
622 rho *=
exp(-Ltotal*lambda_abs_inv ) ;
647 pshort.resize( pshort.size()-1 ) ;
663 for(
int i=0 ; i<nscat ; ++i ) {
675 for(
int i=0; i<nscat; ++i ) {
677 rho *= lambda_scat_inv / (lengths[i]*lengths[i]) *
exp(-lengths[i]*lambda_scat_inv) ;
685 rho *=
exp(-lengths[nscat]*lambda_scat_inv) / ( lengths[nscat]*lengths[nscat] ) ;
691 cerr <<
"NaN for probability density" << endl ;
692 cerr <<
"rho = " << rho << endl ;
Data structure for angles in three dimensions.
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.
const JPosition3D & getPosition() const
JVersor3D generateDirection()
Return a randomly generated direction according to the scattering probability distribution.
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.
void setDirection(JVersor3D &_dir)
This include file serves the purpose of hiding ROOT dependencies and circumphere namespace problems w...
JPMTTarget(JAxis3D _axis, double _alpha, double _r=0.2159)
Constructor.
void setPosition(JPosition3D &_pos)
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 JAcoustics.sh -- typeset -A TRIPODS get_tripods $WORKDIR/tripod.txt TRIPODS XMEAN
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...
bool hitsSphere(const JPosition3D &pos, double r)
Returns whether the photon path intersects a sphere of radius r at position pos.
virtual ~JScatteringModel()
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)
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 .
JHenyeyGreensteinScattering(double _lambda_scat, double _g)
Constructor.
const JPosition3D & getPosition() const
Get position.
then usage $script[distance] fi case set_variable R
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)
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.
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 ("photon target").
double getDot(const JVersor3D &versor) const
Get dot product.
Data structure for position in three dimensions.
virtual double getScatteringLength()
Data structure for normalised vector in three dimensions.
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.
JVersor3D generateDirection()
Return a randomly generated direction according to the emission distribution.