Flip’s LHCO tools

Flip Tanedo (pt267@cornell.edu)
6 November 2012

Based on the Harvard Chameleon notebook and the Cornell MadEvent Analysis notebook. 
What it does: takes the output LHCO file from the MG5/Pythia-PGS pipeline and provides tools for basic data analysis.

How it works: run this notebook. In a separate notebook, do your analysis.

(* Input: lhco file from PGS *)
(* Output: list of row data, not yet processed by event *)
(* USAGE: eventList = ReadLHCO[Import["FILENAME.lhco","Table"]]; *)
ReadLHCO[rawInput_ ] := 
(
(* Define the column title *)
columnTitle={"#","typ","eta","phi","pt","jmas","ntrk","btag","had/em","dum1","dum2"};

(* Define position of column titles *)
pos=Position[rawInput,columnTitle][[1,1]]; 

(* Delete headers using the 'drop' command *)
noheader=Drop[rawInput,pos];

(* Split by event, use length of each line to demarcate each event *)
(* This works because the events are labelled by a 'trigger' line of length 3 *)
byevent=SplitBy[noheader,Length];

(* Remove trigger lines *)
(* These spoil the list *)
DeleteCases[byevent,{{0,_,_}}]
);


(* PARTICLE LIST *)
(* For Event Print *)

particle[0] = "γ";
particle[1] = "e";
particle[2] = "μ";
particle[3] = "(Subscript[tau, had])";
particle[4] = "jet";
particle[6] = "MET";

(* PRINT EVENT *)
(* Input: an event (list of n-tuples ending with MET) *)
(* Output: human-readable output of the event data *)
EventPrint[ event_ ] := ( Print[Join[{{"#","Type", "η", "ϕ", "(Subscript[p, T])", "jet mass", "ntrack", "btag", "had/em", "dummy", "dummy"}}, event/.{x1_,x2_,x3_,x4_,x5_,x6_,x7_,x8_,x9_,x10_,x11_}->{x1,particle[x2],x3,x4,x5,x6,x7,x8,x9,x10,x11}]//MatrixForm]); 

(* Get MET *)
(* Input: an event (list of n-tuples ending with MET) *)
(* Output: the scalar MET for the event *)
MET[event_] :=event[[Length[event]]][[5]];

(* Get list of leptons *)
(* Input: an event (list of n-tuples ending with MET) *)
(* Output: an 'event' containing only leptons *)
LeptonList[event_] := (
DeleteCases[event,Except[
{_,2,_,_,_,_,_,_,_,_,_}| {_,1,_,_,_,_,_,_,_,_,_}]]
);

(* Number of leptons *)
(* Input: an event (list of n-tuples ending with MET) *)
(* Output: number of leptons in the event *)
NLeptons[event_]:=(
Length[LeptonList[event]]
);

(* Sort by Pt *)
(* Input: an event (list of n-tuples ending with MET) *)
(* Output: event listed in order of pT *)
SortByPt[event_]:=Sort[event,#1[[5]]>#2[[5]]&]

(* Get list of jets (including tau jets) *)
(* Input: an event (list of n-tuples ending with MET) *)
(* Output: an 'event' containing only jets and hadronically decaying τs *)
HadronList[event_] := (
DeleteCases[event,Except[
{_,3,_,_,_,_,_,_,_,_,_}| {_,4,_,_,_,_,_,_,_,_,_}]]
);

(* Get HT *)
(* Input: an event *)
(* Output: the HT for the event, including hadronic τs *)
HT[event_]:=(
Sum[Abs[HadronList[event][[i]][[ 5]]],{i,1,Length[HadronList[event]]}]
);

(* Get Lepton ID *)
(* Input: an n-tuple *)
(* Output: the lepton id (+2 = muon, -1 = positron, etc.) *)
LeptonID[ntuple_]:=ntuple[[2]] ntuple[[7]];

(* Get number of Drell Yan pairs *)
(* Input: an event *)
(* Output: number of OS, SF pairs *)
DY[event_]:=
(
(* List of signed lepton IDs*)
IDleps=LeptonID/@LeptonList[event];

(* Return number of pairs *)
Min[Count[IDleps,-1.],Count[IDleps,1.]]+Min[Count[IDleps,-2.],Count[IDleps,2.]]
);