/////////////////////////////////////////////////////////////////////// // // Energy fit based on prompt hits using a seed position, direction and // initial energy estimate. // // Author: name John Walker // // REVISION HISTORY: // 28/09/2015 : J Walker - New file // 01/10/2015 : J Walker - Small changes to ensure positive number // of Cerenkov photons // 18/01/2016 : J Walker - Only perform RSP if inside AV // 02/03/2016 : J Walker - Set default energy errors // 09/03/2016 : M Mottram - Fix compiler warnings // 13/03/2016 : J Walker - Set iterations and tolerance in ratdb // and optimisation // 30/08/2016 : M Mottram - fixed rad2degree conversion bug, added // diagnostics FOMs // 16/12/2016 : I Coulter - Fixed bug in probability of late // scattered photons at request of JW // 20/12/2017 : Logan L - Add several figures of merit // // This method calculates the optical response of each PMT (taking into // account PMT efficiency, solid angle, Cerenkov angular distribution, // transmission probability, attenuation etc.) to estimate the number // of Cerenkov photons from the prompt hits. A lookup table is used to // map from estimated Cerenkov photons to energy. This fitter is only // for a detector geometry with light water in the inner av. // /////////////////////////////////////////////////////////////////////// #ifndef __RAT_Method_EnergyRSP__ #define __RAT_Method_EnergyRSP__ #include #include #include #include #include #include #include #include #include namespace RAT { namespace Methods { class EnergyRSP : public SeededMethod { public: // Returns the method's name virtual std::string GetName() const { return EnergyRSP::Name(); } // Returns the method's name static std::string Name() { return std::string( "energyRSP" ); } // Initialise the method // // "param" is an optional index in the db to use void Initialise( const std::string& param ); void SetI( const std::string& param, const int value ); // This will load the tables from the database void BeginOfRun( DS::Run& run ); void EndOfRun( DS::Run& run ); // Calculate and return the best fit for the seed virtual DS::FitResult GetBestFit(); // Set the seed to the default void DefaultSeed(); private: // Return the wavelength independent PMT response // // "pmtPos" is the PMT position // "pmtDir" is the PMT direction // "pmtID" is the PMT logical channel number double wiPMTResponse( const TVector3& pmtPos, const TVector3& pmtDir, const int pmtID ); // Return the PMT efficiency // // "energy" is the energy (MeV) of the photon // "CosThetan" is the angle with the PMT direction of the photon // "pmtID" is the PMT logical channel number double PMTEfficiency( const double energy, const double CosThetan, const int pmtID ); // Return the Cerenkov angular distribution // // "pmtID" is the PMT logical channel number double CerenkovAngularDist( const int pmtID ); // Fill attenuation values for each region along with rayleigh scattering // that results in late light // // "energy" is the energy (MeV) of the photon // "pmtID" is the PMT logical channel number double Attenuation( const double energy, const int pmtID, double& absInnerAV, double& absAV, double& absCavity, double& rayleigh); // Corrects for multiple photons triggering the PMT // // "nPhotons" is the estimated number of photons incident on an individual PMT double MultiPhotonCorrector( const double nPhotons ); // Performs a 1D interpolation using the TSpline3 method // // "inValue" is the input value // "inTable" is a table of values corresponding to the input // "outTable" is a table of corresponding coordinates to "inTable" values double Interpolate1D( const double inValue, const std::vector& inTable, const std::vector& outTable ); // Returns probability of rayleigh scattering resulting in late (out of window) // light from a 2D lookup table double RayleighLookup( const double radius, const double eventDotLightVec); // Generate a vector of interpolated Cerenkov angular distribution std::vector GenerateCerenkovSpline( const double energy ); // Load absorption pathlength and refractive index vectors for a given material void LoadOpticsVectors(const std::string& materialName, G4MaterialPropertyVector** abslength, G4MaterialPropertyVector** rindex); // Compare values for sorting static int Compare( const void * a, const void * b ) { float fa = *(const float*) a; float fb = *(const float*) b; return (fa > fb) - (fa < fb); } std::string fIndex; // Optional index override for the db table - unused at the moment bool fUseRefracted; // Use the refracted light path for all calculations bool fStoreDiagnostics; // Store a breakdown of the contributions for different components of RSP unsigned int fMaxIterations; // Maximum iterations of method double fTolerance; // Maximum percentage differance between effective and predicted PMT hits double fRMax; // Maximum radius for valid fits std::vector fMeVValues; // The lookup table for energies std::vector fEventDirDotInitialLightVecValues; // The lookup table for the dot product between event direction and initial light vector size_t fDirectionDotVecSize; double fDirectionDotVecWidth; double fDirectionDotVecMin; double fDirectionDotVecMax; std::vector fRadialValues; // The lookup table for the centre of the radial bins size_t fRadialSize; double fRadialWidth; double fRadialMin; double fRadialMax; std::vector fAngularValues; // The lookup table for angular values std::vector fNphotonValues; // The lookup table for Cerenkov photons per MeV std::vector fCerenkovAngularDist; // The lookup table for the Cerenkov angular distribution std::vector fRayleighLateProb; // The lookup table for the probability that a Rayleigh scattered photon will be late std::vector fAngularResponse; // The lookup table for the PMT angular response std::vector fPredictedNphotonValues; // The lookup table for predicted Cerenkov photons per MeV TSpline3* fAngularResponseSpline; double fPhotonScalingFactor; // Scales the number of photons in the predicted table to match the actual number of Cherenkov photons in an event bool fPredictedTableExists; // Flag true if predicted Cerenkov photons per MeV table exists double fCollectionEfficiency; // PMT collection efficiency G4MaterialPropertyVector* fPEEscapeProbability; // Photoelectron escape probability G4MaterialPropertyVector* fInnerAV_ABSLENGTH; // Water absorption length G4MaterialPropertyVector* fInnerAV_RSLENGTH; // Water Rayleigh scattering length G4MaterialPropertyVector* fAV_ABSLENGTH; // Acrylic absorption length G4MaterialPropertyVector* fAV_RSLENGTH; // Acrylic Rayleigh scattering length G4MaterialPropertyVector* fCavity_ABSLENGTH; // Water absorption length G4MaterialPropertyVector* fCavity_RSLENGTH; // Water Rayleigh scattering length double fInnerAVRadius; // Inner AV radius size_t fNumChannels; // Total number of channels double fNoiseRate; // PMT noise rate PMTSelectors::PMTSelector* fTimeResidualCut; // Time residual cut PMT selector PMTSelectors::PMTSelector* fPMTCalSelector; // PMTCalSelector ChannelEfficiency* fChannelEfficiency; // Channel Efficiency Calculator PMTVariation* fPMTVariation; // PMT relative efficiency DU::LightPathCalculator fLightPath; // Light path std::vector fDistInInnerAV; // Vectors of distances for each PMT std::vector fDistInAV; std::vector fDistInWater; std::vector fCosThetan; // Vector of costheta to PMT direction for each PMT std::vector fInitialLightVec; // Vector of initial light vector for each PMT std::vector fCerenkovSpline; // Interpolation of Cerenkov angular distribution TVector3 fEventPosition; // Event position seed TVector3 fEventDirection; // Event direction seed double fSeedEnergy; // Event energy seed double twHighLimit; double twLowLimit; }; } // namespace Methods } // namespace RAT #endif // __RAT_Method_EnergyRSP__