SND@LHC Software
Loading...
Searching...
No Matches
FillingScheme.py
Go to the documentation of this file.
1import ROOT,os
2import rootUtils as ut
3from urllib.request import urlopen
4import numpy
5import time,calendar
6import subprocess
7from XRootD import client
8import pickle
9from rootpyPickler import Pickler
10from rootpyPickler import Unpickler
11import atexit
12import requests
13import json
14import re
15
16def pyExit():
17 print("Make suicide until solution found for freezing")
18 os.system('kill '+str(os.getpid()))
19atexit.register(pyExit)
20
21fromElog = {4361:7902,4362: 7921, 4363: 7921, 4364: 7921, 4365: 7921, 4366: 7922, 4367: 7923, 4368: 7923, 4369: 7923, 4370: 7924, 4371: 7924, 4372: 7925, 4373: 7926, 4374: 7927, 4375: 7928, 4376: 7929, 4377: 7930, 4378: 7931, 4379: 7931, 4380: 7932, 4381: 7933, 4382: 7934, 4383: 7935, 4384: 7936, 4385: 7937, 4386: 7938, 4387: 7939, 4388: 7940, 4389: 7941, 4390: 7942, 4391: 7943, 4392: 7944, 4393: 7945, 4394: 7946, 4395: 7947, 4396: 7948, 4397: 7949, 4398: 7950, 4399: 7951, 4400: 7952, 4401: 7953, 4402: 7954, 4403: 7955, 4404: 7956, 4405: 7957, 4406: 7958, 4407: 7959, 4408: 7959, 4409: 7960, 4410: 7960, 4411: 7961, 4412: 7961, 4413: 7962, 4414: 7963, 4415: 7963, 4416: 7964, 4418: 7964, 4419: 7965, 4420: 7966, 4421: 7966, 4422: 7967, 4423: 7967, 4424: 7967, 4425: 7967, 4426: 7968, 4427: 7969, 4428: 7969, 4429: 7970, 4430: 7971, 4431: 7971, 4432: 7972, 4433: 7973, 4434: 7974, 4435: 7974, 4436: 7975, 4437: 7976, 4438: 7976, 4440: 7977, 4441: 7977, 4443: 7977, 4446: 7977, 4447: 7977, 4448: 7977, 4449: 7978, 4450: 7979, 4451: 7979, 4452: 7979, 4461: 7979, 4462: 7982, 4463: 7983, 4464: 7984, 4465: 7985, 4466: 7986, 4467: 7987, 4468: 7988, 4469: 7989, 4470: 7990, 4471: 7991, 4472: 7992, 4473: 7993, 4474: 7994, 4475: 7995, 4476: 7996, 4477: 7997, 4478: 7998, 4479: 7999, 4480: 7999, 4481: 8000, 4482: 8001, 4483: 8002, 4484: 8003, 4485: 8004, 4486: 8005, 4487: 8006, 4488: 8007, 4489: 8008, 4493: 8008, 4494: 8008, 4495: 8009, 4496: 8010, 4497: 8010, 4498: 8011, 4499: 8012, 4500: 8013, 4501: 8014, 4503: 8016, 4504: 8017, 4505: 8018, 4506: 8018, 4507: 8019, 4508: 8019, 4509: 8020, 4510: 8021, 4511: 8021, 4512: 8022, 4513: 8022, 4514: 8023, 4515: 8023, 4516: 8024, 4523: 8025, 4524: 8025, 4525: 8025, 4526: 8026, 4527: 8027, 4528: 8028, 4529: 8028, 4530: 8029, 4531: 8029, 4532: 8030, 4533: 8031, 4534: 8032, 4535: 8032, 4536: 8033, 4537: 8033, 4539: 8033, 4540: 8033, 4541: 8033, 4542: 8034, 4543: 8034, 4544: 8034, 4557: 8035, 4558: 8035, 4559: 8036, 4561: 8038, 4562: 8039, 4563: 8040, 4564: 8041, 4568: 8043, 4569: 8044, 4570: 8044, 4571: 8045, 4572: 8046, 4573: 8047, 4574: 8047, 4575: 8047, 4578: 8047, 4579: 8047, 4580: 8050, 4581: 8051, 4582: 8052, 4583: 8052, 4585: 8053, 4586: 8053, 4587: 8054, 4588: 8055, 4589: 8056, 4590: 8056, 4591: 8057, 4592: 8058, 4593: 8058, 4594: 8059, 4595: 8059, 4598: 8061, 4601: 8061, 4602: 8061, 4603: 8061, 4604: 8062, 4606: 8063, 4612: 8063, 4613: 8064, 4614: 8064, 4615: 8065, 4616: 8066, 4617: 8067, 4619: 8068, 4620: 8069, 4621: 8069, 4622: 8070, 4623: 8071, 4624: 8071, 4625: 8072, 4626: 8072, 4627: 8073, 4628: 8073, 4629: 8073, 4630: 8073, 4631: 8073, 4632: 8073, 4635: 8073, 4636: 8074, 4637: 8075, 4638: 8076, 4639: 8076, 4641: 8076, 4642: 8076, 4643: 8076, 4644: 8076, 4645: 8076, 4646: 8076, 4647: 8076, 4648: 8077, 4649: 8078, 4650: 8079, 4653: 8080, 4654: 8081, 4657: 8082, 4658: 8082, 4659: 8082, 4660: 8082, 4662: 8083, 4665: 8083, 4666: 8083, 4667: 8083, 4669: 8084, 4670: 8084, 4672: 8084, 4673: 8084, 4675: 8084, 4676: 8084, 4677: 8084, 4678: 8084, 4679: 8084, 4680: 8084, 4681: 8084, 4682: 8084, 4684: 8084, 4686: 8084, 4687: 8084, 4688: 8084, 4689: 8084, 4690: 8084, 4693: 8084, 4560: 8037, 4661: 8083}
22
23# other remarks
24# run 4986 fill 8234 Acc Message : - 300b fill for VELO insertion - issue with MKI - B2, we ramp with 218b
25
26
27emulsionReplacements = {0:1,4573:2,4859:3,5172:4,5485:5,6446:6} # first runs with new emulsion
28
30
31 def Init(self,options):
32 self.options = options
33 self.h={}
34 self.path = options.path
35 self.content = ''
36 self.phaseShift1 = 3564-404 # IP1
37 self.phaseShift2 = 129 # IP2
38 self.lpcFillingscheme = False
39 self.FSdict = {}
40 self.LumiInt = {}
41 self.runInfo = {}
42 self.beamCurrent = {}
43 if "FSdict.root" in os.listdir(options.path):
44 fg = ROOT.TFile.Open(options.path+'FSdict.root')
45 pkl = Unpickler(fg)
46 self.FSdict = pkl.load('FSdict')
47 fg.Close()
48 if "RunInfodict.root" in os.listdir(options.path):
49 fg = ROOT.TFile.Open(options.path+"RunInfodict.root")
50 pkl = Unpickler(fg)
51 self.runInfo = pkl.load('runInfo')
52 fg.Close()
53 if "Lumidict.root" in os.listdir(options.path):
54 fg = ROOT.TFile.Open(options.path+"Lumidict.root")
55 pkl = Unpickler(fg)
56 self.LumiInt = pkl.load('LumiInt')
57 fg.Close()
58
59 self.startTimes = { }
60 self.date = {}
61
62 def myPrint(self,tname,oname):
63 for t in ['.root','.pdf','.png']:
64 self.h[tname].Print(self.options.path+oname+t)
65
66 def readStartTimes(self):
67 months = {'Jan.':'01-','Feb.':'02-','Mar.':'03-','Apr.':'04-','May.':'05-','Jun.':'06-','Jul.':'07-','Aug.':'08-','Sep-':'09-','Oct.':'10-','Nov.':'11-','Dec.':'12-'}
68 s = open(self.path+'/StartTimes.log')
69 L = s.readlines()
70 for i in range(len(L)):
71 l = L[i]
72 if not l.find("(Info,Run)")>0: continue
73 offset = 0
74 if l.find('PM')>0: offset = 12*3600
75 offset -= 2*3600 # stupid time zone
76 for x in ['AM','PM']:
77 l = l.replace(x,'')
78 tmp = l.split("(Info,Run)")
79 tmp[0] = tmp[0].replace(' ','')
80 for m in months: tmp[0] = tmp[0].replace(m,months[m])
81 time_obj = time.strptime(tmp[0], '%m-%d,%Y-%H:%M:%S')
82 r = tmp[1].rfind('Run')
83 if r > 4570: continue #not needed, extracted from json file in run directory
84 run = int(tmp[1][r+3:].split('Started')[0])
85 self.startTimes[run] = calendar.timegm(time_obj)+offset
86
88 # starts with 4362 and ends with 4693
89 # from 4663 - 4687 are not valid
90 e = open(self.path+'/ELOG DAQ Test.html')
91
92 self.fromElogX = { }
93 L = e.readlines()
94 for i in range(len(L)):
95 l = L[i]
96 k = l.find('- New Fill')
97 if k<0: continue
98 fillNr = int(l[k+10:].split(' ')[1])
99 for j in range(i+1,len(L)):
100 k = L[j].find('Started with')
101 if k>0:
102 m = L[j][:k].rfind('Run')
103 runNr = int(L[j][m+3:k])
104 self.fromElogX[runNr] = fillNr
105 if L[j].find('- New Fill')>0: break
106
107# by hand manipulations:
108# 4560, only stops no start
109# 4661, only stops no start
110 self.fromElogX[4560]=8037
111 self.fromElogX[4661]=8083
112
113 def getNameOfFillingscheme(self,fillnr):
114 tmp = fillnr
115 alternative = self.alternativeFill(str(fillnr))
116 if alternative: tmp = alternative
117 Y = "2022"
118 if fillnr > 8500 : Y = "2023"
119 if fillnr >= 9323 : Y="2024"
120 if fillnr >= 10406 : Y="2025"
121 if not self.lpcFillingscheme:
122 with urlopen('https://lpc.web.cern.ch/cgi-bin/fillTable.py?year='+Y) as webpage:
123 self.lpcFillingscheme = webpage.read().decode()
124 self.tagi = '<td>XXXX</td>'
125 self.tagj = 'fillingSchemes/'+Y+'/candidates/'
126 # end of line
127 self.tagl = '</td></tr>'
128 fs = 0
129 i = self.lpcFillingscheme.find(self.tagi.replace('XXXX',str(tmp)))
130 if int(Y)>=2024:
131 if i>0:
132 end_of_line = self.lpcFillingscheme[i:].find(self.tagl)
133 # In rare cases there is a link to download the FS file
134 # One needs to skip this link data
135 end_link_data = self.lpcFillingscheme[i:i+end_of_line].find('</a>')
136 if end_link_data!=-1:
137 end_of_line = self.lpcFillingscheme[i:].find('</a>'+self.tagl)
138 line_of_interest = self.lpcFillingscheme[i:i+end_of_line]
139 last_separator = line_of_interest.rfind("<td>")
140 # Skip download link data, if it exists
141 last_separator_link_check = line_of_interest.rfind('.csv">')
142 if last_separator_link_check!=-1:
143 fs = line_of_interest[last_separator_link_check+6:]
144 else:
145 fs = line_of_interest[last_separator+4:]
146 else:
147 if i>0:
148 j = self.lpcFillingscheme[i:].find(self.tagj)+i+len(self.tagj)
149 if j>0:
150 k = self.lpcFillingscheme[j:].find('.csv')
151 if k>0:
152 fs = self.lpcFillingscheme[j:j+k]
153 return fs
154
155 def getFillNrFromRunNr(self,runNumber):
156 if runNumber in fromElog:
157 return fromElog[runNumber]
158 FILL_NUMBER_MASK = 0x000000000000FFFF
159 R = ROOT.TFile.Open(os.environ['EOSSHIP']+\
160 options.rawData+"/run_"+str(runNumber).zfill(6)+"/data_0000.root")
161 try:
162 if R.Get('event'):
163 rc = R.Get("event").GetEvent(R.Get("event").GetEntries()-1)
164 flags = R.Get("event").flags
165 else:
166 event = R.Get("data")
167 event.GetEvent(0)
168 flags = event.evt_flags
169 fillNumber = numpy.bitwise_and(FILL_NUMBER_MASK,flags)
170 if fillNumber<1: fillNumber = False
171 print('fill number =',fillNumber )
172 except:
173 fillNumber = False
174 if runNumber == 6296 : fillNumber = 8897 # LHC mistake, twice the same fill number
175 return fillNumber
176
177 def getLumiAtIP1(self,fillnr=None,fromnxcals=False, fromAtlas=False):
178 Y = "2022"
179 if fillnr>8500: Y = "2023"
180 if fillnr>=9323: Y="2024"
181 if fillnr>=10406 : Y="2025"
182 if not fromnxcals and not fromAtlas:
183 try:
184 with urlopen('https://lpc.web.cern.ch/cgi-bin/fillAnalysis.py?year='+Y+'&action=fillData&exp=ATLAS&fillnr='+str(fillnr)) as webpage:
185 tmp = webpage.read().decode()
186 except:
187 print('lumi info not avaible from lpc.',fillnr)
188 return -1
189 exec("self.content = "+tmp)
190 """ Each lumi file contains:
191 time stab l dl sl dsl
192 where
193 time = UNIX time, i.e. in seconds since UTC Jan 1, 1970, 00:00:00, not counting leap seconds.
194 stab = stable beam flag: a float value between 0 and 1 which corresponds to the fraction of time spent in stable beams for this time bin.
195 l = luminosity in Hz/ub
196 dl = point-to-point error on luminosity in Hz/ub
197 sl = specific luminosity in Hz/ub (see below)
198 dsl = point-to-point error on specific luminosity in Hz/ub
199 self.content['data']['fillData']['data'][0]
200 [1660892004.2119756, 1.0, 9907.5947265625, 0.0, 3.2100153151839246e-22, 0.0]
201 time.ctime(1660892004.2119756) -> 'Fri Aug 19 08:53:24 2022'
202 """
203 fs = self.getNameOfFillingscheme(fillnr)
204 self.lumiAtIP1 = {'startTime':self.content['data']['fillData']['data'][0][0],'lumiTime':ROOT.TGraph(),'fillingScheme':fs}
205 X = self.content['data']['fillData']['data']
206 t0 = self.lumiAtIP1['startTime']
207 for n in range(len(X)):
208 self.lumiAtIP1['lumiTime'].AddPoint(X[n][0]-t0,X[n][2])
209 return 0
210
211 elif fromnxcals:
212 # nxcals old: /eos/experiment/sndlhc/nxcals_data/
213 fileLoc = os.environ['EOSSHIP']+"/eos/experiment/sndlhc/nxcals_data/fill_"+str(fillnr).zfill(6)+".root"
214 try:
215 fill = ROOT.TFile.Open(fileLoc)
216 except:
217 print('Fill information not found in nxcals ',fillnr)
218 return -1
219
220 if not fill.LuminosityIP1.Get('ATLAS_OFFLINE_LUMI_TOT_INST'):
221 if not fill.LuminosityIP1.Get('ATLAS_LUMI_TOT_INST'):
222 print('Lumi information not found in nxcals ',fillnr)
223 return -1
224 LtreeOff = fill.LuminosityIP1.ATLAS_LUMI_TOT_INST
225 else:
226 LtreeOff = fill.LuminosityIP1.ATLAS_OFFLINE_LUMI_TOT_INST
227# other useful info
228 fillScheme = " "
229 if fill.LHC.FindObjectAny('LHC_STATS_LHC_INJECTION_SCHEME'):
230 rc = fill.LHC.LHC_STATS_LHC_INJECTION_SCHEME.GetEvent(0)
231 fillScheme = str(fill.LHC.LHC_STATS_LHC_INJECTION_SCHEME.var)
232 else:
233 fillScheme = self.getNameOfFillingscheme(fillnr)
234
235 rc = LtreeOff.GetEvent(0)
236 self.lumiAtIP1 = {'startTime':LtreeOff.unix_timestamp,'lumiTime':ROOT.TGraph(),'fillingScheme':fillScheme}
237 t0 = self.lumiAtIP1['startTime']
238 for e in LtreeOff:
239 self.lumiAtIP1['lumiTime'].AddPoint(e.unix_timestamp-t0,e.var)
240
241 return 0
242
243 elif fromAtlas:
244 fileLoc = os.environ['EOSSHIP']+"/eos/experiment/sndlhc/atlas_lumi/fill_"+str(fillnr).zfill(6)+".root"
245 try:
246 fill = ROOT.TFile.Open(fileLoc)
247 except:
248 print('Fill information not found in atlas_lumi ',fillnr)
249 return -1
250
251 LtreeOff = fill.atlas_lumi
252 fillScheme = self.getNameOfFillingscheme(fillnr)
253
254 rc = LtreeOff.GetEvent(0)
255 self.lumiAtIP1 = {'startTime':LtreeOff.unix_timestamp,'lumiTime':ROOT.TGraph(),'fillingScheme':fillScheme}
256 t0 = self.lumiAtIP1['startTime']
257 for e in LtreeOff:
258 self.lumiAtIP1['lumiTime'].AddPoint(e.unix_timestamp-t0,e.var)
259
260 return 0
261 def drawLumi(self,runNumber):
262 R = ROOT.TFile.Open(www+"offline/run"+str(runNumber).zfill(6)+".root")
263 ROOT.gROOT.cd()
264 bCanvas = R.Get("daq").Get('T')
265 Xt = {'time':None,'timeWtDS':None,'timeWt':None}
266 for x in Xt:
267 Xt[x] = bCanvas.FindObject(x)
268 if Xt[x]: self.h[x] = Xt[x].Clone(x)
269 else: self.h[x].Reset()
270 self.h['c1'].cd()
271 ROOT.gROOT.cd()
272 self.options.fillNumbers = self.getFillNrFromRunNr(runNumber)
273 # for overlay, need SND@LHC startTime to adjust with lumiTime
274 runDir = options.rawData+"/run_"+str(runNumber).zfill(6)
275 jname = "run_timestamps.json"
276 dirlist = str( subprocess.check_output("xrdfs "+os.environ['EOSSHIP']+" ls "+runDir,shell=True) )
277 startTime=0
278 if jname in dirlist:
279 with client.File() as f:
280 f.open(os.environ['EOSSHIP']+runDir+"/run_timestamps.json")
281 status, jsonStr = f.read()
282 exec("self.date = "+jsonStr.decode())
283 time_str = self.date['start_time'].replace('Z','')
284 time_obj = time.strptime(time_str, '%Y-%m-%dT%H:%M:%S')
285 self.startTime = calendar.timegm(time_obj)
286 else: # try reading ecs log file
287 if len(self.startTimes)==0: self.readStartTimes()
288 if runNumber in self.startTimes:
289 self.startTime = self.startTimes[runNumber]
290 else: return
291 self.date['start_time'] = time.ctime(self.startTime)
292
293 for x in Xt:
294 if not Xt[x]: continue
295 self.h[x+'10'] = self.h[x].Clone(x+'10')
296 self.h[x+'10'].Rebin(10)
297 self.h[x+'10'].SetMinimum(0)
298 self.h[x+'10'].Scale(self.options.postScale/10.)
299 if self.date['start_time'].find('Tu')<0 and self.date['start_time'].find('Th')<0:
300 tmp = self.date['start_time'].replace('T',' ').replace('Z','')
301 else: tmp = self.date['start_time']
302 self.h['time10'].SetTitle('Run '+str(runNumber)+' Fill '+str(self.options.fillNumbers)+' '+tmp)
303# what to do with spikes?
304 mx = self.h['time10'].GetMaximumBin()
305 side = self.h['time10'].GetBinContent(mx-1)+self.h['time10'].GetBinContent(mx+1)
306 if self.h['time10'].GetBinContent(mx+1)<0.2*self.h['time10'].GetBinContent(mx-1): side = 2*self.h['time10'].GetBinContent(mx-1) # happens at end of fill
307 newMx = max(10,0.75*side)
308# special runs
309 if runNumber==4423: newMx = 100
310 if runNumber==4415: newMx = 150
311 if runNumber==4362: newMx = 7
312 if runNumber==5003: newMx = 4100
313 if self.h['time10'].GetBinContent(mx) > side: self.h['time10'].SetMaximum(newMx)
314 self.h['time10'].SetMinimum(0)
315 self.h['time10'].Draw('hist')
316 self.h['timeWt10'].Draw('histsame')
317 self.h['timeWtDS10'].Draw('histsame')
318 self.h['c1'].Update()
319 self.myPrint('c1','noLumi-run'+str(runNumber).zfill(6))
320
321 nbins = self.h['time'].GetNbinsX()
322 endTime = self.h['time'].GetBinCenter(nbins) # in seconds
323 rc = self.getLumiAtIP1(fillnr=options.fillNumbers,fromnxcals=True,fromAtlas=False)
324 if not rc<0:
325 if FS.lumiAtIP1['lumiTime'].GetN()<2:
326 rc = self.getLumiAtIP1(fillnr=options.fillNumbers,fromnxcals=False,fromAtlas=True)
327 if rc < 0:
328 rc = self.getLumiAtIP1(fillnr=options.fillNumbers,fromnxcals=False,fromAtlas=False)
329 if rc<0: return
330
331 self.lumiAtlas = ROOT.TGraph()
332 deltaT = self.lumiAtIP1['startTime'] - self.startTime # account for timezone/summertime
333 self.Lmax = 0
334 self.Lint = 0
335 self.Lsnd = 0
336 tprev = [-1,0]
337 for n in range(self.lumiAtIP1['lumiTime'].GetN()):
338 t = self.lumiAtIP1['lumiTime'].GetPointX(n) + deltaT
339 l = self.lumiAtIP1['lumiTime'].GetPointY(n)
340 self.lumiAtlas.AddPoint(t,l)
341 if l>self.Lmax: self.Lmax = l
342 if tprev[0] < 0:
343 tprev = [t,l]
344 else:
345 dt = t - tprev[0]
346 X = (tprev[1]+l)/2.
347 self.Lint +=X*dt
348 if t>0 and t<endTime: self.Lsnd += X*dt
349 tprev[0] = t
350 tprev[1] = l
351 self.LumiInt[runNumber] = [self.Lint,self.Lsnd]
352
353 if not self.Lmax>0:
354 print('no lumi for run ',runNumber)
355 return
356# work with second axis
357 rightmax = 1.1*self.Lmax/1000
358 self.scale = ROOT.gPad.GetUymax()/rightmax
359 self.lumiAtlasS = ROOT.TGraph()
360 self.lumiAtlasS.SetLineColor(ROOT.kMagenta)
361 self.lumiAtlasS.SetLineWidth(3)
362
363 for n in range(self.lumiAtIP1['lumiTime'].GetN()):
364 t = self.lumiAtIP1['lumiTime'].GetPointX(n) + deltaT
365 l = self.lumiAtIP1['lumiTime'].GetPointY(n)
366 self.lumiAtlasS.AddPoint(t,l*self.scale/1000)
367 self.lumiAtlasS.Draw('same')
368 self.h['ax1'] = ROOT.TGaxis(ROOT.gPad.GetUxmax(), ROOT.gPad.GetUymin(),
369 ROOT.gPad.GetUxmax(), ROOT.gPad.GetUymax(),
370 0, rightmax, 510, "+L")
371 l = self.Lsnd/1E9
372 ul = 'fb'
373 if l < 0.01:
374 l = self.Lsnd/1E9
375 ul = 'pb'
376 self.h['ax1'].SetTitle('L [Hz/nb] Integral '+"%5.2F %s^{-1}"%(l,ul))
377 self.h['ax1'].SetTextFont(42)
378 self.h['ax1'].SetLabelFont(42)
379 self.h['ax1'].SetTextColor(ROOT.kMagenta)
380 self.h['ax1'].Draw()
381 self.h['time10'].Draw('histsame')
382 self.myPrint('c1','Lumi-run'+str(runNumber).zfill(6))
383 def alternativeFill(self,fillNr):
384 alternative = None
385 if fillNr=='8470':
386 # 25ns_156b_144_90_96_48bpi_4inj_MD7003 from ELOG
387 alternative = '8471'
388 if fillNr=='8461':
389 # 2nominals_10pilots_lossmaps_coll_allIPs from ELOG
390 alternative = '8184'
391 if fillNr=='8256':
392 # 8256 = 25ns_2461b_2448_1737_1733_180bpi_16inj_1INDIV
393 alternative = '8253'
394 if fillNr=='8140':
395 # 8140 = 25ns_2413b_2400_1836_1845_240bpi_12inj_1INDIV, same as 8142
396 alternative = '8142'
397 elif fillNr=='8105':
398 # from Cris: '25ns_2173b_2160_1804_1737_240bpi_11inj_1INDIV' wrong
399 # try: '25ns_1935b_1922_1758_1842_240bpi_12inj_3INDIV'
400 # from elog : 25ns_1935b_1922_1602_1672_192bpi_14inj_3INDIVs
401 alternative = '8106'
402 elif fillNr=='8074':
403 # 8073 25ns_1227b_1214_1054_1102_144bpi_14inj
404 # 8076 25ns_1551b_1538_1404_1467_144bpi_16inj
405 alternative = '8073'
406 elif fillNr=='8070':
407 # 8068 and 8072: 25ns_1227b_1214_1054_1102_144bpi_14inj
408 alternative = '8068'
409 elif fillNr=='8056':
410 # from Cris: '25ns_987b_974_876_912_96bpi_17inj'
411 alternative = '8057'
412 elif fillNr=='8045':
413 # 8043 25ns_987b_974_878_917_144bpi_13inj
414 alternative = '8043'
415 elif fillNr=='8025':
416 # 8023 25ns_603b_590_524_542_48bpi_17inj
417 # 8027 25ns_603b_590_526_547_96bpi_13inj
418 alternative = '8023' # ???
419 elif fillNr=='8011':
420 # 8016 rampup 25ns_603b_590_524_542_48bpi_17inj UFO true false First 600b ramp-up fill
421 # 8007 rampup 25ns_315b_302_237_240_48bpi_11inj RF false false Third 300b ramp-up fill
422 alternative = '8007' # ???
423 elif fillNr=='8388':
424 # 25ns_2462b_2450_1737_1735_180bpi_17inj_2INDIV
425 alternative = '8387' #
426 elif fillNr=='8478':
427 # Single_12b_9_1_3_BSRT_2018_pilot
428 alternative = '8479' #
429 elif fillNr=='8294':
430 # 25ns_315b_302_237_240_48bpi_11inj,
431 alternative = '8295' #? 8293' #
432 return alternative
433
434 def extractFillingScheme(self,fillNr):
435 alternative = self.alternativeFill(str(fillNr))
436 # Get the FS name from the all-year LPC table
437 fs_name_table = self.getNameOfFillingscheme(int(fillNr))
438 # this 2024 ion run FS has its attributes swapped btw LPC table and json
439 # this is the most elegant fix since many fills are affected
440 if fs_name_table=="50ns_119b_58_51_58_56bpi_9inj_3INDIV_4NC_PbPb":
441 fs_name_table="50ns_119b_58_51_58_56bpi_9inj_4NC_3INDIV_PbPb"
442 print('Name of filling scheme: ',fs_name_table)
443 if fs_name_table==0:
444 return -1
445 # since 2024 there is new storage of FS in json files, no csv
446 fs_url="https://gitlab.cern.ch/lhc-injection-scheme/injection-schemes/-/raw/master/"+\
447 fs_name_table+".json"
448 F = ROOT.TFile(self.path+'fillingScheme-'+fillNr+'.root','recreate')
449 nt = ROOT.TNtuple('fill'+fillNr,'b1 IP1 IP2','B1:IP1:IP2:IsB2')
450 # Get the data for year>=2024
451 # 9323 is first fillNr for 2024
452 if int(fillNr) >= 9323:
453 # Get the JSON file content from the web
454 response = requests.get(fs_url)
455 # If JSON file is not found, try adding char 's' to FS name
456 if response.status_code == 404:
457 fs_url="https://gitlab.cern.ch/lhc-injection-scheme/injection-schemes/-/raw/master/"+\
458 fs_name_table+"s.json"
459 response = requests.get(fs_url)
460 response.raise_for_status()
461 fs_data_json = response.json()
462 # check if collision data exists - for 2024 ion runs it doesn't
463 # and SND files from eos are to be used.
464 # These files were generated using the LPC Filling Scheme Editor.
465 if 'collsPatternB1' not in fs_data_json:
466 with open('/eos/experiment/sndlhc/filling_schemes/2024/ion_run/'+
467 fs_name_table+'_snd_generated.json', 'r') as custom_json:
468 fs_data_json = json.load(custom_json)
469 print('Using FS generated using LPC Filling Scheme Editor!\nJSON:'\
470 '/eos/experiment/sndlhc/filling_schemes/2024/ion_run/'+
471 fs_name_table+'_snd_generated.json')
472 # check name of file and fs_name in file are the same
473 if 'schemeName' in fs_data_json:
474 fs_name_json = fs_data_json['schemeName']
475 if fs_name_json!=fs_name_table:
476 print('FS name differs btw the LPC JSON and the all-year LPC table, check!', '\n', \
477 'JSON:', fillNr, fs_name_json, '\n', \
478 'all-year table:', fs_name_table)
479 # Accept if difference is char 's'
480 if fs_name_json!=fs_name_table+'s':
481 return -1
482
483 nB1 = fs_data_json['collsPatternB1']#B1 bucket number,IP1,IP2,IP5,IP8
484 nB2 = fs_data_json['collsPatternB2']#B2 bucket number,IP1,IP2,IP5,IP8
485 # Sanity checks.
486 if not len(nB1)==5 or not len(nB2)==5:
487 print("missing collsPattern data in the FS JSON file")
488 return -1
489
490 # Get N colliding bunches from name of the FS json file
491 # Convention is {spacing}_{bunches}_{IP1/5}_{IP2}_{IP8}_{trainlength}_{injections}_{special info}
492 fs_n_bunches = re.findall(r'(?<=_)\d+(?=_)', fs_name_table)
493 n_ip1_in_title = int(fs_n_bunches[0])
494 n_ip2_in_title = int(fs_n_bunches[1])
495 n_ip8_in_title = int(fs_n_bunches[2])
496
497 summary_collisions = {'B1':[], 'B2':[]}
498 collsPatterns = {'B1':nB1, 'B2':nB2}
499 for beam in collsPatterns.keys():
500 for i in range(1,5):
501 summary_collisions[beam].append(numpy.count_nonzero(numpy.array(collsPatterns[beam][i])))
502 if i==1 or i==3: #IP1/5
503 if not summary_collisions[beam][i-1]== n_ip1_in_title:
504 print("For the number of IP1/5 bunches for {0} got {1} expected {2}. Check the FS!".format(beam,summary_collisions[beam][i-1],n_ip1_in_title))
505 if beam=='B1': # asymmetric collisions in IP2 and IP8
506 if i==2: #IP2
507 if not summary_collisions[beam][i-1]== n_ip2_in_title:
508 print("For the number of IP2 bunches for {0} got {1} expected {2}. Check the FS!".format(beam,summary_collisions[beam][i-1],n_ip2_in_title))
509 if i==4: #IP8
510 if not summary_collisions[beam][i-1]== n_ip8_in_title:
511 print("For the number of IP8 bunches for {0} got {1} expected {2}. Check the FS!".format(beam,summary_collisions[beam][i-1],n_ip8_in_title))
512 print("From the json file, we got:\nBeam1:\nIP1 {0}, IP2 {1}, IP5 {2}, IP8 {3}".format(*summary_collisions['B1']))
513 print("Beam2:\nIP1 {0}, IP2 {1}, IP5 {2}, IP8 {3}".format(*summary_collisions['B2']))
514
515 for i in range(len(nB1[0])):
516 if nB1[1][i]==1: # 1 if headon in IP1
517 b1_ip1_id = int(nB1[0][i])
518 else: b1_ip1_id = -1
519 if nB1[2][i]==1: # 1 if headon in IP2
520 b1_ip2_id = int(nB1[0][i])
521 else: b1_ip2_id = -1
522 rc = nt.Fill(int(nB1[0][i]),b1_ip1_id,b1_ip2_id,0)
523 for i in range(len(nB2[0])):
524 if nB2[1][i]==1: # 1 if headon in IP1
525 b2_ip1_id = int(nB2[0][i])
526 else: b2_ip1_id = -1
527 if nB2[2][i]==1: # 1 if headon in IP2
528 b2_ip2_id = int(nB2[0][i])
529 else: b2_ip2_id = -1
530 rc = nt.Fill(int(nB2[0][i]),b2_ip1_id, b2_ip2_id,1)
531
532 else:
533 # 2022-2023 FS storage in csv files
534 # cases where only a binary csv file exists or the json file is wrong
535 if fillNr in ['9231', '9232']:
536 File = urlopen('https://lpc.web.cern.ch/fillingSchemes/2023/candidates/50ns_1123b_1010_1010_320_56bpi_23inj_3INDIV_PbPb.csv')
537 X = File.read()
538 File.close()
539 csv = X.decode().split('\n')
540 else:
541 fs_url="https://lpc.web.cern.ch/cgi-bin/schemeInfo.py?fill="
542 if alternative:
543 fs_url=fs_url+alternative+'&fmt=json'
544 else:
545 fs_url = fs_url+fillNr+'&fmt=json'
546 with urlopen(fs_url) as webpage:
547 tmp = webpage.read().decode()
548
549 exec("self.content = "+tmp)
550 if len(self.content['fills']) < 1:
551 print('Filling scheme not yet known',fillNr,self.options.runNumbers)
552 return -1
553 if alternative:
554 self.content['fills'][fillNr] = self.content['fills'][alternative]
555 if (self.content['fills'][fillNr]['name'] != fs_name_table ):
556 print('FS data differs btw the LPC JSON and the all-year LPC table, check!', '\n', \
557 'JSON:', fillNr, self.content['fills'][fillNr]['name'], '\n', \
558 'all-year table:', fs_name_table, '\n', \
559 'One can look for the candidates csv files here\n' \
560 "https://lpc.web.cern.ch/fillingSchemes/202X/candidates")
561 return -1
562 csv = self.content['fills'][fillNr]['csv'].split('\n')
563
564 nB1 = csv.index('B1 bucket number,IP1,IP2,IP5,IP8')
565 while nB1>0:
566 tmp = csv[nB1+1].split(',')
567 if len(tmp)!=5: break
568 nB1+=1
569 rc = nt.Fill(int(tmp[0]),int(tmp[1].replace('-','-1')),int(tmp[2].replace('-','-1')),0)
570 nB2 = csv.index('B2 bucket number,IP1,IP2,IP5,IP8')
571 while nB2>0:
572 tmp = csv[nB2+1].split(',')
573 if len(tmp)!=5: break
574 nB2+=1
575 rc = nt.Fill(int(tmp[0]),int(tmp[1].replace('-','-1')),int(tmp[2].replace('-','-1')),1)
576
577 nt.Write()
578 F.Close()
579 return 0
580
581 def extractPhaseShift(self,fillNr,runNumber):
582 # Check if the offline monitoring file exists
583 try:
584 R = ROOT.TFile.Open(www+"offline/run"+str(runNumber).zfill(6)+".root")
585 ROOT.gROOT.cd()
586 try:
587 self.h['bnr'] = R.Get("daq").Get('bunchNumber').FindObject('bnr').Clone('bnr')
588 except:
589 self.h['bnr'] = R.Get("daq").Get('shifter/bunchNumber').FindObject('bnr').Clone('bnr')
590 R.Close()
591 # create the bunch number plot if offline monitoring file is missing
592 except:
593 try:
594 self.h['bnr']=self.h['bnr_from_data']
595 except:
596 self.h['bnr'] = self.BunchNumberPlotFromData(runNumber)
597 Nbunches = self.h['bnr'].GetNbinsX()
598#Filling scheme
599 self.F = ROOT.TFile(self.path+'fillingScheme-'+fillNr+'.root')
600 self.fs = self.F.Get('fill'+fillNr)
601# convert to dictionary
602 self.FSdict[runNumber] = {"fillNumber":fillNr,'phaseShift1':0,'phaseShift2':0,"B1":{},"B2":{}}
603 fsdict = self.FSdict[runNumber]
604 for x in self.fs:
605 if x.IsB2>0:
606 if Nbunches==3564:
607 fsdict['B2'][(x.B1-1)/10]={'IP1':x.IP1>0,'IP2':x.IP2>0}
608 if Nbunches==1782:
609 fsdict['B2'][(x.B1-1)//20]={'IP1':x.IP1>0,'IP2':x.IP2>0}
610 else:
611 if Nbunches==3564:
612 fsdict['B1'][(x.B1-1)/10]={'IP1':x.IP1>0,'IP2':x.IP2>0}
613 if Nbunches==1782:
614 fsdict['B1'][(x.B1-1)//20]={'IP1':x.IP1>0,'IP2':x.IP2>0}
615 self.matches = {}
616 for phase1 in range(0,Nbunches):
617 self.matches[phase1]=0
618 for n in range(0,Nbunches):
619 if not n in fsdict['B1']: continue
620 j = (n+phase1)%Nbunches + 1
621 if fsdict['B1'][n]['IP1']: self.matches[phase1]+=self.h['bnr'].GetBinContent(j)
622 self.phaseShift1 = max(self.matches,key=self.matches.get)
623 print('phaseShift1 found:',self.phaseShift1,Nbunches-self.phaseShift1)
624 self.matches = {}
625 if Nbunches == 1782: self.phaseShift2 = Nbunches-64
626 else:
627 for phase2 in range(0,Nbunches):
628 self.matches[phase2]=0
629 for n in range(0,Nbunches):
630 if not n in fsdict['B2']:
631 continue
632 # calculate the bunch number of SND events
633 j = (n+self.phaseShift1+phase2)%Nbunches + 1 # bin number
634 # adjust it with the found B1 phase to determine(and exclude)
635 # bunches associated with IP1 collisions
636 ip1 = (j-1+Nbunches-self.phaseShift1)%Nbunches
637 if ip1 in fsdict['B1']:
638 if fsdict['B1'][ip1]['IP1']: continue
639 if fsdict['B2'][n]['IP2'] or 1>0:
640 self.matches[phase2]+=self.h['bnr'].GetBinContent(j)
641 self.phaseShift2 = max(self.matches,key=self.matches.get)
642 print('phaseShift2 found:',self.phaseShift2,Nbunches-self.phaseShift2)
643 if not (Nbunches-self.phaseShift2) == 129:
644 print('There is a problem with phaseshift2 for run',runNumber,Nbunches-self.phaseShift2)
645 if self.phaseShift2 == 129:
646 print("!!! Probably beam 1 and beam 2 are interchanged. Try reverse")
647 self.phaseShift1 = FS.phaseShift1 +129
648 if fillNr=="8178":
649 print('special LHCf run. Phaseshift determined by hand: 430-1017+3564')
650 self.phaseShift1 = 430-1017+3564
651 if fillNr=="8056":
652 print('run with very low lumi at IP1, beam 2 background dominates')
653 self.phaseShift1 = 3564-1732+129
654 if fillNr=="8056":
655 print('run with very low lumi at IP1, beam 2 background dominates')
656 self.phaseShift1 = 3564-1603
657 if fillNr=="8140" or fillNr=="8070" or fillNr=="8045":
658 print('run with very low lumi at IP1, beam 2 background dominates')
659 self.phaseShift1 = 3564-1456
660 if fillNr == "8383":
661 print('run with very low lumi at IP1, fit does not converge correctly')
662 self.phaseShift1 = 1978 + 130
663 if fillNr == "8342":
664 print('run with very low lumi at IP1, beam 2 background dominates')
665 self.phaseShift1 = 2107
666 if fillNr == "8256":
667 print('run with very low lumi at IP1, beam 2 background dominates')
668 self.phaseShift1 = 2108
669 if fillNr == "8294":
670 print('run with very low lumi at IP1, beam 2 background dominates')
671 self.phaseShift1 = 2598 +129
672 # force 129 for proton beams
673 # It is determined using ~480m distance to IP1 and 25-ns distance between bunches:
674 # 482.5m/speed of light/bunch distance = 64.33
675 # This is multiplied by 2 to arrive at 129.
676 self.phaseShift2 = Nbunches - 129
677 print('phaseShift2 is set to:',self.phaseShift2,Nbunches-self.phaseShift2)
678 fsdict['phaseShift2'] = self.phaseShift2
679 fsdict['phaseShift1'] = self.phaseShift1
680
681 def plotBunchStructure(self,fillNr,runNumber):
682 h=self.h
683 self.F = ROOT.TFile(self.path+'fillingScheme-'+fillNr+'.root')
684 self.fs = self.F.Get('fill'+fillNr)
685
686 try:
687 R = ROOT.TFile.Open(www+"offline/run"+str(runNumber).zfill(6)+".root")
688 ROOT.gROOT.cd()
689 bCanvas = R.Get("daq").Get('bunchNumber')
690 if not bCanvas:
691 bCanvas = R.Get("daq").Get("shifter").Get('bunchNumber')
692 h['bnr']= bCanvas.FindObject('bnr').Clone('bnr')
693 # create the bunch number plot if offline monitoring file is missing
694 except:
695 try:
696 h['bnr'] = h["bnr_from_data"]
697 except:
698 h['bnr'] = self.BunchNumberPlotFromData(runNumber)
699 Nbunches = h['bnr'].GetNbinsX()
700 ROOT.gROOT.cd()
701
702 ut.bookHist(h,'b1','b1',35640,-0.5,35639.5)
703 ut.bookHist(h,'IP1','IP1',35640,-0.5,35639.5)
704 ut.bookHist(h,'IP2','IP2',35640,-0.5,35639.5)
705 ut.bookHist(h,'b1z','b1',Nbunches,-0.5,Nbunches-0.5)
706 ut.bookHist(h,'b2z','b2',Nbunches,-0.5,Nbunches-0.5)
707 ut.bookHist(h,'IP1z','IP1',Nbunches,-0.5,Nbunches-0.5)
708 ut.bookHist(h,'IP2z','IP2',Nbunches,-0.5,Nbunches-0.5)
709 h['b1'].Draw()
710 h['b1z'].SetLineColor(ROOT.kBlue)
711 h['b2z'].SetLineColor(ROOT.kCyan)
712 h['IP1z'].SetLineColor(ROOT.kRed)
713 h['IP2z'].SetLineColor(ROOT.kOrange)
714 h['b1z'].SetStats(0)
715 h['b2z'].SetStats(0)
716 h['IP1z'].SetStats(0)
717 h['IP2z'].SetStats(0)
718
719 self.Draw()
720
721 def Draw(self):
722 h = self.h
723 Nbunches = h['bnr'].GetNbinsX()
724 eq = ''
725 if Nbunches==1782:
726 equation_str ='floor((B1-1)/20)'
727 if Nbunches==3564:
728 equation_str = '(B1-1)/10'
729 h['c1'].cd()
730 self.fs.Draw('(( '+equation_str+'+'+str(self.phaseShift1)+')%'+str(Nbunches)+')>>b1z','!(IsB2>0)','hist')
731 self.fs.Draw('(( '+equation_str+'+'+str(self.phaseShift1)+')%'+str(Nbunches)+')>>IP1z','IP1>-0.6&&(!(IsB2>0))','hist')
732 self.fs.Draw('(( '+equation_str+'+'+str(self.phaseShift1+ self.phaseShift2)+')%'+str(Nbunches)+')>>IP2z','IP2>-0.6&&IsB2>0','hist')
733 self.fs.Draw('(( '+equation_str+'+'+str(self.phaseShift1+ self.phaseShift2)+')%'+str(Nbunches)+')>>b2z','IsB2>0','hist')
734 norm = h['bnr'].GetBinContent(h['bnr'].GetMaximumBin())
735 h['b1z'].Scale(norm*1.5)
736 h['IP1z'].Scale(norm*1.0)
737 h['b2z'].Scale(norm*0.5)
738 h['IP2z'].Scale(norm*0.3)
739 h['bnr'].SetStats(0)
740 h['bnr'].SetFillColor(17)
741 h['bnr'].SetLineColor(ROOT.kBlack)
742 h['b1z'].Draw('hist')
743 h['bnr'].Draw('histsame')
744 h['b1z'].Draw('histsame')
745 txt = 'phase shift B1, B2: '+str(Nbunches-self.phaseShift1)+','+str(Nbunches-self.phaseShift2)+' for run '+str(options.runNumbers)
746 txt += " fill nr "+options.fillNumbers
747 h['b1z'].SetTitle(txt)
748 if self.options.withIP2:
749 h['IP2z'].Draw('histsame')
750 h['b2z'].Draw('histsame')
751 h['IP1z'].Draw('histsame')
752 h['bnr'].Draw('histsame')
753
754 def Xbunch(self):
755 h = self.h
756 ut.bookHist(h,'Xb1z','b1',3564*4,-0.5,3564*4-0.5)
757 ut.bookHist(h,'Xb2z','b2',3564*4,-0.5,3564*4-0.5)
758 ut.bookHist(h,'XIP1z','IP1',3564*4,-0.5,3564*4-0.5)
759 ut.bookHist(h,'XIP2z','IP2',3564*4,-0.5,3564*4-0.5)
760 h['Xb1z'].SetLineColor(ROOT.kBlue)
761 h['Xb2z'].SetLineColor(ROOT.kCyan)
762 h['XIP1z'].SetLineColor(ROOT.kRed)
763 h['XIP2z'].SetLineColor(ROOT.kOrange)
764 h['Xb1z'].SetStats(0)
765 h['Xb2z'].SetStats(0)
766 h['XIP1z'].SetStats(0)
767 h['XIP2z'].SetStats(0)
768 R = ROOT.TFile.Open(www+"offline/run"+str(runNumber).zfill(6)+".root")
769 ROOT.gROOT.cd()
770 bCanvas = R.Get("daq").Get('sndclock')
771 h['Xbnr']= bCanvas.FindObject('Xbnr').Clone('Xbnr')
772 ROOT.gROOT.cd()
773 h['c1'].cd()
774 for j in range(4):
775 self.fs.Draw('(( (B1-1)/2.5+'+str(self.phaseShift1*4-2)+')%(3564*4))+'+str(j)+'>>+Xb1z','!(IsB2>0)','hist')
776 self.fs.Draw('(( (B1-1)/2.5+'+str(self.phaseShift1*4-2)+')%(3564*4))+'+str(j)+'>>+XIP1z','IP1>-0.6&&(!(IsB2>0))','hist')
777 self.fs.Draw('(( (B1-1)/2.5+'+str(self.phaseShift1*4-2+ self.phaseShift2*4)+')%(3564*4))+'+str(j)+'>>+XIP2z','IP2>-0.6&&IsB2>0','hist')
778 self.fs.Draw('(( (B1-1)/2.5+'+str(self.phaseShift1*4-2+ self.phaseShift2*4)+')%(3564*4))+'+str(j)+'>>+Xb2z','IsB2>0','hist')
779 norm = h['Xbnr'].GetBinContent(h['Xbnr'].GetMaximumBin())
780 h['Xb1z'].Scale(norm*1.5)
781 h['XIP1z'].Scale(norm*1.0)
782 h['Xb2z'].Scale(norm*0.5)
783 h['XIP2z'].Scale(norm*0.3)
784 h['Xbnr'].SetStats(0)
785 h['Xbnr'].SetFillColor(17)
786 h['Xbnr'].SetLineColor(ROOT.kBlack)
787 h['Xb1z'].Draw('hist')
788 h['Xbnr'].Draw('histsame')
789 h['Xb1z'].Draw('histsame')
790 Nbunches = h['bnr'].GetNbinsX()
791 txt = 'phase shift B1, B2: '+str(Nbunches-self.phaseShift1)+','+str(Nbunches-self.phaseShift2)+' for run '+str(options.runNumbers)
792 txt += " fill nr "+options.fillNumbers
793 h['Xb1z'].SetTitle(txt)
794 h['XIP2z'].Draw('histsame')
795 h['Xb2z'].Draw('histsame')
796 h['XIP1z'].Draw('histsame')
797 # overlay all snd subcycles
798 ut.bookHist(h,'scycle','overlay',4,-0.5,3.5)
799 if Nbunches == 3564: div=4
800 if Nbunches == 1782: div=8
801 for n in range(h['XIP1z'].GetNbinsX()):
802 if h['XIP1z'].GetBinContent(n+1)>0:
803 rc = h['scycle'].Fill(n%div, h['Xbnr'].GetBinContent(n+1))
804
805 def Extract(self):
806 if self.options.fillNumbers=='':
807 fillNumber = self.getFillNrFromRunNr(int(options.runNumbers))
808 if not fillNumber:
809 print('Fill number not found')
810 else:
811 rc = self.extractFillingScheme(str(fillNumber))
812 if not rc<0:
813 self.options.fillNumbers = str(fillNumber)
814 self.extractPhaseShift(self.options.fillNumbers,int(self.options.runNumbers))
815 r = int(self.options.runNumbers)
816 self.plotBunchStructure(self.options.fillNumbers,r)
817 self.myPrint('c1','FS-run'+str(r).zfill(6))
818 # add the FS to the file without running all other modules
819 self.merge()
820 self.storeDict(self.FSdict,'FSdict','FSdict')
821
822 else:
823 for r in options.fillNumbers.split(','):
825
826 def test(self,runnr,I=True):
827 h = self.h
828 Nbunches = h['bnr'].GetNbinsX()
829 p = open("FSdict.pkl",'rb')
830 self.FSdict = pickle.load(p)
831 if runnr in self.FSdict: fsdict = self.FSdict[runnr]
832 else:
833 print('run number not known. Exit.')
834 return()
835 options.runNumbers = str(runnr)
836 fillNr = fsdict['fillNumber']
837 if I:
838 f=open('fill'+str(fillNr)+'.csv')
839 X = f.readlines()
840 keys = {'A6R4.B1':0,'A6R4.B2':0,'B6R4.B1':0,'B6R4.B2':0}
841 for k in keys:
842 n=0
843 for l in X:
844 if l.find(k)>0: keys[k]=n
845 n+=1
846 print(keys)
847 A = X[keys['A6R4.B2']//2].replace('\n','').split(',')
848 print(len(A))
849 for bunchNumber in range(1,Nbunches+1):
850 b1 = (bunchNumber-1) in fsdict['B1']
851 if b1 and float(A[bunchNumber])>1: continue
852 if not b1 and float(A[bunchNumber])<1: continue
853 print(bunchNumber,A[bunchNumber],b1)
854 return
855 FS.phaseShift1 = fsdict['phaseShift1']
856 self.plotBunchStructure(fsdict['fillNumber'],runnr)
857 w = {}
858 for x in ['b1z','IP1z','b2z','IP2z']:
859 w[x]=h[x].GetMaximum()
860 h[x].Reset()
861 for bunchNumber in range(0,Nbunches):
862 nb1 = ( Nbunches+bunchNumber - fsdict['phaseShift1'])%Nbunches
863 if nb1 in fsdict['B1']:
864 rc = h['b1z'].Fill(bunchNumber,w['b1z'])
865 if fsdict['B1'][nb1]['IP1']: rc = h['IP1z'].Fill(bunchNumber,w['IP1z'])
866 nb2 = ( Nbunches + bunchNumber - fsdict['phaseShift1'] - fsdict['phaseShift2'])%Nbunches
867 if nb2 in fsdict['B2']:
868 rc = h['b2z'].Fill(bunchNumber,w['b2z'])
869 if fsdict['B2'][nb2]['IP2']: rc = h['IP2z'].Fill(bunchNumber,w['IP2z'])
870
871 def b1b2(self,runnr,b):
872 Nbunches = self.h['bnr'].GetNbinsX()
873 if runnr in self.FSdict:
874 fsdict = self.FSdict[runnr]
875 nb1 = ( Nbunches + b - fsdict['phaseShift1'])%Nbunches
876 nb2 = ( Nbunches + b - fsdict['phaseShift1'] - fsdict['phaseShift2'])%Nbunches
877 print('b1 bunch number',nb1,nb2)
878
879 def calcMu(self):
880 sigma = 80E6 # 80mb
881 Nbunches = self.h['bnr'].GetNbinsX()
882 self.L = ROOT.TFile.Open("Lumi.root")
883 L = self.L
884 h = self.h
885 runInfo = self.runInfo
886 h['muAv']={}
887 for k in L.GetListOfKeys():
888 runNumber = k.GetName()[3:]
889 tc = L.Get(k.GetName())
890 h[runNumber+'_Mu']=ROOT.TGraph()
891 # need to know the number of colliding bunches. Take from filling scheme
892 fs = runInfo[int(runNumber)]['FillingScheme']
893 if fs==' ' or fs=='' or fs==0:
894 print('filling scheme not in runinfo, try ',runNumber)
895 fillnr = runInfo[int(runNumber)]['Fillnumber']
896 fs = self.getNameOfFillingscheme(fillnr)
897 if fs==' ' or fs=='' or fs==0:
898 print('filling scheme not found')
899 continue
900 else:print('filling scheme found!')
901 runInfo[int(runNumber)]['FillingScheme'] = fs
902 tag = 2
903 if fs.find('Multi')==0: tag=3
904 IP1 = int(fs.split('_')[tag])
905 collPerTurn = IP1/Nbunches
906 i=-1
907 for x in tc.GetListOfPrimitives():
908 i+=1
909 if x.ClassName() == "TGraph": gn = i
910 if x.ClassName() == "TGaxis": an = i
911 if x.ClassName() == "TFrame": fn = i
912 if x.GetName().find("timeWt")==0:
913 h[x.GetName()] = x.Clone(x.GetName())
914 name = runNumber+'_'+x.GetName()+'_Mu'
915 h[name] = x.Clone(name)
916 g = tc.GetListOfPrimitives()[gn]
917 axis = tc.GetListOfPrimitives()[an]
918 frame = tc.GetListOfPrimitives()[fn]
919 scale = (axis.GetWmax()-axis.GetWmin())/(frame.GetY2()-frame.GetY1())
920 MuMax = 0
921 for n in range(g.GetN()):
922 mu = g.GetPointY(n) * scale * sigma / collPerTurn * 25E-9
923 h[runNumber+'_Mu'].AddPoint(g.GetPointX(n),mu)
924 if mu>MuMax: MuMax=mu
925 for t in ['timeWt10','timeWtDS10']:
926 hmu = runNumber+'_'+t+'_Mu'
927 for i in range(1,h[t].GetNbinsX()+1):
928 h[hmu].SetBinContent(i,h[t].GetBinContent(i)/collPerTurn * 25E-9)
929 h[hmu].SetBinError(i,h[t].GetBinError(i)/collPerTurn * 25E-9)
930# find mean values start and end time = 10% of max
931 for n in range(h[runNumber+'_Mu'].GetN()):
932 if h[runNumber+'_Mu'].GetPointY(n)>0.1*MuMax:
933 startT = n
934 break
935 for n in range(h[runNumber+'_Mu'].GetN(),0,-1):
936 if h[runNumber+'_Mu'].GetPointY(n)>0.1*MuMax:
937 endT = n
938 break
939 rc = h[runNumber+'_Mu'].Fit('pol0','SQ','',h[runNumber+'_Mu'].GetPointX(startT),h[runNumber+'_Mu'].GetPointX(endT))
940 res = rc.Get()
941 if not res: print('calcMu: something went wrong',runNumber)
942 muAv = res.Parameter(0)
943 h['muAv'][runNumber] = {'':muAv,'Scifi':0,'DS':0}
944 for t in ['timeWt10','timeWtDS10']:
945 hmu = runNumber+'_'+t+'_Mu'
946 if h[hmu].GetSumOfWeights()==0:continue
947 for n in range(1,h[hmu].GetNbinsX()):
948 if h[hmu].GetBinContent(n)>0.1*h[hmu].GetMaximum():
949 startT = h[hmu].GetBinCenter(n)
950 break
951 for n in range(h[hmu].GetNbinsX(),1,-1):
952 if h[hmu].GetBinContent(n)>0.1*h[hmu].GetMaximum():
953 endT = h[hmu].GetBinCenter(n)
954 break
955 rc = h[hmu].Fit('pol0','SQ','',startT,endT)
956 res = rc.Get()
957 if not res: print('calcMu '+t+': something went wrong',runNumber)
958 muAv = res.Parameter(0)
959 tag ='Scifi'
960 if t.find('DS')>0: tag='DS'
961 h['muAv'][runNumber][tag] = muAv
962 fout = ROOT.TFile('Mu.root','recreate')
963 for x in h:
964 if x.find('Mu')>0: h[x].Write()
965 fout.Close()
966# update runInfo
967 for runNumber in h['muAv']:
968 r = int(runNumber)
969 if not r in runInfo:
970 print('calcMu: run not in runInfo',r)
971 continue
972 runInfo[r]['muAv'] = {'':h['muAv'][runNumber][''],'Scifi':h['muAv'][runNumber]['Scifi'],'DS':h['muAv'][runNumber]['DS']}
973
974 def FwBw(self,runNumber):
975 options.runNumbers = runNumber
976 options.fillNumbers = self.runInfo[options.runNumbers]['Fillnumber']
977# analyze run for forward / backward tracks per bunch type
978 h = self.h
979 self.F = ROOT.TFile.Open(www+"offline/run"+str(runNumber).zfill(6)+".root")
980 self.B = ROOT.TFile.Open(www+"offline/BunchStructure.root")
981 self.L = ROOT.TFile.Open(www+"offline/Lumi.root")
982 xing = {'all':False,'B1only':False,'B2noB1':False,'noBeam':False}
983
984 for T in ['Txing','TD','T']:
985 h[T] = self.F.Get("daq").Get(T).Clone(T)
986 for X in ['time','timeWt','timeWtDS']:
987 h[X] = h['T'].FindObject(X).Clone(X)
988 for t in ['time','timeWt','timeWtDS','bnr']:
989 for x in xing:
990 if x=='all': continue
991 X = t+x
992 h[X] = h['Txing'].FindObject(X).Clone(X)
993 for x in xing:
994 X = 'trackDir'+x
995 h[X] = h['TD'].FindObject(X).Clone(X)
996 h[X].SetTitle('Track Velocity; 1/v - 1/c [ns/cm]')
997 h[X].SetName('Track Velocity '+x)
998 h['BunchStructure'] = self.B.Get("run"+str(runNumber).zfill(6)).Clone("run"+str(runNumber).zfill(6))
999 h['Lumi'] = self.L.Get("run"+str(runNumber).zfill(6)).Clone("run"+str(runNumber).zfill(6))
1000 for x in h['Lumi'].GetListOfPrimitives():
1001 if x.ClassName() == "TGraph": break
1002 h['LumiT'] = x.Clone()
1003 Lmax = 0
1004 for n in range(h['LumiT'].GetN()):
1005 l = h['LumiT'].GetPointY(n)
1006 if l>Lmax: Lmax = l
1007
1008 for x in ['IP1z','b1z','b2z']:
1009 h[x] = h['BunchStructure'].FindObject(x)
1010 ROOT.gROOT.cd()
1011 ut.bookCanvas(h,'dirRes','',1600,900,2,1)
1012 tc = h['dirRes'].cd(1)
1013 rc = h['trackDirB2noB1'].Fit('gaus','QS')
1014 rc = h['trackDirB1only'].Fit('gaus','QS')
1015 h['trackDirB2noB1'].Draw()
1016 tc.Update()
1017 stats = h['trackDirB2noB1'].FindObject('stats')
1018 stats.SetOptFit(1111111)
1019 stats.SetX1NDC(0.15)
1020 stats.SetY1NDC(0.5)
1021 stats.SetX2NDC(0.52)
1022 stats.SetY2NDC(0.86)
1023 tc.Update()
1024 tc = h['dirRes'].cd(2)
1025 h['trackDirB1only'].Draw()
1026 tc.Update()
1027 stats = h['trackDirB1only'].FindObject('stats')
1028 stats.SetOptFit(1111111)
1029 stats.SetX1NDC(0.15)
1030 stats.SetY1NDC(0.5)
1031 stats.SetX2NDC(0.52)
1032 stats.SetY2NDC(0.86)
1033 tc.Update()
1034 h['dirRes'].Update()
1035 self.myPrint('dirRes','dirRes-'+str(runNumber).zfill(6))
1036
1037# stats:
1038 norm = {}
1039 for x in xing:
1040 if x=='all': continue
1041 y = h['bnr'+x]
1042 norm[x]=0
1043 for i in range(y.GetNbinsX()+1):
1044 if y.GetBinContent(i)>0: norm[x]+=1
1045 self.stats = {}
1047 bunches = {}
1048 for x in ['IP1z','b1z','b2z']:
1049 bunches[x]=0
1050 for b in range(h[x].GetNbinsX()):
1051 if h[x].GetBinContent(b+1)>0: bunches[x]+=1
1052 print(bunches)
1053 print(norm)
1054 txt = {'timeWt':'scifi tracks','timeWtDS':'DS tracks'}
1055 self.frac = {}
1056 for t in txt:
1057 self.stats[t]={}
1058 self.statsPerBunch[t] = {}
1059 for x in xing:
1060 if x=='all': continue
1061 self.stats[t][x] = h[t+x].GetSumOfWeights()
1062 c = 'noBeam'
1063 self.statsPerBunch[t][c] = self.stats[t][c]/norm[c]
1064 c = 'B2noB1'
1065 self.statsPerBunch[t][c] = (self.stats[t][c] - self.statsPerBunch[t]['noBeam']*norm[c])/norm[c]
1066 c = 'B1only'
1067 self.statsPerBunch[t][c] = (self.stats[t][c] - self.statsPerBunch[t]['noBeam']*norm[c])/norm[c]
1068 print('events with '+t+'/bunch, noBeam subtracted for run ',runNumber)
1069 print('B1 = %5.4F'%(self.statsPerBunch[t]['B1only']))
1070 print('B2 = %5.4F'%(self.statsPerBunch[t]['B2noB1']))
1071 NIP1 = h[t].GetSumOfWeights() - self.stats[t]['noBeam'] - self.stats[t]['B1only'] - self.stats[t]['B2noB1']
1072 self.frac[t] = {}
1073 for c in xing:
1074 if c=='all': continue
1075 self.frac[t][c] = self.statsPerBunch[t][c]*bunches['IP1z']/NIP1
1076
1077 print(txt[t]+' expected B1, B2, noBeam events in IP1: %5.2F%% %5.2F%% %5.2F%%'%(
1078 self.frac[t]['B1only']*100,self.frac[t]['B2noB1']*100,self.frac[t]['noBeam']*100))
1079
1080# plot of B2, B1 and no beam with lumi
1081 self.addBunchCurrent(options.fillNumbers,b=2)
1082 self.addBunchCurrent(options.fillNumbers,b=1)
1083 rebin = {'B2noB1':100,'B1only':1000,'noBeam':1000}
1084 rescale = {'B2noB1':[100,20],'B1only':[5,5],'noBeam':[100,100]}
1085 for B in ['B2noB1','B1only','noBeam']:
1086 ut.bookCanvas(h,B,'',1200,900,1,1)
1087 h[B].cd()
1088 for j in ['time'+B,'timeWtDS'+B,'timeWt'+B]:
1089 h[j+'_100'] = h[j].Clone()
1090 h[j+'_100'].Rebin(rebin[B])
1091 h[j+'_100'].Scale(1/(rebin[B]*norm[B]))
1092
1093 h['time'+B+'_100'].SetTitle(';[s];events/s per bunch')
1094 h['time'+B+'_100'].Draw()
1095 h['timeWt'+B+'_100'].Scale(rescale[B][0])
1096 h['timeWtDS'+B+'_100'].Scale(rescale[B][1])
1097 h['timeWtDS'+B+'_100'].Draw('same')
1098 h['timeWt'+B+'_100'].Draw('same')
1099 h[B].Update()
1100# work with second axis
1101 rightmax = 1.1*Lmax/1000
1102 scale = ROOT.gPad.GetUymax()/rightmax
1103 h['LumiT'+B] = h['LumiT'].Clone('LumiT'+B)
1104 h['LumiT'+B].Scale(scale/1000)
1105 h['LumiT'+B].Draw('same')
1106 h['ax1'+B] = ROOT.TGaxis(ROOT.gPad.GetUxmax(), ROOT.gPad.GetUymin(),
1107 ROOT.gPad.GetUxmax(), ROOT.gPad.GetUymax(),
1108 0, rightmax, 510, "+L")
1109 h['ax1'+B].SetTitle('L [Hz/nb] ')
1110 h['ax1'+B].SetTextFont(42)
1111 h['ax1'+B].SetLabelFont(42)
1112 h['ax1'+B].SetTextColor(ROOT.kMagenta)
1113 h['ax1'+B].Draw()
1114 h['l'+B]=ROOT.TLegend(0.44,0.86,0.88,0.98)
1115 h['l'+B].AddEntry(h['timeB2noB1_100'],'triggered event rate ',"PL")
1116 h['l'+B].AddEntry(h['timeWtB2noB1_100'],'event rate#times'+str(rescale[B][0])+' with Scifi tracks',"PL")
1117 h['l'+B].AddEntry(h['timeWtDSB2noB1_100'],'event rate#times'+str(rescale[B][1])+' with DS tracks',"PL")
1118 h['l'+B].AddEntry(h['LumiT'],'IP1 instanteous luminosity',"PL")
1119 h['l'+B].Draw()
1120 h[B].Update()
1121 self.myPrint(B,B+'-'+str(runNumber).zfill(6))
1122 beam = -1
1123 if B=='B2noB1': beam=2
1124 elif B=='B1only': beam=1
1125 if beam>0:
1126 cur = 'b'+str(beam)+'A'
1127 scale = 0.8*h['time'+B+'_100'].GetMaximum()/self.beamCurrent[cur][0]
1128 self.beamCurrent[cur][1].Scale(scale)
1129 self.beamCurrent[cur][1].SetLineColor(ROOT.kRed)
1130 self.beamCurrent[cur][1].SetLineWidth=2
1131 self.beamCurrent[cur][1].Draw('same')
1132 cur = 'b'+str(beam)+'B'
1133 self.beamCurrent[cur][1].Scale(scale)
1134 self.beamCurrent[cur][1].SetLineColor(ROOT.kBlue)
1135 self.beamCurrent[cur][1].SetLineWidth=2
1136 self.beamCurrent[cur][1].Draw('same')
1137 self.myPrint(B,B+'-'+str(runNumber).zfill(6)+'_withCurrent')
1138
1139# plot track direction / velocity for different xings normalized
1140 for B in ['B2noB1','B1only','noBeam']:
1141 h['trackDir'+B+'_norm'] = h['trackDir'+B].Clone('trackDir'+B+'_norm')
1142 h['trackDirB2noB1_norm'].Scale(bunches['b1z']/norm['B2noB1'])
1143 h['trackDirB1only_norm'].Scale(bunches['b1z']/norm['B1only'])
1144 ut.bookCanvas(h,'V','',1200,900,1,1)
1145 tc = h['V'].cd()
1146 tc.SetLogy(1)
1147 h['trackDirall'].Draw()
1148 h['trackDirall'].SetStats(0)
1149 self.myPrint('V','trackDirAll-'+str(runNumber).zfill(6))
1150 h['trackDirB2noB1_norm'].SetLineColor(ROOT.kCyan)
1151 h['trackDirB1only_norm'].SetLineColor(ROOT.kBlue)
1152 h['trackDirnoBeam_norm'].SetLineColor(ROOT.kGreen)
1153 h['trackDirB2noB1_norm'].Draw('Histsame')
1154 h['trackDirB1only_norm'].Draw('Histsame')
1155 # h['trackDirnoBeam_norm'].Draw('same')
1156 self.myPrint('V','trackDirXing-'+str(runNumber).zfill(6))
1157# plot track angles and positions
1158 ut.bookCanvas(h,'2dslopes','',1200,1200,1,1)
1159 ut.bookCanvas(h,'1dslopes','',1200,900,1,1)
1160 for B in ['B2noB1','B1only','noBeam']:
1161 for x in ['scifi-trackDir','scifi-TtrackPos']:
1162 T = x+B
1163 tmp = self.F.Get("scifi").Get(T)
1164 if not tmp:
1165 tmp = self.F.Get("scifi").Get(B+'/'+T)
1166 h[T] = tmp.Clone(T)
1167 if x.find('Dir')>0:
1168 for y in ['scifi-trackSlopesXL','slopeXL','slopeYL']:
1169 h[y+B] = self.h[T].FindObject(y+B).Clone(y+B)
1170 h[y+B].SetStats(0)
1171 if y.find('track')>0:
1172 cv = '2dslopes'
1173 h[cv].cd()
1174 txt = y+B
1175 h[y+B].Draw('colz')
1176 else:
1177 txt = 'scifi-'+y+B
1178 cv = '1dslopes'
1179 h[cv].cd()
1180 h[y+B].Rebin(4)
1181 h[y+B].Draw()
1182 self.myPrint(cv,txt+'-'+str(runNumber).zfill(6))
1183 elif x.find('Pos')>0:
1184 cv = '2dslopes'
1185 h[cv].cd()
1186 y = 'scifi-trackPos'
1187 h[y+B] = self.h[T].FindObject(y+B).Clone(y+B)
1188 h[y+B].SetStats(0)
1189 h[y+B].Draw('colz')
1190 self.myPrint(cv,y+B+'-'+str(runNumber).zfill(6))
1191
1192 for x in ['muonDSTracks','mufi-TtrackPos']:
1193 T = x+B
1194 tmp = self.F.Get("mufilter").Get(T)
1195 if not tmp:
1196 tmp = self.F.Get("mufilter").Get(B+'/'+T)
1197 h[T] = tmp.Clone(T)
1198 if x.find('DSTrack')>0:
1199 for y in ['mufi-slopes','slopeX','slopeY']:
1200 h[y+B] = self.h[T].FindObject(y+B).Clone(y+B)
1201 h[y+B].SetStats(0)
1202 if not y.find('mufi')<0:
1203 cv = '2dslopes'
1204 h[cv].cd()
1205 txt = y+B
1206 h[y+B].Draw('colz')
1207 else:
1208 cv = '1dslopes'
1209 h[cv].cd()
1210 txt = 'mufi-'+y+B
1211 h[y+B].Draw()
1212 self.myPrint(cv,txt+'-'+str(runNumber).zfill(6))
1213 elif x.find('Pos')>0:
1214 y = 'mufi-trackPos'
1215 cv = '2dslopes'
1216 h[cv].cd()
1217 h[y+B] = self.h[T].FindObject(y+B).Clone(y+B)
1218 h[y+B].SetStats(0)
1219 h[y+B].Draw('colz')
1220 self.myPrint(cv,y+B+'-'+str(runNumber).zfill(6))
1221
1222
1223 def addBunchCurrent(self,fillNr,b=2):
1224 F = ROOT.TFile.Open('root://eospublic.cern.ch//eos/experiment/sndlhc/nxcals_data/fill_'+str(fillNr).zfill(6)+'.root')
1225 LHC = F.LHC
1226 injection_scheme = LHC.LHC_STATS_LHC_INJECTION_SCHEME
1227 rc = injection_scheme.GetEvent(0)
1228 print(injection_scheme.var)
1229 betastar = LHC.HX_BETASTAR_IP1
1230 rc = betastar.GetEvent(0)
1231 print('beta* = ',betastar.var)
1232 beam = {}
1233 beam['b1A'] = LHC.LHC_BCTFR_A6R4_B1_BEAM_INTENSITY
1234 beam['b2A'] = LHC.LHC_BCTFR_A6R4_B2_BEAM_INTENSITY
1235 beam['b1B'] = LHC.LHC_BCTFR_B6R4_B1_BEAM_INTENSITY
1236 beam['b2B'] = LHC.LHC_BCTFR_B6R4_B2_BEAM_INTENSITY
1237 t0 = self.runInfo[options.runNumbers]['StartTime']
1238 X = 'b'+str(b)+'A'
1239 self.beamCurrent[X] = [0,ROOT.TGraph()]
1240 mx = 0
1241 for e in beam[X]:
1242 self.beamCurrent[X][1].AddPoint(e.unix_timestamp-t0,e.var)
1243 if e.var > mx: mx=e.var
1244 self.beamCurrent[X][0] = mx
1245
1246 X = 'b'+str(b)+'B'
1247 self.beamCurrent[X] = [0,ROOT.TGraph()]
1248 mx = 0
1249 for e in beam[X]:
1250 self.beamCurrent[X][1].AddPoint(e.unix_timestamp-t0,e.var)
1251 if e.var > mx: mx=e.var
1252 self.beamCurrent[X][0] = mx
1253
1254# work with second axis
1255 if 1<0:
1256 rightmax = 1.1*Lmax/1000
1257 scale = ROOT.gPad.GetUymax()/rightmax
1258 h['LumiT'+B] = h['LumiT'].Clone('LumiT'+B)
1259 h['LumiT'+B].Scale(scale/1000)
1260 h['LumiT'+B].Draw('same')
1261 h['ax1'+B] = ROOT.TGaxis(ROOT.gPad.GetUxmax(), ROOT.gPad.GetUymin(),
1262 ROOT.gPad.GetUxmax(), ROOT.gPad.GetUymax(),
1263 0, rightmax, 510, "+L")
1264 h['ax1'+B].SetTitle('L [Hz/nb] ')
1265 h['ax1'+B].SetTextFont(42)
1266 h['ax1'+B].SetLabelFont(42)
1267 h['ax1'+B].SetTextColor(ROOT.kMagenta)
1268 h['ax1'+B].Draw()
1269 h['l'+B]=ROOT.TLegend(0.44,0.86,0.91,0.98)
1270 h['l'+B].AddEntry(h['timeB2noB1_100'],'triggered event rate ',"PL")
1271 h['l'+B].AddEntry(h['timeWtB2noB1_100'],'event rate#times'+str(rescale[B][0])+' with Scifi tracks',"PL")
1272 h['l'+B].AddEntry(h['timeWtDSB2noB1_100'],'event rate#times'+str(rescale[B][1])+' with DS tracks',"PL")
1273 h['l'+B].AddEntry(h['LumiT'],'IP1 instanteous luminosity',"PL")
1274 h['l'+B].Draw()
1275 h[B].Update()
1276 self.myPrint(B,B+'-'+str(runNumber).zfill(6))
1277
1278
1279
1280 def merge(self):
1281 h = self.h
1282 for fname in os.listdir():
1283 if fname.find('FS')==0 and fname.find('.root')>0 and fname.find('dict')<0:
1284 rname = fname.split('-')[1].split('.')[0]
1285 F = ROOT.TFile(fname)
1286 h[rname] = F.c1.Clone(rname)
1287 h[rname].SetName(rname)
1288 h[rname].SetTitle(rname)
1289 F = ROOT.TFile('BunchStructure.root','recreate')
1290 keys = list(h.keys())
1291 keys.sort(reverse=True)
1292 for r in keys:
1293 if r.find('run')==0:
1294 if int(r.split('run')[1])<options.rmin: continue
1295 h[r].Write()
1296 F.Close()
1297
1298 def mergeLumi(self):
1299 h = self.h
1300 Llist = []
1301 for fname in os.listdir():
1302 if fname.find('Lumi-run')==0 and fname.find('.root')>0:
1303 rname = fname.split('-')[1].split('.')[0]
1304 F = ROOT.TFile(fname)
1305 if not F.Get('c1'):
1306 print('error file',fname)
1307 continue
1308 h[rname] = F.c1.Clone(rname)
1309 h[rname].SetName(rname)
1310 h[rname].SetTitle(rname)
1311 Llist.append(rname)
1312 F = ROOT.TFile('Lumi.root','recreate')
1313 Llist.sort(reverse=True)
1314 for r in Llist:
1315 print('write 0',r,h[r])
1316 if r.find('run')==0: h[r].Write()
1317 print('write 1',r,h[r])
1318 F.Close()
1319
1320 def lhcNumbering(self):
1321 h = self.h
1322 Nbunches = h['bnr'].GetNbinsX()
1323 F = ROOT.TFile('BunchStructure.root')
1324 p = open("FSdict.pkl",'rb')
1325 self.FSdict = pickle.load(p)
1326 for k in F.GetListOfKeys():
1327 rname = k.GetName()
1328 newname = 'lhc-'+rname
1329 h[newname] = F.Get(rname).Clone(newname)
1330 h[newname].SetName(rname)
1331 h[newname].SetTitle(rname)
1332 F.Close()
1333 F = ROOT.TFile('BunchStructureLHC.root','recreate')
1334 for newname in h:
1335 if not newname.find('lhc-') == 0:continue
1336 X = {}
1337 for p in h[newname].GetListOfPrimitives():
1338 X[p.GetName()] = p
1339 a = X['b1z'].GetTitle().replace(' ','')
1340 X['bnr'].SetLineColor(ROOT.kBlack)
1341 runnr = int(a.split('run')[1].split('fill')[0])
1342 phaseShift1 = self.FSdict[runnr]['phaseShift1']
1343 phaseShift2 = self.FSdict[runnr]['phaseShift2']
1344 for hname in X:
1345 histo=X[hname]
1346 if not hasattr(histo,'GetBinContent'): continue
1347 tmp = {}
1348 for i in range(Nbunches):
1349 tmp[i+1] = histo.GetBinContent(i+1)
1350 for i in range(Nbunches):
1351 newBin = (Nbunches - phaseShift1 + i)%Nbunches
1352 histo.SetBinContent(newBin,tmp[i+1])
1353 h[newname].Update()
1354 h[newname].Write()
1355
1357# check for partitions
1358 runNr = str(r).zfill(6)
1359 partitions = []
1360 conv_data_path = options.rawData.replace("raw_data", "convertedData")
1361 eventChain = ROOT.TChain('rawConv')
1362 eventChain.Add(os.environ['EOSSHIP']+conv_data_path+'run_'+runNr+'/*.root')
1363 nEvents = eventChain.GetEntries()
1364# make the plot
1365 # figure out the number of bunches in the LHC
1366 rc = eventChain.GetEvent(0)
1367 if eventChain.EventHeader.GetAccMode()==12: # ion runs
1368 Nbunches = 1782
1369 div = 8
1370 else: # proton runs
1371 Nbunches = 3564
1372 div = 4
1373 ut.bookHist(self.h,'bnr_from_data','bunch number; LHC bunch number', Nbunches,-0.5,Nbunches-0.5)
1374 # use postscale, same logic as in the monitoring task
1375 postScale = 0
1376 if nEvents>10E6: postScale = 10
1377 if nEvents>100E6: postScale = 100
1378 print('using postScale ',postScale,' for run ',r)
1379 for event in eventChain:
1380 if postScale>0:
1381 if ROOT.gRandom.Rndm()>1./postScale: continue
1382 self.h['bnr_from_data'].Fill(int((event.EventHeader.GetEventTime()%(div*Nbunches))/div+0.5))
1383
1384 return self.h['bnr_from_data']
1385
1386 def getEntriesPerRun(self,r):
1387# check for partitions
1388 runNr = str(r).zfill(6)
1389 partitions = []
1390 dirlist = str( subprocess.check_output("xrdfs "+os.environ['EOSSHIP']+" ls "+options.convpath+"run_"+runNr,shell=True) )
1391 for x in dirlist.split('\\n'):
1392 ix = x.find('sndsw_raw-')
1393 if ix<0: continue
1394 partitions.append(x[ix:])
1395 eventChain = ROOT.TChain('rawConv')
1396 for p in partitions:
1397 eventChain.Add(os.environ['EOSSHIP']+options.convpath+'run_'+runNr+'/'+p)
1398 return eventChain.GetEntries(),partitions
1399
1400 def getTotalStat(self):
1401 L = 0
1402 N = 0
1403 for r in self.runInfo:
1404 L+=self.runInfo[r]['lumiAtIP1withSNDLHC']
1405 N+=self.runInfo[r]['Entries']
1406
1407 def makeLatex(self):
1408 # make latex code
1409 # filling schemes: FS-run004626.pdf and Lumi-run004626.pdf
1410 L,N= 0,0
1411 for r in self.runInfo:
1412 if r>options.rmin:
1413 L+=self.runInfo[r]['lumiAtIP1withSNDLHC']
1414 N+=self.runInfo[r]['Entries']
1415 lines = []
1416 lines.append("\documentclass{beamer}")
1417 lines.append("\mode<presentation>")
1418 lines.append("{\\usetheme{Singapore}}")
1419 lines.append("\\usepackage{graphicx}")
1420 lines.append("\\usepackage[space]{grffile}")
1421 lines.append("\\usepackage[english]{babel}")
1422 lines.append("\\usepackage[latin1]{inputenc}")
1423 lines.append("\\usepackage[T1]{fontenc}")
1424 lines.append("\\title[Short Paper Title] % (optional, use only with long paper titles)")
1425 if options.convpath.find('2022')>0:
1426 lines.append("{SND@LHC Run Summary July - November 2022}")
1427 lines.append("\date[Short Occasion] % (optional)")
1428 lines.append("{ 17 November 2022}")
1429 lines.append("\\begin{document}")
1430 lines.append("\\begin{frame}{}")
1431 lines.append("17 November 2022")
1432 lines.append("\\newline ")
1433 lines.append("\\newline ")
1434 lines.append("Run Summary for July - November 2022")
1435 else:
1436 lines.append("{SND@LHC Run Summary March - July 2023}")
1437 lines.append("\date[Short Occasion] % (optional)")
1438 lines.append("{ 16 July 2023}")
1439 lines.append("\\begin{document}")
1440 lines.append("\\begin{frame}{}")
1441 lines.append("16 July 2023")
1442 lines.append("\\newline ")
1443 lines.append("\\newline ")
1444 lines.append("Run Summary for March - July 2023")
1445 nTXT = "$%5.2F\\times 10^9 $"%(N/1E9)
1446 lines.append("\\begin{itemize}")
1447 lines.append("\item total number of events: "+nTXT)
1448 lines.append("\item integrated luminosity (lower limit): $%5.2F\mathrm{fb}^{-1}$"%(L/1E9))
1449 lines.append("\end{itemize}")
1450 lines.append("\\begin{center}")
1451 lines.append("\includegraphics[width = 0.7\\textwidth]{Lumi-time.pdf}")
1452 lines.append("\end{center}")
1453 lines.append("\end{frame}")
1454 lines.append("\\begin{frame}{}")
1455 lines.append("\\begin{center}")
1456 elist = list(emulsionReplacements.values())
1457 elist.sort(reverse=True)
1458 k=0
1459 for emulsionNr in elist:
1460 if k==6:
1461 k=0
1462 lines.append("\end{center}")
1463 lines.append("\end{frame}")
1464 lines.append("\\begin{frame}{}")
1465 lines.append("\\begin{center}")
1466 if "ScifitrackDens"+str(emulsionNr)+".pdf" in os.listdir('.'):
1467 lines.append("\includegraphics[width = 0.4\\textwidth]{ScifitrackDens"+str(emulsionNr)+".pdf}")
1468 k+=1
1469 lines.append("\end{center}")
1470 lines.append("\end{frame}")
1471
1472 R = list(self.runInfo.keys())
1473 R.sort(reverse=True)
1474
1475 lines.append("\\begin{frame}{}")
1476 lines.append("Overview of runs \\newline")
1477 lines.append("\\scriptsize")
1478 lines.append("\\begin{scriptsize}")
1479 lines.append("\\begin{tabular}{lcrrrl}")
1480 lines.append(" Run & Fill & events & mu & Lumi $\mathrm{pb}^{-1}$ & start time \\\\ ")
1481 ilines = 0
1482 for i in range(len(R)):
1483 r = R[i]
1484 if r < options.rmin: continue
1485 lumi = self.runInfo[r]['lumiAtIP1withSNDLHC']/1E6
1486 N = self.runInfo[r]['Entries']
1487 fill = str(self.runInfo[r]['Fillnumber'])
1488 if fill == '': fill = ' -- '
1489 mu = ''
1490 if 'muAv' in self.runInfo[r]: mu = "%4.1F"%(self.runInfo[r]['muAv'][''])
1491 lines.append(" %i & %6s & %10i & %s & $%5.1F$ & %s \\\\"%(r,fill,N,mu,lumi,self.runInfo[r]['StartTimeC']))
1492 ilines+=1
1493 if ilines%19==0:
1494 lines.append("\end{tabular}")
1495 lines.append("\\end{scriptsize}")
1496 lines.append("\end{frame}")
1497 lines.append("\\begin{frame}{}")
1498 lines.append("\\begin{scriptsize}")
1499 lines.append("\\begin{tabular}{lcrrrl}")
1500 lines.append(" Run & Fill & events & mu & Lumi $\mathrm{pb}^{-1}$ & start time \\\\ ")
1501 lines.append("\end{tabular}")
1502 lines.append("\\end{scriptsize}")
1503 lines.append("\end{frame}")
1504
1505#
1506# runs with beam present measured
1507 RwL = []
1508 self.runsWithBeam()
1509 for r in self.runs:
1510 withBeam = False
1511 if r not in self.runInfo: continue
1512 if 'timeWtDS' in self.runs[r]:
1513 if self.runs[r]['timeWtDS']/self.runs[r]['time'] > 0.25: withBeam = True
1514 if self.runInfo[r]['lumiAtIP1withSNDLHC']>0 or withBeam: RwL.append(r)
1515 RwL.sort(reverse=True)
1516 R = RwL
1517 lines.append("\\begin{frame}{}")
1518 k=0
1519 for i in range(len(R)):
1520 print('at run ',R[i])
1521 r = str(R[i]).zfill(6)
1522 if not "FS-run"+r+".pdf" in os.listdir('.'): continue
1523 if "Lumi-run"+r+".pdf" in os.listdir('.'):
1524 lines.append("\includegraphics[width = 0.3\\textwidth]{Lumi-run"+r+".pdf}")
1525 else:
1526 lines.append("\includegraphics[width = 0.3\\textwidth]{noLumi-run"+r+".pdf}")
1527 if "FS-run"+r+".pdf" in os.listdir('.'):
1528 lines.append("\includegraphics[width = 0.3\\textwidth]{FS-run"+r+".pdf}")
1529 if "run"+r+"Lumi-tracks.pdf" in os.listdir('.'):
1530 lines.append("\includegraphics[width = 0.3\\textwidth]{run"+r+"Lumi-tracks.pdf}")
1531 lines.append("\\newline")
1532 if (3*k+3)%12==0:
1533 lines.append("\end{frame}")
1534 lines.append("\\begin{frame}{}")
1535 k+=1
1536 lines.append("\end{frame}")
1537 lines.append(" ")
1538#
1539 lines.append("\end{document}")
1540 outFile = open(os.environ["HOME"]+'/dummy','w')
1541 for l in lines:
1542 rc = outFile.write(l+"\n")
1543 outFile.close()
1544 os.system('cp $HOME/dummy LumiSummary.tex')
1545
1547 time_objA = time.strptime(dateA,'%m-%d,%Y-%H')
1548 time_objB = time.strptime(dateB,'%m-%d,%Y-%H')
1549 TA = calendar.timegm(time_objA)
1550 TB = calendar.timegm(time_objB)
1551 if not 'cL' in h: ut.readHists(h,'Lumi-time.root')
1552 n=0
1553 for x in h['cL'].GetListOfPrimitives():
1554 if x.ClassName() == 'TGraph':
1555 h['Lumi'+str(n)] = x.Clone('Lumi'+str(n))
1556 n+=1
1557 if x.ClassName() == 'TGaxis':
1558 h['IntLumiAxis'] = x.Clone('IntLumiAxis')
1559 d = h['Lumi1'].Eval(TB)-h['Lumi1'].Eval(TA)
1560 scale = h['IntLumiAxis'].GetWmax() / h['IntLumiAxis'].GetY2()
1561 print('Lumi between ',dateA,dateB,' = ',d*scale)
1562
1564 h = self.h
1565 runInfo = self.runInfo
1566 h['LumiT']=ROOT.TGraph()
1567 h['ILumiT']=ROOT.TGraph()
1568 Lint = 0
1569 runList = list(runInfo.keys())
1570 runList.sort()
1571 lmax = 0
1572 for r in runList:
1573 if r<options.rmin: continue
1574 h['LumiT'].AddPoint(runInfo[r]['StartTime'],0)
1575 X = runInfo[r]['lumiAtIP1withSNDLHC']/1E6
1576 h['LumiT'].AddPoint(runInfo[r]['StartTime']+1,X)
1577 if X>lmax : lmax = X
1578 h['LumiT'].AddPoint(runInfo[r]['StartTime']+3600,0)
1579 Lint+=X
1580 h['ILumiT'].AddPoint(runInfo[r]['StartTime'],Lint/1000)
1581 tstart = h['LumiT'].GetPointX(0)
1582 tend = h['LumiT'].GetPointX(h['LumiT'].GetN()-1)
1583 ut.bookHist(h,'LumiTime','; ; pb^{-1}',100,tstart,tend)
1584 ut.bookCanvas(self.h,'cL','',2400,800,1,1)
1585 h['cL'].cd()
1586 h['LumiTime'].GetXaxis().SetTimeFormat("%d-%m")
1587 h['LumiTime'].GetXaxis().SetTimeOffset(0,'gmt')
1588 h['LumiTime'].SetMaximum(lmax*1.2)
1589 h['LumiTime'].SetStats(0)
1590 h['LumiTime'].Draw()
1591 h['cL'].Update()
1592# work with second axis
1593 rightmax = 1.1*Lint/1000.
1594 scale = ROOT.gPad.GetUymax()/rightmax
1595 h['ILumiT'].Scale(scale)
1596 h['ax2'] = ROOT.TGaxis(ROOT.gPad.GetUxmax(), ROOT.gPad.GetUymin(),
1597 ROOT.gPad.GetUxmax(), ROOT.gPad.GetUymax(),
1598 0, rightmax, 510, "+L")
1599 h['ax2'].SetTitle('L Integral; fb^{-1}')
1600 h['ax2'].SetTextFont(42)
1601 h['ax2'].SetLabelFont(42)
1602 h['ax2'].SetTextColor(ROOT.kRed)
1603 h['ax2'].Draw()
1604 h['LumiT'].SetLineColor(ROOT.kBlue)
1605 h['LumiT'].SetLineWidth(2)
1606 h['LumiT'].Draw('Lsame')
1607 h['ILumiT'].SetLineColor(ROOT.kRed)
1608 h['ILumiT'].SetLineWidth(3)
1609 h['ILumiT'].Draw('same')
1610 self.myPrint('cL','Lumi-time')
1611
1612 def LumiIntegral(self,rmin,rmax):
1613 L = 0
1614 for r in self.runInfo:
1615 if r>= rmin and r<=rmax:
1616 L+=self.runInfo[r]['lumiAtIP1withSNDLHC']
1617 print('Lumi for run ',rmin,' to run ',rmax,' = ',L)
1618 def LumiPerFill(self):
1619 lumiPerFill={}
1620 for r in self.runInfo:
1621 lumiPerFill[ self.runInfo[r]['Fillnumber'] ]=self.runInfo[r]['lumiAtIP1']
1622 I = 0
1623 for x in lumiPerFill:
1624 I+=lumiPerFill[x]
1625 print('total:',l)
1626
1627 def runsWithBeam(self):
1628# potential runs with beam selected by looking for daq rate > cosmics
1629 self.listOfRuns = {}
1630 if options.convpath.find('2022')>0:
1631 offline =www + "offline.html"
1632 else:
1633 offline =www + "offline2023.html"
1634 with client.File() as f:
1635 f.open(offline)
1636 status, L = f.read()
1637 Lhtml = L.decode().split('\n')
1638 f.close()
1639 self.runs = {}
1640 for x in Lhtml:
1641 if x.find('#event')<0: continue
1642 ir = x.find("run ")
1643 if ir<0: continue
1644 k = x[ir:].find(" ")
1645 runNumber = int(x[ir+4:ir+4+k+1])
1646 if runNumber<options.rmin: continue
1647 R = ROOT.TFile.Open(www+"offline/run"+str(runNumber).zfill(6)+".root")
1648 bCanvas = R.Get("daq").Get('T')
1649 if not bCanvas:
1650 print('Error with root file',runNumber)
1651 continue
1652 Xt = {'time':None,'timeWtDS':None,'timeWt':None}
1653 self.runs[runNumber] = {}
1654 for x in Xt:
1655 Xt[x] = bCanvas.FindObject(x)
1656 if Xt[x]:
1657 Ntot = Xt[x].GetSumOfWeights()
1658 Ttot = Xt[x].GetBinCenter(Xt[x].GetNbinsX()-1)
1659 self.runs[runNumber][x] = Ntot/Ttot
1660 for r in self.runs:
1661 if 'timeWtDS' in self.runs[r]:
1662 if self.runs[r]['timeWtDS']/self.runs[r]['time'] > 0.25:
1663 if not r in self.runInfo: print('run with beam but no filling scheme:',r,'Fillnumber',self.runs[r])
1664 else:
1665 print('run not uptodate',r)
1666
1667 def tracksPerLumi(self,aRun=False):
1668 h=self.h
1669 pol1 = ROOT.TF1('pol1','[0]+[1]*x',0,1E7)
1670 self.F = ROOT.TFile.Open(options.path+"Lumi.root")
1671 if aRun:
1672 keys = ['run'+str(aRun).zfill(6)]
1673 Frun = ROOT.TFile.Open(www+"offline/"+keys[0]+'.root')
1674 else:
1675 keys = self.F.GetListOfKeys()
1676 ROOT.gROOT.cd()
1677 for k in keys:
1678 if aRun: nm = k
1679 else: nm = k.GetName()
1680 runNumber = int(nm.split('run')[1])
1681 if runNumber < options.rmin: continue
1682 tc = self.F.Get(nm)
1683 for x in tc.GetListOfPrimitives():
1684 if x.GetTitle()=='' and x.GetName()=='': lumi = x
1685 if x.GetTitle().find('Hz/nb')>0: yaxis = x
1686 if not lumi:
1687 print('lumi not available',runNumber)
1688 continue
1689 calib = yaxis.GetWmax()/yaxis.GetY2()
1690 lumi.Scale(calib)
1691 mean = lumi.GetMean(axis=2)
1692 for iMin in range(lumi.GetN()):
1693 tMin = lumi.GetPointX(iMin)
1694 if lumi.GetPointY(iMin)>mean*0.001: break
1695 for i in range(lumi.GetN()-1,0,-1):
1696 tMax = lumi.GetPointX(i)
1697 if lumi.GetPointY(i)/lumi.GetPointY(iMin)>0.001: break
1698 # special cases
1699 if runNumber == 5157: tMin = 4000
1700 if runNumber == 5122: tMin = 8000
1701 if runNumber == 5120: tMin = 4000
1702 if runNumber == 4449: tMin = 5100
1703 if runNumber == 4504:
1704 tMin = 4900
1705 tMax = 15000
1706 for x in ["timeWtDS","timeWt"]:
1707 hx = nm+x+'100'
1708 if aRun:
1709 h[hx] = Frun.Get("daq").T.FindObject(x).Clone(hx)
1710 h[hx].Rebin(100)
1711 h[hx].Scale(0.01)
1712 else:
1713 h[hx] = tc.FindObject(x+'10').Clone(hx)
1714 h[hx].Rebin(10)
1715 h[hx].Scale(0.1)
1716 h[hx].GetYaxis().SetTitle('dt = 100s [evts/nb^{-1}]')
1717 h[hx].SetTitle('DS(cyan) and Scifi(red) tracks run '+str(runNumber))
1718 if tMin < h[hx].GetBinCenter(1): tMin = h[hx].GetBinCenter(3)
1719 if tMax > h[hx].GetBinCenter(h[hx].GetNbinsX()-1): tMax = h[hx].GetBinCenter(h[hx].GetNbinsX()-1)
1720 if h[hx].GetEntries()>0:
1721 for i in range(h[hx].GetNbinsX()):
1722 t1 = h[hx].GetBinLowEdge(i+1)
1723 t2 = t1 + h[hx].GetBinWidth(i+1)
1724 tl = [0,0]
1725 rMean,rK =0,0
1726 if t1>=tMin and t2<=tMax:
1727 rc = lumi.Fit(pol1,'q','',t1,t2)
1728 L = pol1.Integral(t1,t2)/(t2-t1) # s-1 nb-1
1729 if not L>0:
1730 print('error with lumi',runNumber,t1,t2)
1731 else:
1732 tl[0] = h[hx].GetBinContent(i+1)/L
1733 tl[1] = h[hx].GetBinError(i+1)/L
1734 h[hx].SetBinContent(i+1,tl[0])
1735 rMean += tl[0]
1736 rK+=1
1737 h[hx].SetBinError(i+1,tl[1])
1738# look for outliers
1739 rMean = rMean/rK
1740 for i in range(h[hx].GetNbinsX()):
1741 if h[hx].GetBinContent(i+1)< 0.4*rMean:
1742 h[hx].SetBinContent(i+1,0)
1743 h[hx].SetBinError(i+1,0)
1744 tc = h['c1'].cd()
1745 h[nm+'timeWtDS100'].SetMaximum(140)
1746 h[nm+'timeWtDS100'].SetMinimum(0.0)
1747 h[nm+'timeWtDS100'].Draw()
1748 textInfo = ROOT.TLatex()
1749 if h[nm+'timeWt100'].GetEntries()>0:
1750 h[nm+'timeWt100'].Draw('same')
1751 rc = h[nm+'timeWt100'].Fit('pol1','qS','',tMin,tMax)
1752 res = rc.Get()
1753 if res:
1754 chi2 = res.Chi2()/(tMax-tMin)
1755 fun = h[nm+'timeWt100'].GetFunction('pol1')
1756 if chi2<10 and abs(fun.GetParameter(1)*3600)<5 and fun.GetParameter(0) < 500:
1757 b,m = fun.GetParameter(0),fun.GetParameter(1)
1758 meanV = m*(tMin+tMax)/2+b
1759 txtScifi = "Scifi tracks per nb^{-1} mean: %5.1F slope: %5.2F per hour %5.1F "%(meanV,fun.GetParameter(1)*3600,chi2)
1760 textInfo.DrawLatexNDC(0.2, 0.85,txtScifi)
1761 else:
1762 print('Fit failed',nm+'timeWt100',tMin,tMax,chi2,abs(fun.GetParameter(1)*3600),fun.GetParameter(0))
1763 if h[nm+'timeWt100'].GetFunction('pol1'): h[nm+'timeWt100'].GetFunction('pol1').Delete()
1764
1765 rc = h[nm+'timeWtDS100'].Fit('pol1','Sq','',tMin,tMax)
1766 res = rc.Get()
1767 if res:
1768 chi2 = res.Chi2()/(tMax-tMin)
1769 fun = h[nm+'timeWtDS100'].GetFunction('pol1')
1770 if chi2<20 and abs(fun.GetParameter(1)*3600)<15 and fun.GetParameter(0) < 500:
1771 b,m = fun.GetParameter(0),fun.GetParameter(1)
1772 meanV = m*(tMin+tMax)/2+b
1773 txtDS = " DS tracks per nb^{-1} mean: %5.1F slope: %5.1F per hour %5.1F "%(meanV,m*3600,chi2)
1774 textInfo.DrawLatexNDC(0.2, 0.8,txtDS)
1775 else:
1776 print('Fit failed',nm+'timeWtDS100',tMin,tMax,chi2,abs(fun.GetParameter(1)*3600),fun.GetParameter(0))
1777 if h[nm+'timeWtDS100'].GetFunction('pol1'): h[nm+'timeWtDS100'].GetFunction('pol1').Delete()
1778 tc.Update()
1779 h['c1'].Print(nm+'Lumi-tracks.root')
1780 h['c1'].Print(nm+'Lumi-tracks.pdf')
1781# integral of tracks per cm2
1782 self.R = ROOT.TFile.Open(www+"offline/run"+str(runNumber).zfill(6)+".root")
1783 Nevts = self.R.Get("daq").Get('T').FindObject('time').GetEntries()
1784 postScale = self.runInfo[runNumber]['Entries']/Nevts
1785 T = 'scifi-TtrackPos'
1786 h[T] = self.R.Get("scifi").Get(T).FindObject('scifi-trackPosBeam').Clone(T)
1787 T = 'mufi-TtrackPos'
1788 h[T] = self.R.Get("mufilter").Get(T).FindObject('mufi-trackPosBeam').Clone(T)
1789 ROOT.gROOT.cd()
1790 rList = list(emulsionReplacements.keys())
1791 rList.sort(reverse=True)
1792 for runNr in rList:
1793 if runNumber >= runNr:
1794 emulsionNr = emulsionReplacements[runNr]
1795 break
1796 for k in [ 'scifi-TtrackPos','mufi-TtrackPos']:
1797 if not 'I'+k+str(emulsionNr) in h:
1798 h['I'+k+str(emulsionNr)] = h[k].Clone('I'+k)
1799 h['I'+k+str(emulsionNr)].Scale(postScale)
1800 if k=='scifi-TtrackPos': h['emulsionILumi'+str(emulsionNr)] = self.runInfo[runNumber]['lumiAtIP1withSNDLHC']
1801 else:
1802 h['I'+k+str(emulsionNr)].Add(h[k],postScale)
1803 if k=='scifi-TtrackPos': h['emulsionILumi'+str(emulsionNr)] += self.runInfo[runNumber]['lumiAtIP1withSNDLHC']
1804
1805# merge
1806 for k in self.F.GetListOfKeys():
1807 nm = k.GetName()
1808 tnm = 't'+nm
1809 if not nm+'Lumi-tracks.root' in os.listdir('.'): continue
1810 X = ROOT.TFile(nm+'Lumi-tracks.root')
1811 ROOT.gROOT.cd()
1812 if not X.Get('c1'):
1813 print('problem with',k,nm)
1814 h[tnm] = X.c1.Clone(tnm)
1815 h[tnm].SetName(tnm)
1816 h[tnm].SetTitle(tnm)
1817 F = ROOT.TFile('Lumi-tracks.root','recreate')
1818 keys = list(h.keys())
1819 keys.sort(reverse=True)
1820 for r in keys:
1821 if r.find('trun')==0: h[r].Write()
1822#
1823 elist = list(emulsionReplacements.values())
1824 for emulsionNr in elist:
1825 if not 'Iscifi-TtrackPos'+str(emulsionNr) in h: continue
1826 ut.bookCanvas(h,'ScifitrackDens'+str(emulsionNr),'',900,900,1,1)
1827 tc = h['ScifitrackDens'+str(emulsionNr)].cd()
1828 histo = h['Iscifi-TtrackPos'+str(emulsionNr)]
1829 histo.SetTitle(histo.GetTitle()+' emulsion Nr '+str(emulsionNr)+" #int L=%5.2Ffb^{-1} "%(h['emulsionILumi'+str(emulsionNr)]/1E9))
1830 histo.GetZaxis().SetTitle(' N/cm^{2} ')
1831 histo.Draw('colz')
1832 ROOT.gPad.SetRightMargin(0.13)
1833 zaxis = histo.GetZaxis()
1834 zaxis.SetMaxDigits(3)
1835 zaxis.SetTitleOffset(1.4)
1836 tc.Update()
1837 pal = histo.FindObject('palette')
1838 pal.SetX1NDC(0.86)
1839 pal.SetX2NDC(0.89)
1840 tc.Update()
1841 stats = histo.FindObject('stats')
1842 stats.SetX1NDC(0.11)
1843 stats.SetY1NDC(0.57)
1844 stats.SetX2NDC(0.31)
1845 stats.SetY2NDC(0.81)
1846 tc.Update()
1847 ut.bookCanvas(h,'MufitrackDens'+str(emulsionNr),'',900,900,1,1)
1848 tc = h['MufitrackDens'+str(emulsionNr)].cd()
1849 histo = h['Imufi-TtrackPos'+str(emulsionNr)]
1850 histo.SetTitle(histo.GetTitle()+' emulsion Nr '+str(emulsionNr)+" #int Ldt=%5.2Ffb^{-1} "%(h['emulsionILumi'+str(emulsionNr)]/1E9))
1851 histo.GetZaxis().SetTitle(' N/cm^{2} ')
1852 histo.Draw('colz')
1853 ROOT.gPad.SetRightMargin(0.13)
1854 zaxis = histo.GetZaxis()
1855 zaxis.SetMaxDigits(3)
1856 zaxis.SetTitleOffset(1.4)
1857 pal = histo.FindObject('palette')
1858 if not pal:
1859 print('error',emulsionNr)
1860 return
1861 pal.SetX1NDC(0.86)
1862 pal.SetX2NDC(0.89)
1863 tc.Update()
1864 stats = histo.FindObject('stats')
1865 stats.SetX1NDC(0.11)
1866 stats.SetY1NDC(0.57)
1867 stats.SetX2NDC(0.31)
1868 stats.SetY2NDC(0.81)
1869 tc.Update()
1870 h['ScifitrackDens'+str(emulsionNr)].Print(options.path+'ScifitrackDens'+str(emulsionNr)+'.pdf')
1871 h['MufitrackDens'+str(emulsionNr)].Print(options.path+'MufitrackDens'+str(emulsionNr)+'.pdf')
1872
1873 def fillStats(self,runNr):
1874 fsdict = self.FSdict[runNr]
1875 Nbunches = self.h['bnr'].GetNbinsX()
1876 self.stats = {'B2noB1':0,'B1only':0,'noBeam':0,'IP1':0,'B2andB1':0,'B1':0,'B2':0}
1877 for bunchNumber in range(Nbunches):
1878 nb1 = (Nbunches + bunchNumber - fsdict['phaseShift1'])%Nbunches
1879 nb2 = (Nbunches + bunchNumber - fsdict['phaseShift1']- fsdict['phaseShift2'])%Nbunches
1880 b1 = nb1 in fsdict['B1']
1881 b2 = nb2 in fsdict['B2']
1882 IP1 = False
1883 IP2 = False
1884 if b1:
1885 IP1 = fsdict['B1'][nb1]['IP1']
1886 if b2:
1887 IP2 = fsdict['B2'][nb2]['IP2']
1888 if IP1: self.stats['IP1']+=1
1889 if b1: self.stats['B1']+=1
1890 if b2: self.stats['B2']+=1
1891 if b1 and not IP1 and not b2: self.stats['B1only']+=1
1892 if b2 and not b1: self.stats['B2noB1']+=1
1893 if not b1 and not b2: self.stats['noBeam']+=1
1894 if b1 and b2 and not IP1: self.stats['B2andB1']+=1
1895
1896 def hitMapsNormalized(self,runNumber,Q12MC=False):
1897 marker = {'B2noB1':22,'B1only':21,'noBeam':20}
1898 colors = {'B2noB1':ROOT.kRed,'B1only':ROOT.kOrange,'noBeam':ROOT.kGreen}
1899 h=self.h
1900 Nbunches = self.h['bnr'].GetNbinsX()
1901 self.fillStats(runNumber)
1902 stats = self.stats
1903 norm = {'B2noB1':stats['B2'],'B1only':stats['B1'],'noBeam':Nbunches}
1904 R = ROOT.TFile.Open(www+"offline/run"+str(runNumber).zfill(6)+".root")
1905 ROOT.gROOT.cd()
1906 for hitmaps in ['mufi-barmapsVeto','mufi-barmapsUS','mufi-barmapsDS']:
1907 h[hitmaps] = R.Get("mufilter").Get(hitmaps).Clone(hitmaps)
1908 for B in marker:
1909 tc = R.Get("mufilter").Get(B+'/'+hitmaps+B).Clone(hitmaps+B)
1910 for pad in tc.GetListOfPrimitives():
1911 for plane in pad.GetListOfPrimitives():
1912 if plane.ClassName().find('TH')==0:
1913 hname = plane.GetName()
1914 h[hname] = plane.Clone(hname)
1915 h[hname].Scale(norm[B]/stats[B])
1916 h[hname].SetStats(0)
1917 h[hname].SetMarkerStyle(marker[B])
1918 h[hname].SetLineColor(colors[B])
1919 h[hname].SetMarkerColor(h[hname].GetLineColor())
1920 h[hitmaps].Draw()
1921 tmp = h[hitmaps].Clone('tmp')
1922 j = 0
1923 for pad in tmp.GetListOfPrimitives():
1924 j+=1
1925 tc = h[hitmaps].cd(j)
1926 # tc.SetLogy(1)
1927 for plane in pad.GetListOfPrimitives():
1928 if plane.ClassName().find('TH')==0:
1929 hname = plane.GetName()
1930 plane.SetMinimum(h[hname+'noBeam'].GetMinimum())
1931 plane.SetLineColor(ROOT.kBlue)
1932 plane.SetStats(0)
1933 plane.Draw()
1934 for B in marker:
1935 h[hname+B].Draw('same')
1936 self.myPrint(hitmaps,hitmaps+'-'+str(runNumber).zfill(6))
1937 hitmaps = 'scifi-hitmaps'
1938 h[hitmaps] = R.Get("scifi").Get(hitmaps).Clone(hitmaps)
1939 for B in marker:
1940 tc = R.Get("scifi").Get(B+'/'+hitmaps+B).Clone(hitmaps+B)
1941 for pad in tc.GetListOfPrimitives():
1942 for plane in pad.GetListOfPrimitives():
1943 if plane.ClassName().find('TH')==0:
1944 hname = plane.GetName()
1945 h[hname] = plane.Clone(hname)
1946 h[hname].Scale(Nbunches/stats[B])
1947 h[hname].SetStats(0)
1948 h[hname].SetMarkerStyle(marker[B])
1949 h[hname].SetLineColor(colors[B])
1950 h[hname].SetMarkerColor(h[hname].GetLineColor())
1951 h[hname].SetMarkerSize(0.5)
1952 h[hitmaps].Draw()
1953 tmp = h[hitmaps].Clone('tmp')
1954 j = 0
1955 for pad in tmp.GetListOfPrimitives():
1956 j+=1
1957 tc = h[hitmaps].cd(j)
1958 for plane in pad.GetListOfPrimitives():
1959 if plane.ClassName().find('TH')==0:
1960 hname = plane.GetName()
1961 plane.SetMinimum(h[hname+'noBeam'].GetMinimum())
1962 plane.SetLineColor(ROOT.kBlue)
1963 plane.SetStats(0)
1964 plane.Draw()
1965 for B in marker:
1966 h[hname+B].Draw('same')
1967 self.myPrint(hitmaps,hitmaps+'-'+str(runNumber).zfill(6))
1968 if Q12MC:
1969 fmc = ROOT.TFile(Q12MC)
1970 ROOT.gROOT.cd()
1971 for hitmaps in ['mufi-barmapsVeto','mufi-barmapsUS','mufi-barmapsDS']:
1972 tc = fmc.Get("mufilter").Get(hitmaps).Clone(hitmaps+'MC')
1973 for pad in tc.GetListOfPrimitives():
1974 for plane in pad.GetListOfPrimitives():
1975 if plane.ClassName().find('TH')==0:
1976 hname = plane.GetName()+'MC'
1977 hnameB2 = plane.GetName()+'B2noB1'
1978 h[hname] = plane.Clone(hname)
1979 h[hname].Scale(h[hnameB2].Integral(1,20)/h[hname].Integral(1,20))
1980 h[hname].SetStats(0)
1981 h[hname].SetMarkerStyle(ROOT.kMagenta)
1982 h[hname].SetLineColor(ROOT.kMagenta)
1983 h[hname].SetMarkerColor(h[hname].GetLineColor())
1984 h[hitmaps].Draw()
1985 tmp = h[hitmaps].Clone('tmp')
1986 j = 0
1987 for pad in tmp.GetListOfPrimitives():
1988 j+=1
1989 tc = h[hitmaps].cd(j)
1990 for plane in pad.GetListOfPrimitives():
1991 if plane.ClassName().find('TH')==0:
1992 hname = plane.GetName()
1993 if len(hname)>11: continue
1994 h[hname+'B2noB1'].SetMaximum( 1.1*max(h[hname+'B2noB1'].GetMaximum(),h[hname+'MC'].GetMaximum()))
1995 h[hname+'B2noB1'].SetMinimum(0)
1996 h[hname+'B2noB1'].Draw()
1997 h[hname+'MC'].Draw('sameHist')
1998 self.myPrint(hitmaps,hitmaps+'-Q12MC')
1999
2000
2001 def storeDict(self,dictPtr,dictName,outFileName):
2002 fp = ROOT.TFile.Open(outFileName+'.root','recreate')
2003 pkl = Pickler(fp)
2004 pkl.dump(dictPtr,dictName)
2005 fp.Close()
2006 fp = open(outFileName+'.pkl','wb')
2007 pickle.dump(dictPtr,fp)
2008 fp.close()
2009
2010 def checkSynch(self):
2011 for r in self.FSdict:
2012 if not r in self.runInfo:
2013 print('run does not exist in runInfo',r)
2014 elif not self.runInfo[r]['phaseShift1'] == self.FSdict[r]['phaseShift1']:
2015 print(r,self.runInfo[r]['phaseShift1'], self.FSdict[r]['phaseShift1'])
2016
2017 def modifyFSdict(self,shift=-1):
2018 # shift = -1 for converting 2022 to reproc2022, otherwise 0
2019 fg = ROOT.TFile.Open(options.path+'FSdict.root')
2020 pkl = Unpickler(fg)
2021 self.FSdict = pkl.load('FSdict')
2022 fg.Close()
2023 for r in self.FSdict:
2024 self.FSdict[r]['phaseShift1'] = self.FSdict[r]['phaseShift1'] - shift
2025 self.storeDict(self.FSdict,'FSdict','FSdict')
2026 fg = ROOT.TFile.Open(options.path+'runInfo.root')
2027 pkl = Unpickler(fg)
2028 self.runInfo = pkl.load('runInfo')
2029 fg.Close()
2030 for r in self.runInfo:
2031 self.runInfo[r]['phaseShift1'] = self.FSdict[r]['phaseShift1']
2032 self.storeDict(self.runInfo,'runInfo','RunInfodict')
2033
2034if __name__ == '__main__':
2035
2036 from argparse import ArgumentParser
2037 parser = ArgumentParser()
2038
2039 parser.add_argument("-r", "--runNumbers", dest="runNumbers", help="list of run numbers",required=False, type=str,default="")
2040 parser.add_argument("-F", "--fillNumbers", dest="fillNumbers", help="corresponding fill number",type=str,required=False,default="")
2041 parser.add_argument("-c", "--command", dest="command", help="command", default=None)
2042 parser.add_argument("-p", dest="path", help="path to filling scheme", default="./")
2043 parser.add_argument("-L", dest="lumiversion", help="offline or online lumi from ATLAS", default="offline")
2044 parser.add_argument("-ip2", dest="withIP2", help="with IP2",default=True)
2045 parser.add_argument("-raw", dest="rawData", help="path to rawData",default="/eos/experiment/sndlhc/raw_data/physics/2023") # before "/eos/experiment/sndlhc/raw_data/commissioning/TI18/data"
2046 parser.add_argument("-www", dest="www", help="path to offline folder",default=os.environ['EOSSHIP']+"/eos/experiment/sndlhc/www/")
2047 parser.add_argument("-nMin", dest="nMin", help="min entries for a run",default=100000)
2048 parser.add_argument("--batch", help="run in batch mode",default=False,action='store_true')
2049
2050 options = parser.parse_args()
2051 if options.batch:
2052 ROOT.gROOT.SetBatch(True)
2053 www = options.www
2054 if options.rawData.find('2022')>0:
2055 options.convpath = "/eos/experiment/sndlhc/convertedData/physics/2022/"
2056 options.rmin = 4361-1
2057 offline =www+"offline.html"
2058 elif options.rawData.find('2023')>0:
2059 options.convpath = "/eos/experiment/sndlhc/convertedData/physics/2023/"
2060 options.rmin = 5413-1
2061 offline =www+"offline.html"
2062 elif options.rawData.find('2024')>0:
2063 # extract the em target run from the rawData path
2064 em_run = options.rawData[options.rawData.find("run_"):]
2065 options.convpath = "/eos/experiment/sndlhc/convertedData/physics/2024/"+em_run
2066 options.rmin = 7649-1
2067 offline =www+"offline.html"
2068 elif options.rawData.find('2025')>0:
2069 # extract the em target run from the rawData path
2070 em_run = options.rawData[options.rawData.find("run_"):]
2071 options.convpath = "/eos/experiment/sndlhc/convertedData/physics/2025/"+em_run
2072 options.rmin = 10919-1 #10587 fillN
2073 offline =www+"offline.html"
2075 FS.Init(options)
2076 ut.bookCanvas(FS.h,'c1','c1',1800,900,1,1)
2077 FS.h['c1'].cd()
2078
2079 if options.command == "extract":
2080 FS.Extract()
2081 elif options.command == "draw":
2082 options.withIP2=True
2083 FS.plotBunchStructure(options.fillNumbers,int(options.runNumbers))
2084 elif options.command == "makeAll" or options.command == "update":
2085 problems = {}
2086 with client.File() as f:
2087 f.open(offline)
2088 status, L = f.read()
2089 Lhtml = L.decode().split('\n')
2090 f.close()
2091 runs = []
2092 for x in Lhtml:
2093 if x.find('#event')<0: continue
2094 if x.find('noise run')>0: continue
2095 ir = x.find("run ")
2096 if ir<0: continue
2097 k = x[ir:].find(" ")
2098 runnr = int(x[ir+4:ir+4+k+1])
2099 j = x.find("#events")
2100 tmp=x[j:].split('=')[1].split(' ')
2101 tmp.sort()
2102 for y in tmp:
2103 if y=='': continue
2104 nev = int( y )
2105 break
2106 if nev < options.nMin: continue
2107 scale = 1
2108 k = x.find('post scaled:')
2109 if k>0: scale = int(x[k+12:])
2110 if runnr == 4425: continue # no beam, but with fill number
2111 if runnr > options.rmin: runs.append([runnr,scale])
2112 for z in runs:
2113 r = z[0]
2114 options.postScale = z[1]
2115 FS.options.fillNumbers = ""
2116 FS.options.runNumbers = r
2117
2118 if options.command == "update" and r in FS.runInfo:
2119 if FS.runInfo[r]['lumiAtIP1withSNDLHC']>0:
2120 continue
2121 FS.Extract()
2122 FS.drawLumi(r)
2123# fill dictionary with useful info
2124 Nevts,partitions = FS.getEntriesPerRun(r)
2125 fillScheme = "unknown"
2126 if hasattr(FS,"lumiAtIP1"):
2127 if 'fillingScheme' in FS.lumiAtIP1:
2128 fillScheme = FS.lumiAtIP1['fillingScheme']
2129# cross check
2130 postScale = Nevts/FS.h['time'].GetEntries()
2131 problems[r]={}
2132 if abs(options.postScale/postScale -1) > 0.1 :
2133 print('!!! inconsistent postscale:',r,'html:',options.postScale,' true:',postScale)
2134 x = 1
2135 if Nevts>5E6: x = 10
2136 if Nevts>5E7: x = 100
2137 print('should be ',Nevts,x)
2138 problems[r]['postScale'] = postScale
2139 if not r in FS.LumiInt:
2140 FS.LumiInt[r]= [0,0]
2141 if not r in FS.FSdict:
2142 FS.FSdict[r]={'phaseShift1':0,'phaseShift2':0}
2143 N_ScifiTracks = FS.h['timeWt'].GetEntries()
2144 N_DSTracks = FS.h['timeWtDS'].GetEntries()
2145 FS.runInfo[r] = {'Fillnumber':int(FS.options.fillNumbers),'phaseShift1':FS.FSdict[r]['phaseShift1'],'phaseShift2':FS.FSdict[r]['phaseShift2'],
2146 'StartTime':FS.startTime,'StartTimeC':time.ctime(FS.startTime),'Entries':Nevts,'partitions':partitions,
2147 'N_scifiTracks':N_ScifiTracks*postScale,'N_DSTracks':N_DSTracks*postScale,
2148 'lumiAtIP1':FS.LumiInt[r][0],'lumiAtIP1withSNDLHC':FS.LumiInt[r][1],
2149 'OfflineMonitoring':"https://snd-lhc-monitoring.web.cern.ch/offline/run.html?file=run"+str(r).zfill(6)+".root&lastcycle",
2150 'FillingScheme':fillScheme}
2151 FS.merge()
2152 FS.storeDict(FS.FSdict,'FSdict','FSdict')
2153
2154 FS.mergeLumi()
2155 FS.storeDict(FS.LumiInt,'LumiInt','Lumidict')
2156
2157 # FS.calcMu()
2158 FS.storeDict(FS.runInfo,'runInfo','RunInfodict')
2159
2160 L,N= 0,0
2161 for r in FS.runInfo:
2162 L+=FS.runInfo[r]['lumiAtIP1withSNDLHC']
2163 N+=FS.runInfo[r]['Entries']
2164 print('total nr of events',N,' integrated luminosity ',L/1E9,'fb')
2165 print('run tracks per Lumi')
2166 FS.tracksPerLumi()
2167 FS.plotLumiPerTime()
2168 print('make Latex', 'requires up to date files on EOS! pdflatex LumiSummary.tex')
2169 FS.makeLatex()
2170 # os.system('pdflatex LumiSummary.tex')
2171
2172 print('problems',problems)
2173 # other functionalities of this manager that will be skipped for now
2174 '''
2175 print('do not forget to copy to EOS:')
2176 print('xrdcp -f Lumi.root $EOSSHIP/eos/experiment/sndlhc/www/offline/Lumi2023.root')
2177 print('xrdcp -f Lumidict.root $EOSSHIP//eos/experiment/sndlhc/convertedData/physics/2023/')
2178 print('xrdcp -f Lumidict.pkl $EOSSHIP//eos/experiment/sndlhc/convertedData/physics/2023/')
2179 print('xrdcp -f FSdict.root $EOSSHIP//eos/experiment/sndlhc/convertedData/physics/2023/')
2180 print('xrdcp -f FSdict.pkl $EOSSHIP//eos/experiment/sndlhc/convertedData/physics/2023/')
2181 print('xrdcp -f RunInfodict.root $EOSSHIP//eos/experiment/sndlhc/convertedData/physics/2023/')
2182 print('xrdcp -f RunInfodict.pkl $EOSSHIP//eos/experiment/sndlhc/convertedData/physics/2023/')
2183 print('xrdcp -f Lumi-tracks.root $EOSSHIP//eos/experiment/sndlhc/www/offline/Lumi-tracks2023.root')
2184 print('xrdcp -f LumiSummary.pdf $EOSSHIP//eos/experiment/sndlhc/www/offline/RunSummary2023.pdf')'''
getFillNrFromRunNr(self, runNumber)
storeDict(self, dictPtr, dictName, outFileName)
getNameOfFillingscheme(self, fillnr)
myPrint(self, tname, oname)
getIntegratedLumiFromPlot(dateA, dateB)
addBunchCurrent(self, fillNr, b=2)
test(self, runnr, I=True)
getLumiAtIP1(self, fillnr=None, fromnxcals=False, fromAtlas=False)
plotBunchStructure(self, fillNr, runNumber)
LumiIntegral(self, rmin, rmax)
hitMapsNormalized(self, runNumber, Q12MC=False)
tracksPerLumi(self, aRun=False)
extractPhaseShift(self, fillNr, runNumber)