Jpp  18.0.0-rc.2
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JMarkovIntegrator.hh
Go to the documentation of this file.
1 #ifndef H_MARKOV_INTEGRATOR
2 #define H_MARKOV_INTEGRATOR
3 
4 /**
5  * \author mjongen
6  */
7 
8 // c++ standard library includes
9 #include<iostream>
10 #include<iomanip>
11 #include<cmath>
12 #include<cstdlib>
13 
14 // root includes
15 #include "TObject.h"
16 #include "TH1F.h"
17 #include "TH3F.h"
18 #include "TRandom3.h"
19 
20 // JPP includes
23 
24 namespace JMARKOV {}
25 namespace JPP { using namespace JMARKOV; }
26 
27 namespace JMARKOV {
28 
29  /**
30  Abstract base class for calculating the total probability
31  (/m^2 target cross-section) for a photon from the source to hit
32  the target (with a given, fixed number of scatterings) by Monte
33  Carlo sampling the available nscat*3 dimensional phase space.
34 
35  The sample distribution is implemented in derived classes.
36 
37  Random numbers are drawn from gRandom.
38  **/
40  public:
41 
42  /// standard constructor
44 
45  /**
46  Integrate with N samples.
47 
48  Returns a vector with the contribution to the integral of each sample.
49  The mean of those values is the estimate of the result of the integral,
50  while the distribution itself can be used to estimate the stability of
51  the result.
52  In this distribution, you want to avoid
53  - long tails (because they make the result unstable)
54  - small contributions (because it means that parts of the parameter
55  space are being oversampled, so it is less efficient)
56  This can be achieved by tuning the sample distribution to the problem
57  at hand.
58  **/
59  vector<double> integrate( int N, int nscat, JSourceModel* src, JScatteringModel* sm, JTargetModel* trg, double lambda_abs ) ;
60 
61  /**
62  Integrate a test function with N samples.
63 
64  This can be used as a sanity check for derived classes of
65  JMarkovIntegrator.
66 
67  The integral should yield 1 when the complete relevant part of the
68  volume is taken into account.
69  If it does not, it may be a sign that the implementation is not
70  correct.
71 
72  Returns a vector with the contribution to the integral of each sample.
73  **/
74  vector<double> dummy_integrate( int N, int nscat, JSourceModel* src, JScatteringModel* sm, JTargetModel* trg, double lambda_abs ) ;
75 
76  /**
77  Return photon paths generated with the generatePath method.
78 
79  This can be used to identify the parts of parameter space that are
80  over- or undersampled in a given problem so that the integrator may
81  be optimized to handle those better.
82  **/
83  vector<JPhotonPath> get_diagnostic_ensemble( int N, int nscat, JSourceModel* src, JScatteringModel* sm, JTargetModel* trg, double lambda_abs ) ;
84 
85  protected:
86 
87  /**
88  Generate a random photon path with a given number of
89  scatterings
90 
91  winv must be set to the inverted probability density
92  to generate this particular path.
93  **/
94  virtual JPhotonPath generatePath( int nscat, double& winv ) {
95  // this default implementation assumes that the vertex positions
96  // are completely uncorrelated
97  JPhotonPath p(nscat) ;
98  double _winv = 1 ;
99  for( int nv=0 ; nv<nscat+2 ; ++nv ) {
100  double part_winv ;
101  p[nv] = generatePosition( nscat, nv, part_winv ) ;
102  _winv *= part_winv ;
103  }
104  winv = _winv ;
105  return p ;
106  }
107 
108 
109  /**
110  Generate a random position for vertex nv
111 
112  nv = 1 is the first scattering vertex,
113  nv = 2 is the second scattering vertex
114  etc.
115 
116  So we require that
117  0 < nv <= nscat
118 
119  winv has to be set to 1 over the probability density
120  for this particular position.
121  NOTE: this only works when the positions of each vertex are
122  UNCORRELATED to each other (as in the default implementation
123  of generatePath. If they are not, this method can
124  be given a dummy implementation and be ignored.
125  **/
126  virtual JPosition3D generatePosition( int nscat, int nv, double& winv ) = 0 ;
127 
128  } ;
129 
130  vector<double> JMarkovIntegrator::integrate( int N, int nscat, JSourceModel* src, JScatteringModel* sm, JTargetModel* trg, double lambda_abs ) {
131  vector<double> contributions(N,-1) ;
132 
133  for( int i=0 ; i<N ; ++i ) {
134  double winv ;
135  JPhotonPath p = generatePath(nscat,winv) ;
136 
137  double rho = getPhotonPathProbabilityDensity(p,src,sm,trg,lambda_abs) ;
138  contributions[i] = rho * winv ;
139  }
140  return contributions ;
141  }
142 
143  vector<double> JMarkovIntegrator::dummy_integrate( int N, int nscat, JSourceModel* src, JScatteringModel* sm, JTargetModel* trg, double lambda_abs ) {
144  vector<double> contributions(N,-1) ;
145 
146  for( int i=0 ; i<N ; ++i ) {
147  const double r = 10 ;
148  double winv ;
149  JPhotonPath p = generatePath(nscat,winv) ;
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 ;
153  }
154  return contributions ;
155  }
156 
158  vector<JPhotonPath> paths ;
159 
160  for( int i=0 ; i<N ; ++i ) {
161  double winv ;
162  JPhotonPath p = generatePath(nscat,winv) ;
163  paths.push_back(p) ;
164  }
165  return paths ;
166  }
167 
168 
169  /**
170  In this implementation of the JMarkovIntegrator interface,
171  the sample distribution is a flat distribution.
172  **/
174  public :
175 
176  /**
177  constructor. Posmin and posmax contain the minimal and maximal
178  values of the coordinate in each direction.
179  **/
180  JMarkovUniformIntegrator( JPosition3D _posmin, JPosition3D _posmax ) : gen(_posmin,_posmax) {}
181 
182  protected:
183 
184  virtual JPosition3D generatePosition( int nscat, int nv, double& winv ) {
185  JPosition3D pos = gen.getPosition() ;
186  winv = 1.0/gen.getWeight(pos) ;
187  return pos ;
188  }
189 
191  } ;
192 
193  /**
194  Abstract base class for implementations of the JMarkovIntegrator interface,
195  where the sample distribution is based on histograms filled from an ensemble
196  of representative paths.
197  **/
199  public :
200 
201  // constructor
202  JMarkovEnsembleIntegrator( const vector<JPhotonPath>& _ensemble, JTargetModel* _trg, int _nbinsx, int _nbinsy, int _nbinsz ) : ensemble(_ensemble), nbinsx(_nbinsx), nbinsy(_nbinsy), nbinsz(_nbinsz), trg(_trg) {
203  // check that the ensemble size is at least one photon path
204  if( ensemble.size()==0 ) {
205  cerr << "FATAL ERROR in JMarkov3DensembleIntegrator constructor. Ensemble size = 0." << endl ;
206  exit(1) ;
207  }
208 
209  // check that the position of the first vertex is the same throughout the ensemble
210  for( vector<JPhotonPath>::iterator it=ensemble.begin(); it!=ensemble.end(); ++it ) {
211  if( (*it)[0].getDistance(ensemble.front()[0]) > 0 ) {
212  cerr << "Fatal error in JMarkovEnsembleIntegrator, first vertex positions are not all the same." << endl ;
213  exit(1) ;
214  }
215  }
216  // set source position based on ensemble
217  source_position = ensemble.front()[0] ;
218  }
219 
221 
222  /// write the histograms filled from the ensemble
223  virtual void writeHistograms() = 0 ;
224 
225  protected :
226 
227  /// generate a position for vertex nv for the given number of scatterings, winv will be set to 1/weight
228  JPosition3D generatePosition( int nscat, int nv, double& winv ) {
229  // source vertex
230  if( nv==0 ) {
231  winv = 1.0 ;
232  return source_position ;
233  }
234  // target vertex
235  if( nv == nscat+1 ) {
236  if( nscat+1>(int)target_gens.size() || target_gens[nscat] == NULL ) {
237  initgenerators(nscat) ;
238  }
239  JPosition3D pos = target_gens[nscat]->getPosition() ;
240  winv = 1.0/target_gens[nscat]->getWeight(pos) ;
241  return pos ;
242  } else {
243  // scattering vertex
244  return generateScatteringVertexPosition(nscat,nv,winv) ;
245  }
246  }
247 
248  virtual JPosition3D generateScatteringVertexPosition(int nscat, int nv, double& winv) = 0 ;
249 
250  /// return the internal index for a given number of scatterings and vertex number
251  int getIndex( int nv, int nscat ) {
252  // indices are distributed as follows
253  // [0]: nscat=1, nv=1
254  // [1]: nscat=2, nv=1
255  // [2]: nscat=2, nv=2
256  // [3]: nscat=3, nv=1
257  // [4]: nscat=3, nv=2
258  // etc.
259  return (nscat*(nscat-1))/2 + nv - 1 ;
260  }
261 
262  /// initialize the position generators for a given number of scatterings (i.e. create histograms and fill them from the ensemble)
263  void initgenerators( int nscat ) {
264  // target vertex
265  init_target_generator(nscat) ;
266 
267  // scattering vertices
269  }
270 
271  /// initialize the generator for a target vertex for a given number of scatterings
272  void init_target_generator( int nscat ) {
273  if( (int)target_gens.size() <= nscat ) target_gens.resize( nscat+1 ) ;
274 
275  if( trg->getRadius()>0 ) {
276  // for a finite target, we fill a "sphere generator" from the ensemble
277  char hname[200] ;
278  sprintf( hname, "htarget_nscat%i", nscat ) ;
279  TH2D* ht = new TH2D(hname,";cos(#theta);#phi",100,-1,1,100,-M_PI,M_PI) ;
280  for( vector<JPhotonPath>::iterator it=ensemble.begin() ; it!=ensemble.end() ; ++it ) {
281  if( it->n-2 != nscat ) continue ;
282  JVersor3D dir( it->back() - trg->getPosition() ) ;
283  ht->Fill( dir.getDZ(), dir.getPhi() ) ;
284  }
285  target_gens[nscat] = new JSphereGenerator( trg->getPosition(), trg->getRadius(), ht ) ;
286  delete ht ;
287  } else {
288  // for an infinitesimal target, we create a "generator" that only generates the single point
289  target_gens[nscat] = new JSphereGenerator( trg->getPosition() ) ;
290  }
291  }
292 
293  virtual void init_scattering_vertex_generators( int nscat ) = 0 ;
294 
296  for( vector<JSphereGenerator*>::iterator it=target_gens.begin() ; it!=target_gens.end() ; ++it ) {
297  if( *it != NULL ) delete *it ;
298  }
299  }
300 
303  int nbinsx ;
304  int nbinsy ;
305  int nbinsz ;
308  } ;
309 
310  /**
311  Implementation of JMarkovEnsembleIntegrator interface with 1D histograms.
312  **/
314  public :
315 
316  // constructor
317  JMarkovEnsembleIntegrator1D( const vector<JPhotonPath>& _ensemble, JTargetModel* _trg, int _nbinsx, int _nbinsy, int _nbinsz ) : JMarkovEnsembleIntegrator(_ensemble,_trg,_nbinsx,_nbinsy,_nbinsz) {}
318 
320  free_target_gens() ;
322  }
323 
325  for( vector<JHistGenerator*>::iterator it=gens.begin() ; it!=gens.end() ; ++it ) {
326  if( *it != NULL ) delete *it ;
327  }
328  }
329 
331  for( vector<JHistGenerator*>::iterator it=gens.begin() ; it!=gens.end() ; ++it ) {
332  if( *it != NULL ) {
333  (*it)->hx->Write() ;
334  (*it)->hy->Write() ;
335  (*it)->hz->Write() ;
336  }
337  }
338  for( vector<JSphereGenerator*>::iterator it=target_gens.begin(); it!=target_gens.end(); ++it ) {
339  if( *it != NULL ) {
340  if( (*it)->h != NULL ) (*it)->h->Write() ;
341  }
342  }
343  }
344 
345  protected :
346 
348  int ifirst = getIndex(1,nscat) ;
349  int ilast = getIndex(1,nscat+1) ;
350  if( (int)gens.size() < ilast ) gens.resize( ilast, NULL ) ;
351  for( int nv=1 ; nv<=nscat ; ++nv ) {
352  int index = ifirst + nv - 1 ;
353  // determine minimal and maximal coordinate values for this vertex
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 ;
360  for( vector<JPhotonPath>::iterator it=ensemble.begin() ; it!=ensemble.end() ; ++it ) {
361  if( it->n-2 != nscat ) continue ;
362  JPosition3D vtx( (*it)[nv] ) ;
363  if( vtx.getX()<xmin ) xmin = vtx.getX() ;
364  if( vtx.getX()>xmax ) xmax = vtx.getX() ;
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() ;
369  }
370  // allocate histograms
371  char hname[200] ;
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) ;
378  // fill histograms from ensemble
379  for( vector<JPhotonPath>::iterator it=ensemble.begin() ; it!=ensemble.end() ; ++it ) {
380  if( it->n-2 != nscat ) continue ;
381  JPosition3D vtx( (*it)[nv] ) ;
382  hx.Fill( vtx.getX() ) ;
383  hy.Fill( vtx.getY() ) ;
384  hz.Fill( vtx.getZ() ) ;
385  }
386  // extra: make sure every field is filled at least once
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) ;
393 
394  gens[index] = new JHistGenerator( &hx, &hy, &hz ) ;
395  }
396  }
397 
398 
399  JPosition3D generateScatteringVertexPosition(int nscat, int nv, double& winv) {
400  int index = getIndex(nv,nscat) ;
401  if( index+1>(int)gens.size() || gens[index] == NULL ) initgenerators(nscat) ;
402 
403  JPosition3D pos = gens[index]->getPosition() ;
404  winv = 1.0/gens[index]->getWeight(pos) ;
405  return pos ;
406  }
407 
409  } ;
410 
411  /**
412  This implementation of the JMarkovIntegrator interface generates
413  'natural' looking paths that might sample the phase space well
414  in some cases
415  **/
416  /*class JMarkovNaturalIntegrator : public JMarkovIntegrator {
417  public :
418 
419  JMarkovNaturalIntegrator( double _l ) {
420  l = _l ; // approximate average path length above the minimum
421  }
422 
423  ~JMarkovNaturalIntegrator() {}
424 
425  protected :
426 
427  /// dummy implementation
428  JPosition3D generatePosition( int nscat, int nv, double& winv ) {
429  return JPosition3D(0,0,0) ;
430  } ;
431 
432  /// Generate a 'natural' path
433  JPhotonPath generatePath( int nscat, double& winv ) {
434  double winvtot = 1 ;
435 
436  // distance from source to target
437  double d = (source_pos-target_pos).getLength() ;
438 
439  double g = d / (l+d) ;
440  JHenyeyGreensteinScattering sm(100,g) ;
441 
442  double lambda = (d+l)/(nscat+1) ;
443 
444  JPhotonPath p(nscat) ;
445  p.front() = source_pos ;
446  p.back() = target_pos ;
447 
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) ;
457  // set vertex
458  p[nv] = p[nv-1] + r*JVector3D(dir) ;
459  // double
460  winvtot /= ( r * r * sm.getScatteringProbability(ct) * 1.0/lambda*exp(-r/lambda) ) ;
461  }
462 
463  return p ;
464  }
465 
466  double l ;
467  } ;*/
468 
469  /**
470  Similar to JMarkovEnsembleIntegrator, but using a 3D histogram for each
471  vertex instead of 3 1D histograms.
472  **/
474  public :
475 
476  // constructor
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) {}
478 
480  free_target_gens() ;
482  }
483 
485  // free memory
486  for( vector<J3DhistGenerator*>::iterator it=gens.begin() ; it!=gens.end() ; ++it ) {
487  if( *it != NULL ) delete *it ;
488  }
489  }
490 
492  for( vector<J3DhistGenerator*>::iterator it=gens.begin() ; it!=gens.end() ; ++it ) {
493  if( *it != NULL ) {
494  (*it)->h->Write() ;
495  TH2D* htmp ;
496  htmp = (TH2D*)(*it)->h->Project3D("yx") ;
497  htmp->SetOption("colz") ;
498  htmp->SetEntries(10) ;
499  htmp->Write() ;
500  htmp = (TH2D*)(*it)->h->Project3D("yz") ;
501  htmp->SetOption("colz") ;
502  htmp->SetEntries(10) ;
503  htmp->Write() ;
504  htmp = (TH2D*)(*it)->h->Project3D("zx") ;
505  htmp->SetOption("colz") ;
506  htmp->SetEntries(10) ;
507  htmp->Write() ;
508  }
509  }
510  }
511 
512  JPosition3D generateScatteringVertexPosition(int nscat, int nv, double& winv) {
513  int index = getIndex(nv,nscat) ;
514  if( index+1>(int)gens.size() || gens[index] == NULL ) initgenerators(nscat) ;
515 
516  JPosition3D pos = gens[index]->getPosition() ;
517  winv = 1.0/gens[index]->getWeight(pos) ;
518  return pos ;
519  }
520 
521  protected :
522 
524  int ifirst = getIndex(1,nscat) ;
525  int ilast = getIndex(1,nscat+1) ;
526  if( (int)gens.size() < ilast ) gens.resize( ilast, NULL ) ;
527 
528  for( int nv=1 ; nv<=nscat ; ++nv ) {
529  int index = ifirst + nv - 1 ;
530  // determine minimal and maximal value
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 ;
537  xmin = posmin.getX() ;
538  ymin = posmin.getY() ;
539  zmin = posmin.getZ() ;
540  xmax = posmax.getX() ;
541  ymax = posmax.getY() ;
542  zmax = posmax.getZ() ;
543  // allocate histogram
544  char hname[200] ;
545  sprintf( hname, "h3_nscat%i_vtx%i", nscat, nv ) ;
546  TH3F h(hname,"Vertex position;X;Y;Z",
547  nbinsx,xmin,xmax,
548  nbinsy,ymin,ymax,
549  nbinsz,zmin,zmax ) ;
550  // fill histogram from ensemble
551  for( vector<JPhotonPath>::iterator it=ensemble.begin() ; it!=ensemble.end() ; ++it ) {
552  if( it->n-2 != nscat ) continue ;
553  JPosition3D vtx( (*it)[nv] ) ;
554  h.Fill( vtx.getX(), vtx.getY(), vtx.getZ() ) ;
555  }
556  // extra: make sure every field is filled at least once
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) ;
563  }
564  }
565  }
566  gens[index] = new J3DhistGenerator( &h ) ;
567  }
568  }
569 
573  } ;
574 
575  /**
576  In this implementation of the JMarkovIntegrator interface,
577  the sample distribution is built up out of three components:
578  - a 1/r^2 distribution centered around the source
579  - a 1/r^2 distribution centered around the target
580  - an exponential distribution around the point directly between the source and the target
581  **/
583  public :
584 
585  /**
586  constructor.
587 
588  Arguments:
589  - position of the source
590  - position of the target
591  - the fraction of points that is drawn from the exponential distribution
592  - the minimal radius of the exp. distr.
593  - the maximal radius of the exp. distr.
594  - the length parameter of the exp. distr.
595  - the maximal radius for the 1/r^2 distributions
596  **/
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) {
598  JPosition3D cnt_pos = 0.5*(src_pos+trg_pos) ;
599  gen_exp_shift = new JShiftedGenerator(&gen_exp,cnt_pos) ;
600  double fraction_sing = 0.5*(1-fraction_exp) ;
601  gen = new JTripleGenerator( fraction_sing, &gen_s1, fraction_sing, &gen_s2, fraction_exp, gen_exp_shift ) ;
602  }
603 
605  delete gen_exp_shift ;
606  delete gen ;
607  }
608 
609  protected:
610 
611  virtual JPosition3D generatePosition( int nscat, int nv, double& winv ) {
612  JPosition3D pos = gen->getPosition() ;
613  winv = 1.0/gen->getWeight(pos) ;
614  return pos ;
615  }
616 
622  } ;
623 
624  /**
625  In this implementation of the JMarkovIntegrator interface,
626  the sample distribution is built up out of correlated path vertices.
627 
628 
629  lambda is the effective length scale defined as:
630 
631  1/lambda := 1/lambda_scat + 1/lambda_abs
632  **/
634  public :
635 
636  /**
637  constructor.
638  **/
639  JExperimentalIntegrator( double _lambda ) : gen(_lambda) {}
640 
642 
643  protected:
644 
645  /// dummy implementation
646  JPosition3D generatePosition( int nscat, int nv, double& winv ) {
647  return JPosition3D(0,0,0) ;
648  } ;
649 
650  JPhotonPath generatePath( int nscat, double& winv, JPosition3D source_pos, JPosition3D target_pos ) {
651  JPhotonPath p(nscat) ;
652  p.front() = source_pos ;
653  p.back() = target_pos ;
654  double _winv = 1 ;
655 
656  for( int nv=1 ; nv<=nscat ; ++nv ) {
657  double __winv ;
658  p[nv] = getFirstScatteringVertex(nscat+1-nv,p[nv-1],p.back(),__winv ) ;
659  _winv *= __winv ;
660  }
661  winv = _winv ;
662  return p ;
663  }
664 
665  protected:
666 
667  /**
668  recursive function to randomly generate the first scattering
669  vertex for isotropic scattering from a source to a target
670  **/
671  JPosition3D getFirstScatteringVertex( int nscat, JPosition3D src_pos, JPosition3D trg_pos, double& winv ) {
672  double _winv = 1 ;
673  if( nscat == 0 ) {
674  winv = 1 ;
675  return trg_pos ;
676  }
677  if( gRandom->Integer(2) == 0 ) {
678  src_pos = getFirstScatteringVertex(nscat-1,src_pos,trg_pos,_winv) ;
679  }
680  JPosition3D pos = gen.getPosition() ;
681  winv = 1.0 / gen.getWeight(pos) ;
682  return src_pos + pos ;
683  }
684 
686  } ;
687 
688 }
689 
690 #endif
691 
const double xmax
Definition: JQuadrature.cc:24
Implementation of the JGenerator interface.
Implementation of the JGenerator interface.
virtual JPosition3D generateScatteringVertexPosition(int nscat, int nv, double &winv)=0
JPosition3D getPosition()
Return a randomly generated position.
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...
In this implementation of the JMarkovIntegrator interface, the sample distribution is a flat distribu...
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
double getRadius() const
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)
A photon path.
Definition: JPhotonPath.hh:38
vector< JSphereGenerator * > target_gens
exit
Definition: JPizza.sh:36
double getWeight(JPosition3D pos)
return the weight (=probability density dP/dV) for the given position.
data_type r[M+1]
Definition: JPolint.hh:779
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.
JExperimentalIntegrator(double _lambda)
constructor.
Abstract interface for the generation of points in 3D space.
skip continue
Definition: JTuneHV.sh:92
JMarkovEnsembleIntegrator(const vector< JPhotonPath > &_ensemble, JTargetModel *_trg, int _nbinsx, int _nbinsy, int _nbinsz)
This implementation of the JMarkovIntegrator interface generates &#39;natural&#39; looking paths that might s...
virtual void writeHistograms()=0
write the histograms filled from the ensemble
const JPosition3D & getPosition() const
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.
Virtual base class for a light source.
Implementation of the JGenerator interface.
Implementation of the JGenerator interface.
double getY() const
Get y position.
Definition: JVector3D.hh:104
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.
const double xmin
Definition: JQuadrature.cc:23
then usage $script< input file >[option[primary[working directory]]] nWhere option can be N
Definition: JMuonPostfit.sh:36
Implementation of the JGenerator interface.
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.
double getX() const
Get x position.
Definition: JVector3D.hh:94
Virtual base class for a scattering model.
Implementation of the JGenerator interface.
Virtual base class for a light detector (&quot;photon target&quot;).
virtual JPosition3D generatePosition(int nscat, int nv, double &winv)
Generate a random position for vertex nv.
Data structure for position in three dimensions.
Definition: JPosition3D.hh:36
JMarkovUniformIntegrator(JPosition3D _posmin, JPosition3D _posmax)
constructor.
Data structure for normalised vector in three dimensions.
Definition: JVersor3D.hh:26
vector< J3DhistGenerator * > gens
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.
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.
Definition: JVector3D.hh:115
Abstract base class for calculating the total probability (/m^2 target cross-section) for a photon fr...
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.