#!/usr/bin/env python

##########################################################################################################################
#
# 2016/12/20 Paul Mattione
#
# Python script for transforming Igor's inconsistent data formats to a standard one
#
# Output format:
# E_gamma (GeV)    cos_theta_meson_CM    dsigma/dOmega (ub/sr)    dsigma/dOmega Stat Error (ub/sr)    dsigma/dOmega Syst Error (ub/sr)    
# dsigma/dOmega Total Angle-Dependent Error (ub/sr)    dsigma/dOmega Energy-Dependent Scale Error (ub/sr)      
# dsigma/dOmega Energy-Independent Scale Error (ub/sr)    dsigma/dOmega Total Scale Error (ub/sr)    dsigma/dOmega Total Error (ub/sr) 
#
##########################################################################################################################

import sys
import os
from collections import defaultdict
import math

#################################################### GLOBAL VARIABLES ####################################################

VERBOSE = 0
MAX_COSTHETADIFF = -1.05 # Negative to accept all (or really big!)

###################################################### PRINT HEADER ######################################################

def Print_Header(locOutputFile):
	locOutputFile.write('E_gamma (GeV)   cos_theta_meson_CM   dsigma/dOmega (ub/sr)   ')
	locOutputFile.write('dsigma/dOmega Stat Error (ub/sr)   dsigma/dOmega Syst Error (ub/sr)   dsigma/dOmega Total Angle-Dependent Error (ub/sr)   ')
	locOutputFile.write('dsigma/dOmega Energy-Dependent Scale Error (ub/sr)   dsigma/dOmega Energy-Independent Scale Error (ub/sr)   ')
	locOutputFile.write('dsigma/dOmega Total Scale Error (ub/sr)   dsigma/dOmega Total Error (ub/sr)\n')

##################################################### TRANSFORM BNGA #####################################################

def Transform_BnGa():
	locOutputFile = open('bnga_14.txt', 'w')
	Print_Header(locOutputFile)

	locInputDir = "/Users/pmatt/Scratch/clas/new_ppm_data/other_predictions/bnga_predictions/";
	for locEBin in range (1, 158):
		locInputFileName = locInputDir + "n" + ("%03d" % locEBin) + ".dat"
		locInputFile = open(locInputFileName, 'r')

		if VERBOSE > 0:
			print "file num = " + str(locEBin)

		# loop over lines
		locBeamEnergy = 0.0
		locLineCount = 0
		for locLine in locInputFile:
			locLineCount += 1
			locLineList = locLine.split()
			if len(locLineList) == 0:
				continue

			# get the beam energy
			if locLineCount == 1:
				if locEBin < 57:
					locBeamEnergy = float(locLineList[1])/1000.0
				else:
					locBeamEnergy = float(locLineList[0][2:])/1000.0
				continue

			# skip the line
			if (locLineCount == 2) or (locLineCount == 3):
				continue

			# read values
			locCosTheta = float(locLineList[0])
			if(Filter_Angle(locCosTheta) == False):
				continue
			locCrossSection = float(locLineList[1])

			# write out data
			locDataPointString = str(locBeamEnergy) + "   " + str(locCosTheta) + "   " + str(locCrossSection)
			locDataPointString += "   0.0   0.0   0.0   0.0   0.0   0.0   0.0\n" # no uncertainties provided
			locOutputFile.write(locDataPointString)

		# close input file
		locInputFile.close()

	# close output file
	locInputFile.close()

##################################################### TRANSFORM MAID #####################################################

