---

AlbedoPSpectrum.cxx

Go to the documentation of this file.

// $Id: AlbedoPSpectrum.cxx,v 1.6 2002/06/24 22:37:15 srobinsn Exp $


#include "AlbedoPSpectrum.h"

#include <cmath>
#include <algorithm>
#include <functional>
#include "CLHEP/Random/Random.h"
#include "GPS.h"
#include "Geomag.h"
#include "SpectrumFactory.h"

static SpectrumFactory<AlbedoPSpectrum> factory;
const ISpectrumFactory& AlbedoPSpectrumFactory = factory;


void AlbedoPSpectrum::init(const std::vector<float>& params) {
    
    float lat =  params.size()>0? params[0]: 0.0f;
    float lon =  params.size()>1? params[1]: 0.0f;
    setPosition(lat, lon);
    
    m_particle_name = "p";
    
    // set callback to be notified when the position changes
    m_observer.setAdapter( new ActionAdapter<AlbedoPSpectrum>
        (this, &AlbedoPSpectrum::askGPS) );
    
    GPS::instance()->notification().attach( &m_observer );
    
}


//-------------------------- constructors---------------------

AlbedoPSpectrum::AlbedoPSpectrum(const std::string& paramstring) {
    std::vector<float> params;
    
    parseParamList(paramstring,params);
    
    init(params);
}


//------------------------ title()

std::string AlbedoPSpectrum::title() const {
    return "AlbedoPSpectrum";
}

//------------------------- particleName

const char * AlbedoPSpectrum::particleName() const {
    return m_particle_name.c_str();
}

//-------------------------- setParticleName

void AlbedoPSpectrum::setParticleName(std::string name)
{
    m_particle_name = name;
}

//-------------------------- flux() (current position)

double AlbedoPSpectrum::flux(double) const {
    return m_flux;
}

//--------------------------

double AlbedoPSpectrum::solidAngle() const
{
    return 4.*M_PI;
}

//-------------------------- flux() (specified position)

float AlbedoPSpectrum::flux(float lat, float lon) const {
    double v1,v2, alf1, alf2, emin, emax, Ejoin;
    fitParams(lat, lon, alf1, alf2, emin, emax, v1, v2, Ejoin);
    // The factor 1.47 is induced by the assumed angular dependence.
    return 1.47 * (v1+v2);
}

//-------------------------- flux() (position given in a pair)

float AlbedoPSpectrum::flux( std::pair<double,double> coords) const {
    return flux(coords.first, coords.second);
}

//-------------------------- operator()  sample an energy value

float AlbedoPSpectrum::operator() (float x) const{
    double v1,v2, alf1, alf2, emin, emax, Ejoin;
    fitParams(m_lat, m_lon, alf1, alf2, emin, emax, v1, v2, Ejoin);
    double split = v1/(v1+v2);
    if (x > split) {
        return pow((pow(Ejoin,1.-alf2)+((1.-x)/(1.-split))*(pow(emax,1.-alf2)
            -pow(Ejoin,1.-alf2))), 1./(1.-alf2));
    } else {
        return pow((pow(emin,1.-alf1)+(x/split)*(pow(Ejoin,1.-alf1)
            -pow(emin,1.-alf1))), 1./(1.-alf1));
    }
}

//-------------------------- askGPS()

int AlbedoPSpectrum::askGPS()
{
    setPosition(GPS::instance()->lat(), GPS::instance()->lon());
    return 0; // can't be void in observer pattern
}

//-------------------------- setPosition(separate coordinates)
void AlbedoPSpectrum::setPosition(float lat, float lon) {
    m_lat = lat;
    m_lon = lon;
    m_flux = flux(lat, lon);
}

//-------------------------- setPosition (coordinates packaged in a pair)
void AlbedoPSpectrum::setPosition( std::pair<double,double> coords) {
    AlbedoPSpectrum::setPosition(coords.first, coords.second);
}

//------------------------- calculate_rate()
double AlbedoPSpectrum::calculate_rate(double old_rate)
{
    return flux(GPS::instance()->lat(), GPS::instance()->lon());
}

//-------------------------- dir()
std::pair<float,float> AlbedoPSpectrum::dir(float energy)const
{
    // Random particle direction
    
    float coszenith, earthazi,dens,v;
    const int try_max = 1000;
    static int max_tried = 0;
    int trial = 0;
    earthazi = 2. * M_PI * HepRandom::getTheGenerator()->flat();
    do {
        coszenith = 2. * HepRandom::getTheGenerator()->flat() - 1.;
        dens = 1. + 0.6 * sqrt(1.-coszenith*coszenith);
        v = 1.6 * HepRandom::getTheGenerator()->flat();
    } while (v > dens && trial++ < try_max);
    
    max_tried = std::max(trial, max_tried);
    
    return std::make_pair<float,float>(coszenith, earthazi);
}

//------------------------- fitParams()
void AlbedoPSpectrum::fitParams(const double lat, const double lon, 
                                double& alf1, double& alf2, double& emin, double& emax, 
                                double& v1, double &v2, double& Ejoin) const
{
    double theta = abs(Geomag::geolat(lat,lon)) * M_PI / 180.;
    double e1 = .2;  // pivot point of lower power law
    double e2 = 1.;  // pivot point of higher power law
    emin = .01;      
    emax = 10.;
    // normalization values of the two power laws
    double a1 = 0.142 + theta * (-.4809 + theta * .5517);
    double a2 = .0499 + theta * (-.3239 + theta * (.8077 + theta * (-.8800 + 
        theta * .3607)));
    alf1 = .4913 + theta * (2.017 + theta * (-.6941 - 1.49 * theta));
    alf2 = 2.85 - .875 * theta;
    Ejoin = pow(a1*pow(e1,alf1) / (a2*pow(e2,alf2)), 1./(alf1-alf2));
    v1 = 1000.*a1*pow(e1,alf1)*(pow(Ejoin,1.-alf1)-pow(emin,1.-alf1))
        /(1.-alf1);  // integral of lower power law
    v2 = 1000.*a2*pow(e2,alf2)*(pow(emax,1.-alf2)-pow(Ejoin,1.-alf2)) 
        /(1.-alf2);  // integral of upper power law
}

---