SND@LHC Software
Loading...
Searching...
No Matches
Public Member Functions | Public Attributes | Static Public Attributes | List of all members
shipVertex.Task Class Reference

Public Member Functions

 __init__ (self, hp, sTree)
 
 execute (self)
 
 chi2 (self, res, Vy)
 
 residuals (self, y_data, a, z0)
 
 fcn (self, npar, gin, f, par, iflag)
 
 TwoTrackVertex (self)
 
 VertexError (self, t1, t2, PosDir, CovMat=None, scalFac=None)
 

Public Attributes

 sTree
 
 fPartArray
 
 Particles
 
 reps
 
 states
 
 newPosDir
 
 LV
 
 h
 
 PDG
 
 fitTrackLoc
 
 goodTracksLoc
 
 z0
 

Static Public Attributes

 y_data = np.array([0.,0.,0.,0.,0.,0.,0.,0.,0.,0.])
 
int z0 = 0
 
 Vy = np.zeros(100)
 

Detailed Description

Definition at line 13 of file shipVertex.py.

Constructor & Destructor Documentation

◆ __init__()

shipVertex.Task.__init__ (   self,
  hp,
  sTree 
)

Definition at line 15 of file shipVertex.py.

15 def __init__(self,hp,sTree):
16 self.sTree = sTree
17 self.fPartArray = ROOT.TClonesArray("ShipParticle")
18 if not self.sTree.GetBranch("Particles"):
19 self.Particles = self.sTree.Branch("Particles", self.fPartArray,32000,-1)
20 else:
21 self.Particles = self.sTree.Particles
22 self.reps,self.states,self.newPosDir = {},{},{}
23 self.LV={1:ROOT.TLorentzVector(),2:ROOT.TLorentzVector()}
24 self.h = hp
25 self.PDG = ROOT.TDatabasePDG.Instance()
26 self.fitTrackLoc = "FitTracks"
27 self.goodTracksLoc = "goodTracks"
28 #ut.bookHist(self.h,'Vzpull','Vz pull',100,-3.,3.)
29 #ut.bookHist(self.h,'Vxpull','Vx pull',100,-3.,3.)
30 #ut.bookHist(self.h,'Vypull','Vy pull',100,-3.,3.)
31
32 #ut.bookHist(self.h,'dVx','vertex X residual;X^{RECO}-X^{MC}, cm',400,-10.,10.)
33 #ut.bookHist(self.h,'dVy','vertex Y residual;Y^{RECO}-Y^{MC}, cm',400,-10.,10.)
34 #ut.bookHist(self.h,'dVz','vertex Z residual;Z^{RECO}-Z^{MC}, cm',500,-50.,50.)
35 #ut.bookHist(self.h,'VzpullFit','Vz pull, chi2 fit',100,-3.,3.)
36 #ut.bookHist(self.h,'VxpullFit','Vx pull, chi2 fit',100,-3.,3.)
37 #ut.bookHist(self.h,'VypullFit','Vy pull, chi2 fit',100,-3.,3.)
38 #ut.bookHist(self.h,'dVxFit','vertex X residual, chi2 fit;X^{RECO}-X^{MC}, cm',400,-10.,10.)
39 #ut.bookHist(self.h,'dVyFit','vertex Y residual, chi2 fit;Y^{RECO}-Y^{MC}, cm',400,-10.,10.)
40 #ut.bookHist(self.h,'dVzFit','vertex Z residual, chi2 fit;Z^{RECO}-Z^{MC}, cm',500,-50.,50.)
41 #ut.bookHist(self.h,'dMFit','HNL mass resolution, chi2 fit;M^{RECO}-M^{MC}, GeV',2000,-1.,1.)
42 #ut.bookHist(self.h,'MpullFit','M pull, chi2 fit',100,-3.,3.)
43
44 #ut.bookHist(self.h,'N_raveVtx','Number of RAVE vtx',5,-0.5,4.5)
45 #ut.bookHist(self.h,'N_Vtx','Number of vtx',5,-0.5,4.5)
46

Member Function Documentation

◆ chi2()

shipVertex.Task.chi2 (   self,
  res,
  Vy 
)

Definition at line 55 of file shipVertex.py.

