Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
RooLandau.cxx
Go to the documentation of this file.
1 /*****************************************************************************
2  * Project: RooFit *
3  * Package: RooFitModels *
4  * @(#)root/roofit:$Id$
5  * Authors: *
6  * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu *
7  * DK, David Kirkby, UC Irvine, dkirkby@uci.edu *
8  * *
9  * Copyright (c) 2000-2005, Regents of the University of California *
10  * and Stanford University. All rights reserved. *
11  * *
12  * Redistribution and use in source and binary forms, *
13  * with or without modification, are permitted according to the terms *
14  * listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
15  *****************************************************************************/
16 
17 /** \class RooLandau
18  \ingroup Roofit
19 
20 Landau distribution p.d.f
21 \image html RF_Landau.png "PDF of the Landau distribution."
22 **/
23 
24 #include "RooLandau.h"
25 #include "RooHelpers.h"
26 #include "RooFit.h"
27 #include "RooRandom.h"
28 #include "BatchHelpers.h"
29 
30 #include "TMath.h"
31 
32 ClassImp(RooLandau);
33 
34 ////////////////////////////////////////////////////////////////////////////////
35 
36 RooLandau::RooLandau(const char *name, const char *title, RooAbsReal& _x, RooAbsReal& _mean, RooAbsReal& _sigma) :
37  RooAbsPdf(name,title),
38  x("x","Dependent",this,_x),
39  mean("mean","Mean",this,_mean),
40  sigma("sigma","Width",this,_sigma)
41 {
42  RooHelpers::checkRangeOfParameters(this, {&_sigma}, 0.0);
43 }
44 
45 ////////////////////////////////////////////////////////////////////////////////
46 
47 RooLandau::RooLandau(const RooLandau& other, const char* name) :
48  RooAbsPdf(other,name),
49  x("x",this,other.x),
50  mean("mean",this,other.mean),
51  sigma("sigma",this,other.sigma)
52 {
53 }
54 
55 ////////////////////////////////////////////////////////////////////////////////
56 
57 Double_t RooLandau::evaluate() const
58 {
59  return TMath::Landau(x, mean, sigma);
60 }
61 
62 ////////////////////////////////////////////////////////////////////////////////
63 
64 namespace {
65 
66 //Author: Emmanouil Michalainas, CERN 26 JULY 2019
67 /* Actual computation of Landau(x,mean,sigma) in a vectorization-friendly way
68  * Code copied from function landau_pdf (math/mathcore/src/PdfFuncMathCore.cxx)
69  * and rewritten to take advantage for the most popular case
70  * which is -1 < (x-mean)/sigma < 1. The rest cases are handled in scalar way
71  */
72 
73 constexpr double p1[5] = {0.4259894875,-0.1249762550, 0.03984243700, -0.006298287635, 0.001511162253};
74 constexpr double q1[5] = {1.0 ,-0.3388260629, 0.09594393323, -0.01608042283, 0.003778942063};
75 
76 constexpr double p2[5] = {0.1788541609, 0.1173957403, 0.01488850518, -0.001394989411, 0.0001283617211};
77 constexpr double q2[5] = {1.0 , 0.7428795082, 0.3153932961, 0.06694219548, 0.008790609714};
78 
79 constexpr double p3[5] = {0.1788544503, 0.09359161662,0.006325387654, 0.00006611667319,-0.000002031049101};
80 constexpr double q3[5] = {1.0 , 0.6097809921, 0.2560616665, 0.04746722384, 0.006957301675};
81 
82 constexpr double p4[5] = {0.9874054407, 118.6723273, 849.2794360, -743.7792444, 427.0262186};
83 constexpr double q4[5] = {1.0 , 106.8615961, 337.6496214, 2016.712389, 1597.063511};
84 
85 constexpr double p5[5] = {1.003675074, 167.5702434, 4789.711289, 21217.86767, -22324.94910};
86 constexpr double q5[5] = {1.0 , 156.9424537, 3745.310488, 9834.698876, 66924.28357};
87 
88 constexpr double p6[5] = {1.000827619, 664.9143136, 62972.92665, 475554.6998, -5743609.109};
89 constexpr double q6[5] = {1.0 , 651.4101098, 56974.73333, 165917.4725, -2815759.939};
90 
91 constexpr double a1[3] = {0.04166666667,-0.01996527778, 0.02709538966};
92 constexpr double a2[2] = {-1.845568670,-4.284640743};
93 
94 template<class Tx, class Tmean, class Tsigma>
95 void compute( size_t batchSize,
96  double * __restrict output,
97  Tx X, Tmean M, Tsigma S)
98 {
99  const double NaN = std::nan("");
100  constexpr size_t block=256;
101  double v[block];
102 
103  for (size_t i=0; i<batchSize; i+=block) { //CHECK_VECTORISE
104  const size_t stop = (i+block < batchSize) ? block : batchSize-i ;
105 
106  for (size_t j=0; j<stop; j++) { //CHECK_VECTORISE
107  v[j] = (X[i+j]-M[i+j]) / S[i+j];
108  output[i+j] = (p2[0]+(p2[1]+(p2[2]+(p2[3]+p2[4]*v[j])*v[j])*v[j])*v[j]) /
109  (q2[0]+(q2[1]+(q2[2]+(q2[3]+q2[4]*v[j])*v[j])*v[j])*v[j]);
110  }
111 
112  for (size_t j=0; j<stop; j++) { //CHECK_VECTORISE
113  const bool mask = S[i+j] > 0;
114  /* comparison with NaN will give result false, so the next
115  * loop won't affect output, for cases where sigma <=0
116  */
117  if (!mask) v[j] = NaN;
118  output[i+j] *= mask;
119  }
120 
121  double u, ue, us;
122  for (size_t j=0; j<stop; j++) { //CHECK_VECTORISE
123  // if branch written in way to quickly process the most popular case -1 <= v[j] < 1
124  if (v[j] >= 1) {
125  if (v[j] < 5) {
126  output[i+j] = (p3[0]+(p3[1]+(p3[2]+(p3[3]+p3[4]*v[j])*v[j])*v[j])*v[j]) /
127  (q3[0]+(q3[1]+(q3[2]+(q3[3]+q3[4]*v[j])*v[j])*v[j])*v[j]);
128  } else if (v[j] < 12) {
129  u = 1/v[j];
130  output[i+j] = u*u*(p4[0]+(p4[1]+(p4[2]+(p4[3]+p4[4]*u)*u)*u)*u) /
131  (q4[0]+(q4[1]+(q4[2]+(q4[3]+q4[4]*u)*u)*u)*u);
132  } else if (v[j] < 50) {
133  u = 1/v[j];
134  output[i+j] = u*u*(p5[0]+(p5[1]+(p5[2]+(p5[3]+p5[4]*u)*u)*u)*u) /
135  (q5[0]+(q5[1]+(q5[2]+(q5[3]+q5[4]*u)*u)*u)*u);
136  } else if (v[j] < 300) {
137  u = 1/v[j];
138  output[i+j] = u*u*(p6[0]+(p6[1]+(p6[2]+(p6[3]+p6[4]*u)*u)*u)*u) /
139  (q6[0]+(q6[1]+(q6[2]+(q6[3]+q6[4]*u)*u)*u)*u);
140  } else {
141  u = 1 / (v[j] -v[j]*std::log(v[j])/(v[j]+1) );
142  output[i+j] = u*u*(1 +(a2[0] +a2[1]*u)*u );
143  }
144  } else if (v[j] < -1) {
145  if (v[j] >= -5.5) {
146  u = std::exp(-v[j]-1);
147  output[i+j] = std::exp(-u)*std::sqrt(u)*
148  (p1[0]+(p1[1]+(p1[2]+(p1[3]+p1[4]*v[j])*v[j])*v[j])*v[j])/
149  (q1[0]+(q1[1]+(q1[2]+(q1[3]+q1[4]*v[j])*v[j])*v[j])*v[j]);
150  } else {
151  u = std::exp(v[j]+1.0);
152  if (u < 1e-10) output[i+j] = 0.0;
153  else {
154  ue = std::exp(-1/u);
155  us = std::sqrt(u);
156  output[i+j] = 0.3989422803*(ue/us)*(1+(a1[0]+(a1[1]+a1[2]*u)*u)*u);
157  }
158  }
159  }
160  }
161  }
162 }
163 };
164 
165 ////////////////////////////////////////////////////////////////////////////////
166 
167 /// Compute \f$ Landau(x,mean,sigma) \f$ in batches.
168 /// The local proxies {x, mean, sigma} will be searched for batch input data,
169 /// and if found, the computation will be batched over their
170 /// values. If batch data are not found for one of the proxies, the proxies value is assumed to
171 /// be constant over the batch.
172 /// \param[in] batchIndex Index of the batch to be computed.
173 /// \param[in] batchSize Size of each batch. The last batch may be smaller.
174 /// \return A span with the computed values.
175 
176 RooSpan<double> RooLandau::evaluateBatch(std::size_t begin, std::size_t batchSize) const {
177  using namespace BatchHelpers;
178  auto xData = x.getValBatch(begin, batchSize);
179  auto meanData = mean.getValBatch(begin, batchSize);
180  auto sigmaData = sigma.getValBatch(begin, batchSize);
181  const bool batchX = !xData.empty();
182  const bool batchMean = !meanData.empty();
183  const bool batchSigma = !sigmaData.empty();
184 
185  if (!batchX && !batchMean && !batchSigma) {
186  return {};
187  }
188  batchSize = findSize({ xData, meanData, sigmaData });
189  auto output = _batchData.makeWritableBatchUnInit(begin, batchSize);
190 
191  if (batchX && !batchMean && !batchSigma ) {
192  compute(batchSize, output.data(), xData, BracketAdapter<double>(mean), BracketAdapter<double>(sigma));
193  }
194  else if (!batchX && batchMean && !batchSigma ) {
195  compute(batchSize, output.data(), BracketAdapter<double>(x), meanData, BracketAdapter<double>(sigma));
196  }
197  else if (batchX && batchMean && !batchSigma ) {
198  compute(batchSize, output.data(), xData, meanData, BracketAdapter<double>(sigma));
199  }
200  else if (!batchX && !batchMean && batchSigma ) {
201  compute(batchSize, output.data(), BracketAdapter<double>(x), BracketAdapter<double>(mean), sigmaData);
202  }
203  else if (batchX && !batchMean && batchSigma ) {
204  compute(batchSize, output.data(), xData, BracketAdapter<double>(mean), sigmaData);
205  }
206  else if (!batchX && batchMean && batchSigma ) {
207  compute(batchSize, output.data(), BracketAdapter<double>(x), meanData, sigmaData);
208  }
209  else if (batchX && batchMean && batchSigma ) {
210  compute(batchSize, output.data(), xData, meanData, sigmaData);
211  }
212  return output;
213 }
214 
215 ////////////////////////////////////////////////////////////////////////////////
216 
217 Int_t RooLandau::getGenerator(const RooArgSet& directVars, RooArgSet &generateVars, Bool_t /*staticInitOK*/) const
218 {
219  if (matchArgs(directVars,generateVars,x)) return 1 ;
220  return 0 ;
221 }
222 
223 ////////////////////////////////////////////////////////////////////////////////
224 
225 void RooLandau::generateEvent(Int_t code)
226 {
227  assert(1 == code); (void)code;
228  Double_t xgen ;
229  while(1) {
230  xgen = RooRandom::randomGenerator()->Landau(mean,sigma);
231  if (xgen<x.max() && xgen>x.min()) {
232  x = xgen ;
233  break;
234  }
235  }
236  return;
237 }