MadEvent analysis routines

(* Based on the Chameleon package by Philip Shuster, Jesse Thaler, Natalia Toro *)
(* Modified by Maxim Perelstein and Andi Weiler, 2008-09 *)

(*format of the compulsory event vector*)
eventInfo={_,_,_,_,_,_};

ClearAll[particle,finout];

(*PDG numbers see e.g. http://home.fnal.gov/~maeshima/alignment/ORCA/PYTHIA_particle_codes.ps*)

particle[2] =u;
particle[-2] =uBar;

particle[1] =d;
particle[-1] =dBar;

particle[3] ="s";
particle[-3] ="sBar";

particle[4] ="c";
particle[-4] ="cBar";

particle[5] =b;
particle[-5] =bBar;

particle[6] =t;
particle[-6] =tBar;

particle[11] =(SuperMinus[e]);
particle[-11] =(SuperPlus[e]);

particle[12] ="(Subscript[⋁, e])";
particle[-12] ="((Subscript[⋁, e])^c)";

particle[13] ="(SuperMinus[μ])";
particle[-13]="(SuperPlus[μ])";

particle[14] ="(Subscript[⋁, μ])";
particle[-14] ="((Subscript[⋁, μ])^c)";

particle[16] ="(Subscript[⋁, τ])";
particle[-16] ="((Subscript[⋁, τ])^c)";

particle[21] =gluon;

particle[23] = Z;

particle[24] =(SuperPlus[W]);
particle[-24] =(SuperMinus[W]);

particle[1000006] =(Subscript[OverTilde[t], 1]);
particle[-1000006] =(SuperStar[(Subscript[OverTilde[t], 1])]);
particle[1000022] =(Subscript[OverTilde[χ], 0]);

finout[-1] = in;
finout[1] = out;
finout[2]  = decayed;

ReadME[rawInput_ ] := 
(
(*Search for beginning of events*)
pos=Position[rawInput,{"</init>"}][[1,1]]; 
(*Throw away crap at start of file, combine events in {} array and remove the XML commands and compulsory eventinfo *)
DeleteCases[Split[Drop[ rawInput, pos],#=!={"</event>"}&],{"<event>"}|{"</LesHouchesEvents>"}| {"</event>"} | eventInfo,2]
);

EventPrint[ event_ ] := ( Print[Join[{{pid, in/out,mother1,mother2,color1,color2, px,py,pz,p0,mass,6,hel}}, event/.{x1_,x2_,x3_,x4_,x5_,x6_,x7_,x8_,x9_,x10_,x11_,x12_,x13_}->{particle[x1],finout[x2],x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13}]//MatrixForm]); 

(* those are input or decayed particles*)
decayedOrIn = {_,-1,_,_,_,_,_,_,_,_,_,_,_}  | {_,2,_,_,_,_,_,_,_,_,_,_,_};

EventOut[ event_ ]  := (DeleteCases[event,decayedOrIn]);

SetOut[ event_ ]  := (DeleteCases[event,decayedOrIn,2]);

(* EffMassAll[objList_]  :=Plus @@ Map[ptOf,objList//Transpose,{2}]; *)

EffMass[event_]  :=Plus @@ Map[ptOf,event,{1}];

pT[event_]  :=Map[ptOf,event,{1}]//Flatten;
eta[event_] := Map[etaOf,event,{1}]//Flatten;
theta[event_] := Map[thetaOf,event,{1}]//Flatten;
threeVector[event_] := Map[ThreeVectorFrom,event,{1}];


ptOf[{_,_,_,_,_,_,px_,py_,___}]:=(Sqrt[px^2+py^2]);
enOf[{_,_,_,_,_,_,_,_,_,En_,___}]:=En;
thetaOf[{_,_,_,_,_,_,px_,py_,pz_,___}]:= ArcCos[pz/(Sqrt[px^2 + py^2+ pz^2])];
etaOf[{_,_,_,_,_,_,px_,py_,pz_,___}]:=- Log[ Abs[Tan[ ArcCos[pz/(Sqrt[px^2 + py^2+ pz^2])]/2]]]; 
FourVectorFrom[{_,_,_,_,_,_,px_,py_,pz_,En_,___}] := {En, px, py,pz };

FourLength[{pe_,pz_,px_,py_}] := Sqrt[Max[pe^2 - pz^2 - px^2 - py^2,0.0]];

ThreeVectorFrom[{_,_,_,_,_,_,px_,py_,pz_,___}] := { px, py,pz };

CosthetaTwoJet[{jet1_,jet2_}] := (a = ThreeVectorFrom[jet1]; b=  ThreeVectorFrom[jet2]; a.b/(Sqrt[a.a b.b]));

GetAll[evt_]:=evt;
all[evt_]:=True;

Freq[evtList_,crit_,objsel_,plotfunc_,{min_,max_,nbins_}]:=BinCounts[ Flatten[ plotfunc/@Flatten4[objsel/@Select[evtList,crit]]],{min,max,((max-min)/nbins)}];

Flatten4[List_]:=Flatten[List,Max[0,Depth[List]-4]];

HGet[patt_,All][evt_]:=Cases[evt,patt];

oJetlhe=  {2,1,_,_,_,_,_,_,_,_,_,_,_}  | {-2,1,_,_,_,_,_,_,_,_,_,_,_}|{21,1,_,_,_,_,_,_,_,_,_,_,_}  ;

oMissingEnergy = {8888,1,_,_,_,_,_,_,_,_,_,_,_}  ;

PTSort[event_]:=Sort[event,Greater[ptOf[#1],ptOf[#2]]&];


and[f_][evt_]:=f[evt];
and[f_,g__][evt_]:=f[evt]&&and[g][evt];
or[f_][evt_]:=f[evt];
or[f_,g__][evt_]:=f[evt]||or[g][evt];

join[f_][evt_]:=f[evt];
join[f_,g__][evt_]:=Join[f[evt],join[g[evt]]];