Proton energy array in GeV
199050
def _Ep(self):
        
        return np.logspace(
            np.log10(self.Epmin.to("GeV").value),
            np.log10(self.Epmax.to("GeV").value),
            int(self.nEpd * (np.log10(self.Epmax / self.Epmin))),
        )
630526
def _spectrum(self, photon_energy):
        

        Eph = _validate_ene(photon_energy)

        spec = self.n0 * (
            self.weight_ee * self._emiss_ee(Eph)
            + self.weight_ep * self._emiss_ep(Eph)
        )

        return spec
52095
def pssm_array2pwm_array(arr, background_probs=DEFAULT_BASE_BACKGROUND):
    
    b = background_probs2array(background_probs)
    b = b.reshape([1, 4, 1])
    return (np.exp(arr) * b).astype(arr.dtype)
223663
def main():
    
    usage = "usage: %(prog)s [options] "
    description = "Merge a set of Fermi-LAT files."

    parser = argparse.ArgumentParser(usage=usage, description=description)

    parser.add_argument('-o', '--output', default=None, type=str,
                        help='Output file.')
    parser.add_argument('--clobber', default=False, action='store_true',
                        help='Overwrite output file.')
    parser.add_argument('--hdu', default=None, type=str,
                        help='HDU name.')
    parser.add_argument('--gzip', action='store_true',
                        help='Compress output file')
    parser.add_argument('--rm', action='store_true',
                        help='Remove input files.')
    parser.add_argument('files', nargs='+', default=None,
                        help='List of input files.')

    args = parser.parse_args()

    hpx_map = merge_utils.stack_energy_planes_hpx(args.files, hdu=args.hdu)

    if args.output:
        hpx_map.hpx.write_fits(hpx_map.counts, args.output,
                               extname=args.hdu, clobber=args.clobber)

        if args.gzip:
            os.system('gzip -9 %s' % args.output)

        if args.rm:
            for farg in args.files:
                flist = glob.glob(farg)
                for ffound in flist:
                    os.path.unlink(ffound)
138519
def energy(self):
        """"""
        result = 0.0
        for index1 in range(self.numc):
            for index2 in range(index1):
                if self.scaling[index1, index2] > 0:
                    for se, ve in self.yield_pair_energies(index1, index2):
                        result += se*ve*self.scaling[index1, index2]
        return result
76587
def _one_body_mapping(a_i, a_j, threshold=0.000001):
        
        pauli_list = []
        for alpha in range(2):
            for beta in range(2):
                pauli_prod = Pauli.sgn_prod(a_i[alpha], a_j[beta])
                coeff = 1.0/4 * pauli_prod[1] * np.power(-1j, alpha) * np.power(1j, beta)
                pauli_term = [coeff, pauli_prod[0]]
                if np.absolute(pauli_term[0]) > threshold:
                    pauli_list.append(pauli_term)
        return Operator(paulis=pauli_list)
220724
def _pdist(p):
        
        index, ref, ampl, cutoff, beta = p[:5]
        if cutoff == 0.0:
            pdist = models.PowerLaw(
                ampl * 1e30 * u.Unit("1/eV"), ref * u.TeV, index
            )
        else:
            pdist = models.ExponentialCutoffPowerLaw(
                ampl * 1e30 * u.Unit("1/eV"),
                ref * u.TeV,
                index,
                cutoff * u.TeV,
                beta=beta,
            )
        return pdist
