|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
|
//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;
}
|