// 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 // 18/05/2018 : W Heintzelman - Extensive revisions when using both // charge & time // 28/01/2019 : W Heintzelman - Use Geant geo info to determine photon // incidences // 07/02/2019 : W Heintzelman - Store probabilities of exceeding charge // thresholds in array to avoid recalculation // 11/10/2019 : W Heintzelman - Mods to allow use with grey-disc model // 13/05/2020 : W Heintzelman - Enhancements for fitting with partial fill // 04/05/2021 : M Lei & B Land - Add capability to use Chroma for internal // optical simulation of PMT response. // // 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+". // // Note: Numerous variables have a dimension of 2 (fLightYield, // fRefIdxAVFill, fAVMaterialName, fApplyEAdj, fApplyRAdj, // fEAdjFnName, fRAdjFnName, fCherFnSelPt, fNPhotonSim, // fEffExtraFactor, fUseDirection), and others have their row // dimension set to 2 within the executable part of the code (e.g., // fAdjFnFitE, fAdjFnAdjE, fAdjFnFitR, fAdjFnRFactor, fCherCoefLo, // fCherCoef, fScintCoef, fHitFactor). This is necessary to provide // for fitting events in both the upper and lower volumes of the AV // during partial fill. When the AV is filled entirely with a single // material, only the zeroth element of these variables is used. In // partial fill, the zeroth element applies to the upper volume, and // the oneth element applies to the lower volume. ///////////////////////////////////////////////////////////////////////////// #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 { #ifdef USE_CHROMA class HitData; #endif 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) ; // Chroma version of accumulate data #ifdef USE_CHROMA void AccumulateData(HitData* hit_photons); #endif // Called by the optimiser to invoke the method algorithm, // optimiser will pass current parameters double CalcLikelihood( const std::vector& params, std::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(bool skipSim); // Calculate the expected number of photons produced given an event energy double ExpectedPhotons( double ke); // Make adjustments to fit result to compensate for residual fit biases double EnergyAdjustment( const double fitE); double RadialAdjustment( TVector3& fitPos, const double fitE); // Set upper limit on number of possible PEs considered, // based on estimated energy from the seed unsigned int GetMaxPE(double estKE); // Extract value from fPThresh vector array (This is an in-line function.) double GetPThresh(unsigned int pmtID, unsigned int nPE); // Following to hold the table of 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). // The dimensions 2, 61, and 46 must be equal to the values of the // const int variables nPol, nWL, and nAng defined in PEnergy.cc static float fPMTRsp[2][61][46]; 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 fRAvg; // Average PMT radius size_t fMaxPE; // maximum number of possible photoelectrons to be // considered in the charge or time probability computations. 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. G4MaterialPropertyVector* fPEEscapeProbability; // Probability of Photoelectron escaping photocathode, f(energy) 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< std::vector > fCherCoefLo, fCherCoef, fScintCoef; double fCherFnSelPt[2]; // Approximate number of hits per MeV of event energy: double fHitFactor[2]; // 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 long fCherPhotCount; // Cherenkov photons in auxiliarly simulation long fScintPhotCount; // Scintillation photons in auxiliarly simulation long fIncdPhotCount; // Photons incident on PMTs in aux sim 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 std::vector< std::vector > fPThresh; // fPThresh[i][j] is the // probability that the total charge is above threshold for PMT i // given j+1 PEs. static std::vector< std::vector > fPFac; // storage of factors for // probability computation // Using charge but not time: // For hit PMT i, 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). // 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. // For a PMT i that was not hit, fPFac[i][j] is the probability of no // hit given j+1 photoelectrons. static unsigned int fXIdx[fNPMTsDim]; // Index from full PMT list into // list of selected (i.e., hit) PMTs DS::FitVertex fFitVertex; G4ThreeVector fEvtPosition; // event position from seed fitter TVector3 fTVecEvtPosition; // event position from seed fitter int fVIdx; // Index of volume of event position: 0 for // uniform fill; in partial fill, 0 for top and 1 for bottom 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[2]; // expected number of scintillation photons // produced per unit energy from database double fRefIdxAVFill[2]; // 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 fPulseWidth; // PE pulse width double fIntTime; // charge-integration time in ns double fLongIntTime; // long charge-integration time in ns double fAVInnerRadius; // inner AV radius static 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. Declared static // so it is carried over if new simulation is not done // for another fit of the same event. std::string fAVMaterialName[2]; // Material(s) inside the AV DU::LightPathCalculator fLightPath; // Light path calculator std::string fSelectorName; PMTSelectors::PMTSelector* fSelector; bool fSameGTID; // Following to provide for adjustments to be applied to fitted energy // to remove residual energy and radial dependence bool fApplyEAdj[2]; // is set true if energy adjustment is to be applied bool fApplyRAdj[2]; // is set true if radial adjustment is to be applied std::string fEAdjFnName[2]; // index of energy adjustment function in ratdb std::string fRAdjFnName[2]; // index of radial adjustment function in ratdb std::vector< std::vector > fAdjFnFitE, fAdjFnAdjE; std::vector< std::vector > fAdjFnFitR, fAdjFnRFactor; // For getting time usage of different PEnergy components Profiler* fProfileGenPhotons; //Profiler* fProfileAccumData; Profiler* fProfileCalcIncidProb; Profiler* fProfilePFactors; Profiler* fProfileOptimize; Profiler* fProfileFinish; Profiler* fProfileTotal; // Flag set true in BeginOfRun if grey-disc PMT model is being used bool fUseGreyDisc; // Options that can be changed by macro file inputs bool fUseChroma; // if the current fitter uses chroma for photon tracking 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[2];// 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[2]; // Minimum number of photons simulated for the fit std::string fRespTableName; // Name of PMT response table to be used double fEffScaleFactor; // overall photon detection efficiency // scaling factor double fEffExtraFactor[2]; // Additional photon detection efficiency // scaling factor that varies with fit type (dimension // equal 2 to provide for partial fill) bool fPartialFill; // normally false; set to true if detector is // filled partially. int fFitRegion; // Used only in partial fill; -1 = fit events only in // lower AV; +1 = fit events in upper AV; 0 = fit all double fFillZ; // Fill level for partial fill (in PSUP coordinates) TVector3 fAVPosition; // AV center in PSUP coordinates int fSetNPrimaries; // Set the number of primary electrons to be // generated in the auxiliary simulation // Note: this overrides setting of nPhotonSim double fSetEPrimary; // Set the energy of the primary electrons // 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 fPrevUseTime; static bool fPrevUseDirection; static int fPrevNPhotonSim; static int fPrevChargeType; static bool fGoodPrevAuxSim; static bool fPrevUseChroma; // Arrays & parameters used in accumulating auxiliary simulation statistics: static TVector3 fDirPMT[fNPMTsDim]; static std::vector< std::vector > fTCross; // Array for saving // times at which photons are incident on PMT's face static std::vector< std::vector > fTWeight; // Array for saving // weights to be applied to each time in calculating the pdfs. // Following used in AccumulateData to save info on track step from // the cavity to a PMT envelope G4int fPrevTrackID; std::string fPrevVolName; G4ThreeVector fPrePosition; G4ThreeVector fPostPosition; G4ThreeVector fPostPolarization; double fPostPtTime; double fEKin; double fMaxInclusionTime, fEstMinPropTime; int fNPhotonIncd; // Total number of photon incidences on PMTs // in auxiliary simulation #ifdef USE_CHROMA Chroma* fPtrChroma; // Pointer to Chroma object that sets up a chroma photon propagation interface #endif }; // class PEnergy } // namespace Methods } // namespace RAT #endif