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