This implementation of the JMarkovIntegrator interface generates 'natural' looking paths that might sample the phase space well in some cases. More...

#include <JMarkovIntegrator.hh>

Inheritance diagram for JMARKOV::JMarkovEnsembleIntegrator3D:

Public Member Functions
	JMarkovEnsembleIntegrator3D (const vector< JPhotonPath > &_ensemble, JTargetModel *_trg, int _nbinsx, int _nbinsy, int _nbinsz, JPosition3D _posmin, JPosition3D _posmax)

	~JMarkovEnsembleIntegrator3D ()

void	free_scattering_vertex_gens ()

void	writeHistograms ()
	write the histograms filled from the ensemble More...

JPosition3D	generateScatteringVertexPosition (int nscat, int nv, double &winv)

vector< double >	integrate (int N, int nscat, JSourceModel src, JScatteringModel sm, JTargetModel *trg, double lambda_abs)
	Integrate with N samples. More...

vector< double >	dummy_integrate (int N, int nscat, JSourceModel src, JScatteringModel sm, JTargetModel *trg, double lambda_abs)
	Integrate a test function with N samples. More...

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. More...

Protected Member Functions
void	init_scattering_vertex_generators (int nscat)

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 More...

int	getIndex (int nv, int nscat)
	return the internal index for a given number of scatterings and vertex number More...

void	initgenerators (int nscat)
	initialize the position generators for a given number of scatterings (i.e. create histograms and fill them from the ensemble) More...

void	init_target_generator (int nscat)
	initialize the generator for a target vertex for a given number of scatterings More...

void	free_target_gens ()

virtual JPhotonPath	generatePath (int nscat, double &winv)
	Generate a random photon path with a given number of scatterings. More...

Protected Attributes
vector< J3DhistGenerator * >	gens

JPosition3D	posmin

JPosition3D	posmax

vector< JSphereGenerator * >	target_gens

vector< JPhotonPath >	ensemble

int	nbinsx

int	nbinsy

int	nbinsz

JTargetModel *	trg

JPosition3D	source_position

Detailed Description

This implementation of the JMarkovIntegrator interface generates 'natural' looking paths that might sample the phase space well in some cases.

Similar to JMarkovEnsembleIntegrator, but using a 3D histogram for each vertex instead of 3 1D histograms.

Definition at line 473 of file JMarkovIntegrator.hh.

Constructor & Destructor Documentation

◆ JMarkovEnsembleIntegrator3D()

JMARKOV::JMarkovEnsembleIntegrator3D::JMarkovEnsembleIntegrator3D	(	const vector< JPhotonPath > &	_ensemble,
		JTargetModel *	_trg,
		int	_nbinsx,
		int	_nbinsy,
		int	_nbinsz,
		JPosition3D	_posmin,
		JPosition3D	_posmax
	)

inline

Definition at line 477 of file JMarkovIntegrator.hh.

477 : JMarkovEnsembleIntegrator(_ensemble,_trg,_nbinsx,_nbinsy,_nbinsz), posmin(_posmin), posmax(_posmax) {}

◆ ~JMarkovEnsembleIntegrator3D()

JMARKOV::JMarkovEnsembleIntegrator3D::~JMarkovEnsembleIntegrator3D ( )

inline

Definition at line 479 of file JMarkovIntegrator.hh.

                                    {
       free_target_gens() ;
       free_scattering_vertex_gens() ;
     }

Member Function Documentation

◆ free_scattering_vertex_gens()

void JMARKOV::JMarkovEnsembleIntegrator3D::free_scattering_vertex_gens ( )

inline

Definition at line 484 of file JMarkovIntegrator.hh.

                                        {
       // free memory
       for( vector<J3DhistGenerator*>::iterator it=gens.begin() ; it!=gens.end() ; ++it ) {
         if( *it != NULL ) delete *it ;
       }
     }

◆ writeHistograms()

void JMARKOV::JMarkovEnsembleIntegrator3D::writeHistograms ( )

inlinevirtual

write the histograms filled from the ensemble

Implements JMARKOV::JMarkovEnsembleIntegrator.

Definition at line 491 of file JMarkovIntegrator.hh.

                            {
       for( vector<J3DhistGenerator*>::iterator it=gens.begin() ; it!=gens.end() ; ++it ) {
         if( *it != NULL ) {
           (*it)->h->Write() ;
           TH2D* htmp ;
           htmp = (TH2D*)(*it)->h->Project3D("yx") ;
           htmp->SetOption("colz") ;
           htmp->SetEntries(10) ;
           htmp->Write() ;
           htmp = (TH2D*)(*it)->h->Project3D("yz") ;
           htmp->SetOption("colz") ;
           htmp->SetEntries(10) ;
           htmp->Write() ;
           htmp = (TH2D*)(*it)->h->Project3D("zx") ;
           htmp->SetOption("colz") ;
           htmp->SetEntries(10) ;
           htmp->Write() ;
         }
       }
     }

