diff --git a/HGCalImagingAlgo.py b/HGCalImagingAlgo.py index 63b32cfa..644927f7 100755 --- a/HGCalImagingAlgo.py +++ b/HGCalImagingAlgo.py @@ -16,7 +16,7 @@ sys.setrecursionlimit(100000) # noise thresholds and MIPs -from RecHitCalibration import RecHitCalibration +from RecHitCalibration import RecHitCalibration, RecHitCalibrationNose # definition of Hexel element @@ -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 @@ -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 diff --git a/NtupleDataFormat.py b/NtupleDataFormat.py index b83259fa..14df9a78 100755 --- a/NtupleDataFormat.py +++ b/NtupleDataFormat.py @@ -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) @@ -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): diff --git a/RecHitCalibration.py b/RecHitCalibration.py index 1772887c..ad4be4e3 100755 --- a/RecHitCalibration.py +++ b/RecHitCalibration.py @@ -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 @@ -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) + + +