import ROOT,os,sys,getopt import rootUtils as ut import shipunit as u class FitTrackInfo(object): def __init__(self, tree, debug=0): self.tree = tree self.debug = debug self.count = 0 self.Momentum = {} self.Direction = {} self.Position = {} self.__info = {} self.Vertex = None self.Doca = None # 0 orig, 1 refit, -1 failed refit self.vertexEFlag = 0 def clean(self): self.count = 0 self.Momentum.clear() self.Direction.clear() self.Position.clear() self.__info.clear() self.Vertex = None self.Doca = None def Print(self): print "FitTrackInfo:" for tid in self.__info: print "\t", tid, "{:10.4f}".format(self.__info[tid]['Ndf']), "{:10.4f}".format(self.__info[tid]['Chi2']), print " pos:", " ".join("{:10.4f}".format(self.Position[tid](ii)) for ii in range(0,3)), print " mom:", " ".join("{:10.4f}".format(self.Direction[tid](ii)*self.Momentum[tid]) for ii in range(0,3)), print " P:", "{:10.4f}".format(self.Momentum[tid]) ## \brief returns list of keys #__info . # \return list of MCtrackIDs of fitted tracks. def getTrIDs(self): trID = [] for tid in self.__info: trID.append(tid) return trID def getNtracks(self) : return len(self.__info) def getChi2Ndf(self, tid): return self.__info[tid]['Chi2']/self.__info[tid]['Ndf'] def getNdf(self, tid): return self.__info[tid]['Ndf'] def compareTracks(self, tid, Pos, Dir, Pval): false2 = (Pos==None or Dir==None or Pval==None) false1 = (not tid in self.__info) if (false2 and false1): return True if (false2 or false1): return False return ( (self.Direction[tid]==Dir) and (self.Position[tid]==Pos) and (self.Momentum[tid]==Pval) ) def getPosDirPval(self, tid): if not tid in self.__info: return None, None, None return self.Position[tid], self.Direction[tid], self.Momentum[tid] def getVertex(self) : return self.Vertex def getDoca(self): return self.Doca def myVertex2(self, t1,t2): deltaPos =(self.Position[t1]-self.Position[t2]) # A1-A2 dotDir = self.Direction[t1].Dot(self.Direction[t2]) # a1.a2 crossDir = self.Direction[t1].Cross(self.Direction[t2]) # a1xa2 uPerpend = crossDir*(1./crossDir.Mag()) # (a1xa2)/|a1xa2| from a1 to a2 minDist = deltaPos.Dot(uPerpend) # (A1-A2).(a1xa2)/|a1xa2| # A1 + a1*t1 + (minDist * uPerpend) is (A1 + a1*t1) projected to the plane: # 1) A2+a2*t2 belons to the plane, # 2) A1+a1*t1 is parallel to the plane # cross at t1,t2: A1+a1*t1+(minDist*uPerpend) = A2+a2*t2 t2X = self.Direction[t2].X() if (t2X == 0) : t2X = 0.00000000001 a2a = self.Direction[t2].Y()/t2X alpha = deltaPos - minDist*uPerpend nomin = alpha.Y() - a2a*alpha.X() denom = a2a*self.Direction[t1].X() - self.Direction[t1].Y() s1 = nomin/denom s2 = ( alpha.X() + self.Direction[t1].X()*s1 ) / t2X#self.Direction[t2].X() vec1 = self.Position[t1]+s1*self.Direction[t1] vec2 = self.Position[t2]+s2*self.Direction[t2] ave = (vec1+vec2)*0.5 dif = vec1-vec2 debugNeed = False if(abs(abs(minDist)-dif.Mag())>0.00000001): print "myVertex2 - problem:" debugNeed = True if(self.debug>2 or debugNeed): for tid in (t1,t2): print str(tid)+": Pos : ", "".join(str(self.Position[tid](ii)) +" " for ii in range(0,3)), print "\t\tMom : ", "".join(str(self.Direction[tid](ii))+" " for ii in range(0,3)) print "uPerpend: ","".join(str(uPerpend(ii))+" " for ii in range(0,3)) if(self.debug>1 or debugNeed): print "fit vertex: -> 1st poing : ", vec1.X(), vec1.Y(), vec1.Z() print "fit vertex: -> 2nd point : ", vec2.X(), vec2.Y(), vec2.Z() print "fit vertex: -> average : ", ave.X(), ave.Y(), ave.Z() print "distance", abs(minDist), dif.Mag() return ave, abs(minDist) def readEvent(self): self.clean() indx = -1 for atrack in self.tree.FitTracks: # kill tracks outside fiducial volume # if not checkFiducialVolume(sTree,key,dy): continue fitStatus = atrack.getFitStatus() if not fitStatus.isFitConverged() : continue indx+=1 mcTrID = self.tree.fitTrack2MC[indx] self.__info[mcTrID] = {} self.__info[mcTrID]['Ndf'] = fitStatus.getNdf() self.__info[mcTrID]['Chi2'] = fitStatus.getChi2() fittedState = atrack.getFittedState() self.Momentum[mcTrID] = fittedState.getMomMag() self.Direction[mcTrID] = fittedState.getDir() self.Position[mcTrID] = fittedState.getPos() if(indx>0): if(indx==1): self.createVertex(self.__info.keys()[0], self.__info.keys()[1]) else: pass #print "More than 2 fitterd tracks" return len(self.__info) def createVertex(self, tid1, tid2, flag=0): if( (tid1 in self.__info) and (tid2 in self.__info) ): self.Vertex, self.Doca = self.myVertex2(tid1, tid2) self.vertexEFlag = flag def addNewTrack(self, mcTrID, position, direction, Pval, ndf, chi2, verb = True): if (verb and mcTrID in self.Momentum): print "FotTrackInfo WARNING - trID ", mcTrID, "already filled! Will rewrite all records for this trID!" self.__info[mcTrID]={} self.__info[mcTrID]['Ndf'] = ndf self.__info[mcTrID]['Chi2'] = chi2 self.Momentum[mcTrID] = Pval self.Direction[mcTrID] = direction self.Position[mcTrID] = position