In [1]:
%pylab inline
import numpy as np
import sys
import os
sys.path.append('/nfsopt/root/v6-06-06/build/lib')
import ROOT
from ipywidgets import interact, interactive
from IPython.display import Javascript


Populating the interactive namespace from numpy and matplotlib
Welcome to ROOTaaS 6.06/06

In [2]:
f  = ROOT.TFile('~/alice1/ali-master/AliRoot/v5/FIT.Hits.root')
f2 = ROOT.TFile('~/alice1/ali-master/AliRoot/v5/Kinematics.root')
f3 = ROOT.TFile('~/alice1/ali-master/AliRoot/v5/galice.root')
f4 = ROOT.TFile('~/alice1/ali-master/AliRoot/v4/someting.root')
f5 = ROOT.TFile('~/alice1/ali-master/AliRoot/v5/someting.root')
f6 = ROOT.TFile('~/alice1/ali-master/AliRoot/v6/someting.root')
f7 = ROOT.TFile('~/alice1/ali-master/AliRoot/v7/someting.root')
f8 = ROOT.TFile('~/alice1/ali-master/AliRoot/v8/someting.root')


TClass::Init:0: RuntimeWarning: no dictionary for class AliHit is available
TClass::Init:0: RuntimeWarning: no dictionary for class AliFITHits is available
TClass::Init:0: RuntimeWarning: no dictionary for class AliHeader is available
TClass::Init:0: RuntimeWarning: no dictionary for class AliStack is available
TClass::Init:0: RuntimeWarning: no dictionary for class AliGenEventHeader is available
TClass::Init:0: RuntimeWarning: no dictionary for class AliRun is available
TClass::Init:0: RuntimeWarning: no dictionary for class AliBODY is available
TClass::Init:0: RuntimeWarning: no dictionary for class AliModule is available
TClass::Init:0: RuntimeWarning: no dictionary for class AliPIPEupgrade is available
TClass::Init:0: RuntimeWarning: no dictionary for class AliPIPE is available
TClass::Init:0: RuntimeWarning: no dictionary for class AliMFT is available
TClass::Init:0: RuntimeWarning: no dictionary for class AliDetector is available
TClass::Init:0: RuntimeWarning: no dictionary for class AliMFTSegmentation is available
TClass::Init:0: RuntimeWarning: no dictionary for class AliMFTGeomTGeo is available
TClass::Init:0: RuntimeWarning: no dictionary for class AliITSMFTGeomTGeo is available
TClass::Init:0: RuntimeWarning: no dictionary for class AliFITv5 is available
TClass::Init:0: RuntimeWarning: no dictionary for class AliFIT is available
TClass::Init:0: RuntimeWarning: no dictionary for class AliMC is available
TClass::Init:0: RuntimeWarning: no dictionary for class AliGenDPMjet is available
TClass::Init:0: RuntimeWarning: no dictionary for class AliGenMC is available
TClass::Init:0: RuntimeWarning: no dictionary for class AliGenerator is available
TClass::Init:0: RuntimeWarning: no dictionary for class AliRndm is available
TClass::Init:0: RuntimeWarning: no dictionary for class TDPMjet is available
TClass::Init:0: RuntimeWarning: no dictionary for class AliGenDPMjetEventHeader is available
TClass::Init:0: RuntimeWarning: no dictionary for class AliCollisionGeometry is available
TClass::Init:0: RuntimeWarning: no dictionary for class AliRunLoader is available
TClass::Init:0: RuntimeWarning: no dictionary for class AliLoader is available
TClass::Init:0: RuntimeWarning: no dictionary for class AliDataLoader is available
TClass::Init:0: RuntimeWarning: no dictionary for class AliTreeLoader is available
TClass::Init:0: RuntimeWarning: no dictionary for class AliObjectLoader is available
TClass::Init:0: RuntimeWarning: no dictionary for class AliBaseLoader is available

In [3]:
events = 100

hdir = list(np.zeros(events))
kdir = list(np.zeros(events))


x = 0
for i in range(events):
    eventname = "Event%d;1" % x
    
    htemp = f.GetDirectory(eventname)
    ktemp = f2.GetDirectory(eventname)
    hdir[i] = htemp
    kdir[i] = ktemp
    x += 1

In [4]:
hTree = list(np.zeros(events))
kTree = list(np.zeros(events))
for i in range(events):
    htemp = hdir[i]
    ktemp = kdir[i]
    hTree[i] = htemp.Get("TreeH;1")
    kTree[i] = ktemp.Get("TreeK;1")

