SND@LHC Software
Loading...
Searching...
No Matches
SndlhcDigi.py
Go to the documentation of this file.
1import ROOT
2import os
3import shipunit as u
4import SciFiMapping
5from array import array
6
7stop = ROOT.TVector3()
8start = ROOT.TVector3()
9
10class SndlhcDigi:
11 " convert MC hits to digitized hits"
12 def __init__(self,fout, makeClusterScifi):
13
14 self.iEvent = 0
15
16 outdir=os.getcwd()
17 outfile=outdir+"/"+fout
18 self.fn = ROOT.TFile(fout,'update')
19 self.sTree = self.fn.Get("cbmsim")
20
21 # event header
22 self.header = ROOT.SNDLHCEventHeader()
23 self.eventHeader = self.sTree.Branch("EventHeader.",self.header,32000,1)
24 #scifi
25 self.digiScifi = ROOT.TClonesArray("sndScifiHit")
26 self.digiScifiBranch = self.sTree.Branch("Digi_ScifiHits",self.digiScifi,32000,1)
27 self.digiScifi2MCPoints = ROOT.TClonesArray("Hit2MCPoints")
28 self.digiScifi2MCPoints.BypassStreamer(ROOT.kTRUE)
29 self.digiScifi2MCPointsBranch = self.sTree.Branch("Digi_ScifiHits2MCPoints",self.digiScifi2MCPoints,32000,1)
30 lsOfGlobals = ROOT.gROOT.GetListOfGlobals()
31 self.scifiDet = lsOfGlobals.FindObject('Scifi')
32 self.mufiDet = lsOfGlobals.FindObject('MuFilter')
33 if not self.scifiDet or not self.mufiDet:
34 exit('detector(s) not found, stop')
35# make sipm to fibre mapping
36 self.fibresSiPM = SciFiMapping.getFibre2SiPMCPP(self.scifiDet)
37 self.siPMFibres = SciFiMapping.getSiPM2FibreCPP(self.scifiDet)
38
39 #scifi clusters
40 self.makeClusterScifi = makeClusterScifi
41 if self.makeClusterScifi:
42 self.clusScifi = ROOT.TClonesArray("sndCluster")
43 self.clusScifiBranch = self.sTree.Branch("Cluster_Scifi",self.clusScifi,32000,1)
44
45 #muonFilter
46 self.digiMuFilter = ROOT.TClonesArray("MuFilterHit")
47 self.digiMuFilterBranch = self.sTree.Branch("Digi_MuFilterHits",self.digiMuFilter,32000,1)
48 self.digiMuFilter2MCPoints = ROOT.TClonesArray("Hit2MCPoints")
49 self.digiMuFilter2MCPoints.BypassStreamer(ROOT.kTRUE)
50 self.digiMuFilter2MCPointsBranch = self.sTree.Branch("Digi_MuFilterHits2MCPoints",self.digiMuFilter2MCPoints,32000,1)
51
52 def digitize(self):
53
54 self.header.SetRunId( self.sTree.MCEventHeader.GetRunID() )
55 self.header.SetEventNumber( self.sTree.MCEventHeader.GetEventID() ) # counts from 1
56 self.header.SetBunchType(101);
57 self.eventHeader.Fill()
58 self.digiScifi.Clear('C')
59 self.digiScifi2MCPoints.Clear('C')
60 self.digitizeScifi()
61 self.digiScifiBranch.Fill()
62 self.digiScifi2MCPointsBranch.Fill()
63 if self.makeClusterScifi:
64 self.clusScifi.Clear('C')
65 self.clusterScifi()
66 self.clusScifiBranch.Fill()
67
68 self.digiMuFilter.Clear('C')
69 self.digiMuFilter2MCPoints.Clear('C')
70 self.digitizeMuFilter()
71 self.digiMuFilterBranch.Fill()
72 self.digiMuFilter2MCPointsBranch.Fill()
73
74 def digitizeScifi(self):
75 hitContainer = {}
76 mcLinks = ROOT.Hit2MCPoints()
77 mcPoints = {}
78 norm = {}
79 k=-1
80 for p in self.sTree.ScifiPoint:
81 k+=1
82 # collect all hits in same SiPM channel
83 detID = p.GetDetectorID()
84 locFibreID = detID%100000
85 if not locFibreID in self.siPMFibres: continue # fibres in dead areas
86 for sipmChan in self.siPMFibres[locFibreID]:
87 globsipmChan = int(detID/100000)*100000+sipmChan
88 if not globsipmChan in hitContainer:
89 hitContainer[globsipmChan]=[]
90 mcPoints[globsipmChan] = {}
91 norm[globsipmChan] = 0
92 w = self.siPMFibres[locFibreID][sipmChan]['weight']
93 hitContainer[globsipmChan].append([p,w])
94 dE = p.GetEnergyLoss()*w
95 mcPoints[globsipmChan][k]=dE
96 norm[globsipmChan]+= dE
97 index = 0
98 for detID in hitContainer:
99 allPoints = ROOT.std.vector('ScifiPoint*')()
100 allWeights = ROOT.std.vector('Float_t')()
101 for p in hitContainer[detID]:
102 allPoints.push_back(p[0])
103 allWeights.push_back(p[1])
104 aHit = ROOT.sndScifiHit(detID,allPoints,allWeights)
105 if self.digiScifi.GetSize() == index: self.digiScifi.Expand(index+100)
106 aHit_TCA = self.digiScifi.ConstructedAt(index)
107 ROOT.std.swap(aHit, aHit_TCA)
108 index+=1
109 for k in mcPoints[detID]:
110 mcLinks.Add(detID,k, mcPoints[detID][k]/norm[detID])
111 mcLinks_TCA = self.digiScifi2MCPoints.ConstructedAt(0)
112 ROOT.std.swap(mcLinks, mcLinks_TCA)
113
114 def digitizeMuFilter(self):
115 hitContainer = {}
116 mcLinks = ROOT.Hit2MCPoints()
117 mcPoints = {}
118 norm = {}
119 k=-1
120 for p in self.sTree.MuFilterPoint:
121 k+=1
122 # collect all hits in same detector element
123 detID = p.GetDetectorID()
124 if not detID in hitContainer:
125 hitContainer[detID]=[]
126 mcPoints[detID] = {}
127 norm[detID] = 0
128 hitContainer[detID].append(p)
129 mcPoints[detID][k]=p.GetEnergyLoss()
130 norm[detID]+= p.GetEnergyLoss()
131 index = 0
132 for detID in hitContainer:
133 allPoints = ROOT.std.vector('MuFilterPoint*')()
134 for p in hitContainer[detID]:
135 allPoints.push_back(p)
136 aHit = ROOT.MuFilterHit(detID,allPoints)
137
138 if self.digiMuFilter.GetSize() == index: self.digiMuFilter.Expand(index+100)
139 aHit_TCA = self.digiMuFilter.ConstructedAt(index)
140 ROOT.std.swap(aHit, aHit_TCA)
141 index+=1
142 for k in mcPoints[detID]:
143 mcLinks.Add(detID,k, mcPoints[detID][k]/norm[detID])
144 mcLinks_TCA = self.digiMuFilter2MCPoints.ConstructedAt(0)
145 ROOT.std.swap(mcLinks, mcLinks_TCA)
146
147 def clusterScifi(self):
148 index = 0
149 hitDict = {}
150 for k in range(self.sTree.Digi_ScifiHits.GetEntries()):
151 d = self.sTree.Digi_ScifiHits[k]
152 if not d.isValid(): continue
153 hitDict[d.GetDetectorID()] = k
154 hitList = list(hitDict.keys())
155 if len(hitList)>0:
156 hitList.sort()
157 tmp = [ hitList[0] ]
158 cprev = hitList[0]
159 ncl = 0
160 last = len(hitList)-1
161 hitlist = ROOT.std.vector("sndScifiHit*")()
162 for i in range(len(hitList)):
163 if i==0 and len(hitList)>1: continue
164 c=hitList[i]
165 if (c-cprev)==1:
166 tmp.append(c)
167 if (c-cprev)!=1 or c==hitList[last]:
168 first = tmp[0]
169 N = len(tmp)
170 hitlist.clear()
171 for aHit in tmp: hitlist.push_back( self.sTree.Digi_ScifiHits[hitDict[aHit]],)
172 aCluster = ROOT.sndCluster(first,N,hitlist,self.scifiDet)
173 if self.clusScifi.GetSize() == index: self.clusScifi.Expand(index+10)
174 aCluster_TCA = self.clusScifi.ConstructedAt(index)
175 ROOT.std.swap(aCluster, aCluster_TCA)
176 index+=1
177 if c!=hitList[last]:
178 ncl+=1
179 tmp = [c]
180 cprev = c
181
182
183 def finish(self):
184 print('finished writing tree')
185 self.sTree.Write()
SndlhcDigi.SndlhcDigi.__init__
__init__(self, fout, makeClusterScifi)
Definition SndlhcDigi.py:12
SciFiMapping.getFibre2SiPMCPP
getFibre2SiPMCPP(scifi)
Definition SciFiMapping.py:31
SciFiMapping.getSiPM2FibreCPP
getSiPM2FibreCPP(scifi)
Definition SciFiMapping.py:42