Go to the documentation of this file. 1 #ifndef H_J_SCATTERING_MODEL
2 #define H_J_SCATTERING_MODEL
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 ;
Implementation of the JSourceModel class that represents an isotropic source.
double getPhiMax(double theta)
For a given angle theta with respect to the beam direction, and a given surface normal,...
double getScatteringProbability(double ct)
Return the probability density as a function of cos(theta)
double getEffectiveArea() const
double getScatteringProbability(double ct)
Return the probability density as a function of cos(theta)
Implementation of the JScatteringModel interface with isotropic scattering.
double getEmissionProbability(JVersor3D dir)
Return the probability density.
const JVersor3D & getDirection() const
Implementation of the JSourceModel class that represents a directed source with a flat intensity dist...
double getOpeningAngle() const
return the opening angle in radians
Implementation of the JTargetModel class that represents a spherically symmetric target.
JVersor3D getDirection() const
return the source direction
Implementation of the JScatteringModel interface with that combines two scattering models into one ef...
double getScatteringProbability(const double x)
Function to describe light scattering in water.
JVersor3D generateDirection()
Return a randomly generated direction according to the scattering probability distribution.
JRayleighScattering(double _lambda_scat, double _a)
Constructor.
Implementation of the JTargetModel class that represents a directed target with a cos(theta) acceptan...
Implementation of the JScatteringModel interface with Rayleigh scattering.
void setPosition(JPosition3D &_pos)
Data structure for normalised vector in three dimensions.
JHenyeyGreensteinScattering(double _lambda_scat, double _g)
Constructor.
void setDirection(JVersor3D &_dir)
JCombinedScattering(JScatteringModel *_sm1, JScatteringModel *_sm2)
Constructor.
virtual double getEmissionProbability(JVersor3D dir)=0
Return the probability density.
JDirectedSource(JVersor3D _source_dir, double _alpha, bool _cutSurface=false, const JVersor3D _surface_normal=JVersor3D())
Constructor.
double getScatteringProbability(double ct)
Return the probability density as a function of cos(theta)
virtual double getScatteringLength()
JVersor3D generateDirection()
Return a randomly generated direction according to the scattering probability distribution.
double getEfficiency(JVersor3D dir) const
Return the efficiency, which is defined as the probability that a photon with the given direction hit...
Data structure for position in three dimensions.
Implementation of the JTargetModel class that represents a single PMT on a DOM.
const JPosition3D & getPosition() const
double getDot(const JVersor3D &versor) const
Get dot product.
double getScatteringProbability(double ct)
Return the probability density as a function of cos(theta)
virtual ~JScatteringModel()
Virtual base class for a light source.
void setPosition(JPosition3D &_pos)
JVersor3D generateDirection()
Return a randomly generated direction according to the scattering probability distribution.
Virtual base class for a scattering model.
double getEmissionProbability(JVersor3D dir)
Return the probability density.
Data structure for angles in three dimensions.
virtual double getEffectiveArea()
Return the effective area, i.e.
virtual double getEfficiency(JVersor3D dir) const =0
Return the efficiency, which is defined as the probability that a photon with the given direction hit...
Virtual base class for a light detector ("photon target").
const JPosition3D & getPosition() const
Get position.
virtual double getScatteringLength()
double getOpeningAngle() const
return the PMT opening angle in radians
JPMTTarget(JAxis3D _axis, double _alpha, double _r=0.2159)
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.
double getEfficiency(JVersor3D dir) const
Return the efficiency, which is defined as the probability that a photon with the given direction hit...
JVersor3D generateDirection()
Return a randomly generated direction according to the emission distribution.
JCosineTarget(JVersor3D _target_dir)
Constructor.
Implementation of the JScatteringModel interface with scattering according to the Henyey-Greenstein f...
Auxiliary classes and methods for 3D geometrical objects and operations.
const JPosition3D & getPosition() const
double getEfficiency(JVersor3D dir) const
Return the efficiency, which is defined as the probability that a photon with the given direction hit...
virtual double getScatteringProbability(double ct)=0
Return the probability density as a function of cos(theta)
double getScatteringLength(const double lambda)
Scattering length.
bool hitsSphere(const JPosition3D &pos, double r)
Returns whether the photon path intersects a sphere of radius r at position pos.
JVersor3D generateDirection()
Return a randomly generated direction according to the scattering probability distribution.
JIsotropicScattering(double _lambda_scat)
Constructor.
JVersor3D generateDirection()
Return a randomly generated direction according to the emission distribution.
void setRadius(double _r)