◆ generateScatteringVertexPosition()

JPosition3D JMARKOV::JMarkovEnsembleIntegrator3D::generateScatteringVertexPosition	(	int	nscat,
		int	nv,
		double &	winv
	)

inlinevirtual

Implements JMARKOV::JMarkovEnsembleIntegrator.

Definition at line 512 of file JMarkovIntegrator.hh.

                                                                                   {
       int index = getIndex(nv,nscat) ;
       if( index+1>(int)gens.size() || gens[index] == NULL ) initgenerators(nscat) ;
  
       JPosition3D pos = gens[index]->getPosition() ;
       winv = 1.0/gens[index]->getWeight(pos) ;
       return pos ;
     }

◆ init_scattering_vertex_generators()

void JMARKOV::JMarkovEnsembleIntegrator3D::init_scattering_vertex_generators ( int nscat )

inlineprotectedvirtual

Implements JMARKOV::JMarkovEnsembleIntegrator.

Definition at line 523 of file JMarkovIntegrator.hh.

                                                         {
       int ifirst = getIndex(1,nscat) ;
       int ilast  = getIndex(1,nscat+1) ;
       if( (int)gens.size() < ilast ) gens.resize( ilast, NULL ) ;
       
       for( int nv=1 ; nv<=nscat ; ++nv ) {
         int index = ifirst + nv - 1 ;
         // determine minimal and maximal value
         double xmin =  1.0/0.0 ;
         double xmax = -1.0/0.0 ;
         double ymin =  1.0/0.0 ;
         double ymax = -1.0/0.0 ;
         double zmin =  1.0/0.0 ;
         double zmax = -1.0/0.0 ;
         xmin = posmin.getX() ;
         ymin = posmin.getY() ;
         zmin = posmin.getZ() ;
         xmax = posmax.getX() ;
         ymax = posmax.getY() ;
         zmax = posmax.getZ() ;
         // allocate histogram
         char hname[200] ;
         sprintf( hname, "h3_nscat%i_vtx%i", nscat, nv ) ;
         TH3F h(hname,"Vertex position;X;Y;Z",
                nbinsx,xmin,xmax,
                nbinsy,ymin,ymax,
                nbinsz,zmin,zmax ) ;
         // fill histogram from ensemble
         for( vector<JPhotonPath>::iterator it=ensemble.begin() ; it!=ensemble.end() ; ++it ) {
           if( it->n-2 != nscat ) continue ;
           JPosition3D vtx( (*it)[nv] ) ;
           h.Fill( vtx.getX(), vtx.getY(), vtx.getZ() ) ;
         }
         // extra: make sure every field is filled at least once
         double excuseval = 0.01*h.Integral() / ( h.GetNbinsX() * h.GetNbinsY() * h.GetNbinsZ() ) ;
         for( Int_t xbin=1 ; xbin<=h.GetNbinsX() ; ++xbin ) { 
           for( Int_t ybin=1 ; ybin<=h.GetNbinsY() ; ++ybin ) { 
             for( Int_t zbin=1 ; zbin<=h.GetNbinsZ() ; ++zbin ) { 
               Int_t bin = h.GetBin(xbin,ybin,zbin) ;
               h.AddBinContent(bin,excuseval) ;
             }
           }
         }
         gens[index] = new J3DhistGenerator( &h ) ;
       }
     }

◆ generatePosition()

JPosition3D JMARKOV::JMarkovEnsembleIntegrator::generatePosition	(	int	nscat,
		int	nv,
		double &	winv
	)

inlineprotectedvirtualinherited

generate a position for vertex nv for the given number of scatterings, winv will be set to 1/weight

Implements JMARKOV::JMarkovIntegrator.

Definition at line 228 of file JMarkovIntegrator.hh.

                                                                     {
       // source vertex
       if( nv==0 ) {
         winv = 1.0 ;
         return source_position ;
       }
       // target vertex
       if( nv == nscat+1 ) {
         if( nscat+1>(int)target_gens.size() || target_gens[nscat] == NULL ) { 
           initgenerators(nscat) ; 
         }
         JPosition3D pos = target_gens[nscat]->getPosition() ;
         winv = 1.0/target_gens[nscat]->getWeight(pos) ;
         return pos ;
       } else {
         // scattering vertex
         return generateScatteringVertexPosition(nscat,nv,winv) ;
       }
     }

◆ getIndex()

int JMARKOV::JMarkovEnsembleIntegrator::getIndex	(	int	nv,
		int	nscat
	)