268253
def _get_next_events(self, material):
        """"""

        f = ROOT.TFile(self.filenames[material])
        try:
            t = f.Get('gst')
            n = t.GetEntries()
        except:
            self.log.critical('Could not open the ROOT file with Genie events')
            raise

        for i in range(n):
            t.GetEntry(i)
            next_events = []

            position = convert_3vector_to_dict([self.particle['position']['x'].get_cache(),
                                                self.particle['position']['y'].get_cache(),
                                                self.particle['position']['z'].get_cache()])

            lepton_event = {}
            if t.El ** 2 - (t.pxl ** 2 + t.pyl ** 2 + t.pzl ** 2) < 1e-7:
                lepton_event['pid'] = self.particle[
                                      'pid'].get()  # Either NC or ES
            else:
                lepton_event['pid'] = lookup_cc_partner(
                    self.particle['pid'].get())

            # units: GeV -> MeV
            momentum_vector = [1000 * x for x in [t.pxl, t.pyl, t.pzl]]

            lepton_event['momentum'] = convert_3vector_to_dict(momentum_vector)

            lepton_event['position'] = position

            next_events.append(lepton_event)

            for j in range(t.nf):  # nf, number final hadronic states
                hadron_event = {}
                hadron_event['pid'] = t.pdgf[j]
                hadron_event['position'] = position

                # units: GeV -> MeV
                momentum_vector = [1000 * x for x in [t.pxf[j], t.pyf[j], t.pzf[j]]]

                hadron_event['momentum'] = convert_3vector_to_dict(momentum_vector)

                next_events.append(hadron_event)

            event_type = {}

            event_type['vertex'] = position

            to_save = {} # maps our names to Genie gst names
            to_save['incoming_neutrino'] = 'neu'
            to_save['neutrino_energy'] = 'Ev'
            to_save['target_material'] = 'tgt'

            for key, value in to_save.iteritems():
                self.log.info('%s : %s' % (key, str(t.__getattr__(value))))
                event_type[key] = t.__getattr__(value)

            self.log.debug('Event type:')
            for my_type in ['qel', 'res', 'dis', 'coh', 'dfr',
                            'imd', 'nuel', 'em']:
                if t.__getattr__(my_type) == 1:
                    self.log.debug('\t%s', my_type)
                event_type[my_type] = t.__getattr__(my_type)

            self.log.debug('Propogator:')
            for prop in ['nc', 'cc']:
                if t.__getattr__(prop) == 1:
                    self.log.debug('\t%s', prop)
                event_type[prop] = t.__getattr__(prop)

            yield next_events, event_type

        f.Close()
        os.remove(self.filenames[material])
651463
def evalHeatKernel(dgm1, dgm2, sigma):
    
    kSigma = 0
    I1 = np.array(dgm1)
    I2 = np.array(dgm2)
    for i in range(I1.shape[0]):
        p = I1[i, 0:2]
        for j in range(I2.shape[0]):
            q = I2[j, 0:2]
            qc = I2[j, 1::-1]
            kSigma += np.exp(-(np.sum((p - q) ** 2)) / (8 * sigma)) - np.exp(
                -(np.sum((p - qc) ** 2)) / (8 * sigma)
            )
    return kSigma / (8 * np.pi * sigma)
30963
def fermi_dist(energy, beta):
    
    exponent = np.asarray(beta*energy).clip(-600, 600)
    return 1./(np.exp(exponent) + 1)
355410
def sum_over_energy(self):
        
        # We sum over axis 0 in the array, and drop the energy binning in the
        # hpx object
        return HpxMap(np.sum(self.counts, axis=0), self.hpx.copy_and_drop_energy())
182159
def pwm_array2pssm_array(arr, background_probs=DEFAULT_BASE_BACKGROUND):
    
    b = background_probs2array(background_probs)
    b = b.reshape([1, 4, 1])
    return np.log(arr / b).astype(arr.dtype)
1112846
def gp_sims_panel(version):
  
  inDir, outDir = getWorkDirs()
  inDir = os.path.join(inDir, version)
  mesons = ['pion', 'eta', 'etap', 'rho', 'omega', 'phi', 'jpsi']
  fstems = ['cocktail_contribs/'+m for m in mesons] + ['cocktail', 'cocktail_contribs/ccbar']
  data = OrderedDict((energy, [
      np.loadtxt(open(os.path.join(inDir, fstem+str(energy)+'.dat'), 'rb'))
      for fstem in fstems
  ]) for energy in energies)
  for v in data.values():
      # keep syserrs for total cocktail
      v[-2][:,(2,3)] = 0
      # all errors zero for cocktail contribs
      v[-1][:,2:] = 0
      for d in v[:-2]: d[:,2:] = 0
  make_panel(
      dpt_dict = OrderedDict((
          ' '.join([getEnergy4Key(str(energy)), 'GeV']),
          [ data[energy], [
              'with %s lc %s lw 5 lt %d' % (
                  'lines' if i!=len(fstems)-2 else 'filledcurves pt 0',
                  default_colors[(-2*i-2) if i!=len(fstems)-2 else 0] \
                  if i!=len(fstems)-1 else default_colors[1], int(i==3)+1
              ) for i in xrange(len(fstems))
          ], [particleLabel4Key(m) for m in mesons] + [
              'Cocktail (w/o {/Symbol \162})', particleLabel4Key('ccbar')
          ]]
      ) for energy in energies),
      name = os.path.join(outDir, 'sims_panel'),
      ylog = True, xr = [0.,3.2], yr = [1e-6,9],
      xlabel = xlabel, ylabel = ylabel,
      layout = '2x2', size = '7in,9in',
      key = ['width -4', 'spacing 1.5', 'nobox', 'at graph 0.9,0.95']
  )
