Chameleon 0.51 BETA -- 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, Philip Schuster, Jesse Thaler, Natalia Toro, Lian-Tao Wang 

Harvard University, Jefferson Laboratory
Oct 25, 2005

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

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

Eliminate soft objects to be excluded from all analysis

$PTCUT = 25;
$PTCUTJET = 50;
$ETACUT = 2.0;
$ETACUTJET = 3.0;
$MISSINGPTCUT = 50;
$TOTALPTCUT = 0;

\(MakeMasterCuts[] := \((
If[\(! ValueQ[ETACUT]\), 
          ETACUT = $ETACUT]; 
If[\(! ValueQ[PTCUT]\), 
          PTCUT = $PTCUT]; 
If[\(! ValueQ[ETACUTJET]\), 
          ETACUTJET = $ETACUTJET]; 
If[\(! ValueQ[
              PTCUTJET]\), 
          PTCUTJET = $PTCUTJET]; 
If[\(! ValueQ[
              MISSINGPTCUT]\), 
          MISSINGPTCUT = $MISSINGPTCUT]; 
If[\(! ValueQ[
              TOTALPTCUT]\), 
          TOTALPTCUT = $TOTALPTCUT]; 

goodPhoton = {_, 
            0, _?\((Abs[#] < ETACUT &)\), _, _?\((# > 
                    PTCUT &)\), _, _, _}; 
\ \ \ \ goodElectron = {_, 
            1, _?\((Abs[#] < ETACUT &)\), _, _?\((# > 
                    PTCUT &)\), _, _, _}; 
\ \ \ \ goodMuon = {_, 
            2, _?\((Abs[#] < ETACUT &)\), _, _?\((# > 
                    PTCUT &)\), _, _, _}; \ 
\ \ \ \ goodTau = {_, 
            3, _?\((Abs[#] < ETACUT &)\), _, _?\((# > 
                    PTCUT &)\), _, _, _}; 
\ \ \ \ goodJet = {_, 
            4, _?\((Abs[#] < ETACUTJET &)\), _, _?\((# > 
                    PTCUTJET &)\), _, _, _}; 
\ \ \ \ goodMissingPT = {_, 
            6, _, _, _?\((# > 
                    MISSINGPTCUT &)\), _, _, _}; 
\ \ \ \ goodParticle\ \  = \ \((goodPhoton\  | \ goodElectron\  | \ goodMuon\  | \ goodTau\  | \ 
              goodJet\  | \ 
              goodMissingPT)\);
)\)
\)
ApplyMasterCuts[
          splitData_] :=  (MakeMasterCuts[]; 
DeleteCases[
            DeleteCases[splitData, 
              Except[goodParticle], {2}] , {} | _?(EffMass[#] < 
                      TOTALPTCUT &)]);

Master importing functions

FCheck[fname_] := 
      If[MatchQ[FileType[fname], File], fname, 
        Experimental`FileBrowse[]]
GetFromFile[fname_] := 
        ApplyMasterCuts[
          SplitData[
            Import[fname, "Table"]]];
GetFromFile[$Failed] := Null;
LoadEvents[fname_] := 
        GetFromFile[FCheck[fname]];
\(\(LoadEvents[fname_, nfiles_] := 
        Module[{fhead}, 
fhead = 
            FCheck[fname]; 
fhead = 
            StringReplace[
              fhead, ".dat.0" -> ".dat"]; 
Flatten[
            Join[Table[
                GetFromFile[fhead <> "." <> ToString[i]], {i, 0, 
                  nfiles - 1}]], 1]];\)\)

Particle-matching Patterns

oElectron  =  oE = {_,  
            1, _,  _,  _,  _,  _,  _};
oElectronPlus  =  oEP  =  {_,  1, _,  _,  _,  
            1.0,  _,  _} ;
oElectronMinus  =  oEM  =  {_,  
            1, _,  _,  _,  -1.0,  _,  _} ;

oMuon  =  oμ = {_,  
            2, _,  _,  _,  _,  _,  _};
oMuonPlus  =  oμP = {_,  2, _,  _,  _,  
            1.0,  _,  _} ;
oMuonMinus  =  oμM = {_,  
            2, _,  _,  _,  -1.0,  _,  _} ;

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

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, _,  _,  _,  _,  _,  _};
oBJet  = oB =  {_,  4, _,  _,  _,  _,  _,  
            1.0};
oNonBJet  =  oNonB = {_,  4, _,  _,  _,  _,  _,  
            0.0};
oMissingPT  =  oMPT =  {_,  
            6, _,  _,  _,  _,  _,  _};

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

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

Hardness-selecting patterns

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

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
etaOf[{a_}] :=  etaOf[a]
phiOf[{a_}] :=  phiOf[a]
ptOf[{a_}] :=  ptOf[a]

    
FourVectorFrom[{_,  _, eta_,  phi_,  
            pt_,  ___}]  :=  {pt Cosh[eta],  pt Sinh[eta],  
          pt Cos[phi], pt Sin[phi] };
FourVectorFrom[{_, 4, eta_,  phi_,  pt_,  
            jetM_, ___}]  :=  {(Sqrt[(pt^2) Cosh[eta]^2 + (jetM^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

InvMass[objList_]  :=  
        FourLength[Apply[Plus, Map[FourVectorFrom, objList]]];
EffMass[objList_]   :=  
        Apply[Plus,  Map[ptOf,  objList]];
MissingPT[objList_] := 
        If[(NumOf[oMissingPT]~eq~1)[objList], 
          ptOf[HGet[oMissingPT, 1][objList]], 
          0];
DeltaR[{obj1_, 
            obj2_}] := (Sqrt[(((etaOf[obj1] - etaOf[obj2]))^2) + (((phiOf[obj1] - phiOf[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;

Display an event record with headings

\(ShowEvent[evtList_, \ num_] := \ 
      MatrixForm[
        Join[{{"#", "id", "eta", "phi", "pt", "jmas", "ntrack", "btag"}}, Identity[evtList[\([num]\)]]]]\)

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

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

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

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

Speedup routines (NOT YET EXHAUSTIVELY TESTED)

(modified Oct. 26)

SetAttributes[fastMake, HoldFirst];
SetAttributes[NumOf,  HoldFirst];

    
SetAttributes[SpeciesCount, HoldAll];
fastReset[] := (
ClearAll[
          fast]; 
SetAttributes[fast, 
          HoldFirst]; 
fast /: LeptonContent[typelist_][
            fast[list_, num_]]   := fast /: LeptonContent[typelist][
              fast[list, num]]  = LeptonContent[typelist][
              list[[num]]]; 
fast /: NumOf[obj_][
            fast[list_, num_]] :=  fast /: NumOf[obj][
              fast[list, num]]  =  NumOf[obj][
              list[[num]]]; 

fast /:  eq[attr_, n_][fast[list_, num_]]  :=  
          attr[fast[list, num]] == 
            n; 
fast /:  gt[attr_, n_][
            fast[list_, num_]]  :=  
          attr[fast[list, num]] > 
            n; 
fast /:  lt[attr_, n_][
            fast[list_, num_]]  :=  
          attr[fast[list, num]] < 
            n; 
fast /:  ge[attr_, n_][
            fast[list_, num_]]  :=  
          attr[fast[list, num]] ≥  
            n; 
fast /:  le[attr_, n_][
            fast[list_, num_]]  :=  
          attr[fast[list, num]] ≤  
            n; 

fast  /:  and[f_][
            fast[list_, num_]] := 
          f[fast[list, num]]; 
fast  /:  and[f_, g__][
            fast[list_, num_]] := 
          f[fast[list, num]] && and[g][
              fast[list, num]]; 
fast  /:  or[f_][
            fast[list_, num_]] := 
          f[fast[list, num]]; 
fast  /:  or[f_, g__][
            fast[list_, num_]] := 
          f[fast[list, num]] || or[g][
              fast[list, num]]; 
fast  /:  join[f_][
            fast[list_, num_]] := 
          f[fast[list, num]]; 
fast  /:  join[f_, 
              g__][fast[list_, num_]] := 
          Join[f[fast[list, num]], 
            join[
              g[fast[list, 
                  num]]]]; 

fast  /:  
          func_[fast[list_, num_]]  :=  
          func[list[[num]]];
)
fastReset[];
fastMake[list_]  :=  
        Table[fast[list, i],  {i, 1,  Length[list]}];

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.