In [5]:
c1= ROOT.TCanvas("c1","c1",800,600);
h1_x = ROOT.TH1D("h1_x", "fMCP", 100, -1, 60)
for i in range(events):
    hTree[i].Draw("fMCP>>h1_x")
c1.Draw()



In [6]:
def conc(r, s): #r = radius of array, s = width of single detector
    dP = s/2.
    theta = 2.*np.arctan(dP/r)
    
    #Translations 
    n = [-3.,-2.,-1.,1.,2.,3.]
    gridpoints = []
    for i in n:
        gridpoints.append(r*np.sin((1-1/(2*abs(i)))*i*theta))
        
    xi = np.array([gridpoints[1],gridpoints[2],gridpoints[3],gridpoints[4],
       gridpoints[0],gridpoints[1],gridpoints[2],gridpoints[3],gridpoints[4],gridpoints[5],
       gridpoints[0],gridpoints[1],gridpoints[4],gridpoints[5],
       gridpoints[0],gridpoints[1],gridpoints[4],gridpoints[5],
       gridpoints[0],gridpoints[1],gridpoints[2],gridpoints[3],gridpoints[4],gridpoints[5],
       gridpoints[1],gridpoints[2],gridpoints[3],gridpoints[4]])
    yi = np.array([gridpoints[5],gridpoints[5],gridpoints[5],gridpoints[5],
       gridpoints[4],gridpoints[4],gridpoints[4],gridpoints[4],gridpoints[4],gridpoints[4],
       gridpoints[3],gridpoints[3],gridpoints[3],gridpoints[3],
       gridpoints[2],gridpoints[2],gridpoints[2],gridpoints[2],
       gridpoints[1],gridpoints[1],gridpoints[1],gridpoints[1],gridpoints[1],gridpoints[1],
       gridpoints[0],gridpoints[0],gridpoints[0],gridpoints[0]])
    
    zi = np.array([np.sqrt(r**2 - x**2 - y**2) for x,y in zip(xi,yi)])
    
    #Rotations
    alpha = np.array([(np.arctan(y/x) - np.pi/2 + 2*np.pi) for x,y in zip(xi,yi)])
    beta = np.array([np.arccos(z/r) if x < 0 else -np.arccos(z/r) for x,z in zip(xi,zi)])
    
    #Translations Revisited
    rcomp = r + 2.5
    xf = rcomp*np.cos(alpha + np.pi/2)*np.sin(-beta)
    yf = rcomp*np.sin(alpha + np.pi/2)*np.sin(-beta)
    zf = rcomp*np.cos(beta)
    
    #Conversions
    alpha = alpha*180/np.pi
    beta = beta*180/np.pi
    gamma = -1*alpha
    
    return alpha, beta, gamma, xf, yf, zf

In [7]:
v4hdist2D = f4.Get("v4hdist2D")
v5hdist2D = f5.Get("v5hdist2D")
v6hdist2D = f6.Get("v6hdist2D")
v7hdist2D = f7.Get("v7hdist2D")
v8hdist2D = f8.Get("v8hdist2D")

hdist2Dlist = [v4hdist2D, v5hdist2D, v6hdist2D, v7hdist2D, v8hdist2D]

In [8]:
"v%d" % 3


Out[8]:
'v3'

In [9]:
Rx = lambda ang: np.array([[1,0,0],[0,np.cos(ang), -np.sin(ang)],[0,np.sin(ang),np.cos(ang)]])
Ry = lambda ang: np.array([[np.cos(ang),0,np.sin(ang)],[0,1,0],[-np.sin(ang),0,np.cos(ang)]])
Rz = lambda ang: np.array([[np.cos(ang), -np.sin(ang),0],[np.sin(ang),np.cos(ang),0],[0,0,1]])

boxes = True
def boxmaker2D( index, x, y, ea, title="XY of Hits"):
    c2= ROOT.TCanvas("c2","c2",1200,1200)
    hs = ROOT.THStack("hs", title)
    for i in index:
        hs.Add(hdist2Dlist[i])
    hs.Draw("colz nostack")
    hs.GetXaxis().SetTitle("X (cm)")
    hs.GetYaxis().SetTitle("Y (cm)")
    boxlistq = list(np.zeros(len(x)))
    boxlistm = list(np.zeros(len(x)))
    wq = 5.3/2.
    wm = 5.9/2.
    if boxes:
        for i, n in enumerate(x):
            aq = np.array([ wq, wq, -2.2])  #-2.2 for the z should place plane of hits on the surface of the mcp,  
            bq = np.array([ wq,-wq, -2.2])  #NOT the surface of the quartz radiator. The dimensions of the detector 
            cq = np.array([-wq,-wq, -2.2])  #modules as defined in the alifitv#.cxx looks to be 5.9 x 5.9 x 5 cm, 
            dq = np.array([-wq, wq, -2.2])  #with the quartz element being 5.296 x 5.296 x 2.2 cm and the mcp being
                                            #5.898 x 5.898 x 2.8 cm. 
            
            am = np.array([ wm, wm, .3])  
            bm = np.array([ wm,-wm, .3]) 
            cm = np.array([-wm,-wm, .3])
            dm = np.array([-wm, wm, .3])
            
            Z1 = Rz(ea[0][i])
            X2 = Rx(ea[1][i])
            Z2 = Rz(ea[2][i])
            
            #gamma
            aq = np.matmul(Z2,aq)
            bq = np.matmul(Z2,bq)
            cq = np.matmul(Z2,cq)
            dq = np.matmul(Z2,dq)
            am = np.matmul(Z2,am)
            bm = np.matmul(Z2,bm)
            cm = np.matmul(Z2,cm)
            dm = np.matmul(Z2,dm)



            #beta
            aq = np.matmul(X2,aq)
            bq = np.matmul(X2,bq)
            cq = np.matmul(X2,cq)
            dq = np.matmul(X2,dq)
            am = np.matmul(X2,am)
            bm = np.matmul(X2,bm)
            cm = np.matmul(X2,cm)
            dm = np.matmul(X2,dm)



            #alpha
            aq = np.matmul(Z1,aq)
            bq = np.matmul(Z1,bq)
            cq = np.matmul(Z1,cq)
            dq = np.matmul(Z1,dq)
            am = np.matmul(Z1,am)
            bm = np.matmul(Z1,bm)
            cm = np.matmul(Z1,cm)
            dm = np.matmul(Z1,dm)
            

            boxlistq[i] = ROOT.TPolyLine(4)
            boxlistq[i].SetPoint(0, aq[0]+x[i], aq[1]+y[i])
            boxlistq[i].SetPoint(1, bq[0]+x[i], bq[1]+y[i])
            boxlistq[i].SetPoint(2, cq[0]+x[i], cq[1]+y[i])
            boxlistq[i].SetPoint(3, dq[0]+x[i], dq[1]+y[i])
            boxlistq[i].SetPoint(4, aq[0]+x[i], aq[1]+y[i])
            boxlistq[i].Draw()
            boxlistm[i] = ROOT.TPolyLine(4)
            boxlistm[i].SetPoint(0, am[0]+x[i], am[1]+y[i])
            boxlistm[i].SetPoint(1, bm[0]+x[i], bm[1]+y[i])
            boxlistm[i].SetPoint(2, cm[0]+x[i], cm[1]+y[i])
            boxlistm[i].SetPoint(3, dm[0]+x[i], dm[1]+y[i])
            boxlistm[i].SetPoint(4, am[0]+x[i], am[1]+y[i])
            boxlistm[i].Draw("same")
            
    c2.Draw()
r = 82
s = 3.31735114408*2
a, b, g, x, y, z = conc(r, s)

ea = [np.pi/180*a, np.pi/180*b, np.pi/180*g]
boxmaker2D([4],x, y, ea, title="XY of Hits v8")



In [10]:
xc = np.array([-8.95, -3.05, 3.05, 8.95, 
                  -14.85, -8.95, -3.05, 3.05,
                  8.95,  14.85,-14.85, -8.95, 
                  8.95, 14.85, -14.85, -8.95, 
                  8.95, 14.85, -14.85,-8.95,
                  -3.05, 3.05, 8.95, 14.85, 
                  -8.95, -3.05, 3.05, 8.95])
yc = np.array([14.85, 14.85,14.85,14.85,
                  8.95,8.95,8.95,8.95,8.95,8.95,
                  3.05, 3.05,3.05,3.05,
                  -3.05,-3.05,-3.05,-3.05,
                  -8.95,-8.95,-8.95,-8.95,-8.95,-8.95,
                  -14.85,-14.85,-14.85,-14.85])
zc = np.ones(28)*r
r = 80
s = 5.9 
a, b, g, x, y, z = conc(r, s)

boxmaker2D([0],xc, yc, [np.zeros(28),np.zeros(28),np.zeros(28)], title="XY of Hits v4")