1178489
private void PowerSpectrum() {
        Power = new double[data.length][data[0].length];
        PowerMax = 0;

        for (int i = 0; i < data.length; i++) {
            for (int j = 0; j < data[0].length; j++) {
                double p = data[i][j];
                if (p < 0) p = -p;
                Power[i][j] = p;
                if (p > PowerMax) PowerMax = p;
            }
        }
    }
970980
def read_gbasis(basis_lines, fname):
    '''
    '''

    skipchars = '!#'
    basis_lines = [l for l in basis_lines if l and not l[0] in skipchars]

    bs_data = create_skel('component')

    i = 0
    bs_name = None
    while i < len(basis_lines):
        line = basis_lines[i]
        lsplt = line.split(':')
        elementsym = lsplt[0]

        if bs_name is None:
            bs_name = lsplt[1]
        elif lsplt[1] != bs_name:
            raise RuntimeError("Multiple basis sets in a file")

        element_Z = lut.element_Z_from_sym(elementsym)
        element_Z = str(element_Z)

        if not element_Z in bs_data['elements']:
            bs_data['elements'][element_Z] = {}

        element_data = bs_data['elements'][element_Z]

        if not 'electron_shells' in element_data:
            element_data['electron_shells'] = []

        i += 1

        max_am = int(basis_lines[i].strip())
        i += 1

        for am in range(0, max_am + 1):
            lsplt = basis_lines[i].split()
            shell_am = lut.amchar_to_int(lsplt[0])
            nprim = int(lsplt[1])
            ngen = int(lsplt[2])

            if shell_am[0] != am:
                raise RuntimeError("AM out of order in gbasis?")

            if max(shell_am) <= 1:
                func_type = 'gto'
            else:
                func_type = 'gto_spherical'


            shell = {
                'function_type': func_type,
                'region': '',
                'angular_momentum': shell_am
            }

            exponents = []
            coefficients = []

            i += 1
            for j in range(nprim):
                line = basis_lines[i].replace('D', 'E')
                line = line.replace('d', 'E')
                lsplt = line.split()

                if len(lsplt) != (ngen + 1):
                    raise RuntimeError("Incorrect number of general contractions in gbasis")

                exponents.append(lsplt[0])
                coefficients.append(lsplt[1:])
                i += 1

            shell['exponents'] = exponents

            # We need to transpose the coefficient matrix
            # (we store a matrix with primitives being the column index and
            # general contraction being the row index)
            shell['coefficients'] = list(map(list, zip(*coefficients)))

            element_data['electron_shells'].append(shell)

    return bs_data
1299508
def yield_pair_energies(self, index1, index2):
        """"""
        A = self.As[index1, index2]
        B = self.Bs[index1, index2]
        distance = self.distances[index1, index2]
        yield A*np.exp(-B*distance), 1