def Transform_MAID07():
	locOutputFile = open('maid_07.txt', 'w')
	Print_Header(locOutputFile)

	locInputDir = "/Users/pmatt/Scratch/clas/new_ppm_data/other_predictions/maid07_predictions/";
	for locEBin in range (1, 115):
		locInputFileName = locInputDir + "m" + ("%03d" % locEBin) + ".dat"
		locInputFile = open(locInputFileName, 'r')

		if VERBOSE > 0:
			print "file num = " + str(locEBin)

		# loop over lines
		locBeamEnergy = 0.0
		locLineCount = 0
		for locLine in locInputFile:
			locLineCount += 1
			locLineList = locLine.split()
			if len(locLineList) == 0:
				continue

			# get the beam energy
			if locLineCount == 2:
				locBeamEnergy = float(locLineList[5])/1000.0
				continue

			# skip the line
			if (locLineCount == 1) or (locLineCount == 3):
				continue

			# read values
			locCosTheta = float(locLineList[0])
			if(Filter_Angle(locCosTheta) == False):
				continue
			locCrossSection = float(locLineList[1])
			locCrossSectionTotalError = float(locLineList[2])

			# write out data
			locDataPointString = str(locBeamEnergy) + "   " + str(locCosTheta) + "   " + str(locCrossSection)
			locDataPointString += "   0.0   0.0   " + str(locCrossSectionTotalError) + "   0.0   0.0   0.0   " + str(locCrossSectionTotalError) + "\n"
			locOutputFile.write(locDataPointString)

		# close input file
		locInputFile.close()

	# close output file
	locInputFile.close()

##################################################### TRANSFORM PR15 #####################################################

def Transform_SAIDPR15():
	locOutputFile = open('said_pr15.txt', 'w')
	Print_Header(locOutputFile)

	locInputDir = "/Users/pmatt/Scratch/clas/new_ppm_data/other_predictions/pr15_predictions/";
	for locEBin in range (1, 158):
		locInputFileName = locInputDir + "p" + ("%03d" % locEBin) + ".dat"
		locInputFile = open(locInputFileName, 'r')

		if VERBOSE > 0:
			print "file num = " + str(locEBin)

		# loop over lines
		locBeamEnergy = 0.0
		locLineCount = 0
		for locLine in locInputFile:
			locLineCount += 1
			locLineList = locLine.split()
			if len(locLineList) == 0:
				continue

			# get the beam energy
			if locLineCount == 2:
				locBeamEnergy = float(locLineList[5])/1000.0
				continue

			# skip the line
			if (locLineCount == 1) or (locLineCount == 3):
				continue

			# read values
			locCosTheta = float(locLineList[0])
			if(Filter_Angle(locCosTheta) == False):
				continue
			locCrossSection = float(locLineList[1])
			locCrossSectionTotalError = float(locLineList[2])

			# write out data
			locDataPointString = str(locBeamEnergy) + "   " + str(locCosTheta) + "   " + str(locCrossSection)
			locDataPointString += "   0.0   0.0   " + str(locCrossSectionTotalError) + "   0.0   0.0   0.0   " + str(locCrossSectionTotalError) + "\n"
			locOutputFile.write(locDataPointString)

		# close input file
		locInputFile.close()

	# close output file
	locInputFile.close()

##################################################### TRANSFORM MA27 #####################################################

def Transform_SAIDMA27():
	locOutputFile = open('said_ma27.txt', 'w')
	Print_Header(locOutputFile)

	locInputDir = "/Users/pmatt/Scratch/clas/new_ppm_data/new_SAID_solution/";
	for locEBin in range (1, 158):
		locInputFileName = locInputDir + "v" + ("%03d" % locEBin) + ".dat"
		locInputFile = open(locInputFileName, 'r')

		if VERBOSE > 0:
			print "file num = " + str(locEBin)

		# loop over lines
		locBeamEnergy = 0.0
		locLineCount = 0
		for locLine in locInputFile:
			locLineCount += 1
			locLineList = locLine.split()
			if len(locLineList) == 0:
				continue

			# get the beam energy
			if locLineCount == 2:
				locBeamEnergy = float(locLineList[5])/1000.0
				continue

			# skip the line
			if (locLineCount == 1) or (locLineCount == 3):
				continue

			# read values
			locCosTheta = float(locLineList[0])
			if(Filter_Angle(locCosTheta) == False):
				continue
			locCrossSection = float(locLineList[1])
			locCrossSectionTotalError = float(locLineList[2])

			# write out data
			locDataPointString = str(locBeamEnergy) + "   " + str(locCosTheta) + "   " + str(locCrossSection)
			locDataPointString += "   0.0   0.0   " + str(locCrossSectionTotalError) + "   0.0   0.0   0.0   " + str(locCrossSectionTotalError) + "\n"
			locOutputFile.write(locDataPointString)

		# close input file
		locInputFile.close()

	# close output file
	locInputFile.close()

