I think the easiest way to do this is to simply implement all rules of the operator arithmetic explicitly, rather than trying to be too clever about it. For example, the + operator in vec is not in fact the standard $+:mathbb{C}timesmathbb{C}rightarrowmathbb{C}$, but rather $+:mathbb{(Crightarrow C)}timesmathbb{(Crightarrow C)}rightarrowmathbb{(Crightarrow C)}$. So in some sense, the rules anyway need to be redefined on a case-by-case basis.

Here's how I would do it: First, choose a head operator that designates that an expression should follow our operator arithmetic. We define what happens if this operator is applied to a list of arguments by implementing special cases for the common operations (using the tricks from the question):

operator[op : (List | Plus | Times)[___]][args___] := Through[(operator /@ op)[args]]
operator[op : (Composition | RightComposition)[___]][args___] := (operator /@ op)[args]

We then define a fallback rule for constants:

operator[const_][_] := const

Finally, we define our custom operators:

operator[f1][x_] := c1 x
operator[f2][x_] := c2 x
operator[f3][x_] := c3 x
operator[f4][x_] := c4 x

Now, the composite operator vec:

vec = operator[{f1@*f3 - f2@*f4, f3 + 2 f4 + 4, const1 + Exp[x] - f1}]
(* operator[{f1@*f3 - f2@*f4, 4 + f3 + 2 f4, const1 + E^x - f1}] *)

Applying it to a:

vec[a]
(* {a c1 c3 - a c2 c4, 4 + a c3 + 2 a c4, -a c1 + const1 + E^x} *)

You'll note that what we're doing here is essentially the multilevel/recursive mapping you were suggesting in your question, just with a helper head for readability and ease of use. Also, you'll probably need to add some more rules for arithmetic operations, e.g. Power.

Further remarks

One potential pitfall that comes to mind is your distinction of multiplication and composition: Usually in QM, operators are always linear, and so a*b is used to unambiguously refer to the composition of operators a and b (since (a*b)[c] must be a[b[c]] and not a[c]*b[c], as the latter would not be linear). If you choose to adopt that notation, constants would instead be operators multiplying their argument with themselves, so e.g.: 2[c] is 2*c. This notation has the advantage that a^3 and by extension Exp[a] can be directly interpreted as a*a*a and 1+a+a*a/2+…, respectively. Finally, (1+a)[c] is then c+a[c] which is again linear, rather than 1+a[c], which is not.

An implementation of this would look something like this:

Clear@operator

operator[op : (List | Plus)[___]][args___] := Through[(operator /@ op)[args]]
operator[op : (Composition | Times | NonCommutativeMultiply)[___]][args___] := (operator /@ Composition @@ op)[args]
operator[const_][x_] := const x

operator[f1][x_] := c1 x
operator[f2][x_] := c2 x
operator[f3][x_] := c3 x
operator[f4][x_] := c4 x

vec = operator[{f1 ** f3 - f2 ** f4, f3 + 2 f4 + 4, const1 + Exp[x] - f1}]
(* operator[{f1 ** f3 - f2 ** f4, 4 + f3 + 2 f4, const1 + E^x - f1}] *)

vec[a]
(* {a c1 c3 - a c2 c4, 4 a + a c3 + 2 a c4, -a c1 + a const1 + a E^x} *)

The end result is of course slightly different in this case, but as noted, this way, the vec operator is linear, whereas in your interpretation it is not (which might or might not be desired of course). Also note the use of ** (NonCommuativeMultiply) instead of *, since of course operator multiplication/composition is in general not commutative. Nevertheless, I implemented the second rule in such a way that it supports *, **, and @* for increased flexibility, while interpreting all of them as @*.

Remarks on edit

Here are some remarks on the edit to the question, and the alternate approach proposed there:

• First of all, the second rule can be made quite a bit more readable by replacing op1_ /; MatchQ[op1, opList] by the more direct op1: opList (for the documentation, see Pattern):

  opApply[head_[args__]] := 
   ReplaceAll[Expand[head], {(op1:opList)@*(op2:opList) :> op1[op2[args]], op1:opList :> op1[args]}]
  

• In general, I prefer the approach shown in this answer, because it is (at least in my experience) often more readable and leverages Mathematica's strength in pattern matching better by separating the different rules into separate downvalues. However, sometimes an approach based on ReplaceAll is also the better solution, especially in situations involving tight evaluation control. And of course, it's also a matter of personal preference.

• Furthermore, I think the approach is not so different in it's core, with the two main differences being:

  • The use of ReplaceAll instead of individual rules for all expressions: This enables easier traversal of deeper expressions, since we don't need to implement rules for all levels. Of course, this gives us less control, but it's also less work if we don't need the control. We could emulate this for the operator approach as follows (untested):

    operator[h_[parts___]][args__] := (operator[h][args]) @@ (operator[#][args] & /@ {parts})
    operator[atom_][__] := atom
    

    The first rule essentially makes sure operator is "applied" to every part of an expression, while the second rule effectively terminates the traversal of the expression once the bottom is reached. All other rules would then correspond to the actual replacement rules

  • The use of a function that is applied to a full expression as opposed to a helper head that is wrapped around the operator part. We can also emulate this for the operator approach:

    opApply[op_[args__]]:=operator[op][args]
    

    Clearly this is cheating in some sense, but it does somehow combine the ease of use of the opApply approach with the ease of definition of the operator approach.