667628
def gp_panel(version, skip):
  
  scale = { # QM14 (19 GeV skip later, factor here only informational)
    '19': 1.0340571932983775, '200': 1.0, '39': 0.7776679085382481,
    '27': 0.6412140408244136, '62': 0.9174700031778402
  } if version == 'QM14' else {
      '19': 1.3410566491548412, '200': 1.1051002240771077,
      '39': 1.2719203877292842, '27': 1.350873678084769,
      '62': 1.2664666321635087
  }
  inDir, outDir = getWorkDirs()
  inDir = os.path.join(inDir, version)
  data = {}
  vacRhoTitle = '{/Symbol \162}/{/Symbol \167} VacSF+FB'
  for infile in os.listdir(inDir):
    if infile == "cocktail_contribs": continue
    if infile == 'mediumDmOnly200.dat': continue
    energy = re.compile('\d+').search(infile).group()
    if skip is not None and energy == skip: continue
    data_type = re.sub('%s\.dat' % energy, '', infile)
    file_url = os.path.join(inDir, infile)
    data_import = np.loadtxt(open(file_url, 'rb'))
    if data_type != 'data' and (
        (version == 'QM14' and energy != '19') or version == 'LatestPatrickJieYi'
    ):
       data_import[:,(1,3,4)] /= scale[energy]
    if data_type == 'cocktail': data_import[:,2:] = 0.
    elif fnmatch(data_type, '*medium*') or data_type == 'vacRho':
       data_import = data_import[data_import[:,0] < 0.9] \
               if energy == '200' and data_type == '+medium' \
               else data_import
       data_import[:,(2,3)] = 0.
    key = getEnergy4Key(energy)
    if key not in data: data[key] = {}
    data_type_mod = data_type
    if data_type_mod == 'cocktail': data_type_mod = 'Cocktail'
    elif data_type_mod == 'mediumMedOnly': data_type_mod = 'HMBT'
    elif data_type_mod == 'mediumQgpOnly': data_type_mod = 'QGP'
    elif data_type_mod == '+medium': data_type_mod = 'Cock. + HMBT + QGP'
    elif data_type_mod == 'vacRho': data_type_mod = vacRhoTitle
    data[key][data_type_mod] = data_import
  plot_order = ['Cocktail', vacRhoTitle, 'QGP', 'HMBT', 'Cock. + HMBT + QGP', 'data']
  plot_opts = {
    vacRhoTitle: 'with lines lt 2 lw 5 lc %s' % default_colors[6],
    'QGP': 'with lines lt 2 lw 5 lc %s' % default_colors[1],
    'HMBT': 'with lines lt 2 lw 5 lc %s' % default_colors[2],
    'Cock. + HMBT + QGP': 'with filledcurves lt 1 lw 5 pt 0 lc %s' % default_colors[16],
    'Cocktail': 'with lines lc %s lw 5 lt 1' % default_colors[8],
    'data': 'lt 1 lw 4 ps 1.5 lc %s pt 18' % default_colors[0]
  }
  panel2D_versions = (version == 'LatestPatrickJieYi' or version == 'QM14')
  make_panel(
    dpt_dict = OrderedDict(
      (' '.join([k, 'GeV %s' % (
          '{/=18 PRL 113 022301}'
          if k == '200' and (
              version == 'QM12Latest200' or version == 'QM14' or
              version == 'LatestPatrickJieYi'
          ) else '{/=18 STAR Preliminary}'
      )]), [
        [ data[k][dt] for dt in plot_order if dt in data[k] ],
        [ plot_opts[dt] for dt in plot_order if dt in data[k] ],
        [ dt for dt in plot_order if dt in data[k] ]
      ]) for k in sorted(data, key=float)
    ),
    name = os.path.join(
      outDir, 'panel%s%s' % (version, 'No'+skip if skip is not None else '')
    ),
    ylabel = '1/N@_{mb}^{evt} dN@_{ee}^{acc.}/dM_{ee} [ (GeV/c^2)^{-1} ]',
    xlabel = 'invariant dielectron mass, M_{ee} (GeV/c^{2})',
    ylog = True, xr = [0.05, 1.1], yr = [1e-4, 0.5],
    #lmargin = 0.12 if panel2D_versions else 0.1,
    #bmargin = 0.11 if panel2D_versions else 0.15,
    arrow_length = 0.4, arrow_bar = 0.002,
    gpcalls = [
        'mxtics 2',
        'object 50 rectangle back fc rgb "#C6E2FF" from 0.4,1e-4 to 0.75,2e-2'
    ],
    layout = '3x2' if panel2D_versions else ('%dx1' % len(data)),
    key = ['width -2', 'at graph 0.95,0.85'],
    key_subplot_id = 5 if version != 'QM12' else 0,
    size = '8in,8in' if version != 'QM12' else '5in,8.5in'
  )
  return 'done'