Introduction

The purpose of this twiki page is to gather useful links and information about a preliminary study of the decay H-->Zgamma. This is not intended to be in use long term.

Using the package

The package assumes that you are using smuhpc.smu.edu. Contact Aidan to setup an account. smuhpc (and smuhpc2) have access to all the data and MC, a condor cluster, and large amounts of disk space available.

Create a working directory and setup ATLAS: asetup 17.0.5.6,here

Then check out the package: svn co svn+ssh://@svn.cern.ch/reps/atlasusr/aidan/HToZGammaAnalysis

The README file contains information about getting started with the package.

To use condor add the following line to your .bashrc: source /grid/condor/condor.sh

Modules

The code is written in python with an interface for the C++ modules. The analysis is steered using main_analysis.py, where the task option determines what type of analysis should be performed. The other main modules are sample.py, event_seletion.py, pyrootmagic.py, style.py, trigger.py, truthmatch.py, and variable.py. In addition to these, each object has its own module.

main_analysis.py

The command line options for main_analysis.py are:

To add a new task to main_analysis.py simply go the end of the list of tasks and create a new block.

Reconstruction modes

To keep track of which reconstruction we use, we have the following modes:

sample.py

Each MC and data sample is handled using the sample object. This object is aware of the source of the (unskimmed) files, the number of events for the sample, the number of files, the cross section (for MC) and whether it is MC or data. The samples are specified in the sample.samples list. In main_analysis.py the user can pass a sample index (using the -s command line option) or sample name (using the -n command line option) to find a sample. main_analysis.py will save this information to the variables sample_index. In the case that a sample name is specified this is treated as a regex and the first sample name that matches is used.

In addition the samples are grouped into types and each type has its own style. The types are:

Each type has a function for applying style to a histogram: SetHistogramStyle( histogram )

There are separate constructors for MC and data samples:

class sample:
    def __init__(self, name, MCID, crossSection, nEvents, type, input_path):
class MCSample(sample):
    def __init__(self, name, MCID, crossSection, nEvents, type, MCTag, input_path):
class dataSample(sample):
    def __init__(self, name, runNumber, nEvents, type, input_path):

The stream of the sample is set to either eGamma, Muons or MC by matching the sample name, using the function setStream().

Each sample object has a TChain associated with it, which is accessed by sample.chain. By default this is set to 0 (null pointer). You can add the preskim files to the TChain by using the AddFiles() function. You can add the skimmed files using the AddSkimmedFiles() function (this assumes the skimmed files have been hadded successfully). When skimming the default output directory is output/skimmed.

Each MC sample has a luminosity associated with it (use the Lumi() function) and can be scaled to a given luminosity using the GetScale( lumi ) function, where lumi is given in pb.

To load a sample and obtain its skimmed TChain use the following:

s = sample.samples[sample_index]
s.AddSkimmedFiles()
ch = s.chain

pyrootmagic.py

pyrootmagic is a package developed by Emanuel Strauss (SLAC) that allows python to easily write TTrees. A list of variables that will be added in the skim is included in pyrootmagic.py. You can add a new variable using the following syntax:

After changing the properties of the list of variables the object must be rebuilt using the command ./main_analysis.py -b -t clean.

The pyrootmagic object variables are accessed directly using either pyrootmagic.var_1 = X or pyrootmagic.var_2.push_back( Y ), where the first syntax is for single valued variables and the second syntax is for a vector.

variable.py

The main variables are stored in variable.py. Each variable has a range, a binning, a unit, x- and y-axes (where the y-axis is automatically set to reflect the binning and units) and a ROOT expression for drawing. Using these parameters the variable object can also create a histogram. In addition to these basic parameters the variable object also has a mode (ee or mm), a number of decimal places for the unit (default=2), and excluded regions for placing cuts. There is also a 2D variable class to allow the plotting of one variable against another.

The syntax for creating a new variable is:

def __init__(self,name,branch_name,xaxis,units,mode):

eg: vars['e1_pt'] = variable.MakeVar('e1_pt', 'el_GSF_pt[e1_index[best_ee_index]]*1e-3', 'p_{T}(e_{1})', 'GeV', 200,    0,  200 , 'ee')

Cuts should be added using the exclusions functions. To add cuts, pass a list of excluded regions.

eg: vars['e1_pt'].AddExclusion(0,10) to exclude the region 0<e1pt<10GeV

You can define a list of variables using the var_collection object. The syntax for creating a var_collection object is:

collection = variable.var_collection([], 'ee', 'mainListOfVars_ee')

where the first argument is a list of variables, the second argument is the mode (either ee or mm) and the final argument is a (unique) name to make sure that histogram names do not clash.

This object has methods for creating a series of cutflow plots (in the order that the variables were added) and (n-1) plots.

For example:

