SetDirectory["/Users/maxim/Cargese_Tutorial/"];
Timing[testp=Import["unweighted_events.lhe","Table"];]

Timing[eventListp=ReadME[testp]; ]

Nev = Length[eventListp]-1

tsp={eventListp[[1]],eventListp[[2]]};
EventPrint /@ %;

eventListOutLp=SetOut[eventListp];

{eventListOutLp[[1]],eventListOutLp[[2]]};
EventPrint /@ %;

(* plot muon pair invariant mass *)

InvMass[event_] := Sqrt[2*(enOf[event[[1]]]*enOf[event[[2]]]-ThreeVectorFrom[event[[1]]].ThreeVectorFrom[event[[2]]])];

This next line might be a little confusing:
We want to apply the function InvMass(event) to a list of particles.  To do this, I use the Map function rather than a for loop.  Map applies its first argument to every element of the vector in its second argument and stores the results in a vector.  The # is replaced by each of the particles (and an & is needed to let Mathematica know to do this).  The particles we want are the first ones in the "out" list for events 1-10000, as indicated by the second argument for Map.  The vector is then plotted using the Histogram function trivially.

IM = Map[InvMass[#] &, eventListOutLp[[1;;Nev]]];

Histogram[IM, {80,100,1.0}]

(* plot rapidity of the Z = sum of muon rapidities *)

EtaZ = Map[etaOf[#]&,eventListOutLp[[1;;Nev,1]]]+Map[etaOf[#]&,eventListOutLp[[1;;Nev,2]]];

data1 = Histogram[EtaZ,{-5,5,0.5} ]

(* rapidity of individual muons *)

EtaMU = Map[etaOf[#]&,eventListOutLp[[1;;Nev,1]]];

data2 = Histogram[EtaMU,{-5,5,0.5} ]