/data/refman/pasoursint/CMSSW_6_1_2_SLHC2/src/ElectroWeakAnalysis/Utilities/python/Averages_3Invpb.py

Go to the documentation of this file.

#!/usr/bin/env python
######################################################################################
## Program to average CMS W and Z cross sections for muons and electrons. 
##    Run it as: "python Averages_3Invpb.py"
##
## Notes:
##
## a) For 3 inverse pb statistics, all likelihood profiles are already 
##    Gaussian to a good approximation. Therefore naive combination procedures
##    work.
## b) OPTION allows to average electron and muon measurements in two ways:
##
##    OPTION = "StatisticalAverage" does the average according to statistical 
##    uncertainties. Pros: This average is more rigurous from a statistical point 
##    of view, since some likelihood ansatz is necessary to interpret systematics.
##    Cons: This procedure leads to larger "overall" uncertainties at the level of the 
##    final combination since, for similar statistical uncertainties, measurements 
##    with larger systematics weight the same as measurements with lower systematics.
##
##    OPTION = ""StatisticalPlusSystematicAverage" does the average according
##    to the overall uncertainties. It assumes that systematic ucnertainties 
##    can be treated in a naive Gaussian way and so added quadratically to 
##    statistical uncertainties in the usual way. Correlations are taken into 
##    account.  A covariancia matrix "V" is built, and the solution X corresponds 
##    to the minimization of the expression " sum_{ij}(X-x_i V_{ij}^{-1} (X-x_j)", 
##    where x_i are the electron and muon measurements. Pros: this leads to minimal
##    uncertainties for the overall uncertainty (if stat. and syst. are added in 
##    quadrature, as people usually do). Cons: most of the systematic soruces are 
##    not statistical in origin, so giving them a 68% CL Gaussian meaning is 
##    an ad-hoc assumption. 
##
######################################################################################

from __future__ import division
from math import *

OPTION = "StatisticalPlusSystematicAverage"
#OPTION= "StatisticalAverage"

print ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>"
print ">>>>> METHOD TO AVERAGE MUONS AND ELECTRONS is: '%s'" % (OPTION)
print ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>\n"

######################################################################################
###  INPUTS FOLLOW (they are not the final version)
######################################################################################

# Relative luminosity error
relSysLumi = 11e-2

# Electron inputs
Wenu = 10.1447
absStatWenu = 0.1123
absCorrWenu = Wenu*sqrt(0.008**2+0.0116**2) # theory uncertainty
absUncWenu = 0.5031
print "\nWenu cross section = %.4f +- %.4f (stat.) +- %.4f (syst.) +- %.4f (theo.) +- %.4f (lumi.) [nb]" % (Wenu, absStatWenu, absUncWenu, absCorrWenu, Wenu*relSysLumi)
print "  Systematics >>>>>>>"
print "\tUncorrelated with muons:           %.2f %%" % (100./Wenu*absUncWenu)
print "\tTheory:                            %.2f %%" % (100./Wenu*absCorrWenu)
print "\tTOTAL (LUMI excluded):             %.2f %%\n" % (100./Wenu*sqrt(absUncWenu**2+absCorrWenu**2))

Wplusenu = 5.9349
absStatWplusenu = 0.0793
absCorrWplusenu = Wplusenu*sqrt(0.009**2+0.0133**2) # theory uncertainty
absUncWplusenu = 0.3446
print "\nWplusenu cross section = %.4f +- %.4f (stat.) +- %.4f (syst.) +- %.4f (theo.) +- %.4f (lumi.) [nb]" % (Wplusenu, absStatWplusenu, absUncWplusenu, absCorrWplusenu, Wplusenu*relSysLumi)
print "  Systematics >>>>>>>"
print "\tUncorrelated with muons:           %.2f %%" % (100./Wplusenu*absUncWplusenu)
print "\tTheory:                            %.2f %%" % (100./Wplusenu*absCorrWplusenu)
print "\tTOTAL (LUMI excluded):             %.2f %%\n" % (100./Wplusenu*sqrt(absUncWplusenu**2+absCorrWplusenu**2))

