garching_flux_to_genie_root.py
Go to the documentation of this file.
1 #!/usr/bin/env python
2 
3 # //////////////////////////////////////////////////////////////////////////////
4 # // CONVERT RAW GARCHING FLUX FILES TO ROOT FILES
5 # //
6 # // \brief a tool to convert Garching flux files into a ROOT file that
7 # // GENIE can read.
8 # // \author Justin Vasel <jvasel@indiana.edu>
9 # // \date 14 September 2018
10 # //
11 # //////////////////////////////////////////////////////////////////////////////
12 
13 import numpy as np
14 import pandas as pd
15 from ROOT import TFile, TTree, TH2F, TH1F, TGraph, TCanvas, gROOT
16 gROOT.SetBatch(True)
17 
18 import math
19 import sys
20 
21 import constants as constant
22 import utilities as utilities
23 
24 # Some options
25 energyBins = {
26  'low': 0.,
27  'high': 100.,
28  'bins': 500
29 }
30 
31 timeBins = {
32  'low': -0.5,
33  'high': 15.,
34  'bins': 3100
35 }
36 
37 
38 # /////////////////////////////////////////////////////////////////////////////
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 # /////////////////////////////////////////////////////////////////////////////
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 # /////////////////////////////////////////////////////////////////////////////
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 # .............................................................................
155 if __name__ == '__main__':
156  main()
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]))