#include <G4Electron.hh>
#include <G4Event.hh>
#include <G4EventManager.hh>
#include <G4VFastSimulationModel.hh>
#include <G4GlobalFastSimulationManager.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/DiscOpticalModel.hh>
#include <RAT/GLG4Scint.hh>
#include <RAT/GLG4StringUtil.hh>
#include <RAT/Log.hh>
#include <RAT/Noise.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/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
{

// Static class variables:
double PEnergy::fPIP[fNPMTsDim];
double PEnergy::fPMTEff[fNPMTsDim];
vector< vector<float> > PEnergy::fPThresh;
vector< vector<double> > PEnergy::fPFac;
vector< vector<double> > PEnergy::fTCross;
vector< vector<double> > PEnergy::fTWeight;
unsigned int  PEnergy::fXIdx[fNPMTsDim];
double PEnergy::fExpNoise[fNPMTsDim];
bool   PEnergy::fOnLine[fNPMTsDim];
bool   PEnergy::fUsePMT[fNPMTsDim];
bool   PEnergy::fWasHit[fNPMTsDim];
unsigned int PEnergy::fNPrimaries;

TVector3 PEnergy::fDirPMT[fNPMTsDim];

Int_t  PEnergy::fPrevGTID;
bool   PEnergy::fPrevSeedActual;
bool   PEnergy::fPrevElectronPrime;
bool   PEnergy::fPrevUseTime;
bool   PEnergy::fPrevUseDirection;
int    PEnergy::fPrevNPhotonSim;
int    PEnergy::fPrevChargeType;
bool   PEnergy::fGoodPrevAuxSim = true;

float  PEnergy::fPMTRsp[2][61][46] ;
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 delWL = 10.;       // step in wavelength (in nm)
const double delAng = 2.;       // step in angle (in degrees)
const double lbWL = 200.;       // lower bound in wavelength
const double ubWL = lbWL+(nWL-1)*delWL;    // upper bound in wavelength
const double lbAng = 0.;        // lower bound in angle

// **********************************************************************
inline double PEnergy::GetPThresh(unsigned int pmtID, unsigned int nPE) {
    size_t nEntry = fPThresh[pmtID].size();
    if( nPE > nEntry ) nPE = nEntry;
    return (double)fPThresh[pmtID][nPE-1];
}
// **********************************************************************
PEnergy::PEnergy() {
}

// **********************************************************************
void PEnergy::Initialise( const std::string& ) {
    fRndmNoSeed = 0;
    // Default values:
    fSeedActual = false;
    fUseCharge = false;
    fUseTime = false;
    fChargeType = 0;    // Current default charge type is QHS
    fElectronPrime = true;
    fNPhotonSim = 150000;  // default for scintillator.  Water default is
                                // 200,000, set in BeginOfRun
    fRespTableName = "Jan2019";
    fEffScaleFactor = 0.562;
    fUseDirection = false;
    fApplyEAdj = false;
    fApplyRAdj = false;
    fPartialFill = "";   // Set to "top" or "bottom" to fit in top or
                         // bottom of partial AV fill
    fSetNPrimaries = 0;
    fSetEPrimary = 0.;
}

// **********************************************************************
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);

    fPrevGTID=-10;
    fPrevChargeType=-10;

    DB* db = DB::Get();
    info << "PEnergy PMT response table: "<<fRespTableName << endl;
    RAT::DBLinkPtr respPtr =  db->GetLink("PMT_RESPONSE_TABLE",fRespTableName);
    vector<double> respTable = respPtr->GetDArray("respTable");
    if( respTable.size() != nPol*nWL*nAng )
                RAT::Log::Die( "\nPEnergy: PMT response table size mismatch" );
    size_t ict = 0;
    for(int i=0; i<nPol; i++) {
        for(int j=0; j<nWL; j++) {
            for(int k=0; k<nAng; k++) {
                fPMTRsp[i][j][k] = respTable[ict];
                ict++;
            }
        }
    }

    // Check if grey disc PMT model is being used for inner PMTs
    RAT::DBLinkPtr ptrpmt =  db->GetLink("GEO","innerPMT");
    std::vector<int> greyDisc=  ptrpmt->GetIArray("grey_disc");
    fUseGreyDisc = (greyDisc[0] == 1) ;
    detail << "PEnergy:  fUseGreyDisc = " << fUseGreyDisc << endl;

    string getName = "material";
    if( fPartialFill == "top" ) {
        getName = "material_top";
    }else if( fPartialFill == "bottom" ) {
        getName = "material_bottom";
    }
    RAT::DBLinkPtr ptrav =  db->GetLink("GEO","inner_av");
    G4String avMaterialName;
    try {
        avMaterialName = ptrav->GetS( getName );
    } 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=="internalwater_snoplus" ||
                                    avMaterialName=="lightwater_sno") {
        fEffExtraFactor = respPtr->GetD("adjfactor_water");
        fLightYield = 0.;
        fNPhotonSim = 200000;  // default value for water
    }else{

        fEffExtraFactor = respPtr->GetD("adjfactor_labppo");
        if(avMaterialName=="te_diol_dda_0p5_labppo_scintillator_bisMSB_Jan2018")
                        fEffExtraFactor = respPtr->GetD("adjfactor_tejan2018");

        // Following settings were made for LAB/PPO;
        // Further adjustments should be made in FIT_PEnergy.ratdb
        if( !fUseCharge && !fUseTime) {
            fEffExtraFactor *= 1.0;
        }
        if( fUseCharge && !fUseTime) {
            if( fChargeType == 0) fEffExtraFactor *= 0.999805637;  // QHS
            if( fChargeType == 1) fEffExtraFactor *= 1.009943666;  // QHL
        }
        if( fUseTime) {
            if( fChargeType == 0) fEffExtraFactor *= 1.022796533;  // QHST
            if( fChargeType == 1) fEffExtraFactor *= 1.025637135;  // QHLT
        }

        // Material assumed to be a scintillator, with a light yield.
        fLightYield = avMaterialPtr->GetMaterialPropertiesTable()
                                        ->GetConstProperty("LIGHT_YIELD");
    }

    if( fApplyEAdj) {
        RAT::DBLinkPtr eAdjPtr =  db->GetLink("E_ADJUST_FN",fEAdjFnName);
        fAdjFnFitE = eAdjPtr->GetDArray("fite");
        fAdjFnAdjE = eAdjPtr->GetDArray("adje");
        info << "PEnergy energy adjustment function: "<<fEAdjFnName << endl;
        if (fAdjFnFitE.size() != fAdjFnAdjE.size() ) RAT::Log::Die
                    ("PEnergy: energy adjustment arrays of unequal length");
    }

    if( fApplyRAdj) {
        RAT::DBLinkPtr rAdjPtr =  db->GetLink("R_ADJUST_FN",fRAdjFnName);
        fAdjFnFitR = rAdjPtr->GetDArray("fitr");
        fAdjFnRFactor = rAdjPtr->GetDArray("radjfactor");
        info << "PEnergy radial adjustment function: "<<fRAdjFnName << endl;
        if (fAdjFnFitR.size() != fAdjFnRFactor.size() ) RAT::Log::Die
                    ("PEnergy: radial adjustment arrays of unequal length");
    }

    // Get coefficients from database for function to return the expected
    // number of photons given event energy:
    fAVMaterialName = avMaterialName;
    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" );

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

    DBLinkPtr ldaq = db->GetLink("DAQ");
    if(fUseTime) {
        fUseCharge = true;
        fPulseWidth = ldaq->GetD("pulse_width");
    }

    // Following copied from UnCal.cc
    int ECApattern = ldaq->GetI("eca_pdst_pattern");
    int ECArate = ldaq->GetI("eca_pdst_rate");
    //Get the pattern and rate
    std::stringstream eca_pattern_rate;
    eca_pattern_rate << ECApattern <<"_"<<ECArate;
    DBLinkPtr pedsbank = db->GetLink("ECA_PDST", eca_pattern_rate.str());
    DBLinkPtr pmtCalib = db->GetLink("PMTCALIB");
    if(fChargeType==0) {
        fPedSig    = pedsbank->GetFArrayFromD("pdst_qhswidth");
//      Currently getting thresholds from ChanHWStatus
//      fThreshold = pmtCalib->GetFArrayFromD("QHS_threshold");
        fHHP       = pmtCalib->GetFArrayFromD("QHS_hhp");
        fIntTime   = ldaq->GetD("shortint_time");
    }else if(fChargeType==1) {
        fPedSig    = pedsbank->GetFArrayFromD("pdst_qhlwidth");
//      fThreshold = pmtCalib->GetFArrayFromD("QHL_threshold");
        fHHP       = pmtCalib->GetFArrayFromD("QHL_hhp");
        fIntTime   = ldaq->GetD("longint_time");
    }else{
        fPedSig    = pedsbank->GetFArrayFromD("pdst_qlxwidth");
//      fThreshold = pmtCalib->GetFArrayFromD("QLX_threshold");
        fHHP       = pmtCalib->GetFArrayFromD("QLX_hhp");
        fIntTime   = ldaq->GetD("longint_time");
    }
    fLongIntTime   = ldaq->GetD("longint_time");
    // Decrease charge integration time by 10 ns because total integration time
    // includes 10 ns before the hit that caused the PMT to fire.
    fIntTime -= 10.;

    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 (the null) selector.");
    fSelectorName = (*fSelectors.begin())->GetName();
    fSelector = *fSelectors.begin();
    fSelector->BeginOfRun( run );
    if( fSelectorName != "null" ) {
        RAT::Log::Die( "PEnergy expects the null selector");
    }
    detail<<"PEnergy:  Using the " << fSelectorName << " selector." << 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 ) {
#ifdef PENERGYDETAIL
        detail << "Using PMT noise rate of zero." << endl;
#endif
        // 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();

    fFirstNormal = 0;
    fLastNormal = 0;
    double sumr = 0.;
    int nNormal = 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.
        if(isNormal) {
            fLastNormal = k;
            if( fFirstNormal == 0) fFirstNormal = k;

            double rPMT = (pmtInfo.GetPosition(k)).Mag();
            sumr += rPMT;
            nNormal++;

            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);
            }
        }
    }
    fRAvg = sumr/(double)nNormal;

    // Fill fPThresh array and save PMT directions
    unsigned int dimPE = fPMTChgpdfPtr->GetDimPE();
    int nOnLine = 0;
    fPThresh.clear();
    for(unsigned int k=0; k<=fLastNormal; k++) {
        fDirPMT[k] = pmtInfo.GetDirection(k);
        vector<float> row;
        if( fOnLine[k] ) {
            double hhp = fHHP[k];
            double threshold = chanHWStatus.GetThresholdAboveNoise(k) ;
            if (threshold<0.5*hhp) {
                for(unsigned int j=0; j<dimPE; j++) {
                    double p = 1. -
                        fPMTChgpdfPtr->ProbBelowThreshold(j+1,hhp,threshold);
                    row.push_back(p);
                    if( p > 0.9999999 ) break;
                }
                nOnLine++;
            }else{
                fOnLine[k]=false;
            }
        }
        fPThresh.push_back(row);
    }

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

    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
    if( fPartialFill == "top" ) {
        fRefIdxAVFill = lightPathCalc.GetUpperTargetRI(photonMeV);
    }else if( fPartialFill == "bottom" ) {
        fRefIdxAVFill = lightPathCalc.GetLowerTargetRI(photonMeV);
    }else{
        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 +
                      (fRAvg- 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( "pmtEffFactor" ) ) {
        fEffScaleFactor = value;
    }else if ( param == string( "energyPrimary" ) ) {
        fSetEPrimary = 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;
    }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 if( param == string( "nPrimaries" ) ) {
        fSetNPrimaries = value;
    }else {
        throw Processor::ParamUnknown( param );
    }
}

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

