Taipei Tutorial 2013, problem with the decay chain

Hi I've been working on your tutorials, which I'll put at the end as reference, my problem is its working perfectly when I only run the process as p,p-> uv uv~ however when I add the decay chain to draw the plot of missing transverse energy by phiM1 it draws nothing. I checked all the decay processes one by one but it seems uv's are not decaying. I also tried to adjust the kinematics however it didn't effect. It basically draws a line at zero but measures perfectly when I just measure properties of uv. Since I dont know which part is wrong Im adding both model and batch files.

references: 1. https://feynrules.irmp.ucl.ac.be/attachment/wiki/Talks/Tutorial.pdf
2. arXiv: 1209.0297v1

FeynRules Model File:

IndexRange[Index[scInd]] = Unfold[Range[2]];
IndexStyle[scInd, s];

M$Parameters = {
      lambda == {
      ParameterType -> External,
      ComplexParameter -> False,
      Indices -> {Index[scInd]},
      Value -> {lambda[1]->1, lambda[2]->1},
      InteractionOrder -> {NP,1}
      },

      lambdap == {
      ParameterType -> External,
      ComplexParameter -> False,
      Indices -> {Index[scInd]},
      Value -> {lambdap[1]->1, lambdap[2]->1},
      InteractionOrder -> {NP,1}
      },

      sina == {
         ParameterType -> External,
    ComplexParameter -> False,
    Value -> 0.35
    },
      cosa == {
         ParameterType -> Internal,
    ComplexParameter -> False,
    Value -> Sqrt[1-sina^2]
    },

      MassM == {
      ParameterType -> Internal,
      ComplexParameter -> False,
      Indices -> {Index[scInd], Index[scInd]},
      Definitions -> {MassM[1, 1] -> cosa^2*MphiM1^2 + MphiM2^2*sina^2, MassM[1, 2] -> -(cosa*(MphiM1^2 - MphiM2^2)*sina), MassM[2, 1] -> -(cosa*(MphiM1^2 - MphiM2^2)*sina), MassM[2, 2] -> cosa^2*MphiM2^2 + MphiM1^2*sina^2}
      }
};

M$ClassesDescription = {
       F[20] == {
           ClassName -> uv,
      SelfConjugate -> False,
      Indices -> {Index[Colour]},
      ParticleName -> "uv",
      AntiParticleName -> "uv~",
      QuantumNumbers-> {Q->2/3, Y->2/3},
      Mass -> {Muv, 500},
      Width -> {Wuv, 1}
      },

        F[21] == {
            ClassName -> ev,
      SelfConjugate -> False,
      Indices -> {},
      ParticleName -> "ev",
      AntiParticleName -> "ev~",
      QuantumNumbers-> {Q->-1, Y->-1, LeptonNumber->1},
      Mass -> {Mev, 300},
      Width -> {Wev,1}
         },

   S[21] == {
         ClassName -> phiM,
         SelfConjugate -> True,
         Indices -> {Index[scInd]},
         FlavorIndex -> scInd,

         ClassMembers -> {phiM1, phiM2},
         Mass -> {MphiM, {MphiM1, 200}, {MphiM2,400}},
         Width -> {WphiM, {WphiM1,1}, {WphiM2,1}}
         },

         S[20] == {
               ClassName -> phiT,
        Unphysical -> True,
        SelfConjugate -> True,
        Indices -> {Index[scInd]},

        FlavorIndex -> scInd,
        Definitions -> {phiT[1]-> cosa phiM[1] + sina phiM[2], phiT[2]-> -sina phiM[1] + cosa phiM[2]}
         }

};

Lsq = 1/2 del[phiT[ii],mu] del[phiT[ii],mu] - 1/2 phiT[ii] MassM[ii,jj] phiT[jj];

Lfq = I uvbar.Ga[mu].DC[uv,mu] - Muv uvbar.uv + I evbar.Ga[mu].DC[ev,mu] - Mev evbar.ev;

Lyuk = lambda[ii] phiT[ii] uvbar.ProjP.u + lambdap[ii] phiT[ii] evbar.ProjP.e;

Lnew = Lsq + Lfq + Lyuk + HC[Lyuk];

CalcHEP batch file (last version modified as given in ref 2)

Process: p,p->uv,uv~
Decay: uv -> u, phiM1
Decay: uv~ -> u~, phiM1
Decay: uv-> u,phiM2
Decay: uv~->u~,phiM2
Decay: phiM2 -> e+, ev
Decay: phiM2 -> e-, ev~
Decay: ev -> e-, phiM1
Decay: ev~ -> e+, phiM1

Composite: p=u,u~,d,d~,s,s~,c,c~,b,b~,g
Composite: jet=u,u~,d,d~,s,s~,c,c~,b,b~,g

pdf1: cteq6l (proton)
pdf2: cteq6l (proton)

p1: 4000
p2: 4000

Dist parameter: Z(phiM1)
Dist min: 0
Dist max: 1000
Dist n bins: 300
Dist title: p,p->uv, uv~
Dist x-title: Missing Trans. Eng. (GeV)

Number of events (per run step): 1000
Filename: MC4BSM
NTuple: False
Cleanup: False

Max number of cpus: 2
sleep time: 1
nSess_1: 5
nCalls_1: 100000
nSess_2: 5
nCalls_2: 100000

Thanks a lot

Jack