Jpp test-rotations-old
the software that should make you happy
Loading...
Searching...
No Matches
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
24namespace JMARKOV {}
25namespace JPP { using namespace JMARKOV; }
26
27namespace 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
268 init_scattering_vertex_generators(nscat) ;
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 JDirection3D 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() ;
321 free_scattering_vertex_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 JDirection3D 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() ;
481 free_scattering_vertex_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
Data structure for direction in three dimensions.
Data structure for position in three dimensions.
const JPosition3D & getPosition() const
Get position.
double getY() const
Get y position.
Definition JVector3D.hh:104
double getZ() const
Get z position.
Definition JVector3D.hh:115
double getX() const
Get x position.
Definition JVector3D.hh:94
double getPhi() const
Get phi angle.
Definition JVersor3D.hh:144
double getDZ() const
Get z direction.
Definition JVersor3D.hh:117
Implementation of the JGenerator interface.
Implementation of the JGenerator interface.
In this implementation of the JMarkovIntegrator interface, the sample distribution is built up out of...
JPhotonPath generatePath(int nscat, double &winv, JPosition3D source_pos, JPosition3D target_pos)
JExperimentalIntegrator(double _lambda)
constructor.
JPosition3D generatePosition(int nscat, int nv, double &winv)
dummy implementation
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...
Implementation of the JGenerator interface.
Abstract interface for the generation of points in 3D space.
Implementation of the JGenerator interface.
Implementation of JMarkovEnsembleIntegrator interface with 1D histograms.
JPosition3D generateScatteringVertexPosition(int nscat, int nv, double &winv)
JMarkovEnsembleIntegrator1D(const vector< JPhotonPath > &_ensemble, JTargetModel *_trg, int _nbinsx, int _nbinsy, int _nbinsz)
void writeHistograms()
write the histograms filled from the ensemble
This implementation of the JMarkovIntegrator interface generates 'natural' looking paths that might s...
JMarkovEnsembleIntegrator3D(const vector< JPhotonPath > &_ensemble, JTargetModel *_trg, int _nbinsx, int _nbinsy, int _nbinsz, JPosition3D _posmin, JPosition3D _posmax)
JPosition3D generateScatteringVertexPosition(int nscat, int nv, double &winv)
void writeHistograms()
write the histograms filled from the ensemble
vector< J3DhistGenerator * > gens
Abstract base class for implementations of the JMarkovIntegrator interface, where the sample distribu...
virtual void writeHistograms()=0
write the histograms filled from the ensemble
virtual JPosition3D generateScatteringVertexPosition(int nscat, int nv, double &winv)=0
virtual void init_scattering_vertex_generators(int nscat)=0
JMarkovEnsembleIntegrator(const vector< JPhotonPath > &_ensemble, JTargetModel *_trg, int _nbinsx, int _nbinsy, int _nbinsz)
void init_target_generator(int nscat)
initialize the generator for a target vertex for a given number of scatterings
void initgenerators(int nscat)
initialize the position generators for a given number of scatterings (i.e. create histograms and fill...
int getIndex(int nv, int nscat)
return the internal index for a given number of scatterings and vertex number
vector< JSphereGenerator * > target_gens
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
Abstract base class for calculating the total probability (/m^2 target cross-section) for a photon fr...
virtual JPosition3D generatePosition(int nscat, int nv, double &winv)=0
Generate a random position for vertex nv.
JMarkovIntegrator()
standard constructor
vector< double > integrate(int N, int nscat, JSourceModel *src, JScatteringModel *sm, JTargetModel *trg, double lambda_abs)
Integrate with N samples.
vector< double > dummy_integrate(int N, int nscat, JSourceModel *src, JScatteringModel *sm, JTargetModel *trg, double lambda_abs)
Integrate a test function with N samples.
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.
virtual JPhotonPath generatePath(int nscat, double &winv)
Generate a random photon path with a given number of scatterings.
In this implementation of the JMarkovIntegrator interface, the sample distribution is a flat distribu...
JMarkovUniformIntegrator(JPosition3D _posmin, JPosition3D _posmax)
constructor.
virtual JPosition3D generatePosition(int nscat, int nv, double &winv)
Generate a random position for vertex nv.
A photon path.
Virtual base class for a scattering model.
Implementation of the JGenerator interface.
Implementation of the JGenerator interface.
Virtual base class for a light source.
In this implementation of the JMarkovIntegrator interface, the sample distribution is built up out of...
JSourceTargetIntegrator(double Rsing, JPosition3D src_pos, JPosition3D trg_pos, double fraction_exp, double Lexp, double Rexp)
constructor.
virtual JPosition3D generatePosition(int nscat, int nv, double &winv)
Generate a random position for vertex nv.
Implementation of the JGenerator interface.
Virtual base class for a light detector ("photon target").
Implementation of the JGenerator interface.
Implementation of the JGenerator interface.
double getWeight(JPosition3D pos)
return the weight (=probability density dP/dV) for the given position.
JPosition3D getPosition()
Return a randomly generated 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.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).