<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
<html xmlns="http://www.w3.org/1999/xhtml">
<head>
<meta http-equiv="Content-Type" content="text/xhtml;charset=UTF-8"/>
<meta http-equiv="X-UA-Compatible" content="IE=9"/>
<meta name="generator" content="Doxygen 1.8.3.1"/>
<title>LHAPDF C++ wrapper: /examples/CCTest2.cc</title>
<link href="tabs.css" rel="stylesheet" type="text/css"/>
<script type="text/javascript" src="jquery.js"></script>
<script type="text/javascript" src="dynsections.js"></script>
<link href="doxygen.css" rel="stylesheet" type="text/css" />
</head>
<body>
<div id="top"><!-- do not remove this div, it is closed by doxygen! -->
<div id="titlearea">
<table cellspacing="0" cellpadding="0">
<tbody>
<tr style="height: 56px;">
<td style="padding-left: 0.5em;">
<div id="projectname">LHAPDF C++ wrapper
 <span id="projectnumber">5.4</span>
</div>
</td>
</tr>
</tbody>
</table>
</div>
<!-- end header part -->
<!-- Generated by Doxygen 1.8.3.1 -->
<div id="navrow1" class="tabs">
<ul class="tablist">
<li><a href="index.html"><span>Main Page</span></a></li>
<li><a href="pages.html"><span>Related Pages</span></a></li>
<li><a href="namespaces.html"><span>Namespaces</span></a></li>
<li><a href="annotated.html"><span>Classes</span></a></li>
<li><a href="files.html"><span>Files</span></a></li>
<li><a href="examples.html"><span>Examples</span></a></li>
</ul>
</div>
</div><!-- top -->
<div class="header">
<div class="headertitle">
<div class="title">/examples/CCTest2.cc</div> </div>
</div><!--header-->
<div class="contents">
<p>An example of a program using the C++ interface to <a class="el" href="namespaceLHAPDF.html" title="Namespace containing all the LHAPDF wrapper functions.">LHAPDF</a> to calculate PDF errors.</p>
<div class="fragment"><div class="line"></div>
<div class="line"><span class="comment">// Program to demonstrate usage of the MRST 2006 NNLO PDFs. //</span></div>
<div class="line"><span class="comment">// to calculate errors. //</span></div>
<div class="line"><span class="comment"></span><span class="preprocessor">#include "LHAPDF/LHAPDF.h"</span></div>
<div class="line"><span class="preprocessor">#include <iostream></span></div>
<div class="line"><span class="preprocessor">#include <cmath></span></div>
<div class="line"><span class="preprocessor">#include <cstdio></span></div>
<div class="line"><span class="preprocessor">#include <cstdlib></span></div>
<div class="line"><span class="keyword">using namespace </span>std;</div>
<div class="line"></div>
<div class="line"><span class="keyword">using namespace </span>LHAPDF;</div>
<div class="line"></div>
<div class="line"></div>
<div class="line"><span class="keywordtype">double</span> logdist_x(<span class="keywordtype">double</span> xmin, <span class="keywordtype">double</span> xmax, <span class="keywordtype">size_t</span> ix, <span class="keywordtype">size_t</span> nx) {</div>
<div class="line"> <span class="keyword">const</span> <span class="keywordtype">double</span> log10xmin = log10(xmin);</div>
<div class="line"> <span class="keyword">const</span> <span class="keywordtype">double</span> log10xmax = log10(xmax);</div>
<div class="line"> <span class="keyword">const</span> <span class="keywordtype">double</span> log10x = log10xmin + (ix/<span class="keyword">static_cast<</span><span class="keywordtype">double</span><span class="keyword">></span>(nx-1))*(log10xmax-log10xmin);</div>
<div class="line"> <span class="keyword">const</span> <span class="keywordtype">double</span> x = pow(10.0, log10x);</div>
<div class="line"> <span class="keywordflow">return</span> x;</div>
<div class="line">}</div>
<div class="line"></div>
<div class="line"></div>
<div class="line"><span class="keywordtype">int</span> main() {</div>
<div class="line"> <span class="comment">// Show initialisation banner only once</span></div>
<div class="line"> <a class="code" href="namespaceLHAPDF.html#a54d0a3e29ef9a2bf17777da61c99ad3e" title="Choose level of noisiness.">setVerbosity</a>(LOWKEY); <span class="comment">// or SILENT, for no banner at all</span></div>
<div class="line"></div>
<div class="line"> <span class="comment">// You could explicitly set the path to the PDFsets directory</span></div>
<div class="line"> <span class="comment">// setPDFPath("/home/whalley/local/share/lhapdf/PDFsets");</span></div>
<div class="line"> </div>
<div class="line"> <span class="comment">// Initialize PDF sets</span></div>
<div class="line"> <span class="keyword">const</span> <span class="keywordtype">string</span> NAME = <span class="stringliteral">"MRST2006nnlo"</span>;</div>
<div class="line"> initPDFSetM(1, NAME, LHGRID);</div>
<div class="line"> initPDFSetM(2, NAME, LHGRID);</div>
<div class="line"> initPDFSetM(3, NAME, LHGRID);</div>
<div class="line"> </div>
<div class="line"> <span class="comment">// Find the number of eigensets from numberPDF()</span></div>
<div class="line"> <span class="keyword">const</span> <span class="keywordtype">int</span> neigen = numberPDFM(1)/2;</div>
<div class="line"> cout << <span class="stringliteral">"Number of eigensets in this fit = "</span> << neigen << endl;</div>
<div class="line"> <span class="comment">// Find the min and max values of x and Q2 </span></div>
<div class="line"> <span class="keyword">const</span> <span class="keywordtype">double</span> xmin = <a class="code" href="namespaceLHAPDF.html#a4b68147f30e452ed7bbe300630b46a3f" title="Minimum value considered valid for this set, as specified by the set authors.">getXmin</a>(0);</div>
<div class="line"> <span class="keyword">const</span> <span class="keywordtype">double</span> xmax = <a class="code" href="namespaceLHAPDF.html#a9674750111bc08eee46291537f5aca71" title="Maximum value considered valid for this set, as specified by the set authors.">getXmax</a>(0);</div>
<div class="line"> cout << <span class="stringliteral">"Valid x-range = ["</span> << xmin << <span class="stringliteral">", "</span> << xmax << <span class="stringliteral">"]"</span> << endl;</div>
<div class="line"> <span class="comment">// Number of x values to sample</span></div>
<div class="line"> <span class="keyword">const</span> <span class="keywordtype">int</span> nx = 10;</div>
<div class="line"> <span class="comment">// Set the Q scale and flavour</span></div>
<div class="line"> <span class="keywordtype">double</span> q = 10;</div>
<div class="line"> <span class="keywordtype">int</span> flav = 4;</div>
<div class="line"></div>
<div class="line"> <span class="comment">// Get x's and central PDF values</span></div>
<div class="line"> initPDFM(1, 0);</div>
<div class="line"> vector<double> fc(nx), x(nx);</div>
<div class="line"> <span class="keywordflow">for</span> (<span class="keywordtype">int</span> ix = 0; ix < nx; ++ix) {</div>
<div class="line"> x[ix] = logdist_x(xmin, 0.9*xmax, ix, nx);</div>
<div class="line"> fc[ix] = xfxM(1, x[ix], q, flav);</div>
<div class="line"> }</div>
<div class="line"></div>
<div class="line"> <span class="comment">// Sum over error contributions (two ways, depending on how LHDPAF was compiled)</span></div>
<div class="line"> vector<double> summax(nx), summin(nx), sum(nx);</div>
<div class="line"><span class="preprocessor"> #ifndef LHAPDF_LOWMEM</span></div>
<div class="line"><span class="preprocessor"></span> <span class="comment">// This is the normal, efficient, way to do this, with the error</span></div>
<div class="line"> <span class="comment">// sets being initialised the minimum number of times</span></div>
<div class="line"> cout << <span class="stringliteral">"Using efficient set looping"</span> << endl;</div>
<div class="line"> <span class="keywordflow">for</span> (<span class="keywordtype">int</span> ieigen = 1; ieigen <= neigen; ++ieigen) {</div>
<div class="line"> initPDFM(2, 2*ieigen-1);</div>
<div class="line"> initPDFM(3, 2*ieigen);</div>
<div class="line"> <span class="keywordflow">for</span> (<span class="keywordtype">int</span> ix = 0; ix < nx; ++ix) {</div>
<div class="line"> <span class="comment">// Find central and plus/minus values</span></div>
<div class="line"> <span class="keyword">const</span> <span class="keywordtype">double</span> fp = xfxM(2, x[ix], q, flav);</div>
<div class="line"> <span class="keyword">const</span> <span class="keywordtype">double</span> fm = xfxM(3, x[ix], q, flav);</div>
<div class="line"> <span class="comment">// Construct shifts</span></div>
<div class="line"> <span class="keyword">const</span> <span class="keywordtype">double</span> plus = max(max(fp-fc[ix], fm-fc[ix]),0.0);</div>
<div class="line"> <span class="keyword">const</span> <span class="keywordtype">double</span> minus = min(min(fp-fc[ix], fm-fc[ix]),0.0);</div>
<div class="line"> <span class="keyword">const</span> <span class="keywordtype">double</span> diff = fp-fm;</div>
<div class="line"> <span class="comment">// Add it together</span></div>
<div class="line"> summax[ix] += plus*plus;</div>
<div class="line"> summin[ix] += minus*minus;</div>
<div class="line"> sum[ix] += diff*diff;</div>
<div class="line"> }</div>
<div class="line"> }</div>
<div class="line"><span class="preprocessor"> #else</span></div>
<div class="line"><span class="preprocessor"></span> <span class="comment">// In low memory mode, the sets need to be re-initialised with every </span></div>
<div class="line"> <span class="comment">// change of member. Using the approach above gives wrong answers, and</span></div>
<div class="line"> <span class="comment">// reinitialising in all the nested loops is sloooooow! The best way is </span></div>
<div class="line"> <span class="comment">// to calculate the values, plus and minus errors separately.</span></div>
<div class="line"> cout << <span class="stringliteral">"Using low-mem mode set looping"</span> << endl;</div>
<div class="line"> <span class="keywordflow">for</span> (<span class="keywordtype">int</span> ieigen = 1; ieigen <= neigen; ++ieigen) {</div>
<div class="line"> vector<double> fp(nx), fm(nx);</div>
<div class="line"> initPDFM(2, 2*ieigen-1);</div>
<div class="line"> <span class="keywordflow">for</span> (<span class="keywordtype">int</span> ix = 0; ix < nx; ++ix) {</div>
<div class="line"> fp[ix] = xfxM(2, x[ix], q, flav);</div>
<div class="line"> }</div>
<div class="line"> initPDFM(3, 2*ieigen);</div>
<div class="line"> <span class="keywordflow">for</span> (<span class="keywordtype">int</span> ix = 0; ix < nx; ++ix) {</div>
<div class="line"> fm[ix] = xfxM(3, x[ix], q, flav);</div>
<div class="line"> }</div>
<div class="line"> <span class="keywordflow">for</span> (<span class="keywordtype">int</span> ix = 0; ix < nx; ++ix) {</div>
<div class="line"> <span class="comment">// Construct shifts</span></div>
<div class="line"> <span class="keyword">const</span> <span class="keywordtype">double</span> plus = max(max(fp[ix]-fc[ix], fm[ix]-fc[ix]), 0.0);</div>
<div class="line"> <span class="keyword">const</span> <span class="keywordtype">double</span> minus = min(min(fp[ix]-fc[ix], fm[ix]-fc[ix]), 0.0);</div>
<div class="line"> <span class="keyword">const</span> <span class="keywordtype">double</span> diff = fp[ix]-fm[ix];</div>
<div class="line"> <span class="comment">// Add it together</span></div>
<div class="line"> summax[ix] += plus*plus;</div>
<div class="line"> summin[ix] += minus*minus;</div>
<div class="line"> sum[ix] += diff*diff;</div>
<div class="line"> }</div>
<div class="line"> }</div>
<div class="line"><span class="preprocessor"> #endif</span></div>
<div class="line"><span class="preprocessor"></span></div>
<div class="line"> <span class="comment">// Print out results</span></div>
<div class="line"> cout << <span class="stringliteral">"flavour = "</span> << flav << <span class="stringliteral">" Asymmetric (%) Symmetric (%)"</span> << endl;</div>
<div class="line"> cout << <span class="stringliteral">" x Q**2 xf(x) plus minus +- "</span> << endl;</div>
<div class="line"> <span class="keywordflow">for</span> (<span class="keywordtype">int</span> ix = 0; ix < nx; ++ix) {</div>
<div class="line"> printf(<span class="stringliteral">"%0.7f %.0f %10.2E %8.2f %8.2f %8.2f \n"</span>,</div>
<div class="line"> x[ix], q*q, fc[ix], </div>
<div class="line"> sqrt(summax[ix])*100/fc[ix],</div>
<div class="line"> sqrt(summin[ix])*100/fc[ix],</div>
<div class="line"> 0.5*sqrt(sum[ix])*100/fc[ix]);</div>
<div class="line"> }</div>
<div class="line"> </div>
<div class="line"> <span class="keywordflow">return</span> EXIT_SUCCESS;</div>
<div class="line">}</div>
<div class="line"></div>
<div class="line"></div>
<div class="line"></div>
<div class="line"><span class="preprocessor">#include "FortranWrappers.h"</span></div>
<div class="line"><span class="preprocessor">#ifdef FC_DUMMY_MAIN</span></div>
<div class="line"><span class="preprocessor"></span><span class="keywordtype">int</span> FC_DUMMY_MAIN() { <span class="keywordflow">return</span> 1; }</div>
<div class="line"><span class="preprocessor">#endif</span></div>
</div><!-- fragment --> </div><!-- contents -->
<!-- start footer part -->
<hr class="footer"/><address class="footer"><small>
Generated on Thu Nov 21 2013 06:02:29 for LHAPDF C++ wrapper by  <a href="http://www.doxygen.org/index.html">
<img class="footer" src="doxygen.png" alt="doxygen"/>
</a> 1.8.3.1
</small></address>
</body>
</html>