void PEnergy::SetS( const string& param, const string& value ) {
    if( param == string( "optionPkg" ) ) {
        // This parameter setting causes several PEnergy-specific fitting
        // options to be read as a group from FIT_PEnergy.ratdb.
        DB* db = DB::Get();
        RAT::DBLinkPtr optionPkg =  db->GetLink("OPTION_PKG",value);
        fRespTableName = optionPkg->GetS("pmtRespTable");
        fApplyEAdj = true; fEAdjFnName = optionPkg->GetS("energyAdjFn");
        fApplyRAdj = true; fRAdjFnName = optionPkg->GetS("radialAdjFn");
        fUseDirection = optionPkg->GetI("useDirection");
    }else 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 if( param == string( "pmtRespTable" ) ) {
        if( value != "" ) fRespTableName = value;
    }else if( param == string( "energyAdjFn" ) ) {
        if( value != "" ) {
            fApplyEAdj = true;
            fEAdjFnName = value;
        }
    }else if( param == string( "radialAdjFn" ) ) {
        if( value != "" ) {
            fApplyRAdj = true;
            fRAdjFnName = value;
        }
    }else if( param == string( "partialFill" ) ) {
        if(value == "top" || value == "bottom") {
            fPartialFill = value;
        }else{
            RAT::Log::Die( "\nPEnergy: partialFill entry error" );
        }
    }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=="internalwater_snoplus" ||
                        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=="internalwater_snoplus" ||
                                 fAVMaterialName=="lightwater_sno") &&
                                    fSeedEnergy<0.8 ) fSeedEnergy=0.8;
            if( fFitVertex.ValidPositiveEnergyError() ) {
                fEnergyError = fFitVertex.GetPositiveEnergyError();
            }
        }

        if( fFitVertex.ContainsTime() ) {
            fSeedTime = fFitVertex.GetTime();
            if( fUseTime && fSeedTime<0. ) {
                fValidSeed=false;
                warn << "PEnergy::SetSeed: seed time < 0; "
                    << "fit not attempted." <<endl;
            }
        }else{
            if(fUseTime) {
                fValidSeed=false;
                if(!fSeedActual) warn <<
                    "PEnergy::SetSeed: seed time invalid; "
                    << "fit not attempted." <<endl;
            }
        }

        if( fFitVertex.ContainsPosition() && fFitVertex.ValidPosition() ) {
            fTVecEvtPosition = fFitVertex.GetPosition();
            if( fTVecEvtPosition.Mag()>8000. && !fSeedActual) {
                fValidSeed = false;
                info<< "PEnergy::GetBestFit -- Seed position at r>800 cm; "
                    << "fit not attempted" <<endl;
            }
        }else{
            fValidSeed = false;
            if(!fSeedActual) warn <<
                "PEnergy::SetSeed: seed position invalid; "
                << "fit not attempted." <<endl;
        }

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

    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());

    if(fSetEPrimary != 0.) fSeedEnergy = fSetEPrimary;

