Belle II Software  release-05-01-25
gainAnalysisPixelToPMT.py
1 #!/usr/bin/env python
2 # -*- coding: utf-8 -*-
3 
4 # ---------------------------------------------------------------------------------------
5 # convert an output root file of pixel-by-pixel gain analysis data to pmt-by-pmt format
6 # Contibutors: Maeda Yosuke (KMI, Nagoya Univ.)
7 #
8 # Usage: basf2 gainAnalysisPixelToPMT.py inputRoot(pixel-by-pixel result) outputRoot [maxHitTiming=-8.] [nMinEntries=100]
9 # ---------------------------------------------------------------------------------------
10 
11 import sys
12 import time
13 import numpy
14 from array import array
15 from basf2 import *
16 from ROOT import gROOT
17 gROOT.SetBatch(True)
18 from ROOT import TFile, TTree, TMath
19 
20 
21 def gainana_pixelToPmt(inputRoot, outputRoot, maxHitTiming=-8., nMinEntries=100):
22 
23  f = TFile(inputRoot)
24  tr = f.Get("tree")
25 
26  f_out = TFile(outputRoot, "RECREATE")
27  tr_out = TTree("tree_pmt", "gain/eff. analysis output in PMT-by-PMT")
28 
29  nPmtCh = 16
30  nTopModules = 16
31  nPmtPerModule = 32
32  slotId = array('i', [0])
33  pmtId = array('i', [0])
34  nValidCh = array('i', [0])
35  meanGain = array('f', [0.])
36  stdDevGain = array('f', [0.])
37  maxMinRatio = array('f', [0.])
38  gain = array('f', nPmtCh * [0.])
39  efficiency = array('f', nPmtCh * [0.])
40 
41  tr_out.Branch('slotId', slotId, 'slotId/I')
42  tr_out.Branch('pmtId', pmtId, 'pmtId/I')
43  tr_out.Branch('nValidCh', nValidCh, 'nValidCh/I')
44  tr_out.Branch('meanGain', meanGain, 'meanGain/F')
45  tr_out.Branch('stdDevGain', stdDevGain, 'stdDevGain/F')
46  tr_out.Branch('maxMinRatio', maxMinRatio, 'maxMinRatio/F')
47  tr_out.Branch('gain', gain, 'gain[' + str(nPmtCh) + ']/F')
48  tr_out.Branch('efficiency', efficiency, 'efficiency[' + str(nPmtCh) + ']/F')
49 
50  for iSlot in range(nTopModules):
51  for iPMT in range(nPmtPerModule):
52 
53  slotId[0] = iSlot + 1
54  pmtId[0] = iPMT + 1
55 
56  cut = "slotId==" + str(slotId[0]) + " && pmtId==" + str(pmtId[0]) + \
57  " && hitTiming<(" + str(maxHitTiming) + ") && nEntries>" + str(nMinEntries)
58  tr.Draw("gain:efficiency:pmtChId", cut)
59  nValidCh[0] = tr.GetEntries(cut)
60  if nValidCh[0] < 2:
61  continue
62  elif nValidCh[0] >= nPmtCh:
63  print("ERROR : too many channels for slot" + str(slotId[0]).zfill(2) +
64  " PMT" + str(pmtId[0]).zfill(2) + "(" + str(nValidCh[0]) + ")")
65 
66  gainArray = []
67  effArray = []
68  chArray = []
69  for iCh in range(nValidCh[0]):
70  gainArray.append(tr.GetV1()[iCh])
71  effArray.append(tr.GetV2()[iCh])
72  chArray.append(int(tr.GetV3()[iCh]))
73 
74  meanGain[0] = numpy.average(gainArray)
75  stdDevGain[0] = numpy.std(gainArray)
76  maxMinRatio[0] = max(gainArray) / min(gainArray)
77  print(cut + " : nValidCh = " + str(nValidCh[0]))
78  print(" --> mean = " + str('%05.2f' % meanGain[0]) + ", StdDev = " + str('%05.2f' % stdDevGain[0]) +
79  ", max/min = " + str('%03.2f' % maxMinRatio[0]))
80 
81  for iCh in range(nPmtCh):
82  for jCh in range(nValidCh[0]):
83  if chArray[jCh] - 1 is iCh:
84  gain[iCh] = gainArray[jCh]
85  efficiency[iCh] = effArray[jCh]
86  break
87  else:
88  gain[iCh] = -1
89  efficiency[iCh] = -1
90 
91  tr_out.Fill()
92 
93  f_out.cd()
94  tr_out.Write()
95  f_out.Close()
96 
97 if __name__ == '__main__':
98 
99  args = sys.argv
100  if len(args) > 4:
101  gainana_pixelToPmt(args[1], args[2], float(args[3]), int(args[4]))
102  elif len(args) > 3:
103  gainana_pixelToPmt(args[1], args[2], float(args[3]))
104  elif len(args) > 2:
105  gainana_pixelToPmt(args[1], args[2])
106  else:
107  print("usage:")
108  print(args[0] + " (root file for pixel-by-pixel gain result) (output file)" +
109  " [max. hit timing = -8 ns] [min. number of entries = 100]")