FairShip
Loading...
Searching...
No Matches
shipVertex.py
Go to the documentation of this file.
1# SPDX-License-Identifier: LGPL-3.0-or-later
2# SPDX-FileCopyrightText: Copyright CERN for the benefit of the SHiP Collaboration
3
4# simple vertex reconstruction with errors
5import ctypes
6import math
7from array import array
8
9import global_variables
10import numpy as np
11import ROOT
12import rootUtils as ut
13
14
15class Task:
16 "initialize"
17
18 def __init__(self, hp, sTree, mcTree=None) -> None:
19 self.sTree = sTree
20 self.mcTree = mcTree if mcTree is not None else sTree
21 self.fPartArray = ROOT.std.vector("ShipParticle")()
22 if not self.sTree.GetBranch("Particles"):
23 self.Particles = self.sTree.Branch("Particles", self.fPartArray)
24 else:
25 self.Particles = self.sTree.Particles
26 self.reps, self.states, self.newPosDir = {}, {}, {}
27 self.LV = {1: ROOT.TLorentzVector(), 2: ROOT.TLorentzVector()}
28 self.h = hp
29 self.PDG = ROOT.TDatabasePDG.Instance()
30 self.fitTrackLoc = "FitTracks"
31 self.goodTracksLoc = "goodTracks"
32 # ut.bookHist(self.h,'Vzpull','Vz pull',100,-3.,3.)
33 # ut.bookHist(self.h,'Vxpull','Vx pull',100,-3.,3.)
34 # ut.bookHist(self.h,'Vypull','Vy pull',100,-3.,3.)
35
36 # ut.bookHist(self.h,'dVx','vertex X residual;X^{RECO}-X^{MC}, cm',400,-10.,10.)
37 # ut.bookHist(self.h,'dVy','vertex Y residual;Y^{RECO}-Y^{MC}, cm',400,-10.,10.)
38 # ut.bookHist(self.h,'dVz','vertex Z residual;Z^{RECO}-Z^{MC}, cm',500,-50.,50.)
39 # ut.bookHist(self.h,'VzpullFit','Vz pull, chi2 fit',100,-3.,3.)
40 # ut.bookHist(self.h,'VxpullFit','Vx pull, chi2 fit',100,-3.,3.)
41 # ut.bookHist(self.h,'VypullFit','Vy pull, chi2 fit',100,-3.,3.)
42 # ut.bookHist(self.h,'dVxFit','vertex X residual, chi2 fit;X^{RECO}-X^{MC}, cm',400,-10.,10.)
43 # ut.bookHist(self.h,'dVyFit','vertex Y residual, chi2 fit;Y^{RECO}-Y^{MC}, cm',400,-10.,10.)
44 # ut.bookHist(self.h,'dVzFit','vertex Z residual, chi2 fit;Z^{RECO}-Z^{MC}, cm',500,-50.,50.)
45 # ut.bookHist(self.h,'dMFit','HNL mass resolution, chi2 fit;M^{RECO}-M^{MC}, GeV',2000,-1.,1.)
46 # ut.bookHist(self.h,'MpullFit','M pull, chi2 fit',100,-3.,3.)
47
48 # ut.bookHist(self.h,'N_raveVtx','Number of RAVE vtx',5,-0.5,4.5)
49 # ut.bookHist(self.h,'N_Vtx','Number of vtx',5,-0.5,4.5)
50
51 def execute(self) -> None:
52 # make particles persistent
53 self.TwoTrackVertex()
54
55 # define global data and functions for vertex fit with TMinuit
56 y_data = np.array([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
57 z0 = 0
58 Vy = np.zeros(100)
59
60 def chi2(self, res, Vy) -> int:
61 s = 0
62 for i in range(100):
63 s += Vy[i] * res[i // 10] * res[i % 10]
64 return s
65
66 def residuals(self, y_data, a, z0: int):
67 res = np.zeros(10)
68 res[0] = abs(y_data[0]) - a[5]
69 res[1] = y_data[1] - a[3]
70 res[2] = y_data[2] - a[4]
71 res[3] = y_data[3] - a[0] - a[3] * (a[2] - z0)
72 res[4] = y_data[4] - a[1] - a[4] * (a[2] - z0)
73 res[5] = abs(y_data[5]) - a[8]
74 res[6] = y_data[6] - a[6]
75 res[7] = y_data[7] - a[7]
76 res[8] = y_data[8] - a[0] - a[6] * (a[2] - z0)
77 res[9] = y_data[9] - a[1] - a[7] * (a[2] - z0)
78 return res
79
80 def fcn(self, npar, gin, f, par, iflag) -> None:
81 res = self.residuals(self.y_datay_data, par, self.z0z0)
82 self.chi2(res, self.VyVy)
83 return
84
85 def TwoTrackVertex(self) -> None:
86 self.fPartArray.clear()
87 fittedTracks = getattr(self.sTree, self.fitTrackLoc)
88 goodTracks = getattr(self.sTree, self.goodTracksLoc)
89 if len(goodTracks) < 2:
90 return
91 PosDirCharge, CovMat, scalFac = {}, {}, {}
92 for tr in goodTracks:
93 fittedTracks[tr].getFitStatus()
94 xx = fittedTracks[tr].getFittedState()
95 pid = xx.getPDG()
96 if not global_variables.pidProton and abs(pid) == 2212:
97 pid = int(math.copysign(211, pid))
98 rep = ROOT.genfit.RKTrackRep(xx.getPDG())
99 state = ROOT.genfit.StateOnPlane(rep)
100 rep.setPosMom(state, xx.getPos(), xx.getMom())
101 PosDirCharge[tr] = {
102 "position": xx.getPos(),
103 "direction": xx.getDir(),
104 "momentum": xx.getMom(),
105 "charge": xx.getCharge(),
106 "pdgCode": pid,
107 "state": xx,
108 "rep": rep,
109 "newstate": state,
110 }
111 CovMat[tr] = xx.get6DCov()
112 #
113 if len(PosDirCharge) < 2:
114 return
115 if len(PosDirCharge) > 4:
116 return # abort too busy events
117 for t1 in PosDirCharge:
118 c1 = PosDirCharge[t1]["charge"]
119 for t2 in PosDirCharge:
120 if not t2 > t1:
121 continue
122 # ignore this for background studies
123 if PosDirCharge[t2]["charge"] == c1:
124 continue
125 newPos, doca = self.VertexError(t1, t2, PosDirCharge)
126 # as we have learned, need iterative procedure
127 dz = 99999.0
128 rc = True
129 step = 0
130 while dz > 0.01:
131 zBefore = newPos[2]
132 # make a new rep for track 1,2
133 for tr in [t1, t2]:
134 try:
135 PosDirCharge[tr]["rep"].extrapolateToPoint(PosDirCharge[tr]["newstate"], newPos, False)
136 except Exception:
137 ut.reportError("shipVertex: extrapolation did not work")
138 rc = False
139 break
140 self.newPosDir[tr] = {
141 "position": PosDirCharge[tr]["rep"].getPos(PosDirCharge[tr]["newstate"]),
142 "direction": PosDirCharge[tr]["rep"].getDir(PosDirCharge[tr]["newstate"]),
143 "momentum": PosDirCharge[tr]["rep"].getMom(PosDirCharge[tr]["newstate"]),
144 }
145 if not rc:
146 break
147 newPos, doca = self.VertexError(t1, t2, self.newPosDir)
148 dz = abs(zBefore - newPos[2])
149 step += 1
150 if step > 10:
151 ut.reportError("shipVertex: abort iteration, too many steps")
152 if global_variables.debug:
153 print(
154 "abort iteration, too many steps, pos=",
155 newPos[0],
156 newPos[1],
157 newPos[2],
158 " doca=",
159 doca,
160 "z before and dz",
161 zBefore,
162 dz,
163 )
164 rc = False
165 break
166 #
167 if not rc:
168 continue # extrapolation failed, makes no sense to continue
169 # now go for the last step and vertex error
170 scalFac[t1] = (
171 (PosDirCharge[t1]["position"][2] - newPos[2])
172 / PosDirCharge[t1]["direction"][2]
173 / PosDirCharge[t1]["momentum"].Mag()
174 )
175 scalFac[t2] = (
176 (PosDirCharge[t2]["position"][2] - newPos[2])
177 / PosDirCharge[t2]["direction"][2]
178 / PosDirCharge[t2]["momentum"].Mag()
179 )
180 HNLPos, _covX, _dist = self.VertexError(t1, t2, self.newPosDir, CovMat, scalFac)
181 # monitor Vx resolution and pulls
182 # print "DEBUG",HNLPos[0],HNLPos[1],HNLPos[2],dist,covX[0][0],covX[1][1],covX[2][2]
183 # print " ",mctrack.GetStartX(),mctrack.GetStartY(),mctrack.GetStartZ()
184 # HNL true
185 if self.sTree.GetBranch("fitTrack2MC"):
186 mctrack = self.mcTree.MCTrack[self.sTree.fitTrack2MC[t1]]
187 self.mcTree.MCTrack[self.sTree.fitTrack2MC[t2]]
188 self.mcTree.MCTrack[mctrack.GetMotherId()]
189 # print "true vtx: ",mctrack.GetStartX(),mctrack.GetStartY(),mctrack.GetStartZ()
190 # print "reco vtx: ",HNLPos[0],HNLPos[1],HNLPos[2]
191 # self.h['Vzpull'].Fill( (mctrack.GetStartZ()-HNLPos[2])/ROOT.TMath.Sqrt(covX[2][2]) )
192 # self.h['Vxpull'].Fill( (mctrack.GetStartX()-HNLPos[0])/ROOT.TMath.Sqrt(covX[0][0]) )
193 # self.h['Vypull'].Fill( (mctrack.GetStartY()-HNLPos[1])/ROOT.TMath.Sqrt(covX[1][1]) )
194 # self.h['dVx'].Fill( (mctrack.GetStartX()-HNLPos[0]) )
195 # self.h['dVy'].Fill( (mctrack.GetStartY()-HNLPos[1]) )
196 # self.h['dVz'].Fill( (mctrack.GetStartZ()-HNLPos[2]) )
197
198 # print "*********************************** vertex fit precise ******************************************** "
199
200 detPlane = ROOT.genfit.DetPlane(
201 ROOT.TVector3(0, 0, HNLPos[2]), ROOT.TVector3(1, 0, 0), ROOT.TVector3(0, 1, 0)
202 )
203 detPlanePtr = ROOT.genfit.SharedPlanePtr(detPlane)
204 st1 = fittedTracks[t1].getFittedState()
205 st2 = fittedTracks[t2].getFittedState()
206 try:
207 st1.extrapolateToPlane(detPlanePtr)
208 except Exception:
209 ut.reportError("shipVertex.TwoTrackVertex: extrapolation did not work")
210 continue
211 try:
212 st2.extrapolateToPlane(detPlanePtr)
213 except Exception:
214 ut.reportError("shipVertex.TwoTrackVertex: extrapolation did not work")
215 continue
216 mom1 = st1.getMom()
217 mom2 = st2.getMom()
218 cov1 = st1.getCov()
219 cov2 = st2.getCov()
220 cov = ROOT.TMatrixDSym(10)
221 for i in range(10):
222 for j in range(10):
223 if i < 5 and j < 5:
224 cov[i][j] = cov1[i][j]
225 if i > 4 and j > 4:
226 cov[i][j] = cov2[i - 5][j - 5]
227 covInv = ROOT.TMatrixDSym()
228 ROOT.genfit.tools.invertMatrix(cov, covInv)
229
230 self.y_datay_data = np.array([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
231 stVal1 = st1.getState()
232 stVal2 = st2.getState()
233 for i in range(5):
234 self.y_datay_data[i] = stVal1[i]
235 self.y_datay_data[i + 5] = stVal2[i]
236 self.z0z0 = HNLPos[2]
237 self.VyVy = np.zeros(100)
238 for i in range(100):
239 self.VyVy[i] = covInv[i // 10][i % 10]
240
241 np.array([0.0])
242 gMinuit = ROOT.TMinuit(9)
243 tempFcn = self.fcn
244 gMinuit.SetFCN(tempFcn)
245 gMinuit.SetPrintLevel(-1) # minute quiet mode
246 rc = gMinuit.DefineParameter(0, "X pos", HNLPos[0], 0.1, 0, 0)
247 rc = gMinuit.DefineParameter(1, "Y pos", HNLPos[1], 0.1, 0, 0)
248 rc = gMinuit.DefineParameter(2, "Z pos", HNLPos[2], 0.1, 0, 0)
249 rc = gMinuit.DefineParameter(3, "tan1X", mom1[0] / mom1[2], 0.1, 0, 0)
250 rc = gMinuit.DefineParameter(4, "tan1Y", mom1[1] / mom1[2], 0.1, 0, 0)
251 rc = gMinuit.DefineParameter(5, "1/mom1", 1.0 / mom1.Mag(), 0.1, 0, 0)
252 rc = gMinuit.DefineParameter(6, "tan2X", mom2[0] / mom2[2], 0.1, 0, 0)
253 rc = gMinuit.DefineParameter(7, "tan2Y", mom2[1] / mom2[2], 0.1, 0, 0)
254 rc = gMinuit.DefineParameter(8, "1/mom2", 1.0 / mom2.Mag(), 0.1, 0, 0)
255 gMinuit.Clear()
256 gMinuit.Migrad()
257 try:
258 tmp = array("d", [0])
259 err = array("i", [0])
260 gMinuit.mnexcm("HESSE", tmp, -1, err)
261 # gMinuit.mnexcm( "MINOS", tmp, -1, err )
262 except Exception:
263 ut.reportError("shipVertex::minos does not work")
264 continue
265 # get results from TMinuit:
266 emat = array(
267 "d",
268 [
269 0,
270 ]
271 * 81,
272 )
273 gMinuit.mnemat(emat, 9)
274 values = array(
275 "d",
276 [
277 0,
278 ]
279 * 9,
280 )
281 errors = array(
282 "d",
283 [
284 0,
285 ]
286 * 9,
287 )
288 dValue = ctypes.c_double()
289 dError = ctypes.c_double()
290 rc = gMinuit.GetParameter(0, dValue, dError)
291 values[0] = dValue.value
292 errors[0] = dError.value
293 rc = gMinuit.GetParameter(1, dValue, dError)
294 values[1] = dValue.value
295 errors[1] = dError.value
296 rc = gMinuit.GetParameter(2, dValue, dError)
297 values[2] = dValue.value
298 errors[2] = dError.value
299 rc = gMinuit.GetParameter(3, dValue, dError)
300 values[3] = dValue.value
301 errors[3] = dError.value
302 rc = gMinuit.GetParameter(4, dValue, dError)
303 values[4] = dValue.value
304 errors[4] = dError.value
305 rc = gMinuit.GetParameter(5, dValue, dError)
306 values[5] = dValue.value
307 errors[5] = dError.value
308 rc = gMinuit.GetParameter(6, dValue, dError)
309 values[6] = dValue.value
310 errors[6] = dError.value
311 rc = gMinuit.GetParameter(7, dValue, dError)
312 values[7] = dValue.value
313 errors[7] = dError.value
314 rc = gMinuit.GetParameter(8, dValue, dError)
315 values[8] = dValue.value
316 errors[8] = dError.value
317
318 xFit = values[0]
319 yFit = values[1]
320 zFit = values[2]
321 HNLPosFit = ROOT.TVector3(xFit, yFit, zFit)
322 errors[0]
323 errors[1]
324 errors[2]
325
326 # fixme: mass from track reconstraction needed
327 m1 = self.PDG.GetParticle(PosDirCharge[t1]["pdgCode"]).Mass()
328 m2 = self.PDG.GetParticle(PosDirCharge[t2]["pdgCode"]).Mass()
329
330 # self.h['VxpullFit'].Fill( (mctrack.GetStartX()-xFit)/xFitErr )
331 # self.h['VypullFit'].Fill( (mctrack.GetStartY()-yFit)/yFitErr )
332 # self.h['VzpullFit'].Fill( (mctrack.GetStartZ()-zFit)/zFitErr )
333 # self.h['dVxFit'].Fill( (mctrack.GetStartX()-xFit) )
334 # self.h['dVyFit'].Fill( (mctrack.GetStartY()-yFit) )
335 # self.h['dVzFit'].Fill( (mctrack.GetStartZ()-zFit) )
336
337 def getP(fitValues: array[float], cov: array[float], m1, m2):
338 a3 = fitValues[3]
339 a4 = fitValues[4]
340 a5 = fitValues[5]
341 a6 = fitValues[6]
342 a7 = fitValues[7]
343 a8 = fitValues[8]
344 px1 = a3 / (a5 * ROOT.TMath.Sqrt(1 + a3**2 + a4**2))
345 py1 = a4 / (a5 * ROOT.TMath.Sqrt(1 + a3**2 + a4**2))
346 pz1 = 1 / (a5 * ROOT.TMath.Sqrt(1 + a3**2 + a4**2))
347 px2 = a6 / (a8 * ROOT.TMath.Sqrt(1 + a6**2 + a7**2))
348 py2 = a7 / (a8 * ROOT.TMath.Sqrt(1 + a6**2 + a7**2))
349 pz2 = 1 / (a8 * ROOT.TMath.Sqrt(1 + a6**2 + a7**2))
350 Px = px1 + px2
351 Py = py1 + py2
352 Pz = pz1 + pz2
353 E1 = ROOT.TMath.Sqrt(px1**2 + py1**2 + pz1**2 + m1**2)
354 E2 = ROOT.TMath.Sqrt(px2**2 + py2**2 + pz2**2 + m2**2)
355 M = ROOT.TMath.Sqrt(2 * E1 * E2 + m1**2 + m2**2 - 2 * pz1 * pz2 * (1 + a3 * a6 + a4 * a7))
356 MM = 2 * M
357 A5 = 1 + a3**2 + a4**2
358 A8 = 1 + a6**2 + a7**2
359 M_AtoP = ROOT.TMatrixD(4, 6)
360 MT_AtoP = ROOT.TMatrixD(6, 4)
361 covA = ROOT.TMatrixD(6, 6)
362 M_AtoP[0][0] = (1.0 - a3 * a3 / A5) / (a5 * ROOT.TMath.Sqrt(A5))
363 M_AtoP[0][1] = (-a3 * a4 / A5) / (a5 * ROOT.TMath.Sqrt(A5))
364 M_AtoP[0][2] = (-a3) / (a5 * a5 * ROOT.TMath.Sqrt(A5))
365 M_AtoP[0][3] = (1.0 - a6 * a6 / A8) / (a8 * ROOT.TMath.Sqrt(A8))
366 M_AtoP[0][4] = (-a6 * a7 / A8) / (a8 * ROOT.TMath.Sqrt(A8))
367 M_AtoP[0][5] = (-a6) / (a8 * a8 * ROOT.TMath.Sqrt(A8))
368
369 M_AtoP[1][0] = (-a3 * a4 / A5) / (a5 * ROOT.TMath.Sqrt(A5))
370 M_AtoP[1][1] = (1.0 - a4 * a4 / A5) / (a5 * ROOT.TMath.Sqrt(A5))
371 M_AtoP[1][2] = (-a4) / (a5 * a5 * ROOT.TMath.Sqrt(A5))
372 M_AtoP[1][3] = (-a6 * a7 / A8) / (a8 * ROOT.TMath.Sqrt(A8))
373 M_AtoP[1][4] = (1.0 - a7 * a7 / A8) / (a8 * ROOT.TMath.Sqrt(A8))
374 M_AtoP[1][5] = (-a7) / (a8 * a8 * ROOT.TMath.Sqrt(A8))
375
376 M_AtoP[2][0] = (-a3 / A5) / (a5 * ROOT.TMath.Sqrt(A5))
377 M_AtoP[2][1] = (-a4 / A5) / (a5 * ROOT.TMath.Sqrt(A5))
378 M_AtoP[2][2] = (-1.0) / (a5 * a5 * ROOT.TMath.Sqrt(A5))
379 M_AtoP[2][3] = (-a6 / A8) / (a8 * ROOT.TMath.Sqrt(A8))
380 M_AtoP[2][4] = (-a7 / A8) / (a8 * ROOT.TMath.Sqrt(A8))
381 M_AtoP[2][5] = (-1.0) / (a8 * a8 * ROOT.TMath.Sqrt(A8))
382
383 a5a8 = a5 * a8 * ROOT.TMath.Sqrt(A5) * ROOT.TMath.Sqrt(A8)
384
385 M_AtoP[3][0] = (-2 * a6 / a5a8 + 2 * a3 * E2 / (a5 * a5 * A5 * E1)) / MM
386 M_AtoP[3][1] = (-2 * a7 / a5a8 + 2 * a4 * E2 / (a5 * a5 * A5 * E1)) / MM
387 M_AtoP[3][2] = (
388 2 * (1 + a3 * a6 + a4 * a7) / (a5 * a5a8)
389 - 2 * (1.0 + a3 * a3 + a4 * a4) * E2 / (a5 * a5 * a5 * A5 * E1)
390 ) / MM
391 M_AtoP[3][3] = (-2 * a3 / a5a8 + 2 * a6 * E1 / (a8 * a8 * A8 * E2)) / MM
392 M_AtoP[3][4] = (-2 * a4 / a5a8 + 2 * a7 * E1 / (a8 * a8 * A8 * E2)) / MM
393 M_AtoP[3][5] = (
394 2 * (1 + a3 * a6 + a4 * a7) / (a8 * a5a8)
395 - 2 * (1.0 + a6 * a6 + a7 * a7) * E1 / (a8 * a8 * a8 * A8 * E2)
396 ) / MM
397
398 for i in range(4):
399 for j in range(6):
400 MT_AtoP[j][i] = M_AtoP[i][j]
401
402 for i in range(36):
403 covA[i // 6][i % 6] = cov[i // 6 + 3 + (i % 6 + 3) * 9]
404
405 tmp = ROOT.TMatrixD(4, 6)
406 tmp.Mult(M_AtoP, covA)
407 covP = ROOT.TMatrixD(4, 4)
408 covP.Mult(tmp, MT_AtoP)
409 P = ROOT.TLorentzVector()
410 P.SetXYZM(Px, Py, Pz, M)
411 return P, covP
412
413 P, covP = getP(values, emat, m1, m2)
414 # print "******************************************************************************* "
415
416 # create ship particle
417 # notes:
418 # P-TLorentzVector of fitted HNL prticle
419 # covP - covariance matrix of HNL four-vector (Px,Py,Pz,M)
420 # HNLPosFit - fited position of HNL
421 # covV - comariance matrix of the vtx position
422 covV = array("d", [emat[0], emat[1], emat[2], emat[1 + 9], emat[2 + 9], emat[2 + 2 * 9]])
423 covP = array(
424 "d",
425 [
426 covP[0][0],
427 covP[0][1],
428 covP[0][2],
429 covP[0][3],
430 covP[1][1],
431 covP[1][2],
432 covP[1][3],
433 covP[2][2],
434 covP[2][3],
435 covP[3][3],
436 ],
437 )
438
439 # try to make it persistent
440 vx = ROOT.TLorentzVector(
441 HNLPosFit, 0
442 ) # time at vertex still needs to be evaluated from time of tracks and time of flight
443 particle = ROOT.ShipParticle(9900015, 0, -1, -1, t1, t2, P, vx)
444 particle.SetCovV(covV)
445 particle.SetCovP(covP)
446 particle.SetDoca(doca)
447 self.fPartArray.push_back(particle)
448
449 # self.h['dMFit'].Fill( (1.-P.M()) )
450 # self.h['MpullFit'].Fill( (1.-P.M())/ROOT.TMath.Sqrt(covP[3][3]) )
451
452 # self.h['N_Vtx'].Fill(hasVertex)
453
454 def VertexError(self, t1, t2, PosDir, CovMat=None, scalFac=None):
455 # with improved Vx x,y resolution
456 a, u = PosDir[t1]["position"], PosDir[t1]["direction"]
457 c, v = PosDir[t2]["position"], PosDir[t2]["direction"]
458 Vsq = v.Dot(v)
459 Usq = u.Dot(u)
460 UV = u.Dot(v)
461 ca = c - a
462 denom = Usq * Vsq - UV**2
463 tmp2 = Vsq * u - UV * v
464 Va = ca.Dot(tmp2) / denom
465 tmp2 = UV * u - Usq * v
466 Vb = ca.Dot(tmp2) / denom
467 X = (a + c + Va * u + Vb * v) * 0.5
468 l1 = a - X + u * Va # l2 = c - X + v*Vb
469 dist = 2.0 * ROOT.TMath.Sqrt(l1.Dot(l1))
470 if not CovMat:
471 return X, dist
472 T = ROOT.TMatrixD(3, 12)
473 for i in range(3):
474 for k in range(4):
475 for j in range(3):
476 KD = 0
477 if i == j:
478 KD = 1
479 if k == 0 or k == 2:
480 # cova and covc
481 temp = (u[j] * Vsq - v[j] * UV) * u[i] + (u[j] * UV - v[j] * Usq) * v[i]
482 sign = -1
483 if k == 2:
484 sign = +1
485 T[i][3 * k + j] = 0.5 * (KD + sign * temp / denom)
486 elif k == 1:
487 # covu
488 aNAZ = denom * (ca[j] * Vsq - v.Dot(ca) * v[j])
489 aZAN = (ca.Dot(u) * Vsq - ca.Dot(v) * UV) * 2 * (u[j] * Vsq - v[j] * UV)
490 bNAZ = denom * (ca[j] * UV + (u.Dot(ca) * v[j]) - 2 * ca.Dot(v) * u[j])
491 bZAN = (ca.Dot(u) * UV - ca.Dot(v) * Usq) * 2 * (u[j] * Vsq - v[j] * UV)
492 T[i][3 * k + j] = 0.5 * (
493 Va * KD + u[i] / denom**2 * (aNAZ - aZAN) + v[i] / denom**2 * (bNAZ - bZAN)
494 )
495 elif k == 3:
496 # covv
497 aNAZ = denom * (2 * ca.Dot(u) * v[j] - ca.Dot(v) * u[j] - ca[j] * UV)
498 aZAN = (ca.Dot(u) * Vsq - ca.Dot(v) * UV) * 2 * (v[j] * Usq - u[j] * UV)
499 bNAZ = denom * (ca.Dot(u) * u[j] - ca[j] * Usq)
500 bZAN = (ca.Dot(u) * UV - ca.Dot(v) * Usq) * 2 * (v[j] * Usq - u[j] * UV)
501 T[i][3 * k + j] = 0.5 * (
502 Vb * KD + u[i] / denom**2 * (aNAZ - aZAN) + v[i] / denom**2 * (bNAZ - bZAN)
503 )
504 transT = ROOT.TMatrixD(12, 3)
505 transT.Transpose(T)
506 CovTracks = ROOT.TMatrixD(12, 12)
507 tlist = [t1, t2]
508 for k in range(2):
509 for i in range(6):
510 for j in range(6):
511 xfac = 1.0
512 if i > 2:
513 xfac = scalFac[tlist[k]]
514 if j > 2:
515 xfac = xfac * scalFac[tlist[k]]
516 CovTracks[i + k * 6][j + k * 6] = CovMat[tlist[k]][i][j] * xfac
517 # if i==5 or j==5 : CovMat[tlist[k]][i][j] = 0 # ignore error on z-direction
518 tmp = ROOT.TMatrixD(3, 12)
519 tmp.Mult(T, CovTracks)
520 covX = ROOT.TMatrixD(3, 3)
521 covX.Mult(tmp, transT)
522 return X, covX, dist
523
524
525# usage
526# import shipVertex
527# VertexError = shipVertex.Task()
528# newPos,doca = myVertexError(t1,t2,PosDirCharge)
shipVertex.Task.TwoTrackVertex
None TwoTrackVertex(self)
Definition: shipVertex.py:85
shipVertex.Task.fcn
None fcn(self, npar, gin, f, par, iflag)
Definition: shipVertex.py:80
shipVertex.Task.residuals
def residuals(self, y_data, a, int z0)
Definition: shipVertex.py:66
shipVertex.Task.VertexError
def VertexError(self, t1, t2, PosDir, CovMat=None, scalFac=None)
Definition: shipVertex.py:454
shipVertex.Task.__init__
None __init__(self, hp, sTree, mcTree=None)
Definition: shipVertex.py:18
shipVertex.Task.execute
None execute(self)
Definition: shipVertex.py:51
shipVertex.Task.chi2
int chi2(self, res, Vy)
Definition: shipVertex.py:60