#ifdef PENERGYSEED
    if(fValidSeed) {
        detail << "PEnergy::SetSeed - Evt kE, r, pos(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 and the fitted event time are 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 (which is always zero).  Here we set the
    // seed time equal to the MC event time adjusted to be 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
    if( !fValidSeed) {
        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);
    bool skipSim = (    fSameGTID &&
                        fPrevSeedActual==fSeedActual &&
                        fPrevElectronPrime==fElectronPrime &&
                        fPrevUseDirection==fUseDirection &&
                        fPrevNPhotonSim==fNPhotonSim &&
                        fSetNPrimaries == 0 &&
                        fSetEPrimary == 0.
                    );
    fPrevGTID=gtid ;
    fPrevSeedActual=fSeedActual ;
    fPrevElectronPrime=fElectronPrime ;
    fPrevUseDirection=fUseDirection ;
    fPrevNPhotonSim=fNPhotonSim ;

    // Get list of selected PMTs (should be all PMTs because the NULL selector
    // is used).  fSelectedPMTData is inherited from SelectorMethod class and
    // 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];
    }

    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, is not cross-talk, and
        // the hit time is no more than 30ns before the seed time:
        // Need to apply the hit time test to avoid problems with time
        // distributions.
        fUsePMT[idx] &= ( pmtCalStatus.GetHitStatus(fSelectedPMTData[iSel])==0
                && (fSelectedPMTData[iSel].GetTime()-fSeedTime)>-30. );
    }

