Changeset View
Changeset View
Standalone View
Standalone View
src/EvtGenBase/EvtLinSample.cpp
- This file was added.
#include "EvtGenBase/EvtLinSample.hh" | |||||
#include <algorithm> | |||||
#include <cmath> | |||||
#include <fstream> | |||||
#include <iostream> | |||||
EvtLinSample::EvtLinSample( const char* fname ) | |||||
{ | |||||
std::ifstream IN( fname ); | |||||
float x, y; | |||||
while ( IN >> x >> y ) { | |||||
v.push_back( { x, y } ); | |||||
} | |||||
std::cout << "Point read -- " << v.size() << std::endl; | |||||
I.push_back( 0 ); | |||||
for ( unsigned i = 0; i < v.size() - 1; i++ ) { | |||||
I.push_back( ( v[i + 1].y + v[i + 0].y ) * ( v[i + 1].x - v[i + 0].x ) + | |||||
I.back() ); | |||||
} | |||||
double norm = 1 / I.back(); | |||||
for ( unsigned i = 0; i < v.size(); i++ ) | |||||
I[i] *= norm; | |||||
} | |||||
void EvtLinSample::init( const std::vector<EvtPointf>& _v ) | |||||
{ | |||||
v = _v; | |||||
std::cout << "Point provided -- " << v.size() << std::endl; | |||||
I.push_back( 0 ); | |||||
for ( unsigned i = 0; i < v.size() - 1; i++ ) { | |||||
I.push_back( ( v[i + 1].y + v[i + 0].y ) * ( v[i + 1].x - v[i + 0].x ) + | |||||
I.back() ); | |||||
} | |||||
double norm = 1 / I.back(); | |||||
for ( unsigned i = 0; i < v.size(); i++ ) | |||||
I[i] *= norm; | |||||
} | |||||
std::pair<double, double> EvtLinSample::operator()( double r ) const | |||||
{ | |||||
int j = upper_bound( I.begin(), I.end(), r ) - I.begin(); | |||||
double dI = I[j] - I[j - 1]; | |||||
r = ( r - I[j - 1] ) / dI; | |||||
double f0 = v[j - 1].y, f1 = v[j].y, x0 = v[j - 1].x, x1 = v[j].x, | |||||
df = f1 - f0, dx = x1 - x0, z; | |||||
if ( fabs( df ) > f0 * 1e-3 ) { | |||||
z = ( f1 * x0 - x1 * f0 + dx * sqrt( df * ( f0 + f1 ) * r + f0 * f0 ) ) / | |||||
df; | |||||
} else { | |||||
if ( f0 > f1 ) { | |||||
z = x0 + ( r * dx ) * ( f1 - r * df ) / f0; | |||||
} else { | |||||
z = x1 - ( r * dx ) * ( f0 + r * df ) / f1; | |||||
} | |||||
} | |||||
return std::make_pair( z, f0 + ( df / dx ) * ( z - x0 ) ); | |||||
} |