1 #ifndef H_MARKOV_INTEGRATOR
2 #define H_MARKOV_INTEGRATOR
25 namespace JPP {
using namespace JMARKOV; }
99 for(
int nv=0 ; nv<nscat+2 ; ++nv ) {
133 for(
int i=0 ; i<
N ; ++i ) {
138 contributions[i] = rho * winv ;
140 return contributions ;
146 for(
int i=0 ; i<
N ; ++i ) {
147 const double r = 10 ;
150 double rho = 1.0/(4.0/3.0*M_PI*r*r*
r) ;
151 if( p[1].getLength()>
r ) rho = 0 ;
152 contributions[i] = rho * winv ;
154 return contributions ;
160 for(
int i=0 ; i<
N ; ++i ) {
205 cerr <<
"FATAL ERROR in JMarkov3DensembleIntegrator constructor. Ensemble size = 0." << endl ;
211 if( (*it)[0].getDistance(
ensemble.front()[0]) > 0 ) {
212 cerr <<
"Fatal error in JMarkovEnsembleIntegrator, first vertex positions are not all the same." << endl ;
235 if( nv == nscat+1 ) {
259 return (nscat*(nscat-1))/2 + nv - 1 ;
278 sprintf( hname,
"htarget_nscat%i", nscat ) ;
279 TH2D* ht =
new TH2D(hname,
";cos(#theta);#phi",100,-1,1,100,-M_PI,M_PI) ;
283 ht->Fill( dir.getDZ(), dir.getPhi() ) ;
297 if( *it != NULL )
delete *it ;
326 if( *it != NULL )
delete *it ;
340 if( (*it)->h != NULL ) (*it)->h->Write() ;
350 if( (
int)
gens.size() < ilast )
gens.resize( ilast, NULL ) ;
351 for(
int nv=1 ; nv<=nscat ; ++nv ) {
352 int index = ifirst + nv - 1 ;
354 double xmin = 1.0/0.0 ;
355 double xmax = -1.0/0.0 ;
356 double ymin = 1.0/0.0 ;
357 double ymax = -1.0/0.0 ;
358 double zmin = 1.0/0.0 ;
359 double zmax = -1.0/0.0 ;
365 if( vtx.
getY()<ymin ) ymin = vtx.
getY() ;
366 if( vtx.
getY()>ymax ) ymax = vtx.
getY() ;
367 if( vtx.
getZ()<zmin ) zmin = vtx.
getZ() ;
368 if( vtx.
getZ()>zmax ) zmax = vtx.
getZ() ;
372 sprintf( hname,
"hx_nscat%i_vtx%i", nscat, nv ) ;
373 TH1F hx(hname,
"Vertex position;X",
nbinsx,xmin,xmax) ;
374 sprintf( hname,
"hy_nscat%i_vtx%i", nscat, nv ) ;
375 TH1F hy(hname,
"Vertex position;Y",
nbinsy,ymin,ymax) ;
376 sprintf( hname,
"hz_nscat%i_vtx%i", nscat, nv ) ;
377 TH1F hz(hname,
"Vertex position;Z",
nbinsz,zmin,zmax) ;
382 hx.Fill( vtx.
getX() ) ;
383 hy.Fill( vtx.
getY() ) ;
384 hz.Fill( vtx.
getZ() ) ;
387 for( Int_t bin=1 ; bin<=hx.GetNbinsX() ; ++bin )
388 hx.AddBinContent(bin,1) ;
389 for( Int_t bin=1 ; bin<=hy.GetNbinsX() ; ++bin )
390 hy.AddBinContent(bin,1) ;
391 for( Int_t bin=1 ; bin<=hz.GetNbinsX() ; ++bin )
392 hz.AddBinContent(bin,1) ;
404 winv = 1.0/
gens[index]->getWeight(pos) ;
428 JPosition3D generatePosition( int nscat, int nv, double& winv ) {
429 return JPosition3D(0,0,0) ;
433 JPhotonPath generatePath( int nscat, double& winv ) {
436 // distance from source to target
437 double d = (source_pos-target_pos).getLength() ;
439 double g = d / (l+d) ;
440 JHenyeyGreensteinScattering sm(100,g) ;
442 double lambda = (d+l)/(nscat+1) ;
444 JPhotonPath p(nscat) ;
445 p.front() = source_pos ;
446 p.back() = target_pos ;
448 for( int nv=1 ; nv<=nscat ; ++nv ) {
449 // generate travel length
450 double r = gRandom->Exp(lambda) ;
451 // generate direction from a HG distribution, relative to THE TARGET.
452 JVersor3D dir( p.back()-p[nv-1] ) ;
453 JRotation3D rot(dir) ;
454 dir = sm.generateDirection() ;
455 double ct = dir.getDZ() ;
456 dir = dir.rotate_back(rot) ;
458 p[nv] = p[nv-1] + r*JVector3D(dir) ;
460 winvtot /= ( r * r * sm.getScatteringProbability(ct) * 1.0/lambda*exp(-r/lambda) ) ;
477 JMarkovEnsembleIntegrator3D(
const vector<JPhotonPath>& _ensemble,
JTargetModel* _trg,
int _nbinsx,
int _nbinsy,
int _nbinsz,
JPosition3D _posmin,
JPosition3D _posmax ) :
JMarkovEnsembleIntegrator(_ensemble,_trg,_nbinsx,_nbinsy,_nbinsz),
posmin(_posmin),
posmax(_posmax) {}
487 if( *it != NULL )
delete *it ;
496 htmp = (TH2D*)(*it)->h->Project3D(
"yx") ;
497 htmp->SetOption(
"colz") ;
498 htmp->SetEntries(10) ;
500 htmp = (TH2D*)(*it)->h->Project3D(
"yz") ;
501 htmp->SetOption(
"colz") ;
502 htmp->SetEntries(10) ;
504 htmp = (TH2D*)(*it)->h->Project3D(
"zx") ;
505 htmp->SetOption(
"colz") ;
506 htmp->SetEntries(10) ;
517 winv = 1.0/
gens[index]->getWeight(pos) ;
526 if( (
int)
gens.size() < ilast )
gens.resize( ilast, NULL ) ;
528 for(
int nv=1 ; nv<=nscat ; ++nv ) {
529 int index = ifirst + nv - 1 ;
531 double xmin = 1.0/0.0 ;
532 double xmax = -1.0/0.0 ;
533 double ymin = 1.0/0.0 ;
534 double ymax = -1.0/0.0 ;
535 double zmin = 1.0/0.0 ;
536 double zmax = -1.0/0.0 ;
545 sprintf( hname,
"h3_nscat%i_vtx%i", nscat, nv ) ;
546 TH3F h(hname,
"Vertex position;X;Y;Z",
557 double excuseval = 0.01*h.Integral() / ( h.GetNbinsX() * h.GetNbinsY() * h.GetNbinsZ() ) ;
558 for( Int_t xbin=1 ; xbin<=h.GetNbinsX() ; ++xbin ) {
559 for( Int_t ybin=1 ; ybin<=h.GetNbinsY() ; ++ybin ) {
560 for( Int_t zbin=1 ; zbin<=h.GetNbinsZ() ; ++zbin ) {
561 Int_t bin = h.GetBin(xbin,ybin,zbin) ;
562 h.AddBinContent(bin,excuseval) ;
597 JSourceTargetIntegrator(
double Rsing,
JPosition3D src_pos,
JPosition3D trg_pos,
double fraction_exp,
double Lexp,
double Rexp ) :
gen_exp(0,Rexp,Lexp),
gen_s1(Rsing,src_pos),
gen_s2(Rsing,trg_pos) {
600 double fraction_sing = 0.5*(1-fraction_exp) ;
652 p.front() = source_pos ;
653 p.back() = target_pos ;
656 for(
int nv=1 ; nv<=nscat ; ++nv ) {
677 if( gRandom->Integer(2) == 0 ) {
682 return src_pos + pos ;
Implementation of the JGenerator interface.
Implementation of the JGenerator interface.
virtual JPosition3D generateScatteringVertexPosition(int nscat, int nv, double &winv)=0
void initgenerators(int nscat)
initialize the position generators for a given number of scatterings (i.e. create histograms and fill...
JPosition3D getFirstScatteringVertex(int nscat, JPosition3D src_pos, JPosition3D trg_pos, double &winv)
recursive function to randomly generate the first scattering vertex for isotropic scattering from a s...
Abstract base class for implementations of the JMarkovIntegrator interface, where the sample distribu...
int getIndex(int nv, int nscat)
return the internal index for a given number of scatterings and vertex number
In this implementation of the JMarkovIntegrator interface, the sample distribution is built up out of...
JMarkovEnsembleIntegrator1D(const vector< JPhotonPath > &_ensemble, JTargetModel *_trg, int _nbinsx, int _nbinsy, int _nbinsz)
vector< JSphereGenerator * > target_gens
Implementation of JMarkovEnsembleIntegrator interface with 1D histograms.
double getWeight(JPosition3D pos)
return the weight (=probability density dP/dV) for the given position.
virtual JPosition3D generatePosition(int nscat, int nv, double &winv)
Generate a random position for vertex nv.
JSingularityGenerator gen_s1
JExperimentalIntegrator(double _lambda)
constructor.
Abstract interface for the generation of points in 3D space.
~JSourceTargetIntegrator()
JMarkovEnsembleIntegrator(const vector< JPhotonPath > &_ensemble, JTargetModel *_trg, int _nbinsx, int _nbinsy, int _nbinsz)
~JMarkovEnsembleIntegrator3D()
void init_scattering_vertex_generators(int nscat)
This implementation of the JMarkovIntegrator interface generates 'natural' looking paths that might s...
JSingularityGenerator gen_s2
virtual void writeHistograms()=0
write the histograms filled from the ensemble
const JPosition3D & getPosition() const
void free_scattering_vertex_gens()
double getWeight(JPosition3D pos)
return the weight (=probability density dP/dV) for the given position.
JPhotonPath generatePath(int nscat, double &winv, JPosition3D source_pos, JPosition3D target_pos)
Implementation of the JGenerator interface.
virtual JPosition3D generatePosition(int nscat, int nv, double &winv)=0
Generate a random position for vertex nv.
JPosition3D generateScatteringVertexPosition(int nscat, int nv, double &winv)
Implementation of the JGenerator interface.
~JMarkovEnsembleIntegrator1D()
void init_scattering_vertex_generators(int nscat)
Virtual base class for a light source.
Implementation of the JGenerator interface.
Implementation of the JGenerator interface.
double getY() const
Get y position.
JMarkovIntegrator()
standard constructor
virtual void init_scattering_vertex_generators(int nscat)=0
JPosition3D generatePosition(int nscat, int nv, double &winv)
generate a position for vertex nv for the given number of scatterings, winv will be set to 1/weight ...
JPosition3D generatePosition(int nscat, int nv, double &winv)
dummy implementation
void init_target_generator(int nscat)
initialize the generator for a target vertex for a given number of scatterings
vector< JPhotonPath > get_diagnostic_ensemble(int N, int nscat, JSourceModel *src, JScatteringModel *sm, JTargetModel *trg, double lambda_abs)
Return photon paths generated with the generatePath method.
JPosition3D getPosition()
Return a randomly generated position.
vector< double > dummy_integrate(int N, int nscat, JSourceModel *src, JScatteringModel *sm, JTargetModel *trg, double lambda_abs)
Integrate a test function with N samples.
then usage $script< input file >[option[primary[working directory]]] nWhere option can be N
virtual JPhotonPath generatePath(int nscat, double &winv)
Generate a random photon path with a given number of scatterings.
JMarkovEnsembleIntegrator3D(const vector< JPhotonPath > &_ensemble, JTargetModel *_trg, int _nbinsx, int _nbinsy, int _nbinsz, JPosition3D _posmin, JPosition3D _posmax)
void writeHistograms()
write the histograms filled from the ensemble
JSourceTargetIntegrator(double Rsing, JPosition3D src_pos, JPosition3D trg_pos, double fraction_exp, double Lexp, double Rexp)
constructor.
vector< JHistGenerator * > gens
JGenerator * gen_exp_shift
virtual ~JMarkovEnsembleIntegrator()
double getX() const
Get x position.
Virtual base class for a scattering model.
Implementation of the JGenerator interface.
Virtual base class for a light detector ("photon target").
Data structure for position in three dimensions.
vector< JPhotonPath > ensemble
Data structure for normalised vector in three dimensions.
vector< J3DhistGenerator * > gens
JPosition3D source_position
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.
~JExperimentalIntegrator()
JPosition3D generateScatteringVertexPosition(int nscat, int nv, double &winv)
In this implementation of the JMarkovIntegrator interface, the sample distribution is built up out of...
double getZ() const
Get z position.
void free_scattering_vertex_gens()
Abstract base class for calculating the total probability (/m^2 target cross-section) for a photon fr...
JExponentialGenerator gen_exp
vector< double > integrate(int N, int nscat, JSourceModel *src, JScatteringModel *sm, JTargetModel *trg, double lambda_abs)
Integrate with N samples.
Implementation of the JGenerator interface.
void writeHistograms()
write the histograms filled from the ensemble
JPosition3D getPosition()
Return a randomly generated position.