#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

    // 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.
    // This is used as an estimate of the time of the PMT hit(s) that produce
    // the trigger.
    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 (relative to the event time)
    //            for a PE that would be included in charge integration for an
    //            event.  Used for initial test if an incident photon in
    //            auxiliary simulation should be counted.
    fMaxInclusionTime = fEstMinPropTime + fTotalTriggerDelay + fLongIntTime;

    fPMTTimeDistnsPtr->SetEventParameters(
                          fEstMinPropTime + fFECToTriggerDelay +10., fMaxPE );

    // 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;
    }
    // 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(skipSim);

    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();
    // First apply radial and energy adjustments if specified
    if(fApplyRAdj) params[0]=RadialAdjustment(fTVecEvtPosition.Mag(),params[0]);
    if(fApplyEAdj) params[0] = EnergyAdjustment( 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() );
    // Save the numbers of Cherenkov and scintillation photons generated in the
    // auxiliary simulation
    fFitResult.SetFOM( "CherPhotCount", (double)fCherPhotCount );
    fFitResult.SetFOM( "ScintPhotCount", (double)fScintPhotCount );
    // Save the total number of photon incidences in the auxiliary simulation
    fFitResult.SetFOM( "IncdPhotCount", (double)fIncdPhotCount );

#ifdef PENERGYFITOUT
    std::stringstream strFitName;
    strFitName << "PEnergy    " ;
    int chgSetting=0;
    if(fUseCharge) chgSetting = fChargeType +1;
    info << strFitName.str().substr(0,15) <<
        fSeedActual <<
        chgSetting<<
        fUseTime <<
        fUseDirection <<" "<<
        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::EnergyAdjustment( const double fitE) {
    // Make piece-wise linear adjustment to compensate for
    // energy-dependent fit bias.
    double f;
    unsigned int n = fAdjFnFitE.size();
    if( fitE<=fAdjFnFitE[0] ) {
        f = fAdjFnAdjE[0]/fAdjFnFitE[0];
    }else if( fitE>=fAdjFnFitE[n-1] ) {
        f = fAdjFnAdjE[n-1]/fAdjFnFitE[n-1];
    }else{
        for(unsigned int i=0; i<=n-2; i++) {
            if( fitE>fAdjFnFitE[i] && fitE<=fAdjFnFitE[i+1] ) {
                double f1 = fAdjFnAdjE[i]/fAdjFnFitE[i];
                double f2 = fAdjFnAdjE[i+1]/fAdjFnFitE[i+1];
                f = f1 + (f2-f1)*(fitE - fAdjFnFitE[i])/
                            (fAdjFnFitE[i+1] - fAdjFnFitE[i]);
                break;
            }
        }
    }
    return f*fitE;
}

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

