DoubleCrystalBallGenerator.cc

Go to the documentation of this file.

//FAMOS headers
#include "FastSimulation/Utilities/interface/DoubleCrystalBallGenerator.h"
#include "FastSimulation/Utilities/interface/RandomEngineAndDistribution.h"

//ROOT headers
#include "TMath.h"
#include <iostream>

using namespace TMath;

double DoubleCrystalBallGenerator::shoot(double mu, double sigma, double aL, double nL, double aR, double nR,
                                         RandomEngineAndDistribution const* random){
    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;
}