55 def chi2(self,res,Vy):
56 s=0
57 for i in range(100):
58 s+=Vy[i]*res[i//10]*res[i%10]
59 return s

◆ execute()

shipVertex.Task.execute (   self)

Definition at line 47 of file shipVertex.py.

47 def execute(self):
48 # make particles persistent
49 self.TwoTrackVertex()
50 self.Particles.Fill()

◆ fcn()

shipVertex.Task.fcn (   self,
  npar,
  gin,
  f,
  par,
  iflag 
)

Definition at line 73 of file shipVertex.py.

73 def fcn(self,npar, gin, f, par, iflag):
74 res = self.residuals(self.y_data,par,self.z0)
75 f = self.chi2(res,self.Vy)
76 return
77

◆ residuals()

shipVertex.Task.residuals (   self,
  y_data,
  a,
  z0 
)

Definition at line 60 of file shipVertex.py.

60 def residuals(self,y_data,a,z0):
61 res = np.zeros(10)
62 res[0] = abs(y_data[0]) - a[5]
63 res[1] = y_data[1] - a[3]
64 res[2] = y_data[2] - a[4]
65 res[3] = y_data[3] - a[0] - a[3]*(a[2] - z0)
66 res[4] = y_data[4] - a[1] - a[4]*(a[2] - z0)
67 res[5] = abs(y_data[5]) - a[8]
68 res[6] = y_data[6] - a[6]
69 res[7] = y_data[7] - a[7]
70 res[8] = y_data[8] - a[0] - a[6]*(a[2] - z0)
71 res[9] = y_data[9] - a[1] - a[7]*(a[2] - z0)
72 return res

◆ TwoTrackVertex()

shipVertex.Task.TwoTrackVertex (   self)

Definition at line 78 of file shipVertex.py.

78 def TwoTrackVertex(self):
79 self.fPartArray.Delete()
80 fittedTracks = getattr(self.sTree,self.fitTrackLoc)
81 goodTracks = getattr(self.sTree,self.goodTracksLoc)
82 if goodTracks.size() < 2: return
83 particles = self.fPartArray
84 PosDirCharge,CovMat,scalFac = {},{},{}
85 for tr in goodTracks:
86 fitStatus = fittedTracks[tr].getFitStatus()
87 xx = fittedTracks[tr].getFittedState()
88 pid = xx.getPDG()
89 if not global_variables.pidProton and abs(pid) == 2212:
90 pid = int(math.copysign(211,pid))
91 rep = ROOT.genfit.RKTrackRep(xx.getPDG())
92 state = ROOT.genfit.StateOnPlane(rep)
93 rep.setPosMom(state,xx.getPos(),xx.getMom())
94 PosDirCharge[tr] = {'position':xx.getPos(),'direction':xx.getDir(),\
95 'momentum':xx.getMom(),'charge':xx.getCharge(),'pdgCode':pid,'state':xx,'rep':rep,'newstate':state}
96 CovMat[tr] = xx.get6DCov()
97#
98 if len(PosDirCharge) < 2: return
99 if len(PosDirCharge) > 4: return # abort too busy events
100 for t1 in PosDirCharge:
101 c1 = PosDirCharge[t1]['charge']
102 for t2 in PosDirCharge:
103 if not t2>t1: continue
104 # ignore this for background studies
105 if PosDirCharge[t2]['charge'] == c1 : continue
106 newPos,doca = self.VertexError(t1,t2,PosDirCharge)
107# as we have learned, need iterative procedure
108 dz = 99999.
109 rc = True
110 step = 0
111 while dz > 0.01:
112 zBefore = newPos[2]
113 # make a new rep for track 1,2
114 for tr in [t1,t2]:
115 try:
116 PosDirCharge[tr]['rep'].extrapolateToPoint(PosDirCharge[tr]['newstate'], newPos, False)
117 except:
118 ut.reportError('shipVertex: extrapolation did not worked')
119 rc = False
120 break
121 self.newPosDir[tr] = {'position':PosDirCharge[tr]['rep'].getPos(PosDirCharge[tr]['newstate']),\
122 'direction':PosDirCharge[tr]['rep'].getDir(PosDirCharge[tr]['newstate']),\
123 'momentum':PosDirCharge[tr]['rep'].getMom(PosDirCharge[tr]['newstate'])}
124 if not rc: break
125 newPos,doca = self.VertexError(t1,t2,self.newPosDir)
126 dz = abs(zBefore-newPos[2])
127 step+=1
128 if step > 10:
129 ut.reportError("shipVertex::abort iteration, too many steps")
130 if global_variables.debug:
131 print('abort iteration, too many steps, pos=',newPos[0],newPos[1],newPos[2],' doca=',doca,'z before and dz',zBefore,dz)
132 rc = False
133 break
134#
135 if not rc: continue # extrapolation failed, makes no sense to continue
136# now go for the last step and vertex error
137 scalFac[t1] = (PosDirCharge[t1]['position'][2]-newPos[2])/PosDirCharge[t1]['direction'][2]/PosDirCharge[t1]['momentum'].Mag()
138 scalFac[t2] = (PosDirCharge[t2]['position'][2]-newPos[2])/PosDirCharge[t2]['direction'][2]/PosDirCharge[t2]['momentum'].Mag()
139 HNLPos,covX,dist = self.VertexError(t1,t2,self.newPosDir,CovMat,scalFac)
140# monitor Vx resolution and pulls
141 #print "DEBUG",HNLPos[0],HNLPos[1],HNLPos[2],dist,covX[0][0],covX[1][1],covX[2][2]
142 #print " ",mctrack.GetStartX(),mctrack.GetStartY(),mctrack.GetStartZ()
143# HNL true
144 if self.sTree.GetBranch("fitTrack2MC"):
145 mctrack = self.sTree.MCTrack[self.sTree.fitTrack2MC[t1]]
146 mctrack2 = self.sTree.MCTrack[self.sTree.fitTrack2MC[t2]]
147 mcHNL = self.sTree.MCTrack[mctrack.GetMotherId()]
148 #print "true vtx: ",mctrack.GetStartX(),mctrack.GetStartY(),mctrack.GetStartZ()
149 #print "reco vtx: ",HNLPos[0],HNLPos[1],HNLPos[2]
150 #self.h['Vzpull'].Fill( (mctrack.GetStartZ()-HNLPos[2])/ROOT.TMath.Sqrt(covX[2][2]) )
151 #self.h['Vxpull'].Fill( (mctrack.GetStartX()-HNLPos[0])/ROOT.TMath.Sqrt(covX[0][0]) )
152 #self.h['Vypull'].Fill( (mctrack.GetStartY()-HNLPos[1])/ROOT.TMath.Sqrt(covX[1][1]) )
153 #self.h['dVx'].Fill( (mctrack.GetStartX()-HNLPos[0]) )
154 #self.h['dVy'].Fill( (mctrack.GetStartY()-HNLPos[1]) )
155 #self.h['dVz'].Fill( (mctrack.GetStartZ()-HNLPos[2]) )
156
157
158 #print "*********************************** vertex fit precise ******************************************** "
159
160 detPlane = ROOT.genfit.DetPlane(ROOT.TVector3(0,0,HNLPos[2]),ROOT.TVector3(1,0,0),ROOT.TVector3(0,1,0))
161 plane = ROOT.genfit.RKTrackRep().makePlane(ROOT.TVector3(0,0,HNLPos[2]),ROOT.TVector3(1,0,0),ROOT.TVector3(0,1,0))
162 st1 = fittedTracks[t1].getFittedState()
163 st2 = fittedTracks[t2].getFittedState()
164 try:
165 st1.extrapolateToPlane(plane)
166 except:
167 ut.reportError("shipVertex::TwoTrackVertex: extrapolation does not worked")
168 continue
169 try:
170 st2.extrapolateToPlane(plane)
171 except:
172 ut.reportError("shipVertex::TwoTrackVertex: extrapolation does not worked")
173 continue
174 mom1 = st1.getMom()
175 mom2 = st2.getMom()
176 cov1 = st1.getCov()
177 cov2 = st2.getCov()
178 cov = ROOT.TMatrixDSym(10)
179 for i in range(10):
180 for j in range(10):
181 if i<5 and j<5: cov[i][j]=cov1[i][j]
182 if i>4 and j>4: cov[i][j]=cov2[i-5][j-5]
183 covInv = ROOT.TMatrixDSym()
184 ROOT.genfit.tools.invertMatrix(cov,covInv)
185
186 self.y_data = np.array([0.,0.,0.,0.,0.,0.,0.,0.,0.,0.])
187 stVal1 = st1.getState()
188 stVal2 = st2.getState()
189 for i in range(5):
190 self.y_data[i]=stVal1[i]
191 self.y_data[i+5]=stVal2[i]
192 self.z0 = HNLPos[2]
193 self.Vy = np.zeros(100)
194 for i in range(100):
195 self.Vy[i] = covInv[i//10][i%10]
196
197 f=np.array([0.])
198 gMinuit = ROOT.TMinuit(9)
199 tempFcn = self.fcn
200 gMinuit.SetFCN(tempFcn)
201 gMinuit.SetPrintLevel(-1)#minute quiet mode
202 rc = gMinuit.DefineParameter(0, 'X pos',HNLPos[0], 0.1,0,0)
203 rc = gMinuit.DefineParameter(1, 'Y pos',HNLPos[1], 0.1,0,0)
204 rc = gMinuit.DefineParameter(2, 'Z pos',HNLPos[2], 0.1,0,0)
205 rc = gMinuit.DefineParameter(3, 'tan1X',mom1[0]/mom1[2], 0.1,0,0)
206 rc = gMinuit.DefineParameter(4, 'tan1Y',mom1[1]/mom1[2], 0.1,0,0)
207 rc = gMinuit.DefineParameter(5, '1/mom1',1./mom1.Mag(), 0.1,0,0)
208 rc = gMinuit.DefineParameter(6, 'tan2X',mom2[0]/mom2[2], 0.1,0,0)
209 rc = gMinuit.DefineParameter(7, 'tan2Y',mom2[1]/mom2[2], 0.1,0,0)
210 rc = gMinuit.DefineParameter(8, '1/mom2',1./mom2.Mag(), 0.1,0,0)
211 gMinuit.Clear()
212 gMinuit.Migrad()
213 try:
214 tmp = array('d',[0])
215 err = array('i',[0])
216 gMinuit.mnexcm( "HESSE", tmp, -1, err )
217 #gMinuit.mnexcm( "MINOS", tmp, -1, err )
218 except:
219 ut.reportError("shipVertex::minos does not work")
220 continue
221 #get results from TMinuit:
222 emat = array('d',[0,]*81)
223 gMinuit.mnemat(emat,9)
224 values = array('d',[0,]*9)
225 errors = array('d',[0,]*9)
226 dValue = ctypes.c_double()
227 dError = ctypes.c_double()
228 rc = gMinuit.GetParameter(0, dValue, dError)
229 values[0]=dValue.value
230 errors[0]=dError.value
231 rc = gMinuit.GetParameter(1, dValue, dError)
232 values[1]=dValue.value
233 errors[1]=dError.value
234 rc = gMinuit.GetParameter(2, dValue, dError)
235 values[2]=dValue.value
236 errors[2]=dError.value
237 rc = gMinuit.GetParameter(3, dValue, dError)
238 values[3]=dValue.value
239 errors[3]=dError.value
240 rc = gMinuit.GetParameter(4, dValue, dError)
241 values[4]=dValue.value
242 errors[4]=dError.value
243 rc = gMinuit.GetParameter(5, dValue, dError)
244 values[5]=dValue.value
245 errors[5]=dError.value
246 rc = gMinuit.GetParameter(6, dValue, dError)
247 values[6]=dValue.value
248 errors[6]=dError.value
249 rc = gMinuit.GetParameter(7, dValue, dError)
250 values[7]=dValue.value
251 errors[7]=dError.value
252 rc = gMinuit.GetParameter(8, dValue, dError)
253 values[8]=dValue.value
254 errors[8]=dError.value
255
256 xFit = values[0]
257 yFit = values[1]
258 zFit = values[2]
259 HNLPosFit = ROOT.TVector3(xFit,yFit,zFit)
260 xFitErr = errors[0]
261 yFitErr = errors[1]
262 zFitErr = errors[2]
263
264 #fixme: mass from track reconstraction needed
265 m1 = self.PDG.GetParticle(PosDirCharge[t1]['pdgCode']).Mass()
266 m2 = self.PDG.GetParticle(PosDirCharge[t2]['pdgCode']).Mass()
267
268 #self.h['VxpullFit'].Fill( (mctrack.GetStartX()-xFit)/xFitErr )
269 #self.h['VypullFit'].Fill( (mctrack.GetStartY()-yFit)/yFitErr )
270 #self.h['VzpullFit'].Fill( (mctrack.GetStartZ()-zFit)/zFitErr )
271 #self.h['dVxFit'].Fill( (mctrack.GetStartX()-xFit) )
272 #self.h['dVyFit'].Fill( (mctrack.GetStartY()-yFit) )
273 #self.h['dVzFit'].Fill( (mctrack.GetStartZ()-zFit) )
274
275 def getP(fitValues,cov,m1,m2):
276 a3=fitValues[3]
277 a4=fitValues[4]
278 a5=fitValues[5]
279 a6=fitValues[6]
280 a7=fitValues[7]
281 a8=fitValues[8]
282 px1 = a3/(a5*ROOT.TMath.Sqrt(1 + a3**2 + a4**2))
283 py1 = a4/(a5*ROOT.TMath.Sqrt(1 + a3**2 + a4**2))
284 pz1 = 1/(a5*ROOT.TMath.Sqrt(1 + a3**2 + a4**2))
285 px2 = a6/(a8*ROOT.TMath.Sqrt(1 + a6**2 + a7**2))
286 py2 = a7/(a8*ROOT.TMath.Sqrt(1 + a6**2 + a7**2))
287 pz2 = 1/(a8*ROOT.TMath.Sqrt(1 + a6**2 + a7**2))
288 Px = px1 + px2
289 Py = py1 + py2
290 Pz = pz1 + pz2
291 E1 = ROOT.TMath.Sqrt(px1**2 + py1**2 + pz1**2 + m1**2)
292 E2 = ROOT.TMath.Sqrt(px2**2 + py2**2 + pz2**2 + m2**2)
293 M = ROOT.TMath.Sqrt(2*E1*E2+m1**2+m2**2-2*pz1*pz2*(1+a3*a6+a4*a7))
294 MM = 2*M
295 A5 = 1 + a3**2 + a4**2
296 A8 = 1 + a6**2 + a7**2
297 M_AtoP = ROOT.TMatrixD(4,6)
298 MT_AtoP = ROOT.TMatrixD(6,4)
299 covA = ROOT.TMatrixD(6,6)
300 M_AtoP[0][0] = (1.-a3*a3/A5)/(a5*ROOT.TMath.Sqrt(A5))
301 M_AtoP[0][1] = (-a3*a4/A5)/(a5*ROOT.TMath.Sqrt(A5))
302 M_AtoP[0][2] = (-a3)/(a5*a5*ROOT.TMath.Sqrt(A5))
303 M_AtoP[0][3] = (1.-a6*a6/A8)/(a8*ROOT.TMath.Sqrt(A8))
304 M_AtoP[0][4] = (-a6*a7/A8)/(a8*ROOT.TMath.Sqrt(A8))
305 M_AtoP[0][5] = (-a6)/(a8*a8*ROOT.TMath.Sqrt(A8))
306
307 M_AtoP[1][0] = (-a3*a4/A5)/(a5*ROOT.TMath.Sqrt(A5))
308 M_AtoP[1][1] = (1.-a4*a4/A5)/(a5*ROOT.TMath.Sqrt(A5))
309 M_AtoP[1][2] = (-a4)/(a5*a5*ROOT.TMath.Sqrt(A5))
310 M_AtoP[1][3] = (-a6*a7/A8)/(a8*ROOT.TMath.Sqrt(A8))
311 M_AtoP[1][4] = (1.-a7*a7/A8)/(a8*ROOT.TMath.Sqrt(A8))
312 M_AtoP[1][5] = (-a7)/(a8*a8*ROOT.TMath.Sqrt(A8))
313
314 M_AtoP[2][0] = (-a3/A5)/(a5*ROOT.TMath.Sqrt(A5))
315 M_AtoP[2][1] = (-a4/A5)/(a5*ROOT.TMath.Sqrt(A5))
316 M_AtoP[2][2] = (-1.)/(a5*a5*ROOT.TMath.Sqrt(A5))
317 M_AtoP[2][3] = (-a6/A8)/(a8*ROOT.TMath.Sqrt(A8))
318 M_AtoP[2][4] = (-a7/A8)/(a8*ROOT.TMath.Sqrt(A8))
319 M_AtoP[2][5] = (-1.)/(a8*a8*ROOT.TMath.Sqrt(A8))
320
321 a5a8 = a5*a8*ROOT.TMath.Sqrt(A5)*ROOT.TMath.Sqrt(A8)
322
323 M_AtoP[3][0] = (-2*a6/a5a8 + 2*a3*E2/(a5*a5*A5*E1))/MM
324 M_AtoP[3][1] = (-2*a7/a5a8 + 2*a4*E2/(a5*a5*A5*E1))/MM
325 M_AtoP[3][2] = (2*(1+a3*a6+a4*a7)/(a5*a5a8) - 2*(1.+a3*a3+a4*a4)*E2/(a5*a5*a5*A5*E1))/MM
326 M_AtoP[3][3] = (-2*a3/a5a8 + 2*a6*E1/(a8*a8*A8*E2))/MM
327 M_AtoP[3][4] = (-2*a4/a5a8 + 2*a7*E1/(a8*a8*A8*E2))/MM
328 M_AtoP[3][5] = (2*(1+a3*a6+a4*a7)/(a8*a5a8) - 2*(1.+a6*a6+a7*a7)*E1/(a8*a8*a8*A8*E2))/MM
329
330 for i in range(4):
331 for j in range(6):
332 MT_AtoP[j][i] = M_AtoP[i][j]
333
334 for i in range(36):
335 covA[i//6][i%6] = cov[i//6+3+(i%6+3)*9]
336
337 tmp = ROOT.TMatrixD(4,6)
338 tmp.Mult(M_AtoP,covA)
339 covP = ROOT.TMatrixD(4,4)
340 covP.Mult(tmp,MT_AtoP)
341 P = ROOT.TLorentzVector()
342 P.SetXYZM(Px,Py,Pz,M)
343 return P,covP
344
345 P,covP = getP(values,emat,m1,m2)
346 #print "******************************************************************************* "
347
348 #create ship particle
349 #notes:
350 #P-TLorentzVector of fitted HNL prticle
351 #covP - covariance matrix of HNL four-vector (Px,Py,Pz,M)
352 #HNLPosFit - fited position of HNL
353 #covV - comariance matrix of the vtx position
354 covV = array('d',[emat[0],emat[1],emat[2],emat[1+9],emat[2+9],emat[2+2*9]])
355 covP = array('d',[covP[0][0],covP[0][1],covP[0][2],covP[0][3],covP[1][1],covP[1][2],covP[1][3],covP[2][2],covP[2][3],covP[3][3]])
356
357# try to make it persistent
358 vx = ROOT.TLorentzVector(HNLPosFit,0) # time at vertex still needs to be evaluated from time of tracks and time of flight
359 particle = ROOT.ShipParticle(9900015,0,-1,-1,t1,t2,P,vx)
360 particle.SetCovV(covV)
361 particle.SetCovP(covP)
362 particle.SetDoca(doca)
363 nParts = particles.GetEntries()
364 particles[nParts] = particle
365
366 #self.h['dMFit'].Fill( (1.-P.M()) )
367 #self.h['MpullFit'].Fill( (1.-P.M())/ROOT.TMath.Sqrt(covP[3][3]) )
368
369 #self.h['N_Vtx'].Fill(hasVertex)
370

◆ VertexError()

shipVertex.Task.VertexError (   self,
  t1,
  t2,
  PosDir,
  CovMat = None,
  scalFac = None 
)

Definition at line 371 of file shipVertex.py.

371 def VertexError(self,t1,t2,PosDir,CovMat=None,scalFac=None):
372# with improved Vx x,y resolution
373 a,u = PosDir[t1]['position'],PosDir[t1]['direction']
374 c,v = PosDir[t2]['position'],PosDir[t2]['direction']
375 Vsq = v.Dot(v)
376 Usq = u.Dot(u)
377 UV = u.Dot(v)
378 ca = c-a
379 denom = Usq*Vsq-UV**2
380 tmp2 = Vsq*u-UV*v
381 Va = ca.Dot(tmp2)/denom
382 tmp2 = UV*u-Usq*v
383 Vb = ca.Dot(tmp2)/denom
384 X = (a+c+Va*u+Vb*v) * 0.5
385 l1 = a - X + u*Va # l2 = c - X + v*Vb
386 dist = 2. * ROOT.TMath.Sqrt( l1.Dot(l1) )
387 if not CovMat: return X,dist
388 T = ROOT.TMatrixD(3,12)
389 for i in range(3):
390 for k in range(4):
391 for j in range(3):
392 KD = 0
393 if i==j: KD = 1
394 if k==0 or k==2:
395 # cova and covc
396 temp = ( u[j]*Vsq - v[j]*UV )*u[i] + (u[j]*UV-v[j]*Usq)*v[i]
397 sign = -1
398 if k==2 : sign = +1
399 T[i][3*k+j] = 0.5*( KD + sign*temp/denom )
400 elif k==1:
401 # covu
402 aNAZ = denom*( ca[j]*Vsq-v.Dot(ca)*v[j] )
403 aZAN = ( ca.Dot(u)*Vsq-ca.Dot(v)*UV )*2*( u[j]*Vsq-v[j]*UV )
404 bNAZ = denom*( ca[j]*UV+(u.Dot(ca)*v[j]) - 2*ca.Dot(v)*u[j] )
405 bZAN = ( ca.Dot(u)*UV-ca.Dot(v)*Usq )*2*( u[j]*Vsq-v[j]*UV )
406 T[i][3*k+j] = 0.5*( Va*KD + u[i]/denom**2*(aNAZ-aZAN) + v[i]/denom**2*(bNAZ-bZAN) )
407 elif k==3:
408 # covv
409 aNAZ = denom*( 2*ca.Dot(u)*v[j] - ca.Dot(v)*u[j] - ca[j]*UV )
410 aZAN = ( ca.Dot(u)*Vsq-ca.Dot(v)*UV )*2*( v[j]*Usq-u[j]*UV )
411 bNAZ = denom*( ca.Dot(u)*u[j]-ca[j]*Usq )
412 bZAN = ( ca.Dot(u)*UV-ca.Dot(v)*Usq )*2*( v[j]*Usq-u[j]*UV )
413 T[i][3*k+j] = 0.5*(Vb*KD + u[i]/denom**2*(aNAZ-aZAN) + v[i]/denom**2*(bNAZ-bZAN) )
414 transT = ROOT.TMatrixD(12,3)
415 transT.Transpose(T)
416 CovTracks = ROOT.TMatrixD(12,12)
417 tlist = [t1,t2]
418 for k in range(2):
419 for i in range(6):
420 for j in range(6):
421 xfac = 1.
422 if i>2: xfac = scalFac[tlist[k]]
423 if j>2: xfac = xfac * scalFac[tlist[k]]
424 CovTracks[i+k*6][j+k*6] = CovMat[tlist[k]][i][j] * xfac
425 # if i==5 or j==5 : CovMat[tlist[k]][i][j] = 0 # ignore error on z-direction
426 tmp = ROOT.TMatrixD(3,12)
427 tmp.Mult(T,CovTracks)
428 covX = ROOT.TMatrixD(3,3)
429 covX.Mult(tmp,transT)
430 return X,covX,dist
431
432#usage
433# import shipVertex
434# VertexError = shipVertex.Task()
435# newPos,doca = myVertexError(t1,t2,PosDirCharge)

Member Data Documentation

◆ fitTrackLoc

shipVertex.Task.fitTrackLoc

Definition at line 26 of file shipVertex.py.

◆ fPartArray

shipVertex.Task.fPartArray

Definition at line 17 of file shipVertex.py.

◆ goodTracksLoc

shipVertex.Task.goodTracksLoc

Definition at line 27 of file shipVertex.py.

◆ h

shipVertex.Task.h

Definition at line 24 of file shipVertex.py.

◆ LV

shipVertex.Task.LV

Definition at line 23 of file shipVertex.py.

◆ newPosDir

shipVertex.Task.newPosDir

Definition at line 22 of file shipVertex.py.

◆ Particles

shipVertex.Task.Particles

Definition at line 19 of file shipVertex.py.

◆ PDG

shipVertex.Task.PDG

Definition at line 25 of file shipVertex.py.

◆ reps

shipVertex.Task.reps

Definition at line 22 of file shipVertex.py.

◆ states

shipVertex.Task.states

Definition at line 22 of file shipVertex.py.

◆ sTree

shipVertex.Task.sTree

Definition at line 16 of file shipVertex.py.

◆ Vy

shipVertex.Task.Vy = np.zeros(100)
static

Definition at line 54 of file shipVertex.py.

◆ y_data

shipVertex.Task.y_data = np.array([0.,0.,0.,0.,0.,0.,0.,0.,0.,0.])
static

Definition at line 52 of file shipVertex.py.

◆ z0 [1/2]

int shipVertex.Task.z0 = 0
static

Definition at line 53 of file shipVertex.py.

◆ z0 [2/2]

shipVertex.Task.z0

Definition at line 192 of file shipVertex.py.

The documentation for this class was generated from the following file: