Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Added Nose updated HGCAL calibration #27

Open
wants to merge 2 commits into
base: master
Choose a base branch
from
Open

Added Nose updated HGCAL calibration #27

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
721677a
updated calibration constants and added Nose calibration
selvaggi Oct 24, 2019
2abd265
added MC particle class
selvaggi Oct 24, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
23 changes: 20 additions & 3 deletions HGCalImagingAlgo.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@

sys.setrecursionlimit(100000)
# noise thresholds and MIPs
from RecHitCalibration import RecHitCalibration
from RecHitCalibration import RecHitCalibration, RecHitCalibrationNose

# definition of Hexel element

Expand Down Expand Up @@ -652,6 +652,8 @@ def recHitAboveThreshold(rHit, ecut, dependSensor=True, usePandas=False):
layer = rHit.layer()
thickness = rHit.thickness()
energy = rHit.energy()

#print(rHit.eta(), layer, thickness, energy, HGCalImagingAlgo.lastLayerFH)

if(dependSensor):
thickIndex = -1
Expand All @@ -660,14 +662,29 @@ def recHitAboveThreshold(rHit, ecut, dependSensor=True, usePandas=False):
if(thickness > 99. and thickness < 101.): thickIndex = 0
elif(thickness > 199. and thickness < 201.): thickIndex = 1
elif(thickness > 299. and thickness < 301.): thickIndex = 2
else: print("ERROR - silicon thickness has a nonsensical value")

## adding this case to avoid warnings for now
elif(thickness > 119. and thickness < 121.): thickIndex = 0

else:
print("ERROR - silicon thickness has a nonsensical value")
#print(rHit.eta(), layer, thickness, energy, HGCalImagingAlgo.lastLayerFH)

if abs(rHit.eta()) > 3.0 and (thickness < 119. or thickness > 121.):
print(rHit.eta(), layer, thickness, energy)

# determine noise for each sensor/subdetector using RecHitCalibration library

RecHitCalib = RecHitCalibration()
sigmaNoise = 0.001 * RecHitCalib.sigmaNoiseMeV(layer, thickIndex) # returns threshold for EE, FH, BH (in case of BH thickIndex does not play a role)

# treat Nose hits
if abs(rHit.z()) > 1000.0:
RecHitCalibNose = RecHitCalibrationNose()
sigmaNoise = 0.001 * RecHitCalibNose.sigmaNoiseMeV(layer, thickIndex)

aboveThreshold = energy >= ecut * sigmaNoise # this checks if rechit energy is above the threshold of ecut (times the sigma noise for the sensor, if that option is set)
return sigmaNoise, aboveThreshold


def getEnergy(item):
return item.energy
30 changes: 30 additions & 0 deletions NtupleDataFormat.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -204,6 +204,10 @@ def genParticles(self, prefix="genpart"):
"""Returns generator particles object."""
return GenParticles(self._tree, prefix)

def mcParticles(self, prefix="gen"):
"""Returns generator particles object."""
return McParticles(self._tree, prefix)

def primaryVertex(self, prefix="vtx"):
"""Returns PrimaryVertex object."""
return PrimaryVertex(self._tree, prefix)
Expand Down Expand Up @@ -360,6 +364,32 @@ def __init__(self, tree, prefix):
# self.prefix = prefix
super(GenParticles, self).__init__(tree, prefix + "_pt", GenParticle, prefix)

##########
class McParticle(_Object):
"""Class representing a McParticle."""

def __init__(self, tree, index, prefix):
"""Constructor.

Arguments:
tree -- TTree object
index -- Index of the McParticle
prefix -- TBranch prefix
"""
super(McParticle, self).__init__(tree, index, prefix)

class McParticles(_Collection):
"""Class presenting a collection of McParticles."""

def __init__(self, tree, prefix):
"""Constructor.

Arguments:
tree -- TTree object
prefix -- TBranch prefix
"""
# self.prefix = prefix
super(McParticles, self).__init__(tree, prefix + "_pt", McParticle, prefix)