################################################## TRANSFORM WORLD DATA ##################################################

def Transform_WorldData():
	locInputDir = "/Users/pmatt/Scratch/clas/new_ppm_data/";
	locInputFileName = locInputDir + "world_data.dat"
	locInputFile = open(locInputFileName, 'r')

	# loop over lines
	locBeamEnergy = 0.0
	locNewSetFlag = False
#	locFileDictionary = defaultdict(file) # file name string, file
	locOutputFile = file('dummy.txt', 'w')
	locFileDictionary = {}
	locFileDictionary['dummy.txt'] = locOutputFile
	locDataPointList = [] # data is in reverse order: must switch
	for locLine in locInputFile:
		locLineList = locLine.split()
		if len(locLineList) == 0:
			continue

		# extract, build dataset string
		if locNewSetFlag:
			locParenthesisIndex = locLineList[2].find('(')
			locFileName = locLineList[0] + "_" + locLineList[1] + "_" + locLineList[2][0:locParenthesisIndex] + ".txt"
			locFileName = locFileName.replace(',', '')
			if VERBOSE > 0:
				print "file name, energy = " + locFileName + ", " + str(locBeamEnergy) + "\n"
			if locFileName in locFileDictionary:
				locOutputFile = locFileDictionary[locFileName]
			else:
				locOutputFile = file(locFileName, 'w')
				Print_Header(locOutputFile)
				locFileDictionary[locFileName] = locOutputFile
			locNewSetFlag = False
			continue

		# see if is new dataset
		if len(locLineList) > 3:
			#beginning of new dataset for this energy: extract beam energy
			locBeamEnergy = float(locLineList[0])/1000.0
			locNewSetFlag = True

			# write out previous set's data
			for locDataPointString in locDataPointList:
				locOutputFile.write(locDataPointString)
			locDataPointList = []
			continue

		# read values
		locCosTheta = math.cos(float(locLineList[0])*math.pi/180.0)
		if(abs(locCosTheta) < 1E-7):
			locCosTheta = 0.0
		if(Filter_Angle(locCosTheta) == False):
			continue
		locCrossSection = float(locLineList[1])
		locCrossSectionTotalError = float(locLineList[2])

		# save data point
		locDataPointString = str(locBeamEnergy) + "   " + str(locCosTheta) + "   " + str(locCrossSection)
		locDataPointString += "   0.0   0.0   " + str(locCrossSectionTotalError) + "   0.0   0.0   0.0   " + str(locCrossSectionTotalError) + "\n"
		locDataPointList.insert(0, locDataPointString)

	# write out last data point
	for locDataPointString in locDataPointList:
		locOutputFile.write(locDataPointString)

	# close input file
	locInputFile.close()

	# close output files
	for locFileName in locFileDictionary:
		locFileDictionary[locFileName].close()

################################################## TRANSFORM g10 NO FSI ##################################################

