// Class: RAT::Methods::PEnergy // // Brief: A simulation-based energy fitter // // Author: W. Heintzelman // // REVISION HISTORY: // 19/11/2014 : W Heintzelman - New file // 29/04/2015 : W Heintzelman - Add capability to use PMT charge in fit // 09/07/2015 : GDOG - New tables for PMT noise parameters // 02/02/2016 : W Heintzelman - Take into account finite event build // interval; effect of photon polarization on PMT response // 10/06/2016 : W Heintzelman - New method for using hit times // 05/12/2016 : W Heintzelman - Mods for fitting events in water // 13/03/2017 : W Heintzelman - Change program flow to reduce execution time // 22/03/2017 : W Heintzelman - Variable number of PEs considered // 17/04/2017 : W Heintzelman - Changes needed to fit events in zdab file // 19/04/2017 : W Heintzelman - New pmtRsp table due to 3-D PMT model // updates // 08/05/2017 : W Heintzelman - Mods to allow use of ModeCut with PEnergy // 09/05/2017 : W Heintzelman - Revise seeding logic; add hit validity // check // 30/10/2017 : W Heintzelman - Improvements to noise treatment // 12/01/2018 : W Heintzelman - Recover from Geant simulation failures // // Detail: Fit for the energy of an event within the AV using // detector response derived from an auxiliary simulation of an // event located at the previously determined event position. Has // options to also take into account the recorded PMT charges and // hit times. // Additional discussion in SNO+-doc-2871, "A Simulation-Based Energy // Fitter for SNO+". // ///////////////////////////////////////////////////////////////////////////// #ifndef __RAT_Method_PEnergy_ #define __RAT_Method_PEnergy_ #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include const size_t fNPMTsDim = 9728; namespace RAT { namespace Methods { class PEnergy : public SeededMethod, public OptimisedMethod, public SelectorMethod { public: PEnergy(); void Initialise( const std::string& ); void BeginOfRun( DS::Run& run ); void SetD( const std::string& param, const double value ); void SetI( const std::string& param, const int value ); void SetS( const std::string& param, const std::string& value ); void EndOfRun( DS::Run& run ); virtual std::string GetName() const { return PEnergy::Name(); } static std::string Name() { return std::string( "pEnergy" ); } void DefaultSeed(); // Set the seed to the default void SetSeed( const DS::FitResult& seed ); // Set the seed from seed DS::FitResult GetBestFit(); // Invokes the method's algorithm // Called by the optimiser to get the starting/seed fit parameter values: virtual std::vector GetParams() const; virtual std::vector GetPositiveErrors() const; virtual std::vector GetNegativeErrors() const; // Converts the vector of params into the fFitResult void SetParams( const std::vector& params ); // Converts the vector of errors into the fFitResult virtual void SetNegativeErrors( const std::vector& errors ); virtual void SetPositiveErrors( const std::vector& errors ); // Called by SteppingAction at every particle step void AccumulateData(const G4Step* step) ; // Called by the optimiser to invoke the method algorithm, // optimiser will pass current parameters double CalcLikelihood(const vector& params, vector& llSum); virtual double operator()( const std::vector& params ); // Utility to dump the likelihood function when debugging void PrintLkhdFunction(const double eCenter) ; // Returns probability a photoelectron is produced in a PMT given photon // energy in MeV, angle of incidence in radians, and cosine of angle // between polarization direction and the senkrecht direction // (perpendicular to plane of incidence) double PMTResponse(const double photonE, const double incAngle, const double cosPol) const; private: // Set the seed vertex and time equal to the MC event position and time void SetMCSeed() ; bool DoAuxiliarySimulation(); // Initialize for accumulating data from auxiliary simulation void CalcInit(); // Determine the probabilities of photon incidence on the various PMTs from // the auxiliary simulation results void CalcIncidenceProbs(); // Calculate the expected number of photons produced given an event energy double ExpectedPhotons( double ke); // Make adjustment to fit result to compensate for energy-dependent fit // bias in water double WaterAdjustment( const double fitE); // Set upper limit on number of possible PEs considered, // based on estimated energy from the seed unsigned int GetMaxPE(double estKE); std::stringstream fMainSeqState; // To save states of main and local std::stringstream fLocalSeqState; // random number sequences Gsim* fGsimPtr; PMTCharge* fPMTChargePtr; PMTChgpdf* fPMTChgpdfPtr; PMTTimeDistns* fPMTTimeDistnsPtr; PMTTransitTime* fPMTTransitTimePtr; size_t fNPMTs; // number of PMTs unsigned int fFirstNormal; // index of first normal PMT unsigned int fLastNormal; // index of last normal PMT double rConc; // radius of PMT concentrator lip double rPMTMax; // largest distance from AV center of normal PMTs double fRTest; // value of radius to test for photon approaching PSUP double fRsqTest; // value of radius^2 to test for photon approaching PSUP size_t fMaxPE; // maximum number of possible photoelectrons to be // considered in the charge or time probability computations. G4MaterialPropertyVector* fPEEscapeProbability; // Probability of // photoelectron escaping photocathode, f(energy) G4double fCollectionEfficiency; // Collection efficiency of the PMTs -- // single value from database. Supposed to be the // probability that a photoelectron results in a signal. // See p. 54 of Phil Jones's thesis. std::vector fPedSig; // PMT charge pedestal spreads std::vector fHHP; // PMT high half-points // Polynomial coefficients for expected photons produced as a function of // event energy (used by function ExpectedPhotons): std::vector fCherCoef, fScintCoef; // Following variables declared static -- saves memory in case more than // one PEnergy fitter variation is used in the run. Also permits skipping // a new simulation if PEnergy fitter variation to be executed has same // simulation options as the one just run on the same event. static double fPIP[fNPMTsDim]; // photon incidence probability on each PMT // from auxiliary simulation static double fPMTEff[fNPMTsDim]; // average probability a photoelectron // is produced by a photon incident on a PMT -- // extracted from auxiliary simulation static double fExpNoise[fNPMTsDim]; // expected noise PEs per PMT per ns double fNoiseWindow; // Window for calculated total noise PEs static bool fOnLine[fNPMTsDim]; // true if a PMT is normal and on-line static bool fUsePMT[fNPMTsDim]; // true if a PMT is to be used in fit static bool fWasHit[fNPMTsDim]; // flag indicating if a PMT was hit in // the event being fitted static vector< vector > fPQFac; // storage of factors for // probability computation when using charge in fit static vector< vector > fPTFac; // storage of factors for // probability computation when using time in fit static unsigned int fXIdx[fNPMTsDim]; // Index from full PMT list into // list of selected (i.e., hit) PMTs static double fPhotonTransitTime[fNPMTsDim]; // transit time from fitted // position to PMT; needed if TimeResidualCut selector is used DS::FitVertex fFitVertex; G4ThreeVector fEvtPosition; // event position from seed fitter TVector3 fTVecEvtPosition; // event position from seed fitter G4ThreeVector fEvtDirection;// event direction from seed fitter TVector3 fTVecEvtDirection;// event direction from seed fitter bool fValidDirection; // true if valid direction from seed exists bool fValidSeed; // true if valid seed position double fSeedTime; // needed if hit times are used in fit double fSeedEnergy, fEnergyError; // fit starting values double fLightYield; // expected number of scintillation photons // produced per unit energy from database double fRefIdxAVFill; // index of refraction of AV fill double fPropTimeAdder; // propagation time from inner AV surface to // PSUP for photon in radial direction double fFECToTriggerDelay; double fTotalTriggerDelay; // total trigger delay double fAVInnerRadius; // inner AV radius unsigned int fNPrimaries; // Number of primary particles to be generated // in auxiliary simulation. This is usually greater // than 1 to ensure that sufficient photons are // generated for statistical reasons. std::string fAVMaterialName; // Material inside the AV bool fMakeWaterAdj; // Flag indicating if fit result is to be // adjusted to compensate for bias in water. This is // set to true within the code automatically if the // default number of photons to be simulated is not // overridden. Otherwise it is set false. std::string fSelectorName; PMTSelectors::PMTSelector* fSelector; PMTSelectors::TimeResidualCut* fTRCut; double fLowCut, fHighCut; // Ends of selector time interval // for either the ModeCut or TimeResidualCut bool fSameGTID, fSameSelector; // For getting time usage of different PEnergy components Profiler* fProfileGenPhotons; // Profiler* fProfileAccumData; Profiler* fProfileCalcIncidProb; Profiler* fProfilePFactors; Profiler* fProfileOptimize; Profiler* fProfileFinish; Profiler* fProfileTotal; // Options that can be changed by macro file inputs long fRndmNoSeed; // to specify a random number seed bool fSeedActual; // if true and if event to be fitted was simulated, // ignore seed position specified (if any) and use // MC event position in fit bool fUseCharge; // if true, fit uses charge in likelihood function bool fUseTime; // if true, fit uses hit times in likelihood function bool fUseDirection; // if true, fits use direction from seed for electrons // in auxiliary simulations -- for fits in water int fChargeType; // 0 - use QHS; 1 - use QHL; 2 - use QLX bool fElectronPrime; // false if primary particles should be photons // rather than electrons. Default is true(electrons). int fNPhotonSim; // Minimum number of photons simulated for the fit double fTabAreaFactor; // The tabulation area is this number // times the area of the PMT concentrator. double fEffScaleFactor; // overall photon detection efficiency // scaling factor int fCallCounter; // Count the number of times GetBestFit is called // Used to save GTID and settings in previous PEnergy fit to allow tests // to possibly avoid recalculations in subsequent fit of the same event static Int_t fPrevGTID; static bool fPrevSeedActual; static bool fPrevElectronPrime; static bool fPrevUseDirection; static int fPrevNPhotonSim; static double fPrevTabAreaFactor; static double fPrevEffScaleFactor; static int fPrevChargeType; static string fPrevSelector; static double fPrevLowCut, fPrevHighCut; static bool fGoodPrevAuxSim; // Arrays & parameters used in accumulating auxiliary simulation statistics: static TVector3 fPosPMT[fNPMTsDim],fDirPMT[fNPMTsDim]; static double fPosPMTUX[fNPMTsDim], fPosPMTUY[fNPMTsDim], fPosPMTUZ[fNPMTsDim]; double fMaxInclusionTime, fEstMinPropTime, fTabLim, fRhoSqLim; vector< vector > fTCross; int fNReachPSUP, fNTabRegions; }; // class PEnergy } // namespace Methods } // namespace RAT #endif