trento/src/rapidity_profile.h

0001 // TRENTO: Reduced Thickness Event-by-event Nuclear Topology
0002 // Copyright 2015 Jonah E. Bernhard, J. Scott Moreland
0003 // TRENTO3D: Three-dimensional extension of TRENTO by Weiyao Ke
0004 // MIT License
0005
0006 #ifndef ETA_H
0007 #define ETA_H
0008
0009 #include <cmath>
0010 #include <vector>
0011 #include <gsl/gsl_fft_complex.h>
0012
0013 /// GSL fft of complex array, z_i = x_i + j*y_i
0014 /// data are stored in one real array A that A[2i] = x_i and A[2i+1] = y_i
0015 #define REAL(z,i) ((z)[2*(i)])
0016 #define IMAG(z,i) ((z)[2*(i)+1])
0017
0018 double const sqrt2 = std::sqrt(2);
0019 constexpr double TINY = 1e-12;
0020 constexpr int relative_skew_switch = 1;
0021 constexpr int absolute_skew_switch = 2;
0022 /// the mean of the rapidity distribution as function of Ta(x,y),
0023 /// Tb(x,y) and exp(y_beam).
0024 /// For the case both Ta and Tb are tiny small, mean = 0.0 is returned.
0025 /// For the case only Ta(Tb) = 0, it returns +(-)y_beam,
0026 /// For the case y_beam >> 1, mean ~ 0.5*log(Ta/Tb)
0027 double inline mean_function(double ta, double tb, double exp_ybeam){
0028   if (ta < TINY && tb < TINY) return 0.;
0029   return 0.5 * std::log((ta*exp_ybeam + tb/exp_ybeam) / std::max(ta/exp_ybeam + tb*exp_ybeam, TINY));
0030 }
0031
0032 /// The coefficient of normalized width of the rapidity distribution
0033 /// Currently, it is simply a constant independent of Ta(x,y) and Tb(x,y)
0034 double inline std_function(double ta, double tb){
0035   (void)ta;
0036   (void)tb;
0037   return 1.;
0038 }
0039
0040 /// The normalized skewness as function of Ta(x,y) and Tb(x,y)
0041
0042 double inline skew_function(double ta, double tb, int type_switch){
0043   if (type_switch == relative_skew_switch)
0044     return (ta - tb)/std::max(ta + tb, TINY);
0045   else if(type_switch == absolute_skew_switch)
0046       return ta-tb;
0047   else return 0.;
0048 }
0049
0050 /// A fast pseudorapidity to rapidity transformer using pretabulated values
0051 /// It returns the transfoamtion y(eta) and the jacobian dy/deta(eta)
0052 class fast_eta2y {
0053 private:
0054   double etamax_;
0055   double deta_;
0056   std::size_t neta_;
0057   std::vector<double> y_;
0058   std::vector<double> dydeta_;
0059 public:
0060   fast_eta2y(double J, double etamax, double deta)
0061       : etamax_(etamax),
0062         deta_(deta),
0063         neta_(std::ceil(2.*etamax_/(deta_+1e-15))+1),
0064         y_(neta_, 0.),
0065         dydeta_(neta_, 0.) {
0066
0067     for (std::size_t ieta = 0; ieta < neta_; ++ieta) {
0068       double eta = -etamax_ + ieta*deta_;
0069       double Jsh = J*std::sinh(eta);
0070       double sq = std::sqrt(1. + Jsh*Jsh);
0071       y_[ieta] = std::log(sq + Jsh);
0072       dydeta_[ieta] = J*std::cosh(eta)/sq;
0073     }
0074   }
0075
0076   double rapidity(double eta){
0077     double steps = (eta + etamax_)/deta_;
0078     double xi = std::fmod(steps, 1.);
0079     std::size_t index = std::floor(steps);
0080     return y_[index]*(1. - xi) + y_[index+1]*xi;
0081   }
0082
0083   double Jacobian(double eta){
0084     double steps = (eta + etamax_)/deta_;
0085     double xi = std::fmod(steps, 1.);
0086     std::size_t index = std::floor(steps);
0087     return dydeta_[index]*(1. - xi) + dydeta_[index+1]*xi;
0088   }
0089
0090 };
0091
0092
0093 /// A class for cumulant generating function inversion
0094 /// This class handles inverse fourier transform the cumulant generating function
0095 /// A direct transformation F^{-1} exp(i*m*k-(std*k)^2/2-skew*(std*k)^3/6) results
0096 /// in an ill-behaved function. Instead, skew is replaced by skew*exp(-(std*k)^2/2).
0097 /// In this way higher order cumulant are included to regulated the function.
0098 /// The reconstruction range is taken to be [-3.33, 3.33]*std, which does not
0099 /// seem to be very large, however, by trail and error, this range elimiates
0100 /// oscillations at large rapidity and leaves the details close to mid-rapidity
0101 /// unchanged. We used GSL FFT library (more specialized library would be FFTW,
0102 /// e.g.), with 256 points. The transformed results are stored and interpolated.
0103 class cumulant_generating{
0104 private:
0105   size_t const N;
0106   double * data, * dsdy;
0107   double eta_max;
0108   double deta;
0109   double center;
0110
0111 public:
0112   cumulant_generating(): N(256), data(new double[2*N]), dsdy(new double[2*N]){};
0113
0114   /// This function set the mean, std and skew of the profile and use FFT to
0115   /// transform cumulant generating function at zero mean.
0116   void calculate_dsdy(double mean, double std, double skew){
0117     double k1, k2, k3, amp, arg;
0118     // adaptive eta_max = 3.33*std;
0119     center = mean;
0120     eta_max = std*3.33;
0121     deta = 2.*eta_max/(N-1.);
0122     double fftmean = eta_max/std;
0123       for(size_t i=0;i<N;i++){
0124           k1 = M_PI*(i-N/2.0)/eta_max*std;
0125       k2 = k1*k1;
0126       k3 = k2*k1;
0127
0128           amp = std::exp(-k2/2.0);
0129           arg = fftmean*k1+skew/6.0*k3*amp;
0130
0131       REAL(data,i) = amp*std::cos(arg);
0132           IMAG(data,i) = amp*std::sin(arg);
0133       }
0134        gsl_fft_complex_radix2_forward(data, 1, N);
0135
0136       for(size_t i=0;i<N;i++){
0137           dsdy[i] = REAL(data,i)*(2.0*static_cast<double>(i%2 == 0)-1.0);
0138       }
0139   }
0140
0141   /// When interpolating the funtion, the mean is put back by simply shifting
0142   /// the function by y = y - mean + dy/2, the last term is correcting for
0143   /// interpolating bin edge instead of bin center
0144   double interp_dsdy(double y){
0145     y = y-center+deta/2.;
0146     if (y < -eta_max || y >= eta_max) return 0.0;
0147     double xy = (y+eta_max)/deta;
0148     size_t iy = std::floor(xy);
0149     double ry = xy-iy;
0150     return dsdy[iy]*(1.-ry) + dsdy[iy+1]*ry;
0151   }
0152 };
0153 #endif