SND@LHC Software
Loading...
Searching...
No Matches
Monitor.py
Go to the documentation of this file.
1#!/usr/bin/env python
2import ROOT
3import os,sys,subprocess
4import time
5import ctypes
6from array import array
7import rootUtils as ut
8import shipunit as u
9import SndlhcGeo
10from XRootD import client
11from XRootD.client.flags import DirListFlags, OpenFlags, MkDirFlags, QueryCode
12from rootpyPickler import Unpickler
13
14A,B=ROOT.TVector3(),ROOT.TVector3()
15# for fixing a root bug, will be solved in the forthcoming 6.26 release.
16
17ROOT.gInterpreter.Declare("""
18#include "MuFilterHit.h"
19#include "AbsMeasurement.h"
20#include "TrackPoint.h"
21
22void fixRoot(MuFilterHit& aHit,std::vector<int>& key,std::vector<float>& value, bool mask) {
23 std::map<int,float> m = aHit.GetAllSignals(false);
24 std::map<int, float>::iterator it = m.begin();
25 while (it != m.end())
26 {
27 key.push_back(it->first);
28 value.push_back(it->second);
29 it++;
30 }
31}
32void fixRootT(MuFilterHit& aHit,std::vector<int>& key,std::vector<float>& value, bool mask) {
33 std::map<int,float> m = aHit.GetAllTimes(false);
34 std::map<int, float>::iterator it = m.begin();
35 while (it != m.end())
36 {
37 key.push_back(it->first);
38 value.push_back(it->second);
39 it++;
40 }
41}
42void fixRoot(MuFilterHit& aHit, std::vector<TString>& key,std::vector<float>& value, bool mask) {
43 std::map<TString, float> m = aHit.SumOfSignals();
44 std::map<TString, float>::iterator it = m.begin();
45 while (it != m.end())
46 {
47 key.push_back(it->first);
48 value.push_back(it->second);
49 it++;
50 }
51}
52
53void fixRoot(std::vector<genfit::TrackPoint*>& points, std::vector<int>& d,std::vector<int>& k, bool mask) {
54 for(std::size_t i = 0; i < points.size(); ++i) {
55 genfit::AbsMeasurement* m = points[i]->getRawMeasurement();
56 d.push_back( m->getDetId() );
57 k.push_back( int(m->getHitId()/1000) );
58 }
59}
60""")
61
62Tkey = ROOT.std.vector('TString')()
63Ikey = ROOT.std.vector('int')()
64Value = ROOT.std.vector('float')()
65
66class Monitoring():
67 " set of monitor histograms "
68 def __init__(self,options,FairTasks):
69 self.options = options
70 self.EventNumber = -1
71 self.TStart = -1
72 self.TEnd = -1
73 self.MonteCarlo = False
74 self.Weight = 1
75
76 path = options.path
77 self.myclient = None
78 if path.find('eos')>0:
79 path = options.server+options.path
80 if options.online:
81 path = path.replace("raw_data","convertedData").replace("data/","")
82 self.myclient = client.FileSystem(options.server)
83# setup geometry
84 if (options.geoFile).find('../')<0: self.snd_geo = SndlhcGeo.GeoInterface(path+options.geoFile)
85 else: self.snd_geo = SndlhcGeo.GeoInterface(options.geoFile[3:])
86 self.MuFilter = self.snd_geo.modules['MuFilter']
87 self.Scifi = self.snd_geo.modules['Scifi']
88
89# MuFilter mapping of planes and bars
90 self.systemAndPlanes = {1:self.MuFilter.GetConfParI("MuFilter/NVetoPlanes"),
91 2:self.MuFilter.GetConfParI("MuFilter/NUpstreamPlanes"),
92 3:2*self.MuFilter.GetConfParI("MuFilter/NDownstreamPlanes")-1} # to arrive at 7 DS planes
93 self.systemAndBars = {1:self.MuFilter.GetConfParI("MuFilter/NVetoBars"),
94 2:self.MuFilter.GetConfParI("MuFilter/NUpstreamBars"),
95 3:self.MuFilter.GetConfParI("MuFilter/NDownstreamBars")}
96 self.systemAndChannels = {1:[self.MuFilter.GetConfParI("MuFilter/VetonSiPMs"),0],
97 2:[self.MuFilter.GetConfParI("MuFilter/UpstreamnSiPMs")-2,2],
98 3:[self.MuFilter.GetConfParI("MuFilter/DownstreamnSiPMs"),0]}
99 self.sdict = {0:'Scifi',1:'Veto',2:'US',3:'DS'}
100
101 self.freq = 160.316E6
102 self.TDC2ns = 1E9/self.freq
103
105
106 self.h = {} # histogram storage
107
108 self.runNr = str(options.runNumber).zfill(6)
109# presenter file
110 if hasattr(self, "saveTo") and options.saveTo!="":
111 name = options.saveTo+'run'+self.runNr+'_'+str(options.nStart//1000000)+'.root'
112 else:
113 name = 'run'+self.runNr+'.root'
114 if options.interactive: name = 'I-'+name
115 self.presenterFile = ROOT.TFile(name,'recreate')
116 for role in ['','shifter', 'expert']:
117 self.presenterFile.mkdir('scifi/'+role)
118 self.presenterFile.mkdir('mufilter/'+role)
119 self.presenterFile.mkdir('daq/'+role)
120 if options.interactive: self.presenterFile.mkdir('eventdisplay/'+role)
121 self.FairTasks = {}
122 for x in FairTasks: # keeps extended methods if from python class
123 self.FairTasks[x.GetName()] = x
124
125# setup input
126 if options.online:
127 import ConvRawData
128 options.chi2Max = 2000.
129 options.saturationLimit = 0.95
130 options.stop = False
131 options.withGeoFile = True
133 self.converter.Init(options)
134 self.options.online = self.converter
135 self.fsdict = False
136 self.hasBunchInfo = False
137 self.eventTree = options.online.fSink.GetOutTree()
138 self.Nkeys = 38 # need to find a way to get this number automatically
139 if self.converter.newFormat: self.Nkeys = 1
140 for t in self.FairTasks:
141 T = self.FairTasks[t]
142 self.converter.run.AddTask(T)
143 if t=='simpleTracking':
144 T.Init(self.converter)
145 self.trackTask = self.FairTasks[t]
146 self.Reco_MuonTracks = T.fittedTracks
147 self.clusMufi = T.clusMufi
148 self.clusScifi = T.clusScifi
149 else: T.Init()
150 self.run = self.converter.run
151 if self.eventTree.EventHeader.GetAccMode()==12: # ion runs
152 self.Nbunches = 1782
153 self.div = 8
154 else: # proton runs
155 self.Nbunches = 3564
156 self.div = 4
157 return
158 else:
159 if options.fname:
160 f=ROOT.TFile.Open(options.fname)
161 eventChain = f.Get('rawConv')
162 if not eventChain:
163 eventChain = f.cbmsim
164 if eventChain.GetBranch('MCTrack'): self.MonteCarlo = True
165 partitions = []
166 else:
167 partitions = 0
168 if options.partition < 0:
169 partitions = []
170 if path.find('eos')>0:
171# check for partitions
172 dirlist = str( subprocess.check_output("xrdfs "+options.server+" ls "+options.path+"run_"+self.runNr,shell=True) )
173 for x in dirlist.split('\\n'):
174 ix = x.find('sndsw_raw-')
175 if ix<0: continue
176 partitions.append(x[ix:])
177 else:
178# check for partitions
179 dirlist = os.listdir(options.path+"run_"+self.runNr)
180 for x in dirlist:
181 if not x.find('sndsw_raw-')<0:
182 partitions.append(x)
183 else:
184 partitions = ["sndsw_raw-"+ str(options.partition).zfill(4)+".root"]
185 if options.runNumber>0:
186 eventChain = ROOT.TChain('rawConv')
187 for p in partitions:
188 eventChain.Add(path+'run_'+self.runNr+'/'+p)
189
190 rc = eventChain.GetEvent(0)
191 self.TStart = eventChain.EventHeader.GetEventTime()
192 if options.nEvents <0:
193 rc = eventChain.GetEvent(eventChain.GetEntries()-1)
194 else:
195 rc = eventChain.GetEvent(options.nEvents-1)
196 self.TEnd = eventChain.EventHeader.GetEventTime()
197 rc = eventChain.GetEvent(0)
198# start FairRunAna
199 self.run = ROOT.FairRunAna()
200 ioman = ROOT.FairRootManager.Instance()
201 ioman.SetTreeName(eventChain.GetName())
202 outFile = ROOT.TMemFile('dummy','CREATE')
203 source = ROOT.FairFileSource(eventChain.GetCurrentFile())
204 for i in range(1,len(partitions)):
205 p = partitions[i]
206 source.AddFile(path+'run_'+self.runNr+'/'+p)
207 self.run.SetSource(source)
208 self.sink = ROOT.FairRootFileSink(outFile)
209 self.run.SetSink(self.sink)
210
211 for t in FairTasks:
212 self.run.AddTask(t)
213
214#avoiding some error messages
215 xrdb = ROOT.FairRuntimeDb.instance()
216 xrdb.getContainer("FairBaseParSet").setStatic()
217 xrdb.getContainer("FairGeoParSet").setStatic()
218
219 self.run.Init()
220 if len(partitions)>0: self.eventTree = ioman.GetInChain()
221 else: self.eventTree = ioman.GetInTree()
222#fitted tracks
223 if "simpleTracking" in self.FairTasks:
224 self.trackTask = self.FairTasks["simpleTracking"]
225 self.Reco_MuonTracks = self.trackTask.fittedTracks
226 self.clusMufi = self.trackTask.clusMufi
227 self.clusScifi = self.trackTask.clusScifi
228 self.trackTask.DSnPlanes = 3
229
230# initialize detector class for access to eventheader
231 rc = eventChain.GetEvent(0)
232 if not self.MonteCarlo:
233 self.snd_geo.modules['Scifi'].InitEvent(eventChain.EventHeader)
234 self.snd_geo.modules['MuFilter'].InitEvent(eventChain.EventHeader)
235
236 # define the number of bunches in the LHC
237 if eventChain.EventHeader.GetAccMode()==12: # ion runs
238 self.Nbunches = 1782
239 self.div = 8
240 else: # proton runs
241 self.Nbunches = 3564
242 self.div = 4
243
244 # get filling scheme, only necessary if not encoded in EventHeader, before 2022 reprocessing
245 self.hasBunchInfo = False
246 self.fsdict = False
247 if hasattr(eventChain.EventHeader,"GetBunchType"):
248 if not eventChain.EventHeader.GetBunchType()<0:
249 self.hasBunchInfo = True
250 print('take bunch info from event header')
251 if not self.hasBunchInfo:
252 try:
253 fg = ROOT.TFile.Open(options.server+options.path+'FSdict.root')
254 pkl = Unpickler(fg)
255 FSdict = pkl.load('FSdict')
256 fg.Close()
257 if options.runNumber in FSdict: self.fsdict = FSdict[options.runNumber]
258 except:
259 print('continue without knowing filling scheme',options.server+options.path)
260 if self.fsdict:
261 print('extract bunch info from filling scheme')
262 if self.fsdict or self.hasBunchInfo:
263 for x in ['B1only','B2noB1','noBeam']:
264 for role in ['shifter', 'expert']:
265 self.presenterFile.mkdir('mufilter/'+role+'/'+x)
266 self.presenterFile.mkdir('scifi/'+role+'/'+x)
267
268 def GetEntries(self):
269 if self.options.online:
270 if self.converter.newFormat: return self.converter.fiN.Get('data').GetEntries()
271 else: return self.converter.fiN.event.GetEntries()
272 else:
273 return self.eventTree.GetEntries()
274
275 def updateSource(self,fname):
276 # only needed in auto mode
277 notOK = True
278 nIter = 0
279 while notOK:
280 # self.converter.fiN.Close()
281 if nIter > 100:
282 print('too many attempts to read the file ',fname,' I am exiting.')
283 quit()
284 nIter+=1
285 if self.converter.fiN.GetName() != fname:
286 self.converter.fiN = ROOT.TFile.Open(fname)
287 else:
288 self.converter.fiN = ROOT.TFile.Open(fname)
289 Nkeys = self.converter.fiN.ReadKeys()
290 if Nkeys != self.Nkeys:
291 time.sleep(10)
292 self.converter.fiN = ROOT.TFile.Open(fname)
293 print('wrong number of keys',Nkeys)
294 continue
295
296 listOfKeys = []
297 for b in self.converter.fiN.GetListOfKeys(): listOfKeys.append(b.GetName())
298 for name in listOfKeys:
299 exec("self.converter.fiN."+name+".Refresh()")
300
301 first = -1
302 for name in listOfKeys:
303 nentries = eval("self.converter.fiN."+name+".GetEntries()")
304 if first<0: first = nentries
305 if first != nentries:
306 print('wrong number of entries',first,nentries)
307 first = -1
308 break
309 if first<0: continue
310
311 for name in listOfKeys:
312 if name.find('board')==0:
313 self.converter.boards[name]=eval("self.converter.fiN."+name)
314
315 notOK = False
316
317 def GetEvent(self,n):
318 if not self.options.online: # offline, FairRoot in charge
319
320 if self.eventTree.GetBranchStatus('Reco_MuonTracks'):
321 for aTrack in self.eventTree.Reco_MuonTracks:
322 if aTrack: aTrack.Delete()
323 self.eventTree.Reco_MuonTracks.Delete()
324 if "simpleTracking" in self.FairTasks:
325 self.Reco_MuonTracks.Delete()
326
327 if self.options.online:
328 online = self.options.online
329 online.executeEvent(n)
330 self.Reco_MuonTracks.Delete()
331 if self.options.FairTask_convRaw:
332 self.options.online.sTree.GetEvent(self.options.online.sTree.GetEntries()-1)
333 for t in self.FairTasks: self.FairTasks[t].ExecuteTask()
334 self.eventTree = self.options.online.sTree
335 else:
336 self.eventTree.GetEvent(n)
337 if not self.MonteCarlo:
338 # initialize detector class for access to eventheader
339 self.snd_geo.modules['Scifi'].InitEvent(self.eventTree.EventHeader)
340 self.snd_geo.modules['MuFilter'].InitEvent(self.eventTree.EventHeader)
341 if self.MonteCarlo: self.Weight = self.eventTree.MCTrack[0].GetWeight()
342 for t in self.FairTasks:
343 if t =='simpleTracking':
344 self.FairTasks[t].ExecuteTask(self.options.trackType)
345 else:
346 self.FairTasks[t].ExecuteTask()
347 self.EventNumber = n
348
349# check for bunch xing type
350 self.xing = {'all':True,'B1only':False,'B2noB1':False,'noBeam':False}
351 if self.hasBunchInfo:
352 binfo = self.eventTree.EventHeader
353 self.xing['IP1'] = binfo.isIP1()
354 self.xing['IP2'] = binfo.isIP2()
355 self.xing['B1'] = binfo.isB1()
356 self.xing['B2'] = binfo.isB2()
357 self.xing['B1only'] = binfo.isB1Only()
358 self.xing['B2noB1'] = binfo.isB2noB1()
359 self.xing['noBeam'] = binfo.isNoBeam()
360 elif self.fsdict:
361 T = self.eventTree.EventHeader.GetEventTime()
362 bunchNumber = int((T%(self.div*self.Nbunches))/self.div+0.5)
363 nb1 = (self.Nbunches + bunchNumber - self.fsdict['phaseShift1'])%self.Nbunches
364 nb2 = (self.Nbunches + bunchNumber - self.fsdict['phaseShift1']- self.fsdict['phaseShift2'])%self.Nbunches
365 b1 = nb1 in self.fsdict['B1']
366 b2 = nb2 in self.fsdict['B2']
367 IP1 = False
368 IP2 = False
369 if b1:
370 IP1 = self.fsdict['B1'][nb1]['IP1']
371 if b2:
372 IP2 = self.fsdict['B2'][nb2]['IP2']
373 self.xing['IP1'] = IP1
374 self.xing['IP2'] = IP2
375 self.xing['B1'] = b1
376 self.xing['B2'] = b2
377 self.xing['B1only'] = b1 and not IP1 and not b2
378 self.xing['B2noB1'] = b2 and not b1
379 self.xing['noBeam'] = not b1 and not b2
380 if self.xing['B1only'] and self.xing['B2noB1'] or self.xing['B1only'] and self.xing['noBeam'] : print('error with b1only assignment',self.xing)
381 if self.xing['B2noB1'] and self.xing['noBeam'] : print('error with b2nob1 assignment',self.xing)
382
383 return self.eventTree
384
386 # try to copy root file with TCanvas to EOS
387 self.presenterFile.Close()
388 if self.options.online:
389 wwwPath = "/eos/experiment/sndlhc/www/online"
390 elif self.options.path.find('2022')>0 :
391 wwwPath = "/eos/experiment/sndlhc/www/reprocessing"
392 else:
393 wwwPath = "/eos/experiment/sndlhc/www/offline"
394 if self.options.sudo:
395 try:
396 rc = os.system("xrdcp -f "+self.presenterFile.GetName()+" "+os.environ['EOSSHIP']+wwwPath)
397 except:
398 print("copy of root file failed. Token expired?")
399 self.presenterFile = ROOT.TFile('run'+self.runNr+'.root','update')
400
402 wwwPath = "/eos/experiment/sndlhc/www/online/"
403 for r in self.options.runNumbers.split(','):
404 if r!= '': runList.append(int(r))
405 for r in runList:
406 runNr = str(r).zfill(6)+'.root'
407 f = ROOT.TFile.Open(options.server+wwwPath+ runNr)
408 g = ROOT.TFile('tmp.root','recreate')
409 for d in f.GetListOfKeys():
410 D = f.Get(d.GetName())
411 D.Purge()
412 for d in f.GetListOfKeys():
413 name = d.GetName()
414 s = f.Get(name)
415 g.mkdir(name)
416 g.cd(name)
417 for x in s.GetListOfKeys():
418 s.Get(x.GetName()).Write()
419 g.Close()
420 f.Close()
421 os.system('xrdcp -f tmp.root '+self.options.server+options.path+rName)
422
423 def updateHtml(self):
424 if self.options.online: destination="online"
425 elif self.options.path.find('2022')>0: destination="reprocessing"
426 else: destination="offline"
427 rc = os.system("xrdcp -f "+os.environ['EOSSHIP']+"/eos/experiment/sndlhc/www/"+destination+".html . ")
428 old = open(destination+".html")
429 oldL = old.readlines()
430 old.close()
431 tmp = open("tmp.html",'w')
432 found = False
433 for L in oldL:
434 if not L.find(self.runNr)<0: return
435 if L.find("https://snd-lhc-monitoring.web.cern.ch/"+destination+"/run.html?file=run")>0 and not found:
436 r = str(self.options.runNumber)
437 Lnew = ' <li> <a href="https://snd-lhc-monitoring.web.cern.ch/'+destination+'/run.html?file=run'
438 Lnew+= self.runNr+'.root&lastcycle">run '+r+' </a>'+self.options.startTime
439 if self.options.postScale>1: Lnew+=" post scaled:"+str(self.options.postScale)
440 Lnew+='\n'
441 tmp.write(Lnew)
442 found = True
443 tmp.write(L)
444 tmp.close()
445 os.system('cp '+destination+'.html '+destination+time.ctime().replace(' ','')+'.html ') # make backup
446 os.system('cp tmp.html '+destination+'.html')
447 try:
448 rc = os.system("xrdcp -f "+destination+".html "+os.environ['EOSSHIP']+"/eos/experiment/sndlhc/www/")
449 except:
450 print("copy of html failed. Token expired?")
451
452 def cleanUpHtml(self,destination="online"):
453 rc = os.system("xrdcp -f "+os.environ['EOSSHIP']+"/eos/experiment/sndlhc/www/"+destination+".html . ")
454 old = open(destination+".html")
455 oldL = old.readlines()
456 old.close()
457 tmp = open("tmp.html",'w')
458 dirlist = str( subprocess.check_output("xrdfs "+os.environ['EOSSHIP']+" ls /eos/experiment/sndlhc/www/"+destination+"/",shell=True) )
459 for L in oldL:
460 OK = True
461 if L.find("https://snd-lhc-monitoring.web.cern.ch/"+destination)>0:
462 k = L.find("file=")+5
463 m = L.find(".root")+5
464 R = L[k:m]
465 if not R in dirlist: OK = False
466 if OK: tmp.write(L)
467 tmp.close()
468 os.system('cp tmp.html '+destination+'.html')
469 try:
470 rc = os.system("xrdcp -f "+destination+".html "+os.environ['EOSSHIP']+"/eos/experiment/sndlhc/www/")
471 except:
472 print("copy of html failed. Token expired?")
473
474 def checkAlarm(self,minEntries=1000):
475 self.alarms = {'scifi':[]}
476 # check for empty histograms in scifi signal
477 entries = {}
478 for mat in range(30):
479 entries[mat] = self.h['scifi-mat_'+str(mat)].GetEntries()
480 res = sorted(entries.items(), key=operator.itemgetter(1),reverse=True)
481 if res[0][1]>minEntries: # choose limit which makes it sensitive to check for empty mats
482 for mat in range(30):
483 if entries[mat] < 1:
484 s = mat//6 + 1
485 if mat%6 < 3 : p='H'
486 else : p='V'
487 e = str(s)+p+str(mat%3)
488 self.alarms['scifi'].append(e)
489 if len(self.alarms['scifi']) > 0: self.sendAlarm()
490
491 self.alarms['Veto']=[]
492 # check for empty histograms in veto signal
493 entries = {}
494 s = 1
495 for p in range(2):
496 for side in ['L','R']:
497 entries[str(10*s+p)+side] = self.h['mufi-sig'+p+str(10*s+p)].GetEntries()
498 res = sorted(entries.items(), key=operator.itemgetter(1),reverse=True)
499 if res[0][1]>minEntries: # choose limit which makes it sensitive to check for empty boards
500 for p in range(2):
501 for side in ['L','R']:
502 if entries[str(10*s+p)+side] < 1:
503 self.alarms['Veto'].append(str(10*s+p)+side())
504 # check for empty histograms in US signal
505 entries = {}
506 s = 2
507 for p in range(5):
508 for side in ['L','R']:
509 entries[str(s*10+p)+side] = self.h['mufi-sig2'+p+str(s*10+p)].GetEntries()
510 res = sorted(entries.items(), key=operator.itemgetter(1),reverse=True)
511 if res[0][1]>1000: # choose limit which makes it sensitive to check for empty boards
512 for p in range(5):
513 for side in ['L','R']:
514 if entries[str(s*10+p)+side] < 1:
515 self.alarms['Veto'].append(str(s*10+p)+side())
516
517
518 def sendAlarm(self):
519 print('ALARM',self.alarms)
520
521 def systemAndOrientation(self,s,plane):
522 if s==1 or s==2: return "horizontal"
523 if plane%2==1 or plane == 6: return "vertical"
524 return "horizontal"
525
527 zPos={'MuFilter':{},'Scifi':{}}
528 for s in self.systemAndPlanes:
529 for plane in range(self.systemAndPlanes[s]):
530 bar = 4
531 p = plane
532 if s==3 and (plane%2==0 or plane==7):
533 bar = 90
534 p = plane//2
535 elif s==3 and plane%2==1:
536 bar = 30
537 p = plane//2
538 self.MuFilter.GetPosition(s*10000+p*1000+bar,A,B)
539 zPos['MuFilter'][s*10+plane] = (A.Z()+B.Z())/2.
540 for s in range(1,self.Scifi.GetConfParI("Scifi/nscifi")+1):
541 mat = 1
542 sipm = 1
543 channel = 64
544 for o in range(2):
545 self.Scifi.GetPosition(channel+1000*sipm+10000*mat+100000*o+1000000*s,A,B)
546 zPos['Scifi'][s*10+o] = (A.Z()+B.Z())/2.
547 return zPos
548
549 def smallSiPMchannel(self,i):
550 if i==2 or i==5 or i==10 or i==13: return True
551 else: return False
552
553# Scifi specific code
554 def Scifi_xPos(self,detID):
555 orientation = (detID//100000)%10
556 nStation = 2*(detID//1000000-1)+orientation
557 mat = (detID%100000)//10000
558 X = detID%1000+(detID%10000)//1000*128
559 return [nStation,mat,X] # even numbers are Y (horizontal plane), odd numbers X (vertical plane)
560
561# decode MuFilter detID
562 def MuFilter_PlaneBars(self,detID):
563 s = detID//10000
564 l = (detID%10000)//1000 # plane number
565 bar = (detID%1000)
566 if s>2:
567 l=2*l
568 if bar>59:
569 bar=bar-60
570 if l<6: l+=1
571 return {'station':s,'plane':l,'bar':bar}
572
573 def map2Dict(self,aHit,T='GetAllSignals',mask=True):
574 if T=="SumOfSignals":
575 key = Tkey
576 elif T=="GetAllSignals" or T=="GetAllTimes":
577 key = Ikey
578 else:
579 print('use case not known',T)
580 1/0
581 key.clear()
582 Value.clear()
583 if T=="GetAllTimes": ROOT.fixRootT(aHit,key,Value,mask)
584 else: ROOT.fixRoot(aHit,key,Value,mask)
585 theDict = {}
586 for k in range(key.size()):
587 if T=="SumOfSignals": theDict[key[k].Data()] = Value[k]
588 else: theDict[key[k]] = Value[k]
589 return theDict
590
591 def fit_langau(self,hist,o,bmin,bmax,opt=''):
592 if opt == 2:
593 params = {0:'Width(scale)',1:'mostProbable',2:'norm',3:'sigma',4:'N2'}
594 F = ROOT.TF1('langau',langaufun,0,200,len(params))
595 else:
596 params = {0:'Width(scale)',1:'mostProbable',2:'norm',3:'sigma'}
597 F = ROOT.TF1('langau',twoLangaufun,0,200,len(params))
598 for p in params: F.SetParName(p,params[p])
599 rc = hist.Fit('landau','S'+o,'',bmin,bmax)
600 res = rc.Get()
601 if not res: return res
602 F.SetParameter(2,res.Parameter(0))
603 F.SetParameter(1,res.Parameter(1))
604 F.SetParameter(0,res.Parameter(2))
605 F.SetParameter(3,res.Parameter(2))
606 F.SetParameter(4,0)
607 F.SetParLimits(0,0,100)
608 F.SetParLimits(1,0,100)
609 F.SetParLimits(3,0,10)
610 rc = hist.Fit(F,'S'+o,'',bmin,bmax)
611 res = rc.Get()
612 return res
613
614 def twoLangaufun(self,x,par):
615 N1 = self.langaufun(x,par)
616 par2 = [par[0],par[1]*2,par[4],par[3]]
617 N2 = self.langaufun(x,par2)
618 return N1+abs(N2)
619
620 def langaufun(self,x,par):
621 #Fit parameters:
622 #par[0]=Width (scale) parameter of Landau density
623 #par[1]=Most Probable (MP, location) parameter of Landau density
624 #par[2]=Total area (integral -inf to inf, normalization constant)
625 #par[3]=Width (sigma) of convoluted Gaussian function
626 #
627 #In the Landau distribution (represented by the CERNLIB approximation),
628 #the maximum is located at x=-0.22278298 with the location parameter=0.
629 #This shift is corrected within this function, so that the actual
630 #maximum is identical to the MP parameter.
631#
632 # Numeric constants
633 invsq2pi = 0.3989422804014 # (2 pi)^(-1/2)
634 mpshift = -0.22278298 # Landau maximum location
635#
636 # Control constants
637 np = 100.0 # number of convolution steps
638 sc = 5.0 # convolution extends to +-sc Gaussian sigmas
639#
640 # Variables
641 summe = 0.0
642#
643 # MP shift correction
644 mpc = par[1] - mpshift * par[0]
645#
646 # Range of convolution integral
647 xlow = max(0,x[0] - sc * par[3])
648 xupp = x[0] + sc * par[3]
649#
650 step = (xupp-xlow) / np
651#
652 # Convolution integral of Landau and Gaussian by sum
653 i=1.0
654 if par[0]==0 or par[3]==0: return 9999
655 while i<=np/2:
656 i+=1
657 xx = xlow + (i-.5) * step
658 fland = ROOT.TMath.Landau(xx,mpc,par[0]) / par[0]
659 summe += fland * ROOT.TMath.Gaus(x[0],xx,par[3])
660#
661 xx = xupp - (i-.5) * step
662 fland = ROOT.TMath.Landau(xx,mpc,par[0]) / par[0]
663 summe += fland * ROOT.TMath.Gaus(x[0],xx,par[3])
664#
665 return (par[2] * step * summe * invsq2pi / par[3])
666
667 def myPrint(self,tc,name,subdir='',withRootFile=True):
668 srun = 'run'+str(self.options.runNumber)
669 if isinstance(tc, ROOT.TCanvas) :
670 tc.Update()
671 if withRootFile:
672 self.presenterFile.cd(subdir)
673 tc.Write()
674 else:
675 if not os.path.isdir(srun): os.system('mkdir '+srun)
676 pname = srun+'/'+name+'-'+srun
677 tc.Print(pname+'.png')
678 tc.Print(pname+'.pdf')
679
680 def fillHist1(self,hname,parx):
681 for x in ['','B1only','B2noB1','noBeam']:
682 if x=='':
683 rc = self.h[hname].Fill(parx,self.Weight)
684 elif self.xing[x]:
685 rc = self.h[hname+x].Fill(parx,self.Weight)
686 def fillHist2(self,hname,parx,pary):
687 for x in ['','B1only','B2noB1','noBeam']:
688 if x=='':
689 rc = self.h[hname].Fill(parx,pary,self.Weight)
690 elif self.xing[x]:
691 rc = self.h[hname+x].Fill(parx,pary,self.Weight)
692
694 " run reconstruction, select events with tracks"
695 def __init__(self,options):
696 self.options = options
697 self.EventNumber = -1
698 self.MonteCarlo = False
699
700 path = options.path
701 if path.find('eos')>0:
702 path = options.server+options.path
703# setup geometry
704 if (options.geoFile).find('../')<0: self.snd_geo = SndlhcGeo.GeoInterface(path+options.geoFile)
705 else: self.snd_geo = SndlhcGeo.GeoInterface(options.geoFile[3:])
706 self.MuFilter = self.snd_geo.modules['MuFilter']
707 self.Scifi = self.snd_geo.modules['Scifi']
708
709 self.runNr = str(options.runNumber).zfill(6)
710 if options.runNumber > 0:
711 if options.partition < 0:
712 partitions = []
713 if path.find('eos')>0:
714# check for partitions
715 print("xrdfs "+options.server+" ls "+options.path+"run_"+self.runNr)
716 dirlist = str( subprocess.check_output("xrdfs "+options.server+" ls "+options.path+"run_"+self.runNr,shell=True) )
717 for x in dirlist.split('\\n'):
718 ix = x.find('sndsw_raw-')
719 if ix<0: continue
720 partitions.append(x[ix:])
721 else:
722# check for partitions
723 dirlist = os.listdir(options.path+"run_"+self.runNr)
724 for x in dirlist:
725 if not x.find("sndsw_raw-")<0:
726 partitions.append(x)
727 else:
728 partitions = ["sndsw_raw-"+ str(options.partition).zfill(4)+".root"]
729 eventChain = ROOT.TChain('rawConv')
730 for p in partitions:
731 eventChain.Add(path+'run_'+self.runNr+'/'+p)
732 else:
733# for MC data
734 eventChain = ROOT.TChain("cbmsim")
735 eventChain.Add(options.fname)
736 self.MonteCarlo = True
737 partitions = []
738 rc = eventChain.GetEvent(0)
739 self.snd_geo.modules['Scifi'].InitEvent(eventChain.EventHeader)
740 self.snd_geo.modules['MuFilter'].InitEvent(eventChain.EventHeader)
741# start FairRunAna
742 self.run = ROOT.FairRunAna()
743 ioman = ROOT.FairRootManager.Instance()
744 ioman.SetTreeName(eventChain.GetName())
745 source = ROOT.FairFileSource(eventChain.GetCurrentFile())
746 ioman.SetInChain(eventChain)
747 first = True
748 for p in partitions:
749 if first:
750 first = False
751 continue
752 source.AddFile(path+'run_'+self.runNr+'/'+p+'.root')
753 self.run.SetSource(source)
754
755#avoiding some error messages
756 xrdb = ROOT.FairRuntimeDb.instance()
757 xrdb.getContainer("FairBaseParSet").setStatic()
758 xrdb.getContainer("FairGeoParSet").setStatic()
759
760# init() tracking tasks
761 if self.options.HoughTracking:
762 if self.options.trackType == 'Scifi' or self.options.trackType == 'ScifiDS':
763 self.muon_reco_task_Sf = options.FairTasks["houghTransform_Sf"]
764 self.muon_reco_task_Sf.Init()
765 self.genfitTrack = self.muon_reco_task_Sf.genfitTrack
766 if self.options.trackType == 'DS' or self.options.trackType == 'ScifiDS':
767 self.muon_reco_task_DS = options.FairTasks["houghTransform_DS"]
768 self.muon_reco_task_DS.Init()
769 self.genfitTrack = self.muon_reco_task_DS.genfitTrack
770 if self.options.simpleTracking:
771 self.trackTask = options.FairTasks["simpleTracking"]
772 if not self.options.HoughTracking:
773 self.genfitTrack = self.options.genfitTrack
774
775 self.trackTask.SetTrackClassType(self.genfitTrack)
776 self.trackTask.Init()
777
778# prepare output tree, same branches as input plus track(s)
779 self.outFile = ROOT.TFile(options.oname,'RECREATE')
780 self.fSink = ROOT.FairRootFileSink(self.outFile)
781
782 self.outTree = eventChain.CloneTree(0)
783 ROOT.gDirectory.pwd()
784
785 # after track tasks init(), output track format is known
786 if self.genfitTrack:
787 self.fittedTracks = ROOT.TClonesArray("genfit::Track")
788 self.fittedTracks.BypassStreamer(ROOT.kFALSE)
789 else:
790 self.fittedTracks = ROOT.TClonesArray("sndRecoTrack")
791 self.MuonTracksBranch = self.outTree.Branch("Reco_MuonTracks",self.fittedTracks,32000,0)
792
793 if self.options.simpleTracking and not self.options.trackType.find('Scifi')<0 and not eventChain.GetBranch("Cluster_Scifi"):
794 self.clusScifi = ROOT.TClonesArray("sndCluster")
795 self.clusScifiBranch = self.outTree.Branch("Cluster_Scifi",self.clusScifi,32000,1)
796 if self.options.simpleTracking and not self.options.trackType.find('DS')<0 and not eventChain.GetBranch("Cluster_Mufi"):
797 self.clusMufi = ROOT.TClonesArray("sndCluster")
798 self.clusMufiBranch = self.outTree.Branch("Cluster_Mufi",self.clusMufi,32000,1)
799
800 B = ROOT.TList()
801 B.SetName('BranchList')
802 B.Add(ROOT.TObjString('Reco_MuonTracks'))
803 B.Add(ROOT.TObjString('Scifi_sndCluster'))
804 B.Add(ROOT.TObjString('Mufi_sndCluster'))
805 B.Add(ROOT.TObjString('sndScifiHit'))
806 B.Add(ROOT.TObjString('MuFilterHit'))
807 B.Add(ROOT.TObjString('FairEventHeader'))
808 self.fSink.WriteObject(B,"BranchList", ROOT.TObject.kSingleKey)
809 self.fSink.SetRunId(options.runNumber)
810 self.fSink.SetOutTree(self.outTree)
811
812 self.eventTree = eventChain
813 self.run.SetSink(self.fSink)
814 self.OT = ioman.GetSink().GetOutTree()
815
816 def ExecuteEvent(self,event):
817 track_container_list = []
818 # Delete SndlhcTracking fitted tracks container
819 if self.options.simpleTracking:
820 self.trackTask.fittedTracks.Delete()
821
822 if self.options.trackType == 'ScifiDS':
823 if self.options.HoughTracking:
824 self.muon_reco_task_Sf.Exec(0)
825 self.muon_reco_task_DS.Exec(0)
826 track_container_list = [self.muon_reco_task_Sf.kalman_tracks,self.muon_reco_task_DS.kalman_tracks]
827 if self.options.simpleTracking:
828 self.trackTask.ExecuteTask(option='ScifiDS')
829 track_container_list.append(self.trackTask.fittedTracks)
830
831 elif self.options.trackType == 'Scifi':
832 if self.options.HoughTracking:
833 self.muon_reco_task_Sf.Exec(0)
834 track_container_list.append(self.muon_reco_task_Sf.kalman_tracks)
835 if self.options.simpleTracking:
836 self.trackTask.ExecuteTask(option='Scifi')
837 track_container_list.append(self.trackTask.fittedTracks)
838
839 elif self.options.trackType == 'DS':
840 if self.options.HoughTracking:
841 self.muon_reco_task_DS.Exec(0)
842 track_container_list.append(self.muon_reco_task_DS.kalman_tracks)
843 if self.options.simpleTracking:
844 self.trackTask.ExecuteTask(option='DS')
845 track_container_list.append(self.trackTask.fittedTracks)
846
847 i_muon = -1
848 for item in track_container_list:
849 for aTrack in item:
850 i_muon += 1
851 self.fittedTracks[i_muon] = aTrack
852
853 def Execute(self):
854 for n in range(self.options.nStart,self.options.nStart+self.options.nEvents):
855
856 if self.options.scaleFactor > 1:
857 if ROOT.gRandom.Rndm() > 1.0/self.options.scaleFactor: continue
858
859 self.eventTree.GetEvent(n)
860 # delete track containers
861 self.fittedTracks.Delete()
862
863 self.ExecuteEvent(self.eventTree)
864 if not self.MonteCarlo and self.fittedTracks.GetEntries() == 0: continue
865 if self.options.simpleTracking and not self.options.trackType.find('Scifi')<0:
866 if not self.eventTree.GetBranch("Cluster_Scifi"):
867 self.clusScifi.Delete()
868 self.clusScifi.Expand(len(self.trackTask.clusScifi))
869 for index, aCl in enumerate(self.trackTask.clusScifi):
870 self.clusScifi[index] = aCl
871 if self.options.simpleTracking and not self.options.trackType.find('DS')<0:
872 if not self.eventTree.GetBranch("Cluster_Mufi"):
873 self.clusMufi.Delete()
874 self.clusMufi.Expand(len(self.trackTask.clusMufi))
875 for index, aCl in enumerate(self.trackTask.clusMufi):
876 self.clusMufi[index] = aCl
877
878 # if using FairEventHeader, i.e. before sndlhc header was introduced
879 if hasattr(self.OT.EventHeader, "SetMCEntryNumber"):
880 self.OT.EventHeader.SetMCEntryNumber(n)
881 self.fSink.Fill()
882
883 def Finalize(self):
884 self.fSink.Write()
885 self.outFile.Close()
__init__(self, options, FairTasks)
Definition Monitor.py:68
publishRootFile(self)
Definition Monitor.py:385
MuFilter_PlaneBars(self, detID)
Definition Monitor.py:562
myPrint(self, tc, name, subdir='', withRootFile=True)
Definition Monitor.py:667
map2Dict(self, aHit, T='GetAllSignals', mask=True)
Definition Monitor.py:573
Scifi_xPos(self, detID)
Definition Monitor.py:554
twoLangaufun(self, x, par)
Definition Monitor.py:614
GetEvent(self, n)
Definition Monitor.py:317
fit_langau(self, hist, o, bmin, bmax, opt='')
Definition Monitor.py:591
getAverageZpositions(self)
Definition Monitor.py:526
purgeMonitorHistos(self)
Definition Monitor.py:401
langaufun(self, x, par)
Definition Monitor.py:620
fillHist2(self, hname, parx, pary)
Definition Monitor.py:686
checkAlarm(self, minEntries=1000)
Definition Monitor.py:474
cleanUpHtml(self, destination="online")
Definition Monitor.py:452
smallSiPMchannel(self, i)
Definition Monitor.py:549
fillHist1(self, hname, parx)
Definition Monitor.py:680
systemAndOrientation(self, s, plane)
Definition Monitor.py:521
updateSource(self, fname)
Definition Monitor.py:275
__init__(self, options)
Definition Monitor.py:695
ExecuteEvent(self, event)
Definition Monitor.py:816