(*******************************************************************
This file was generated automatically by the Mathematica front end.
It contains Initialization cells from a Notebook file, which
typically will have the same name as this file except ending in
".nb" instead of ".m".

This file is intended to be loaded into the Mathematica kernel using
the package loading commands Get or Needs.  Doing so is equivalent
to using the Evaluate Initialization Cells menu command in the front
end.

DO NOT EDIT THIS FILE.  This entire file is regenerated
automatically each time the parent Notebook file is saved in the
Mathematica front end.  Any changes you make to this file will be
overwritten.
***********************************************************************)









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



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



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



If[Or[$PGSVersion\[Equal]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\[Equal]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\[Equal]4,fdef,Clear[f]];





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



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"\[Rule]".dat"];
      Flatten[
        Join[Table[GetFromFile[fhead<>"."<>ToString[i]],{i,0,nfiles-1}]],1]];
Options[LoadEvents]={Criteria\[Rule](True&),MaxEvents\[Rule]\[Infinity]};
TakeUpTo[list_,n_]:=Take[list,Min[Length[list],n]];
LoadEvents[namespec_,LoadEventsOptions___]:=
    TakeUpTo[Select[LoadEvents[namespec],
        Criteria/.{LoadEventsOptions}/.Options[LoadEvents]],
      MaxEvents/.{LoadEventsOptions}/.Options[LoadEvents]];



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

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

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

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

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

If[$PGSVersion\[Equal]3,
  oTau = o\[Tau]={_, 3,_, _, _, _, _, __};
  oTauPlus= o\[Tau]P={_, 3,_, _, _, 1.0,_, __} ;
  oTauMinus = o\[Tau]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\[Gamma]={_, 0,_, _, _, _, _, __};
oJet = {_, 4,_, _, _, _, _, __};
oBJetTight=oBT=oBJet2=oB2={_, 4,_, _, _, _, _, 2.0,___};
oBJetLoose=oBL=oBJet1 =oB1={_, 4,_, _, _, _, _, 1.0,___}|oB2;
oB=oBJet=If[$PGSVersion\[Equal]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];



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

oTrigger={0,_,_};





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\[Rule] 4, 
              ObjectCounter\[Rule] thejet[[1]] ,NTrack\[Rule]thejet[[7]]+0.1,
              Btag\[Rule]thejet[[8]],HadEmFraction\[Rule]thejet[[9]]];
          tempevent =DeleteCases[tempevent,themuon],
          tempevent =DeleteCases[tempevent,themuon]
          ],
        {i,Length[muonlist],1,-1}];
      tempevent
      ];

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

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









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





\!\(etaOf[{_, \ _, eta_, \ phi_, \ pt_, \ ___}] := \ eta\[IndentingNewLine]
  phiOf[{_, \ _, eta_, \ phi_, \ pt_, \ ___}] := \ phi\[IndentingNewLine]
  ptOf[{_, \ _, eta_, \ phi_, \ pt_, \ ___}] := \ pt\[IndentingNewLine]
  \(nTracksOf[{_, 4, _, _, _, _, nTrack_, ___}] := Floor[nTrack];\)\[IndentingNewLine]
  \(nMergedMuonsOf[{_, 4, _, _, _, _, nTrack_, ___}] := \((nTrack - Floor[nTrack])\)*10;\)\[IndentingNewLine]
  \(\(hadEmFractionOf[{_, _, _, _, _, _, _, _, fraction_, ___}] := fraction;\)\(\[IndentingNewLine]\)
  \)\[IndentingNewLine]
  etaOf[{a_}] := \ etaOf[a]\[IndentingNewLine]
  phiOf[{a_}] := \ phiOf[a]\[IndentingNewLine]
  \(\(ptOf[{a_}] := \ ptOf[a]\)\(\[IndentingNewLine]\)
  \)\[IndentingNewLine]
  \(FourVectorFrom[{_, \ _, eta_, \ phi_, \ pt_, \ mass_, ___}]\  := \ {\@\(pt\^2\ Cosh[eta]\^2 + mass\^2\), \ pt\ Sinh[eta], \ pt\ Cos[phi], pt\ Sin[phi]\ };\)\[IndentingNewLine]
  \(FourLength[{pe_, pz_, px_, py_}]\  := \ Sqrt[Max[pe^2\  - \ pz^2\  - \ px^2\  - \ py^2, 0.0]];\)\)







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



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



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



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



