DoubleCrystalBallGenerator Class Reference

#include <DoubleCrystalBallGenerator.h>

Public Member Functions
	DoubleCrystalBallGenerator ()
double	shoot (double mu, double sigma, double aL, double nL, double aR, double nR, RandomEngineAndDistribution const *random)
virtual	~DoubleCrystalBallGenerator ()

Detailed Description

Definition at line 12 of file DoubleCrystalBallGenerator.h.

Constructor & Destructor Documentation

DoubleCrystalBallGenerator::DoubleCrystalBallGenerator ( )

inline

Definition at line 16 of file DoubleCrystalBallGenerator.h.

{ }

virtual DoubleCrystalBallGenerator::~DoubleCrystalBallGenerator ( )

inlinevirtual

Definition at line 18 of file DoubleCrystalBallGenerator.h.

References RPCpg::mu, random, and shoot().

{}

Member Function Documentation

double DoubleCrystalBallGenerator::shoot	(	double	mu,
		double	sigma,
		double	aL,
		double	nL,
		double	aR,
		double	nR,
		RandomEngineAndDistribution const *	random
	)

Definition at line 11 of file DoubleCrystalBallGenerator.cc.

References RandomEngineAndDistribution::flatShoot(), N, Pi, x, and y.

Referenced by ~DoubleCrystalBallGenerator().

                                                                                   {
    if(nL<=1 || nR<=1) return 0; //n>1 required

    double dL = nL/aL;
    double dR = nR/aR;
    double N = 1/(sigma*(nL/aL*1/(nL-1)*Exp(-aL*aL/2) + Sqrt(Pi()/2)*(Erf(aL/Sqrt(2))+Erf(aR/Sqrt(2))) + nR/aR*1/(nR-1)*Exp(-aR*aR/2)));

    //start x as NaN
    double x = 0./0.;
    while(!Finite(x)){
        //shoot a flat random number
        double y = random->flatShoot();

        //check range
        //crystal ball CDF changes at (x-mu)/sigma = -aL && (x-mu)/sigma = aR
        //compute this in pieces because it is awful
        double AL = dL/(nL-1)*Exp(-aL*aL/2);
        double CL = Sqrt(Pi()/2)*Erf(aL/Sqrt(2));
        double CR = Sqrt(Pi()/2)*Erf(aR/Sqrt(2));
        if(y < sigma*N*AL){//below lower boundary, use inverted power law CDF (left side)
            double BL = dL/(nL-1)*Exp(-aL*aL/2);
            x = mu + sigma*(-dL*Power(y/(sigma*N*BL),1/(-nL+1)) - aL + dL);
        }
        else if(y > sigma*N*(AL+CL+CR)){//above lower boundary, use inverted power law CDF (right side)
            double AR = dR/(nR-1)*Exp(-aR*aR/2);
            double BR = dR/(1-nR)*Exp(-aR*aR/2);
            double D = (y/(sigma*N)-AL-CL-CR-AR)/BR;
            x = mu + sigma*(dR*Power(D,1/(-nR+1)) + aR - dR);
        }
        else{//between boundaries, use gaussian CDF with proper normalization (in terms of erfc)
            double D = 1 - Sqrt(2/Pi())*(y/(sigma*N)-AL-CL);
            x = mu + sigma*Sqrt(2)*ErfcInverse(D);
        }
    }
    return x;
}