##########
class LayerCluster(_Object):
Expand Down
158 changes: 104 additions & 54 deletions RecHitCalibration.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,62 +8,64 @@ class RecHitCalibration:
def __init__(self):
"""set variables used in the functions"""
# https://github.com/cms-sw/cmssw/blob/CMSSW_9_3_X/RecoLocalCalo/HGCalRecProducers/python/HGCalRecHit_cfi.py#L5
self.dEdX_weights = (0.0, # there is no layer zero
8.603, # Mev
8.0675,
8.0675,
8.0675,
8.0675,
8.0675,
8.0675,
8.0675,
8.0675,
8.9515,
10.135,
10.135,
10.135,
10.135,
10.135,
10.135,
10.135,
10.135,
10.135,
11.682,
13.654,
13.654,
13.654,
13.654,
13.654,
13.654,
13.654,
38.2005,
55.0265,
49.871,
49.871,
49.871,
49.871,
49.871,
49.871,
49.871,
49.871,
49.871,
49.871,
62.005,
83.1675,
92.196,
92.196,
92.196,
92.196,
92.196,
92.196,
92.196,
92.196,
92.196,
92.196,
46.098)
self.dEdX_weights = (0.0, # there is no layer zero
8.366557, # Mev
10.425456,
10.425456,
10.425456,
10.425456,
10.425456,
10.425456,
10.425456,
10.425456,
10.425456,
10.425456,
10.425456,
10.425456,
10.425456,
10.425456,
10.425456,
10.425456,
10.425456,
10.425456,
10.425456,
10.425456,
10.425456,
10.425456,
10.425456,
10.425456,
10.425456,
10.425456,
31.497849,
51.205434,
52.030486,
52.030486,
52.030486,
52.030486,
52.030486,
52.030486,
52.030486,
52.030486,
52.030486,
52.030486,
71.265149,
90.499812,
90.894274,
90.537470,
89.786205,
89.786205,
89.786205,
89.786205,
89.786205,
89.786205,
89.786205,
89.786205,
89.786205)

# https://github.com/cms-sw/cmssw/blob/CMSSW_9_3_X/RecoLocalCalo/HGCalRecProducers/python/HGCalRecHit_cfi.py#L86
self.thicknessCorrection = (1.132,1.092,1.084) # 100, 200, 300 um
self.thicknessCorrection = (0.759,0.760,0.773) #120um, 200um, 300um

#### # understand what is should be for 120 um

# Base configurations for HGCal digitizers
# https://github.com/cms-sw/cmssw/blob/CMSSW_9_3_X/SimCalorimetry/HGCalSimProducers/python/hgcalDigitizer_cfi.py#L5
Expand Down Expand Up @@ -96,3 +98,51 @@ def sigmaNoiseMIP(self, layer, thicknessIndex):

def sigmaNoiseMeV(self, layer, thicknessIndex):
return self.sigmaNoiseMIP(layer, thicknessIndex) * self.MeVperMIP(layer, thicknessIndex)


class RecHitCalibrationNose:

def __init__(self):
"""set variables used in the functions"""
# https://github.com/cms-sw/cmssw/blob/CMSSW_9_3_X/RecoLocalCalo/HGCalRecProducers/python/HGCalRecHit_cfi.py#L5
self.dEdX_weights = (0.0, # there is no layer zero
39.500245, # Mev
39.756638,
39.756638,
39.756638,
39.756638,
66.020266,
92.283895,
92.283895)

# https://github.com/cms-sw/cmssw/blob/CMSSW_9_3_X/RecoLocalCalo/HGCalRecProducers/python/HGCalRecHit_cfi.py#L86

### ALL sensors seem to have thickness = 120 um in NOSE
self.thicknessCorrection = 0.759 # understand what is should be for 120 um

# Base configurations for HGCal digitizers
# https://github.com/cms-sw/cmssw/blob/CMSSW_9_3_X/SimCalorimetry/HGCalSimProducers/python/hgcalDigitizer_cfi.py#L5
# self.eV_per_eh_pair = 3.62
self.fC_per_ele = 1.6020506e-4
self.nonAgedNoises = 2100.0 # 100,200,300 um (in electrons)

# https://github.com/cms-sw/cmssw/blob/CMSSW_9_3_X/RecoLocalCalo/HGCalRecProducers/python/HGCalUncalibRecHit_cfi.py#L25
self.fCPerMIP = 1.25 # 100um, 200um, 300um

# https://github.com/cms-sw/cmssw/blob/CMSSW_9_3_X/SimCalorimetry/HGCalSimProducers/python/hgcalDigitizer_cfi.py#L127
self.noise_MIP = 1.0/7.0 #expectation based on latest SiPM performance

def MeVperMIP(self, layer, thicknessIndex):
return self.dEdX_weights[layer]/self.thicknessCorrection

def MIPperGeV(self, layer, thicknessIndex):
return 1000./MeVperMIP(layer, thicknessIndex)

def sigmaNoiseMIP(self, layer, thicknessIndex):
return self.fC_per_ele * self.nonAgedNoises / self.fCPerMIP

def sigmaNoiseMeV(self, layer, thicknessIndex):
return self.sigmaNoiseMIP(layer, thicknessIndex) * self.MeVperMIP(layer, thicknessIndex)