It seems that the modeling of the step input for the current is rather not physical. The invariant is delayed perturbed in a non-neglectable manner.

First, take into account that NDSolve is proven in practice and not a built-in that is made for step size control. That was a wish by the Mathematica customers and a consequence of the test suits on which it is verified for the highest possible versatility. It saves time and effort for the users.

There are still of course workarounds to control [step size][1] for tedious and error-prone partial differential equations. This case is not close to such problems.

I did put the precision options out of the command. Automatic means for the Projection option too the use of [`Runge-Kutta`][2]. A proper choice for it is to match the order of the coefficients favored. So with `Automatic` these whole lot of optimization are done by Mathematica in NDSolve for the best solution.

With this in mind the results for 

 Remove["Global`*"]
 
 (*PARAMETERS*)
 \[Gamma] = 
 1.76*10.^7*10.^4 (*T^-1s^-1*); \[Alpha] = 0.06; \[Theta]sh = 0.3; Ms \
 = 7.*10.^5;(*A/m*)tFM = 2.*10.^-9;(*m*)currAMP = 
 30.*10.^-3 (*A*); sHM = 10.*10.^-6*3.*10.^-9(*m^2,tHM=3 nm*);
 Hamp = 0.03; Hk = 0.1;
 a = 0. \[Degree]; b = 0. \[Degree]; c = 0. \[Degree]; Ang = 
 0. \[Degree];
 temp = 3; rad = 20.*10.^-9;(*Thermal effect*)kTL = 2.;
 tSt = 0.; tFin = 1.5*10.0^-8; tStep = 1.*10.0^-12;
 
 (*CURRENT PULSE*)
 tPulse = 5.*10^-9; tRise = 0.5*10^-9;
 jcurr[t_] := Piecewise[{{1., -1 <= t <= tPulse}}, 0.];
 Cjconst = (3.29375`*^-16) (\[Theta]sh/Ms/tFM)*currAMP/sHM;
 Cj[t_] := Cj[t] = jcurr[t]*Cjconst;
 
 (*FIELD COMPONENTS*)
 \[Sigma]f[temp_, rad_] := 
 Sqrt[(2.*\[Alpha]*1.38*10.^-23*temp)/(\[Gamma]*Ms*tFM*3.1415*rad^2*
 tStep)]; \[Sigma] = \[Sigma]f[temp, rad];(*stand.dev*)
 RandComp[] := 
 RandomVariate[NormalDistribution[0., \[Sigma]]]/
 Sqrt[3.];(*random component*)Happ = ({{Cos[c], -Sin[c], 
 0.}, {Sin[c], Cos[c], 0.}, {0., 0., 1.}}).({{Cos[b], 0., 
 Sin[b]}, {0., 1., 0.}, {-Sin[b], 0., Cos[b]}}).({{1., 0., 
 0.}, {0., Cos[a], -Sin[a]}, {0., Sin[a], 
 Cos[a]}}).({{Hamp*
 Cos[Ang]}, {Hamp*(-Sin[
 Ang])}, {0.}});(*permanent applied component*)
 Hrand[t_] := Hrand[t] = {RandComp[], RandComp[], RandComp[]};
 H[t_, mz_] := Flatten[Happ + Hrand[t] + {0., -kTL*Cj[t], Hk*mz}]
 
 (*ODE SYSTEM*)
 t0 = 0;
 sol = NDSolve[({{mX'[
 t] == \[Gamma] (-mY[t]*H[t, mZ[t]][[3]] + 
 mZ[t]*H[t, 
 mZ[t]][[2]] - \[Alpha] (mY[
 t] (mX[t]*H[t, mZ[t]][[2]] - mY[t]*H[t, mZ[t]][[1]]) - 
 mZ[t] (mZ[t]*H[t, mZ[t]][[1]] - mX[t]*H[t, mZ[t]][[3]])) +
 Cj[t]*mY[t]*mX[t] + \[Alpha]*Cj[t]*mZ[t]), 
 mX[0] == 
 0.}, {mY'[
 t] == \[Gamma] (-mZ[t]*H[t, mZ[t]][[1]] + 
 mX[t]*H[t, 
 mZ[t]][[3]] - \[Alpha] (mZ[
 t] (mY[t]*H[t, mZ[t]][[3]] - mZ[t]*H[t, mZ[t]][[2]]) - 
 mX[t] (mX[t]*H[t, mZ[t]][[2]] - mY[t]*H[t, mZ[t]][[1]])) -
 Cj[t] (mX[t]^2 + mZ[t]^2)), 
 mY[0] == 
 0.}, {mZ'[
 t] == \[Gamma] (-mX[t]*H[t, mZ[t]][[2]] + 
 mY[t]*H[t, 
 mZ[t]][[1]] - \[Alpha] (mX[
 t] (mZ[t]*H[t, mZ[t]][[1]] - mX[t]*H[t, mZ[t]][[3]]) - 
 mY[t] (mY[t]*H[t, mZ[t]][[3]] - mZ[t]*H[t, mZ[t]][[2]])) +
 Cj[t]*mY[t]*mZ[t] - \[Alpha]*Cj[t]*mX[t]), 
 mZ[0] == 1.}}), {mX, mY, mZ}, {t, tSt, tFin}, 
 StartingStepSize -> 1, 
 Method -> {"Projection", Method -> "Automatic", 
 "Invariants" -> (mX[t]^2 + mY[t]^2 + mZ[t]^2 == 1)}(*,
 StepMonitor\[RuleDelayed]Block[{},\[Delta]t=t-t0;t0=t;
 Print["Current time: ",t];
 Print["Current step: ",\[Delta]t];
 Print["mX: ",mX[t]];
 Print["mY: ",mY[t]];
 Print["mZ: ",mZ[t],"
 "];
 AbsoluteTiming[Pause[1]]]*)]

are

[![InterpolationFunction][3]][3]

[![Graphics][4]][4]

The peak at the delayed time t about `5.5 ns` for the invariant is a violation of the foundational physics for the problem. Since the problem is reduced from the mathematical standpoint this is the solution.

The best path for a more physical solution will be to model that the peak does not appear. For example model some counter current from the other side of the electric system or some immanent current sources that delay in the calculated or near calculated manner.

Another interpretation might be, that there are times of unphysical solutions unavoidable in this modeling. Another one is a sign is set wrong. Nevertheless, the system recovers. 

I change the method especially for the question to:

 t0 = 0;
 sol = NDSolve[{{mX'[
 t] == \[Gamma] (-mY[t]*H[t, mZ[t]][[3]] + 
 mZ[t]*H[t, 
 mZ[t]][[2]] - \[Alpha] (mY[
 t] (mX[t]*H[t, mZ[t]][[2]] - mY[t]*H[t, mZ[t]][[1]]) - 
 mZ[t] (mZ[t]*H[t, mZ[t]][[1]] - mX[t]*H[t, mZ[t]][[3]])) +
 Cj[t]*mY[t]*mX[t] + \[Alpha]*Cj[t]*mZ[t]), 
 mX[0] == 
 0.}, {mY'[
 t] == \[Gamma] (-mZ[t]*H[t, mZ[t]][[1]] + 
 mX[t]*H[t, 
 mZ[t]][[3]] - \[Alpha] (mZ[
 t] (mY[t]*H[t, mZ[t]][[3]] - mZ[t]*H[t, mZ[t]][[2]]) - 
 mX[t] (mX[t]*H[t, mZ[t]][[2]] - mY[t]*H[t, mZ[t]][[1]])) -
 Cj[t] (mX[t]^2 + mZ[t]^2)), 
 mY[0] == 
 0.}, {mZ'[
 t] == \[Gamma] (-mX[t]*H[t, mZ[t]][[2]] + 
 mY[t]*H[t, 
 mZ[t]][[1]] - \[Alpha] (mX[
 t] (mZ[t]*H[t, mZ[t]][[1]] - mX[t]*H[t, mZ[t]][[3]]) - 
 mY[t] (mY[t]*H[t, mZ[t]][[3]] - mZ[t]*H[t, mZ[t]][[2]])) +
 Cj[t]*mY[t]*mZ[t] - \[Alpha]*Cj[t]*mX[t]), 
 mZ[0] == 1.}}, {mX, mY, mZ}, {t, tSt, tFin}, 
 StartingStepSize -> (tFin - tSt)/1000, 
 Method -> {"FixedStep", Method -> "ExplicitRungeKutta"}];
 StepDataPlot[sol, PlotRange -> {tSt, tFin}]
 Plot[{mX[t], mY[t], mZ[t], mX[t]^2 + mY[t]^2 + mZ[t]^2} /. sol, {t, 
 tSt, tFin}]

and got:

[![StepDataPlot][5]][5]

[![Better results in the invariant][6]][6]

This still not 1 anymore after the current is set on. The peak is gone!

Another attempt at this question is rather different but still with Mathematica possible. Use the SystemModelExamples new in newer versions of Mathematica.


 [1]: https://reference.wolfram.com/language/tutorial/NDSolveFixedStep.html
 [2]: https://reference.wolfram.com/language/tutorial/NDSolveExplicitRungeKutta.html
 [3]: https://i.sstatic.net/5eQrf.png
 [4]: https://i.sstatic.net/PDhjG.png
 [5]: https://i.sstatic.net/JGEzk.png
 [6]: https://i.sstatic.net/XpMJr.png