inlineprotectedinherited

return the internal index for a given number of scatterings and vertex number

Definition at line 251 of file JMarkovIntegrator.hh.

                                       {
       // indices are distributed as follows
       // [0]: nscat=1, nv=1
       // [1]: nscat=2, nv=1
       // [2]: nscat=2, nv=2
       // [3]: nscat=3, nv=1
       // [4]: nscat=3, nv=2
       // etc.
       return (nscat*(nscat-1))/2 + nv - 1 ;
     }

◆ initgenerators()

void JMARKOV::JMarkovEnsembleIntegrator::initgenerators ( int nscat )

inlineprotectedinherited

initialize the position generators for a given number of scatterings (i.e. create histograms and fill them from the ensemble)

Definition at line 263 of file JMarkovIntegrator.hh.

                                      {
       // target vertex
       init_target_generator(nscat) ;
  
       // scattering vertices
       init_scattering_vertex_generators(nscat) ;
     }

◆ init_target_generator()

void JMARKOV::JMarkovEnsembleIntegrator::init_target_generator ( int nscat )

inlineprotectedinherited

initialize the generator for a target vertex for a given number of scatterings

Definition at line 272 of file JMarkovIntegrator.hh.

                                             {
       if( (int)target_gens.size() <= nscat ) target_gens.resize( nscat+1 ) ;
  
       if( trg->getRadius()>0 ) {
         // for a finite target, we fill a "sphere generator" from the ensemble
         char hname[200] ;
         sprintf( hname, "htarget_nscat%i", nscat ) ;
         TH2D* ht = new TH2D(hname,";cos(#theta);#phi",100,-1,1,100,-M_PI,M_PI) ;
         for( vector<JPhotonPath>::iterator it=ensemble.begin() ; it!=ensemble.end() ; ++it ) {
           if( it->n-2 != nscat ) continue ;
           JVersor3D dir( it->back() - trg->getPosition() ) ;
           ht->Fill( dir.getDZ(), dir.getPhi() ) ;
         }
         target_gens[nscat] = new JSphereGenerator( trg->getPosition(), trg->getRadius(), ht ) ;
         delete ht ;
       } else {
         // for an infinitesimal target, we create a "generator" that only generates the single point
         target_gens[nscat] = new JSphereGenerator( trg->getPosition() ) ;
       }
     }

◆ free_target_gens()

void JMARKOV::JMarkovEnsembleIntegrator::free_target_gens ( )

inlineprotectedinherited

Definition at line 295 of file JMarkovIntegrator.hh.

                             {
       for( vector<JSphereGenerator*>::iterator it=target_gens.begin() ; it!=target_gens.end() ; ++it ) {
         if( *it != NULL ) delete *it ;
       }
     }

◆ integrate()

vector< double > JMARKOV::JMarkovIntegrator::integrate	(	int	N,
		int	nscat,
		JSourceModel *	src,
		JScatteringModel *	sm,
		JTargetModel *	trg,
		double	lambda_abs
	)

inherited

Integrate with N samples.

Returns a vector with the contribution to the integral of each sample. The mean of those values is the estimate of the result of the integral, while the distribution itself can be used to estimate the stability of the result. In this distribution, you want to avoid

long tails (because they make the result unstable)
small contributions (because it means that parts of the parameter space are being oversampled, so it is less efficient) This can be achieved by tuning the sample distribution to the problem at hand.

Definition at line 130 of file JMarkovIntegrator.hh.

                                                                                                                                                {
     vector<double> contributions(N,-1) ;
  
     for( int i=0 ; i<N ; ++i ) {
       double winv ;
       JPhotonPath p = generatePath(nscat,winv) ;
  
       double rho = getPhotonPathProbabilityDensity(p,src,sm,trg,lambda_abs) ;
       contributions[i] = rho * winv ; 
     }
     return contributions ;
   }

◆ dummy_integrate()

vector< double > JMARKOV::JMarkovIntegrator::dummy_integrate	(	int	N,
		int	nscat,
		JSourceModel *	src,
		JScatteringModel *	sm,
		JTargetModel *	trg,
		double	lambda_abs
	)

inherited

Integrate a test function with N samples.

This can be used as a sanity check for derived classes of JMarkovIntegrator.

The integral should yield 1 when the complete relevant part of the volume is taken into account. If it does not, it may be a sign that the implementation is not correct.

Returns a vector with the contribution to the integral of each sample.

