JMARKOV::JMarkovEnsembleIntegrator1D Class Reference

Implementation of JMarkovEnsembleIntegrator interface with 1D histograms. More...

#include <JMarkovIntegrator.hh>

Inheritance diagram for JMARKOV::JMarkovEnsembleIntegrator1D:

Public Member Functions
	JMarkovEnsembleIntegrator1D (const vector< JPhotonPath > &_ensemble, JTargetModel *_trg, int _nbinsx, int _nbinsy, int _nbinsz)
	~JMarkovEnsembleIntegrator1D ()
void	free_scattering_vertex_gens ()
void	writeHistograms ()
	write the histograms filled from the ensemble More...
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	generateScatteringVertexPosition (int nscat, int nv, double &winv)
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< JHistGenerator * >	gens
vector< JSphereGenerator * >	target_gens
vector< JPhotonPath >	ensemble
int	nbinsx
int	nbinsy
int	nbinsz
JTargetModel *	trg
JPosition3D	source_position

Detailed Description

Implementation of JMarkovEnsembleIntegrator interface with 1D histograms.

Definition at line 313 of file JMarkovIntegrator.hh.

Constructor & Destructor Documentation

JMARKOV::JMarkovEnsembleIntegrator1D::JMarkovEnsembleIntegrator1D ( const vector< JPhotonPath > &  _ensemble, JTargetModel *  _trg, int  _nbinsx, int  _nbinsy, int  _nbinsz  )

inline

Definition at line 317 of file JMarkovIntegrator.hh.

: JMarkovEnsembleIntegrator(_ensemble,_trg,_nbinsx,_nbinsy,_nbinsz) {}

JMARKOV::JMarkovEnsembleIntegrator1D::~JMarkovEnsembleIntegrator1D ( )

inline

Definition at line 319 of file JMarkovIntegrator.hh.

                                   {
      free_target_gens() ;
      free_scattering_vertex_gens() ;
    }

Member Function Documentation

void JMARKOV::JMarkovEnsembleIntegrator1D::free_scattering_vertex_gens ( )

inline

Definition at line 324 of file JMarkovIntegrator.hh.

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

void JMARKOV::JMarkovEnsembleIntegrator1D::writeHistograms ( )

inlinevirtual

write the histograms filled from the ensemble

Implements JMARKOV::JMarkovEnsembleIntegrator.

Definition at line 330 of file JMarkovIntegrator.hh.

                           {
      for( vector<JHistGenerator*>::iterator it=gens.begin() ; it!=gens.end() ; ++it ) {
        if( *it != NULL ) {
          (*it)->hx->Write() ;
          (*it)->hy->Write() ;
          (*it)->hz->Write() ;
        }
      }
      for( vector<JSphereGenerator*>::iterator it=target_gens.begin(); it!=target_gens.end(); ++it ) {
        if( *it != NULL ) {
          if( (*it)->h != NULL ) (*it)->h->Write() ;
        }
      }
    }

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

inlineprotectedvirtual

Implements JMARKOV::JMarkovEnsembleIntegrator.

Definition at line 347 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 coordinate values for this vertex
        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 ;
        for( vector<JPhotonPath>::iterator it=ensemble.begin() ; it!=ensemble.end() ; ++it ) {
          if( it->n-2 != nscat ) continue ;
          JPosition3D vtx( (*it)[nv] ) ;
          if( vtx.getX()<xmin ) xmin = vtx.getX() ;
          if( vtx.getX()>xmax ) xmax = vtx.getX() ;
          if( vtx.getY()<ymin ) ymin = vtx.getY() ;
          if( vtx.getY()>ymax ) ymax = vtx.getY() ;
          if( vtx.getZ()<zmin ) zmin = vtx.getZ() ;
          if( vtx.getZ()>zmax ) zmax = vtx.getZ() ;
        }
        // allocate histograms
        char hname[200] ;
        sprintf( hname, "hx_nscat%i_vtx%i", nscat, nv ) ;
        TH1F hx(hname,"Vertex position;X",nbinsx,xmin,xmax) ;
        sprintf( hname, "hy_nscat%i_vtx%i", nscat, nv ) ;
        TH1F hy(hname,"Vertex position;Y",nbinsy,ymin,ymax) ;
        sprintf( hname, "hz_nscat%i_vtx%i", nscat, nv ) ;
        TH1F hz(hname,"Vertex position;Z",nbinsz,zmin,zmax) ;
        // fill histograms from ensemble
        for( vector<JPhotonPath>::iterator it=ensemble.begin() ; it!=ensemble.end() ; ++it ) {
          if( it->n-2 != nscat ) continue ;
          JPosition3D vtx( (*it)[nv] ) ;
          hx.Fill( vtx.getX() ) ;
          hy.Fill( vtx.getY() ) ;
          hz.Fill( vtx.getZ() ) ;
        }
        // extra: make sure every field is filled at least once
        for( Int_t bin=1 ; bin<=hx.GetNbinsX() ; ++bin ) 
          hx.AddBinContent(bin,1) ;
        for( Int_t bin=1 ; bin<=hy.GetNbinsX() ; ++bin ) 
          hy.AddBinContent(bin,1) ;
        for( Int_t bin=1 ; bin<=hz.GetNbinsX() ; ++bin ) 
          hz.AddBinContent(bin,1) ;
        
        gens[index] = new JHistGenerator( &hx, &hy, &hz ) ;
      }
    }

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

inlineprotectedvirtual

Implements JMARKOV::JMarkovEnsembleIntegrator.

Definition at line 399 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 ;
    }

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) ;
      }
    }

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 ;
    }

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) ;
    }

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() ) ;
      }
    }

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 ;
      }
    }

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 ;
  }

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 ;
  }

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 ;
  }

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

vector<JHistGenerator*> JMARKOV::JMarkovEnsembleIntegrator1D::gens

protected

Definition at line 408 of file JMarkovIntegrator.hh.

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

protectedinherited

Definition at line 301 of file JMarkovIntegrator.hh.

vector<JPhotonPath> JMARKOV::JMarkovEnsembleIntegrator::ensemble

protectedinherited

Definition at line 302 of file JMarkovIntegrator.hh.

int JMARKOV::JMarkovEnsembleIntegrator::nbinsx

protectedinherited

Definition at line 303 of file JMarkovIntegrator.hh.

int JMARKOV::JMarkovEnsembleIntegrator::nbinsy

protectedinherited

Definition at line 304 of file JMarkovIntegrator.hh.

int JMARKOV::JMarkovEnsembleIntegrator::nbinsz

protectedinherited

Definition at line 305 of file JMarkovIntegrator.hh.

JTargetModel* JMARKOV::JMarkovEnsembleIntegrator::trg

protectedinherited

Definition at line 306 of file JMarkovIntegrator.hh.

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