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
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
fZenithDistanceMax0.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
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
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()
}}}