Utilities/src/DoubleCrystalBallGenerator.cc

0001 //FAMOS headers
0002 #include "FastSimulation/Utilities/interface/DoubleCrystalBallGenerator.h"
0003 #include "FastSimulation/Utilities/interface/RandomEngineAndDistribution.h"
0004
0005 //ROOT headers
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;  //n>1 required
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   //start x as NaN
0023   double x = 0. / 0.;
0024   while (!Finite(x)) {
0025     //shoot a flat random number
0026     double y = random->flatShoot();
0027
0028     //check range
0029     //crystal ball CDF changes at (x-mu)/sigma = -aL && (x-mu)/sigma = aR
0030     //compute this in pieces because it is awful
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) {  //below lower boundary, use inverted power law CDF (left side)
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)) {  //above lower boundary, use inverted power law CDF (right side)
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 {  //between boundaries, use gaussian CDF with proper normalization (in terms of erfc)
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 }