I don’t know much on processes involving Z decays (I am doing simulations with fast neutrons), but you should see in your vector of secondaries a mu+ and mu- . For these, you can then get kinetic energy, momentum, creation process etc so that you should be able to calculate the invariant mass.
Let’s say you call your vector of secondaries “secondary”, then you can do something like:
const std::vector<const G4Track*>* secondary = step->GetSecondaryInCurrentStep();
G4String name = (*secondary)[i]->GetDefinition()->GetParticleName();
G4String process = (*secondary)[i]-> GetCreatorProcess()-> GetProcessName();
etc. Look at the G4Track class reference to see more quantities that you can retrieve.
I hope this helps. Maybe there are other possibilities that I don’t know that are more elegant.
I recommend again the Hadr03 example.