collection = variable.var_collection([], 'ee', 'mainListOfVars_ee')
collection.add_var_with_exclusions(variable.vars['e1ph_DR'           ], [[0.0 , 0.7]])
collection.add_var_with_exclusions(variable.vars['e2ph_DR'           ], [[0.0 , 1.0]])
collection.add_var_with_exclusions(variable.vars['e1_pt'             ], [[0,40]])
collection.add_var_with_exclusions(variable.vars['e2_pt'             ], [[0,20]])
collection.add_var_with_exclusions(variable.vars['Zee_m'             ], [[0.0,91.2-10.0],[91.2+10.0,1e30]])
collection.add_var_with_exclusions(variable.vars['Hee_m'             ], [])

import style
legend = style.legend(0.6, 0.5, 0.8, 0.9)

legend.AddEntry( types['signal'     ].histogram, types['signal'     ].title, 'f')
legend.AddEntry( types['ZllGamma'   ].histogram, types['ZllGamma'   ].title, 'f')
legend.AddEntry( types['ttbar'      ].histogram, types['ttbar'      ].title, 'f')
legend.AddEntry( types['qcd'        ].histogram, types['qcd'        ].title, 'f')
legend.AddEntry( types['WZ'         ].histogram, types['WZ'         ].title, 'f')
legend.AddEntry( types['ZttJets'    ].histogram, types['ZttJets'    ].title, 'f')
legend.AddEntry( types['ZllJets'    ].histogram, types['ZllJets'    ].title, 'f')

histogram_signal = ROOT.TH1F('h_signal','',100,0,1)
histogram_signal.SetLineColor(ROOT.kRed)
histogram_signal.SetLineStyle(ROOT.kDashed)
legend.AddEntry(histogram_signal, 'Signal template', 'l')

for s in sample.samples:
    if 'data' in s.type.name:
        break
    collection.perform_cutflow        (s, 1.01e3, 'cutflow_test_', file_out)
collection.make_cutflow_plots  (file_out,legend)

for s in sample.samples:
    if 'data' in s.type.name:
        break
    collection.perform_nMinusOne_plots(s, lumi, 'cutflow_test_', file_out)
    collection.make_nMinusOne_plots(file_out, legend)
    collection.clean_nMinusOne()

truthmatchh.py

truthmatchh.py is used to truthmatch the candidates in an event. It defines the following variables:

When truthmatching the mode must be passed to the truthmatching module, where the muon blocks must be specified for muons.

truthmatch.py also finds the indices for the real leptons, Z, photon and Higgs, if they exist. (Indices are set to -1 otherwise.)

The module also contains some high level functions for navigating the MC block.

Particle modules

Each particle (electron, muon, photon, Z boson, Higgs boson) has its own module that manages that particle's properties. Each particle has a p4, and where the particle has children the children are also stored. Cuts relating to the particle should be stored in the relevant particle module. For example the Zboson module reads:

DeltaRCut = 0.2

class Zboson:
    def __init__(self,l1,l2):
      self.l1 = l1
      self.l2 = l2
      self.p4 = self.l1.p4 + self.l2.p4
    
    def pass_cuts(self, options):
        if options.nocut:
            return True
        if self.p4.M() < 10e3:
            return False
        if self.l1.p4.DeltaR( self.l2.p4) < DeltaRCut:
            return False
        return True

Each pass_cuts function must be aware of the variable options.nocut so that it can return True when the user requests that no cuts be applied.

The muon module contains more information than other modules because it must store the kinematic and isolation variables for ease of acces, and choose the correct block on a candidate by candidate basis.

You can define different kinds of selection criteria for the particles, depending on where you took recommendations from. You can change the definitions under the appropriate block, for exampe:

# Set definitions for different objects
electron_definition = 'HSG2'
muon_definition     = 'HSG2'
photon_definition   = 'SM_ZGamma'
jet_definition      = 'SM_ZGamma'
tau_definition      = 'tau_WG'

if nocut:
    electron_definition = 'No_Cut'
    muon_definition     = 'No_Cut'
    photon_definition   = 'No_Cut'
    jet_definition      = 'No_Cut'
    tau_definition      = 'No_Cut'

counter.py

The counter.py module contains the counter object that keeps track of the numbers of candidates on a per candidate and per event basis. Each counter has its own unique name each counter contains a histogram which can be filled on a per event basis.

At the start of a task you can initialize the counter with:

counters['electron_candidates'] = counter.counter('electron_candidates',  20) (This counter has a histogram going from 0 to 20, so it can store information about electron multiplicity up to 20.)

At the start of each event loop the counter should be initialized for that event counters['electron_candidates'].startEvent()

Increment the number of candidates with counters['electron_candidates'].incrementCount()

Increment the number of events with counters['electron_candidates'].incrementEvent()

Prevent further incrementing per event with counters['electron_candidates'].finishEvent()

