Sorting function for non commuting bosons

I am trying to write a sorting function which will sort expressions involving products of bosonic objects which do not commute. For example, I can have objects like $a,\ a^\dagger,\ b,\ b^\dagger$ where $[a,a^\dagger] = 1$, $[b,b^\dagger] = 1$ and all other commutators vanish. So the sorting function should for example sort a term like $a a^\dagger b a b^\dagger$ to $a + a^\dagger a^2 + ab^\dagger b+a^\dagger a^2 b^\dagger b$. So I want to sort them first lexicographically and then such that dagger objects will be to the left.

I have been using a function from the notebook provided by Simon Tyler at arXiv. So I have created a bosonic object using

Clear[Boson, BosonC, BosonA]
Boson /: MakeBoxes[
 Boson[cr : (True | False), sym_], fmt_] := 
  With[{sbox = If[StringQ[sym], sym, ToBoxes[sym]]}, 
   With[{abox = If[cr, SuperscriptBox[#, FormBox["\[Dagger]", Bold]], #] &@sbox},

    InterpretationBox[abox, Boson[cr, sym]]
   ]
  ]
BosonA[sym_:"a"] := Boson[False, sym]
BosonC[sym_:"a"] := Boson[True, sym]

Next following the code, I want to turn the atomic expressions into lists and use Mathematica's ordering after that. So I have functions like

NCOrder[Boson[cr : True | False, sym_]] := 
NCOrder[Boson[cr, sym]] = {{sym, Boole[cr]}}


NCSort[expr_] := Module[{temp, NCM},
 temp = expr /. NonCommutativeMultiply -> NCM;   
 temp = temp /. (NCM[a : _Boson ..] :> (NCM[a][[Ordering[#]]] &[NCOrder /@ {a}]));
 temp = temp /. NCM[singleArg_] :> singleArg;   
 temp = temp /. NCM -> NonCommutativeMultiply
]

This NCSort function right now is just ordering it lexicographically. I want to modify it so that it would sort and order the expressions as I explained above.

Answer

I suggest the following function, which looks rather ugly and is not really efficient, but seems to do the job:

expand[expr_] :=
 Module[{times},
   expr /.
     NonCommutativeMultiply[b___] :> times[b] //.
     {
       times[left___, cnum_ /; FreeQ[cnum, _Boson], right___] :> 
             cnum*times[left, right],
       times[left___, fst : Boson[_, s_], middle__, sec : Boson[_, s_],right___] /;
          FreeQ[{middle}, Boson[_, s]] &&
              OrderedQ[{Last@Union[Cases[{middle}, Boson[_, sym_] :> sym]], s}] :>

                  times[left, middle, fst, sec, right],
       times[left___, Boson[False, s_], Boson[True, s_], right___] :>
          times[left, right] + times[left, Boson[True, s], Boson[False, s], right],
       times[b_Boson] :> b
     } //.
    times[left___, p : Longest[PatternSequence[(b : Boson[type_, s_]) ..]],right___] :>
      times[left, b^Length[{p}], right] /. 
    times[arg__] :> NonCommutativeMultiply[arg] /. times[] -> OverHat[1]
];

You can use it as

expand[NonCommutativeMultiply[BosonA["a"], BosonC["a"], BosonA["b"], BosonA["a"], BosonC["b"]]]

(*
    a+a^\[Dagger]**a^2+a**b^\[Dagger]**b+a^\[Dagger]**a^2**b^\[Dagger]**b
*)