Wminusenu = 4.1401
absStatWminusenu = 0.0690
absCorrWminusenu = Wminusenu*sqrt(0.015**2+0.0090**2) # theory uncertainty
absUncWminusenu = 0.2432
print "\nWminusenu cross section = %.4f +- %.4f (stat.) +- %.4f (syst.) +- %.4f (theo.) +- %.4f (lumi.) [nb]" % (Wminusenu, absStatWminusenu, absUncWminusenu, absCorrWminusenu, Wminusenu*relSysLumi)
print "  Systematics >>>>>>>"
print "\tUncorrelated with muons:           %.2f %%" % (100./Wminusenu*absUncWminusenu)
print "\tTheory:                            %.2f %%" % (100./Wminusenu*absCorrWminusenu)
print "\tTOTAL (LUMI excluded):             %.2f %%\n" % (100./Wminusenu*sqrt(absUncWminusenu**2+absCorrWminusenu**2))

Zee = 0.9595
absStatZee = 0.0370
absCorrZee = Zee*sqrt(0.011**2+0.0134**2) # theory uncertainty
absUncZee = 0.0565
print "\nZee cross section = %.4f +- %.4f (stat.) +- %.4f (syst.) +- %.4f (theo.) +- %.4f (lumi.) [nb]" % (Zee, absStatZee, absUncZee, absCorrZee, Zee*relSysLumi)
print "  Systematics >>>>>>>"
print "\tUncorrelated with muons:           %.2f %%" % (100./Zee*absUncZee)
print "\tTheory:                            %.2f %%" % (100./Zee*absCorrZee)
print "\tTOTAL (LUMI excluded):             %.2f %%\n" % (100./Zee*sqrt(absUncZee**2+absCorrZee**2))

Ratioenu = Wplusenu/Wminusenu
absStatRatioenu = 0.02883
absCorrRatioenu = Ratioenu*sqrt(0.017**2+0.0127**2) # theory uncertainty
absUncRatioenu = 0.076
print "\nRatioenu cross section = %.4f +- %.4f (stat.) +- %.4f (syst.) +- %.4f (theo.)" % (Ratioenu, absStatRatioenu, absUncRatioenu, absCorrRatioenu)
print "  Systematics >>>>>>>"
print "\tUncorrelated with muons:           %.2f %%" % (100./Ratioenu*absUncRatioenu)
print "\tTheory:                            %.2f %%" % (100./Ratioenu*absCorrRatioenu)
print "\tTOTAL:                             %.2f %%\n" % (100./Ratioenu*sqrt(absUncRatioenu**2+absCorrRatioenu**2))

WZe = Wenu/Zee
absStatWZe = 0.425
absCorrWZe = WZe*sqrt(0.009**2+0.0103**2) # theory uncertainty
absUncWZe = 0.524
print "\nWZe cross section = %.4f +- %.4f (stat.) +- %.4f (syst.) +- %.4f (theo.)" % (WZe, absStatWZe, absUncWZe, absCorrWZe)
print "  Systematics >>>>>>>"
print "\tUncorrelated with muons:           %.2f %%" % (100./WZe*absUncWZe)
print "\tTheory:                            %.2f %%" % (100./WZe*absCorrWZe)
print "\tTOTAL:                             %.2f %%\n" % (100./WZe*sqrt(absUncWZe**2+absCorrWZe**2))

# Muon inputs 
Wmunu = 9.922
absStatWmunu = 0.090
absCorrWmunu = Wmunu*sqrt(0.011**2+0.0136**2) # theory uncertainty
relUncFit = 1.3e-2
relUncPreTrig = 0.5e-2
relUncSysEff = 0.5e-2
relUncEff = sqrt(relUncFit**2+relUncPreTrig**2+relUncSysEff**2)
relUncMomRes = 0.3e-2
relUncRecoil = 0.4e-2
relUncMCStat = 1.4e-3/sqrt(2)
relUncBkg = sqrt(2.0e-2**2+0.2e-2**2)
absUncWmunu=Wmunu*sqrt(relUncEff**2+relUncMomRes**2+relUncRecoil**2+relUncMCStat**2+relUncBkg**2)
print "\nWmunu cross section = %.4f +- %.4f (stat.) +- %.4f (syst.) +- %.4f (theo.) +- %.4f (lumi.) [nb]" % (Wmunu, absStatWmunu, absUncWmunu, absCorrWmunu, Wmunu*relSysLumi)
print "  Systematics >>>>>>>"
print "\tEfficiency(Zfit,Pretrig,Z->W):     %.2f %%" % (100*relUncEff)
print "\tMomentum scale/resolution:         %.2f %%" % (100*relUncMomRes)
print "\tBackground subtraction:            %.2f %%" % (100*relUncBkg)
print "\tSignal Recoil modeling:            %.2f %%" % (100*relUncRecoil)
print "\tMC statistics (acceptance):        %.2f %%" % (100*relUncMCStat)
print "\tTheory:                            %.2f %%" % (100./Wmunu*absCorrWmunu)
print "\tTOTAL (LUMI excluded):             %.2f %%\n" % (100./Wmunu*sqrt(absUncWmunu**2+absCorrWmunu**2))

