code: def _density_par(self,dangle,tdisrupt=None,approx=True,
                       higherorder=None):
          
          if higherorder is None: higherorder= self._higherorderTrack
          if tdisrupt is None: tdisrupt= self._tdisrupt
          if approx:
              return self._density_par_approx(dangle,tdisrupt,
                                              higherorder=higherorder)
          else:
              return integrate.quad(lambda T: numpy.sqrt(self._sortedSigOEig[2])\
                                        *(1+T*T)/(1-T*T)**2.\
                                        *self.pOparapar(T/(1-T*T)\
                                                            *numpy.sqrt(self._sortedSigOEig[2])\
                                                            +self._meandO,dangle),
                                    -1.,1.)[0]
docs: The raw density as a function of parallel angle,
          approx= use faster method that directly integrates the spline
          representation