// #include <CLHEP/Units/SystemOfUnits.h>
#include <G4Electron.hh>
#include <G4Event.hh>
#include <G4EventManager.hh>
#include <G4Material.hh>
#include <G4MaterialPropertiesTable.hh>
#include <G4OpticalPhoton.hh>
#include <G4PhysicalConstants.hh>
#include <G4VPhysicalVolume.hh>
#include <G4PrimaryVertex.hh>
#include <G4RunManager.hh>
#include <G4StateManager.hh>
#include <RAT/AuxSimInfo.hh>
#include <RAT/Cerenkov.hh>
#include <RAT/ChannelEfficiency.hh>
#include <RAT/GLG4Scint.hh>
#include <RAT/HitTimeProbabilityFactor.hh>
#include <RAT/Log.hh>
#include <RAT/ModeCut.hh>
#include <RAT/Noise.hh>
#include <RAT/Optics.hh>
#include <RAT/Optimiser.hh>
#include <RAT/PEnergy.hh>
#include <RAT/PhotonThinning.hh>
#include <RAT/PMTSelector.hh>
#include <RAT/PMTVariation.hh>
#include <RAT/Profiler.hh>
#include <RAT/Processor.hh>
#include <RAT/RunManager.hh>
#include <RAT/SphereCrossPoint.hh>
#include <RAT/Trajectory.hh>
#include <RAT/SteppingAction.hh>
#include <RAT/TrackingAction.hh>
#include <RAT/DB.hh>
#include <RAT/DS/MC.hh>
#include <RAT/DS/MCPMT.hh>
#include <RAT/DS/MCTrack.hh>
#include <RAT/DU/ChanHWStatus.hh>
#include <RAT/DU/PMTInfo.hh>
#include <RAT/DU/Utility.hh>
#include <cmath>

using std::vector;
using std::string;

using CLHEP::second;
using CLHEP::electronvolt;
using CLHEP::nanometer;

extern RAT::RunManager* gTheRatRunManager;
const double radToDeg = 180./pi;