Wplusmunu = 5.844
absStatWplusmunu = 0.069
absCorrWplusmunu = Wplusmunu*sqrt(0.013**2+0.0142**2) # theory uncertainty
relUncFit = 1.3e-2
relUncPreTrig = 0.5e-2
relUncSysEff = 0.5e-2
relUncEff = sqrt(relUncFit**2+relUncPreTrig**2+relUncSysEff**2)
relUncMomRes = 0.3e-2
relUncRecoil = 0.4e-2
relUncMCStat = 1.4e-3
relUncBkg = sqrt(1.7e-2**2+0.2e-2**2)
absUncWplusmunu=Wplusmunu*sqrt(relUncEff**2+relUncMomRes**2+relUncRecoil**2+relUncMCStat**2+relUncBkg**2)
print "\nWplusmunu cross section = %.4f +- %.4f (stat.) +- %.4f (syst.) +- %.4f (theo.) +- %.4f (lumi.) [nb]" % (Wplusmunu, absStatWplusmunu, absUncWplusmunu, absCorrWplusmunu, Wplusmunu*relSysLumi)
print "  Systematics >>>>>>>"
print "\tEfficiency(Zfit,Pretrig,Z->W):     %.2f %%" % (100*relUncEff)
print "\tMomentum scale/resolution:         %.2f %%" % (100*relUncMomRes)
print "\tBackground subtraction:            %.2f %%" % (100*relUncBkg)
print "\tSignal Recoil modeling:            %.2f %%" % (100*relUncRecoil)
print "\tMC statistics (acceptance):        %.2f %%" % (100*relUncMCStat)
print "\tTheory:                            %.2f %%" % (100./Wplusmunu*absCorrWplusmunu)
print "\tTOTAL (LUMI excluded):             %.2f %%\n" % (100./Wplusmunu*sqrt(absUncWplusmunu**2+absCorrWplusmunu**2))

Wminusmunu = 4.078
absStatWminusmunu = 0.057
absCorrWminusmunu = Wminusmunu*sqrt(0.019**2+0.0126**2) # theory uncertainty
relUncFit = 1.3e-2
relUncPreTrig = 0.5e-2
relUncSysEff = 0.5e-2
relUncEff = sqrt(relUncFit**2+relUncPreTrig**2+relUncSysEff**2)
relUncMomRes = 0.3e-2
relUncRecoil = 0.4e-2
relUncMCStat = 1.4e-3
relUncBkg = sqrt(2.3e-2**2+0.2e-2**2)
absUncWminusmunu=Wminusmunu*sqrt(relUncEff**2+relUncMomRes**2+relUncRecoil**2+relUncMCStat**2+relUncBkg**2)
print "\nWminusmunu cross section = %.4f +- %.4f (stat.) +- %.4f (syst.) +- %.4f (theo.) +- %.4f (lumi.) [nb]" % (Wminusmunu, absStatWminusmunu, absUncWminusmunu, absCorrWminusmunu, Wminusmunu*relSysLumi)
print "  Systematics >>>>>>>"
print "\tEfficiency(Zfit,Pretrig,Z->W):     %.2f %%" % (100*relUncEff)
print "\tMomentum scale/resolution:         %.2f %%" % (100*relUncMomRes)
print "\tBackground subtraction:            %.2f %%" % (100*relUncBkg)
print "\tSignal Recoil modeling:            %.2f %%" % (100*relUncRecoil)
print "\tMC statistics (acceptance):        %.2f %%" % (100*relUncMCStat)
print "\tTheory:                            %.2f %%" % (100./Wminusmunu*absCorrWminusmunu)
print "\tTOTAL (LUMI excluded):             %.2f %%\n" % (100./Wminusmunu*sqrt(absUncWminusmunu**2+absCorrWminusmunu**2))

