#!/bin/env python3 import math import numpy as np from array import array #import matplotlib.pyplot as plt from readMDA import readMDA import ROOT import ctypes class DistributionModel(): def __call__(self, xList, pList): x=xList[0] y=xList[1] Amp = pList[0] y0 = pList[1] sigmaY = pList[2] Bkg = pList[3] dY = (y - y0 ) / sigmaY value = Bkg + Amp * math.exp(-0.5 * dY * dY ) return value def fitGauss2D( histo, fitFun ): binX = ctypes.c_int(0) binY = ctypes.c_int(0) binZ = ctypes.c_int(0) histo.GetMaximumBin( binX, binY, binZ ) print ("Bin Center is {0} or {1}".format(binY, binY.value)) fitFun.SetParameter(0, histo.GetMaximum() ) # Amplitude fitFun.SetParameter(1, histo.GetYaxis().GetBinCenter( binY.value ) ) # Mean fitFun.SetParameter(2, histo.GetRMS(2) ) # Sigma fitFun.SetParameter(3, 0 ) # Background histo.Fit( fitFun, "S V" ) # 1. Load the MDA file # d will be a list containing the dimension structures (0 = headers, 1 = outer scan, 2 = inner scan) mda_file = '/gluex/data/MicroscopeScans/sscanData/tagm_scan_0004.mda' mdaData = readMDA(mda_file, maxdim=2) innerSettingsScan = mdaData[2] outerMotorScan = mdaData[1] upstreamMotor = outerMotorScan.p[0] middleMotor = outerMotorScan.p[1] downstreamMotor = outerMotorScan.p[2] tagmSettings = innerSettingsScan.p[0] detectorArray = innerSettingsScan.d #print ( "Upstream motor name is ", upstreamMotor.name ) #print ( "Middle motor name is ", middleMotor.name ) #print ( "Downstream motor name is ", downstreamMotor.name ) # 2. Extract Positioners (Motors) for 2D # d[1] is the outer loop (e.g., Y-axis), d[2] is the inner loop (e.g., X-axis) # x_axis_positions = mdaData[1].p[0].data # u_axis_positions = mdaData[1].p[1].data # y_axis_positions = mdaData[2].p[0].data # detector_index = 0 # z_axis_data = detectorArray[detector_index].data # print( "This is X \n") # print( x_axis_positions ) # print( "This is Y \n") # print( y_axis_positions ) # print( "This is U \n") # print( u_axis_positions ) # 3. Extract Detector Data # Detectors are generally 2D arrays shaped as (Outer points, Inner points) #det_data = detectorArray[detector_index].data #motorError = 0.010 #numberOfGraphPoints = len(tagmSettings.data) * len(tagmSettings.data[0]) #xPoints = np.tile(upstreamMotor.data, (len(detectorSettings.data), 1) ) #xPoints2D = [list(upstreamMotor.data) for _ in range(len(detectorSettings.data))] # xPoints2D = [[item] * len(tagmSettings.data[0]) for item in upstreamMotor.data] # xPoints = array('d', [element for sublist in xPoints2D for element in sublist] ) # xPointErrors = array('d', [motorError for element in xPoints] ) # yPoints2D = tagmSettings.data # yPoints = array('d', [element for sublist in yPoints2D for element in sublist] ) # yPointErrors = array('d', [0 for element in yPoints] ) # zPoints2D = det_data # zPoints = array('d', [element for sublist in zPoints2D for element in sublist] ) # zPointErrors = array('d', [math.sqrt(element) for element in zPoints] ) # print( "This is X-points \n " ) # print( xPoints ) # print( "This is Y-points \n " ) # print( yPoints ) # print( "This is Z-points \n " ) # print( zPoints ) #print( "Number of graph points is " , numberOfGraphPoints ) numberofRowsTAGM = len(tagmSettings.data[0] ) numberOfColumnsTAGM = len(detectorArray) nMotorPositions = len(outerMotorScan.p[0].data) minMotor = min( outerMotorScan.p[0].data ) - ( max(outerMotorScan.p[0].data) - min(outerMotorScan.p[0].data) ) / (2. * (nMotorPositions-1)) maxMotor = max( outerMotorScan.p[0].data ) + ( max(outerMotorScan.p[0].data) - min(outerMotorScan.p[0].data) ) / (2. * (nMotorPositions-1)) ROOT.gRandom = ROOT.TRandom3() meanHeight = 10.0 sigmaHeight = 3.0 rowHisto = [None] * numberofRowsTAGM for tagmColumn in range( 0, numberOfColumnsTAGM ): simHisto = ROOT.TH1D( "SIM", "SIM", nMotorPositions, minMotor, maxMotor ) for iEvt in range(10000): simData = ROOT.gRandom.Gaus( meanHeight, sigmaHeight ) simHisto.Fill( simData ) for motorIndex in range(0, nMotorPositions ) : motorPosition = upstreamMotor.data[motorIndex] for tagmRowIndex in range( 0, len(tagmSettings.data[motorIndex]) ) : tagmRow = tagmSettings.data[motorIndex][tagmRowIndex] if( rowHisto[tagmRowIndex] == None ) : histoName = "Row_{0}".format(tagmRow) histoTitle = "TAGM Row # {0}".format(tagmRow) rowHisto[tagmRowIndex] = ROOT.TH2D( histoName, histoTitle, numberOfColumnsTAGM, 0.5, numberOfColumnsTAGM+0.5, nMotorPositions, minMotor, maxMotor ) # rowHisto[tagmRowIndex].Print() # rowHisto[tagmRowIndex].Dump() # measuredData = detectorArray[tagmColumn].data[motorIndex][tagmRowIndex] measuredData = simHisto.GetBinContent(motorIndex+1) # print( "Data in [[{0}][{1}] is {2}".format(tagmColumn, motorPosition, measuredData ) ) rowHisto[tagmRowIndex].Fill(tagmColumn, motorPosition, measuredData ) ROOT.gStyle.SetPalette(ROOT.kRainBow) ROOT.gStyle.SetPalette(ROOT.kTemperatureMap) ROOT.gStyle.SetPalette(ROOT.kSunset) ROOT.gStyle.SetPalette(ROOT.kColorPrintableOnGrey) ROOT.gStyle.SetPalette(ROOT.kBlueRedYellow) ROOT.gStyle.SetOptStat(1) ROOT.gStyle.SetOptFit(1) ROOT.gStyle.SetOptTitle(1) model = DistributionModel() canvas = [None] * numberofRowsTAGM for rowIndex in range( 0, numberofRowsTAGM ): fitFun = ROOT.TF2("gauss2d", model, rowHisto[rowIndex].GetXaxis().GetXmin(), rowHisto[rowIndex].GetXaxis().GetXmax(), rowHisto[rowIndex].GetYaxis().GetXmin(), rowHisto[rowIndex].GetYaxis().GetXmax(), 4) fitFun.SetParNames("Amplitude", "Mean_Y", "Sigma_Y", "Bkg" ) canvas[rowIndex] = ROOT.TCanvas() # rightPad.cd( rowIndex+1 ) rowHisto[rowIndex].SetContour(50) rowHisto[rowIndex].Draw("LEGO2") fitGauss2D( rowHisto[rowIndex], fitFun ) canvas[rowIndex] .Update() # # We pass the NumPy arrays directly; PyROOT automatically handles the memory buffers. # graph = ROOT.TGraph2DErrors( numberOfGraphPoints, xPoints, yPoints, zPoints, xPointErrors, yPointErrors, zPointErrors ) # # Set some visual properties # graphTitle = "Result of the TAGM scan " + mdaData[1].scan_name.decode('utf-8') + " on " + (mdaData[1].time.decode('utf-8') ) # graph.SetTitle(graphTitle) # graph.SetMarkerStyle(20) # graph.SetMarkerColor(ROOT.kBlue) # graph.SetLineColor(ROOT.kRed) # graph.Draw("TRI1 P0") input("Press Enter to continue...")