FreePathSolver Class Reference

Finds photon paths from a nanobeacon source to a target PMT that have a non-zero probability. More...

Public Member Functions
	FreePathSolver ()
JPhotonPath	solve (int nscat, const JDirectedSource *src, const JPMTTarget *trg)

Protected Member Functions
double	sourceScore (const JPhotonPath &p, const JDirectedSource *src)
double	targetScore (const JPhotonPath &p, const JPMTTarget *trg)
double	intersectScore (const JPhotonPath &p, const JPMTTarget *trg)
double	getScore (const JPhotonPath &p, const JDirectedSource *src, const JPMTTarget *trg)

Detailed Description

Finds photon paths from a nanobeacon source to a target PMT that have a non-zero probability.

This means that

• the emission angle has to be within the nanobeacon opening angle
• the photon has to hit the PMT at an angle less than 90 degrees
• the photon may not intersect the target

It works by maximizing a score function. If we find a score of 0, the path meets all criteria.

Definition at line 162 of file JBeaconSimulator.cc.

Constructor & Destructor Documentation

FreePathSolver::FreePathSolver ( )

inline

Definition at line 165 of file JBeaconSimulator.cc.

{}

Member Function Documentation

JPhotonPath FreePathSolver::solve ( int  nscat, const JDirectedSource *  src, const JPMTTarget *  trg  )

inline

Definition at line 167 of file JBeaconSimulator.cc.

                                                                                   {
    //cout << "Running FreePathSolver::solve with nscat = " << nscat << endl ;

    JPhotonPath p(min(2,nscat)) ;

    // start at source
    p.front() = src->getPosition() ; 
    // end on PMT surface
    p.back()  = trg->getPosition() ;
    if( trg->getRadius()>0 ) p.back() = trg->getPosition() + trg->getRadius() * JVector3D(trg->getDirection()) ;

    if( nscat == 0 ) return p ;

    // create intermediate points
    for( int nv=1 ; nv<(int)p.size()-1 ; ++nv ) {
      double factor = 1.0 * nv / (p.size()-1) ;
      p[nv] = (1-factor)*p.front() + factor*p.back() ;
    }

    int maxnsteps = -1 ;

    if( nscat == 1 ) {
      double score = getScore(p,src,trg) ;
      double dl = 1 ; // m
      const double dlmin = 0.0001 ;
      JPosition3D dx(1,0,0) ;
      JPosition3D dy(0,1,0) ;
      JPosition3D dz(0,0,1) ;
      int istep = 0 ;
      while( score<0 ) {
        if( ++istep == maxnsteps ) break ;
        double initscore = score ;

        // try step in x-direction
        p[1] = p[1] + dl*dx ;
        if( getScore(p,src,trg) < score ) {
          // try moving in the opposite direction
          p[1] = p[1] - 2*dl*dx ;
          if( getScore(p,src,trg) < score ) p[1] = p[1] + dl*dx ;
        }

        // try step in y-direction
        p[1] = p[1] + dl*dy ;
        if( getScore(p,src,trg) < score ) {
          // try moving in the opposite direction
          p[1] = p[1] - 2*dl*dy ;
          if( getScore(p,src,trg) < score ) p[1] = p[1] + dl*dy ; 
        }

        // try step in z-direction
        p[1] = p[1] + dl*dz ;
        if( getScore(p,src,trg) < score ) {
          // try moving in the opposite direction
          p[1] = p[1] - 2*dl*dz ;
          if( getScore(p,src,trg) < score ) p[1] = p[1] + dl*dz ;
        }

        // see if it worked
        score = getScore(p,src,trg) ;
        // check if we are stuck
        if( score == initscore ) {
          dl = 0.5 * dl ;
          //cout << "dl = " << dl << endl ;
          if( dl<dlmin ) break ;
        }
        if( p.getLength() > 10000 ) break ;
      }
    }

    if( nscat>1) {
      //cout << "nscat more than one" << endl ;
      //cout << "starting with " << endl ;
      //for( vector<JPosition3D>::iterator it=p.begin() ; it!= p.end() ; ++it ) cout << *it << endl ;
      //cout << endl ;

      double score = getScore(p,src,trg) ;
      double dl = 1 ; // m
      const double dlmin = 1e-10 ;

      const int nsm = 6 ;
      JPosition3D single_moves[nsm] = { JPosition3D(1,0,0), JPosition3D(0,1,0), JPosition3D(0,0,1),
                                        JPosition3D(-1,0,0), JPosition3D(0,-1,0), JPosition3D(0,0,-1) } ;
      /*ctor< pair<JPosition3D,JPosition3D> > moves ;
        for( int i=0 ; i<nsm ; ++i ) {
        for( int j=0 ; j<nsm ; ++j ) {
        moves.push_back( pair<JPosition3D,JPosition3D>(single_moves[i],single_moves[j]) ) ;
        }
        }*/

      int i = 0 ;
      
      while( score<0 ) {
        if( ++i>100 ) break ;
        //cout << "@ " << JPosition3D(p[1]-src->getPosition()) << ", " << JPosition3D(p[2]-src->getPosition()) << ", score = " << score << ", dl = " << dl << endl ;
        double initscore = score ;

        for( int i=0; i<nsm; ++i ) {
          // try move
          p[1] = p[1] + dl * single_moves[i]  ;
          double newscore = getScore(p,src,trg) ;
          if( newscore < score ) {
            // if the score is reduced, undo the move
            p[1] = p[1] - dl * single_moves[i]  ;
          } else {
            //cout << "First vertex adds " << single_moves[i] 
            //   << ", score difference = " << newscore - score << endl ;
            score = newscore ;
          }
        }

        for( int i=0; i<nsm; ++i ) {
          // try move
          p[2] = p[2] + dl * single_moves[i]  ;
          double newscore = getScore(p,src,trg) ;
          if( newscore < score ) {
            // if the score is reduced, undo the move
            p[2] = p[2] - dl * single_moves[i]  ;
          } else {
            //cout << "Second vertex adds " << single_moves[i] 
            //   << ", score difference = " << newscore - score << endl ;
            score = newscore ;
          }
        }

        // see if it worked
        score = getScore(p,src,trg) ;
        // check if we are stuck
        if( score == initscore ) {
          //cout << "Stuck" << endl ;
          dl = 0.5*dl ;
          if( dl<dlmin ) break ;
        }
        if( p.getLength() > 10000 ) {
          //cout << "Too long" << endl ;
          break ;
        }
      }

      // fill in the gaps
      JPhotonPath p2(nscat) ;
      p2[0] = p[0] ;
      p2[1] = p[1] ;
      p2[2] = p[2] ;
      p2.back() = p.back() ;
      for( int nv=3 ; nv<nscat+1 ; ++nv ) {
        double factor = 1.0 * (nv-2) / (nscat-1) ;
        p2[nv] = (1-factor)*p2[2] + factor*p2.back() ;
      }
      return p2 ;
    }
  
    return p ;
  
  }