Zmumu = 0.924 # 0.893*1.025*1.01
absStatZmumu = 0.031 # 0.030*1.025*1.01
absCorrZmumu = Zmumu*sqrt(0.012**2+0.0158**2) # theory uncertainty
relUncEff = 0.5e-2 # pre-triggering
relUncFit= 1.0e-2
relUncMomRes = 0.2e-2
relUncMCStat = 0.4e-2
relUncBkg = sqrt(relUncFit**2+0.2e-2**2)
absUncZmumu=Zmumu*sqrt(relUncEff**2+relUncMomRes**2+relUncBkg**2+relUncMCStat**2)
print "\nZmumu cross section = %.4f +- %.4f (stat.) +- %.4f (syst.) +- %.4f (theo.) +- %.4f (lumi.) [nb]" % (Zmumu, absStatZmumu, absUncZmumu, absCorrZmumu, Zmumu*relSysLumi)
print "  Systematics >>>>>>>"
print "\tEfficiency(Pretriggering):         %.2f %%" % (100*relUncEff)
print "\tMomentum scale/resolution:         %.2f %%" % (100*relUncMomRes)
print "\tBackground subtraction/fit:        %.2f %%" % (100*relUncBkg)
print "\tMC statistics (acceptance):        %.2f %%" % (100*relUncMCStat)
print "\tTheory:                            %.2f %%" % (100./Zmumu*absCorrZmumu)
print "\tTOTAL (LUMI excluded):             %.2f %%\n" % (100./Zmumu*sqrt(absUncZmumu**2+absCorrZmumu**2))

Ratiomunu = Wplusmunu/Wminusmunu
absStatRatiomunu = 0.026
absCorrRatiomunu = Ratiomunu*sqrt(0.021**2+0.0119**2) # theory uncertainty
relUncEff = 2.8e-2
relUncMomRes = 0.3e-2
relUncMCStat = sqrt(2)*1.4e-3
relUncBkg = 0.7e-2
absUncRatiomunu = Ratiomunu*sqrt(relUncEff**2+relUncMomRes**2+relUncMCStat**2+relUncBkg**2)
print "\nRatiomunu cross section = %.4f +- %.4f (stat.) +- %.4f (syst.) +- %.4f (theo.)" % (Ratiomunu, absStatRatiomunu, absUncRatiomunu, absCorrRatiomunu)
print "  Systematics >>>>>>>"
print "\tEfficiency(W+ versus W- tests):    %.2f %%" % (100*relUncEff)
print "\tMomentum scale/resolution:         %.2f %%" % (100*relUncMomRes)
print "\tBackground subtraction:            %.2f %%" % (100*relUncBkg)
print "\tMC statistics (acceptance):        %.2f %%" % (100*relUncMCStat)
print "\tTheory:                            %.2f %%" % (100./Ratiomunu*absCorrRatiomunu)
print "\tTOTAL:                             %.2f %%\n" % (100./Ratiomunu*sqrt(absUncRatiomunu**2+absCorrRatiomunu**2))

WZmu = Wmunu/Zmumu
absStatWZmu = 0.3675
absCorrWZmu = WZmu*sqrt(0.011**2+0.0135**2) # theory uncertainty
relUncEffW = 1.3e-2
relUncMomResW = 0.3e-2
relUncMomResZ = 0.2e-2
relSysSubtract = sqrt(relUncEffW**2+relUncMomResW**2+relUncMomResZ**2)
relSysAdd = abs(relUncMomResW-relUncMomResZ)
absUncWZmu = WZmu*sqrt((absUncWmunu/Wmunu)**2 + (absUncZmumu/Zmumu)**2 - relSysSubtract**2 + relSysAdd**2)
print "\nWZmu cross section = %.4f +- %.4f (stat.) +- %.4f (syst.) +- %.4f (theo.)" % (WZmu, absStatWZmu, absUncWZmu, absCorrWZmu)
print "  STATISTICAL UNCERTAINTY INCLUDES EFFICIENCY (via Zmumu fit)"
print "  Systematics >>>>>>>"
print "\tUncorrelated with electrons:       %.2f %%" % (100./WZmu*absUncWZmu)
print "\tTheory:                            %.2f %%" % (100./WZmu*absCorrWZmu)
print "\tTOTAL:                             %.2f %%\n" % (100./WZmu*sqrt(absUncWZmu**2+absCorrWZmu**2))

######################################################################################
###  Utility functions
######################################################################################

###
def vbtfXSectionCheck(title, xsection, exsection, sysUnc, sysCor, relSysLumi):
      absSysLumi = xsection*relSysLumi
      print "VBTF inputs: %s = %.4f +- %.4f (stat.) +- %.4f (exp.) +- %.4f (the.) +- %.4f (lumi) [nb]" % (title, xsection, exsection, sysUnc, sysCor, absSysLumi) 

###
def vbtfXSectionAverage(title, xsection1, xsection2, exsection1, exsection2, sysUnc1, sysUnc2, sysCor1, sysCor2, relSysLumi):
      if OPTION== "StatisticalAverage":
            V11 = exsection1**2
            V22 = exsection2**2
            V12 = 0
      else:
            V11 = exsection1**2+sysUnc1**2+sysCor1**2
            V22 = exsection2**2+sysUnc2**2+sysCor2**2
            V12 = sysCor1*sysCor2

      a1 = (V22-V12)/(V11+V22-2*V12)
      a2 = (V11-V12)/(V11+V22-2*V12)
      average = a1*xsection1 + a2*xsection2
      errstat = sqrt(a1**2*exsection1**2+a2**2*exsection2**2)
      errunco = sqrt(a1**2*sysUnc1**2 + a2**2*sysUnc2**2)
      errtheo = sqrt(a1**2*sysCor1**2 + a2**2*sysCor2**2 + 2*a1*a2*sysCor1*sysCor2)
      errsyst = sqrt(errunco**2+errtheo**2)

      print "VBTF average: %s = %.4f +- %.4f (stat.) [nb]" % (title, average, errstat) 

      absSysLumi = average*relSysLumi
      print "\tVBTF systematics (1): +- %.4f (exp) +- %.4f (the) +- %.4f (lumi) [nb]" % (errunco, errtheo, absSysLumi) 
      print "\tVBTF systematics (2): +- %.4f (exp+the) +- %.4f (lumi) [nb]" % (errsyst, absSysLumi) 

###
def vbtfRatioCheck(title, ratio, absStat, sysUnc, sysCor):
      print "VBTF inputs: %s = %.4f +- %.4f (stat.) +- %.4f (exp.) +- %.4f (the.)" % (title, ratio, absStat, sysUnc, sysCor) 

###
def vbtfRatioAverage(title, ratio1, ratio2, eratio1, eratio2, sysUnc1, sysUnc2, sysCor1, sysCor2):
      if OPTION== "StatisticalAverage":
            V11 = eratio1**2
            V22 = eratio2**2
            V12 = 0
      else:
            V11 = eratio1**2+(sysUnc1**2+sysCor1**2)
            V22 = eratio2**2+(sysUnc2**2+sysCor2**2)
            V12 = sysCor1*sysCor2

      a1 = (V22-V12)/(V11+V22-2*V12)
      a2 = (V11-V12)/(V11+V22-2*V12)
      average = a1*ratio1 + a2*ratio2
      errstat = sqrt(a1**2*eratio1**2+a2**2*eratio2**2)
      errunco = sqrt(a1**2*sysUnc1**2 + a2**2*sysUnc2**2)
      errtheo = sqrt(a1**2*sysCor1**2 + a2**2*sysCor2**2 + 2*a1*a2*sysCor1*sysCor2)
      errsyst = sqrt(errunco**2+errtheo**2)

      print "VBTF average: %s = %.4f +- %.4f (stat.)" % (title, average, errstat) 

      print "\tVBTF systematics (1): +- %.4f (exp) +- %.4f (the)" % (errunco, errtheo) 
      print "\tVBTF systematics (2): +- %.4f (exp+the)" % (errsyst) 