double PEnergy::RadialAdjustment( const double fitR, const double fitE) {
    // Make piece-wise linear adjustment to compensate for
    // radially dependent fit bias.
    double f;
    unsigned int n = fAdjFnFitR.size();
    if( fitR>=fAdjFnFitR[n-1] ) {
        f = fAdjFnRFactor[n-1];
    }else{
        for(unsigned int i=0; i<=n-2; i++) {
            if( fitR>=fAdjFnFitR[i] && fitR<=fAdjFnFitR[i+1] ) {
                double f1 = fAdjFnRFactor[i];
                double f2 = fAdjFnRFactor[i+1];
                f = f1 + (f2-f1)*(fitR - fAdjFnFitR[i])/
                            (fAdjFnFitR[i+1] - fAdjFnFitR[i]);
                break;
            }
        }
    }
    return f*fitE;
}

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

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));
            if(fSetNPrimaries != 0) fNPrimaries = fSetNPrimaries;
#ifdef PENERGYDETAIL
            detail << "PEnergy::DoAuxiliarySimulation -- "<< endl;
            detail << "    fNPhotonSim, fSeedEnergy, ExpectedPhotons = " <<
                    fNPhotonSim<<" "<< fSeedEnergy <<" "<<
                    ExpectedPhotons(fSeedEnergy) << endl;
            detail << "    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:
        fPrevTrackID = -1;      // Initialize for use in AccumData
        fProfileGenPhotons->Start();
        evtManager->ProcessOneEvent(auxEvent);
        fProfileGenPhotons->Stop();
        haveGoodAuxSim = asiPtr->IsSimEvtGood();

        fCherPhotCount  = 0.;
        fScintPhotCount = 0.;
        if( haveGoodAuxSim ) {
            fCherPhotCount  = Cerenkov::GetPhotonCount();
            fScintPhotCount = GLG4Scint::GetScintillatedCount();
#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

    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.;
    }

    // Arrays for saving times at which photons intersect each PMT's tabulation
    // area, and the corresponding weights:
    fTCross.clear();
    fTCross.resize(fNPMTsDim);
    fTWeight.clear();
    fTWeight.resize(fNPMTsDim);
    fNPhotonIncd = 0 ; // to count number of tabulation areas incremented

#ifdef PENERGYDBG
    detail << "Leaving PEnergy::CalcInit" << endl;
#endif
}

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

void PEnergy::AccumulateData(const G4Step* step)  {
//     fProfileAccumData->Start();

    G4Track* g4TrackPtr = step->GetTrack();
    if( g4TrackPtr->GetParticleDefinition()->GetParticleName()
                                                != "opticalphoton" ) return;
    G4int trackID = g4TrackPtr->GetTrackID();
    if(trackID != fPrevTrackID) {               // new track
        fPrevVolName = "None";
        fPrevTrackID = trackID;
    }

    G4StepPoint* preStepPt = step->GetPreStepPoint();
    string preStepVolName = preStepPt->GetPhysicalVolume()->GetName();
    G4StepPoint* postStepPt = step->GetPostStepPoint();

    // Save step info if step began in the cavity:
    if (preStepVolName=="cavity") {
        fPrePosition  = preStepPt->GetPosition();
        fPostPosition = postStepPt->GetPosition();
        fPostPolarization = postStepPt->GetPolarization();
        fPostPtTime = postStepPt->GetGlobalTime();
        fEKin = postStepPt->GetKineticEnergy();
    }

    // Return if previous step did not begin in cavity:
    if( fPrevVolName != "cavity") {
        fPrevVolName = preStepVolName;
        return;
    }
    fPrevVolName = preStepVolName;

    // Return if this step does not begin in an innerPMT_pmtenv:
    if( preStepVolName.find("innerPMT_pmtenv") == string::npos) return;

    string postStepVolName = postStepPt->GetPhysicalVolume()->GetName();
    bool keepflag = true;
    if(fUseGreyDisc) {
        // When grey-disc model is in use, photon steps tested for disc
        // intersection should have the same pre- and post-step volume
        if (preStepVolName != postStepVolName) {
            string msg = "\nPEnergy: Pre- and post-step volumes different "
                + preStepVolName + "  " + postStepVolName ;
            RAT::Log::Die( msg );
        }
        G4VFastSimulationModel* fsm =
            G4GlobalFastSimulationManager::GetGlobalFastSimulationManager()->
            GetFastSimulationModel("innerPMT_optical_model");
        DiscOpticalModel* discOptModel = dynamic_cast<DiscOpticalModel*> (fsm);
        keepflag = discOptModel->GetLastTestResult();
    }else{
        // Require that this step ends at the concentrator, the PMT, or
        // is absorbed or scattered within the pmtenv.
        if ( postStepVolName != "innerPMT_concentrator" &&
             postStepVolName != "innerPMT_pmt" ) {
            keepflag = false;
            // Retain cases in which photon is absorbed or scattered in water
            // within pmtenv after entering:
            if( postStepVolName.find("innerPMT_pmtenv") != string::npos ) {
                string processName = postStepPt->GetProcessDefinedStep()
                                                            ->GetProcessName();
                keepflag = ( processName == "OpAbsorption" ||
                             processName == "OpRayleigh"   ||
                             processName == "OpMie"        );
            }
        }
    }
    if( !keepflag ) return;

    // Extract the pmt ID from the preStepVolName:
    string strpmtID = preStepVolName.substr(15,4);
    int  pmtID = util_to_int(strpmtID);
    if( !fUsePMT[pmtID]) return;

    // We want to use the track segment from a point in the cavity to a point
    // at the pmtenv edge to calculate the incidence angle, polarization, and
    // time.

    TVector3 prevPos, ptPos, polarization;
    for(int i=0; i<=2; i++) {
        ptPos(i) = fPostPosition(i);
        prevPos(i) = fPrePosition(i);
        polarization(i) = fPostPolarization(i);
    }

    TVector3 segmentVect = ptPos - prevPos;
    TVector3 photonDirection = segmentVect.Unit();

    double cosIncdAngle=photonDirection.Dot(fDirPMT[pmtID]);
    if(cosIncdAngle>1.) cosIncdAngle = 1.;
    if(cosIncdAngle>0.){

        // Require that the time is within the interval for a hit
        // generated by the photon to be included in the event:
        if( fPostPtTime<fMaxInclusionTime ) {
#ifdef PENERGYDETAIL
            fNPhotonIncd++;
#endif
            fPIP[pmtID] += 1;
            double incdAngle = acos(cosIncdAngle);
            // vector perpendicular to the plane of incidence:
            TVector3 senkrecht =
                    (fDirPMT[pmtID].Cross(photonDirection)).Unit();
            double cosPol = senkrecht.Dot(polarization);
            double pmtResp = PMTResponse(fEKin,incdAngle,cosPol);

            fPMTEff[pmtID] += pmtResp;
            if(fUseTime) {
                fTCross[pmtID].push_back(fPostPtTime);
                fTWeight[pmtID].push_back(pmtResp);
            }
        }    //if( t<fMaxInclusionTime )
    }      // if(cosIncdAngle>0. )


//     fProfileAccumData->Stop();
    return;
}       //  PEnergy::AccumulateData

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