def Transform_g10NoFSI():
	locOutputFile = open('CLAS_g10_noFSI.txt', 'w')
	Print_Header(locOutputFile)

	locInputDir = "/Users/pmatt/Scratch/clas/prc/rootfiles/crossfiles/";
	locInputFileName = locInputDir + "CLAS_PPM_Cross_WeiChen.txt"
	locInputFile = open(locInputFileName, 'r')

	# loop over lines
	locPreviousBeamEnergy = 0.0
	locDataPointList = [] # data is in reverse order: must switch
	for locLine in locInputFile:
		locLineList = locLine.split()
		if len(locLineList) == 0:
			continue

		# read values
		locBeamEnergy = float(locLineList[0])/1000.0
		locCosTheta = math.cos(float(locLineList[1])*math.pi/180.0)
		if(abs(locCosTheta) < 1E-7):
			locCosTheta = 0.0
		if(Filter_Angle(locCosTheta) == False):
			continue
		locCrossSection = float(locLineList[2])
		locCrossSectionStatError = float(locLineList[3])

		# fix errors
		if locBeamEnergy <= 1.2: 
			locEDependentScaleError = locCrossSection*11.54/100.0;
		elif locBeamEnergy <= 1.4: 
			locEDependentScaleError = locCrossSection*11.59/100.0;
		elif locBeamEnergy <= 1.6: 
			locEDependentScaleError = locCrossSection*11.65/100.0;
		elif locBeamEnergy <= 1.8: 
			locEDependentScaleError = locCrossSection*11.73/100.0;
		elif locBeamEnergy <= 2.0: 
			locEDependentScaleError = locCrossSection*11.82/100.0;
		elif locBeamEnergy <= 2.2: 
			locEDependentScaleError = locCrossSection*12.03/100.0;
		elif locBeamEnergy <= 2.4: 
			locEDependentScaleError = locCrossSection*12.23/100.0;
		elif locBeamEnergy <= 2.6: 
			locEDependentScaleError = locCrossSection*12.49/100.0;
		elif locBeamEnergy <= 2.8: 
			locEDependentScaleError = locCrossSection*12.67/100.0;
		elif locBeamEnergy <= 3.0: 
			locEDependentScaleError = locCrossSection*12.84/100.0;
		elif locBeamEnergy <= 3.2: 
			locEDependentScaleError = locCrossSection*13.09/100.0;
		elif locBeamEnergy <= 3.4: 
			locEDependentScaleError = locCrossSection*13.28/100.0;
		else:
			locEDependentScaleError = locCrossSection*13.70/100.0;

		# see if is new beam energy
		if abs(locBeamEnergy - locPreviousBeamEnergy) > 0.001:
			#beginning of new dataset for this energy: extract beam energy
			locPreviousBeamEnergy = locBeamEnergy

			# write out previous set's data
			for locDataPointString in locDataPointList:
				locOutputFile.write(locDataPointString)
			locDataPointList = []

		# save data point
		locCrossSectionTotalError = math.sqrt(locCrossSectionStatError*locCrossSectionStatError + locEDependentScaleError*locEDependentScaleError)
		locDataPointString = str(locBeamEnergy) + "   " + str(locCosTheta) + "   " + str(locCrossSection) + "   "
		locDataPointString += str(locCrossSectionStatError) + "   0.0   " + str(locCrossSectionStatError) + "   "
		locDataPointString += str(locEDependentScaleError) + "   0.0   " + str(locEDependentScaleError) + "   "
		locDataPointString += str(locCrossSectionTotalError) + "\n"
		locDataPointList.insert(0, locDataPointString)

	# write out last data point
	for locDataPointString in locDataPointList:
		locOutputFile.write(locDataPointString)

	# close files
	locInputFile.close()
	locOutputFile.close()

#################################################### TRANSFORM g13 FSI ###################################################

