PMNS_Base.cxx
//.....................................................................
///
/// Compute the sample points for a bin of L/E with width dLoE
///
/// This is used for averaging the probability over a bin of L/E.
/// It should be a private function, but I'm keeping it public for
/// now for debugging purposes. The number of sample points seems
/// too high for most purposes. The number of subdivisions needs
/// to be optimized.
///
/// @param LoE  - The neutrino  L/E value in the bin center in km/GeV
/// @param dLoE   - The L/E bin width in km/GeV
///
Compute the sample points for a bin of L/E with width dLoE

This is used for averaging the probability over a bin of L/E.
It should be a private function, but I'm keeping it public for
now for debugging purposes. The number of sample points seems
too high for most purposes. The number of subdivisions needs
to be optimized.
Compute the sample points for a bin of L / E with width dLoE This is used for averaging the probability over a bin of L / E . It should be a private function but I ' m keeping it public for now for debugging purposes . The number of sample points seems too high for most purposes . The number of subdivisions needs to be optimized .
vector<double> PMNS_Base::GetSamplePoints(double LoE, double dLoE)
{
  SolveHam();
  const double k1267 = kKm2eV / (4 * kGeV2eV);
  vector<double> effDm;
  for(int i=0; i<fNumNus-1; i++){
    for(int j=i+1; j<fNumNus; j++){
      effDm.push_back( 2 * kGeV2eV * fEnergy * fabs(fEval[j] - fEval[i]) );
    }
  }
  int numDm = effDm.size();
  sort(effDm.begin(), effDm.end());
  int n_div = ceil( 20 * pow(dLoE/LoE,0.8) );
  vector<double> allSamples;
  for(int k=0; k<n_div; k++){
    double bctr = LoE - dLoE/2 + (k+0.5)*dLoE/n_div;
    double bwdt = dLoE/n_div;
    vector<double> samples;
    samples.push_back(bctr);
    for(int i=0; i<numDm; i++){
      vector<double> prev = samples;
      double Width = 2*min(prev[0] - (bctr - bwdt/2), (bctr + bwdt/2) - prev[0]);
      const double arg = k1267 * effDm[i] * Width;
      if(i < numDm-1){
        if(arg<0.9) continue;
      }
      else{
        if(arg<0.1) continue;
      }
      samples.clear();
      for(int j=0; j<int(prev.size()); j++){
        double sample = (1/sqrt(3)) * (Width/2);
        if(arg!=0) sample = acos(sin(arg)/arg)/arg * (Width/2);
        samples.push_back(prev[j]-sample);
        samples.push_back(prev[j]+sample);
      }
    }
    allSamples.insert(allSamples.end(), samples.begin(), samples.end());
  }
  return allSamples;
}
vector < double > PMNS_Base :: GetSamplePoints ( double LoE , double dLoE ) { SolveHam ( ) ; const double k1267 = kKm2eV / ( 4 * kGeV2eV ) ; vector < double > effDm ; for ( int i = 0 ; i < fNumNus - 1 ; i ++ ) { for ( int j = i + 1 ; j < fNumNus ; j ++ ) { effDm . push_back ( 2 * kGeV2eV * fEnergy * fabs ( fEval [ j ] - fEval [ i ] ) ) ; } } int numDm = effDm . size ( ) ; sort ( effDm . begin ( ) , effDm . end ( ) ) ; int n_div = ceil ( 20 * pow ( dLoE / LoE , 0.8 ) ) ; vector < double > allSamples ; for ( int k = 0 ; k < n_div ; k ++ ) { double bctr = LoE - dLoE / 2 + ( k + 0.5 ) * dLoE / n_div ; double bwdt = dLoE / n_div ; vector < double > samples ; samples . push_back ( bctr ) ; for ( int i = 0 ; i < numDm ; i ++ ) { vector < double > prev = samples ; double Width = 2 * min ( prev [ 0 ] - ( bctr - bwdt / 2 ) , ( bctr + bwdt / 2 ) - prev [ 0 ] ) ; const double arg = k1267 * effDm [ i ] * Width ; if ( i < numDm - 1 ) { if ( arg < 0.9 ) continue ; } else { if ( arg < 0.1 ) continue ; } samples . clear ( ) ; for ( int j = 0 ; j < int ( prev . size ( ) ) ; j ++ ) { double sample = ( 1 / sqrt ( 3 ) ) * ( Width / 2 ) ; if ( arg != 0 ) sample = acos ( sin ( arg ) / arg ) / arg * ( Width / 2 ) ; samples . push_back ( prev [ j ] - sample ) ; samples . push_back ( prev [ j ] + sample ) ; } } allSamples . insert ( allSamples . end ( ) , samples . begin ( ) , samples . end ( ) ) ; } return allSamples ; }
.....................................................................
Compute the sample points for a bin of L/E with width dLoE
..................................................................... Compute the sample points for a bin of L / E with width dLoE
// Solve Hamiltonian to get eigenvalues // Define conversion factor [km/GeV -> 1/(4 eV^2)] // Get list of all effective Dm^2 // Sort the effective Dm^2 list // Set a number of sub-divisions to achieve "good" accuracy // This needs to be studied better //int n_div = 1; // A vector to store sample points // Loop over sub-divisions // Define sub-division center and width // Make a vector of L/E sample values // Initialized in the sub-division center // Loop over all Dm^2 to average each frequency // This will recursively sample points in smaller // bins so that all relevant frequencies are used // Copy the list of sample L/E values // Redefine bin width to lie within full sub-division // Compute oscillation argument sorted from lowest  to highest // Skip small oscillation values. // If it's the last one, lower the tolerance // Reset samples to redefine them // Loop over previous samples // Compute new sample points around old samples // This is based on a oscillatory quadrature rule // Add samples above and below center // End of loop over Dm^2 // Add sub-division samples to the end of allSamples vector // End of loop over sub-divisions // Return all sample points