include/IRT2/CherenkovRadiator.h

0001 #pragma once
0002
0003 #include <map>
0004 #include <vector>
0005
0006 #include <TRef.h>
0007 #include <TObject.h>
0008 #include <TVector3.h>
0009 #include <TString.h>
0010 class TH1D;
0011
0012 #include "ParametricSurface.h"
0013 class G4LogicalVolume;
0014 class G4RadiatorMaterial;
0015
0016 class G4DataInterpolation;
0017
0018 namespace IRT2 {
0019
0020 struct CherenkovRadiatorCalibration: public TObject {
0021   CherenkovRadiatorCalibration(): m_Stat(0),
0022     m_AverageZvtx(0.0), m_hcalib(0), m_Coffset(0.0), m_Csigma(0.0) {};
0023   ~CherenkovRadiatorCalibration() {};
0024
0025   unsigned m_Stat;
0026   double /*m_AverageRefractiveIndex,*/ m_AverageZvtx;
0027   // Averaged across photons produced in this particular radiator (and for parent particles
0028   // in a given theta bin), calculated for ALL registered radiators (will be needed for
0029   // IRT algorithm to function properly);
0030   std::vector<double> m_AverageRefractiveIndices;
0031
0032   TH1D *m_hcalib;
0033   double m_Coffset, m_Csigma;
0034
0035 #ifndef DISABLE_ROOT_IO
0036   ClassDef(CherenkovRadiatorCalibration, 1);
0037 #endif
0038 };
0039
0040 struct CherenkovRadiatorPlots {
0041   CherenkovRadiatorPlots(const char *tag);
0042   ~CherenkovRadiatorPlots() {};
0043
0044   void SetRefractiveIndexRange(double min, double max);
0045   void SetPhotonVertexRange(double min, double max);
0046   void SetCherenkovAngleRange(double min, double max);
0047
0048   // Monte-Carlo plots;
0049   TH1D *hvtx()              const { return m_hvtx;  };
0050   TH1D *hnpe()              const { return m_hnpe;  };
0051   TH1D *hwl()               const { return m_hwl;  };
0052   TH1D *hri()               const { return m_hri;  };
0053
0054   // Reconstruction plots;
0055   TH1D *hnhits()            const { return m_hnhits;  };
0056   TH1D *hthph()             const { return m_hthph;  };
0057   TH1D *hccdfph()           const { return m_hccdfph;  };
0058   TH1D *hthtr()             const { return m_hthtr;  };
0059
0060 private:
0061   std::string m_Tag;
0062   TH1D *m_hvtx, *m_hnpe, *m_hwl, *m_hri, *m_hnhits, *m_hthph, *m_hccdfph, *m_hthtr;
0063 };
0064
0065 class CherenkovRadiator: public TObject {
0066  public:
0067   // NB: do not want to use physical volume here because a particle can cross more than one of them
0068   // (at the sector boundary), while there is no good reason to separate these contributions;
0069  CherenkovRadiator(const G4LogicalVolume *volume = 0, const G4RadiatorMaterial *material = 0):
0070    /*m_LogicalVolume(volume),*/
0071    m_Material(material), m_OpticalPhotonGenerationEnabled(true),
0072    m_ReferenceRefractiveIndex(0.0), m_ReferenceAttenuationLength(0.0),
0073    m_ID(0), m_TrajectoryBinCount(1), m_Smearing(0.0),
0074    m_GaussianSmearing(false), m_CalibrationPhotonCount(0), m_DetectedPhotonCount(0), m_YieldStat(0),
0075    m_YieldCff(0.0), m_DetectedToCalibrationPhotonRatio(0.0),
0076    m_UsedInRingImaging(false), m_Plots(0) {
0077     m_LogicalVolumes.push_back(volume);
0078   };
0079   ~CherenkovRadiator() {};
0080
0081   double n( void )                               const { return m_ReferenceRefractiveIndex; };
0082   double GetReferenceAttenuationLength( void )   const { return m_ReferenceAttenuationLength; };
0083
0084   void SetReferenceRefractiveIndex(double n)   { m_ReferenceRefractiveIndex   = n; };
0085   double GetReferenceRefractiveIndex(void)       const { return m_ReferenceRefractiveIndex; };
0086   void SetReferenceAttenuationLength(double l) { m_ReferenceAttenuationLength = l; };
0087
0088   void AddLogicalVolume(const G4LogicalVolume *volume) { m_LogicalVolumes.push_back(volume); };
0089
0090   const G4RadiatorMaterial *GetMaterial( void )  const { return m_Material; };
0091
0092   ParametricSurface *GetFrontSide(unsigned path) {
0093     if (m_Borders.find(path) == m_Borders.end()) return 0;
0094
0095     return dynamic_cast<ParametricSurface*>(m_Borders[path].first.GetObject());
0096   };
0097   ParametricSurface *GetRearSide(unsigned path) {
0098     if (m_Borders.find(path) == m_Borders.end()) return 0;
0099
0100     return dynamic_cast<ParametricSurface*>(m_Borders[path].second.GetObject());
0101   };
0102   std::map<unsigned, std::pair<TRef, TRef>> m_Borders;
0103
0104   // Material name in dd4hep world;
0105   void SetAlternativeMaterialName(const char *name) { m_AlternativeMaterialName = name; };
0106   const char *GetAlternativeMaterialName( void ) const {
0107     return m_AlternativeMaterialName.Data();
0108   };
0109
0110   void DisableOpticalPhotonGeneration( void ) { m_OpticalPhotonGenerationEnabled = false; };
0111   bool OpticalPhotonGenerationEnabled( void ) const { return m_OpticalPhotonGenerationEnabled; };
0112
0113   CherenkovRadiator *UseInRingImaging( void ) {
0114     m_UsedInRingImaging = true;
0115
0116     return this;
0117   };
0118   bool UsedInRingImaging( void ) const { return m_UsedInRingImaging; };
0119
0120  protected:
0121   // Run-time variables for the GEANT pass;
0122   //+const G4LogicalVolume *m_LogicalVolume;          //!
0123   std::vector<const G4LogicalVolume *> m_LogicalVolumes;          //!
0124   const G4RadiatorMaterial *m_Material;            //!
0125
0126  private:
0127   bool m_OpticalPhotonGenerationEnabled;
0128
0129   // Refractive index calculated for some fixed reference wave length (supposedly the average
0130   // one as seen on the detected photon wave length plot);
0131   double m_ReferenceRefractiveIndex;
0132   double m_ReferenceAttenuationLength;
0133
0134   TString m_AlternativeMaterialName;
0135
0136  public:
0137   void SetTrajectoryBinCount(unsigned bins) { m_TrajectoryBinCount = bins; };
0138   double GetTrajectoryBinCount( void) const { return m_TrajectoryBinCount; };
0139
0140   inline double GetSmearing( void ) const { return m_Smearing; };
0141   void SetGaussianSmearing(double sigma) { m_GaussianSmearing = true;  m_Smearing = sigma; }
0142   void SetUniformSmearing (double range) { m_GaussianSmearing = false; m_Smearing = range; }
0143   bool UseGaussianSmearing( void )  const { return m_GaussianSmearing; };
0144
0145   // FIXME: memory leak;
0146   void ResetLocations( void ) { m_Locations.clear(); }
0147   void ResetTimes( void ) { m_Times.clear(); }
0148   void AddLocation(/*unsigned sector,*/ const TVector3 &x, const TVector3 &n) {
0149     m_Locations/*[sector]*/.push_back(std::make_pair(x, n));
0150   };
0151   void AddTime(double value) { m_Times.push_back(value); };
0152
0153   // Transient variables for the ReconstructionFactory convenience;
0154   unsigned m_ID;                                            //!
0155   unsigned m_TrajectoryBinCount;                            //!
0156
0157   // This is a hack for now;
0158   double m_Smearing;                                        //!
0159   bool m_GaussianSmearing;                                  //!
0160   std::vector<std::pair<TVector3, TVector3>> m_Locations;   //!
0161   std::vector<double> m_Times;                              //!
0162
0163   std::vector<std::pair<double, double>> m_ri_lookup_table; //!
0164
0165   // Overall counts of "calibration" (did not pass QE check) and "real" (passed)
0166   // photons; since they originate from the same parent distribution, one can choose
0167   // basically any way to calculate their effective ratio for Poisson statistics purposes;
0168   unsigned m_CalibrationPhotonCount;                        //!
0169   unsigned m_DetectedPhotonCount;                           //!
0170   unsigned m_YieldStat;                                     //!
0171
0172   double m_YieldCff;
0173   double m_DetectedToCalibrationPhotonRatio;
0174   std::vector<CherenkovRadiatorCalibration> m_Calibrations;
0175
0176   bool m_UsedInRingImaging;                                 //!
0177
0178   CherenkovRadiator *InitializePlots(const char *tag) {
0179     m_Plots = new CherenkovRadiatorPlots(tag);
0180
0181     // Simplify scripting;
0182     return this;
0183   };
0184   CherenkovRadiatorPlots *Plots( void ) const { return m_Plots; };
0185   void DisplayStandardPlots(const char *cname, int wtopx, unsigned wtopy, unsigned wx, unsigned wy) const;
0186
0187   CherenkovRadiatorPlots  *m_Plots;                         //!
0188
0189   G4DataInterpolation *m_RefractiveIndex;                   //!
0190
0191 #ifndef DISABLE_ROOT_IO
0192   ClassDef(CherenkovRadiator, 9);
0193 #endif
0194 };
0195
0196 } // namespace IRT2