eq[attr_,n_][objList_]:=attr[objList]\[Equal]n;
gt[attr_,n_][objList_]:=attr[objList]>n;
lt[attr_,n_][objList_]:=attr[objList]<n;
ge[attr_,n_][objList_]:=attr[objList]\[GreaterEqual]n;
le[attr_,n_][objList_]:=attr[objList]\[LessEqual]n;



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



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



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





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

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



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





all[evt_]:=True;



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]]\[Equal]0,True,
        False]];

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



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



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;





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;





(*NiceTake[List_,n_]=
      Take[List,If[n>0,Min[Length[List],n],-Min[Length[List],-n]]];*)

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





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





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



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-","\[Mu]-","\[Tau]-","e+","\[Mu]+","\[Tau]+"},
          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]]]\[LessEqual]F1[l[[2]]];
      MatrixForm[
        Transpose[
          Join[{{"-","e-","\[Mu]-","\[Tau]-","e+","\[Mu]+","\[Tau]+"}},
            Transpose[
              Join[{{"e-","\[Mu]-","\[Tau]-","e+","\[Mu]+","\[Tau]+"}},
                
                Map[If[UTfnc[#],CriterionCount[TMPLIST,LeptonContent[#]],0]&,
                  Outer[List,LEPLIST,LEPLIST],{2}]]]]]]];



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



\!\(\(\(Flatten4[List_] := Flatten[List, Max[0, Depth[List] - 4]];\)\(\[IndentingNewLine]\)
  \)\[IndentingNewLine]
  \(Freq[evtList_, crit_, objsel_, plotfunc_, {min_, max_, nbins_}] := BinCounts[plotfunc /@ Flatten4[objsel /@ Select[evtList, crit]], {min, max, \(max - min\)\/nbins}];\)\[IndentingNewLine]
  \(Hist[freqcombo__, {min_, max_, nbins_}, HistogramOptions___] := Histogram[freqcombo /. Bin[args__] \[Rule] Freq[args, {min, max, nbins}], FrequencyData \[Rule] True, HistogramCategories \[Rule] Table[i, {i, min, max, \(max - min\)\/nbins}], HistogramOptions];\)\[IndentingNewLine]
  \(MakeFreq[freqcombo__, {min_, max_, nbins_}] := freqcombo /. Bin[args__] \[Rule] Freq[args, {min, max, nbins}];\)\)





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

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

\!\(\(NewObject[Type_, {p0_, pz_, px_, py_}, NewObjectOptions___] := Module[{mass, eta, phi, pt, forcemass}, \[IndentingNewLine]forcemass = \((\((ForceMass /. {NewObjectOptions})\) /. Options[NewObject])\); \[IndentingNewLine]phi = Mod[ArcTan[px, py], 2\ \[Pi]]; \[IndentingNewLine]pt = \@\(py\^2 + px\^2\); \[IndentingNewLine]eta = ArcSinh[pz/pt]; \[IndentingNewLine]mass = \@\(p0\^2 - pz\^2 - px\^2 - py\^2\); \[IndentingNewLine]If[\(! MatchQ[forcemass, False]\), {mass, pt, eta, phi} = ForceMassAdjust[{mass, pt, eta, phi}, forcemass]]; \[IndentingNewLine]\((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]\[IndentingNewLine]];\)\)

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

Options[ReconstructParent]=Options[NewObject];







RunInclusiveCounts[evtList_]:=
  Print["\nPhotons: ",SpeciesCount[evtList,o\[Gamma],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]]]