def Transform_g13FSI():
	locOutputFile = open('CLAS_g13_FSI.txt', 'w')
	Print_Header(locOutputFile)

	locInputFile_noFSI = open("/Users/pmatt/Scratch/clas/prc/data_files/CLAS_g13_noFSI.txt", 'r')
	locInputFile_noFSI.readline() #skip header line
	locLastFilePosition = locInputFile_noFSI.tell();
	locInputDir = "/Users/pmatt/Scratch/clas/new_ppm_data/my_data/temp/";

	for locEBin in range (1, 158):
		locInputFileName_FSI = locInputDir + "a" + ("%03d" % locEBin) + ".dat"
		locInputFile_FSI = open(locInputFileName_FSI, 'r')

		if VERBOSE > 0:
			print "file num = " + str(locEBin)

		# loop over igor e-bins, collecting FSI data
		locLineCount = 0
		locDataPointList = [] # Igor data is in reverse order: must switch
		for locLine in locInputFile_FSI:
			locLineCount += 1
			locLineList = locLine.split()
			if len(locLineList) == 0:
				continue

			# skip the line
			if (locLineCount == 2) or (locLineCount == 1):
				continue

			# read and save values
			locCrossSection = float(locLineList[1])
			locCrossSectionAngleError = float(locLineList[2])
			locTuple = (locCrossSection, locCrossSectionAngleError)
			locDataPointList.append(locTuple) #CHANGE: SAVE PAIR OF FLOATS

		# close input file
		locInputFile_FSI.close()

		# now, loop over noFSI data, correcting them and saving them
		locPreviousBeamEnergy = 0.0
		locLineCount = 0
		while True:
			locLastFilePosition = locInputFile_noFSI.tell();
			locLine = locInputFile_noFSI.readline()
			locLineList = locLine.split()
			if len(locLineList) == 0:
				break

			# read in the beam energy
			locBeamEnergy = float(locLineList[0])
			if (abs(locBeamEnergy - locPreviousBeamEnergy) > 0.001) and (locLineCount != 0) :
				#beginning of new energy bin: stop
				locPreviousBeamEnergy = locBeamEnergy
				locInputFile_noFSI.seek(locLastFilePosition) #backtrack a line
				break

			locPreviousBeamEnergy = locBeamEnergy

			# read in the rest of the values
			locCosTheta = float(locLineList[1])
			if(Filter_Angle(locCosTheta) == False):
				continue
			locCrossSection = float(locLineList[2])
			locCrossSectionError_Stat = float(locLineList[3])
			locCrossSectionError_Syst = float(locLineList[4])
			locCrossSectionError_Angle = float(locLineList[5])
			locCrossSectionError_Scale_Dependent = float(locLineList[6])
			locCrossSectionError_Scale_Independent = float(locLineList[7])
			locCrossSectionError_Scale = float(locLineList[8])
			locCrossSectionError_Total = float(locLineList[9])

			# Apply FSI corrections # FSI_cross = no_FSI_cross * scale_factor
			locScaleFactor = locDataPointList[locLineCount][0] / locCrossSection # scale_factor = FSI_cross / no_FSI_cross
			# Update value
			locCrossSection = locDataPointList[locLineCount][0]
			# Update errors: Fixed %
			locCrossSectionError_Stat *= locScaleFactor
			# new_syst_error = sqrt(new_angle_error*new_angle_error - new_stat_error*new_stat_error)
			locCrossSectionError_Syst = math.sqrt(locDataPointList[locLineCount][1]*locDataPointList[locLineCount][1] - locCrossSectionError_Stat*locCrossSectionError_Stat)
			locCrossSectionError_Angle = locDataPointList[locLineCount][1]
			locCrossSectionError_Scale_Dependent *= locScaleFactor
			locCrossSectionError_Scale_Independent *= locScaleFactor
			locCrossSectionError_Scale *= locScaleFactor

			# combined errors
			locCrossSectionError_Total = math.sqrt(locCrossSectionError_Angle*locCrossSectionError_Angle + locCrossSectionError_Scale*locCrossSectionError_Scale)

			# Save the new data
			locDataPointString = str(locBeamEnergy) + "   " + str(locCosTheta) + "   " + str(locCrossSection) + "   "
			locDataPointString += str(locCrossSectionError_Stat) + "   " + str(locCrossSectionError_Syst) + "   " + str(locCrossSectionError_Angle) + "   "
			locDataPointString += str(locCrossSectionError_Scale_Dependent) + "   " + str(locCrossSectionError_Scale_Independent) + "   " + str(locCrossSectionError_Scale) + "   "
			locDataPointString += str(locCrossSectionError_Total) + "\n"
			locOutputFile.write(locDataPointString)

			locLineCount += 1

	# close files
	locInputFile_noFSI.close()
	locOutputFile.close()

########################################################## MAIN ##########################################################

def Filter_Angle(locCosTheta):
	if(MAX_COSTHETADIFF < 0.0):
		return True;
	if(abs(locCosTheta - 0.6) < MAX_COSTHETADIFF):
		return True;
	elif(abs(locCosTheta - 0.3) < MAX_COSTHETADIFF):
		return True;
	elif(abs(locCosTheta - 0.0) < MAX_COSTHETADIFF):
		return True;
	elif(abs(locCosTheta + 0.35) < MAX_COSTHETADIFF):
		return True;
	return False;

########################################################## MAIN ##########################################################

def main(argv):
	Transform_BnGa()
	Transform_MAID07()
	Transform_SAIDPR15()
	Transform_SAIDMA27()
	Transform_WorldData()
	Transform_g10NoFSI()
	Transform_g13FSI()

if __name__ == "__main__":
	main(sys.argv[1:])