In [11]:
def ThreeDModel(RunID = 0, Theta=45,Phi=45,Zoom = 1):
    #Draw the canvas with the chamber view
    c3= ROOT.TCanvas("c3","c3",800,600);
    h3 = ROOT.TH3F("h3","XYZ of Hits", 10, -90*Zoom, -80*Zoom, 10, -25*Zoom, 25*Zoom, 10, -25*Zoom, 25*Zoom);
    h3.SetStats(0)
    h3.SetXTitle("  Z(cm)")
    h3.SetYTitle("  X(cm)") #Note, Y and Z are flipped to get the beam perspective right
    h3.SetZTitle("  Y(cm)")
    c1.SetPhi(Phi)
    c1.SetTheta(Theta)
    for i in range(events):
        hTree[i].Draw("fY:fX:fZ>>h3")
    h3.Draw()
    c3.Draw()

In [12]:
interact(ThreeDModel,RunID=(0,3,1),Theta=(0,90,1),Phi=(0,90,1),Zoom = (0,2,0.1))


Widget Javascript not detected.  It may not be installed or enabled properly.
Out[12]:
<function __main__.ThreeDModel>
[IPKernelApp] WARNING | Widget Javascript not detected.  It may not be installed or enabled properly.

In [13]:
"""
for i in range(10):
    FIT1 = hTree[i].GetBranch('b_FIT_')
    for entry in hTree[i]:
        x = entry.FIT
        for j in x:
            
            #print(j)
            """


Out[13]:
"\nfor i in range(10):\n    FIT1 = hTree[i].GetBranch('b_FIT_')\n    for entry in hTree[i]:\n        x = entry.FIT\n        for j in x:\n            \n            #print(j)\n            "

In [14]:
c4 = ROOT.TCanvas()
#hTree1.Draw("fX:fY")
for i in range(100):
    kTree[i].Draw("fPdgCode", "")
c4.Draw()
#c1.Print("XY.png")



In [15]:
print('FIT.Hits.root HTree1 has ' + str(hTree[0].GetEntries()) + ' entries')
print('FIT.Hits.root HTree2 has ' + str(hTree[1].GetEntries()) + ' entries')
print('Kinematics.root KTree1 has ' + str(kTree[0].GetEntries()) + ' entries')
print('Kinematics.root KTree2 has ' + str(kTree[1].GetEntries()) + ' entries')
print('\n???')


FIT.Hits.root HTree1 has 73 entries
FIT.Hits.root HTree2 has 66 entries
Kinematics.root KTree1 has 175 entries
Kinematics.root KTree2 has 217 entries

???

In [16]:
particles1 = kTree[0].GetBranch("Particles")
particles2 = kTree[1].GetBranch("Particles")

particles2.GetNleaves()


Out[16]:
1

In [17]:
f5.Print()
f6.Print()
f7.Print()


TFile: name=/nfshome/nmille16/alice1/ali-master/AliRoot/v5/someting.root, title=, option=READ
TH1.Print Name  = v5hdist2D, Entries= 50698648, Total sum= 5.06986e+07
TFile: name=/nfshome/nmille16/alice1/ali-master/AliRoot/v6/someting.root, title=, option=READ
TH1.Print Name  = v6hdist2D, Entries= 44237713, Total sum= 4.42377e+07
TFile: name=/nfshome/nmille16/alice1/ali-master/AliRoot/v7/someting.root, title=, option=READ
TH1.Print Name  = v7hdist2D, Entries= 242055724, Total sum= 2.42056e+08

In [18]:
f5.ls()


