File indexing completed on 2023-10-25 09:46:13
0001
0002 #include "FastSimulation/Utilities/interface/DoubleCrystalBallGenerator.h"
0003 #include "FastSimulation/Utilities/interface/RandomEngineAndDistribution.h"
0004
0005
0006 #include "TMath.h"
0007 #include <iostream>
0008
0009 using namespace TMath;
0010
0011 double DoubleCrystalBallGenerator::shoot(
0012 double mu, double sigma, double aL, double nL, double aR, double nR, RandomEngineAndDistribution const* random) {
0013 if (nL <= 1 || nR <= 1)
0014 return 0;
0015
0016 double dL = nL / aL;
0017 double dR = nR / aR;
0018 double N = 1 / (sigma * (nL / aL * 1 / (nL - 1) * Exp(-aL * aL / 2) +
0019 Sqrt(Pi() / 2) * (Erf(aL / Sqrt(2)) + Erf(aR / Sqrt(2))) +
0020 nR / aR * 1 / (nR - 1) * Exp(-aR * aR / 2)));
0021
0022
0023 double x = 0. / 0.;
0024 while (!Finite(x)) {
0025
0026 double y = random->flatShoot();
0027
0028
0029
0030
0031 double AL = dL / (nL - 1) * Exp(-aL * aL / 2);
0032 double CL = Sqrt(Pi() / 2) * Erf(aL / Sqrt(2));
0033 double CR = Sqrt(Pi() / 2) * Erf(aR / Sqrt(2));
0034 if (y < sigma * N * AL) {
0035 double BL = dL / (nL - 1) * Exp(-aL * aL / 2);
0036 x = mu + sigma * (-dL * Power(y / (sigma * N * BL), 1 / (-nL + 1)) - aL + dL);
0037 } else if (y > sigma * N * (AL + CL + CR)) {
0038 double AR = dR / (nR - 1) * Exp(-aR * aR / 2);
0039 double BR = dR / (1 - nR) * Exp(-aR * aR / 2);
0040 double D = (y / (sigma * N) - AL - CL - CR - AR) / BR;
0041 x = mu + sigma * (dR * Power(D, 1 / (-nR + 1)) + aR - dR);
0042 } else {
0043 double D = 1 - Sqrt(2 / Pi()) * (y / (sigma * N) - AL - CL);
0044 x = mu + sigma * Sqrt(2) * ErfcInverse(D);
0045 }
0046 }
0047 return x;
0048 }