Python Example (Test Page) » History » Version 2
Version 1 (Thomas Bretz, 2018-09-14 15:42) → Version 2/6 (Thomas Bretz, 2018-09-14 15:43)
[[TOC]]
Here is an example how to do a simple database based analysis in Python. This examples illustrates how to fill theta-square plots and are the Python counterparts to the example given on the main page.
The following Python code was kindly provided by Julian Kemp MSc.
# Query
Both examples on this page read the following query from an ascii file `query.txt`:
~~~ sql {{{#!sql
SELECT
Images.*,
Position.X,
Position.Y
FROM
Images
LEFT JOIN Position USING (FileId, EvtNumber)
LEFT JOIN RunInfo USING (FileId)
WHERE
fSourceKey=5
AND
fRunTypeKey=1
AND
FileId BETWEEN 180100000 AND 180200000
AND
fZenithDistanceMax<35
AND
fR750Cor>0.9*fR750Ref
AND
NumUsedPixels>5.5
AND
NumIslands<3.5
AND
Leakage1<0.1
AND
Width*Length*PI() < LOG10(Size)*898-1535;
~~~ }}}
# One shot example
The first example runs the query and filling the on- and off-histograms in a single python program:
~~~ python {{{#!python
import numpy as np
import matplotlib.pyplot as plt
import mysql.connector
def FillHistogram(bins,n,value,weight=1):
idx = np.searchsorted(bins,value)
if not idx==0 and not idx>n.size:
n[idx-1]+=weight
mm2deg = 0.0117193246260285378
# Create bins for the histograms
bins = np.linspace(0,1,56)
non = np.zeros(55)
noff = np.zeros(55)
# Read the MySQL query from a text file
queryFile = open('query.txt')
query = queryFile.read()
queryFile.close()
# Open a connection to the MySQL database
conn = mysql.connector.connect(user='fact',password='[...]',host='ihp-pc45.ethz.ch',database='factdata',port='3306')
cursor = conn.cursor(dictionary=True)
cursor.execute(query)
# Loop over all events received from the database
for Event in cursor:
# Intialize all variables needed in the calculations
Length = Event['Length']
Width = Event['Width']
NumIslands = Event['NumIslands']
NumUsedPixels = Event['NumUsedPixels']
Leakage1 = Event['Leakage1']
Size = Event['Size']
X = Event['X']
Y = Event['Y']
CosDelta = Event['CosDelta']
SinDelta = Event['SinDelta']
M3Long = Event['M3Long']
SlopeLong = Event['SlopeLong']
MeanX = Event['MeanX']
MeanY = Event['MeanY']
# First calculate all cuts to speedup the analysis
area = np.pi*Width*Length
# The abberation correction does increase also Width and Length by 1.02
cutq = (NumIslands<3.5 and NumUsedPixels>5.5 and Leakage1<0.1)
if not cutq:
continue
cut0 = area<(np.log10(Size)*898-1535)
if not cut0:
continue
# Loop over all wobble positions in the camera
for angle in range(0,360,60):
# ----------- Source dependent parameter calculation ----------
cr = np.cos(np.deg2rad(angle))
sr = np.sin(np.deg2rad(angle))
px = cr*X-sr*Y
py = cr*Y+sr*X
dx = MeanX-px*1.02
dy = MeanY-py*1.02
norm = np.sqrt(dx*dx+dy*dy)
dist = norm*mm2deg
lx = np.min([np.max([(CosDelta*dy-SinDelta*dx)/norm,-1.]),1.])
ly = np.min([np.max([(CosDelta*dx+SinDelta*dy)/norm,-1.]),1.])
alpha = np.arcsin(lx)
sgn = np.sign(ly)
# ------------------------ Application ---------------------
m3l = M3Long*sgn*mm2deg
slope = SlopeLong*sgn/mm2deg
# -------------------------- Analysis ----------------------
xi = 1.39252+0.154247*slope+1.67972*(1-1/(1+4.86232*Leakage1))
sign1 = m3l+0.07
sign2 = (dist-0.5)*7.2-slope
disp = (-xi if (sign1<0 or sign2<0) else xi)*(1-Width/Length)
thetasq = disp*disp+dist*dist-2*disp*dist*np.sqrt(1-lx*lx)
if angle==0:
FillHistogram(bins,non,thetasq)
else:
FillHistogram(bins,noff,thetasq,weight=1./5.)
conn.close()
plt.hist(bins[:-1],bins=bins,histtype='step',weights=non)
plt.hist(bins[:-1],bins=bins,histtype='step',weights=noff)
plt.show()
~~~ }}}
# == Two step example
==
This example works similar to `rootifysql` and a ROOT macro. It first requests the data from a database an stores it into a file and then processed the file with a second program.
## === Requesting and storing the data
~~~ python ===
{{{#!python
import numpy as np
import pickle
import mysql.connector
# Read the MySQL query from a text file
queryFile = open('query.txt')
query = queryFile.read()
queryFile.close()
# Open a connection to the MySQL database
conn = mysql.connector.connect(user='fact',password='[...]',host='ihp-pc45.ethz.ch',database='factdata',port='3306')
cursor = conn.cursor()
cursor.execute(query)
with open('crab-data-only.p','wb') as outFile:
description = np.array(cursor.description).T[0]
pickle.dump(description,outFile,pickle.HIGHEST_PROTOCOL)
for Event in cursor:
pickle.dump(Event,outFile,pickle.HIGHEST_PROTOCOL)
conn.close()
}}}
=== Reading and processing the data ===
{{{#!python
import numpy as np
import pickle
import matplotlib.pyplot as plt
def FillHistogram(bins,n,value,weight=1):
idx = np.searchsorted(bins,value)
if not idx==0 and not idx>n.size:
n[idx-1]+=weight
inFile = open('crab-data-only.p','rb')
description = pickle.load(inFile)
# Find indices corresponding to the different variables
Length_idx = np.nonzero(description=='Length')[0][0]
Width_idx = np.nonzero(description=='Width')[0][0]
NumIslands_idx = np.nonzero(description=='NumIslands')[0][0]
NumUsedPixels_idx = np.nonzero(description=='NumUsedPixels')[0][0]
Leakage1_idx = np.nonzero(description=='Leakage1')[0][0]
Size_idx = np.nonzero(description=='Size')[0][0]
X_idx = np.nonzero(description=='X')[0][0]
Y_idx = np.nonzero(description=='Y')[0][0]
CosDelta_idx = np.nonzero(description=='CosDelta')[0][0]
SinDelta_idx = np.nonzero(description=='SinDelta')[0][0]
M3Long_idx = np.nonzero(description=='M3Long')[0][0]
SlopeLong_idx = np.nonzero(description=='SlopeLong')[0][0]
MeanX_idx = np.nonzero(description=='MeanX')[0][0]
MeanY_idx = np.nonzero(description=='MeanY')[0][0]
mm2deg = 0.0117193246260285378
# Create bins for the histogram
bins = np.linspace(0,1,56)
non = np.zeros(55)
noff = np.zeros(55)
while True:
try:
Event = pickle.load(inFile)
except:
break
# Intialize all variables needed in the calculations
Length = Event[Length_idx]
Width = Event[Width_idx]
NumIslands = Event[NumIslands_idx]
NumUsedPixels = Event[NumUsedPixels_idx]
Leakage1 = Event[Leakage1_idx]
Size = Event[Size_idx]
X = Event[X_idx]
Y = Event[Y_idx]
CosDelta = Event[CosDelta_idx]
SinDelta = Event[SinDelta_idx]
M3Long = Event[M3Long_idx]
SlopeLong = Event[SlopeLong_idx]
MeanX = Event[MeanX_idx]
MeanY = Event[MeanY_idx]
# First calculate all cuts to speedup the analysis
area = np.pi*Width*Length
# The abberation correction does increase also Width and Length by 1.02
cutq = (NumIslands<3.5 and NumUsedPixels>5.5 and Leakage1<0.1)
if not cutq:
continue
cut0 = area<(np.log10(Size)*898-1535)
if not cut0:
continue
#print(area,cutq,cut0)
# Loop over all wobble positions in the camera
for angle in range(0,360,60):
# ----------- Source dependent parameter calculation ----------
cr = np.cos(np.deg2rad(angle))
sr = np.sin(np.deg2rad(angle))
px = cr*X-sr*Y
py = cr*Y+sr*X
dx = MeanX-px*1.02
dy = MeanY-py*1.02
norm = np.sqrt(dx*dx+dy*dy)
dist = norm*mm2deg
lx = np.min([np.max([(CosDelta*dy-SinDelta*dx)/norm,-1.]),1.])
ly = np.min([np.max([(CosDelta*dx+SinDelta*dy)/norm,-1.]),1.])
alpha = np.arcsin(lx)
sgn = np.sign(ly)
# ------------------------ Application ---------------------
m3l = M3Long*sgn*mm2deg
slope = SlopeLong*sgn/mm2deg
# -------------------------- Analysis ----------------------
xi = 1.39252+0.154247*slope+1.67972*(1-1/(1+4.86232*Leakage1))
sign1 = m3l+0.07
sign2 = (dist-0.5)*7.2-slope
disp = (-xi if (sign1<0 or sign2<0) else xi)*(1-Width/Length)
thetasq = disp*disp+dist*dist-2*disp*dist*np.sqrt(1-lx*lx)
if angle==0:
FillHistogram(bins,non,thetasq)
else:
FillHistogram(bins,noff,thetasq,weight=1./5.)
plt.hist(bins[:-1],bins=bins,histtype='step',weights=non)
plt.hist(bins[:-1],bins=bins,histtype='step',weights=noff)
plt.show()
~~~ }}}
Here is an example how to do a simple database based analysis in Python. This examples illustrates how to fill theta-square plots and are the Python counterparts to the example given on the main page.
The following Python code was kindly provided by Julian Kemp MSc.
# Query
Both examples on this page read the following query from an ascii file `query.txt`:
~~~ sql {{{#!sql
SELECT
Images.*,
Position.X,
Position.Y
FROM
Images
LEFT JOIN Position USING (FileId, EvtNumber)
LEFT JOIN RunInfo USING (FileId)
WHERE
fSourceKey=5
AND
fRunTypeKey=1
AND
FileId BETWEEN 180100000 AND 180200000
AND
fZenithDistanceMax<35
AND
fR750Cor>0.9*fR750Ref
AND
NumUsedPixels>5.5
AND
NumIslands<3.5
AND
Leakage1<0.1
AND
Width*Length*PI() < LOG10(Size)*898-1535;
~~~ }}}
# One shot example
The first example runs the query and filling the on- and off-histograms in a single python program:
~~~ python {{{#!python
import numpy as np
import matplotlib.pyplot as plt
import mysql.connector
def FillHistogram(bins,n,value,weight=1):
idx = np.searchsorted(bins,value)
if not idx==0 and not idx>n.size:
n[idx-1]+=weight
mm2deg = 0.0117193246260285378
# Create bins for the histograms
bins = np.linspace(0,1,56)
non = np.zeros(55)
noff = np.zeros(55)
# Read the MySQL query from a text file
queryFile = open('query.txt')
query = queryFile.read()
queryFile.close()
# Open a connection to the MySQL database
conn = mysql.connector.connect(user='fact',password='[...]',host='ihp-pc45.ethz.ch',database='factdata',port='3306')
cursor = conn.cursor(dictionary=True)
cursor.execute(query)
# Loop over all events received from the database
for Event in cursor:
# Intialize all variables needed in the calculations
Length = Event['Length']
Width = Event['Width']
NumIslands = Event['NumIslands']
NumUsedPixels = Event['NumUsedPixels']
Leakage1 = Event['Leakage1']
Size = Event['Size']
X = Event['X']
Y = Event['Y']
CosDelta = Event['CosDelta']
SinDelta = Event['SinDelta']
M3Long = Event['M3Long']
SlopeLong = Event['SlopeLong']
MeanX = Event['MeanX']
MeanY = Event['MeanY']
# First calculate all cuts to speedup the analysis
area = np.pi*Width*Length
# The abberation correction does increase also Width and Length by 1.02
cutq = (NumIslands<3.5 and NumUsedPixels>5.5 and Leakage1<0.1)
if not cutq:
continue
cut0 = area<(np.log10(Size)*898-1535)
if not cut0:
continue
# Loop over all wobble positions in the camera
for angle in range(0,360,60):
# ----------- Source dependent parameter calculation ----------
cr = np.cos(np.deg2rad(angle))
sr = np.sin(np.deg2rad(angle))
px = cr*X-sr*Y
py = cr*Y+sr*X
dx = MeanX-px*1.02
dy = MeanY-py*1.02
norm = np.sqrt(dx*dx+dy*dy)
dist = norm*mm2deg
lx = np.min([np.max([(CosDelta*dy-SinDelta*dx)/norm,-1.]),1.])
ly = np.min([np.max([(CosDelta*dx+SinDelta*dy)/norm,-1.]),1.])
alpha = np.arcsin(lx)
sgn = np.sign(ly)
# ------------------------ Application ---------------------
m3l = M3Long*sgn*mm2deg
slope = SlopeLong*sgn/mm2deg
# -------------------------- Analysis ----------------------
xi = 1.39252+0.154247*slope+1.67972*(1-1/(1+4.86232*Leakage1))
sign1 = m3l+0.07
sign2 = (dist-0.5)*7.2-slope
disp = (-xi if (sign1<0 or sign2<0) else xi)*(1-Width/Length)
thetasq = disp*disp+dist*dist-2*disp*dist*np.sqrt(1-lx*lx)
if angle==0:
FillHistogram(bins,non,thetasq)
else:
FillHistogram(bins,noff,thetasq,weight=1./5.)
conn.close()
plt.hist(bins[:-1],bins=bins,histtype='step',weights=non)
plt.hist(bins[:-1],bins=bins,histtype='step',weights=noff)
plt.show()
~~~ }}}
# == Two step example
==
This example works similar to `rootifysql` and a ROOT macro. It first requests the data from a database an stores it into a file and then processed the file with a second program.
## === Requesting and storing the data
~~~ python ===
{{{#!python
import numpy as np
import pickle
import mysql.connector
# Read the MySQL query from a text file
queryFile = open('query.txt')
query = queryFile.read()
queryFile.close()
# Open a connection to the MySQL database
conn = mysql.connector.connect(user='fact',password='[...]',host='ihp-pc45.ethz.ch',database='factdata',port='3306')
cursor = conn.cursor()
cursor.execute(query)
with open('crab-data-only.p','wb') as outFile:
description = np.array(cursor.description).T[0]
pickle.dump(description,outFile,pickle.HIGHEST_PROTOCOL)
for Event in cursor:
pickle.dump(Event,outFile,pickle.HIGHEST_PROTOCOL)
conn.close()
}}}
=== Reading and processing the data ===
{{{#!python
import numpy as np
import pickle
import matplotlib.pyplot as plt
def FillHistogram(bins,n,value,weight=1):
idx = np.searchsorted(bins,value)
if not idx==0 and not idx>n.size:
n[idx-1]+=weight
inFile = open('crab-data-only.p','rb')
description = pickle.load(inFile)
# Find indices corresponding to the different variables
Length_idx = np.nonzero(description=='Length')[0][0]
Width_idx = np.nonzero(description=='Width')[0][0]
NumIslands_idx = np.nonzero(description=='NumIslands')[0][0]
NumUsedPixels_idx = np.nonzero(description=='NumUsedPixels')[0][0]
Leakage1_idx = np.nonzero(description=='Leakage1')[0][0]
Size_idx = np.nonzero(description=='Size')[0][0]
X_idx = np.nonzero(description=='X')[0][0]
Y_idx = np.nonzero(description=='Y')[0][0]
CosDelta_idx = np.nonzero(description=='CosDelta')[0][0]
SinDelta_idx = np.nonzero(description=='SinDelta')[0][0]
M3Long_idx = np.nonzero(description=='M3Long')[0][0]
SlopeLong_idx = np.nonzero(description=='SlopeLong')[0][0]
MeanX_idx = np.nonzero(description=='MeanX')[0][0]
MeanY_idx = np.nonzero(description=='MeanY')[0][0]
mm2deg = 0.0117193246260285378
# Create bins for the histogram
bins = np.linspace(0,1,56)
non = np.zeros(55)
noff = np.zeros(55)
while True:
try:
Event = pickle.load(inFile)
except:
break
# Intialize all variables needed in the calculations
Length = Event[Length_idx]
Width = Event[Width_idx]
NumIslands = Event[NumIslands_idx]
NumUsedPixels = Event[NumUsedPixels_idx]
Leakage1 = Event[Leakage1_idx]
Size = Event[Size_idx]
X = Event[X_idx]
Y = Event[Y_idx]
CosDelta = Event[CosDelta_idx]
SinDelta = Event[SinDelta_idx]
M3Long = Event[M3Long_idx]
SlopeLong = Event[SlopeLong_idx]
MeanX = Event[MeanX_idx]
MeanY = Event[MeanY_idx]
# First calculate all cuts to speedup the analysis
area = np.pi*Width*Length
# The abberation correction does increase also Width and Length by 1.02
cutq = (NumIslands<3.5 and NumUsedPixels>5.5 and Leakage1<0.1)
if not cutq:
continue
cut0 = area<(np.log10(Size)*898-1535)
if not cut0:
continue
#print(area,cutq,cut0)
# Loop over all wobble positions in the camera
for angle in range(0,360,60):
# ----------- Source dependent parameter calculation ----------
cr = np.cos(np.deg2rad(angle))
sr = np.sin(np.deg2rad(angle))
px = cr*X-sr*Y
py = cr*Y+sr*X
dx = MeanX-px*1.02
dy = MeanY-py*1.02
norm = np.sqrt(dx*dx+dy*dy)
dist = norm*mm2deg
lx = np.min([np.max([(CosDelta*dy-SinDelta*dx)/norm,-1.]),1.])
ly = np.min([np.max([(CosDelta*dx+SinDelta*dy)/norm,-1.]),1.])
alpha = np.arcsin(lx)
sgn = np.sign(ly)
# ------------------------ Application ---------------------
m3l = M3Long*sgn*mm2deg
slope = SlopeLong*sgn/mm2deg
# -------------------------- Analysis ----------------------
xi = 1.39252+0.154247*slope+1.67972*(1-1/(1+4.86232*Leakage1))
sign1 = m3l+0.07
sign2 = (dist-0.5)*7.2-slope
disp = (-xi if (sign1<0 or sign2<0) else xi)*(1-Width/Length)
thetasq = disp*disp+dist*dist-2*disp*dist*np.sqrt(1-lx*lx)
if angle==0:
FillHistogram(bins,non,thetasq)
else:
FillHistogram(bins,noff,thetasq,weight=1./5.)
plt.hist(bins[:-1],bins=bins,histtype='step',weights=non)
plt.hist(bins[:-1],bins=bins,histtype='step',weights=noff)
plt.show()
~~~ }}}