Functions | Variables
garching_flux_to_genie_root Namespace Reference

Functions

def Alpha (ebar, ebar2)
 
def SpectralFunction (N, ebar, alpha, elo=energyBins['low'], ehi=energyBins['high'], binsize=(energyBins['high']-energyBins['low'])/energyBins['bins'])
 
def main ()
 MAIN FUNCTION. More...
 

Variables

dictionary energyBins
 
dictionary timeBins
 

Function Documentation

def garching_flux_to_genie_root.Alpha (   ebar,
  ebar2 
)
Compute the alpha parameter given the first two energy moments.

Definition at line 39 of file garching_flux_to_genie_root.py.

Referenced by main().

39 def Alpha(ebar, ebar2):
40  """Compute the alpha parameter given the first two energy moments."""
41 
42  if (ebar**2 == ebar2):
43  return 0
44  else:
45  return (ebar2 - 2 * ebar**2) / (ebar**2 - ebar2)
46 
47 
48 # /////////////////////////////////////////////////////////////////////////////
def garching_flux_to_genie_root.main ( )

MAIN FUNCTION.

1D plots hits vs w //

Definition at line 82 of file garching_flux_to_genie_root.py.

References Alpha(), g4zmq.apply(), Divide(), check_time_usage.float, novadaq::HexUtils.format(), print, utilities.progressBar(), PandAna.Demos.demo1.range, SetBinContent(), and SpectralFunction().

