JMARKOV::JMarkovEnsembleIntegrator3D Class Reference

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

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.

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

JMARKOV::JMarkovEnsembleIntegrator3D::~JMarkovEnsembleIntegrator3D ( )

inline

Definition at line 479 of file JMarkovIntegrator.hh.

                                   {
      free_target_gens() ;
      free_scattering_vertex_gens() ;
    }

Member Function Documentation

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

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

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

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

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 ;
          JDirection3D 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<J3DhistGenerator*> JMARKOV::JMarkovEnsembleIntegrator3D::gens

protected

Definition at line 570 of file JMarkovIntegrator.hh.

JPosition3D JMARKOV::JMarkovEnsembleIntegrator3D::posmin

protected

Definition at line 571 of file JMarkovIntegrator.hh.

JPosition3D JMARKOV::JMarkovEnsembleIntegrator3D::posmax

protected

Definition at line 572 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