laserPixelIlumination.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
 
 
 
# ------------------------------------------------------------------------
# example of making histograms of pixel hits coming from individual fibers
# needs as input the file produced by top/examples/simLaserCalibSystem.py
# ------------------------------------------------------------------------
 
import basf2 as b2
import os
import sys
from ROOT import Belle2
from ROOT import TH2F, TFile
 
inputFile = 'opticalGun.root'
if not os.path.exists(inputFile):
    b2.B2ERROR(inputFile + ': file not found')
    b2.B2INFO('File can be generated with top/examples/simLaserCalibSystem.py')
    sys.exit()
 
 
class Histogrammer(b2.Module):
 
    ''' A module to histogram pixel hits from individual fibers'''
 
    
    hist = [TH2F('fiber' + str(k + 1), 'Pixel hit distribution from fiber No.' + str(k + 1),
                 64, 0.5, 64.5, 8, 0.5, 8.5) for k in range(9)]
    
    pixelID = 482
    
    propTime = TH2F('propTime', 'Photon propagation times for pixel ' + str(pixelID) +
                    ' vs. fiber number', 9, 0.5, 9.5, 200, 0.0, 1.0)
    propTime.GetXaxis().SetTitle('fiber number')
    propTime.GetYaxis().SetTitle('propagation time [ns]')
 
    def event(self):
        ''' Event processor: fill histograms '''
 
        digits = Belle2.PyStoreArray('TOPDigits')
        for digit in digits:
            simhits = digit.getRelationsWith('TOPSimHits')  # several simhits can contribute
            for simhit in simhits:
                photon = simhit.getRelated('TOPSimPhotons')  # one or none possible
                if photon:
                    k = int((photon.getEmissionPoint().X() + 22.5) / 5.5)
                    if k >= 0 and k < 9:
                        self.histhist[k].Fill(digit.getPixelCol(), digit.getPixelRow())
                    else:
                        b2.B2ERROR('wrong decoding of fiber number: ' + str(k + 1))
                    if digit.getPixelID() == self.pixelIDpixelID:
                        t = photon.getDetectionTime() - photon.getEmissionTime()
                        self.propTimepropTime.Fill(k + 1, t)
 
    def terminate(self):
        ''' Write histograms to file '''
 
        tfile = TFile('laserPixelIlumination.root', 'recreate')
        for k in range(9):
            self.histhist[k].Write()
        self.propTimepropTime.Write()
        tfile.Close()
 
 
# Create path
main = b2.create_path()
 
# Input
roinput = b2.register_module('RootInput')
roinput.param('inputFileName', inputFile)
main.add_module(roinput)
 
# Histogrammer
main.add_module(Histogrammer())
 
# Show progress of processing
progress = b2.register_module('Progress')
main.add_module(progress)
 
# Process events
b2.process(main)
 
# Print call statistics
print(b2.statistics)