82 def main():
83 
84  channels = {
85  'fluence_e': 'data/fluxes/neutrino_signal_nue-LS220-z9.6co.dat',
86  'fluence_ae': 'data/fluxes/neutrino_signal_nuebar-LS220-z9.6co.dat',
87  'fluence_x': 'data/fluxes/neutrino_signal_nux-LS220-z9.6co.dat',
88  'fluence_ax': 'data/fluxes/neutrino_signal_nuxbar-LS220-z9.6co.dat'
89  }
90 
91  channelTitles = {
92  'fluence_e': 'nue flux, {} nu/MeV/s'.format(constant.LUM_UNITS),
93  'fluence_ae': 'nuebar flux, {} nu/MeV/s'.format(constant.LUM_UNITS),
94  'fluence_x': 'nux flux, {} nu/MeV/s'.format(constant.LUM_UNITS),
95  'fluence_ax': 'nuxbar flux, {} nu/MeV/s'.format(constant.LUM_UNITS)
96  }
97 
98  rootOutFilePath = 'test.root'
99  rootFile = TFile(rootOutFilePath, 'RECREATE')
100 
101  # Loop over time and energy and fill the histograms
102  model, hists = {}, {}
103  for channel, chanFile in channels.items():
104  hTime = TH1F("time", "time", timeBins['bins'], timeBins['low'], timeBins['high'])
105  hEnergy = TH1F("energy", "energy", energyBins['bins'], energyBins['low'], energyBins['high'])
106  hCounts = TH2F("counts", "counts", timeBins['bins'], timeBins['low'], timeBins['high'], energyBins['bins'], energyBins['low'], energyBins['high'])
107 
108  hists[channel] = TH2F(channel, channelTitles[channel], timeBins['bins'], timeBins['low'], timeBins['high'], energyBins['bins'], energyBins['low'], energyBins['high'])
109 
110  # UNITS: [time= s, luminosity_foe= foe/s, <e>= MeV, <e^2>= MeV^2]
111  model[channel] = pd.read_csv(chanFile, sep='\s+', comment='#')
112  model[channel].columns = ['time', 'luminosity_foe', '<e>', '<e^2>']
113 
114  # Custom columns
115  model[channel]['luminosity_MeV'] = model[channel]['luminosity_foe'] * constant.MEV_PER_FOE / constant.LUM_UNITS
116  model[channel]['alpha'] = model[channel].apply(lambda row: Alpha(row['<e>'], row['<e^2>']), axis=1)
117  model[channel]['N'] = model[channel]['luminosity_MeV'] / model[channel]['<e>']
118 
119  # Fill any NaNs with 0
120  model[channel].fillna(0, inplace=True)
121 
122  N_timesteps = len(model[channel]['time'])
123  energyRange = np.arange(energyBins['low'], energyBins['high'], (energyBins['high']-energyBins['low'])/energyBins['bins'])
124  for timeIdx in range(1, N_timesteps):
125  thisTime = model[channel]['time'][timeIdx]
126  timeBin = hTime.FindBin(thisTime)
127  (timeIdx % 100 == 0 or timeIdx == N_timesteps - 1) and utilities.progressBar(timeIdx, N_timesteps, msg=channel)
128 
129  spectrum = SpectralFunction(model[channel]['N'][timeIdx], model[channel]['<e>'][timeIdx], model[channel]['alpha'][timeIdx])
130 
131  for energyIdx in energyRange:
132  energyBin = hEnergy.FindBin(energyIdx)
133  currentValue = hists[channel].GetBinContent(timeBin, energyBin)
134  newValue = currentValue + spectrum[energyBin - 1]
135  hists[channel].SetBinContent(timeBin, energyBin, float(newValue))
136  hCounts.SetBinContent(timeBin, energyBin, hCounts.GetBinContent(timeBin, energyBin) + 1)
137 
138  hists[channel].Divide(hCounts)
139 
140  # We don't want to write these to the file, so clear them
141  hTime = None
142  hEnergy = None
143  hCounts = None
144 
145  print('')
146 
147 
148  rootFile.Write()
149  rootFile.Close()
150 
151  return
152 
153 
154 # .............................................................................
def apply(command)
Definition: g4zmq.py:38
def progressBar(value, endvalue, bar_length=40, msg='Progress')
Definition: utilities.py:13
ratio_hxv Divide(hxv, goal_hxv)
def SpectralFunction(N, ebar, alpha, elo=energyBins['low'], ehi=energyBins['high'], binsize=(energyBins['high']-energyBins['low'])/energyBins['bins'])
bool print
std::string format(const int32_t &value, const int &ndigits=8)
Definition: HexUtils.cpp:14
cout<< "--"<< endl;for(Int_t iP=1;iP<=hyz->GetNbinsX();iP++){for(Int_t iC=1;iC<=hyz->GetNbinsY();iC++){if(hyv->GetBinContent(iP, iC)>-999){goal_hyv-> SetBinContent(iP, iC,-(dy[iP-1][iC-1]))
def garching_flux_to_genie_root.SpectralFunction (   N,
  ebar,
  alpha,
  elo = energyBins['low'],
  ehi = energyBins['high'],
  binsize = (energyBins['high']-energyBins['low'])/energyBins['bins'] 
)
Compute the spectral shape, given the necessary parameters.

  returns a numpy array of values of length binsize.

  :param N:       Number of neutrinos, constant scaling factor.
  :param ebar:    Average neutrino energy.
  :param alpha:   Alpha parameter
  :param elo:     Min energy (default taken from energyBins)
  :param ehi:     Max energy (default taken from timeBins)
  :param binsize: Energy binning

Definition at line 49 of file garching_flux_to_genie_root.py.

Referenced by main().

49 def SpectralFunction(N, ebar, alpha, elo=energyBins['low'], ehi=energyBins['high'], binsize=(energyBins['high']-energyBins['low'])/energyBins['bins']):
50  """Compute the spectral shape, given the necessary parameters.
51 
52  returns a numpy array of values of length binsize.
53 
54  :param N: Number of neutrinos, constant scaling factor.
55  :param ebar: Average neutrino energy.
56  :param alpha: Alpha parameter
57  :param elo: Min energy (default taken from energyBins)
58  :param ehi: Max energy (default taken from timeBins)
59  :param binsize: Energy binning
60  """
61 
62  if (ebar == 0):
63  return np.zeros(energyBins['bins'])
64 
65  energies = np.arange(elo, ehi, binsize)
66  spectrum = np.array([])
67 
68  # Nomalization factors (energy independent)
69  # norm1 -- area = 1
70  # norm2 -- instantaneous neutrino rate (luminosity/average energy)
71  norm1 = (alpha + 1)**(alpha + 1) / (ebar * math.gamma(alpha + 1))
72  norm2 = N
73 
74  for energy in energies:
75  shape = (energy / ebar)**alpha * math.exp(-(alpha + 1) * energy / ebar)
76  spectrum = np.append(spectrum, norm1 * norm2 * shape)
77 
78  return spectrum
79 
80 
81 # /////////////////////////////////////////////////////////////////////////////
def SpectralFunction(N, ebar, alpha, elo=energyBins['low'], ehi=energyBins['high'], binsize=(energyBins['high']-energyBins['low'])/energyBins['bins'])

Variable Documentation

dictionary garching_flux_to_genie_root.energyBins
Initial value:
1 = {
2  'low': 0.,
3  'high': 100.,
4  'bins': 500
5 }

Definition at line 25 of file garching_flux_to_genie_root.py.

dictionary garching_flux_to_genie_root.timeBins
Initial value:
1 = {
2  'low': -0.5,
3  'high': 15.,
4  'bins': 3100
5 }

Definition at line 31 of file garching_flux_to_genie_root.py.