DoubleCrystalBallGenerator Class Reference

#include <DoubleCrystalBallGenerator.h>

List of all members.

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

---

Detailed Description

Definition at line 12 of file DoubleCrystalBallGenerator.h.

---

Constructor & Destructor Documentation

DoubleCrystalBallGenerator::DoubleCrystalBallGenerator

(

const RandomEngine *

engine

)

[inline]

Definition at line 16 of file DoubleCrystalBallGenerator.h.

: random(engine) {}

virtual DoubleCrystalBallGenerator::~DoubleCrystalBallGenerator

(

)

[inline, virtual]

Definition at line 18 of file DoubleCrystalBallGenerator.h.

{}

---

Member Function Documentation

double DoubleCrystalBallGenerator::shoot	(	double	mu,
		double	sigma,
		double	aL,
		double	nL,
		double	aR,
		double	nR
	)

Definition at line 11 of file DoubleCrystalBallGenerator.cc.

References funct::D, PFRecoTauDiscriminationAgainstElectronDeadECAL_cfi::dR, N, Pi, random, x, and detailsBasic3DVector::y.

Referenced by HCALResponse::cballShootNoNegative().

                                                                                                           {
        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;
        
}

---

Member Data Documentation

const RandomEngine* DoubleCrystalBallGenerator::random [private]

Definition at line 25 of file DoubleCrystalBallGenerator.h.