double FreePathSolver::sourceScore ( const JPhotonPath &  p, const JDirectedSource *  src  )

inlineprotected

Definition at line 324 of file JBeaconSimulator.cc.

                                                                         {
    double ct = JVersor3D(p[1]-src->getPosition()).getDot(src->getDirection()) ;
    const double ctmin = cos(0.5*src->getOpeningAngle()) ;
    double val =  max(ctmin-ct,0.0) / (1-ctmin) ;
    return -val*val ;
  }

double FreePathSolver::targetScore ( const JPhotonPath &  p, const JPMTTarget *  trg  )

inlineprotected

Definition at line 331 of file JBeaconSimulator.cc.

                                                                    {
    double ct = JVersor3D(p[p.size()-2]-p[p.size()-1]).getDot(trg->getDirection()) ;
    const double ctmin = cos(0.5*trg->getOpeningAngle()) ;
    double val =  max(ctmin-ct,0.0) / (1-ctmin) ;
    return -val*val ;
  }

double FreePathSolver::intersectScore ( const JPhotonPath &  p, const JPMTTarget *  trg  )

inlineprotected

Definition at line 338 of file JBeaconSimulator.cc.

                                                                       {
    double score = 0 ;

    // check whether the photon path goes through the target sphere
    // we do not need to consider the last vertex, because if that intersects the sphere
    // the target score would be 0.
    if( p.n == 2 ) return score ;

    for( int nv=1 ; nv<(int)p.size()-2 ; ++nv ) {
      JPhotonPath seg(0) ;
      seg[0] = p[nv] ;
      seg[1] = p[nv+1] ;
      if( seg.hitsSphere(trg->getPosition(),trg->getRadius()) ) {
        double d = JAxis3D( seg[0], JVersor3D(seg[1]-seg[0]) ).getDistance(trg->getPosition()) ;
        double R = trg->getRadius() ;
        double val = (R-d)/R ;
        score -= val*val ;
        //cout << "Segment from " << nv << " to " << nv+1 << "intersects sphere"
        //     << ", distance = " << d << ", val = " << val << endl ;
      }
    }

    /*JPhotonPath pshort(p) ;
    pshort.resize( pshort.size()-1 ) ;
    --pshort.n ;
    if( pshort.hitsSphere(trg->getPosition(),trg->getRadius()) ) {
      JVersor3D hit_dir( pshort.getSphereHitPosition(trg->getPosition(),trg->getRadius()) - trg->getPosition() ) ;
      double val = hit_dir.getDot(trg->getDirection()) ;
      score -= val * val ; 
    }
    
    // extra penalty when a vertex is inside the target
    for( int nv=1 ; nv<p.n-1 ; ++nv ) {
      double r =  p[nv].getDistance(trg->getPosition()) ;
      if( r < trg->getRadius() ) {
        double val = ( trg->getRadius() - r ) / trg->getRadius() ;
        score -= val * val ;
      } 
      }*/

    return score ;
  }

double FreePathSolver::getScore ( const JPhotonPath &  p, const JDirectedSource *  src, const JPMTTarget *  trg  )

inlineprotected

Definition at line 382 of file JBeaconSimulator.cc.

                                                                                              {
    double score = sourceScore(p,src) + targetScore(p,trg) + intersectScore(p,trg) ;
    return score ;
  }

---

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

• JBeaconSimulator.cc