Print the summary of the counter at the end of the analysis with counters['electron_candidates'].printSummary()

Example analysis

As an example consider the multiplicity of Z and Higgs boson candidates per event per sample. To add this task we can do the following:

elif task=='count_candidates':
    # Create a dictionary of histograms to store information
    histogram_Hee_candidates_MC = {}
    histogram_Hmm_candidates_MC = {}
    histogram_Hee_candidates_data = ROOT.TH1F('h_Hee_candidates_data', '', 20, -0.5, 19.5)
    histogram_Hmm_candidates_data = ROOT.TH1F('h_Hmm_candidates_data', '', 20, -0.5, 19.5)
    
    histogram_Hee_candidates_data.GetXaxis().SetTitle('Number of Higgs (ee#gamma) candidates')
    histogram_Hee_candidates_data.GetYaxis().SetTitle('Number of events')
            
    histogram_Hmm_candidates_data.GetXaxis().SetTitle('Number of Higgs (#mu#mu#gamma) candidates')
    histogram_Hmm_candidates_data.GetYaxis().SetTitle('Number of events')
    
    # Load and loop over the samples
    import sample as sample_module
    for s in sample_module.samples:
        s.AddSkimmedFiles()
        ch = s.chain
        
        # Create histogram if we're using MC
        # The stream is already set, so use that to see if we have MC or data
        if s.stream == 'MC':
            histogram_Hee_candidates_MC[s.name] = ROOT.TH1F('h_Hee_candidates_%s'%(s.name), '', 20, -0.5, 19.5)
            histogram_Hmm_candidates_MC[s.name] = ROOT.TH1F('h_Hmm_candidates_%s'%(s.name), '', 20, -0.5, 19.5)
            
            histogram_Hee_candidates_MC[s.name].GetXaxis().SetTitle('Number of Higgs (ee#gamma) candidates')
            histogram_Hee_candidates_MC[s.name].GetYaxis().SetTitle('Number of events')
            
            histogram_Hmm_candidates_MC[s.name].GetXaxis().SetTitle('Number of Higgs (#mu#mu#gamma) candidates')
            histogram_Hmm_candidates_MC[s.name].GetYaxis().SetTitle('Number of events')
        
        # Loop over the events
        for j in range(0,ch.GetEntries()):
            ch.GetEntry(j)
            
            # Since we are using pyroot we already have direct access to all the branches in the TTree:
            n_Hee = len(ch.Hee_m)
            n_Hmm = len(ch.Hmm_m)
            if s.stream == 'MC':
                histogram_Hee_candidates_MC[s.name].Fill(n_Hee)
                histogram_Hmm_candidates_MC[s.name].Fill(n_Hmm)
            else:
                histogram_Hee_candidates_data.Fill(n_Hee)
                histogram_Hmm_candidates_data.Fill(n_Hmm)
    
    # Now all the histograms have been filled with the correct values
    # Set the styles of the histograms:
    
    for s in sample_module.samples:
        s.type.SetHistogramStyle(histogram_Hee_candidates_MC[s.name])
        s.type.SetHistogramStyle(histogram_Hmm_candidates_MC[s.name])
        s.type.SetHistogramStyle(histogram_Hee_candidates_data)
        s.type.SetHistogramStyle(histogram_Hmm_candidates_data)
    
    # Create a legend
    import style
    legend = style.TLegend(05, 0.9, 0.8, 0.5) # Style of this TLegend is nice
    for t in sample_module.types:
        if 'data' in t:
            legend.AddEntry(sample_module.types[t].histogram,t,'p')
        else:
            legend.AddEntry(sample_module.types[t].histogram,t,'f')
    
    # Combine results.  Only add MC to the THStack
    stack_ee = ROOT.THStack()
    stack_mm = ROOT.THStack()
    for s in sample_module.samples:
        if 'data' in s.name:
            break
        stack_ee.Add(histogram_Hee_candidates_MC[s.name])
        stack_mm.Add(histogram_Hmm_candidates_MC[s.name])
    
    # Grab canvas from the style, already set to the recommended size
    c = style.canvas
    
    # Prints results to file
    histogram_Hee_candidates_MC.Draw('pe')
    stack_ee.Draw('same')
    histogram_Hee_candidates_MC.Draw('pe:same') # Redraw to get axes titles
    c.Print('plots/histogram_Hee_candidates.eps')
    
    histogram_Hmm_candidates_MC.Draw('pe')
    stack_mm.Draw('same')
    histogram_Hmm_candidates_MC.Draw('pe:same')
    c.Print('plots/histogram_Hmm_candidates.eps')
    
    # No need to garbage collect, python does that for us!

Then run this from the command line like this:

./main_analysis.py -b - count_candidates

Running the skim

When running the skim the output is by default written to output/skimmed and later hadded to output/total. To run the skim:

• Create the necessary directories
  • mkdir output
  • mkdir output/skimmed
  • mkdir output/total
• Create the pyrootmagic object
  • ./main_analysis.py -b -t clean
• Create condor submission scripts for the skim
  • ./main_analysis.py -b -t condor
• Submit the condor jobs (no more than 400 at a time, please!)
  • head -n 400 condor.sh > condor_1.sh
  • head -n 800 condor.sh | tail -n 400 > condor_2.sh etc
  • source condor_1.sh etc
• hadd the output to save time when loading files
  • source hadd.sh

This should allow you to create your own skim. Skimming takes about 5 hours for all the MC and the first fb^-1 from 2011, allowing you to perform all the necessary analyses and crosscheckeds with data and MC.

Analysis

Event selection

Trigger

For MC the trigger responses are luminosity weighted. See trigger.py for details.

GRL

Use GoodRunsLists-00-00-96 and data11_7TeV.periodAllYear_HEAD_CoolRunQuery-00-04-08_All_Good.xml.

Primary vertices

At least one primary vertex with at least 3 tracks.

LAr error

larError!=2

Objects

EITHER

• Two muons where:
  • Either mu_staco or mu_muid
    • If a mu_staco candidate is within ΔR of 0.05 of a mu_muid candidate, remove the mu_muid candidate.
  • Both pass loose requirements (release 17)
  • Muon pT>=7GeV

OR

• Two electrons where:
  • Both pass loose requirements (release 17)
  • Electron pT>=7GeV

AND

• One photon where:
  • Passes loose requirements

Muons

Hit requirements:

   Int_t n = (*mu_staco_nTRTHits)[mu_i] + (*mu_staco_nTRTOutliers)[mu_i];
   Double_t mu_eta = fabs((*mu_staco_eta)[mu_i]);
   Bool_t case1 = ((n > 5) && (*mu_staco_nTRTOutliers)[mu_i] < .9*n) ; //require n>5 && the = of trt outliers < .9n
   Bool_t case2 = (n > 5) ? ((*mu_staco_nTRTOutliers)[mu_i] < .9*n) : true ; //if n>5, we require the = of trt outliers < .9n
   Bool_t trt_ext = ( 0.1 < mu_eta < 1.9) ? case1 : case2 ; // if |eta| < 1.9, then case 1; if |eta|>=1.9 or |eta| <0.1 , then case 2 

See the HSG2 page for the latest recommendations.

Reconstruction

The decay chain H-->Zgamma, Z-->ll is reconstructed.

All possible candidates are considered.

In the skim Z candidates are required to have m(ll)>10GeV

MC reweighting

Cross section reweighting

Reference: (pdf) SM cross sections

Reference for diboson processes: ATLAS Z gamma etc cross sections

LHC Yellow Book for Higgs cross sections (CDS)

Pileup weighting

Pileup weighting is implemented using PhysicsAnalysis/AnalysisCommon/PileupReweighting PileupReweighting-00-02-05

Smearing and scaling

The following packages are used for smearing/scaling MC (to be implemented):

Tag	svn Path
ggFReweighting-00-00-08	svn+ssh://svn.cern.ch/reps/atlasoff/PhysicsAnalysis/HiggsPhys/ggFReweighting/tags/
MuonMomentumCorrections-00-06-08	svn+ssh://svn.cern.ch/reps/atlasoff/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonMomentumCorrections/tags/
MuonEfficiencyCorrections-02-01-01	svn+ssh://svn.cern.ch/reps/atlasoff/PhysicsAnalysis/MuonID/MuonIDAnalysis/MuonEfficiencyCorrections/tags/
MuonIsolationCorrection-00-08	svn+ssh://svn.cern.ch/reps/penn/schae/MuonIsolationCorrection/tags/
egammaAnalysisUtils-00-03-28	svn+ssh://svn.cern.ch/reps/atlasoff/Reconstruction/egamma/egammaAnalysis/egammaAnalysisUtils/tags/
TrigMuonEfficiency-00-02-05	svn+ssh://svn.cern.ch/reps/atlasoff/Trigger/TrigAnalysis/TrigMuonEfficiency/tags/
ZMassConstraint-00-00-06	svn+ssh://svn.cern.ch/reps/atlasphys/Physics/Higgs/HSG2/Code/ZMassConstraint/tags/

MC samples

Signal samples:

Background SM samples:

Cross sections can be found at HSG6 background MC page

The following samples are located on smuhpc:

To be added: single top WW ZZ

General links and references

Search for excited leptons with llgamma final state:

Pubic note

Internal note

MCTruthClassifier

Muon combined performance recommendations

egamma recommendations

H-->4l selection (Winter 2012)

MCTruthClassifier (SVNWeb)

-- AidanRandleConde - 09-Jan-2012