TFluka/TFlukaCerenkov.cxx

   1 #include "TFlukaCerenkov.h"
   2
   3 Double_t TFlukaCerenkov::fgGlobalMaximumEfficiency = 0.;
   4
   5 ClassImp(TFlukaCerenkov);
   6
   7
   8 TFlukaCerenkov::TFlukaCerenkov()
   9     : fSamples(0),
  10       fIsMetal(kFALSE),
  11       fIsSensitive(kFALSE),
  12       fEnergy(0),
  13       fWaveLength(0),
  14       fAbsorptionCoefficient(0),
  15       fQuantumEfficiency(0),
  16       fRefractionIndex(0),
  17       fMaximumEfficiency(0.)
  18 {
  19 // Default constructor
  20 }
  21
  22
  23 TFlukaCerenkov::TFlukaCerenkov(Int_t npckov, Float_t *ppckov, Float_t *absco, Float_t *effic, Float_t *rindex)
  24 {
  25 // Constructor
  26     fSamples = npckov;
  27     fEnergy                = new Float_t[fSamples];
  28     fWaveLength            = new Float_t[fSamples];
  29     fAbsorptionCoefficient = new Float_t[fSamples];
  30     fRefractionIndex       = new Float_t[fSamples];
  31     fQuantumEfficiency     = new Float_t[fSamples];
  32
  33     for (Int_t i = 0; i < fSamples; i++) {
  34         fEnergy[i]             = ppckov[i];
  35         fWaveLength[i]         = khc / ppckov[i];
  36         if (absco[i] > 0.) {
  37             fAbsorptionCoefficient[i]   = 1./ absco[i];
  38         } else {
  39             fAbsorptionCoefficient[i]   = 1.e15;
  40         }
  41         fRefractionIndex[i]    = rindex[i];
  42         fQuantumEfficiency[i]  = effic[i];
  43         //
  44         // Find local maximum quantum efficiency
  45         if (effic[i] > fMaximumEfficiency) fMaximumEfficiency = effic[i];
  46         //
  47         // Flag is sensitive if quantum efficiency 0 < eff < 1 for at least one value.
  48         if (effic[i] < 1. && effic[i] > 0.) fIsSensitive = 1;
  49     }
  50     // Find global  maximum quantum efficiency
  51     if (fMaximumEfficiency > GetGlobalMaximumEfficiency()) {
  52         SetGlobalMaximumEfficiency(fMaximumEfficiency);
  53     }
  54     printf("Maximum eff. %f\n",  GetGlobalMaximumEfficiency());
  55
  56
  57 }
  58
  59 Float_t TFlukaCerenkov::GetAbsorptionCoefficient(Float_t energy)
  60 {
  61 //
  62 // Get AbsorptionCoefficient for given energy
  63 //
  64     return Interpolate(energy, fEnergy, fAbsorptionCoefficient);
  65
  66 }
  67
  68 Float_t TFlukaCerenkov::GetRefractionIndex(Float_t energy)
  69 {
  70 //
  71 // Get RefractionIndex for given energy
  72 //
  73     return Interpolate(energy, fEnergy, fRefractionIndex);
  74
  75 }
  76
  77 Float_t TFlukaCerenkov::GetQuantumEfficiency(Float_t energy)
  78 {
  79 //
  80 // Get QuantumEfficiency for given energy
  81 //
  82     return Interpolate(energy, fEnergy, fQuantumEfficiency);
  83
  84 }
  85
  86
  87 Float_t TFlukaCerenkov::GetAbsorptionCoefficientByWaveLength(Float_t wavelength)
  88 {
  89 //
  90 // Get AbsorptionCoefficient for given wavelength
  91 //
  92     Float_t energy = khc / wavelength;
  93     return Interpolate(energy, fEnergy, fAbsorptionCoefficient);
  94
  95 }
  96
  97 Float_t TFlukaCerenkov::GetRefractionIndexByWaveLength(Float_t wavelength)
  98 {
  99 //
 100 // Get RefractionIndex for given wavelenth
 101 //
 102     Float_t energy = khc / wavelength;
 103     return Interpolate(energy, fEnergy, fRefractionIndex);
 104 }
 105
 106 Float_t TFlukaCerenkov::GetQuantumEfficiencyByWaveLength(Float_t wavelength)
 107 {
 108 //
 109 // Get QuantumEfficiency for given wavelength
 110 //
 111     Float_t energy = khc / wavelength;
 112     return Interpolate(energy, fEnergy, fQuantumEfficiency);
 113 }
 114
 115 Float_t TFlukaCerenkov::Interpolate(Float_t value, Float_t* array1, Float_t* array2)
 116 {
 117 //
 118 // Interpolate array values
 119 //
 120     if (value < array1[0] && value >= array1[fSamples - 1]) {
 121         Warning("Interpolate()", "Photon energy out of range. Returning 0.");
 122         return (0.);
 123     }
 124
 125     Int_t i = TMath::BinarySearch(fSamples, array1, value);
 126     if (i == (fSamples-1)) {
 127         return (array2[i]);
 128     } else {
 129         return (array2[i] + (array2[i+1] - array2[i]) / (array1[i+1] - array1[i]) * (value - array1[i]));
 130     }
 131 }
 132