void PEnergy::CalcIncidenceProbs(bool skipSim) {
#ifdef PENERGYDBG
    detail << "Entering PEnergy::CalcIncidenceProbs" << endl;
#endif
double noiseMult = fNPrimaries;
if (!skipSim) {
    // 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./ (double)(fCherPhotCount+fScintPhotCount);
        // Denominator is total number of photons created in the auxiliary
        // simulation
    }else{
        adjFactor = 1./( (double)fNPhotonSim );
    }

    PMTVariation* pmtVariationPtr = PMTVariation::Get();
    fIncdPhotCount =0;
    for(unsigned int k=fFirstNormal; k<=fLastNormal; k++) {
        if( fUsePMT[k] ) {
            // Note: At this point, fPIP[k] is the number of photons incident
            // on PMT k in the auxiliary simulation, and fPMTEff[k] is the sum
            // of the PMTResponse values for the incident photons.

            if( fPIP[k]>0. ) {
                fIncdPhotCount += fPIP[k];
                // Calculate average PMT response and multiply by various other
                // factors
                fPMTEff[k] = fCollectionEfficiency*fPMTEff[k]
                               *pmtVariationPtr->GetVariation(k)/fPIP[k];
            }else{  // No simulated photon incident on PMT; use response for
                    // wavelength of 370 nm and 30 deg incidence angle.
                fPMTEff[k] = fCollectionEfficiency
                            *PMTResponse(3.35*electronvolt, 30./radToDeg, 0.71)
                            *pmtVariationPtr->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 were 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();
}

    // Test if we can avoid redoing the fPFac computations
    bool skipPFacComp = ( skipSim && fPrevChargeType==fChargeType
            && fPrevUseTime==fUseTime );
    fPrevChargeType = fChargeType;
    fPrevUseTime = fUseTime;

    // Not using charge:
    //   For hit PMT i, fPFac[i][j] is the probability that the total
    //   charge is above the threshold given j+1 photoelectrons, i.e., the
    //   probability that j+1 photoelectrons will produce a hit.
    if ( !fUseCharge && !skipPFacComp ) {
        fProfilePFactors->Start();
        fPFac.clear();
        fPFac.resize( fNPMTsDim, vector<double>(fMaxPE,0.) );
        for(unsigned int i=fFirstNormal; i<=fLastNormal; i++) {
            if( fUsePMT[i] ) {
                for(unsigned int j=0; j<fMaxPE; j++) {
                    fPFac[i][j] = GetPThresh(i,j+1);
                }
            }
        }
        fProfilePFactors->Stop();
    }

    if ( fUseCharge && !skipPFacComp ) {
        fProfilePFactors->Start();
        const DU::ChanHWStatus&
                        chanHWStatus= DU::Utility::Get()->GetChanHWStatus();
        // 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.
        fPFac.clear();
        fPFac.resize( fNPMTsDim, vector<double>(fMaxPE,0.) );

        // Calculate fPFac[][].  For hit PMT i, using charge but not time,
        // fPFac[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, fPFac[i][j] is the probability of no hit given j+1
        // photoelectrons.
        // Using both charge and time, for hit PMT i, fPFac[i][j] is the
        // probability of a hit with the observed charge and time given j+1
        // photoelectrons.

        // Calculate probability factors for PMTs that were hit:
        size_t nHitPMTs= fSelectedPMTData.size();
        for(size_t i=0; i<nHitPMTs; i++) {
            UInt_t hitPMTID = fSelectedPMTData[i].GetID();
            if( hitPMTID>=fNPMTsDim ) {
                std::stringstream msg;
                msg << "PEnergy::CalcIncidenceProbs -- dimension error: "
                        << hitPMTID <<" "<< fNPMTsDim ;
                RAT::Log::Die(msg.str());
            }
            double hhp = fHHP[ hitPMTID ];
            double threshold= chanHWStatus.GetThresholdAboveNoise(hitPMTID);

            if( fUsePMT[hitPMTID] ) {
                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();
                }

                vector<double> pQ;
                pQ.resize(fMaxPE+1);
                for(size_t  j=1; j<=fMaxPE; j++) {
                    // Probability of total recorded charge q from j PEs
                    // (i.e., given that charge is above threshold):
                    pQ[j] = fPMTChgpdfPtr->
                                PdfRecordedVal(j,q,hhp,threshold,pedSig);
                }

                vector< vector<double> > g, Gwid, Gprm, Glim;
                vector< vector<double> > pIntDelta, pPulseDelta;

                if(fUseTime ) {
                    double relHitTime =fSelectedPMTData[i].GetTime()- fSeedTime;
                    double nExpNoisePEs =
                                    noiseMult*fNoiseWindow*fExpNoise[hitPMTID];
                    fPMTTimeDistnsPtr->CalcOrderStats(
                        nExpNoisePEs, fTCross[hitPMTID], fTWeight[hitPMTID],
                        relHitTime, hitPMTID, fPulseWidth, fIntTime,
                        fMaxInclusionTime , g, Gwid, Gprm, Glim);

                    fPMTTimeDistnsPtr->CalcIntervalProbs(nExpNoisePEs,
                        fTCross[hitPMTID], fTWeight[hitPMTID],
                        relHitTime, hitPMTID, fIntTime, pIntDelta);

                    fPMTTimeDistnsPtr->CalcIntervalProbs(nExpNoisePEs,
                        fTCross[hitPMTID], fTWeight[hitPMTID],
                        relHitTime, hitPMTID, fPulseWidth, pPulseDelta);

                }

                // Loop over possible total numbers of PEs
                for(size_t j=0; j<fMaxPE; j++) {
                    size_t jPE = j+1;              // number of PEs
                    if(hitPMTID>=fNPMTsDim || jPE>fMaxPE)
                        warn << "fPFac index outside of range: " <<
                        hitPMTID<<" "<<fNPMTsDim <<" "<< jPE<<" "<<fMaxPE<<endl;

                    if(fUseTime ) {             // using both charge and time
                        double sumk = 0.;
                        double OneMHProd = 1.;
                        // loop over total number being considered
                        for(size_t k=1; k<=jPE; k++) {
                            // Calculate probability of observed charge given
                            // total number of PEs = k
                            double sumii = 0.;
                            for(size_t ii=1; ii<=jPE-k+1; ii++) {
                                sumii += pIntDelta[jPE-k][ii-1]*pQ[ii];
                            }
                            double summ = 0.;
                            for(size_t m=1; m<=jPE-k+1; m++) {
                                summ += pPulseDelta[jPE-k][m-1]*
                                                GetPThresh(hitPMTID,m);
                            }
                            // Limited precision can cause summ to exceed 1 by
                            // a small amount
                            if(summ>1.) summ = 1.;
                            double H = OneMHProd*summ;
                            sumk += H*g[jPE][k]*sumii;
                            OneMHProd *= (1. - H);
                        }       // for(size_t k=1; k<=jPE; k++)
                        fPFac[hitPMTID][j] = sumk;
                    }else{                      // using charge but not time
                        // hit probability factor (probability total charge is
                        // above the threshold):
                        // Probability of a hit with the observed charge if
                        // there are jPE=j+1 PEs:
                        fPFac[hitPMTID][j] = GetPThresh(hitPMTID,jPE)*pQ[jPE];
                    }
                    // Following needed to avoid prob=0 in operator()
                    // function, which produces an error in log(prob).
                    if(fPFac[hitPMTID][j]== 0.) fPFac[hitPMTID][j]=1.e-33;
                }
            }   //  if( fUsePMT[hitPMTID] )
        }       // for(size_t i=0; i<nHitPMTs; i++)

        // Calculate probability factors for PMTs that were not hit
        for(unsigned int i=fFirstNormal; i<=fLastNormal; i++) {
            if(fUsePMT[i] && !fWasHit[i] ) {
                if(fUseTime) {
                    // When using time, we assume that each PE is an
                    // independent trial in firing the PMT.
                    double qF= 1. - GetPThresh(i,1);
                    double qFexp = qF;
                    for(unsigned int j=0; j<fMaxPE; j++) {
                        if( qFexp > 0. ) {
                            fPFac[i][j] = qFexp;
                        }else{
                            fPFac[i][j] = 1.e-33;
                        }
                        qFexp *= qF;
                    }
                }else{      // if(fUseTime)
                    // When using charge without time, we assume that the
                    // PE charges add in determining whether the PMT fires.
                    for(unsigned int j=0; j<fMaxPE; j++) {
                        fPFac[i][j] = 1. - GetPThresh(i,j+1);
                    }
                }      // if(fUseTime)
            }           // if(fUsePMT[i] && !fWasHit[i] )
        }       // for(int i=fFirstNormal; i<=fLastNormal; i++)

        fProfilePFactors->Stop();
    }    // if ( fUseCharge && !skipPFacComp )

#ifdef PENERGYDETAIL
  detail<<"Number of photon incidences on PMTs = "<<fNPhotonIncd <<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>ubWL ) {
        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;
        float omt = 1.-t;
        float omu = 1.-u;
        double pmtRespP = omu*( omt*fPMTRsp[0][ia][ja] + t*fPMTRsp[0][ia+1][ja])
                    +u*( t*fPMTRsp[0][ia+1][ja+1] + omt*fPMTRsp[0][ia][ja+1] );
        double pmtRespS = omu*( omt*fPMTRsp[1][ia][ja] + t*fPMTRsp[1][ia+1][ja])
                    +u*( t*fPMTRsp[1][ia+1][ja+1] + omt*fPMTRsp[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 are obtained from fits to RAT simulation
    // outputs of electron events in the medium in question.

    // 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]);

    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;

    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]*
                            fEffExtraFactor*fEffScaleFactor*fPMTEff[i] ;
            // Expected number of noise photoelectrons:
            double expNoisePEs = fNoiseWindow*fExpNoise[i];
            bool mostlyNoise = ( lambda<expNoisePEs );
            // Add the noise contribution:
            lambda += expNoisePEs;

            if(fUseCharge) {    // Note: if using time, than by default also
                                // using charge

                double probPE = exp(-lambda); // probability of zero PEs
                double sumProbPE = probPE ;
                double sumProb ;
                int jLim = fMaxPE; if(mostlyNoise) jLim = 1;
                if( fWasHit[i] ) {
                    sumProb = 0.;
                    for(int j=0; j<jLim; j++) {
                        probPE = probPE*lambda/(double)(j+1);
                        sumProbPE += probPE ;
                        sumProb += probPE*fPFac[i][j];
                    }
                }else{
                    sumProb = probPE;
                    for(int j=0; j<jLim; j++) {
                        probPE = probPE*lambda/(double)(j+1);
                        sumProbPE += probPE ;
                        sumProb += probPE*fPFac[i][j];
                    }
                }
                if(sumProbPE > 0. && sumProbPE < 1.)
                            sumProb += (1.-sumProbPE)*fPFac[i][jLim-1];

                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;
                    }
                }

            }else{   // if(!fUseCharge)
                double pPoisson  =  exp(-lambda);  // probability of no PEs
                double sumPPois  = pPoisson;
                double pHit = 0.;

                for(unsigned int j=1; j<=fMaxPE; j++) {
                    pPoisson = pPoisson*lambda/(double)j;
                    sumPPois  += pPoisson;
                    // pPoisson is now the probability that j PEs are produced
                    // in the PMT.
                    pHit = pHit + pPoisson*fPFac[i][j-1];
                }
                // Add in residual probability for numbers of PEs above
                // fMaxPE, with assumption that with that many PEs, a hit is
                // certain:
                if(sumPPois < 1.) pHit += (1.-sumPPois);

                if( fWasHit[i] ) {
                    if(pHit==0.) pHit = 1.e-33;
                    logLike +=  log(pHit) ;
                }else{
                    if(pHit>=1.) pHit = 1. - 1.e-33;
                    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