import ROOT,argparse,random,math #parser=argparse.ArgumentParser(description='test the covarianz matrix') #parser.add_argument( ROOT.gROOT.ProcessLine(".x rootlogon.C") ROOT.gROOT.ProcessLine('gROOT->SetStyle("Plain")') ROOT.gROOT.LoadMacro("stlPYROOT.h+") #set_palette() #pcl=ROOT.TpcPrelimCluster() digis=[] digiamp=[] #for i in range(5): # digis.append(ROOT.TVector3(random.randrange(0,5),random.randrange(0,5),0))#random.randrange(0,5))) # digiamp.append(random.randrange(1,2)) digis.append(ROOT.TVector3(0,0,1)) digiamp.append(1) #digis.append(ROOT.TVector3(0,0,1)) #digiamp.append(1) #digis.append(ROOT.TVector3(0,0,-1)) #digiamp.append(1) digis.append(ROOT.TVector3(0,0,-1)) digiamp.append(1) digis.append(ROOT.TVector3(0,1,0)) digiamp.append(1) #digis.append(ROOT.TVector3(1,0,0)) #digiamp.append(1) digis.append(ROOT.TVector3(0,-1,0)) digiamp.append(1) #digis.append(ROOT.TVector3(-1,0,0)) #digiamp.append(1) #digis.append(ROOT.TVector3(1,1,0)) #digiamp.append(1) #digis.append(ROOT.TVector3(-1,-1,0)) #digiamp.append(1) digis.append(ROOT.TVector3(0,1,1)) digiamp.append(1) #digis.append(ROOT.TVector3(-1,-1,0)) #digiamp.append(1) #digis.append(ROOT.TVector3(1,1,1)) #digiamp.append(1) #digis.append(ROOT.TVector3(-1,-1,-1)) #digiamp.append(1) clpos=ROOT.TVector3() clamp=0 dc=-1 for dig in digis: dc+=1 clpos+=ROOT.TVector3(dig.X()*digiamp[dc],dig.Y()*digiamp[dc],dig.Z()*digiamp[dc]) clamp+=digiamp[dc] #clpos*=1./len(digis) clpos*=1./clamp #clpos=ROOT.TVector3(0,0,0) sig=ROOT.TVector3(0,0,0) dc=-1 for dig in digis: dc+=1 df=dig-clpos sigx=df.X()**2*digiamp[dc] sigy=df.Y()**2*digiamp[dc] sigz=df.Z()**2*digiamp[dc] thissig=ROOT.TVector3(sigx,sigy,sigz) sig+=thissig #sig*=1./len(digis) sig*=1./clamp #clpos=ROOT.TVector3(0,0,0) cl=ROOT.TpcCluster(clpos,sig,len(digis),0) print'clusterpos:' clpos.Print() print 'cluster sig:' sig.Print() cov=ROOT.TMatrixD(3,3) dc=-1 for dig in digis: dc+=1 c=ROOT.TMatrixD(3,1) for i in range(3): c[i][0]=(dig[i]-clpos[i]) c_t=ROOT.TMatrixD(ROOT.TMatrixD.kTransposed,c) acov=ROOT.TMatrixD(c,ROOT.TMatrixD.kMult,c_t) acov*=float(digiamp[dc]) # c.Print() # c_t.Print() # acov.Print() cov+=acov print '***********************' #cov*=1./len(digis) cov*=1./clamp #test cov: for dig in digis: dc+=1 c=ROOT.TMatrixD(3,1) for i in range(3): c[i][0]=(dig[i]-clpos[i]) c_t=ROOT.TMatrixD(ROOT.TMatrixD.kTransposed,c) mult1=ROOT.TMatrixD(cov,ROOT.TMatrixD.kMult,c) num=ROOT.TMatrixD(c_t,ROOT.TMatrixD.kMult,mult1) print 'cov test:',num[0][0] eigenproblem=ROOT.TMatrixDEigen(cov) eigenwert=eigenproblem.GetEigenValues() eigenvector=eigenproblem.GetEigenVectors() eigenvector_i=ROOT.TMatrixD(ROOT.TMatrixD.kInverted,eigenvector) print 'eigenvalues:' eigenwert.Print() #calculate ellipse: S=ROOT.TMatrixD(3,3) for i in range(0,3): for j in range(0,3): if i==j: if eigenwert[i][j]!=0: S[i][j]=1./(eigenwert[i][j]**2) S[i][j]=eigenwert[i][j] else: S[i][j]=0 else: S[i][j]=0 print 'S matrix:' S.Print() temp=ROOT.TMatrixD(S,ROOT.TMatrixD.kMult,eigenvector_i) ellipse=ROOT.TMatrixD(eigenvector,ROOT.TMatrixD.kMult,temp) print 'ellipse:' ellipse.Print() print 'cof:' cov.Print() ellipse_eigenprob=ROOT.TMatrixDEigen(ellipse) ellipse_eigenwert=ellipse_eigenprob.GetEigenValues() ellipse_eigenvect=ellipse_eigenprob.GetEigenVectors() print 'ellipse eigenwert:' ellipse_eigenwert.Print() print 'ellipse eigenvect:' ellipse_eigenvect.Print() print '1/e^2:' S.Print() print 'eigenvect:' eigenvector.Print() #xaxis=ROOT.TVector3(1,0,0) print 'intersection calculation: start ************************' for i in range(0,4): c=ROOT.TMatrixD(3,1) if i<3: for j in range(0,3): print i,j if i==j: c[j][0]=1 else: c[j][0]=0 else: c[0][0]=0 c[1][0]=1 c[2][0]=1 c.Print() c_t=ROOT.TMatrixD(ROOT.TMatrixD.kTransposed,c) mult1=ROOT.TMatrixD(ellipse,ROOT.TMatrixD.kMult,c) num=ROOT.TMatrixD(c_t,ROOT.TMatrixD.kMult,mult1) if num[0][0]!=0: print 'intersection:',i,1/math.sqrt(num[0][0]) print 'intersection calculation: end **************************' print 'the covarianz:' cov.Print() cl.SetCov(cov) axis=[] for i in range(3): ax=cl.axis(i) print 'axis',i,'of cov:' cl.axis(i).Print() thisaxis=ROOT.TMatrixD(3,1) for j in range(3): thisaxis[j][0]=ax[j] theaxis=ROOT.TMatrixD(cov,ROOT.TMatrixD.kMult,thisaxis) axis.append(theaxis) print 'scaled axis:' axis[-1].Print() #eigenwert=ROOT.TVector3() #(axis[0][0][0]/cl.axis(0)[0],axis[1][0][0]/cl.axis(1)[0],axis[2][0][0]/cl.axis(2)[0]) hsize=5 hist=ROOT.TH3D('hist','hist',100,-hsize,hsize,100,-hsize,hsize,100,-hsize,hsize) #points=[] points=ROOT.TPolyMarker3D() for dig in digis: points.SetNextPoint(dig.X(),dig.Y(),dig.Z()) cluster=ROOT.TPolyMarker3D() cluster.SetNextPoint(clpos.X(),clpos.Y(),clpos.Z()) cluster.SetMarkerStyle(20) cluster.SetMarkerSize(2) cluster.SetMarkerColor(2) axisl=[] for i in range(3): axisl.append(ROOT.TPolyLine3D()) # axisl[-1].SetNextPoint(clpos.X()-axis[i][0][0]/2,clpos.Y()-axis[i][1][0]/2,clpos.Z()-axis[i][2][0]/2) # axisl[-1].SetNextPoint(clpos.X()+axis[i][0][0]/2,clpos.Y()+axis[i][1][0]/2,clpos.Z()+axis[i][2][0]/2) axisl[-1].SetNextPoint(clpos.X(),clpos.Y(),clpos.Z()) axisl[-1].SetNextPoint(clpos.X()+axis[i][0][0],clpos.Y()+axis[i][1][0],clpos.Z()+axis[i][2][0]) points.SetMarkerStyle(20) hist.Draw() hist.GetXaxis().SetTitle("X") hist.GetYaxis().SetTitle("Y") hist.GetZaxis().SetTitle("Z") points.Draw('same') cluster.Draw('same') count=1 sphere=ROOT.TGeoSphere(0,1) trafo=ROOT.TGeoScale() trafo.SetScale(eigenwert[0][0],eigenwert[1][1],eigenwert[2][2]) sphere.SetTransform(trafo) #sphere.Draw('same') #volume=ROOT.TGeoVolume() #volume.AddNode(sphere,0,trafo) #volume.Draw('same') for ax in axisl: count+=1 ax.SetLineColor(count) ax.Draw('same') u=raw_input('done?')