######################################################################################
###  MAIN CALLS ...
######################################################################################

#############################################################
########## Wlnu total cross section
#############################################################
print "\n>>>>>>>>>>>>>>>"
vbtfXSectionCheck("W -> munu cross section",Wmunu,absStatWmunu,absUncWmunu,absCorrWmunu,relSysLumi)
vbtfXSectionCheck("W -> enu cross section",Wenu,absStatWenu,absUncWenu,absCorrWenu,relSysLumi)
vbtfXSectionAverage("W -> lnu cross section",Wmunu,Wenu,absStatWmunu,absStatWenu,absUncWmunu,absUncWenu,absCorrWmunu,absCorrWenu,relSysLumi)

#############################################################
########## Wplus -> lnu cross section
#############################################################
print "\n>>>>>>>>>>>>>>>"
vbtfXSectionCheck("W+ -> munu cross section",Wplusmunu,absStatWplusmunu,absUncWplusmunu,absCorrWplusmunu,relSysLumi)
vbtfXSectionCheck("W+ -> enu cross section",Wplusenu,absStatWplusenu,absUncWplusenu,absCorrWplusenu,relSysLumi)
vbtfXSectionAverage("W+ -> lnu cross section",Wplusmunu,Wplusenu,absStatWplusmunu,absStatWplusenu,absUncWplusmunu,absUncWplusenu,absCorrWplusmunu,absCorrWplusenu,relSysLumi)

#############################################################
########## Wminus -> lnu cross section
#############################################################
print "\n>>>>>>>>>>>>>>>"
vbtfXSectionCheck("W- -> munu cross section",Wminusmunu,absStatWminusmunu,absUncWminusmunu,absCorrWminusmunu,relSysLumi)
vbtfXSectionCheck("W- -> enu cross section",Wminusenu,absStatWminusenu,absUncWminusenu,absCorrWminusenu,relSysLumi)
vbtfXSectionAverage("W- -> lnu cross section",Wminusmunu,Wminusenu,absStatWminusmunu,absStatWminusenu,absUncWminusmunu,absUncWminusenu,absCorrWminusmunu,absCorrWminusenu,relSysLumi)

#############################################################
########## W+/W- ratio
#############################################################
#
print "\n>>>>>>>>>>>>>>>"
vbtfRatioCheck("W+ / W- cross section ratio, muon channel",Ratiomunu,absStatRatiomunu,absUncRatiomunu,absCorrRatiomunu)
vbtfRatioCheck("W+ / W- cross section ratio, electron channel",Ratioenu,absStatRatioenu,absUncRatioenu,absCorrRatioenu)
vbtfRatioAverage("W+ / W- cross section ratio",Ratiomunu,Ratioenu,absStatRatiomunu,absStatRatioenu,absUncRatiomunu,absUncRatioenu,absCorrRatiomunu,absCorrRatioenu)

#############################################################
########## Z > ll cross section (in 60 < Mll < 120 GeV)
#############################################################
#
print "\n>>>>>>>>>>>>>>>"
vbtfXSectionCheck("Z -> mumu cross section",Zmumu,absStatZmumu,absUncZmumu,absCorrZmumu,relSysLumi)
vbtfXSectionCheck("Z -> ee cross section",Zee,absStatZee,absUncZee,absCorrZee,relSysLumi)
vbtfXSectionAverage("Z -> ll cross section",Zmumu,Zee,absStatZmumu,absStatZee,absUncZmumu,absUncZee,absCorrZmumu,absCorrZee,relSysLumi)

#############################################################
########## W/Z ratio
#############################################################
#
print "\n>>>>>>>>>>>>>>>"
vbtfRatioCheck("W/Z ratio muons",WZmu,absStatWZmu,absUncWZmu,absCorrWZmu)
vbtfRatioCheck("W/Z ratio electrons",WZe,absStatWZe,absUncWZe,absCorrWZe)
vbtfRatioAverage("W/Z ratio",WZmu,WZe,absStatWZmu,absStatWZe,absUncWZmu,absUncWZe,absCorrWZmu,absCorrWZe)