When converting a LiDAR dataset to a DEM, you are taking a set of discrete data points and converting them into a single, continuous dataset. Let's say that your .las file contains X (latitude), Y (longitude) and Z (elevation) values with an average resolution of ~1 meters. The resolution here is really important- we're only talking about an average and so we'll hardly find that ~1 meter resolution throughout the dataset. Instead we'll find values that fall within a 'ballpark' estimate of that resolution. So you take these points and convert them to a raster DEM, or maybe a TIN. X and Y values should see minimum distortion, but you'll notice that your Z values might not be what you expected. This is because the computer doesn't know what the correct Z values are in cells that do not fall on one of your LiDAR points. In between LiDAR points, an interpolation algorithm has been applied to estimate what reasonable Z values might be. Choosing the correct interpolation method relative to the goals of your analysis is a crucial part of going from LiDAR to DEM. Setting the proper resolution on that output DEM is important- always set a lower resolution than the resolution of your LiDAR dataset. So for a ~1 meter resolution I would set a 3 meter resolution for the DEM, in an effort to minimize distortion.
I have experience studying landslides and debris flows with DEMs derived from LiDAR. Landslides and debris flows are very linear features which occur near other linear features in the topography. So when I convert from LiDAR to DEM, I want an interpolation method that best accentuates linear features. This happens to be a TIN (Triangulated Irregular Network). You say that you aim to do a hydrological analysis. Maybe you should try a spline interpolation method to build your DEM. Spline interpolation works be drawing continuous, overlapping lines through all of your data points in order to create a very smooth raster surface. Identify your sinks, fill them, draw countours, repeat.
This is a bit of a ramble, but what I'm trying to get at here is that it seems to me like you're asking the wrong question. Instead of asking for a software workflow you should use to build a hydrologically correct DEM, you should be asking what interpolation method to use. If I were you, I'd try a spline interpolation method.
In terms of software, processing LiDAR data is waaay CPU/RAM intensive. If you have >6 GB of RAM I would recommend GRASS GIS. They've got the best LiDAR processing software I've ever used (it's FOSS), but you have to do some memory allocation. Otherwise, I would recommend sticking with ArcGIS. They've got some great documentation on how to do what you want to do on their website.
I just finished writing a script to accomplish this task using the free and open-source GIS Whitebox Geospatial Analysis Tools (download here), for which I am the lead developer. The source code of the script can be found here. Although the script is not yet part of the current official Whitebox release (v. 3.2.1) you can get an early version of it by selecting Update Scripts From Repository from the Tools menu. You will need to relaunch Whitebox after doing this. You should then find a new script in the LiDAR toolbox called Filter LAS Scan Angles. Unfortunately you won't have an associated Help file until the public release of the tool with the next version of Whitebox.
The tool takes any number of LAS file inputs and will filter them in parallel. Thus, the performance should scale with the number of cores you have on your system. If you want to integrate the tool within a workflow involving multiple steps or run it in batch mode, you can also call the tool from a Python script in the Whitebox Scripter using the following as an example:
wd = pluginHost.getWorkingDirectory()
# You may have multiple input files but they must # be separated by semicolons in the string.
inputs = wd + "input1.las" + ";" + wd + "input2.las" + ";" + wd + "input3.las"
threshold = "15.0"
args = [inputs, threshold]
pluginHost.runPlugin("FilterLasScanAngles", args, False)
The tool will output a ShapeFile of MULTIPOINTZ shapetype for each LAS input. The elevation and intensity values are stored in the Z and M properties of the points respectively. I hope that suffices, if not, I can work on it to output LAS files, but it'll take more time.
Here's an example, where I've filtered a LAS file to remove all points with a scan angle greater than 4 degrees. The resulting ShapeFile (in green palette) has been overlayed overtop the original LAS file (in blue palette) and you can see the flight line pattern associated with varying point scan angle.
Best Answer
FUSION LiDAR Toolkit has clipping capabilities (PolyClipData tool) and unlike LASTools, its usage is unrestricted. However, despite the fact that some SVN repository on SourceForge exists, the source code published is incomplete and very old. If you can proceed without knowing the code and just use the compiled binary, then FUSION should be fine for this task.