helmart変換を行ったうえでパターンマッチを行うpythonコードです。importにあるgrid_affineとaffineはJyoshidaのgistから入手してください。gorgデータは54番サーバー内にあるのでそこから入手してください。

helmart.py

# -*- coding: utf-8 -*-
"""
Created on Wed Aug 17 15:27:51 2016

@author: ryousuke
"""

import affine
import numpy as np
#from array import array
#import re
import ROOT
import pandas as pd
from os.path import join, relpath
from glob import glob
import itertools

import grid_affine

import gc



ROOT.std.__file__ = 'dummy_for_old_pyastro'

ROOT.gROOT.Reset()
c1 = ROOT.TCanvas( 'c1', 'Canvas', 700, 700 )

factorx = 0.267
factory = 0.267

gridfile = ['c_coordinate','e_coordinate','n_coordinate','s_coordinate','w_coordinate']
s_l_d = ['c_mt_d','e_mt_d','n_mt_d','s_mt_d','w_mt_d']
s_l_u = ['c_mt_u','e_mt_u','n_mt_u','s_mt_u','w_mt_u']

mod_n = '004'
d_p_n = '01'
u_p_n = '02'

path_d = 'C:\\Users\\ryousuke\\Desktop\\study_log\\code\\python\\JYoshida\\patternmatch\\mod{0}pl{1}tracks\\'.format(mod_n,d_p_n)
files_d = [relpath(x, path_d) for x in glob(join(path_d, '*'))]

path_u = 'C:\\Users\\ryousuke\\Desktop\\study_log\\code\\python\\JYoshida\\patternmatch\\mod{0}pl{1}tracks\\'.format(mod_n,u_p_n)
files_u = [relpath(x, path_u) for x in glob(join(path_u, '*'))]

gorg = grid_affine.GridMark('.\\grid.gorg')


for q,l in itertools.product(range(len(s_l_d)), range(len(files_d))):
    
#    print '{0}_{1}'.format(i,l)
    
    if files_d[l].find(s_l_d[q]) > -1:
        
        filename_d = path_d + files_d[l]
        
#        print filename_d
        
        
        for m in range(len(files_u)):
            
            if files_u[m].find(s_l_u[q]) > -1:
                filename_u = path_u + files_u[m]
                
#                print filename_u
                
                
                filename_coor_d = path_d + gridfile[q] + '.txt'
                af_d = affine.Affine()                
                for text in open(filename_coor_d).readlines():
                    column = np.array(text.split())
        
                    x_p = float(column[0])
                    y_p = float(column[1])
                    z_p = float(column[2])
                    
                    nega_p = gorg.getTheNearest(x_p,y_p)                    
                    
                    print '{0}_{1} {2}_{3}'.format(x_p,y_p,nega_p[0],nega_p[1])
                    
                    af_d.add(x_p,y_p, nega_p[0],nega_p[1])
                    
                prams = af_d.calc()
                af_d.printParams()
                
                
                filename_coor_u = path_u + gridfile[q] + '.txt'
                af_u = affine.Affine()                
                for text in open(filename_coor_u).readlines():
                    column = np.array(text.split())
        
                    x_p = float(column[0])
                    y_p = float(column[1])
                    z_p = float(column[2])
                    
                    nega_p = gorg.getTheNearest(x_p,y_p)
                    
                    print '{0}_{1} {2}_{3}'.format(x_p,y_p,nega_p[0],nega_p[1])
                    
                    af_u.add(x_p,y_p, nega_p[0],nega_p[1])
                    
                prams = af_u.calc()
                af_u.printParams()
              
                ROOT.gROOT.Reset()

                col_names = ['ph','pv','iax','iay','px','py','vx','vy']
                    
                mtd = pd.read_csv(filename_d, header=None, sep=r"\s+", names=col_names)
                mtu = pd.read_csv(filename_u, header=None, sep=r"\s+", names=col_names)                
    
                tnd = ROOT.TNtuple('tnd','track','ph:pv:iax:iay:px:py:vx:vy:tx:ty')
                for i, t in mtd.iterrows():
                    
                    tx = t.vx*1000.0 + t.px*factorx
                    ty = t.vy*1000.0 - t.py*factory
                    
                    atra = af_d.affTrans(tx,ty)
                    
                    tnd.Fill(t.ph, t.pv, t.iax, t.iay, t.px, t.py, t.vx, t.vy, atra[0], atra[1])
                
                tnu = ROOT.TNtuple('tnu','track','ph:pv:iax:iay:px:py:vx:vy:tx:ty')
                for i, t in mtu.iterrows():
                    tx = t.vx*1000.0 + t.px*factorx
                    ty = t.vy*1000.0 - t.py*factory
                
                    atra = af_u.affTrans(tx,ty)    
                    
                    tnu.Fill(t.ph, t.pv, t.iax, t.iay, t.px, t.py, t.vx, t.vy, atra[0],atra[1])
                
                
                h2dydx = ROOT.TH2D("h2dydx","up;x[mm]",360, -1800, 1800, 360, -1800, 1800)
                 
                    #combination
                for di, dt in mtd.iterrows():
                            
                    if dt.pv < 30:
                        continue
                    
                    tx_d = dt.vx*1000.0 + dt.px*factorx
                    ty_d = dt.vy*1000.0 - dt.py*factory
                    
                    atra_d = af_d.affTrans(tx_d,ty_d)
                    
                    dox = atra_d[0]
                    doy = atra_d[1]
                        
                        
                    utarr = mtu[ np.logical_and(abs(mtu.vx - dt.vx )<0.5, abs(mtu.vy - dt.vy)<0.5)]
                
                    for ui, ut in utarr.iterrows():
                        
                        if ut.pv < 30:
                            continue
                        
                        tx_u = ut.vx*1000 + ut.px*factorx
                        ty_u = ut.vy*1000 - ut.py*factory
                        
                        atra_u = af_u.affTrans(tx_u,ty_u)
                        
                        uox = atra_u[0]
                        uoy = atra_u[1]
                        
                        dx = uox - dox
                        dy = uoy - doy
                        
                        h2dydx.Fill(dy,dx)
                        
                    if di%500==0:
                        print di, " / ", len(mtd) 
                      
                    
                print 'finish'
                h2dydx.Draw('colz')
                c1.Print('h2dydx_mod{0}pl{1}-{2}_{3}.png'.format(mod_n, d_p_n, u_p_n, gridfile[q]))
                gc.collect()