namespace RAT
{
namespace Methods
{

// PMT reponse table
const int nPol = 2;             // number of polarization entries
const int nWL = 61;             // number of wavelength entries in table
const int nAng = 46;            // number of angle entries in table
const double lbWL = 200.;       // lower bound in wavelength
const double lbAng = 0.;        // lower bound in angle
const double delWL = 10.;       // step in wavelength (in nm)
const double delAng = 2.;       // step in angle (in degrees)

// Following are probabilities that a photon incident on a PMT produces a
// photoelectron, as a function of photon polarization, wavelength (in nm), and
// angle of incidence (in degrees).
const float pmtRsp[nPol][nWL][nAng] =
{
{
{ 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07},
{ 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07},
{ 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07},
{ 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07},
{ 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07},
{0.0040643,0.0040000,0.0039904,0.0040162,0.0039695,0.0039348,0.0039994,0.0039659,0.0039363,0.0040126,0.0039842,0.0039555,0.0039269,0.0038166,0.0038486,0.0038143,0.0037465,0.0035957,0.0035421,0.0034831,0.0033657,0.0033508,0.0032224,0.0030984,0.0029095,0.0027135,0.0025157,0.0021687,0.0018128,0.0014285,0.0010424,0.0007923,0.0006047,0.0004407,0.0002907,0.0002002,0.0001868,0.0001727,0.0001585,0.0001629,0.0001414,0.0001372,0.0001182,0.0001020,0.0000803, 1.49e-06},
{0.0088993,0.0088769,0.0088215,0.0088526,0.0087977,0.0088271,0.0088417,0.0088256,0.0088843,0.0088439,0.0088603,0.0088058,0.0087269,0.0086266,0.0084958,0.0084273,0.0081921,0.0080667,0.0079171,0.0077050,0.0075853,0.0073608,0.0070665,0.0068126,0.0064797,0.0060829,0.0055433,0.0048901,0.0040070,0.0031014,0.0023508,0.0017896,0.0013527,0.0009666,0.0006641,0.0004569,0.0004305,0.0003849,0.0003695,0.0003401,0.0003393,0.0002895,0.0002624,0.0002262,0.0001819, 3.41e-06},
{0.0232539,0.0232104,0.0231934,0.0230260,0.0231105,0.0231930,0.0233340,0.0232839,0.0232554,0.0231900,0.0230643,0.0231144,0.0229401,0.0227282,0.0223834,0.0219265,0.0215348,0.0212307,0.0205676,0.0202637,0.0197674,0.0193284,0.0187069,0.0180469,0.0171389,0.0159915,0.0145171,0.0126315,0.0104360,0.0083548,0.0061832,0.0046958,0.0035193,0.0025531,0.0017067,0.0011800,0.0011376,0.0009879,0.0009761,0.0009174,0.0008488,0.0007753,0.0006773,0.0006051,0.0004520, 6.39e-06},
{0.0482784,0.0484824,0.0482083,0.0480897,0.0482124,0.0483274,0.0483659,0.0483065,0.0485312,0.0483395,0.0484175,0.0480273,0.0477436,0.0472901,0.0466053,0.0458641,0.0447866,0.0442179,0.0429994,0.0418993,0.0409282,0.0400622,0.0388328,0.0374818,0.0357184,0.0334209,0.0302737,0.0266605,0.0220454,0.0173555,0.0130711,0.0097919,0.0072976,0.0053301,0.0035442,0.0024343,0.0023465,0.0021423,0.0019810,0.0019105,0.0017333,0.0016183,0.0014413,0.0012460,0.0009103,0.0000124},
{0.0760496,0.0758471,0.0756364,0.0756865,0.0755108,0.0758603,0.0762245,0.0761293,0.0761778,0.0762375,0.0759827,0.0757114,0.0751276,0.0744671,0.0732886,0.0723495,0.0707523,0.0693625,0.0678665,0.0661864,0.0648850,0.0630694,0.0612413,0.0590869,0.0562852,0.0526895,0.0479908,0.0419676,0.0349412,0.0271789,0.0203996,0.0156113,0.0115693,0.0083274,0.0055274,0.0039086,0.0036000,0.0034151,0.0031744,0.0029845,0.0027872,0.0025540,0.0022824,0.0019225,0.0014556,0.0000136},
{0.1012799,0.1011061,0.1009353,0.1008616,0.1008520,0.1009091,0.1010222,0.1010409,0.1012735,0.1017372,0.1010473,0.1006452,0.1003387,0.0993344,0.0980328,0.0964309,0.0943329,0.0924173,0.0903725,0.0881864,0.0863376,0.0840435,0.0816015,0.0786752,0.0750790,0.0703785,0.0638963,0.0559964,0.0466305,0.0363727,0.0274119,0.0207879,0.0155406,0.0112763,0.0074728,0.0052505,0.0048207,0.0045987,0.0042475,0.0039691,0.0037286,0.0033721,0.0029855,0.0025551,0.0019240,0.0000185},
{0.1248988,0.1248288,0.1247625,0.1247201,0.1249832,0.1252058,0.1255738,0.1254952,0.1260302,0.1261634,0.1261512,0.1256860,0.1255508,0.1247137,0.1233265,0.1221318,0.1200356,0.1178158,0.1156294,0.1132441,0.1112618,0.1086496,0.1063806,0.1037774,0.1003717,0.0950388,0.0866658,0.0766613,0.0631504,0.0484909,0.0356326,0.0261487,0.0195408,0.0141917,0.0095063,0.0067159,0.0061446,0.0057378,0.0053936,0.0051717,0.0047376,0.0042652,0.0036707,0.0031301,0.0023093,0.0000222},
{0.1443542,0.1442780,0.1441039,0.1447977,0.1452711,0.1456239,0.1458703,0.1461091,0.1463951,0.1475555,0.1477050,0.1478611,0.1478967,0.1476532,0.1470824,0.1456106,0.1438507,0.1418141,0.1391883,0.1376396,0.1363849,0.1342739,0.1323323,0.1306354,0.1276426,0.1225781,0.1127683,0.1000650,0.0819091,0.0614660,0.0434688,0.0311382,0.0228922,0.0166087,0.0113977,0.0081080,0.0075715,0.0069919,0.0066032,0.0061408,0.0056801,0.0048948,0.0042490,0.0035614,0.0026711,0.0000279},
{0.1526571,0.1528128,0.1527659,0.1533666,0.1538767,0.1548394,0.1553497,0.1561534,0.1571061,0.1584209,0.1584654,0.1593938,0.1595583,0.1599606,0.1592479,0.1582707,0.1567131,0.1547242,0.1522812,0.1512417,0.1497227,0.1481729,0.1462696,0.1458283,0.1430948,0.1374573,0.1272298,0.1126019,0.0917917,0.0682626,0.0476049,0.0335772,0.0245267,0.0179911,0.0123906,0.0089252,0.0082481,0.0076741,0.0072153,0.0067461,0.0061487,0.0052679,0.0045350,0.0037798,0.0027563,0.0000215},
{0.1608426,0.1607495,0.1612765,0.1622671,0.1630253,0.1636532,0.1644053,0.1657099,0.1661174,0.1675578,0.1686021,0.1697178,0.1700976,0.1708368,0.1705088,0.1696639,0.1677966,0.1658662,0.1640922,0.1627108,0.1617299,0.1598443,0.1592392,0.1585695,0.1562577,0.1506019,0.1398630,0.1239255,0.1009757,0.0746965,0.0515580,0.0360098,0.0261076,0.0193416,0.0132882,0.0095987,0.0088748,0.0082473,0.0077233,0.0072569,0.0066501,0.0056654,0.0047672,0.0039548,0.0029216,0.0000273},
{0.1629406,0.1632332,0.1634383,0.1645974,0.1651629,0.1666011,0.1676812,0.1684253,0.1699237,0.1715597,0.1726102,0.1738948,0.1749495,0.1754873,0.1754690,0.1748432,0.1731225,0.1708735,0.1692415,0.1683527,0.1670288,0.1660532,0.1651930,0.1652556,0.1635301,0.1578922,0.1468835,0.1303691,0.1063401,0.0782268,0.0536437,0.0371770,0.0266907,0.0198552,0.0137261,0.0100337,0.0092767,0.0085378,0.0080262,0.0074897,0.0069003,0.0057783,0.0048253,0.0040481,0.0030309,0.0000211},
{0.1660799,0.1663133,0.1665755,0.1674892,0.1691985,0.1699259,0.1712466,0.1725414,0.1737987,0.1759951,0.1774083,0.1785186,0.1798068,0.1804275,0.1808064,0.1797787,0.1789108,0.1769148,0.1750100,0.1736462,0.1730986,0.1720609,0.1718694,0.1721525,0.1704294,0.1652292,0.1538756,0.1366343,0.1113627,0.0817908,0.0558818,0.0384714,0.0276179,0.0206030,0.0141257,0.0103323,0.0096420,0.0088818,0.0083186,0.0078377,0.0070679,0.0059144,0.0049453,0.0040997,0.0030530,0.0000215},
{0.1679685,0.1677785,0.1685949,0.1696831,0.1712354,0.1723321,0.1735821,0.1747776,0.1769482,0.1787998,0.1801122,0.1816145,0.1833295,0.1838026,0.1840398,0.1836577,0.1824852,0.1803787,0.1785587,0.1774322,0.1769095,0.1761689,0.1759629,0.1765570,0.1754347,0.1703335,0.1588886,0.1415896,0.1151925,0.0844652,0.0573517,0.0393512,0.0283600,0.0209595,0.0144860,0.0106976,0.0097949,0.0091208,0.0085510,0.0079734,0.0072814,0.0060794,0.0050837,0.0042875,0.0030958,0.0000181},
{0.1680281,0.1687069,0.1686194,0.1699685,0.1714858,0.1726179,0.1739395,0.1755342,0.1772509,0.1790075,0.1807876,0.1824705,0.1835182,0.1846933,0.1845510,0.1843046,0.1830230,0.1812552,0.1789137,0.1784032,0.1776554,0.1771587,0.1763706,0.1776217,0.1762596,0.1712153,0.1598946,0.1423885,0.1160204,0.0850552,0.0576415,0.0398530,0.0285423,0.0211957,0.0145623,0.0107474,0.0098620,0.0091945,0.0085457,0.0079502,0.0073103,0.0060941,0.0051342,0.0042479,0.0030952,0.0000200},
{0.1698055,0.1700848,0.1704251,0.1719267,0.1730000,0.1740863,0.1756245,0.1768543,0.1785252,0.1809237,0.1824999,0.1843461,0.1854124,0.1865674,0.1869071,0.1865127,0.1845974,0.1827120,0.1807468,0.1792931,0.1790937,0.1782200,0.1774826,0.1787923,0.1774809,0.1724360,0.1612752,0.1440637,0.1170173,0.0857163,0.0582501,0.0403201,0.0288554,0.0212854,0.0148541,0.0108492,0.0098458,0.0092409,0.0086087,0.0081329,0.0073514,0.0061651,0.0051331,0.0042354,0.0031908,0.0000179},
{0.1722801,0.1724890,0.1726058,0.1739118,0.1756673,0.1767778,0.1779574,0.1788143,0.1804931,0.1828149,0.1844262,0.1857702,0.1875372,0.1883542,0.1885406,0.1884775,0.1862993,0.1843321,0.1817521,0.1802703,0.1795263,0.1781883,0.1769583,0.1781589,0.1767171,0.1719874,0.1609029,0.1438905,0.1171525,0.0864053,0.0583184,0.0400812,0.0286890,0.0215291,0.0148812,0.0108833,0.0100733,0.0092558,0.0087233,0.0080880,0.0074656,0.0061380,0.0051612,0.0042882,0.0031293,0.0000153},
{0.1745478,0.1749802,0.1753289,0.1767653,0.1780113,0.1789069,0.1802809,0.1811346,0.1826567,0.1843512,0.1864707,0.1878520,0.1895557,0.1904377,0.1904769,0.1899736,0.1880285,0.1857186,0.1831678,0.1812087,0.1801174,0.1781806,0.1772618,0.1771376,0.1763908,0.1717565,0.1606892,0.1439683,0.1176019,0.0864068,0.0583259,0.0402217,0.0288658,0.0215840,0.0150351,0.0109853,0.0099783,0.0093018,0.0087774,0.0082247,0.0074809,0.0061991,0.0051521,0.0043282,0.0031178,0.0000149},
{0.1757819,0.1757067,0.1758204,0.1775246,0.1784415,0.1795919,0.1807129,0.1811876,0.1830009,0.1847499,0.1864364,0.1880562,0.1892882,0.1899344,0.1908036,0.1901578,0.1883678,0.1858179,0.1830441,0.1805416,0.1791399,0.1771329,0.1754279,0.1754319,0.1744362,0.1700571,0.1592673,0.1425634,0.1168456,0.0856055,0.0582100,0.0398285,0.0285925,0.0214015,0.0149276,0.0109523,0.0099825,0.0093027,0.0087186,0.0081474,0.0073919,0.0062247,0.0051142,0.0043201,0.0031073,0.0000124},
{0.1757069,0.1757291,0.1761299,0.1774034,0.1787994,0.1796567,0.1808997,0.1818032,0.1827331,0.1846933,0.1859568,0.1875830,0.1890094,0.1901493,0.1904126,0.1900677,0.1881157,0.1852875,0.1821472,0.1796650,0.1778153,0.1756990,0.1733472,0.1731301,0.1721404,0.1678897,0.1573182,0.1408775,0.1155570,0.0849968,0.0574850,0.0394995,0.0284993,0.0212475,0.0147994,0.0109321,0.0100302,0.0092760,0.0086927,0.0081003,0.0073161,0.0061470,0.0051378,0.0042528,0.0031314,0.0000117},
{0.1732812,0.1731796,0.1735472,0.1750085,0.1761455,0.1773056,0.1779261,0.1787770,0.1797949,0.1817063,0.1830364,0.1845646,0.1858647,0.1871802,0.1876641,0.1871757,0.1849977,0.1822066,0.1787088,0.1751482,0.1741248,0.1713621,0.1685188,0.1687154,0.1674839,0.1634984,0.1531182,0.1373593,0.1125370,0.0830153,0.0560878,0.0386513,0.0279425,0.0209646,0.0145904,0.0106581,0.0098317,0.0091093,0.0085821,0.0079137,0.0072675,0.0060209,0.0050939,0.0041101,0.0031237,0.0000126},
{0.1686971,0.1687059,0.1691787,0.1703269,0.1716709,0.1724021,0.1732588,0.1740420,0.1749991,0.1770239,0.1783129,0.1795804,0.1813427,0.1820343,0.1827681,0.1818788,0.1797544,0.1767180,0.1731216,0.1703044,0.1681495,0.1651468,0.1625636,0.1620292,0.1612731,0.1573796,0.1474854,0.1322375,0.1087700,0.0800688,0.0540669,0.0371862,0.0269463,0.0201715,0.0141789,0.0103165,0.0095480,0.0087968,0.0082696,0.0076754,0.0070460,0.0058275,0.0048675,0.0040492,0.0030149, 7.24e-06},
{0.1639214,0.1634982,0.1639419,0.1654109,0.1667018,0.1676807,0.1683450,0.1691810,0.1697636,0.1717001,0.1727793,0.1740801,0.1752556,0.1764292,0.1769725,0.1761557,0.1744028,0.1712168,0.1678079,0.1649014,0.1625806,0.1595433,0.1563936,0.1558149,0.1549378,0.1511708,0.1414479,0.1271979,0.1048961,0.0771257,0.0521395,0.0360132,0.0261091,0.0195830,0.0137406,0.0100175,0.0092354,0.0085461,0.0079951,0.0074562,0.0068213,0.0056909,0.0046820,0.0039523,0.0028782, 6.18e-06},
{0.1547595,0.1551186,0.1552460,0.1571131,0.1578856,0.1591027,0.1596890,0.1600364,0.1610373,0.1626978,0.1636294,0.1652575,0.1665063,0.1675152,0.1678339,0.1669118,0.1652241,0.1624626,0.1581778,0.1555348,0.1530735,0.1499208,0.1471312,0.1460415,0.1448716,0.1413563,0.1326185,0.1194671,0.0983076,0.0728298,0.0491461,0.0339817,0.0245325,0.0185787,0.0129587,0.0095278,0.0087297,0.0080890,0.0075700,0.0071142,0.0064432,0.0053325,0.0045282,0.0037167,0.0027301, 6.39e-06},
{0.1466692,0.1464294,0.1470109,0.1484190,0.1494530,0.1500241,0.1511900,0.1514599,0.1523506,0.1538449,0.1551048,0.1563365,0.1573178,0.1587866,0.1587770,0.1583410,0.1560443,0.1529529,0.1496100,0.1460875,0.1439499,0.1411031,0.1376277,0.1368560,0.1354758,0.1326503,0.1240741,0.1118322,0.0921082,0.0682575,0.0459808,0.0318403,0.0231559,0.0175774,0.0122682,0.0089359,0.0081485,0.0076403,0.0071734,0.0066658,0.0060535,0.0050288,0.0041876,0.0034969,0.0025789, 4.90e-06},
{0.1399651,0.1397945,0.1399438,0.1416703,0.1423042,0.1436292,0.1443644,0.1448215,0.1453642,0.1466072,0.1480085,0.1493003,0.1505186,0.1512475,0.1519570,0.1512777,0.1490858,0.1462809,0.1427112,0.1393429,0.1370004,0.1339012,0.1302403,0.1294447,0.1283672,0.1254899,0.1176362,0.1058856,0.0874954,0.0646950,0.0437255,0.0302884,0.0219676,0.0165478,0.0115853,0.0086375,0.0077983,0.0072850,0.0068756,0.0063634,0.0057610,0.0048183,0.0039883,0.0033255,0.0024049, 4.47e-06},
{0.1328613,0.1330673,0.1330660,0.1344302,0.1357649,0.1367448,0.1376422,0.1378763,0.1381957,0.1393587,0.1409804,0.1422765,0.1433793,0.1442539,0.1449845,0.1438251,0.1419568,0.1390699,0.1354298,0.1324867,0.1302403,0.1267894,0.1231110,0.1220400,0.1212948,0.1182430,0.1110317,0.0998226,0.0829606,0.0613188,0.0412545,0.0287955,0.0208398,0.0158200,0.0112109,0.0081879,0.0074224,0.0069327,0.0065489,0.0060522,0.0055555,0.0046641,0.0038479,0.0031348,0.0023010, 4.69e-06},
{0.1256013,0.1256773,0.1258599,0.1272543,0.1279868,0.1292032,0.1302920,0.1302875,0.1307979,0.1323050,0.1333393,0.1345932,0.1355108,0.1364317,0.1368417,0.1364417,0.1341727,0.1315832,0.1279468,0.1245344,0.1223827,0.1193461,0.1159657,0.1145987,0.1135031,0.1109738,0.1039898,0.0939740,0.0778899,0.0577095,0.0389261,0.0270637,0.0196102,0.0149879,0.0105065,0.0077082,0.0070236,0.0065538,0.0061387,0.0057987,0.0052665,0.0043764,0.0035732,0.0030053,0.0021681, 4.47e-06},
{0.1098300,0.1099340,0.1102375,0.1113323,0.1122884,0.1127764,0.1136995,0.1138000,0.1144786,0.1157806,0.1164407,0.1176835,0.1188656,0.1194596,0.1199674,0.1192956,0.1171561,0.1147989,0.1116679,0.1088300,0.1068313,0.1036658,0.1005906,0.1000155,0.0986618,0.0962251,0.0903446,0.0818070,0.0677476,0.0498158,0.0339563,0.0234805,0.0170633,0.0129578,0.0091661,0.0067700,0.0061353,0.0057367,0.0053561,0.0049934,0.0046209,0.0037399,0.0031674,0.0026362,0.0019080, 1.49e-06},
{0.0868403,0.0869936,0.0872155,0.0882230,0.0890092,0.0895953,0.0903176,0.0905619,0.0908021,0.0917401,0.0930437,0.0935063,0.0945261,0.0950469,0.0954117,0.0948373,0.0932256,0.0911849,0.0884181,0.0862011,0.0844241,0.0819859,0.0794535,0.0784302,0.0775576,0.0756379,0.0709393,0.0641237,0.0533042,0.0396002,0.0268356,0.0185928,0.0134969,0.0103182,0.0072677,0.0053116,0.0049197,0.0045772,0.0042886,0.0039838,0.0036153,0.0029655,0.0025086,0.0020863,0.0015022, 1.70e-06},
{0.0684820,0.0688799,0.0686330,0.0697749,0.0702166,0.0706756,0.0710884,0.0714662,0.0717940,0.0726952,0.0734488,0.0739250,0.0746892,0.0750599,0.0750422,0.0747429,0.0736505,0.0717340,0.0696652,0.0678452,0.0665199,0.0641629,0.0620831,0.0614620,0.0608656,0.0593105,0.0553199,0.0503012,0.0416760,0.0310390,0.0210675,0.0146319,0.0106618,0.0081540,0.0057521,0.0042767,0.0038662,0.0036000,0.0033646,0.0031668,0.0028567,0.0023783,0.0019870,0.0016609,0.0012104, 8.52e-07},
{0.0613621,0.0610492,0.0610509,0.0621870,0.0628965,0.0630592,0.0633683,0.0638128,0.0640004,0.0647478,0.0655091,0.0659052,0.0669568,0.0669365,0.0673397,0.0668703,0.0653655,0.0641080,0.0621478,0.0604927,0.0589623,0.0571576,0.0552204,0.0544848,0.0537800,0.0525853,0.0491169,0.0445849,0.0370526,0.0276799,0.0188413,0.0131355,0.0095188,0.0072824,0.0051623,0.0037457,0.0034756,0.0032192,0.0030179,0.0028505,0.0025644,0.0021265,0.0017898,0.0014750,0.0010596, 1.06e-06},
{0.0491736,0.0491453,0.0493778,0.0498302,0.0503419,0.0506113,0.0511774,0.0513753,0.0516599,0.0521689,0.0527289,0.0533112,0.0537278,0.0541329,0.0542650,0.0538483,0.0528897,0.0516115,0.0500407,0.0485406,0.0473295,0.0458354,0.0440733,0.0436873,0.0428918,0.0422034,0.0393133,0.0357836,0.0296763,0.0222592,0.0150164,0.0105740,0.0076993,0.0058520,0.0041683,0.0030132,0.0028245,0.0025785,0.0024066,0.0022944,0.0020645,0.0016997,0.0014166,0.0011619,0.0008654, 8.52e-07},
{0.0419504,0.0423989,0.0425401,0.0429947,0.0433666,0.0436524,0.0439199,0.0440019,0.0445046,0.0450897,0.0453442,0.0458933,0.0463698,0.0464869,0.0465163,0.0462988,0.0454848,0.0442999,0.0430750,0.0419808,0.0407078,0.0394984,0.0380611,0.0374136,0.0367749,0.0359625,0.0336637,0.0306826,0.0254888,0.0190581,0.0129310,0.0090330,0.0066213,0.0050607,0.0035559,0.0026645,0.0024006,0.0022328,0.0021144,0.0019376,0.0017687,0.0015033,0.0012415,0.0010371,0.0007565, 6.39e-07},
{0.0355708,0.0355033,0.0357923,0.0358707,0.0364920,0.0367678,0.0369114,0.0370383,0.0373312,0.0377157,0.0383231,0.0383625,0.0390620,0.0390852,0.0391231,0.0387793,0.0382818,0.0371879,0.0360809,0.0350015,0.0341512,0.0329945,0.0317408,0.0312124,0.0306194,0.0300577,0.0281308,0.0256530,0.0212437,0.0160724,0.0108362,0.0076136,0.0055774,0.0042807,0.0030515,0.0022475,0.0020049,0.0019246,0.0017793,0.0017016,0.0014956,0.0012782,0.0010624,0.0008400,0.0006552, 8.52e-07},
{0.0283900,0.0286398,0.0285306,0.0289715,0.0292226,0.0293866,0.0295037,0.0297604,0.0299463,0.0304345,0.0306277,0.0308818,0.0312458,0.0314850,0.0314611,0.0311154,0.0307059,0.0300021,0.0289312,0.0280677,0.0272592,0.0263815,0.0254788,0.0249972,0.0247212,0.0240285,0.0223131,0.0204011,0.0171073,0.0128351,0.0087099,0.0061150,0.0044635,0.0034290,0.0024760,0.0017712,0.0016273,0.0014999,0.0014360,0.0013406,0.0012190,0.0010185,0.0008398,0.0007095,0.0005174, 4.26e-07},
{0.0217804,0.0218347,0.0218599,0.0222405,0.0224192,0.0225122,0.0226995,0.0228260,0.0229401,0.0232739,0.0234258,0.0238030,0.0239378,0.0240452,0.0240241,0.0239519,0.0234635,0.0230869,0.0223001,0.0215702,0.0209610,0.0202581,0.0195161,0.0190660,0.0187625,0.0182390,0.0171225,0.0156413,0.0130539,0.0097642,0.0067157,0.0046886,0.0034411,0.0026569,0.0019323,0.0013991,0.0013033,0.0011632,0.0011039,0.0010051,0.0009346,0.0007676,0.0006588,0.0005406,0.0003798, 1.00e-07},
{0.0153957,0.0154424,0.0154128,0.0157470,0.0158803,0.0159157,0.0160307,0.0161723,0.0162072,0.0163534,0.0165974,0.0167706,0.0169474,0.0169663,0.0169604,0.0168922,0.0166963,0.0160873,0.0156941,0.0151834,0.0147923,0.0142356,0.0135534,0.0135414,0.0133527,0.0128818,0.0121144,0.0110720,0.0091876,0.0069048,0.0047830,0.0033455,0.0024490,0.0019105,0.0013719,0.0009602,0.0009112,0.0008483,0.0007902,0.0007387,0.0006824,0.0005555,0.0004630,0.0003808,0.0002982, 2.13e-07},
{0.0111453,0.0110560,0.0110518,0.0112652,0.0113883,0.0113802,0.0114703,0.0116545,0.0116509,0.0117678,0.0118984,0.0120662,0.0121397,0.0122579,0.0122665,0.0121911,0.0118539,0.0115495,0.0113169,0.0109632,0.0106550,0.0102066,0.0098790,0.0096771,0.0094462,0.0091766,0.0086128,0.0079289,0.0065693,0.0049838,0.0034771,0.0024075,0.0017836,0.0013606,0.0009625,0.0007284,0.0006656,0.0006013,0.0005883,0.0005406,0.0004846,0.0004179,0.0003361,0.0002895,0.0002166, 2.13e-07},
{0.0072209,0.0071738,0.0072558,0.0073842,0.0073167,0.0074552,0.0074569,0.0074882,0.0075404,0.0076379,0.0077499,0.0077791,0.0078417,0.0079655,0.0079112,0.0079365,0.0077983,0.0075517,0.0073904,0.0070705,0.0069223,0.0066605,0.0063700,0.0062360,0.0061169,0.0060036,0.0056138,0.0051450,0.0043237,0.0032756,0.0022302,0.0015693,0.0011593,0.0009059,0.0006445,0.0004630,0.0004319,0.0003989,0.0003683,0.0003632,0.0003265,0.0002552,0.0002147,0.0001840,0.0001318, 2.13e-07},
{0.0044837,0.0044575,0.0045014,0.0045446,0.0045889,0.0046273,0.0046464,0.0046918,0.0047544,0.0047661,0.0047915,0.0048707,0.0048797,0.0048965,0.0049363,0.0048750,0.0047144,0.0046816,0.0045600,0.0044330,0.0042852,0.0041069,0.0039389,0.0038726,0.0038173,0.0037195,0.0034867,0.0031894,0.0026679,0.0020079,0.0013930,0.0009942,0.0007184,0.0005476,0.0004100,0.0002967,0.0002624,0.0002330,0.0002230,0.0002209,0.0002051,0.0001608,0.0001363,0.0001163,0.0000846, 1.00e-07},
{0.0024400,0.0024243,0.0023951,0.0025001,0.0025208,0.0025018,0.0025563,0.0025874,0.0025632,0.0025674,0.0026528,0.0026106,0.0026726,0.0027112,0.0026863,0.0027131,0.0026839,0.0025600,0.0024532,0.0024051,0.0023450,0.0022253,0.0021391,0.0021295,0.0020863,0.0020422,0.0018939,0.0017318,0.0014371,0.0011295,0.0007653,0.0005438,0.0004124,0.0003024,0.0002106,0.0001595,0.0001387,0.0001331,0.0001359,0.0001203,0.0001080,0.0000950,0.0000699,0.0000626,0.0000443, 1.00e-07},
{0.0013020,0.0012939,0.0012995,0.0013078,0.0013295,0.0013316,0.0013585,0.0013691,0.0013527,0.0013834,0.0013996,0.0014160,0.0014298,0.0014584,0.0014064,0.0014396,0.0013911,0.0013444,0.0013527,0.0012728,0.0012285,0.0012028,0.0011657,0.0011551,0.0011165,0.0010952,0.0010249,0.0009235,0.0007994,0.0006049,0.0003919,0.0002803,0.0002162,0.0001685,0.0001161,0.0000918,0.0000835,0.0000735,0.0000756,0.0000630,0.0000624,0.0000441,0.0000449,0.0000330,0.0000279, 1.00e-07},
{0.0006155,0.0006479,0.0006251,0.0006456,0.0006867,0.0006532,0.0006622,0.0006669,0.0006641,0.0006633,0.0006703,0.0006914,0.0007229,0.0007095,0.0007171,0.0007001,0.0006769,0.0006556,0.0006543,0.0006305,0.0006055,0.0005868,0.0005653,0.0005393,0.0005482,0.0005327,0.0004916,0.0004641,0.0003981,0.0002931,0.0002030,0.0001361,0.0001027,0.0000871,0.0000584,0.0000454,0.0000407,0.0000343,0.0000296,0.0000364,0.0000277,0.0000239,0.0000194,0.0000145,0.0000145, 1.00e-07},
{0.0003180,0.0003235,0.0003076,0.0003342,0.0003301,0.0003174,0.0003395,0.0003304,0.0003335,0.0003576,0.0003431,0.0003574,0.0003712,0.0003476,0.0003542,0.0003606,0.0003410,0.0003465,0.0003384,0.0003091,0.0003007,0.0002837,0.0002824,0.0002809,0.0002654,0.0002599,0.0002575,0.0002428,0.0001938,0.0001393,0.0001065,0.0000773,0.0000588,0.0000403,0.0000264,0.0000241,0.0000202,0.0000181,0.0000175,0.0000149,0.0000153,0.0000119,0.0000104, 8.95e-06, 5.54e-06, 1.00e-07},
{0.0001361,0.0001346,0.0001297,0.0001431,0.0001465,0.0001399,0.0001440,0.0001384,0.0001382,0.0001444,0.0001421,0.0001499,0.0001495,0.0001453,0.0001538,0.0001495,0.0001529,0.0001350,0.0001404,0.0001314,0.0001342,0.0001229,0.0001240,0.0001206,0.0001135,0.0001103,0.0001137,0.0000892,0.0000824,0.0000637,0.0000445,0.0000332,0.0000224,0.0000187,0.0000102,0.0000109, 8.73e-06, 6.60e-06, 9.58e-06, 6.60e-06, 7.03e-06, 5.11e-06, 3.19e-06, 4.69e-06, 3.62e-06, 1.00e-07},
{0.0000765,0.0000786,0.0000784,0.0000869,0.0000794,0.0000858,0.0000835,0.0000848,0.0000860,0.0000841,0.0000890,0.0000888,0.0000890,0.0000852,0.0000826,0.0000848,0.0000743,0.0000816,0.0000812,0.0000748,0.0000720,0.0000758,0.0000767,0.0000707,0.0000628,0.0000588,0.0000562,0.0000599,0.0000479,0.0000324,0.0000209,0.0000185,0.0000164, 9.37e-06, 8.31e-06, 5.54e-06, 4.90e-06, 4.26e-06, 3.19e-06, 4.05e-06, 4.26e-06, 4.26e-06, 2.34e-06, 2.98e-06, 1.28e-06, 1.00e-07},
{0.0000665,0.0000671,0.0000667,0.0000654,0.0000671,0.0000660,0.0000603,0.0000656,0.0000630,0.0000669,0.0000703,0.0000688,0.0000724,0.0000756,0.0000665,0.0000630,0.0000675,0.0000707,0.0000601,0.0000611,0.0000596,0.0000567,0.0000564,0.0000558,0.0000488,0.0000526,0.0000532,0.0000490,0.0000379,0.0000300,0.0000213,0.0000138, 9.80e-06, 7.67e-06, 5.32e-06, 4.26e-06, 4.90e-06, 2.13e-06, 2.98e-06, 2.56e-06, 2.56e-06, 1.70e-06, 2.56e-06, 2.34e-06, 6.39e-07, 1.00e-07},
{0.0000498,0.0000452,0.0000530,0.0000460,0.0000513,0.0000490,0.0000575,0.0000545,0.0000524,0.0000530,0.0000560,0.0000558,0.0000518,0.0000569,0.0000498,0.0000560,0.0000522,0.0000515,0.0000490,0.0000494,0.0000494,0.0000447,0.0000426,0.0000413,0.0000428,0.0000417,0.0000381,0.0000426,0.0000324,0.0000277,0.0000149, 9.58e-06, 6.60e-06, 7.24e-06, 4.69e-06, 1.92e-06, 1.92e-06, 2.98e-06, 2.77e-06, 2.13e-06, 1.92e-06, 1.28e-06, 1.92e-06, 1.28e-06, 1.49e-06, 1.00e-07},
{0.0000409,0.0000435,0.0000360,0.0000392,0.0000366,0.0000394,0.0000411,0.0000396,0.0000407,0.0000407,0.0000464,0.0000424,0.0000449,0.0000458,0.0000415,0.0000392,0.0000432,0.0000354,0.0000413,0.0000328,0.0000411,0.0000375,0.0000351,0.0000341,0.0000343,0.0000341,0.0000302,0.0000298,0.0000222,0.0000175, 9.58e-06, 8.95e-06, 5.96e-06, 4.90e-06, 5.75e-06, 1.70e-06, 2.56e-06, 2.56e-06, 2.98e-06, 2.77e-06, 6.39e-07, 1.06e-06, 1.28e-06, 4.26e-07, 4.26e-07, 1.00e-07},
{0.0000281,0.0000305,0.0000288,0.0000334,0.0000292,0.0000390,0.0000305,0.0000313,0.0000271,0.0000337,0.0000300,0.0000328,0.0000351,0.0000328,0.0000273,0.0000330,0.0000339,0.0000292,0.0000330,0.0000305,0.0000281,0.0000309,0.0000260,0.0000226,0.0000256,0.0000234,0.0000234,0.0000190,0.0000200,0.0000134, 9.37e-06, 7.03e-06, 3.83e-06, 4.47e-06, 2.98e-06, 1.28e-06, 3.62e-06, 1.92e-06, 2.56e-06, 2.13e-06, 6.39e-07, 1.49e-06, 1.28e-06, 8.52e-07, 6.39e-07, 1.00e-07},
{0.0000194,0.0000249,0.0000230,0.0000228,0.0000200,0.0000239,0.0000209,0.0000268,0.0000196,0.0000224,0.0000211,0.0000264,0.0000300,0.0000258,0.0000253,0.0000249,0.0000236,0.0000258,0.0000247,0.0000236,0.0000251,0.0000219,0.0000187,0.0000170,0.0000164,0.0000204,0.0000211,0.0000168,0.0000141,0.0000132, 7.88e-06, 5.96e-06, 3.19e-06, 5.11e-06, 1.92e-06, 1.70e-06, 2.34e-06, 1.06e-06, 1.28e-06, 4.26e-07, 6.39e-07, 1.06e-06, 4.26e-07, 2.13e-07, 8.52e-07, 1.00e-07},
{0.0000194,0.0000211,0.0000160,0.0000153,0.0000160,0.0000219,0.0000217,0.0000196,0.0000151,0.0000173,0.0000183,0.0000187,0.0000211,0.0000183,0.0000194,0.0000179,0.0000183,0.0000190,0.0000158,0.0000192,0.0000162,0.0000145,0.0000207,0.0000160,0.0000115,0.0000134,0.0000138,0.0000141, 7.45e-06, 9.16e-06, 5.11e-06, 2.77e-06, 3.62e-06, 3.41e-06, 2.13e-06, 8.52e-07, 1.49e-06, 1.49e-06, 6.39e-07, 8.52e-07, 1.49e-06, 1.06e-06, 4.26e-07, 6.39e-07, 6.39e-07, 1.00e-07},
{0.0000109,0.0000121,0.0000121,0.0000121,0.0000132,0.0000138,0.0000155,0.0000155,0.0000126,0.0000143,0.0000147,0.0000153,0.0000132,0.0000175,0.0000147,0.0000109,0.0000132,0.0000121,0.0000151,0.0000104,0.0000138,0.0000104,0.0000100, 9.37e-06,0.0000121,0.0000138,0.0000117, 9.37e-06, 6.39e-06, 7.24e-06, 4.47e-06, 3.83e-06, 1.92e-06, 1.92e-06, 2.13e-06, 1.49e-06, 1.70e-06, 4.26e-07, 4.26e-07, 2.13e-07, 6.39e-07, 4.26e-07, 4.26e-07, 4.26e-07, 4.26e-07, 1.00e-07},
{ 7.88e-06,0.0000104, 8.73e-06,0.0000104, 7.88e-06, 8.52e-06, 9.37e-06, 7.03e-06, 9.16e-06,0.0000102,0.0000136,0.0000102,0.0000100,0.0000117, 6.82e-06, 7.67e-06,0.0000113, 9.58e-06, 6.60e-06, 8.09e-06, 9.58e-06, 5.32e-06, 7.24e-06, 5.54e-06, 7.24e-06, 5.96e-06, 4.47e-06, 7.24e-06, 4.26e-06, 3.19e-06, 3.41e-06, 2.34e-06, 1.49e-06, 1.70e-06, 6.39e-07, 1.06e-06, 6.39e-07, 6.39e-07, 2.13e-07, 1.06e-06, 4.26e-07, 2.13e-07, 6.39e-07, 1.06e-07, 1.06e-07, 1.00e-07},
{ 2.56e-06, 3.41e-06, 5.11e-06, 4.69e-06, 4.47e-06, 3.41e-06, 4.26e-06, 4.26e-06, 5.96e-06, 4.90e-06, 4.47e-06, 5.11e-06, 4.90e-06, 5.11e-06, 3.41e-06, 4.05e-06, 4.69e-06, 3.83e-06, 3.62e-06, 3.19e-06, 3.19e-06, 2.98e-06, 4.05e-06, 2.77e-06, 2.98e-06, 4.05e-06, 3.62e-06, 3.62e-06, 2.56e-06, 1.70e-06, 1.06e-06, 8.52e-07, 8.52e-07, 1.06e-06, 2.13e-07, 6.39e-07, 2.13e-07, 3.19e-07, 3.19e-07, 1.06e-07, 1.06e-07, 1.06e-07, 1.06e-07, 2.13e-07, 1.00e-07, 1.00e-07},
{ 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07}
},

{
{ 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07},
{ 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07},
{ 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07},
{ 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07},
{ 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07},
{0.0040652,0.0039629,0.0039321,0.0039606,0.0039589,0.0040494,0.0040077,0.0040292,0.0039840,0.0040673,0.0039623,0.0039459,0.0039365,0.0038356,0.0038138,0.0037702,0.0036062,0.0035785,0.0035071,0.0034051,0.0032590,0.0031700,0.0029555,0.0028422,0.0026337,0.0024190,0.0021442,0.0018509,0.0014982,0.0011968,0.0008995,0.0006922,0.0005061,0.0003789,0.0002833,0.0002230,0.0001979,0.0001985,0.0001791,0.0001502,0.0001623,0.0001340,0.0001308,0.0001042,0.0000754, 3.62e-06},
{0.0088903,0.0088532,0.0088814,0.0088400,0.0088635,0.0088181,0.0088852,0.0089610,0.0088547,0.0088673,0.0088890,0.0087985,0.0087817,0.0086072,0.0084835,0.0083161,0.0081152,0.0080228,0.0077529,0.0075318,0.0072873,0.0069983,0.0066939,0.0063478,0.0059078,0.0053186,0.0047695,0.0040916,0.0033338,0.0026590,0.0019751,0.0015308,0.0011212,0.0008132,0.0005987,0.0004622,0.0004251,0.0004183,0.0003983,0.0003764,0.0003563,0.0003233,0.0002831,0.0002513,0.0001678, 5.32e-06},
{0.0233468,0.0230965,0.0231052,0.0231506,0.0231945,0.0232064,0.0232386,0.0232850,0.0233791,0.0231972,0.0231768,0.0230775,0.0229738,0.0226324,0.0222922,0.0218222,0.0213806,0.0208939,0.0203770,0.0197421,0.0190386,0.0184143,0.0174236,0.0167037,0.0156117,0.0141438,0.0125827,0.0106954,0.0088224,0.0068976,0.0052739,0.0039657,0.0029580,0.0022441,0.0015653,0.0012667,0.0011438,0.0010624,0.0010281,0.0009401,0.0008767,0.0008034,0.0006922,0.0006098,0.0004492,0.0000166},
{0.0483921,0.0482347,0.0481649,0.0480141,0.0480556,0.0483987,0.0483949,0.0482479,0.0483465,0.0483167,0.0482343,0.0478592,0.0473706,0.0469146,0.0463625,0.0456703,0.0445259,0.0433140,0.0424036,0.0411171,0.0399770,0.0382726,0.0365783,0.0347798,0.0326243,0.0295485,0.0261005,0.0224963,0.0185397,0.0144132,0.0109423,0.0082104,0.0061879,0.0046494,0.0033442,0.0025591,0.0023625,0.0022104,0.0021180,0.0019480,0.0018239,0.0016903,0.0014556,0.0012483,0.0009651,0.0000309},
{0.0758801,0.0756903,0.0760622,0.0756026,0.0757827,0.0758213,0.0760040,0.0761423,0.0760854,0.0761600,0.0755491,0.0750763,0.0747804,0.0740596,0.0728222,0.0714937,0.0701440,0.0685231,0.0666026,0.0650271,0.0632030,0.0607410,0.0580390,0.0547453,0.0511176,0.0467401,0.0413606,0.0356473,0.0292405,0.0228939,0.0171928,0.0130000,0.0099808,0.0072867,0.0051514,0.0040341,0.0037497,0.0034741,0.0033174,0.0030812,0.0029389,0.0026596,0.0022939,0.0019768,0.0014705,0.0000403},
{0.1012224,0.1010571,0.1006720,0.1006277,0.1006637,0.1010958,0.1012799,0.1015099,0.1014230,0.1012462,0.1010452,0.1004141,0.0998530,0.0986650,0.0973793,0.0951996,0.0935913,0.0913048,0.0891585,0.0867357,0.0841698,0.0810530,0.0772831,0.0733146,0.0685050,0.0625429,0.0557629,0.0477436,0.0390603,0.0305663,0.0230865,0.0176311,0.0132426,0.0097485,0.0069101,0.0054085,0.0049487,0.0046592,0.0044360,0.0040888,0.0038786,0.0035335,0.0030780,0.0026360,0.0019248,0.0000464},
{0.1246922,0.1247666,0.1243063,0.1244053,0.1249370,0.1251482,0.1257678,0.1256978,0.1260283,0.1260232,0.1262004,0.1253433,0.1246224,0.1239608,0.1224645,0.1209010,0.1186317,0.1162754,0.1139857,0.1112816,0.1084952,0.1052662,0.1009542,0.0963840,0.0908405,0.0838752,0.0750049,0.0643318,0.0522909,0.0401838,0.0297540,0.0221668,0.0165923,0.0123399,0.0090147,0.0069012,0.0064275,0.0059512,0.0056535,0.0052675,0.0048867,0.0043766,0.0038021,0.0032324,0.0023836,0.0000505},
{0.1442271,0.1437736,0.1437649,0.1441738,0.1447983,0.1457840,0.1459512,0.1465725,0.1468403,0.1475161,0.1474820,0.1472417,0.1468982,0.1465519,0.1455847,0.1444660,0.1423928,0.1402028,0.1382575,0.1357448,0.1331195,0.1296724,0.1259095,0.1213144,0.1153866,0.1074797,0.0965323,0.0829821,0.0670155,0.0503672,0.0362452,0.0263674,0.0197093,0.0147729,0.0107453,0.0084895,0.0078013,0.0072407,0.0067742,0.0063431,0.0059367,0.0051176,0.0043416,0.0036692,0.0027069,0.0000560},
{0.1528675,0.1526281,0.1529979,0.1529208,0.1538694,0.1547747,0.1553186,0.1560060,0.1568375,0.1573928,0.1580652,0.1581966,0.1584616,0.1580784,0.1573604,0.1567289,0.1549101,0.1525468,0.1508609,0.1488294,0.1462271,0.1430437,0.1393923,0.1350927,0.1290066,0.1202656,0.1082907,0.0932601,0.0749565,0.0558449,0.0397450,0.0286978,0.0211627,0.0159561,0.0116062,0.0092360,0.0084939,0.0079306,0.0074181,0.0069457,0.0063278,0.0054532,0.0046383,0.0037855,0.0028379,0.0000571},
{0.1608656,0.1607325,0.1608718,0.1618449,0.1621695,0.1633889,0.1644501,0.1646745,0.1652473,0.1668422,0.1674055,0.1677331,0.1681491,0.1684379,0.1681184,0.1669766,0.1655934,0.1637998,0.1617614,0.1598496,0.1577674,0.1546477,0.1512236,0.1466349,0.1409708,0.1317433,0.1190973,0.1024302,0.0821261,0.0610375,0.0427338,0.0304852,0.0225819,0.0170428,0.0125921,0.0099416,0.0091048,0.0084004,0.0079176,0.0074224,0.0068164,0.0057881,0.0048245,0.0040473,0.0030486,0.0000571},
{0.1632814,0.1631382,0.1631301,0.1638763,0.1647497,0.1659796,0.1670049,0.1674660,0.1686948,0.1697457,0.1703425,0.1712479,0.1717823,0.1723608,0.1722075,0.1714204,0.1701474,0.1687945,0.1670102,0.1651853,0.1632763,0.1605244,0.1568780,0.1532890,0.1480234,0.1378986,0.1251269,0.1078858,0.0861429,0.0635977,0.0443698,0.0316801,0.0231883,0.0174428,0.0128833,0.0102490,0.0094264,0.0087470,0.0082017,0.0076302,0.0070471,0.0059067,0.0049555,0.0040858,0.0030494,0.0000494},
{0.1659979,0.1658575,0.1660901,0.1669401,0.1682111,0.1691097,0.1699093,0.1708971,0.1719042,0.1732699,0.1742688,0.1751389,0.1763702,0.1767642,0.1770060,0.1765972,0.1747653,0.1737459,0.1720852,0.1709069,0.1690215,0.1666569,0.1631928,0.1596166,0.1537821,0.1446469,0.1311761,0.1131827,0.0905382,0.0666023,0.0459167,0.0326190,0.0238895,0.0182204,0.0133233,0.0105847,0.0097489,0.0090424,0.0084362,0.0079708,0.0072388,0.0061016,0.0050115,0.0041627,0.0031286,0.0000407},
{0.1679459,0.1680690,0.1681781,0.1691167,0.1700899,0.1711928,0.1722009,0.1730473,0.1737169,0.1760435,0.1767614,0.1778273,0.1791676,0.1793163,0.1797325,0.1791974,0.1783114,0.1768019,0.1750386,0.1737687,0.1722388,0.1697380,0.1670456,0.1641593,0.1581930,0.1489091,0.1353712,0.1169874,0.0932432,0.0682311,0.0474873,0.0332354,0.0245235,0.0183461,0.0136452,0.0109095,0.0099374,0.0092728,0.0086562,0.0081197,0.0074100,0.0061957,0.0050822,0.0042417,0.0031395,0.0000441},
{0.1682379,0.1683580,0.1682941,0.1690784,0.1703542,0.1713316,0.1720106,0.1730809,0.1742709,0.1760213,0.1767063,0.1783595,0.1794918,0.1798812,0.1801363,0.1800330,0.1788477,0.1773114,0.1760539,0.1744026,0.1732268,0.1711020,0.1681227,0.1650694,0.1594550,0.1500577,0.1366773,0.1183120,0.0943574,0.0691214,0.0477284,0.0336117,0.0246009,0.0185007,0.0137197,0.0109517,0.0100192,0.0092464,0.0088081,0.0080775,0.0073561,0.0061357,0.0050954,0.0042556,0.0031595,0.0000383},
{0.1697463,0.1699891,0.1701738,0.1707088,0.1722102,0.1726993,0.1739983,0.1748175,0.1760479,0.1774452,0.1786797,0.1800087,0.1809010,0.1812247,0.1817427,0.1811214,0.1803723,0.1790106,0.1772053,0.1758209,0.1745945,0.1723029,0.1696594,0.1665212,0.1612226,0.1523710,0.1384494,0.1197063,0.0959186,0.0704332,0.0486166,0.0341031,0.0250982,0.0188858,0.0138878,0.0110051,0.0101346,0.0093712,0.0088317,0.0081476,0.0074170,0.0062136,0.0050799,0.0043261,0.0031859,0.0000360},
{0.1723578,0.1718526,0.1725093,0.1732871,0.1742138,0.1748967,0.1760609,0.1771698,0.1778829,0.1794743,0.1807425,0.1816128,0.1829970,0.1837378,0.1838264,0.1837401,0.1825555,0.1813429,0.1792852,0.1781355,0.1773210,0.1744104,0.1718426,0.1686594,0.1639050,0.1546809,0.1413118,0.1221039,0.0979602,0.0717557,0.0493972,0.0347293,0.0252993,0.0190984,0.0139644,0.0111046,0.0101691,0.0094645,0.0088011,0.0081695,0.0074471,0.0062187,0.0052360,0.0043223,0.0032066,0.0000332},
{0.1745932,0.1747632,0.1747502,0.1757468,0.1770573,0.1776513,0.1786560,0.1791700,0.1803876,0.1818454,0.1827257,0.1840349,0.1850722,0.1858241,0.1863105,0.1862379,0.1848526,0.1837029,0.1818599,0.1800396,0.1783167,0.1762341,0.1735621,0.1704850,0.1656707,0.1574733,0.1433653,0.1243608,0.0998119,0.0730083,0.0501274,0.0350315,0.0255796,0.0193604,0.0140360,0.0110673,0.0101925,0.0094935,0.0088579,0.0082072,0.0075359,0.0062954,0.0051487,0.0042997,0.0032158,0.0000279},
{0.1756867,0.1753934,0.1760311,0.1769308,0.1776187,0.1785934,0.1795595,0.1799868,0.1810400,0.1822364,0.1832771,0.1844072,0.1855310,0.1864903,0.1869536,0.1870354,0.1853740,0.1836939,0.1823674,0.1805612,0.1791091,0.1767764,0.1735559,0.1712398,0.1660820,0.1576939,0.1441408,0.1252839,0.1003815,0.0735027,0.0505657,0.0351233,0.0256941,0.0193140,0.0140718,0.0110865,0.0101022,0.0094160,0.0088605,0.0083193,0.0074671,0.0062745,0.0051042,0.0042871,0.0031768,0.0000260},
{0.1756735,0.1761872,0.1757410,0.1769425,0.1778718,0.1787178,0.1799657,0.1802660,0.1815727,0.1824358,0.1837421,0.1846247,0.1856390,0.1866145,0.1870141,0.1868386,0.1856079,0.1841585,0.1820645,0.1801604,0.1784288,0.1765197,0.1735698,0.1704588,0.1659759,0.1576624,0.1446256,0.1257061,0.1008462,0.0734183,0.0503710,0.0354345,0.0256786,0.0193574,0.0140132,0.0110132,0.0100841,0.0094049,0.0087257,0.0081736,0.0075229,0.0062319,0.0052109,0.0043340,0.0031380,0.0000234},
{0.1731203,0.1732296,0.1732185,0.1740528,0.1751212,0.1762285,0.1773744,0.1779097,0.1788100,0.1794630,0.1810688,0.1819570,0.1829881,0.1838758,0.1843229,0.1841186,0.1829740,0.1812630,0.1789915,0.1772356,0.1762677,0.1738526,0.1703091,0.1678328,0.1637242,0.1556437,0.1425619,0.1237529,0.0995759,0.0728100,0.0496927,0.0346931,0.0253947,0.0190051,0.0137666,0.0108198,0.0098775,0.0091082,0.0085174,0.0080062,0.0072739,0.0060705,0.0050707,0.0042181,0.0030829,0.0000219},
{0.1686611,0.1686307,0.1687291,0.1698162,0.1707297,0.1716183,0.1726458,0.1730775,0.1738211,0.1750758,0.1761367,0.1774117,0.1783086,0.1792999,0.1798754,0.1792731,0.1780671,0.1765125,0.1748718,0.1728483,0.1711297,0.1688692,0.1661186,0.1634424,0.1594343,0.1518865,0.1392860,0.1211887,0.0976096,0.0714707,0.0486354,0.0338903,0.0246818,0.0184856,0.0135223,0.0103881,0.0096690,0.0089171,0.0083685,0.0077974,0.0071091,0.0058929,0.0049184,0.0040784,0.0029981,0.0000213},
{0.1639619,0.1637563,0.1639640,0.1649286,0.1654371,0.1665137,0.1676155,0.1681063,0.1687936,0.1704878,0.1714880,0.1720956,0.1737693,0.1744366,0.1751693,0.1746481,0.1734392,0.1716841,0.1698537,0.1681855,0.1667444,0.1642386,0.1612460,0.1590173,0.1556799,0.1482829,0.1357991,0.1188386,0.0954452,0.0697683,0.0475519,0.0331849,0.0240332,0.0181361,0.0130931,0.0101044,0.0093353,0.0086145,0.0080743,0.0075514,0.0068914,0.0056912,0.0047876,0.0039210,0.0029020,0.0000141},
{0.1550875,0.1551112,0.1552895,0.1560993,0.1572948,0.1582253,0.1591395,0.1595091,0.1602714,0.1616524,0.1624328,0.1635719,0.1646522,0.1659378,0.1664224,0.1660034,0.1648952,0.1633793,0.1609546,0.1599734,0.1581297,0.1558375,0.1533900,0.1513146,0.1480143,0.1410513,0.1294662,0.1134349,0.0912126,0.0663900,0.0453791,0.0315338,0.0229768,0.0172577,0.0124880,0.0095640,0.0088053,0.0081061,0.0076709,0.0071325,0.0065304,0.0053960,0.0044997,0.0036899,0.0027338,0.0000141},
{0.1465917,0.1465845,0.1464735,0.1477723,0.1488705,0.1498814,0.1506650,0.1509674,0.1514752,0.1527135,0.1537395,0.1553154,0.1558028,0.1572679,0.1575225,0.1574023,0.1560656,0.1547932,0.1531732,0.1509870,0.1501235,0.1480786,0.1456085,0.1439408,0.1402535,0.1342731,0.1236072,0.1080486,0.0866479,0.0635291,0.0431069,0.0299725,0.0219791,0.0164447,0.0118051,0.0090413,0.0083161,0.0076204,0.0072449,0.0067565,0.0061372,0.0051078,0.0042479,0.0035272,0.0025957,0.0000134},
{0.1398758,0.1397174,0.1399519,0.1410773,0.1420347,0.1433076,0.1438023,0.1442616,0.1454026,0.1460430,0.1471923,0.1481908,0.1496775,0.1504471,0.1510528,0.1510239,0.1497640,0.1478711,0.1469150,0.1449438,0.1441977,0.1419012,0.1395553,0.1381086,0.1353468,0.1297693,0.1191789,0.1043751,0.0841612,0.0613442,0.0417376,0.0289431,0.0210743,0.0157898,0.0112909,0.0085651,0.0079246,0.0074051,0.0069024,0.0064556,0.0058552,0.0048545,0.0041024,0.0032786,0.0024464, 7.67e-06},
{0.1328705,0.1330843,0.1329380,0.1339523,0.1353367,0.1365708,0.1369721,0.1375022,0.1382760,0.1397974,0.1402239,0.1413949,0.1428219,0.1437110,0.1440771,0.1441012,0.1431791,0.1418098,0.1399178,0.1384675,0.1377361,0.1359259,0.1338539,0.1322532,0.1301676,0.1245706,0.1145559,0.1008931,0.0813261,0.0592445,0.0401404,0.0278462,0.0201110,0.0151983,0.0108202,0.0081808,0.0075184,0.0069661,0.0065621,0.0061289,0.0055830,0.0046526,0.0038716,0.0031659,0.0023105,0.0000119},
{0.1254322,0.1256026,0.1257018,0.1268492,0.1278326,0.1289140,0.1296903,0.1302752,0.1306735,0.1319101,0.1328914,0.1340094,0.1349203,0.1362483,0.1370204,0.1366843,0.1357977,0.1347925,0.1328745,0.1318777,0.1312596,0.1294036,0.1277105,0.1261712,0.1240993,0.1192803,0.1101534,0.0967919,0.0779265,0.0568194,0.0384965,0.0266599,0.0193144,0.0145783,0.0102432,0.0077870,0.0071810,0.0067299,0.0061849,0.0057749,0.0052590,0.0043457,0.0036302,0.0030277,0.0021698,0.0000100},
{0.1097035,0.1098417,0.1098818,0.1109474,0.1117078,0.1126428,0.1134982,0.1137953,0.1145742,0.1154901,0.1164573,0.1173742,0.1184460,0.1195397,0.1201140,0.1200782,0.1188219,0.1179749,0.1164993,0.1159727,0.1151440,0.1138980,0.1121063,0.1112111,0.1095225,0.1055593,0.0977555,0.0858575,0.0690728,0.0505080,0.0341608,0.0235744,0.0171169,0.0127180,0.0090626,0.0067578,0.0061802,0.0057781,0.0054963,0.0050786,0.0045548,0.0037964,0.0032083,0.0026200,0.0019574, 6.18e-06},
{0.0870754,0.0872909,0.0871497,0.0878895,0.0887625,0.0894892,0.0900234,0.0904880,0.0908545,0.0918292,0.0922784,0.0934914,0.0942526,0.0950532,0.0952909,0.0955387,0.0950319,0.0939717,0.0929293,0.0923962,0.0920164,0.0908066,0.0897636,0.0891940,0.0879444,0.0848204,0.0785472,0.0693970,0.0558068,0.0408520,0.0275182,0.0190204,0.0137465,0.0102882,0.0071764,0.0053712,0.0049289,0.0045778,0.0042618,0.0040075,0.0036488,0.0029917,0.0025506,0.0020816,0.0014858, 4.69e-06},
{0.0685612,0.0685480,0.0687099,0.0694961,0.0699936,0.0704011,0.0707915,0.0711184,0.0715838,0.0723321,0.0730556,0.0737137,0.0743987,0.0752518,0.0754705,0.0754990,0.0750895,0.0746388,0.0737830,0.0733291,0.0731872,0.0724720,0.0713823,0.0714094,0.0705018,0.0679474,0.0628956,0.0556701,0.0450153,0.0329934,0.0222878,0.0153491,0.0110400,0.0082520,0.0056869,0.0042204,0.0039140,0.0036126,0.0033870,0.0030794,0.0028592,0.0023908,0.0019906,0.0016541,0.0012015, 3.83e-06},
{0.0612449,0.0612085,0.0611747,0.0619687,0.0624784,0.0628151,0.0635001,0.0639018,0.0640920,0.0648954,0.0655870,0.0659005,0.0667289,0.0672464,0.0676882,0.0678009,0.0673617,0.0669184,0.0662530,0.0655855,0.0657272,0.0650995,0.0646006,0.0642577,0.0635902,0.0618799,0.0575048,0.0507325,0.0410807,0.0301536,0.0203370,0.0139770,0.0101246,0.0074584,0.0051702,0.0037774,0.0034882,0.0032477,0.0030479,0.0028147,0.0025668,0.0021761,0.0017721,0.0014703,0.0011022, 2.77e-06},
{0.0490143,0.0492684,0.0494984,0.0497559,0.0502839,0.0506605,0.0510539,0.0511114,0.0515983,0.0520933,0.0526260,0.0531864,0.0539442,0.0541968,0.0546822,0.0547913,0.0544946,0.0538360,0.0535076,0.0533450,0.0535116,0.0530034,0.0526055,0.0526087,0.0521340,0.0506351,0.0471472,0.0419463,0.0339636,0.0247874,0.0166654,0.0114860,0.0081766,0.0060524,0.0042341,0.0030820,0.0028300,0.0025793,0.0024392,0.0022658,0.0020550,0.0017504,0.0014626,0.0012215,0.0008801, 1.92e-06},
{0.0422028,0.0422643,0.0424373,0.0428705,0.0434581,0.0437348,0.0440045,0.0441433,0.0445250,0.0450611,0.0454912,0.0458019,0.0464435,0.0468473,0.0471674,0.0475054,0.0471798,0.0466620,0.0466933,0.0464253,0.0462763,0.0459834,0.0457725,0.0460170,0.0457563,0.0445870,0.0416560,0.0369218,0.0299702,0.0219050,0.0147751,0.0101276,0.0072354,0.0053263,0.0036432,0.0026760,0.0024192,0.0022426,0.0021120,0.0019474,0.0017919,0.0015203,0.0012686,0.0010515,0.0007335, 1.06e-06},
{0.0352654,0.0353476,0.0355823,0.0360315,0.0362128,0.0365717,0.0368053,0.0370362,0.0373712,0.0377448,0.0379600,0.0383983,0.0390288,0.0394543,0.0396950,0.0398400,0.0396213,0.0396933,0.0393257,0.0391265,0.0393024,0.0389977,0.0390400,0.0391461,0.0390403,0.0381446,0.0356920,0.0319338,0.0259740,0.0189048,0.0128064,0.0086696,0.0061974,0.0045674,0.0030503,0.0022647,0.0020341,0.0018718,0.0018049,0.0016754,0.0015007,0.0012658,0.0010398,0.0008758,0.0006230, 8.52e-07},
{0.0282807,0.0284937,0.0284948,0.0288028,0.0291148,0.0293902,0.0296085,0.0298211,0.0300226,0.0303327,0.0306332,0.0308769,0.0312245,0.0317314,0.0318686,0.0319602,0.0319766,0.0319133,0.0315864,0.0317570,0.0316264,0.0316886,0.0316058,0.0319195,0.0319544,0.0312656,0.0293791,0.0260658,0.0213862,0.0156503,0.0104175,0.0072153,0.0050520,0.0037836,0.0025229,0.0017555,0.0016688,0.0014941,0.0014437,0.0013284,0.0012075,0.0010155,0.0008541,0.0006854,0.0004990, 1.49e-06},
{0.0218501,0.0218946,0.0218982,0.0219904,0.0223723,0.0225227,0.0228890,0.0228298,0.0230062,0.0233001,0.0234496,0.0237425,0.0241996,0.0243740,0.0245678,0.0245945,0.0245632,0.0245171,0.0243499,0.0243378,0.0245163,0.0245127,0.0245817,0.0247989,0.0249544,0.0245335,0.0231404,0.0207465,0.0169323,0.0123817,0.0083359,0.0057101,0.0040383,0.0029555,0.0020083,0.0014085,0.0012456,0.0011676,0.0011118,0.0010236,0.0009401,0.0008038,0.0006447,0.0005414,0.0003951, 1.28e-06},
{0.0153338,0.0155044,0.0154160,0.0155546,0.0157519,0.0159853,0.0160026,0.0161802,0.0162569,0.0165493,0.0165740,0.0168381,0.0171052,0.0173376,0.0173242,0.0176192,0.0175152,0.0174070,0.0172777,0.0174967,0.0175898,0.0174918,0.0175695,0.0180175,0.0179540,0.0177870,0.0168317,0.0151870,0.0122856,0.0090645,0.0060901,0.0040922,0.0029267,0.0021670,0.0014379,0.0010032,0.0008765,0.0008594,0.0007810,0.0007502,0.0006526,0.0005559,0.0004584,0.0003930,0.0002767, 4.26e-07},
{0.0110686,0.0111229,0.0111421,0.0111793,0.0113440,0.0115046,0.0115091,0.0115632,0.0116170,0.0118460,0.0120298,0.0120899,0.0123054,0.0123267,0.0125670,0.0126181,0.0125715,0.0126319,0.0126017,0.0125736,0.0127161,0.0126707,0.0128912,0.0130524,0.0132070,0.0130818,0.0123555,0.0111968,0.0092251,0.0067195,0.0045257,0.0030809,0.0021936,0.0015755,0.0010415,0.0007056,0.0006564,0.0006060,0.0005681,0.0005299,0.0004731,0.0004222,0.0003346,0.0002722,0.0002179, 1.00e-07},
{0.0071135,0.0073078,0.0072505,0.0073018,0.0073659,0.0073651,0.0074892,0.0075007,0.0075657,0.0076854,0.0077188,0.0077719,0.0079604,0.0081027,0.0081764,0.0082486,0.0082420,0.0081193,0.0082539,0.0083744,0.0083120,0.0083885,0.0084933,0.0086443,0.0088309,0.0086833,0.0082524,0.0074703,0.0062207,0.0044794,0.0030688,0.0020556,0.0014309,0.0010477,0.0006888,0.0004728,0.0004307,0.0003838,0.0003810,0.0003376,0.0003146,0.0002573,0.0002124,0.0001813,0.0001378, 1.00e-07},
{0.0044420,0.0044505,0.0045314,0.0044933,0.0045670,0.0045979,0.0046394,0.0046839,0.0046990,0.0047842,0.0048219,0.0048775,0.0049785,0.0049947,0.0051078,0.0051278,0.0051755,0.0051949,0.0051103,0.0051529,0.0052038,0.0052952,0.0052854,0.0054439,0.0055534,0.0054941,0.0052522,0.0048228,0.0039898,0.0029033,0.0019683,0.0013318,0.0009338,0.0006950,0.0004573,0.0002927,0.0002709,0.0002560,0.0002458,0.0002100,0.0002085,0.0001638,0.0001346,0.0001199,0.0000818, 1.00e-07},
{0.0024513,0.0024194,0.0024288,0.0024886,0.0024841,0.0025148,0.0025188,0.0025357,0.0025966,0.0026243,0.0026356,0.0026777,0.0027450,0.0027619,0.0027655,0.0028019,0.0028292,0.0028628,0.0028228,0.0028520,0.0028305,0.0028843,0.0029551,0.0030639,0.0030771,0.0030507,0.0029331,0.0026871,0.0022026,0.0016226,0.0011099,0.0007146,0.0005220,0.0003817,0.0002381,0.0001597,0.0001561,0.0001272,0.0001252,0.0001105,0.0001161,0.0000933,0.0000748,0.0000609,0.0000526, 1.00e-07},
{0.0013180,0.0012848,0.0012956,0.0012963,0.0013286,0.0013248,0.0013638,0.0013482,0.0013789,0.0014072,0.0014375,0.0014072,0.0014432,0.0014592,0.0014843,0.0014916,0.0014731,0.0015229,0.0015108,0.0014892,0.0015685,0.0015612,0.0015987,0.0016505,0.0016699,0.0016696,0.0015977,0.0015105,0.0012332,0.0009063,0.0006002,0.0003902,0.0002699,0.0002089,0.0001350,0.0000822,0.0000797,0.0000728,0.0000724,0.0000654,0.0000599,0.0000439,0.0000411,0.0000347,0.0000245, 1.00e-07},
{0.0006441,0.0006545,0.0006349,0.0006403,0.0006658,0.0006777,0.0006620,0.0006628,0.0006790,0.0006961,0.0007146,0.0007046,0.0007125,0.0007069,0.0007229,0.0007404,0.0007448,0.0007431,0.0007306,0.0007284,0.0007493,0.0007736,0.0007974,0.0007998,0.0008128,0.0008230,0.0008273,0.0007174,0.0006530,0.0004564,0.0003035,0.0002241,0.0001508,0.0001059,0.0000754,0.0000426,0.0000371,0.0000334,0.0000341,0.0000351,0.0000332,0.0000262,0.0000224,0.0000158,0.0000126, 1.00e-07},
{0.0003301,0.0003152,0.0003093,0.0003267,0.0003346,0.0003112,0.0003378,0.0003365,0.0003331,0.0003419,0.0003531,0.0003465,0.0003614,0.0003525,0.0003723,0.0003951,0.0003983,0.0003736,0.0003727,0.0003894,0.0003885,0.0003985,0.0004209,0.0004347,0.0004449,0.0004241,0.0004315,0.0004034,0.0003225,0.0002394,0.0001653,0.0001131,0.0000756,0.0000560,0.0000354,0.0000215,0.0000239,0.0000166,0.0000181,0.0000155,0.0000132,0.0000145, 8.95e-06,0.0000104, 5.75e-06, 1.00e-07},
{0.0001485,0.0001412,0.0001333,0.0001427,0.0001399,0.0001438,0.0001480,0.0001397,0.0001602,0.0001436,0.0001468,0.0001387,0.0001463,0.0001536,0.0001499,0.0001442,0.0001642,0.0001583,0.0001736,0.0001546,0.0001651,0.0001738,0.0001834,0.0001940,0.0001904,0.0001840,0.0001772,0.0001783,0.0001308,0.0001108,0.0000667,0.0000486,0.0000305,0.0000266,0.0000158,0.0000100, 8.52e-06, 7.45e-06, 6.39e-06, 5.75e-06, 7.24e-06, 5.54e-06, 3.83e-06, 4.26e-06, 3.41e-06, 1.00e-07},
{0.0000775,0.0000837,0.0000767,0.0000852,0.0000788,0.0000831,0.0000880,0.0000837,0.0000807,0.0000818,0.0000871,0.0000858,0.0000912,0.0000878,0.0000884,0.0000927,0.0000988,0.0000920,0.0000958,0.0000963,0.0000967,0.0000916,0.0001103,0.0001039,0.0001146,0.0001086,0.0001061,0.0000971,0.0000807,0.0000667,0.0000411,0.0000281,0.0000192,0.0000151, 9.37e-06, 5.75e-06, 6.60e-06, 5.32e-06, 4.90e-06, 5.11e-06, 3.19e-06, 2.56e-06, 1.70e-06, 1.49e-06, 2.34e-06, 1.00e-07},
{0.0000637,0.0000573,0.0000633,0.0000675,0.0000652,0.0000579,0.0000590,0.0000645,0.0000662,0.0000656,0.0000739,0.0000741,0.0000692,0.0000756,0.0000728,0.0000724,0.0000758,0.0000703,0.0000737,0.0000812,0.0000786,0.0000733,0.0000856,0.0000822,0.0000875,0.0000939,0.0000854,0.0000814,0.0000669,0.0000599,0.0000351,0.0000228,0.0000168,0.0000128, 6.39e-06, 5.54e-06, 5.32e-06, 3.19e-06, 3.83e-06, 3.19e-06, 2.77e-06, 2.13e-06, 2.77e-06, 2.13e-06, 2.56e-06, 1.00e-07},
{0.0000515,0.0000498,0.0000477,0.0000469,0.0000520,0.0000486,0.0000501,0.0000526,0.0000507,0.0000481,0.0000518,0.0000520,0.0000652,0.0000532,0.0000556,0.0000620,0.0000635,0.0000630,0.0000643,0.0000609,0.0000571,0.0000720,0.0000611,0.0000699,0.0000786,0.0000703,0.0000741,0.0000694,0.0000547,0.0000371,0.0000315,0.0000185,0.0000136, 7.03e-06, 5.32e-06, 3.83e-06, 4.69e-06, 2.98e-06, 2.77e-06, 1.06e-06, 2.34e-06, 2.13e-06, 1.70e-06, 1.70e-06, 6.39e-07, 1.00e-07},
{0.0000383,0.0000377,0.0000411,0.0000356,0.0000417,0.0000452,0.0000392,0.0000443,0.0000426,0.0000424,0.0000409,0.0000462,0.0000458,0.0000494,0.0000447,0.0000458,0.0000469,0.0000456,0.0000475,0.0000456,0.0000524,0.0000458,0.0000528,0.0000532,0.0000575,0.0000605,0.0000539,0.0000543,0.0000435,0.0000283,0.0000198,0.0000143,0.0000121, 8.09e-06, 4.26e-06, 2.77e-06, 2.56e-06, 2.13e-06, 1.70e-06, 2.34e-06, 2.77e-06, 1.28e-06, 1.06e-06, 1.49e-06, 6.39e-07, 1.00e-07},
{0.0000307,0.0000311,0.0000292,0.0000309,0.0000283,0.0000290,0.0000328,0.0000315,0.0000298,0.0000292,0.0000339,0.0000317,0.0000354,0.0000313,0.0000315,0.0000343,0.0000328,0.0000366,0.0000351,0.0000394,0.0000322,0.0000371,0.0000364,0.0000392,0.0000441,0.0000447,0.0000403,0.0000456,0.0000302,0.0000279,0.0000175,0.0000100, 7.88e-06, 4.05e-06, 4.26e-06, 3.62e-06, 3.41e-06, 8.52e-07, 1.28e-06, 2.13e-06, 8.52e-07, 1.49e-06, 1.06e-06, 1.70e-06, 4.26e-07, 1.00e-07},
{0.0000213,0.0000209,0.0000230,0.0000251,0.0000194,0.0000247,0.0000222,0.0000232,0.0000226,0.0000194,0.0000288,0.0000266,0.0000268,0.0000288,0.0000260,0.0000264,0.0000298,0.0000294,0.0000264,0.0000296,0.0000232,0.0000313,0.0000311,0.0000277,0.0000324,0.0000302,0.0000337,0.0000300,0.0000290,0.0000241,0.0000153, 8.52e-06, 5.54e-06, 1.70e-06, 2.13e-06, 1.06e-06, 1.06e-06, 2.13e-06, 2.13e-06, 6.39e-07, 1.06e-06, 1.49e-06, 8.52e-07, 8.52e-07, 1.00e-07, 1.00e-07},
{0.0000190,0.0000181,0.0000158,0.0000192,0.0000196,0.0000170,0.0000173,0.0000198,0.0000175,0.0000151,0.0000181,0.0000179,0.0000190,0.0000198,0.0000183,0.0000213,0.0000213,0.0000192,0.0000224,0.0000236,0.0000204,0.0000232,0.0000239,0.0000279,0.0000228,0.0000211,0.0000277,0.0000277,0.0000204,0.0000147, 9.80e-06, 8.95e-06, 6.18e-06, 5.54e-06, 2.34e-06, 1.49e-06, 1.28e-06, 1.28e-06, 8.52e-07, 8.52e-07, 1.28e-06, 1.06e-06, 6.39e-07, 6.39e-07, 1.00e-07, 1.00e-07},
{0.0000119,0.0000136,0.0000134,0.0000147,0.0000126,0.0000134,0.0000119,0.0000128,0.0000141,0.0000126,0.0000106,0.0000138,0.0000168,0.0000134,0.0000124,0.0000149,0.0000164,0.0000132,0.0000136,0.0000179,0.0000190,0.0000126,0.0000162,0.0000175,0.0000185,0.0000209,0.0000158,0.0000200,0.0000158,0.0000104, 7.88e-06, 3.41e-06, 2.98e-06, 3.19e-06, 1.28e-06, 1.28e-06, 1.49e-06, 1.06e-06, 4.26e-07, 6.39e-07, 8.52e-07, 6.39e-07, 8.52e-07, 8.52e-07, 2.13e-07, 1.00e-07},
{ 8.95e-06, 7.88e-06, 7.67e-06, 9.80e-06, 8.31e-06, 9.16e-06, 6.60e-06, 8.52e-06,0.0000115, 7.03e-06,0.0000106,0.0000109,0.0000102, 7.67e-06, 9.58e-06,0.0000104,0.0000104,0.0000102,0.0000106, 8.09e-06,0.0000134,0.0000143,0.0000117,0.0000117,0.0000115,0.0000134,0.0000104, 9.80e-06, 8.52e-06,0.0000106, 3.62e-06, 2.34e-06, 3.41e-06, 1.92e-06, 4.26e-07, 6.39e-07, 2.13e-06, 2.13e-07, 2.13e-07, 6.39e-07, 2.13e-07, 6.39e-07, 2.13e-07, 1.00e-07, 1.00e-07, 1.00e-07},
{ 4.47e-06, 4.90e-06, 3.62e-06, 3.19e-06, 4.47e-06, 2.13e-06, 3.83e-06, 4.26e-06, 4.05e-06, 2.77e-06, 5.11e-06, 2.98e-06, 6.18e-06, 4.26e-06, 3.62e-06, 4.26e-06, 4.90e-06, 5.32e-06, 5.96e-06, 5.11e-06, 5.96e-06, 7.45e-06, 7.88e-06, 4.90e-06, 5.32e-06, 6.82e-06, 7.03e-06, 6.18e-06, 4.47e-06, 2.56e-06, 2.98e-06, 1.06e-06, 1.28e-06, 6.39e-07, 8.52e-07, 1.00e-07, 1.00e-07, 1.00e-07, 2.13e-07, 2.13e-07, 2.13e-07, 2.13e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07},
{ 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07, 1.00e-07}
}
};


// Static class variables:
double PEnergy::fPIP[fNPMTsDim];
double PEnergy::fPMTEff[fNPMTsDim];
vector< vector<double> > PEnergy::fPQFac;
vector< vector<double> > PEnergy::fPTFac;
unsigned int  PEnergy::fXIdx[fNPMTsDim];
double PEnergy::fExpNoise[fNPMTsDim];
bool   PEnergy::fOnLine[fNPMTsDim];
bool   PEnergy::fUsePMT[fNPMTsDim];
bool   PEnergy::fWasHit[fNPMTsDim];

TVector3 PEnergy::fPosPMT[fNPMTsDim];
TVector3 PEnergy::fDirPMT[fNPMTsDim];
double PEnergy::fPosPMTUX[fNPMTsDim];
double PEnergy::fPosPMTUY[fNPMTsDim];
double PEnergy::fPosPMTUZ[fNPMTsDim];
double PEnergy::fPhotonTransitTime[fNPMTsDim];

Int_t  PEnergy::fPrevGTID;
bool   PEnergy::fPrevSeedActual;
bool   PEnergy::fPrevElectronPrime;
bool   PEnergy::fPrevUseDirection;
int    PEnergy::fPrevNPhotonSim;
double PEnergy::fPrevTabAreaFactor;
double PEnergy::fPrevEffScaleFactor;
int    PEnergy::fPrevChargeType;
string PEnergy::fPrevSelector;
double PEnergy::fPrevLowCut, PEnergy::fPrevHighCut;
bool   PEnergy::fGoodPrevAuxSim = true;

// **********************************************************************
PEnergy::PEnergy() {
}

// **********************************************************************
void PEnergy::Initialise( const std::string& ) {
    fRndmNoSeed = 0;
    // Default values:
    fSeedActual = false;
    fUseCharge = true;
    fUseTime = false;
    fChargeType = 1;    // Current default charge type is QHL
    fElectronPrime = true;
    fNPhotonSim = 100000;  // default for scintillator.  Water default is
                                // 200,000, set in BeginOfRun
    fTabAreaFactor = 2.;
//     fEffScaleFactor = 0.609646*(9.91524/10.)*(6.74992/7.)*(6.96544/7.);
    fEffScaleFactor = 0.5800054*0.99187362*1.003;
    fUseDirection = false;      // reset to true in BeginOfRun for water
    fMakeWaterAdj = true;       // reset to false if non-default setting of
                                // nPhotonSim is used in water
}

// **********************************************************************
void PEnergy::BeginOfRun( DS::Run& run ){
    fPMTChgpdfPtr = PMTChgpdf::Instance();
    fPMTTimeDistnsPtr = PMTTimeDistns::Instance();
    fPMTTimeDistnsPtr->BeginOfRun();
    fPMTTransitTimePtr = PMTTransitTime::Get();

    if(fRndmNoSeed==0) {
        time_t start_time = time(NULL);
        pid_t pid = getpid();
        fRndmNoSeed = start_time ^ (pid << 16);
    }
    detail << "PEnergy starting random number seed: " << fRndmNoSeed <<endl;
    // Save the state of the main random number sequence, initialize and save
    // a sequence for PEnergy, and restore the main sequence:
    CLHEP::HepRandom* hepRandom = CLHEP::HepRandom::getTheGenerator();
    hepRandom->saveStaticRandomStates(fMainSeqState);
    hepRandom->setTheSeed(fRndmNoSeed,0);
    hepRandom->saveStaticRandomStates(fLocalSeqState);
    hepRandom->restoreStaticRandomStates(fMainSeqState);

    fCallCounter = 0;
    fPrevGTID=-10;
    DB* db = DB::Get();

    RAT::DBLinkPtr ptrav =  db->GetLink("GEO","inner_av");
    G4String avMaterialName;
    try {
        avMaterialName = ptrav->GetS( "material" );
    } catch ( RAT::DBNotFoundError& e) {
        RAT::Log::Die(
                "\nPEnergy::BeginOfRun: Cannot find inner_av material" );
    }
    G4Material* avMaterialPtr = G4Material::GetMaterial(avMaterialName );
    if (avMaterialPtr == NULL) RAT::Log::Die
            ( "\nPEnergy::BeginOfRun: Cannot find inner_av material pointer" );
    if(avMaterialName=="lightwater_sno") {
        fLightYield = 0.;
        fUseDirection = true;
        if(fMakeWaterAdj) fNPhotonSim = 200000;  // default value for water
    }else{
        fMakeWaterAdj = false;
        // Material assumed to be a scintillator, with a light yield.
        fLightYield = avMaterialPtr->GetMaterialPropertiesTable()
                                        ->GetConstProperty("LIGHT_YIELD");
    }

    // Get coefficients from database for function to return the expected
    // number of photons given event energy:
    // If material is not one of the two for which coefficients are provided,
    // use labppo_scintillator coefficients (scaled in the function
    // ExpectedPhotons by scintillation light yield).
    fAVMaterialName = avMaterialName;
    if(fAVMaterialName!="labppo_scintillator" &&
                                           fAVMaterialName!="lightwater_sno") {
        fAVMaterialName="labppo_scintillator";
        detail << "\nPEnergy:  energy function coefficients not provided for "
            << avMaterialName <<".  \nUsing " << fAVMaterialName <<
            " coefficients.\n\n";
    }
    DBLinkPtr lPhotFns = db->GetLink( "PHOTON_PROD_FNS", fAVMaterialName );
    fCherCoef = lPhotFns->GetDArray("cherfn");
    fScintCoef = lPhotFns->GetDArray("scintfn");

    // Copied from PMTOpticalModel.cc
    DBLinkPtr pmtModelTable =
                    db->GetLink( "PMT_OPTICAL_PARAMETERS", "PMTOptics0_0" );
    fCollectionEfficiency = pmtModelTable->GetD( "collection_efficiency" );
//     detail << "PEnergy::BeginOfRun,  fCollectionEfficiency = " <<
//                                               fCollectionEfficiency << endl;
    fPEEscapeProbability = Optics::LoadWavelengthPropertyVector
                                    ( "pe_escape_probability", pmtModelTable );

    const DU::PMTInfo& pmtInfo = DU::Utility::Get()->GetPMTInfo();
    fNPMTs = pmtInfo.GetCount();
    if( fNPMTs > fNPMTsDim) Log::Die("PEnergy: fNPMTsDim insufficiently large");

    if(fUseTime) fUseCharge = true;

    if(fUseCharge) {
        // This section copied from UnCal.cc
        DBLinkPtr fUnCalbank = db->GetLink("DAQ");
        int ECApattern = fUnCalbank->GetI("eca_pdst_pattern");
        int ECArate = fUnCalbank->GetI("eca_pdst_rate");
       //Get the pattern and rate
        std::stringstream eca_pattern_rate;
        eca_pattern_rate << ECApattern <<"_"<<ECArate;
        DBLinkPtr fPedsbank = db->GetLink("ECA_PDST", eca_pattern_rate.str());
        DBLinkPtr fPmtCalib = db->GetLink("PMTCALIB");
        if(fChargeType==0) {
            fPedSig    = fPedsbank->GetFArrayFromD("pdst_qhswidth");
//          Currently getting thresholds from ChanHWStatus
//          fThreshold = fPmtCalib->GetFArrayFromD("QHS_threshold");
            fHHP       = fPmtCalib->GetFArrayFromD("QHS_hhp");
        }else if(fChargeType==1) {
            fPedSig    = fPedsbank->GetFArrayFromD("pdst_qhlwidth");
//          fThreshold = fPmtCalib->GetFArrayFromD("QHL_threshold");
            fHHP       = fPmtCalib->GetFArrayFromD("QHL_hhp");
        }else{
            fPedSig    = fPedsbank->GetFArrayFromD("pdst_qlxwidth");
//          fThreshold = fPmtCalib->GetFArrayFromD("QLX_threshold");
            fHHP       = fPmtCalib->GetFArrayFromD("QLX_hhp");
        }
    }

    for(unsigned int i=0; i<fPedSig.size(); i++) {
        // For FFT to work, can't have PedSig=0 or too large
        if( fPedSig[i] < 0.5 ) fPedSig[i] = 0.5;
        if( fPedSig[i] > 200. ) fPedSig[i] = 200.;
    }

    fPMTChargePtr = PMTCharge::Get();

    DBLinkPtr fLdaq = db->GetLink( "DAQ" );
    // Delay of channel info from FEC (front-end card) to trigger
    fFECToTriggerDelay = fLdaq->GetD("triggerfecdelay") ;
    // Add delay of global trigger back to FEC
    fTotalTriggerDelay = fFECToTriggerDelay+ fLdaq->GetD("gtriggerdelay");
#ifdef PENERGYDETAIL
    detail<<"PEnergy::BeginOfRun: fFECToTriggerDelay, FECfromTriggerdelay = "<<
        fFECToTriggerDelay<<" "<< fLdaq->GetD("gtriggerdelay") << endl;
#endif
    fNoiseWindow = fLdaq->GetD( "triggergate" );
    // fSelectors is inherited from SelectorMethod class.
    if( ++(fSelectors.begin()) != fSelectors.end() ) RAT::Log::Die(
                                        "PEnergy expects only one selector.");
    // If a non-null selector is being applied, use selector interval to
    // calculate the expected number of noise hits.
    fSelectorName = (*fSelectors.begin())->GetName();
    fSelector = *fSelectors.begin();
    fSelector->BeginOfRun( run );
    fLowCut  = -30.;
    fHighCut = +300.;
    if( fSelectorName=="timeResidualCut" ) {
        fTRCut = dynamic_cast<PMTSelectors::TimeResidualCut*> (fSelector);
        fLowCut  = fTRCut->GetLowCut();
        fHighCut = fTRCut->GetHighCut();
        fNoiseWindow = fHighCut - fLowCut;
    }else if( fSelectorName=="modeCut" ) {
        PMTSelectors::ModeCut* fModeCut =
                            dynamic_cast<PMTSelectors::ModeCut*> (fSelector);
        fLowCut  = fModeCut->GetLowCut();
        fHighCut = fModeCut->GetHighCut();
        fNoiseWindow = fHighCut - fLowCut;
    }else if( fSelectorName != "null" ) {
        RAT::Log::Die(
            "PEnergy expects the null, modeCut, or timeResidualCut selector");
    }
    detail<<"PEnergy:  Using the " << fSelectorName << " selector." << endl;
    if(fSelectorName != "null")
            detail<<"          Interval Limits = "<<fLowCut<<" "<<fHighCut<<endl;

    DBLinkPtr noiseMC = db->GetLink( "NOISE_MC" );
    DBLinkPtr noiseDB = db->GetLink( "NOISE_RUN_LEVEL" );
    Noise::ENoiseModel fNoiseFlag =
            static_cast<Noise::ENoiseModel>( noiseMC->GetI( "noise_flag" ) );

    if( fNoiseFlag == Noise::eNoNoise ) {
        // Use a small, but non-zero noise value
        for(unsigned int i=0; i<fNPMTs; i++) {
            fExpNoise[i] = 1.e-33;
        }
    }else if( fNoiseFlag == Noise::eAverageNoise ) {
        // Use average noise rate.  The rate from the database in is Hz;
        // divide by "second" to convert to inverse ns.  Divide by channel
        // efficiency because values in database are observed values.
        // Channel efficiency is reapplied later.

#ifdef PENERGYDETAIL
        detail << "Noise window = "<<fNoiseWindow <<" ns" << endl;
        detail << "expected noise PEs per PMT= " <<" "<<
            fNoiseWindow*noiseDB->GetD( "average_noiserate" )/second << endl;
#endif

        for(unsigned int i=0; i<fNPMTs; i++) {
            fExpNoise[i] = noiseDB->GetD( "average_noiserate" )/
                  (second*ChannelEfficiency::Get()->GetChannelEfficiency(i));
        }
    }else if( fNoiseFlag == Noise::ePerChannelNoise ) {
        // Get tube status and noise rate
        unsigned int fNumChannels =
                          db->GetLink("PMTINFO")->GetIArray("pmt_type").size();
        if( fNPMTs != fNumChannels) {
                    Log::Die( "PEnergy: fNPMTs and fNumChannels not equal" );
        }
        // Noise rates per channel (lcn):
        vector<double> fChannelRate = noiseDB->GetDArray( "perpmt_noiserate" );
        for(unsigned int i=0; i<fNPMTs; i++) {
            fExpNoise[i] = fChannelRate[i]/
                    (second*ChannelEfficiency::Get()->GetChannelEfficiency(i));
        }
    }

#ifdef PENERGYDETAIL
    detail << "fNoiseFlag = " << fNoiseFlag << endl;
    detail << "Samples of fExpNoise:"<<endl;
    int icount = 0;
    for(unsigned int i=0; i<fNPMTs; i+=100) {
        detail << fExpNoise[i] <<" ";
        icount++;
        if(icount==10) {
            detail << endl;
            icount=0;
        }
    }
    detail << endl;
#endif

    fGsimPtr = gTheRatRunManager->GetGsim();

    DBLinkPtr ptrInnerPMT = db->GetLink( "GEO", "innerPMT" );
    vector<string> conType = ptrInnerPMT->GetSArray("concentrator_type");
    DBLinkPtr ptrConcentrator = db->GetLink( "CONCENTRATOR", conType[0]);
    rConc = ptrConcentrator->GetD("lip_radius");

    fFirstNormal = 0;
    fLastNormal = 0;
    double minR = 10000.;
    rPMTMax = 0.;
    const DU::ChanHWStatus&
                chanHWStatus = DU::Utility::Get()->GetChanHWStatus();
    for(size_t k=0; k<fNPMTs; k++) {
        bool isNormal = (pmtInfo.GetType(k) == pmtInfo.NORMAL) ;
        fOnLine[k] = isNormal;
        // Determine indices of first and last normal PMTs in the arrays, so
        // subsequent loops can be shortened.  Also determine the least and
        // greatest distances of any PMT from the center of the detector.
        if(isNormal) {
            fLastNormal = k;
            if( fFirstNormal == 0) fFirstNormal = k;

            double rPMT = (pmtInfo.GetPosition(k)).Mag();
            if(rPMT<minR) minR = rPMT;
            if(rPMT>rPMTMax) rPMTMax = rPMT;

            if( chanHWStatus.IsEnabled() ) {    //  Skip off-line PMTs
#ifdef PENERGYDETAIL
                info << "PEnergy, Status flags for tube "<<k<<" = "
                     << fOnLine[k]
                     <<" "<< chanHWStatus.IsTubeOnline(k)
                     <<" "<< chanHWStatus.IsTubeAtHighVoltage(k)
                     <<" "<< chanHWStatus.IsChannelOnline(k)
                     <<" "<< chanHWStatus.IsDAQEnabled(k)
                     <<" "<< PMTVariation::Get()->GetVariation(k)<<endl;
#endif
                // Some of the following tests may be redundant??
                fOnLine[k] = fOnLine[k] && chanHWStatus.IsTubeOnline(k)
                                        && chanHWStatus.IsTubeAtHighVoltage(k)
                                        && chanHWStatus.IsChannelOnline(k)
                                        && chanHWStatus.IsDAQEnabled(k);
            }
        }
    }

    int nOnLine = 0;
    for(unsigned int k=fFirstNormal; k<=fLastNormal; k++) {
        if( fOnLine[k] ) {
            nOnLine++;
            fPosPMT[k] = pmtInfo.GetPosition(k);
            TVector3 posPMTUnit = fPosPMT[k].Unit();
            fPosPMTUX[k] = posPMTUnit.X();
            fPosPMTUY[k] = posPMTUnit.Y();
            fPosPMTUZ[k] = posPMTUnit.Z();

            fDirPMT[k] = pmtInfo.GetDirection(k);
        }
    }

#ifdef PENERGYDETAIL
    detail << "PEnergy::BeginOfRun:  fLightYield, tubes on-line = " <<
                                        fLightYield <<" "<< nOnLine << endl;
#endif

    fRTest = minR - 1.;   // Slightly smaller than smallest distance of PMT
                          // positions from detector center
    fRsqTest = fRTest*fRTest;
    double epsilon = acos(fRTest/rPMTMax);
    double delta = atan(rConc*sqrt(fTabAreaFactor)/fRTest) ;
    // Dot product limit used in first screen to determine if photon is inside
    // a PMT's tabulation area:
    fTabLim = cos(epsilon+delta);
    // Square of distance used in final test to determine if photon is incident
    // on tabulation area:
    fRhoSqLim = fTabAreaFactor*rConc*rConc;

    DU::LightPathCalculator lightPathCalc =
                                DU::Utility::Get()->GetLightPathCalculator();
    double photonMeV = 2.8e-6;  // typical energy of photon reaching PSUP
                                // Corresponds to wavelength of 442.8 nm
    fRefIdxAVFill = lightPathCalc.GetInnerAVRI(photonMeV);
    double refIdxAV = lightPathCalc.GetAVRI(photonMeV);
    double refIdxWater = lightPathCalc.GetWaterRI(photonMeV);
    // The inner and outer radii of the acrylic vessel
    fAVInnerRadius = db->GetLink( "SOLID", "acrylic_vessel_inner" )->
                                                      GetD( "r_sphere" );
    double AVOuterRadius = db->GetLink( "SOLID", "acrylic_vessel_outer" )->
                                                      GetD( "r_sphere" );
    fPropTimeAdder = ( (AVOuterRadius-fAVInnerRadius)*refIdxAV +
              (0.5*(minR+rPMTMax) - AVOuterRadius)*refIdxWater )/c_light;

    // Create profilers for accumulating time spent in various computations
    fProfileGenPhotons = new Profiler("PEnergy", "GenPhotons");
//     fProfileAccumData = new Profiler("PEnergy", "AccumData");
    fProfileCalcIncidProb = new Profiler("PEnergy", "CalcIncdPrb");
    fProfilePFactors = new Profiler("PEnergy", "PFactors");
    fProfileOptimize = new Profiler("PEnergy", "Optimize");
    fProfileFinish = new Profiler("PEnergy", "Finish");
    fProfileTotal = new Profiler("PEnergy", "Total");
}       // void PEnergy::BeginOfRun( DS::Run& )

// **********************************************************************

void PEnergy::EndOfRun( DS::Run& ){
    info << "\nTime usage by different parts of PEnergy:\n";
    fProfileGenPhotons->OutputTiming();
//     fProfileAccumData->OutputTiming();
    fProfileCalcIncidProb->OutputTiming();
    fProfilePFactors->OutputTiming();
    fProfileOptimize->OutputTiming();
    fProfileFinish->OutputTiming();
    fProfileTotal->OutputTiming();
    info << endl;

    fProfileGenPhotons->~Profiler();    fProfileGenPhotons = NULL;
//     fProfileAccumData->~Profiler();     fProfileAccumData = NULL;
    fProfileCalcIncidProb->~Profiler(); fProfileCalcIncidProb = NULL;
    fProfilePFactors->~Profiler();      fProfilePFactors = NULL;
    fProfileOptimize->~Profiler();      fProfileOptimize = NULL;
    fProfileFinish->~Profiler();        fProfileFinish = NULL;
    fProfileTotal->~Profiler();         fProfileTotal = NULL;
}

// **********************************************************************

void PEnergy::SetD( const string& param, const double value ) {
    if( param == string( "tabAreaFactor" ) ) {
        fTabAreaFactor = value;
    }else if ( param == string( "pmtEffFactor" ) ) {
        fEffScaleFactor = value;
    }else{
        throw Processor::ParamUnknown( param );
    }
}

// **********************************************************************

void PEnergy::SetI( const string& param, const int value ) {
    if( param == string( "seedActual" ) ) {
        fSeedActual = (value==1);
    }else if( param == string( "electronPrime" ) ) {
        fElectronPrime = (value==1);
    }else if( param == string( "useCharge" ) ) {
        fUseCharge = (value==1);
    }else if( param == string( "useTime" ) ) {
        fUseTime = (value==1);
    }else if( param == string( "nPhotonSim" ) ) {
        fNPhotonSim = value;
        fMakeWaterAdj = false;  // Do not apply fit adjustment if non-default
                                // value of fNPhotonSim is used in water
    }else if( param == string( "randomSeed1" ) ) {
        // higher-order digits of random number seed
        fRndmNoSeed = 1000000*value;
    }else if( param == string( "randomSeed2" ) ) {
        // lower-order digits of random number seed
        fRndmNoSeed = fRndmNoSeed + value;
    }else if( param == string( "useDirection" ) ) {
        fUseDirection = (value==1);
    }else {
        throw Processor::ParamUnknown( param );
    }
}

// **********************************************************************

void PEnergy::SetS( const string& param, const string& value ) {
    if( param == string( "chargeType" ) ) {
        if( value=="qhs" ) {
            fChargeType = 0;
        }else if( value=="qhl" ) {
            fChargeType = 1;
//      For now, not allowing QLX; there are no threshold values in
//      PMTCALIB.ratdb for QLX.
//      }else if( value=="qlx" ) {
//          fChargeType = 2;
        }else{
            throw Processor::ParamUnknown( param+" "+value );
        }
    }else{
        throw Processor::ParamUnknown( param );
    }
}

// **********************************************************************

void PEnergy::DefaultSeed() {
    // No provision for a default seed
    fValidSeed = false;
}

// **********************************************************************

void PEnergy::SetSeed( const DS::FitResult& seed ) {
    fFitResult = seed;
    fValidSeed = true;
    bool haveSeedVertex = true;
    try {
        fFitVertex = fFitResult.GetVertex(0);
    } catch( DS::FitResult::NoVertexError& error ) {
        haveSeedVertex = false;
        fValidSeed = false;
        if(!fSeedActual) warn <<
                        "PEnergy::SetSeed: no seed vertex provided. "  << endl;
    }

    fSeedEnergy = 2.5; // Defaults if valid energy not provided by seed
    if( fAVMaterialName=="lightwater_sno")  fSeedEnergy = 7.;
    fEnergyError = 2.;

    if(haveSeedVertex) {
        if( fFitVertex.ContainsEnergy() && fFitVertex.ValidEnergy() ) {
            double temp = fFitVertex.GetEnergy();
            if( !std::isnan(temp) ) fSeedEnergy = temp;
            if( fAVMaterialName=="lightwater_sno" && fSeedEnergy<0.8 )
                                                            fSeedEnergy=0.8;
            if( fFitVertex.ValidPositiveEnergyError() ) {
                fEnergyError = fFitVertex.GetPositiveEnergyError();
            }
        }

        if( fFitVertex.ContainsTime() ) {
            fSeedTime = fFitVertex.GetTime();
        }else{
            if(fUseTime) fValidSeed=false;
        }

        if( fFitVertex.ContainsPosition() && fFitVertex.ValidPosition() ) {
            fTVecEvtPosition = fFitVertex.GetPosition();
        }else{
            fValidSeed = false;
//          if(!fSeedActual) warn <<
//                         "PEnergy::SetSeed: seed position invalid. " << endl;
        }

        if( fFitVertex.ValidDirection() ) {
            fValidDirection = true;
            fTVecEvtDirection = fFitVertex.GetDirection();
        }else{
            fValidDirection = false;
            if(fUseDirection) {
                fValidSeed=false;
//              if(!fSeedActual) warn <<
//                         "PEnergy::SetSeed: seed direction invalid. " << endl;
            }
        }
    }

    if(fSeedActual) {           // seed with MC position, direction, & time
        if( !haveSeedVertex) fFitVertex = DS::FitVertex();
        fValidSeed = true;
        SetMCSeed();
        if( !haveSeedVertex) fFitResult.SetVertex( 0, fFitVertex );
    }

    fEvtPosition.set(
            fTVecEvtPosition.X(),fTVecEvtPosition.Y(),fTVecEvtPosition.Z());
    fEvtDirection.set(
            fTVecEvtDirection.X(),fTVecEvtDirection.Y(),fTVecEvtDirection.Z());

#ifdef PENERGYSEED
    if(fValidSeed) {
        detail << "PEnergy::SetSeed - Evt kE,r,pos (in cm),dir,t = "<<
            fSeedEnergy <<" "<<
            0.1*fTVecEvtPosition.Mag() <<" "<<
            0.1*fTVecEvtPosition.X() <<" "<<
            0.1*fTVecEvtPosition.Y() <<" "<<
            0.1*fTVecEvtPosition.Z() <<" "<<
            fTVecEvtDirection.X() <<" "<<
            fTVecEvtDirection.Y() <<" "<<
            fTVecEvtDirection.Z() <<" "<<
            fSeedTime << endl ;
    }
#endif

}

// **********************************************************************

void PEnergy::SetMCSeed() {
    if( !fRun->GetMCFlag() ) RAT::Log::Die(
      "\nPEnergy::SetMCSeed: seedActual can be used only with MC data");
    // Note:  fDS is defined in Method.hh
    if(fDS==NULL) RAT::Log::Die( "\nPEnergy::SetMCSeed: fDS is NULL" );
    RAT::DS::MC mc = fDS->GetMC();
    RAT::DS::MCEV mcev = fDS->GetMCEV(0);
    if( mc.GetMCParticleCount()==0 )
                RAT::Log::Die( "\nPEnergy::SetMCSeed: no particles" );
    RAT::DS::MCParticle mcpart = mc.GetMCParticle(0);

    // PMT times are recorded relative to an origin that is equal to the global
    // trigger time + gtriggerdelay(110 ns) -500 ns.  For a MC event, GetGTTime
    // returns the time of the global trigger relative to the MC event time.
    // Here we set the seed time equal to the MC event time relative to the
    // same origin as the PMT times.
    DBLinkPtr fLdaq = DB::Get()->GetLink( "DAQ" );
    fSeedTime = 500. - mcev.GetGTTime() - fLdaq->GetD("gtriggerdelay");
    fFitVertex.SetTime(fSeedTime, true, true);

    fTVecEvtPosition  = mcpart.GetPosition();
    fFitVertex.SetPosition(fTVecEvtPosition, true, true);
    fFitVertex.SetPositionErrors(
                    ROOT::Math::XYZVectorF(0.,0.,0.), true, true);

    fTVecEvtDirection  = (mcpart.GetMomentum()).Unit();
    fFitVertex.SetDirection(fTVecEvtDirection, true, true);
    fFitVertex.SetDirectionErrors( TVector3(0.,0.,0.), true, true);
}

// **********************************************************************

DS::FitResult PEnergy::GetBestFit() {
#ifdef PENERGYDBG
    detail << "Entering PEnergy::GetBestFit"<< endl;
#endif
    fCallCounter++;
    if( !fValidSeed) {
        detail<< "PEnergy::GetBestFit -- seed invalid, no fit attempted"<<endl;
        fFitVertex.SetEnergyValid(false);
        return fFitResult;
    }
    if( fTVecEvtPosition.Mag() > 8000.) {
        detail<<
           "PEnergy::GetBestFit -- Seed position at r>800 cm; no fit attempted"
           <<endl;
        fFitVertex.SetEnergyValid(false);
        return fFitResult;
    }
    // fEvent is inherited from Method class
    fProfileTotal->Start();

    fMaxPE =  GetMaxPE(fSeedEnergy);
    fMaxPE = fPMTChgpdfPtr->SetMaxPE(fMaxPE);   // fMaxPE may be reduced
                                                // in this function call

    // Check if we can avoid re-doing the auxiliary simulation if event is the
    // same one (same GTID) already fitted by another PEnergy fitter:
    Int_t gtid = fEvent->GetGTID();
    fSameGTID = (fPrevGTID==gtid);
    fSameSelector = (   fPrevSelector==fSelectorName &&
                        fPrevLowCut==fLowCut &&
                        fPrevHighCut==fHighCut );
    bool skipSim = (    fSameGTID &&
                        fPrevSeedActual==fSeedActual &&
                        fPrevElectronPrime==fElectronPrime &&
                        fPrevUseDirection==fUseDirection &&
                        fPrevNPhotonSim==fNPhotonSim &&
                        fPrevTabAreaFactor==fTabAreaFactor &&
                        fPrevEffScaleFactor==fEffScaleFactor &&
                        fSameSelector );
    fPrevGTID=gtid ;
    fPrevSeedActual=fSeedActual ;
    fPrevElectronPrime=fElectronPrime ;
    fPrevUseDirection=fUseDirection ;
    fPrevNPhotonSim=fNPhotonSim ;
    fPrevTabAreaFactor=fTabAreaFactor ;
    fPrevEffScaleFactor=fEffScaleFactor ;
    fPrevSelector=fSelectorName;
    fPrevLowCut=fLowCut;
    fPrevHighCut=fHighCut;

    // Get list of selected PMTs (all PMTs if the NULL selector is used).
    // fSelectedPMTData is inherited from SelectorMethod class.
    // fSelectedPMTData is defined as "std::vector<FitterPMT> fSelectedPMTData;"
    fSelectedPMTData = fSelector->
                        GetSelectedPMTs( fPMTData, fFitResult.GetVertex( 0 ) );
    if( fSelectedPMTData.empty() ) {
        string msg = "PEnergy: No hits to fit";
        detail << msg << endl;
        throw MethodFailError(msg);
    }
    if( fSelectedPMTData.size()<6 ) {
        string msg = "PEnergy: Fit abandoned, Nhits<6";
        detail << msg << endl;
        throw MethodFailError(msg);
    }
    for(unsigned int k=fFirstNormal; k<=fLastNormal; k++) {
        fWasHit[k] = false;
        fXIdx[k] = -1;
        fUsePMT[k] = fOnLine[k];
        if(fOnLine[k] && fSelectorName=="timeResidualCut" ) {
             fPhotonTransitTime[k] =
                     fTRCut->CalcTransitTime(fTVecEvtPosition,fPosPMT[k]);
        }else{
             fPhotonTransitTime[k] = -999.;
        }
    }

    const DU::PMTCalStatus& pmtCalStatus =DU::Utility::Get()->GetPMTCalStatus();
    for(size_t iSel=0; iSel<fSelectedPMTData.size(); iSel++) {
        int idx = fSelectedPMTData[iSel].GetID();
        fXIdx[idx] = iSel;
        fWasHit[ idx ] = true;
        // Use a hit only if it is has valid ECA, PCA, and is not cross-talk:
        fUsePMT[idx] &= (pmtCalStatus.GetHitStatus(fSelectedPMTData[iSel])==0);

//      double relHitTime =fSelectedPMTData[iSel].GetTime()- fSeedTime;
//      detail<<"GetBestFit: iSel, relHitTime = "<<iSel <<" "<<relHitTime<<endl;
    }

#ifdef PENERGYDETAIL
        int nUsePMT = 0;
        for(unsigned int i=fFirstNormal; i<=fLastNormal; i++) {
            if (fUsePMT[i]) nUsePMT++;
        }
        info <<"PEnergy: number of tubes used in fit = " << nUsePMT << endl;
#endif

    // Run auxiliary Geant simulation
    bool goodAuxSim = fGoodPrevAuxSim;
    if(!skipSim) goodAuxSim = DoAuxiliarySimulation();
    fGoodPrevAuxSim = goodAuxSim;
    if (!goodAuxSim) {
        info << "PEnergy:  failed auxiliary simulation" << endl;
        fFitVertex.SetEnergyValid(false);
        return fFitResult;
    }

    fProfileOptimize->Start();
    fOptimiser->SetComponent( this );
    // Note: fOptimiser is defined in OptimisedComponent.hh.
    // Call the designated optimiser to fit for energy and save the fom
    const double logL = fOptimiser->Maximise();
    fFitResult.SetFOM( "ELogL", logL );
    fProfileOptimize->Stop();
    fProfileFinish->Start();

    // Now save the optimised values
    vector<double> params =  fOptimiser->GetParams();
    if(fMakeWaterAdj) params[0] = WaterAdjustment( params[0] );
    SetParams( params );
    SetPositiveErrors( fOptimiser->GetPositiveErrors() );
    SetNegativeErrors( fOptimiser->GetNegativeErrors() );
    // Save the number of hits used in the fit as a FOM
    fFitResult.SetFOM( "fitNHits", (double)fSelectedPMTData.size() );
#ifdef PENERGYFITOUT
    std::stringstream strFitName;
    strFitName << "PEnergy" ;
    if (fSelectorName == "null") {
        strFitName << "Null     ";
    }else if(fSelectorName == "timeResidualCut") {
        strFitName << "TRC" << fHighCut << "      ";
    }else if(fSelectorName == "modeCut") {
        strFitName << "Mode" << fHighCut << "      ";
    }
    int chgSetting=0;
    if(fUseCharge) chgSetting = fChargeType +1;
    info << strFitName.str().substr(0,15) <<
        fSeedActual <<
        chgSetting<<
        fUseTime <<
        fUseDirection <<" "<<
//         fMakeWaterAdj <<" "<<
        fNPhotonSim/1000 <<" "<<
        "Hits, seed, result: ";

    char outbuf[30];
    int nHits = fSelectedPMTData.size();
    sprintf(outbuf,"%4d %6.3f %6.3f +/- %6.3f",
        nHits, fSeedEnergy, params[0],
        (fOptimiser->GetPositiveErrors())[0] );
    info << outbuf << endl;

//     PrintLkhdFunction( params[0] );
#endif

    fProfileFinish->Stop();
    fProfileTotal->Stop();
#ifdef PENERGYDBG
    detail << "Leaving PEnergy::GetBestFit" << endl;
#endif
    return fFitResult;
}       // PEnergy::GetBestFit

// **********************************************************************

unsigned int PEnergy::GetMaxPE( const double estKE) {
    // Adjust number of PEs considered, to save execution time at
    // lower energy levels
    double maxPE;
    if( fLightYield == 0.) {
        // Material in AV assumed to be water
        if(estKE>=40.) {
            maxPE =  8. + (15.-8.)*(estKE-40.)/(80.-40.) ;
        }else if(estKE>=20.){
            maxPE =  6. + ( 8.-6.)*(estKE-20.)/(40.-20.) ;
        }else if(estKE>=10.){
            maxPE =  5. + ( 6.-5.)*(estKE-10.)/(20.-10.) ;
        }else if(estKE>= 5.){
            maxPE =  4. + ( 5.-4.)*(estKE- 5.)/(10.- 5.) ;
        }else if(estKE>= 3.){
            maxPE =  3. + ( 4.-3.)*(estKE- 3.)/( 5.- 3.) ;
        }else{
            maxPE = 3.;
        }
    }else{
        // Scintillator
        maxPE = 3.+ estKE*fLightYield/11900.;
    }
    return (int)round(maxPE);
}
// **********************************************************************

double PEnergy::WaterAdjustment( const double fitE) {
    // Make adjustment to compensate for energy-dependent fit bias in water.
    // Note:  The values of f10 and f80 used here are valid only for 200,000
    // photons in the auxiliary simulation, which is the default for water.
    double adjVal;
    double f10 = 1.02925;
    if(fitE < f10*10.) {
        adjVal = fitE/f10;
        return adjVal;
    }
    double f80 = 0.983533;
    double slope = (f80 - f10)/(log10(80.) - 1.);
    double lastVal = 0.01;
    adjVal = fitE;
    while ( (adjVal-lastVal)>0.001 ) {  // solve by successive approximation
        lastVal = adjVal;
        adjVal = fitE/(f10 + slope*(log10(adjVal) - 1.));
    }
    return adjVal;
}

// **********************************************************************

bool PEnergy::DoAuxiliarySimulation() {
#ifdef PENERGYDBG
    detail << "Entering PEnergy::DoAuxSim" << endl;
#endif
    G4RunManager* g4RunManager = G4RunManager::GetRunManager();
    G4EventManager* evtManager = G4EventManager::GetEventManager();
    G4TrackingManager* fTrackingManager = evtManager->GetTrackingManager();

    // Save the run setting of the "StoreTrajectory" control variable
    G4int fMainStoreTrajectory = fTrackingManager->GetStoreTrajectory();

    // Save the photon thinning factor used in the main simulation
    double mainThinningFactor = PhotonThinning::GetFactor();
    PhotonThinning::SetFactor( 1. );

    // Save the state of the main random number sequence and restore to the
    // state of the local (PEnergy) sequence:
    CLHEP::HepRandom* hepRandom = CLHEP::HepRandom::getTheGenerator();
    fMainSeqState.clear(); fMainSeqState.seekp(0); fMainSeqState.seekg(0);
    hepRandom->saveStaticRandomStates(fMainSeqState);
    hepRandom->restoreStaticRandomStates(fLocalSeqState);

    // Store existing configurations for later restore
    const G4UserTrackingAction* mainTrackingAction =
                                    g4RunManager->GetUserTrackingAction();
    const G4UserSteppingAction* mainSteppingAction =
                                    g4RunManager->GetUserSteppingAction();
    TrackingAction* fTrackingAction = new TrackingAction();
    SteppingAction* fSteppingAction = new SteppingAction(this);
    g4RunManager->SetUserAction( fTrackingAction );
    g4RunManager->SetUserAction( fSteppingAction );

    // Turn off storing of track trajectories to save time:
    fTrackingManager->SetStoreTrajectory(0);

    G4StateManager* gStateMgr = G4StateManager::GetStateManager();
    // Save the current Geant state
    G4ApplicationState mainAppState = gStateMgr->GetCurrentState();
    // Set flag indicating that an auxiliary simulation is being done.
    // The processing in various Gsim functions (BeginOfRunAction,
    // BeginOfEventAction, EndOfEventAction) is skipped for an auxiliary
    // simulation.
    fGsimPtr->SetAuxiliarySim(true);
    AuxSimInfo* asiPtr =  AuxSimInfo::GetAuxSimInfo();
    asiPtr->SetSimType(true);
    // If this is a MC run, Geant will be in the event processing (EventProc)
    // state
    bool isMCRun = (mainAppState==G4State_EventProc);
    if(!isMCRun) {
        // If this is the first event in a run to process events from a zdab or
        // event root file rather than newly generated MC events, then Geant
        // will be in the G4State_Idle state and the Geant run will not have
        // been initialized, which needs to be done before simulating any
        // PEnergy auxiliary events.  However, the setting of fAuxiliarySim in
        // Gsim causes Gsim::BeginOfRunAction not to be executed in the call
        // from the G4RunManager.
        if(mainAppState==G4State_Idle){
            g4RunManager->RunInitialization();
         }
    }else{
        // If it is a MC run, set the state to GeomClosed
        if( !gStateMgr->SetNewState(G4State_GeomClosed) )
                            RAT::Log::Die("PEnergy: state change A failed");
    }

    // Save existing values from the main simulation so they can be restored
    // later, and initialize for the auxiliary simulation:
    unsigned int mainCherPhotCount  = Cerenkov::GetPhotonCount();
    unsigned int mainScintPhotCount = GLG4Scint::GetScintillatedCount();
    unsigned int mainReemPhotCount  = GLG4Scint::GetReemittedCount();
    G4double mainTotEdep, mainTotEdepQuenched, mainTotEdepTime;
    G4ThreeVector  mainScintCentroidSum;
    GLG4Scint::GetTotEdepAll( mainTotEdep, mainTotEdepQuenched,
                                mainTotEdepTime, mainScintCentroidSum );

    bool haveGoodAuxSim = false;
    int trycount = 0;
    // This loop causes the auxiliary simulation to be
    // attempted up to three times.
    while( !haveGoodAuxSim && trycount<3) {
        asiPtr->SetSimEvtStatus(true);
        trycount++;

        Cerenkov::ResetPhotonCount();
        GLG4Scint::ResetPhotonCount();
        GLG4Scint::ResetTotEdep( );

        G4PrimaryVertex* fPrimaryVertex = new G4PrimaryVertex(fEvtPosition,0.);
        G4Event* auxEvent = new G4Event();
        asiPtr->SetEventPointer(auxEvent);
        auxEvent->AddPrimaryVertex(fPrimaryVertex);
        if(fElectronPrime) {
            fNPrimaries = 1 + (int)(fNPhotonSim/ExpectedPhotons(fSeedEnergy));
#ifdef PENERGYDETAIL
            detail << "PEnergy::DoAuxiliarySimulation -- fNPrimaries, "<<
                "fValidDirection = "<< fNPrimaries<<" "<< fValidDirection<<endl;
            detail << "total expected photons in aux simulation = " <<
                            fNPrimaries*ExpectedPhotons(fSeedEnergy) << endl;
#endif
        }else{
            // Probably should allow the fElectronPrime=false  option only if
            // fUseTime=false.
            fNPrimaries = fNPhotonSim;
        }

        //  Create primary particles
        G4PrimaryParticle* primary;
        for(unsigned int i=0; i<fNPrimaries; i++) {
            G4ThreeVector momentum;
            if(fUseDirection && fValidDirection) {
                momentum = fEvtDirection;
            }else{
                // Generate direction vector from isotropic distribution:
                double cosTheta = -1. + 2.*G4UniformRand();
                double phi = 2.*pi* G4UniformRand();
                double sinTheta = sqrt( 1.-cosTheta*cosTheta);
                momentum = G4ThreeVector( sinTheta*cos(phi),
                                      sinTheta*sin(phi), cosTheta);
            }

            if(fElectronPrime) {            // generate electrons
                double E = fSeedEnergy + electron_mass_c2;
                momentum=sqrt(E*E - electron_mass_c2*electron_mass_c2)*momentum;
                primary = new G4PrimaryParticle( G4Electron::Definition(),
                      momentum.x(), momentum.y(), momentum.z(), E);
            }else{                          // generate photons
                double E = 3.35*electronvolt;
                momentum = E*momentum;
                primary = new G4PrimaryParticle( G4OpticalPhoton::Definition(),
                      momentum.x(), momentum.y(), momentum.z(), E);
                // Generate random polarization:
                double phi = twopi* G4UniformRand();
                G4ThreeVector e1 = momentum.orthogonal().unit();
                G4ThreeVector e2 = momentum.unit().cross(e1);
                G4ThreeVector pol = e1*cos(phi)+e2*sin(phi);
                primary->SetPolarization(pol);
            }
            fPrimaryVertex->SetPrimary(primary);
        }
        // Initialize arrays to accumulate auxiliary simulation data:
        CalcInit();

        // Simulate the auxiliary event and accumulate data:
        fProfileGenPhotons->Start();
        evtManager->ProcessOneEvent(auxEvent);
        fProfileGenPhotons->Stop();
        haveGoodAuxSim = asiPtr->IsSimEvtGood();
        if( haveGoodAuxSim ) {
            // G4ApplicationState.hh   defines the states

            // For each PMT, calculate the probability that any single photon
            // is incident on it, and the average probability that an incident
            // photon produces a signal in the PMT:
            CalcIncidenceProbs();

#ifdef PENERGYDETAIL
            // Extract some stuff from the Geant event for display:
            G4PrimaryVertex* pv = auxEvent->GetPrimaryVertex();
            double x = pv->GetX0();
            double y = pv->GetY0();
            double z = pv->GetZ0();
            double t = pv->GetT0();
            double r = sqrt(x*x + y*y + z*z);
            detail<<"DoAuxiliarySimulation - Geant auxiliary event data (mm):" ;
            detail << "  Vertex x, y, z, r, t = "<< x <<" "<< y <<" "<< z
                <<" "<< r <<" "<< t << endl;
#endif
        }
        // Destruction of the primary vertex and particles is taken care of
        // when the event is destroyed.
        delete auxEvent;  auxEvent = NULL;
        if( !haveGoodAuxSim ) {
            detail<<"PEnergy:  Failed auxiliary simulation attempt "
                                                            << trycount <<endl;
        }
    }

    if( !haveGoodAuxSim ) {
        detail<<"PEnergy:  Multiple failures of auxiliary simulation."<<endl;
        detail<<"          Fit abandoned"<<endl;
    }

    fGsimPtr->SetAuxiliarySim(false);
    asiPtr->SetSimType(false);

    // Restore the main setting of StoreTrajectory:
    fTrackingManager->SetStoreTrajectory(fMainStoreTrajectory);
    // Restore the photon thinning factor:
    PhotonThinning::SetFactor( mainThinningFactor );

    // Save the state of the local random number sequence and restore the main
    // sequence:
    fLocalSeqState.clear(); fLocalSeqState.seekp(0); fLocalSeqState.seekg(0);
    hepRandom->saveStaticRandomStates(fLocalSeqState);
    hepRandom->restoreStaticRandomStates(fMainSeqState);

    // If it's a MC run, restore the state of the primary simulation.
    // Otherwise leave Geant in the GeomClosed state so that RunInitialization
    // is not done for any event except the first in a run.
    if(isMCRun){
        if( !gStateMgr->SetNewState(mainAppState) )
                            RAT::Log::Die("PEnergy: state change B failed");
    }

    // Restore the main TrackingAction and SteppingAction objects
    if(fTrackingAction!=NULL) delete fTrackingAction;
    fTrackingAction = NULL;
    if(fSteppingAction!=NULL) delete fSteppingAction;
    fSteppingAction = NULL;
    g4RunManager->
        SetUserAction(const_cast<G4UserTrackingAction*>(mainTrackingAction));
    g4RunManager->
        SetUserAction(const_cast<G4UserSteppingAction*>(mainSteppingAction));

#ifdef PENERGYDETAIL
    if (fElectronPrime) {
        detail << "mainScintPhotCount, mainReemPhotCount = " <<
                            GLG4Scint::GetScintillatedCount() <<" "<<
                            GLG4Scint::GetReemittedCount() << endl;
    }
#endif
    Cerenkov::ResetPhotonCount( mainCherPhotCount );
    GLG4Scint::ResetPhotonCount( mainScintPhotCount, mainReemPhotCount);
    GLG4Scint::SetTotEdep( mainTotEdep, mainTotEdepQuenched,
                                    mainTotEdepTime, mainScintCentroidSum );
#ifdef PENERGYDBG
    detail << "Leaving PEnergy::DoAuxSim" << endl;
#endif
    return haveGoodAuxSim;
}       // PEnergy::DoAuxiliarySimulation

// **********************************************************************

void PEnergy::CalcInit() {
#ifdef PENERGYDBG
    detail << "Entering PEnergy::CalcInit" << endl;
#endif
  // Estimated minimum photon propagation time -- i.e., the time it takes
  // photons by the shortest path to reach the PSUP from the event location.
  // Note:  There could be a problem here for Cherenkov events, which do
  // not have an isotropic photon distribution wrt. the event location.
  fEstMinPropTime= ( fAVInnerRadius - fTVecEvtPosition.Mag() )*
                                        fRefIdxAVFill/c_light + fPropTimeAdder;
  // Triggers are set on ticks of the 50 MHz clock.
  // In following, 10 ns is 1/2 the interval between ticks.
  // fMaxInclusionTime - estimated latest time for a hit that would be built,
  //                     i.e., included in the event.  Used to test if an
  //                     incident photon in auxiliary simulation should be used.

  fMaxInclusionTime = fEstMinPropTime + fTotalTriggerDelay + 10.;
  if(fSelectorName == "modeCut") fMaxInclusionTime = fEstMinPropTime + fHighCut;
  fPMTTimeDistnsPtr->SetEventParameters(
          fEstMinPropTime + fFECToTriggerDelay +10. , fLowCut, fHighCut );

  for(unsigned int k=fFirstNormal; k<=fLastNormal; k++) {
      // Initially use fPIP[k] to sum the number of photons incident on the
      // PMT's tabulation area and fPMTEff[k] to sum the PMTResponse values for
      // the incident photons.
      fPIP[k] = 0.;
      fPMTEff[k] = 0.;
  }

  // Array for saving times at which photons intersect each PMT's tabulation
  // area:
  fTCross.resize(fNPMTsDim);
  fNReachPSUP = 0 ;  // to count number of photons that reach the PSUP
  fNTabRegions = 0 ; // to count number of tabulation areas incremented
#ifdef PENERGYDBG
  detail << "Leaving PEnergy::CalcInit" << endl;
#endif
}

// **********************************************************************

void PEnergy::AccumulateData(const G4Step* step)  {
//     fProfileAccumData->Start();
    if( step->GetTrack()->GetParticleDefinition()->GetParticleName()
                                                != "opticalphoton" ) return;

    G4StepPoint* postStepPt = step->GetPostStepPoint();
    G4ThreeVector  postPosition = postStepPt->GetPosition();
    double   ptTime   = postStepPt->GetGlobalTime();   // Time at current point

    TVector3 ptPos;        // Current point on the track
    for(int i=0; i<=2; i++) {
        ptPos(i) = postPosition(i);
    }

    // In the following logic, we set checkIncdt=true if ptPos is outside the
    // test radius and the photon track is directed outward at that point.
    // Also require that the time at the point is within the interval
    // for a hit generated by the photon to be included in the event.
    if( !( ptPos.Mag2()>fRsqTest && ptTime<fMaxInclusionTime) ) {
//         fProfileAccumData->Stop();
        return;
    }

    G4StepPoint* preStepPt  = step->GetPreStepPoint();
    double prevTime = preStepPt->GetGlobalTime();   // Time at previous point
    G4ThreeVector g4Polarization = postStepPt->GetPolarization();
    double eKin = postStepPt->GetKineticEnergy();
    G4ThreeVector prePosition = preStepPt->GetPosition();
    TVector3 prevPos, polarization;      // Previous point on the trajectory
    for(int i=0; i<=2; i++) {
        prevPos(i) = prePosition(i);
        polarization(i) = g4Polarization(i);
    }

    double rsqPrev = prevPos.Mag2();

    bool checkIncdt = false;
    if( rsqPrev<=fRsqTest ) {  // Is prevPos inside test radius?
        checkIncdt = true;
    }else{              // prevPos outside test radius
        // Check if photon is travelling inward at point prevPos
        TVector3 segmentVect = ptPos - prevPos;
        double dotprod = segmentVect.Dot(prevPos);
        if ( dotprod < 0.) {
            double rho = -dotprod/segmentVect.Mag2();
            // rho is the fraction of the distance between prevPos and
            // ptPos at which the line segment joining those two points
            // is closest to the center of the detector.

            // Is location of the closest point to the detector center
            // between the previous and this position on the segment?
            // If so, set checkIncdt true and reset prevPos & prevTime
            // to be the position & time of the closest point.
            if (rho>=0. && rho<1.) {
                prevPos = prevPos + rho*segmentVect;
                rsqPrev = prevPos.Mag2();
                prevTime = prevTime + rho*(ptTime-prevTime);
                checkIncdt = (rsqPrev<=fRsqTest) ;
            }
        }
    }

    if( checkIncdt) {   // Now check on which (if any) PMT tabulation
                        // areas the photon is incident
        fNReachPSUP++;
        TVector3 photonDirection = (ptPos - prevPos).Unit();

        // Point at which photon crosses test radius:
        TVector3 crossPt = SphereCrossPoint(fRTest, prevPos, ptPos);
        TVector3 crossUnit = crossPt.Unit();
        double crossUX = crossUnit.X();
        double crossUY = crossUnit.Y();
        double crossUZ = crossUnit.Z();

        for(unsigned int k=fFirstNormal; k<=fLastNormal; k++){ // Loop over PMTs
            if( fUsePMT[k]) {
                if( (fPosPMTUX[k]*crossUX + fPosPMTUY[k]*crossUY +
                        fPosPMTUZ[k]*crossUZ) > fTabLim ) {
                    // Calculate where photon crosses PMT face or
                    // tabulation area.  d is the distance along the photon
                    // track from the previous photon position to the
                    // crossing point.
                    double cosIncdAngle=photonDirection.Dot(fDirPMT[k]);
                    if(cosIncdAngle>1.) cosIncdAngle = 1.;
                    double d= (fPosPMT[k]-prevPos).Dot(fDirPMT[k])/cosIncdAngle;
                    // Crossing point = prevPos + d*photonDirection;
                    // Square of distance from PMT position to where photon
                    // intersects (crosses) the "PMT face":
                    double rhoSq=(prevPos+d*photonDirection -fPosPMT[k]).Mag2();

                    // In rare cases there may be a photon that, although
                    // directed outward, is almost perpendicular to a vector
                    // from the AV center, and a PMT that is not pointed
                    // directly inward such that the photon is moving away from
                    // the PMT face.  We eliminate those cases by requiring
                    // that the cosine of the incident angle be positive.
                    if(rhoSq<fRhoSqLim && cosIncdAngle>0.) {
                       // Time the photon intersects the tabulation area:
                        double t = prevTime +
                                (ptTime-prevTime)*d/(ptPos-prevPos).Mag();
                        double weight = 1.;
                        if(fSelectorName=="timeResidualCut") {
                            // tResid is time residual for photon arrival at
                            // PMT.  Total time residual is photon time
                            // residual plus PMT transit time.  PMT transit
                            // time distribution is apparently offset by amount
                            // such that the prompt peak has nominally a zero
                            // PMT transit time.
                            double tResid = t - fPhotonTransitTime[k];
                            // Weight the photon by the probability that the
                            // time of a hit generated by it (i.e., tResid + PMT
                            // transit time) would be within the time residual
                            // cut acceptance interval:
                            weight = fPMTTransitTimePtr->GetIntervalProbability
                                     ( fLowCut-tResid, fHighCut-tResid, false );
//                             detail<<"times: "<<t<<" "<< fPMTTransitTime[k]
//                                     <<" "<<tResid<<" "<<acceptPhoton<< endl;
                        }else if(fSelectorName=="modeCut") {
                            // tTest is the time of the photon arrival relative
                            // to the estimated time of peak arrivals (the
                            // mode).
                            double tTest = t - fEstMinPropTime;
                            weight = fPMTTransitTimePtr->GetIntervalProbability
                                     ( fLowCut-tTest, fHighCut-tTest, false );
                        }
                        if(weight>0.) {
                            fTCross[k].push_back(t);
                            fNTabRegions++;
                            fPIP[k] += weight;
                            double incdAngle = acos(cosIncdAngle);
                            // vector perpendicular to the plane of incidence:
                            TVector3 senkrecht =
                                    (fDirPMT[k].Cross(photonDirection)).Unit();
                            double cosPol = senkrecht.Dot(polarization);
                            fPMTEff[k] +=
                                    weight*PMTResponse(eKin, incdAngle, cosPol);
                        }
                    }   // if(rhoSq<fRhoSqLim)
                }
            }    // if( fUsePMT[k])
        }        // Loop over PMTs
    }            // if( checkIncdt)
//     fProfileAccumData->Stop();
}

// **********************************************************************

void PEnergy::CalcIncidenceProbs() {
#ifdef PENERGYDBG
    detail << "Entering PEnergy::CalcIncidenceProbs" << endl;
#endif
// For each PMT, calculate the probability that any single photon is
// incident on it, and the average probability that a photon produces a
// hit in the PMT, including the electronics channel

    fProfileCalcIncidProb->Start();
    double adjFactor;
    if(fElectronPrime) {
        adjFactor= 1./ (fTabAreaFactor*(double)
            ( Cerenkov::GetPhotonCount() + GLG4Scint::GetScintillatedCount() ));
            // Above line is sum of photons created in the auxiliary simulation
    }else{
        adjFactor = 1./( fTabAreaFactor*(double)fNPhotonSim );
    }

    double noiseMult = fTabAreaFactor*fNPrimaries;
    for(unsigned int k=fFirstNormal; k<=fLastNormal; k++) {
        if( fUsePMT[k] ) {
            // Note: At this point, fPIP[k] is the number of photons incident
            // on the tabulation area of PMT k in the auxiliary simulation, and
            // fPMTEff[k] is the sum of the PMTResponse values for the incident
            // photons.
            if( fPIP[k]>0. ) {
                // Calculate average PMT response and multiply by various other
                // factors
                fPMTEff[k] = fEffScaleFactor*fCollectionEfficiency*fPMTEff[k]
                               *PMTVariation::Get()->GetVariation(k)/fPIP[k];
            }else{        // No simulated photon incident on PMT; use escape
                          // probability for wavelength of 370 nm.
                fPMTEff[k] = fEffScaleFactor*fCollectionEfficiency
                    *fPEEscapeProbability->Value (3.35*electronvolt)
                                *PMTVariation::Get()->GetVariation(k);
            }

            // Now fPIP[k] becomes the probability any single photon is incident
            // on PMT k.
            fPIP[k] = adjFactor*fPIP[k];
            // What probability to assign if no photons were incident in the
            // auxiliary simulation?  Can't use zero.  Geometric PMT coverage
            // is about 70%, and there are about 10,000 PMTs.  If photons were
            // isotropically distributed and there was no attenuation, would
            // have incidence probability of about 0.7/10000 = 0.00007.  Use
            // value equal to half of the smallest possible non-zero value, but
            // no greater than 0.00001
            if( fPIP[k]==0. ) {
                fPIP[k]= 0.5*adjFactor;
                if(fPIP[k]>0.00001) fPIP[k]= 0.00001;
            }
        }else{
            fPIP[k] = 0.;
        }         // if( fUsePMT[k] )
    }             //  for(int k=fFirstNormal; k<=fLastNormal; k++)
    fProfileCalcIncidProb->Stop();
    fProfilePFactors->Start();

    // Test if we can avoid redoing the fPQFac computations
    bool skipPFacComp = ( fSameGTID &&
                                fPrevChargeType==fChargeType && fSameSelector);
    const DU::ChanHWStatus& chanHWStatus= DU::Utility::Get()->GetChanHWStatus();
    if ( fUseCharge && !skipPFacComp ) {
        // Warning:  the following clear() statement is essential.  Otherwise
        // the compiled code does not expand the size of the inner vectors even
        // if the new value of fMaxPE is larger than the value for the previous
        // event.
        fPQFac.clear();
        fPQFac.resize( fNPMTsDim, vector<double>(fMaxPE,0.) );

        fPrevChargeType = fChargeType;
        // Calculate fPQFac[][].  For hit PMT i, fPQFac[i][j] is the probability
        // of a hit with the observed charge given j+1 photoelectrons (i.e.,
        // the probability that the charge is above the threshold times the
        // probability of the observed charge given that the charge is above
        // the threshold).  For a PMT i that was not hit, fPQFac[i][j] is the
        // probability of no hit given j+1 photoelectrons.

        // Calculate probability factors for PMTs that were hit:
        double nExpct =  ExpectedPhotons(fSeedEnergy);
        double timePdfPtVal, timeCdfPtVal;
        size_t nHitPMTs= fSelectedPMTData.size();

        HitTimeProbabilityFactor htpf(fMaxPE);
        if(fUseTime) {
            fPTFac.clear();
            fPTFac.resize( nHitPMTs, vector<double>(fMaxPE,0.) );
        }

        for(size_t i=0; i<nHitPMTs; i++) {
            UInt_t hitPMTID = fSelectedPMTData[i].GetID();
            if( hitPMTID>=fNPMTsDim ) info <<
                        "PEnergy::CalcIncidenceProbs -- dimension error: "
                        << hitPMTID <<" "<< fNPMTsDim << endl;
            double hhp = fHHP[ hitPMTID ];
            // Charge threshold:
            double threshold= chanHWStatus.GetThresholdAboveNoise(hitPMTID);
            if( fUsePMT[hitPMTID] && threshold<0.5*hhp) {  // Note: PMTs
                // with threshold>=0.5*hhp  are ignored in function operator().
                double pedSig = fPedSig[ fSelectedPMTData[i].GetCCCC() ];
                double q;
                if(fChargeType==0){
                    q = fSelectedPMTData[i].GetQHS();
                }else if(fChargeType==1){
                    q = fSelectedPMTData[i].GetQHL();
                }else{
                    q = fSelectedPMTData[i].GetQLX();
                }

                if(fUseTime ) {
                    double relHitTime =fSelectedPMTData[i].GetTime()- fSeedTime;
                    fPMTTimeDistnsPtr->GetRelHitTimeProbs(
                        noiseMult*fNoiseWindow*fExpNoise[hitPMTID],relHitTime,
                        fTCross[hitPMTID], timePdfPtVal, timeCdfPtVal);
                    htpf.SetChargeParameters(threshold, hhp, q);
                    htpf.SetTimeParameters( timePdfPtVal, timeCdfPtVal);
                }
                // Calculate the expected number of photoelectrons produced
                // in the PMT by incident photons, at the seed energy:
                double lambda = nExpct*fPIP[hitPMTID]*fPMTEff[hitPMTID];

                double hFac, qFac;
                double probPE = exp(-lambda);    // probability of zero PEs
                double probMax = 0.;
                for(size_t  j=0; j<fMaxPE; j++) {
                    int jP1 = j+1;
                    // hit probability factor (probability total charge is
                    // above the threshold):
                    hFac =  1. -
                       fPMTChgpdfPtr->ProbBelowThreshold(jP1,hhp,threshold);
                    // charge probability factor (probability of observed
                    // charge given that charge is above the threshold)):
                    qFac = fPMTChgpdfPtr->
                               PdfRecordedVal(jP1,q,hhp, threshold, pedSig);
                    if(hitPMTID>=fNPMTsDim || j>=fMaxPE)
                        warn << "fPQFac index outside of range: " <<
                        hitPMTID<<" "<<fNPMTsDim <<" "<< j<<" "<<fMaxPE<<endl;

                    fPQFac[hitPMTID][j] = hFac*qFac;
                    // Following needed to avoid prob=0 in operator()
                    // function, which produces an error in log(prob).
                    if(fPQFac[hitPMTID][j]== 0.) fPQFac[hitPMTID][j]=1.e-33;
                    if(fUseTime ) {
                        // probPE = probability of (j+1) photoelectrons
                        probPE = probPE*lambda/(double)(j+1);
                        double probFac = probPE*fPQFac[hitPMTID][j] ;
                        if(probFac>probMax) probMax = probFac;
                        if(probFac > 0.001*probMax) {
                            fPTFac[i][j] = htpf.GetTimeProbabilityFactor( jP1 );
                        }else{
                            if(j>0) {
                                fPTFac[i][j] = fPTFac[i][j-1];
                            }else{
                                fPTFac[i][j] = 1.e-33;
                            }
                        }
                        if(fPTFac[i][j]<= 0.) fPTFac[i][j]=1.e-33;
                    }
                }
            }   //  if( fUsePMT[hitPMTID] && threshold<0.5*hhp)
        }       // for(size_t i=0; i<fSelectedPMTData.size(); i++)

        // Calculate probability factors for PMTs that were not hit
        for(unsigned int i=fFirstNormal; i<=fLastNormal; i++) {
            if(fUsePMT[i] && !fWasHit[i] ) {
                double hhp = fHHP[i];
                double threshold = chanHWStatus.GetThresholdAboveNoise(i) ;
                if (threshold<0.5*hhp) {
                    for(unsigned int j=0; j<fMaxPE; j++) {
                        fPQFac[i][j] =
                           fPMTChgpdfPtr->ProbBelowThreshold(j+1,hhp,threshold);
                    }
                }
            }
        }       // for(int i=fFirstNormal; i<=fLastNormal; i++)

    }    // if ( fUseCharge && !skipPFacComp )

  fProfilePFactors->Stop();

#ifdef PENERGYDETAIL
  detail<<"number of photons reaching the PSUP = "<< fNReachPSUP << endl;
  detail<<"number of increments to tabulation area sums  "<<fNTabRegions <<endl;
  detail<<"Average tabulation areas intersected per photon = " <<
                              (double)fNTabRegions/(double)fNReachPSUP << endl;
#endif
#ifdef PENERGYDBG
    detail << "Leaving PEnergy::CalcIncidenceProbs" << endl;
#endif
}
// **********************************************************************
double PEnergy::PMTResponse
    (const double photonE, const double incAngle, const double cosPol) const {
    // Get PMT response from the table, given photon energy, incidence angle,
    // and cosine of polarization angle
    // photonE = photon energy in MeV
    // incAngle = photon incidence angle on concentrator face, in radians
    // cosPol = cosine of angle between polarization direction and the
    //             senkrecht direction (perpendicular to plane of incidence)

    double degAngle = radToDeg*incAngle; // incidence angle in degrees
    if(degAngle<0.) {
        degAngle = 0.;
    }
    if(degAngle>90.){
        degAngle = 90.;
    }
    double wl = twopi*hbarc/(photonE*nanometer);    // wavelength in nm
    double pmtResponse;
    if(wl<lbWL || wl>lbWL+(nWL-1)*delWL ) {
        // detail<<"PEnergy::PMTResponse: wavelength out of bounds: "<<wl<<"\n";
        pmtResponse = 1.0e-7;
    }else{
        // Bilinear interpolation using formula from Numerical Recipes:
        int ia = (int)( (wl-lbWL)/delWL );
        int ja = (int)( (degAngle-lbAng)/delAng );
        float t = (wl-lbWL-ia*delWL)/delWL;
        float u = (degAngle-lbAng-ja*delAng)/delAng;
        double pmtRespP =
                    (1.-u)*( (1.-t)*pmtRsp[0][ia][ja] + t*pmtRsp[0][ia+1][ja])
                    +u*( t*pmtRsp[0][ia+1][ja+1] + (1.-t)*pmtRsp[0][ia][ja+1] );
        double pmtRespS =
                    (1.-u)*( (1.-t)*pmtRsp[1][ia][ja] + t*pmtRsp[1][ia+1][ja])
                    +u*( t*pmtRsp[1][ia+1][ja+1] + (1.-t)*pmtRsp[1][ia][ja+1] );
        // Calculate weighted average of response in the two perpendicular
        // polarization directions relative to the plane of incidence:
        double cossq = cosPol*cosPol;
        pmtResponse = (1.-cossq)*pmtRespP + cossq*pmtRespS;
    }
return pmtResponse;
}

// **********************************************************************

double PEnergy::ExpectedPhotons( double ke) {
    // The coefficients for the expected number of photons produced as a
    // function of event energy were obtained from fits to RAT simulation
    // outputs of electrons in light water and in LAB/PPO, for which the
    // nominal scintillation yield was 10800 photons/MeV.  For other
    // scintillation cocktails new fits should be done, but for now I assume
    // that the number of Cherenkov photons does not change in other cocktails
    // and the number of scintillation photons changes in proportion to the
    // scintillation yield.

    // Expected number of Cherenkov photons:
    double nCher =
      ((fCherCoef[3]*ke +fCherCoef[2])*ke +fCherCoef[1])*ke +fCherCoef[0];
    if(nCher<0.) nCher = 0.;
    // Expected number of scintillation photons:
    double nScint =
      (((fScintCoef[3]*ke +fScintCoef[2])*ke +fScintCoef[1])*ke +fScintCoef[0])
                                                          *fLightYield/10800.;
    if(nScint<0.) nScint = 0.;
    double nExpct = nCher + nScint; // Expected total photons
    return nExpct;
}

// **********************************************************************
void PEnergy::PrintLkhdFunction(const double eCenter) {
    vector<double> params(1);
    vector<double> llSums(5);
    int iMid = 10*(int)eCenter;
    int ia = iMid-100;
    if(ia<0) ia=1;
    detail << "#                   *--- More than noise -------*" <<
                  "*------- Noise only --------* " << endl;
    detail << "#ke      Total          "
            "Hit            Not hit         Hit           Not hit" << endl;
    for(int i=ia; i<iMid+100; i++) {
        for(int j=0; j<5; j++) llSums[j]=0.;
        double ke = 0.1*(double)i;
        params[0] = ke;
        CalcLikelihood ( params, llSums );
        char outbuf[80];
        sprintf(outbuf, "%4.1f %14.7e %14.7e %14.7e %14.7e %14.7e",
                ke, llSums[0], llSums[1], llSums[2], llSums[3], llSums[4] );
        detail << outbuf << endl;
    }
}

// **********************************************************************
double PEnergy::operator()( const vector<double>& params ) {
    vector<double> llSums(5, 0.);
    return CalcLikelihood ( params, llSums );
}

// **********************************************************************
// Following function calculates the log-likelihood value given the current fit
// values in params.  Only params[0] (energy) is used.  Various contributions
// to the log(Likelihood) are returned in llSums if charge was use in the fit.
double PEnergy::CalcLikelihood( const vector<double>& params,
                                                vector<double>& llSums  ) {
    double ke = params[0];
    // Need to avoid energy values less than or equal zero in calculation below
    bool negFlag = false;
    if(ke<0.01) {
        negFlag = true;
        ke = 0.01;
    }

    double nExpct =  ExpectedPhotons(ke); // expected total photons in event
    double logLike = 0.0;

    double extraFactor = 1.0;
    if( fUseCharge && !fUseTime) extraFactor = 1.0;
    if( fUseCharge &&  fUseTime) extraFactor = 0.99851441;
    const DU::ChanHWStatus&
                chanHWStatus = DU::Utility::Get()->GetChanHWStatus();

    for(unsigned int i=fFirstNormal; i<=fLastNormal; i++) {
        if( fUsePMT[i] ) {
            // Calculate the expected number of photoelectrons
            // produced in the PMT by incident photons:
            double lambda= nExpct*fPIP[i]*extraFactor*fPMTEff[i] ;
            // Expected number of noise photoelectrons:
            double expNoisePEs = fNoiseWindow*fExpNoise[i];
            bool mostlyNoise = ( lambda<expNoisePEs );
            // Add the noise contribution:
            lambda += expNoisePEs;

            if(fUseCharge) {
                double hhp = fHHP[i];
                double threshold = chanHWStatus.GetThresholdAboveNoise( i ) ;
                if (threshold<0.5*hhp) {
                    double probPE = exp(-lambda); // probability of zero PEs
                    double sumProb ;
                    if( fWasHit[i] ) {
                        sumProb = 0.;
                        if(fUseTime) {
                            int jLim = fMaxPE; if(mostlyNoise) jLim = 1;
                            for(int j=0; j<jLim; j++) {
                                // probPE = probability of (j+1) photoelectrons
                                probPE = probPE*lambda/(double)(j+1);
                                sumProb += probPE*fPQFac[i][j]*
                                                        fPTFac[fXIdx[i]][j];
                            }
                        }else{
                            int jLim = fMaxPE; if(mostlyNoise) jLim = 1;
                            for(int j=0; j<jLim; j++) {
                                probPE = probPE*lambda/(double)(j+1);
                                sumProb += probPE*fPQFac[i][j];
                            }
                        }
                    }else{
                        sumProb = probPE;
                        int jLim = fMaxPE; if(mostlyNoise) jLim = 1;
                        for(int j=0; j<jLim; j++) {
                            probPE = probPE*lambda/(double)(j+1);
                            sumProb += probPE*fPQFac[i][j];
                        }
                    }

                    if(sumProb <= 0.) {
                        sumProb = 1.e-33;
                    }
                    double deltaLogLike = log( sumProb );
                    if(!mostlyNoise) {
                        if( fWasHit[i] ) {
                            llSums[1] += deltaLogLike;
                        }else{
                            llSums[2] += deltaLogLike;
                        }
                    }else{
                        if( fWasHit[i] ) {
                            llSums[3] += deltaLogLike;
                        }else{
                            llSums[4] += deltaLogLike;
                        }
                    }
                }       // if (threshold<0.5*hhp)

            }else{   // if(!fUseCharge)
                double pPoisson  =  exp(-lambda);
                double pHit = 0.;
                double pNoSignal = 1.;    // Initially probability of no hit
                                          // given no PEs
                double chanEff =
                            ChannelEfficiency::Get()->GetChannelEfficiency(i);
                if(chanEff==0.) chanEff = 0.00000001;
                double factor = 1. - chanEff;  // probability a single
                                               // PE does not result in a hit
                for(unsigned int j=1; j<=fMaxPE; j++) {
                    pPoisson = pPoisson*lambda/(double)j;
                    // pPoisson is now the probability that j PEs are produced
                    // in the PMT.
                    pNoSignal = pNoSignal*factor;
                    // pNoSignal now the probability of no hit given j PEs
                    pHit = pHit + pPoisson*(1. - pNoSignal);
                }
                if( fWasHit[i] ) {
                    logLike +=  log(pHit) ;
                }else{
                    logLike += log(1.-pHit);
                }
            }   // if(fUseCharge)
        }       // if( fUsePMT[i]
    }           // for(int i=fFirstNormal; i<=fLastNormal; i++)

    if(fUseCharge) {
        llSums[0] = llSums[1] + llSums[2] + llSums[3] + llSums[4] ;
        logLike = llSums[0];
    }

    // Keep us out of trouble if optimizer is trying an energy less than zero:
    if (negFlag) logLike = logLike - 1000.*(0.01-params[0]);
    return logLike;
}       // PEnergy::operator()


// **********************************************************************
vector<double> PEnergy::GetParams() const {
  vector<double> params;
  params.push_back( fSeedEnergy );
  return params;
}

// The following two functions do the same thing, but both are required
// to match with the other Fit Methods
// Note that Minuit assigns the same value to the positive and negative errors.
vector<double> PEnergy::GetPositiveErrors() const {
  vector<double> errors;
  errors.push_back( fEnergyError );
  return errors;
}

vector<double> PEnergy::GetNegativeErrors() const {
  return GetPositiveErrors();
}

void PEnergy::SetParams( const std::vector<double>& params ) {
  DS::FitVertex vertex = fFitResult.GetVertex(0);
  vertex.SetEnergy( params[0], fOptimiser->GetValid() );
  fFitResult.SetVertex( 0, vertex );
}

void PEnergy::SetPositiveErrors( const std::vector<double>& errors ) {
    DS::FitVertex vertex = fFitResult.GetVertex(0);
    vertex.SetPositiveEnergyError( errors[0], fOptimiser->GetValid() );
    fFitResult.SetVertex( 0, vertex );
}

void PEnergy::SetNegativeErrors( const std::vector<double>& errors ) {
    DS::FitVertex vertex = fFitResult.GetVertex(0);
    vertex.SetNegativeEnergyError( errors[0], fOptimiser->GetValid() );
    fFitResult.SetVertex( 0, vertex );
}

}       // namespace Methods
}       // namespace RAT