def get_analysis_and_structure(self,
                                   structure,
                                   calculate_valences=True,
                                   guesstimate_spin=False,
                                   op_threshold=0.1):
        """
        Obtain an analysis of a given structure and if it may be Jahn-Teller
        active or not. This is a heuristic, and may give false positives and
        false negatives (false positives are preferred).

        :param structure: input structure
        :param calculate_valences (bool): whether to attempt to calculate valences or not, structure
        should have oxidation states to perform analysis
        :param guesstimate_spin (bool): whether to guesstimate spin state from magnetic moments
        or not, use with caution
        :param op_threshold (float): threshold for order parameter above which to consider site
        to match an octahedral or tetrahedral motif, since Jahn-Teller structures can often be
        quite distorted, this threshold is smaller than one might expect
        :return (dict): analysis of structure, with key 'strength' which may be 'none', 'strong',
        'weak', or 'unknown'
        """

        structure = structure.get_primitive_structure()

        if calculate_valences:
            bva = BVAnalyzer()
            structure = bva.get_oxi_state_decorated_structure(structure)

        # no point testing multiple equivalent sites, doesn't make any difference to analysis
        # but makes returned
        symmetrized_structure = SpacegroupAnalyzer(structure).get_symmetrized_structure()

        # to detect structural motifs of a given site
        op = LocalStructOrderParams(['oct', 'tet'])

        # dict of site index to the Jahn-Teller analysis of that site
        jt_sites = []
        non_jt_sites = []

        for indices in symmetrized_structure.equivalent_indices:

            idx = indices[0]
            site = symmetrized_structure[idx]

            # only interested in sites with oxidation states
            if isinstance(site.specie, Specie) and site.specie.element.is_transition_metal:

                # get motif around site
                order_params = op.get_order_parameters(symmetrized_structure, idx)

                if order_params[0] > order_params[1] and order_params[0] > op_threshold:
                    motif = 'oct'
                    motif_order_parameter = order_params[0]
                elif order_params[1] > op_threshold:
                    motif = 'tet'
                    motif_order_parameter = order_params[1]
                else:
                    motif = 'unknown'
                    motif_order_parameter = None

                if motif == "oct" or motif == "tet":

                    # guess spin of metal ion
                    if guesstimate_spin and 'magmom' in site.properties:

                        # estimate if high spin or low spin
                        magmom = site.properties['magmom']
                        spin_state = self._estimate_spin_state(site.specie, motif, magmom)
                    else:
                        spin_state = "unknown"

                    magnitude = self.get_magnitude_of_effect_from_species(site.specie,
                                                                          spin_state,
                                                                          motif)

                    if magnitude != "none":

                        ligands = get_neighbors_of_site_with_index(structure, idx,
                                                                   approach="min_dist",
                                                                   delta=0.15)
                        ligand_bond_lengths = [ligand.distance(structure[idx])
                                               for ligand in ligands]
                        ligands_species = list(set([str(ligand.specie) for ligand in ligands]))
                        ligand_bond_length_spread = max(ligand_bond_lengths) - \
                                                    min(ligand_bond_lengths)

                        def trim(f):
                            # avoid storing to unreasonable precision, hurts readability
                            return float("{:.4f}".format(f))

                        # to be Jahn-Teller active, all ligands have to be the same
                        if len(ligands_species) == 1:
                            jt_sites.append({'strength': magnitude,
                                             'motif': motif,
                                             'motif_order_parameter': trim(motif_order_parameter),
                                             'spin_state': spin_state,
                                             'species': str(site.specie),
                                             'ligand': ligands_species[0],
                                             'ligand_bond_lengths': [trim(length) for length in
                                                                     ligand_bond_lengths],
                                             'ligand_bond_length_spread':
                                                 trim(ligand_bond_length_spread),
                                             'site_indices': indices})

                    # store reasons for not being J-T active
                    else:
                        non_jt_sites.append({'site_indices': indices,
                                             'strength': "none",
                                             'reason': "Not Jahn-Teller active for this "
                                                       "electronic configuration."})
                else:
                    non_jt_sites.append({'site_indices': indices,
                                         'strength': "none",
                                         'reason': "motif is {}".format(motif)})

        # perform aggregation of all sites
        if jt_sites:
            analysis = {'active': True}
            # if any site could exhibit 'strong' Jahn-Teller effect
            # then mark whole structure as strong
            strong_magnitudes = [site['strength'] == "strong" for site in jt_sites]
            if any(strong_magnitudes):
                analysis['strength'] = "strong"
            else:
                analysis['strength'] = "weak"
            analysis['sites'] = jt_sites
            return analysis, structure
        else:
            return {'active': False, 'sites': non_jt_sites}, structure