Chameleon 1.02 -- A black-box analysis package in Mathematica for the LHC Olympics

Harvard Blackbox Team:
Nima Arkani-Hamed, Matthew Baumgart, Clifford Cheung, Liam Fitzpatrick, Tom Hartman, Can Kilic, Aaron Pierce, Philip Schuster, Jesse Thaler, Natalia Toro, Lian-Tao Wang 

Harvard University, Jefferson Laboratory
July 9, 2006

Setup functions (file-loading)

Suppress warnings

Off[General::spell1];
Off[General::spell1];

Load packages

<< Graphics`Graphics`
<< Graphics`MultipleListPlot`
<< Statistics`DataManipulation`
<< DiscreteMath`Combinatorica`

Function definition for version < 5.1 (for backwards-compatibility) added Oct 20

If[$VersionNumber < 5.1, 
      Except[p1_, 
          p2_] := (p2?(Not[
                  MatchQ[#, p1]] &)); 
Except[
          p1_] := _?(Not[MatchQ[#, p1]] &);]

Check Chameleon Version Settings

\(\(If[Or[$PGSVersion == 4, \(! ValueQ[$PGSVersion]\)], 
        Print["Reading data as PGS Version 4...set $PGSVersion=3 for compatibility with older (6-column, no-trigger) data format"]; $PGSVersion = 4, 
If[$PGSVersion == 3, 
          Print["Reading data as PGS Version 3...set $PGSVersion=4 for compatibility with new (11-column, triggered) data format"]]];\)\)

DefPGS4[f_, fdef_] := 
        If[$PGSVersion == 4, fdef, Clear[f]];

Pre-Processing Function Definitions

Split into events (automatically delete comment line, if there is one)

SplitData[rawObjList_] := 
        Split[Drop[rawObjList, If[NumberQ[rawObjList[[1, 1]]], 0, 1]], 
          Less[First[#1], First[#2]] &];

Master importing functions

FCheck[fname_] := 
      If[MatchQ[FileType[fname], File], fname, 
        Experimental`FileBrowse[]]
GetFromFile[fname_] := 
        SplitData[Import[fname, "Table"]];
GetFromFile[$Failed] := Null;
LoadEvents[fname_] := 
        GetFromFile[FCheck[fname]];
\(\(LoadEvents[fname_, nfiles_?NumberQ] := 
        Module[{fhead}, 
fhead = 
            FCheck[fname]; 
fhead = 
            StringReplace[
              fhead, ".dat.0" -> ".dat"]; 
Flatten[
            Join[Table[
                GetFromFile[fhead <> "." <> ToString[i]], {i, 0, 
                  nfiles - 1}]], 1]];\)\)
Options[LoadEvents] = {Criteria -> (True &), 
          MaxEvents -> ∞};
TakeUpTo[list_, n_] := 
        Take[list, Min[Length[list], n]];
LoadEvents[namespec_, LoadEventsOptions___] := 
        TakeUpTo[
          Select[LoadEvents[namespec], Criteria /. {LoadEventsOptions} /. 
              Options[LoadEvents]], MaxEvents /. {LoadEventsOptions} /. 
            Options[LoadEvents]];

Particle-matching Patterns

If[$PGSVersion == 
        4, 
oElectron  =  oE = {_,  
            1, _,  _,  _,  _,  _,  __}; 
oElectronPlus  =  oEP  =  {_,  
            1, _,  _,  _,  _,  
            1.0,  __} ; 
oElectronMinus  =  oEM  =  {_,  1, _,  _,  _,  _,  -1.0,  __} ;]

If[$PGSVersion == 
        3, 
oElectron  =  oE = {_,  
            1, _,  _,  _,  _,  _,  __}; 
oElectronPlus  =  oEP  =  {_,  1, _,  _,  _,  
            1.0, _,  __} ; 
oElectronMinus  =  oEM  =  {_,  1, _,  _,  _, -1.0, _,  __} ;]

If[$PGSVersion == 
        4, 
oMuon  =  oμ = {_,  
            2, _,  _,  _,  _,  _,  __}; 
oMuonPlus  =  oμP = {_,  2, _,  _,  _,  _,  
            1.0,  __} ; 
oMuonMinus  =  oμM = {_,  2, _,  _,  _, _,  -1.0,  __} ;]

If[$PGSVersion == 
        3, 
oMuon  =  oμ = {_,  
            2, _,  _,  _,  _,  _,  __}; 
oMuonPlus  =  oμP = {_,  2, _,  _,  _,  
            1.0, _, __} ; 
oMuonMinus  =  oμM = {_,  2, _,  _,  _, -1.0, _, __} ;]

If[$PGSVersion == 
        4, 
oTau  =  oτ = {_,  
            3, _,  _,  _,  _,  _,  __}; 
oTauPlus3Prong  =  oτP3 = {_,  
            3, _,  _,  _, _, 
            3.0,  __} ; 
oTauMinus3Prong  =  oτM3 = {_,  
            3, _,  _,  _, _,  -3.0,  __} ; 
oTauPlus1Prong  =  oτP1 = {_,  
            3, _,  _,  _,  _, 
            1.0,  __} ; 
oTauMinus1Prong  =  oτM1 = {_,  
            3, _,  _,  _, _, -1.0,  __} ; 
oTauPlus  =  oτP = 
          oτP1 | 
            oτP3 ; 
oTauMinus  =  oτM = 
          oτM1 | 
            oτM3;]

If[$PGSVersion == 
        3, 
oTau  =  oτ = {_,  
            3, _,  _,  _,  _,  _,  __}; 
oTauPlus =  oτP = {_,  3, _,  _,  _,  
            1.0, _,  __} ; 
oTauMinus  =  oτM = {_,  3, _,  _,  _, -1.0, _, __} ;
]

oLepton   =  oL = (oElectron  |  oMuon  |  
              oTau);
oLeptonPlus   =  oLP = (oElectronPlus  |  oMuonPlus  |  
              oTauPlus);
oLeptonMinus  =  oLM = (oElectronMinus  |  oMuonMinus  |  
              oTauMinus);
oLightLepton   =  oLL = (oElectron  |  
              oMuon );
oLightLeptonPlus   =  oLLP = (oElectronPlus  |  
              oMuonPlus );
oLightLeptonMinus  =  oLLM = (oElectronMinus  |  
              oMuonMinus );

oLeptonType[-3]  =   oTauMinus;
oLeptonType[-2]  =   oMuonMinus;
oLeptonType[-1]  =   oElectronMinus;
oLeptonType[0]  =    object_  /;  
          False;
oLeptonType[1]  =   oElectronPlus;
oLeptonType[2]  =   oMuonPlus;
oLeptonType[3]  =   oTauPlus;

oPhoton  =  oγ = {_,  
            0, _,  _,  _,  _,  _,  __};
oJet  =  {_,  
          4, _,  _,  _,  _,  _,  __};
oBJetTight = oBT = oBJet2 = oB2 = {_,  
                4, _,  _,  _,  _,  _,  
                2.0, ___};
oBJetLoose = oBL = oBJet1  = oB1 = {_,  
                  4, _,  _,  _,  _,  _,  1.0, ___} | 
                oB2;
oB = oBJet = 
          If[$PGSVersion == 4, oB1, 
            oB1];
oNonBJet  =  oNonB = {_,  4, _,  _,  _,  _,  _,  
            0.0, ___};
oNonLooseBJet  =  oNonB = {_,  4, _,  _,  _,  _,  _,  
            0.0, ___};
oNonTightBJet = oNonB | oB1;
oMissingPT  =  oMPT =  {_,  
            6, _,  _,  _,  _,  _,  __};

oReconstructed = oRec = oR = {_, 
              7, _, _, _, _, _, __};

oAny  =  {_, _, _, _, _, _, _, __};

oNot[oType_] := Except[oType, oAny];

Hardness-selecting patterns

oHard[oType_, pthard_]  = 
        oType?((ptOf[ #] ≥ 
                  pthard) &);
oLessHard[oType_, pthard_]  = 
        oType?((ptOf[#] < pthard) &);

oTrigger = {0, _, __};

Cleaning functions

Enforcing muon isolation

ptIsoOf[{x_}] :=  ptIsoOf[x]
etRatioOf[{x_}] :=  etRatioOf[x]
ptIsoOf[{_, 2, _, _, _, _, _, _, ptiso_, __}]  :=  
        N[IntegerPart[ptiso]];
etRatioOf[{_, 2, _, _, _, _, _, _, etratio_, __}]  :=  
        N[FractionalPart[etratio]];
oIsolatedMuon  =  {_, 
            2, _, _, _, _, _, _, _, __}?((ptIsoOf[#] < 5.0 &&  
                  etRatioOf[#] <  .1*9/
                        8) &);
oMergeableMuon  =  {_, 
            2, _, _, _, _, _, _, _, __}?((ptIsoOf[#] ≥ 
                    5.0 ||  etRatioOf[#] ≥  .1*9/
                        8) &);

FindListMax[list_,  function_]  :=  
      Module[{tempmaxindex,  tempmaxvalue, 
          tempfunctionvalue}, 
tempmaxindex  =  
          1; 
tempmaxvalue  =  
          function[
            list[[1]]]; 
Do[
tempfunctionvalue  =  
            function[list[[i]]]; 
If[
            tempfunctionvalue  >  tempmaxvalue,  tempmaxindex  =  i;  
            tempmaxvalue  =  tempfunctionvalue,  
            Null]
, {i, 2, 
            Length[list]}]; 
{tempmaxindex, tempmaxvalue, 
          list[[tempmaxindex]]}
]


DeltaRToMuon[muon_, object_, Options___]   := 
        If[ MatchQ[object, oJet], DeltaR[muon, object], 
          99];

CheckMuonsInEvent[event_, Options___]  :=  
        Module[{muonlist, tempevent, themuon, thejetinfo, thejet, 
            thejetindex}, 
muonlist  =  
            PTSort[HGet[oMergeableMuon, All][
                event]]; 
tempevent  =  
            event; 
Do[
themuon =  
              muonlist[[i]]; 
thejetinfo  =  
              FindListMax[
                tempevent,  -DeltaRToMuon[themuon, #, 
                      Options] &]; 
thejetindex  =  
              thejetinfo[[1]]; 
thejet  =  
              thejetinfo[[3]]; 
If[
              MatchQ[thejet, 
                oJet] , 
tempevent[[thejetindex]]  =  ReconstructParent[{themuon, thejet}, Type ->  4,  
                  ObjectCounter ->  thejet[[1]] , 
                  NTrack -> thejet[[7]] + 0.1, 
                  Btag -> thejet[[8]], 
                  HadEmFraction -> 
                    thejet[[9]]]; 
tempevent  = 
                DeleteCases[tempevent, 
                  themuon], 
tempevent  = 
                DeleteCases[tempevent, 
                  themuon]
], 
{i, 
              Length[muonlist], 
              1, -1}]; 
tempevent
];

ApplyMuonIsolation[evtList_,  
        Options___] :=  CheckMuonsInEvent[#, Options] &  /@  
        evtList 

RemoveTriggerObject[event_] := 
DeleteCases[event, 
        oTrigger]

ShowEvent[x_]  :=  ShowEvent[{x}, 1]

Stripping soft objects (TO BE ADDED)

Analysis Functions

Object level

Query function (boolean from pattern-testing)

Q[patt_][obj_] := MatchQ[obj, patt];

Kinematics

FourVectorFrom definition for jets changed Oct. 26

etaOf[{_,  _, eta_,  phi_,  pt_,  ___}] :=  
      eta
phiOf[{_,  _, eta_,  phi_,  pt_,  ___}] :=  
      phi
ptOf[{_,  _, eta_,  phi_,  pt_,  ___}] :=  
      pt
nTracksOf[{_, 4, _, _, _, _, nTrack_, ___}] := 
        Floor[nTrack];
nMergedMuonsOf[{_, 4, _, _, _, _, 
            nTrack_, ___}] := (nTrack - 
              Floor[nTrack])*10;
hadEmFractionOf[{_, _, _, _, _, _, _, _, fraction_, ___}] := 
        fraction;

    
etaOf[{a_}] :=  etaOf[a]
phiOf[{a_}] :=  phiOf[a]
ptOf[{a_}] :=  ptOf[a]

    
FourVectorFrom[{_,  _, eta_,  phi_,  pt_,  
            mass_, ___}]  :=  {(Sqrt[(pt^2) Cosh[eta]^2 + (mass^2)]),  
          pt Sinh[eta],  pt Cos[phi], 
          pt Sin[phi] };
FourLength[{pe_, pz_, px_, py_}]  :=  
        Sqrt[Max[pe^2  -  pz^2  -  px^2  -  py^2, 0.0]];

Multi-Object level

Sorting objects

by pt:

PTSort[objList_] := 
      Sort[objList, Greater[ptOf[#1], ptOf[#2]] &]

by any function (can also be used to sort ktuples)

FSort[func_][objList_] := 
      Sort[objList, Greater[func[#1], func[#2]] &]

Kinematic quantities for a list of objects

FourVectorSum[objList_] := 
        Apply[Plus, 
          Map[FourVectorFrom, objList]];
InvMass[objList_]  :=  
        FourLength[FourVectorSum[RemoveTriggerObject[objList]]];
EffMass[objList_]   :=  
        Apply[Plus,  
          Map[ptOf,  
            RemoveTriggerObject[objList]]];
MissingPT[objList_] := 
        If[(NumOf[oMissingPT]~eq~1)[objList], 
          ptOf[HGet[oMissingPT, 1][objList]], 
          0];
DeltaPhi[obj1_, obj2_] := 
        Module[{dϕ}, 
          dϕ = phiOf[obj1] - phiOf[obj2]; 
Min[
            Abs[dϕ], 
            2 π - Abs[dϕ]]];
DeltaPhi[{a_}] := DeltaPhi[a];
DeltaR[obj1_, 
          obj2_] := (Sqrt[(((etaOf[obj1] - etaOf[obj2]))^2) + DeltaPhi[obj1, obj2]^2]);
DeltaR[{a_}] := DeltaR[a];
InvMass[{}] = 0;
EffMass[{}] = 0;

Counting objects

NumOf[objType_][objList_] := 
        Count[objList, objType];

Demand that a numerical function of an object list satisfy (in)equality

eq[attr_, n_][objList_] := attr[objList] == n;
gt[attr_, n_][objList_] := attr[objList] > n;
lt[attr_, n_][objList_] := attr[objList] < n;
ge[attr_, n_][objList_] := 
        attr[objList] ≥ n;
le[attr_, n_][objList_] := 
        attr[objList] ≤ n;

Ktuple testing and processing

ListAlternative[p_] := 
      Alternatives @@ 
        If[MatchQ[Head[p], PatternTest], p[[1]], p]

Select kTuples from an object list that satisfy a criterion

KTuples[k_, crit__][objList_] := 
        Select[KSubsets[objList, k], crit];

Select ordered kTuples matching a specific pattern and of given length.  Note you can also implement criteria with patt?crit (see Mathematica help)

OrdTuples[k_, patt_][objList_] := 
        Cases[Flatten[Permutations /@ KSubsets[objList, k], 1], 
          patt];

Note: if necessary for speed, it's easy to add a criterion argument that's tested on the subset level so most subsets aren't even permuted.

Does any kTuple satisfy my criterion?

KTest[k_, crit_][objList_] := 
        MemberQ[KSubsets[objList, k], _?crit];

OrdTest[k_, patt_][objList_] := 
        MemberQ[Flatten[Permutations /@ KSubsets[objList, k], 1], 
          patt];

Select a tuple by sorting in order of sortf (like HGet)

PickATuple[TupleFunc_, sortf_, n_][evt_] := 
        Take[FSort[sortf][TupleFunc[evt]], 
          n];
PickATuple[TupleFunc_, sortf_, {n__}][evt_] := 
        Part[FSort[sortf][TupleFunc[evt]], {n}];

Event level

Trivial criterion

all[evt_] := True;

Trigger information selection

DefPGS4[TriggerWord, 
        TriggerWord[evt_] := evt[[1, 3]]];

BinaryWord[list_?ListQ] := 
        Plus @@ ((2^(# - 1) &) /@ 
              list);
BinaryWord[word_] := word;
BinaryExpand[
          word_] := Position[Reverse[IntegerDigits[word, 2]], 1][[All, 
            1]];

DefPGS4[
        PassesTriggers, PassesTriggers[testWord_][evt_] := 
          If[BitAnd[BitNot[TriggerWord[evt]], BinaryWord[testWord]] == 
              0, True, False]];

DefPGS4[
        ExactTriggers, ExactTriggers[testWord_][
            evt_] := (TriggerWord[evt] == 
              BinaryWord[testWord])];

Display an event record with headings

\(ShowEvent[evtList_, \ num_, 
        ShowEventOptions___] := \ \((If[$PGSVersion == 3, 
          MatrixForm[
            Join[{{"#", "id", "eta", "phi", "pt", "jmas", "ntrack", "btag"}}, 
              Identity[
                evtList[\([num]\)]]]], 
If[$PGSVersion == 4, \ 
            Module[{tempObjects}, 
              tempObjects\  = \ 
                DeleteCases[
                  evtList[[num]], \ 
                  oTrigger]; 
              Print["Triggers Passed:", 
                BinaryExpand[
                  TriggerWord[
                    evtList[[num]]]]]; 
              MatrixForm[
                Join[{{"#", "typ", "eta", "phi", "pt", "jmas", "ntrk", "btag", \ "had/em", \ "dum1", \ "dum2", \ "E", \ "pz", \ "px", \ "py"}}, 
                  MapThread[
                    Join, {tempObjects, \ 
                      FourVectorFrom\  /@ \ tempObjects}]]]]]])\)\)

Demand a particular set of leptons (other particles unconstrained)

LeptonIndexCount[typelist_, n_][
        evt_] := (NumOf[oLeptonType[n]]~eq~Count[typelist, n])[
        evt]

LeptonContent[typelist_][
          evt_] := LeptonIndexCount[typelist, -3][
            evt] && 
LeptonIndexCount[typelist, -2][
            evt] && 
LeptonIndexCount[typelist, -1][
            evt] && 
LeptonIndexCount[typelist, 1][
            evt] && 
LeptonIndexCount[typelist, 2][
            evt] && 
LeptonIndexCount[typelist, 3][
            evt];

OSDiLightLepton[
          evt_] := or[LeptonContent[{1, -1}], 
            LeptonContent[{2, -2}]][evt];
SSDiLightLepton[
          evt_] := or[LeptonContent[{1, 1}], LeptonContent[{2, 2}], 
            LeptonContent[{-1, -1}], 
            LeptonContent[{-2, -2}]][
          evt];
OSDiLepton[
          evt_] := or[LeptonContent[{1, -1}], 
            LeptonContent[{2, -2}], LeptonContent[{3, -3}]][
          evt];
SSDiLepton[
          evt_] := or[LeptonContent[{1, 1}], LeptonContent[{2, 2}], 
            LeptonContent[{3, 3}], LeptonContent[{-1, -1}], 
            LeptonContent[{-2, -2}], 
            LeptonContent[{-3, -3}]][
          evt];
OSDL = OSDiLepton;
SSDL = SSDiLepton;
OSDLL = OSDiLightLepton; SSDLL = SSDiLightLepton;

Choose hardest n objects of given type, or a list of objects of given type, indexed by hardness

Usage:
	HGet[patt,n] returns hardest n objects matching patt
	HGet[patt,-n]    "      softest n     "          "           "    
	HGet[patt,All] returns all objects in an event matching patt
	HGet[patt,{k1,k2...,km}] returns k1'th,k2'th,..km'th hardest objects matching patt
	HGet[{{patt1,n},{patt2,{k}},{patt3,{-n}}}] joins what you'd get from acting on each pair separately

HGet[patt_, n_][evt_] := 
        Take[PTSort[Cases[evt, patt]], n];
HGet[patt_, {n__}][evt_] := 
        Part[PTSort[Cases[evt, patt]], {n}];
HGet[patt_, All][evt_] := 
        Cases[evt, patt];
HGet[{{patt1_, sel1_}, MoreArgs___}][evt_] := 
        Join[HGet[patt1, sel1][evt], HGet[{MoreArgs}][
            evt]];
HGet[patt_, UpTo[n_]][evt_] := 
        If[Abs[n] < Length[Cases[evt, patt]], 
          Take[PTSort[Cases[evt, patt]], n], 
          Cases[evt, patt]];
HGet[{}][evt_] := List[];

GetAll[evt_] := evt;

Selection as in HGet but with arbitrary function used for sorting (can also be used on lists of ktuples, with function acting on tuples)

Helper function (Take that can handle requests for more elements than there are in a list) NOT USED OR DOCUMENTED

FGet[func_, patt_, n_][evt_] := 
        Take[FSort[func][Cases[evt, patt]], 
          n];
FGet[func_, patt_, {n__}][evt_] := 
        Part[FSort[func][
            Cases[evt, patt]], {n}];
FGet[func_, patt_, All][evt_] := 
        Cases[evt, patt];
FGet[func_, patt_, UpTo[n_]][evt_] := 
        If[Abs[n] < Length[Cases[evt, patt]], 
          Take[FSort[func][Cases[evt, patt]], n], 
          Cases[evt, patt]];
FGet[{{func1_, patt1_, sel1_}, MoreArgs___}][evt_] := 
        Join[FGet[func1, patt1, sel1][evt], HGet[{MoreArgs}][
            evt]];
FGet[{}][evt_] := List[];

Logical operators on functions (defined recursively)

These functions "distribute" the And, Or, and Join operations over function calls--i.e. and[f,g] is a function that returns true on x if f and g both return true.  
For Boolean functions f and g,
and[f,g][evt]=And[f[evt],g[evt]] = f[evt] && g[evt]
or[f,g][evt]=Or[f[evt],g[evt]]=f[evt] || g[evt]
join acts on functions that return lists, and returns a function that returns the union of the two lists f and g would have produced.
join[f,g][evt]=Join[f[evt],g[evt]]

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]]];

Mulit-Event Level (high-level routines)

Counting helper function

CriterionCount[evtList_, crit_] := 
Count[
          evtList, _?crit];

Count number of events with k particles of type patt, for k from 0 to Nmax

SpeciesCount[evtList_, patt_, Nmax_] := 
MatrixForm[
          Transpose[
            Table[{i, CriterionCount[evtList, NumOf[patt]~eq~i]}, {i, 0, 
                Nmax}]]];

LEPLIST = {-1, -2, -3, 1, 2, 
          3};

    
\(\(\(SingleLeptonCount[evtList_] := 
        Module[{TMPLIST}, 
TMPLIST = 
            Select[evtList, \((NumOf[oLepton]~eq~
                  1)\)]; 
MatrixForm[{{"e-", "μ-", "τ-", "e+", "μ+", "τ+"}, 
Map[CriterionCount[TMPLIST, LeptonContent[{#}]] &, 
                LEPLIST]}]];\)\(
\)
    \)\)
\(\(DiLeptonCount[evtList_] := 
        Module[{TMPLIST, F1, UTfnc}, 
TMPLIST = 
            Select[evtList, \((NumOf[oLepton]~eq~2)\)]; 
F1[
              a_] := Abs[a] + 2\ Sign[a]; 
UTfnc[l_] := 
            F1[l[\([1]\)]] ≤ 
              F1[l[\([2]\)]]; 
MatrixForm[
            Transpose[
              Join[{{"-", "e-", "μ-", "τ-", "e+", "μ+", "τ+"}}, 
                Transpose[
                  Join[{{"e-", "μ-", "τ-", "e+", "μ+", "τ+"}}, 
Map[
                      If[UTfnc[#], CriterionCount[TMPLIST, LeptonContent[#]], 
                          0] &, Outer[List, LEPLIST, 
                        LEPLIST], {2}]]]]]]];\)\)

Advanced counting (MAY BECOME OBSOLETE, NOT DOCUMENTED)

AdvancedCount[evtList_, crit_, objsel_, boolFunc_] := 
        Count[objsel /@ Select[evtList, crit], 
          ev_ /; boolFunc[ev]];

Histogramming

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

    
Freq[evtList_, crit_, objsel_, plotfunc_, {min_, max_, nbins_}] := 
        BinCounts[
          plotfunc /@ Flatten4[objsel /@ Select[evtList, crit]], {min, 
            max, ((max - min)/nbins)}];
Hist[freqcombo__, {min_, max_, nbins_}, HistogramOptions___] := 
        Histogram[
          freqcombo /. Bin[args__] -> Freq[args, {min, max, nbins}], 
          FrequencyData -> True, 
          HistogramCategories -> 
            Table[i, {i, min, max, ((max - min)/nbins)}], 
          HistogramOptions];
MakeFreq[freqcombo__, {min_, max_, nbins_}] := 
        freqcombo /. 
          Bin[args__] -> Freq[args, {min, max, nbins}];

Object/event modification functions

Create object

Options[NewObject] = {ForceMass -> False, NTrack -> 0, 
          Btag -> 0, HadEmFraction -> 0, Dummy1 -> 0, 
          Dummy2 -> 0, ObjectCounter -> -1};

ForceMassAdjust[{m_, pt_, eta_, phi_}, forcemass_] := {forcemass, pt, 
          eta, phi};

NewObject[Type_, {p0_, pz_, px_, py_}, NewObjectOptions___] := 
        Module[{mass, eta, phi, pt, 
            forcemass}, 
forcemass = ((ForceMass /. {NewObjectOptions}) /. Options[NewObject]); 
phi = 
            Mod[ArcTan[px, py], 
              2 π]; 
pt = (Sqrt[(py^2) + (px^2)]); 
eta = 
            ArcSinh[pz/
                pt]; 
mass = (Sqrt[(p0^2) - (pz^2) - (px^2) - (py^2)]); 
If[! MatchQ[forcemass, False], {mass, pt, eta, 
                phi} = ForceMassAdjust[{mass, pt, eta, phi}, 
                forcemass]]; 
(If[$PGSVersion == 
                    4, {ObjectCounter, Type, eta, phi, pt, mass, NTrack, 
                    Btag, HadEmFraction, Dummy1, Dummy2}, {ObjectCounter, 
                    Type, eta, phi, pt, mass, NTrack, 
                    Btag}] /. {NewObjectOptions}) /. 
            Options[NewObject]
];

ReconstructParent[{objects__}, Type -> type_:  7, Options___] := 
      NewObject[type, FourVectorSum[{objects}], Options]

Options[ReconstructParent] = Options[NewObject];

Reconstruction

Summary output routines (highest-level)

Inclusive count set

RunInclusiveCounts[evtList_] := 
      Print["\nPhotons: ", 
        SpeciesCount[evtList, oγ, 4], "   Jets: ", 
        SpeciesCount[evtList, oJet, 6], "   BJets: ", 
        SpeciesCount[evtList, oBJet, 4], "   Lepton: ", 
        SpeciesCount[evtList, oL, 4] , "  (OSDL: ", 
        CriterionCount[evtList, OSDiLepton], ", SSDL: ", 
        CriterionCount[evtList, 
          SSDiLepton], ")\nSingle- and Di- Lepton Counts", 
        SingleLeptonCount[evtList], 
        DiLeptonCount[
          evtList], "\n Single- and Di-Lepton Counts (events with no jets): \
", SingleLeptonCount[Select[evtList, NumOf[oJet]~eq~0]], 
        DiLeptonCount[Select[evtList, NumOf[oJet]~eq~0]]]

Changelog

v. 0.5 / Oct 20:

added fake definition of Except for compatibility with Mathematica 5.0.  Time required for loading data in earlier versions of Mathematica is greater than in 5.0 by a factor of 2.  (Natalia Toro)

v. 0.51 / Nov 1:

corrected FourVectorFrom to use jet invariant mass in calculating invariant mass for jet objects. (Natalia Toro)
changes to "fast" routines: "fast" information (results of previous NumOf evaluations that are stored to memory) is associated with the name of the raw data array it was called on, but is not erased when this array is cleared.  A partial fix, in the current version, is the fastReset[] function.  It must be run whenever you load new data into an old array name, so that "fast" results on the old data are not mis-interpreted as results for the new data.  Unfortunately this must be run manually, and destroys ALL "fast" data (no event data is lost irretrievably, only pre-calculated function calls are lost), which we hope to change soon.
corrected description of master cuts from "triggers" and changed default $TOTALPTCUT to 0 (which was used in the tutorial and Harvard BB results).
added null definitions of InvMass,EffMass: InvMass[{}]=0, etc.
HGet[patt,UpTo[n]] takes the n hardest objects of type patt, or all such objects if there are fewer than n.  HGet[patt,UpTo[-n]] works similarly, picking the softest n objects, or all all objects.  The analogous syntax for FGet has also been added.

v. 1.0/April 6, 2006:

Modified LoadEvents to take two options: MaxEvents (max. number of events to keep) and Criteria (function that must return true if event is to be kept).  This should be sped up in the future, but functionality should be compatible. 
Added $PGSVersion variable (default = 4) that knows what kind of PGS output format to expect.
New object types: BJet2 and BJet1 (BJet=BJet2 in v4-mode, BJet1 in v3-mode), 1- and 3-pronged taus,oLL/oLLP/oLLM for light leptons (not tau)
When $PGSVersion=4: Trigger operations
  -TriggerWord returns trigger word
  -BinaryWord converts list to binary word (BinaryExpand does the reverse)
  -ExactTriggers[word][event] returns true if event's trigger word is exactly word.
  -PassesTriggers[word][event] returns true if all of the triggers in word are passed by event.
  -Both PassesTriggers and ExactTriggers can take either a trigger word or a list of triggers to be converted to a word by BinaryWord (e.g. {2,1} -> 3)
Multi-Object: Fixed DeltaR, added DeltaPhi
Single-Object: Added CounterOf,ChargeOf,EMHadFractionOf,NTrackOf,NMuonOf (?)---
InclusiveCount looseB...
goodbye fastMake, goodbye ApplyMasterCuts, etc...
SumFourVectors, ReconstructParent, NewObject

1.02: July 9, 2006

Fixed muon isolation--now a mergeable muon is a muon that is not isolated, as it should be!  (Thanks to Steve Ellis)
Generalized oTrigger definition to an object starting with 0 with 3 OR MORE columns
RemoveTriggerObject[event] returns the event with no trigger object.  InvMass and EffMass automatically call RemoveTriggerObject, so they can act on a whole event. (Thanks to Brent Nelson)