TFile**		/nfshome/nmille16/alice1/ali-master/AliRoot/v5/someting.root	
 TFile*		/nfshome/nmille16/alice1/ali-master/AliRoot/v5/someting.root	
  OBJ: TH2F	v5hdist2D	v5 : 0 at: 0xa4dbb50
  KEY: TH1F	v5hpdgcode;1	v5
  KEY: TH2F	v5htof;1	tofv5
  KEY: TH2F	v5hdist2D;1	v5
  KEY: TH1F	v5ehist;1	Eloss summed v5, E = 0 Excluded
  KEY: TH1F	v5ehist2;1	Etot summed
  KEY: TH1F	v5ehist3;1	Etot
  KEY: TH1F	v5ehist4;1	Eloss
  KEY: TH1F	v5hdist;1	v5
  KEY: TH1F	v5nPrim;1	v5nPrimaries per event
  KEY: TH2F	v5effichist;1	& Efficiency vs. number of Primaries for v5
  KEY: TH1F	v5nptrue;1	Number of true primaries for v5
  KEY: TH1F	v5nptrue2;1	Number of true primaries for v5
  KEY: TH1F	v5npmeas;1	Number of & measured primaries for v5
  KEY: TH1F	v5npmeas2;1	Number of || measured primaries for v5
  KEY: TH1F	v5effich;1	& Efficiency for v5
  KEY: TH1F	v5effich2;1	|| Efficiency for v5
  KEY: TH1F	v5TracksH;1	v5TracksH per event
  KEY: THStack	hs2;1	v5 nPrimaries vs. TracksH
  KEY: THStack	hs3;1	v5 Measured vs. Expected Primaries
  KEY: TH3F	v5hits3d;1	XYZ of hits for v5
  KEY: TH2F	v5PtEta;1	Hit Eta vs. Transverse Momentum for v5
  KEY: TH2F	v5pPtEta;1	Primary Hit Eta vs. Transverse Momentum for v5
  KEY: TH2F	v5kPtEta;1	Primary Eta vs. Transverse Momentum for v5

In [19]:
hpdgcode4 = f4.Get("v4hpdgcode")
hpdgcode5 = f5.Get("v5hpdgcode")
hpdgcode6 = f6.Get("v6hpdgcode")
hpdgcode7 = f7.Get("v7hpdgcode")
hpdgcode8 = f8.Get("v8hpdgcode")

In [20]:
c5 = ROOT.TCanvas()

hs = ROOT.THStack("hs", "PDG Codes")

hpdgcode4.SetLineColorAlpha(6, 1)
hs.Add(hpdgcode4)

hpdgcode5.SetLineColorAlpha(4, 1)
hs.Add(hpdgcode5)

hpdgcode6.SetLineColorAlpha(2, 1)
hs.Add(hpdgcode6)

hpdgcode7.SetLineColorAlpha(8, 1)
hs.Add(hpdgcode7)

hpdgcode8.SetLineColorAlpha(9, 1)
hs.Add(hpdgcode8)


hs.Draw()
hs.GetXaxis().SetRangeUser(-100,200)
hs.GetXaxis().SetTitle("pdgcode")
c5.BuildLegend(0.78,0.7,0.98, 0.95)
c5.SetPhi(0)
c5.SetTheta(0)
c5.Draw()



In [ ]:


In [21]:
htof4 = f4.Get("v4htof")
htof5 = f5.Get("v5htof")
htof6 = f6.Get("v6htof")
htof7 = f7.Get("v7htof")
htof8 = f8.Get("v8htof")

In [22]:
c6 = ROOT.TCanvas("c6","c6",1200,1200)
c6.Divide(2,2)
xmin = 0
xmax = 6500

htof4.GetXaxis().SetRangeUser(xmin, xmax)
htof4.SetLineColorAlpha(1, 1)
htof4.SetMarkerColorAlpha(1, 1)
htof4.SetFillColor(1)
htof4.SetFillStyle(3001)

htof5.GetXaxis().SetRangeUser(xmin, xmax)
htof5.SetLineColorAlpha(4, 1)
htof5.SetMarkerColorAlpha(4, 1)
htof5.SetFillColor(4)
htof5.SetFillStyle(3001)

htof6.GetXaxis().SetRangeUser(xmin, xmax)
htof6.SetLineColorAlpha(2, 1)
htof6.SetMarkerColorAlpha(2, 1)
htof6.SetFillColor(2)
htof6.SetFillStyle(3001)

htof7.GetXaxis().SetRangeUser(xmin, xmax)
htof7.SetLineColorAlpha(8, 1)
htof7.SetMarkerColorAlpha(8, 1)
htof7.SetFillColor(8)
htof7.SetFillStyle(3001)


c6.cd(1)
htof4.Draw()
c6.cd(2)
htof5.Draw()
c6.cd(3)
htof6.Draw()
c6.cd(4)
htof7.Draw()


c6.cd(0)
legend = ROOT.TLegend(0.40,0.459,0.60, 0.54)
#legend.SetHeader("The Legend Title")
legend.AddEntry(htof4,"v4, mean: %7.2f" % (htof4.GetMean(2)),"f")
legend.AddEntry(htof5,"v5, mean: %7.2f" % (htof5.GetMean(2)),"f")
legend.AddEntry(htof6,"v6, mean: %7.2f" % (htof6.GetMean(2)),"f")
legend.AddEntry(htof7,"v7, mean: %7.2f" % (htof7.GetMean(2)),"f")
legend.Draw()

c6.Draw()