Definition at line 143 of file JMarkovIntegrator.hh.

                                                                                                                                                      {
     vector<double> contributions(N,-1) ;
  
     for( int i=0 ; i<N ; ++i ) {
       const double r = 10 ;
       double winv ;
       JPhotonPath p = generatePath(nscat,winv) ;
       double rho = 1.0/(4.0/3.0*M_PI*r*r*r) ;
       if( p[1].getLength()>r ) rho = 0 ;
       contributions[i] = rho * winv ; 
     }
     return contributions ;
   }

◆ get_diagnostic_ensemble()

vector< JPhotonPath > JMARKOV::JMarkovIntegrator::get_diagnostic_ensemble	(	int	N,
		int	nscat,
		JSourceModel *	src,
		JScatteringModel *	sm,
		JTargetModel *	trg,
		double	lambda_abs
	)

inherited

Return photon paths generated with the generatePath method.

This can be used to identify the parts of parameter space that are over- or undersampled in a given problem so that the integrator may be optimized to handle those better.

Definition at line 157 of file JMarkovIntegrator.hh.

                                                                                                                                                                   {
     vector<JPhotonPath> paths ;
  
     for( int i=0 ; i<N ; ++i ) {
       double winv ;
       JPhotonPath p = generatePath(nscat,winv) ;
       paths.push_back(p) ;
     }
     return paths ;
   }

◆ generatePath()

virtual JPhotonPath JMARKOV::JMarkovIntegrator::generatePath	(	int	nscat,
		double &	winv
	)

inlineprotectedvirtualinherited

Generate a random photon path with a given number of scatterings.

winv must be set to the inverted probability density to generate this particular path.

Definition at line 94 of file JMarkovIntegrator.hh.

                                                                 {
       // this default implementation assumes that the vertex positions
       // are completely uncorrelated
       JPhotonPath p(nscat) ;
       double _winv = 1 ;
       for( int nv=0 ; nv<nscat+2 ; ++nv ) {
         double part_winv ;
         p[nv] = generatePosition( nscat, nv, part_winv ) ;
         _winv *= part_winv ;
       }
       winv = _winv ;
       return p ;
     }

Member Data Documentation

◆ gens

vector<J3DhistGenerator*> JMARKOV::JMarkovEnsembleIntegrator3D::gens

protected

Definition at line 570 of file JMarkovIntegrator.hh.

◆ posmin

JPosition3D JMARKOV::JMarkovEnsembleIntegrator3D::posmin

protected

Definition at line 571 of file JMarkovIntegrator.hh.

◆ posmax

JPosition3D JMARKOV::JMarkovEnsembleIntegrator3D::posmax

protected

Definition at line 572 of file JMarkovIntegrator.hh.

◆ target_gens

vector<JSphereGenerator*> JMARKOV::JMarkovEnsembleIntegrator::target_gens

protectedinherited

Definition at line 301 of file JMarkovIntegrator.hh.

◆ ensemble

vector<JPhotonPath> JMARKOV::JMarkovEnsembleIntegrator::ensemble

protectedinherited

Definition at line 302 of file JMarkovIntegrator.hh.

◆ nbinsx

int JMARKOV::JMarkovEnsembleIntegrator::nbinsx

protectedinherited

Definition at line 303 of file JMarkovIntegrator.hh.

◆ nbinsy

int JMARKOV::JMarkovEnsembleIntegrator::nbinsy

protectedinherited

Definition at line 304 of file JMarkovIntegrator.hh.

◆ nbinsz

int JMARKOV::JMarkovEnsembleIntegrator::nbinsz

protectedinherited

Definition at line 305 of file JMarkovIntegrator.hh.

◆ trg

JTargetModel* JMARKOV::JMarkovEnsembleIntegrator::trg

protectedinherited

Definition at line 306 of file JMarkovIntegrator.hh.

◆ source_position

JPosition3D JMARKOV::JMarkovEnsembleIntegrator::source_position

protectedinherited

Definition at line 307 of file JMarkovIntegrator.hh.

The documentation for this class was generated from the following file:

JMarkovIntegrator.hh

Public Member Functions

Protected Member Functions

Protected Attributes

Detailed Description

Constructor & Destructor Documentation

◆ JMarkovEnsembleIntegrator3D()

◆ ~JMarkovEnsembleIntegrator3D()

Member Function Documentation

◆ free_scattering_vertex_gens()

◆ writeHistograms()

◆ generateScatteringVertexPosition()

◆ init_scattering_vertex_generators()

◆ generatePosition()

◆ getIndex()

◆ initgenerators()

◆ init_target_generator()

◆ free_target_gens()

◆ integrate()

◆ dummy_integrate()

◆ get_diagnostic_ensemble()

◆ generatePath()

Member Data Documentation

◆ gens

◆ posmin

◆ posmax

◆ target_gens

◆ ensemble

◆ nbinsx

◆ nbinsy

◆ nbinsz

◆ trg

◆ source_position