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 class G4DataInterpolation;
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
0083   double n( void )                               const { return m_ReferenceRefractiveIndex; };
0084   double GetReferenceAttenuationLength( void )   const { return m_ReferenceAttenuationLength; };
0085
0086   void SetReferenceRefractiveIndex(double n)   { m_ReferenceRefractiveIndex   = n; };
0087   double GetReferenceRefractiveIndex(void)       const { return m_ReferenceRefractiveIndex; };
0088   void SetReferenceAttenuationLength(double l) { m_ReferenceAttenuationLength = l; };
0089
0090   void AddLogicalVolume(const G4LogicalVolume *volume) { m_LogicalVolumes.push_back(volume); };
0091
0092   const G4RadiatorMaterial *GetMaterial( void )  const { return m_Material; };
0093
0094   ParametricSurface *GetFrontSide(unsigned path) {
0095     if (m_Borders.find(path) == m_Borders.end()) return 0;
0096
0097     return dynamic_cast<ParametricSurface*>(m_Borders[path].first.GetObject());
0098   };
0099   ParametricSurface *GetRearSide(unsigned path) {
0100     if (m_Borders.find(path) == m_Borders.end()) return 0;
0101
0102     return dynamic_cast<ParametricSurface*>(m_Borders[path].second.GetObject());
0103   };
0104   std::map<unsigned, std::pair<TRef, TRef>> m_Borders;
0105
0106   // Material name in dd4hep world;
0107   void SetAlternativeMaterialName(const char *name) { m_AlternativeMaterialName = name; };
0108   const char *GetAlternativeMaterialName( void ) const {
0109     return m_AlternativeMaterialName.Data();
0110   };
0111
0112   void DisableOpticalPhotonGeneration( void ) { m_OpticalPhotonGenerationEnabled = false; };
0113   bool OpticalPhotonGenerationEnabled( void ) const { return m_OpticalPhotonGenerationEnabled; };
0114
0115   CherenkovRadiator *UseInRingImaging( void ) {
0116     m_UsedInRingImaging = true;
0117
0118     return this;
0119   };
0120   bool UsedInRingImaging( void ) const { return m_UsedInRingImaging; };
0121
0122  protected:
0123   // Run-time variables for the GEANT pass;
0124   //+const G4LogicalVolume *m_LogicalVolume;          //!
0125   std::vector<const G4LogicalVolume *> m_LogicalVolumes;          //!
0126   const G4RadiatorMaterial *m_Material;            //!
0127
0128  private:
0129   bool m_OpticalPhotonGenerationEnabled;
0130
0131   // Refractive index calculated for some fixed reference wave length (supposedly the average
0132   // one as seen on the detected photon wave length plot);
0133   double m_ReferenceRefractiveIndex;
0134   double m_ReferenceAttenuationLength;
0135
0136   TString m_AlternativeMaterialName;
0137
0138  public:
0139   void SetTrajectoryBinCount(unsigned bins) { m_TrajectoryBinCount = bins; };
0140   double GetTrajectoryBinCount( void) const { return m_TrajectoryBinCount; };
0141
0142   inline double GetSmearing( void ) const { return m_Smearing; };
0143   void SetGaussianSmearing(double sigma) { m_GaussianSmearing = true;  m_Smearing = sigma; }
0144   void SetUniformSmearing (double range) { m_GaussianSmearing = false; m_Smearing = range; }
0145   bool UseGaussianSmearing( void )  const { return m_GaussianSmearing; };
0146
0147   // FIXME: memory leak;
0148   void ResetLocations( void ) { m_Locations.clear(); }
0149   void ResetTimes( void ) { m_Times.clear(); }
0150   void AddLocation(/*unsigned sector,*/ const TVector3 &x, const TVector3 &n) {
0151     m_Locations/*[sector]*/.push_back(std::make_pair(x, n));
0152   };
0153   void AddTime(double value) { m_Times.push_back(value); };
0154
0155   // Transient variables for the ReconstructionFactory convenience;
0156   unsigned m_ID;                                            //!
0157   unsigned m_TrajectoryBinCount;                            //!
0158
0159   // This is a hack for now;
0160   double m_Smearing;                                        //!
0161   bool m_GaussianSmearing;                                  //!
0162   std::vector<std::pair<TVector3, TVector3>> m_Locations;   //!
0163   std::vector<double> m_Times;                              //!
0164
0165   std::vector<std::pair<double, double>> m_ri_lookup_table; //!
0166
0167   // Overall counts of "calibration" (did not pass QE check) and "real" (passed)
0168   // photons; since they originate from the same parent distribution, one can choose
0169   // basically any way to calculate their effective ratio for Poisson statistics purposes;
0170   unsigned m_CalibrationPhotonCount;                        //!
0171   unsigned m_DetectedPhotonCount;                           //!
0172   unsigned m_YieldStat;                                     //!
0173
0174   double m_YieldCff;
0175   double m_DetectedToCalibrationPhotonRatio;
0176   std::vector<CherenkovRadiatorCalibration> m_Calibrations;
0177
0178   bool m_UsedInRingImaging;                                 //!
0179
0180   CherenkovRadiator *InitializePlots(const char *tag) {
0181     m_Plots = new CherenkovRadiatorPlots(tag);
0182
0183     // Simplify scripting;
0184     return this;
0185   };
0186   CherenkovRadiatorPlots *Plots( void ) const { return m_Plots; };
0187   TCanvas *DisplayStandardPlots(const char *cname, const char *wname,
0188                 int wtopx, unsigned wtopy, unsigned wx, unsigned wy) const;
0189
0190   CherenkovRadiatorPlots  *m_Plots;                         //!
0191
0192   G4DataInterpolation *m_RefractiveIndex;                   //!
0193
0194   bool m_OutputPlotVisualizationEnabled;                    //!
0195   int m_wtopx;                                              //!
0196   unsigned m_wtopy;                                         //!
0197   unsigned m_wx;                                            //!
0198   unsigned m_wy;                                            //!
0199   //TCanvas *m_cv;                                            //!
0200
0201 #ifndef DISABLE_ROOT_IO
0202   ClassDef(CherenkovRadiator, 9);
0203 #endif
0204 };
0205
0206 } // namespace IRT2