#!/bin/env python3 import math import numpy as np from array import array from readMDA import readMDA import ROOT import ctypes import sys, re from optparse import OptionParser class DistributionModel(): def __call__(self, xList, pList): x=xList[0] y=xList[1] Bkg = pList[0] Amp = pList[1] sigmaY = pList[2] y0 = pList[3] + pList[4] * x 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, 0 ) # Background fitFun.SetParameter(1, histo.GetMaximum() ) # Amplitude fitFun.SetParameter(2, histo.GetRMS(2) ) # Sigma fitFun.SetParameter(3, histo.GetYaxis().GetBinCenter( binY.value ) ) # Mean fitFun.SetParameter(4, 0 ) # Slope of the peak histo.Fit( fitFun, "S" ) return def getScanData( mdaFileName ): mdaData = readMDA(mdaFileName, 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 returnMap = {"uMotor" : upstreamMotor, "mMotor" : middleMotor, "dMotor" : downstreamMotor, "tagm" : tagmSettings, "detectors": detectorArray } return returnMap def openRootFile( fileName = "ROOTFile.root"): rootFile = ROOT.TFile.Open(fileName, "RECREATE") return rootFile def closeRootFile( rootFile ): if rootFile.IsOpen(): rootFile.Write() rootFile.Close() return def createNtuple( scanMap ): upstreamMotor = scanMap["uMotor"] middleMotor = scanMap["mMotor"] downstreamMotor = scanMap["dMotor"] tagmSettings = scanMap["tagm"] detectorArray = scanMap["detectors"] numberofRowsTAGM = len(tagmSettings.data[0] ) numberOfDetectors = len(detectorArray) nMotorPositions = len(upstreamMotor.data) numberOfColumnsTAGM = numberOfDetectors - 1 ROOT.gROOT.ProcessLine(f""" struct MultiLeafData {{ Int_t id[{numberOfColumnsTAGM}]; Double_t rate[{numberOfColumnsTAGM}]; }}; """) tagmDetPattern = r"^TAGM_T_(\d+)_scaler_latch$" rootTree = ROOT.TTree("ScanTree", "TAGM Scan Tree") branchVecDict = {} for positionerKey in [ "uMotor", "mMotor", "dMotor", "tagm" ]: branchName = scanMap[positionerKey].readback_name.decode("utf-8").replace(":", "_") branchVecDict[branchName] = ROOT.std.vector('double')() rootTree.Branch(branchName, branchVecDict[branchName]) print( "Created branch {0}".format(branchName)) leafList = f"id[{numberOfColumnsTAGM}]/I:rate[{numberOfColumnsTAGM}]/D" tagmDetectorData = ROOT.MultiLeafData() rootTree.Branch("DetRates", tagmDetectorData, leafList ) print( "Created branch {0}".format("DetRates") ) for detector in detectorArray: detName = detector.name.decode("utf-8").replace(":", "_") patternMatch = re.match(tagmDetPattern, detName) if not patternMatch : branchVecDict[detName] = ROOT.std.vector('double')() # Pass the vector object directly to the branch rootTree.Branch(detName, branchVecDict[detName]) print( "Created branch {0}".format(detName)) for motorIndex in range(0, nMotorPositions ) : uMotorPosition = upstreamMotor.data[motorIndex] mMotorPosition = middleMotor.data[motorIndex] dMotorPosition = downstreamMotor.data[motorIndex] branchName = upstreamMotor.readback_name.decode("utf-8").replace(":", "_") branchVecDict[branchName].clear() branchVecDict[branchName].push_back( uMotorPosition ) branchName = middleMotor.readback_name.decode("utf-8").replace(":", "_") branchVecDict[branchName].clear() branchVecDict[branchName].push_back( mMotorPosition ) branchName = downstreamMotor.readback_name.decode("utf-8").replace(":", "_") branchVecDict[branchName].clear() branchVecDict[branchName].push_back( dMotorPosition ) for tagmRowIndex in range( 0, len(tagmSettings.data[motorIndex]) ) : tagmRow = tagmSettings.data[motorIndex][tagmRowIndex] branchName = tagmSettings.readback_name.decode("utf-8").replace(":", "_") branchVecDict[branchName].clear() branchVecDict[branchName].push_back( tagmRow ) for tagmIndex in range( numberOfColumnsTAGM ): tagmDetectorData.id[tagmIndex] = 0 tagmDetectorData.rate[tagmIndex] = 0 tagmIndex = 0 for detectorIndex in range( 0, numberOfDetectors ): detectorName = detectorArray[detectorIndex].name.decode("utf-8").replace(":", "_") detectorValue = detectorArray[detectorIndex].data[motorIndex][tagmRowIndex] patternMatch = re.match(tagmDetPattern, detectorName) if patternMatch : tagmNumber = int( patternMatch.group(1) ) # print ( f"Detector name is {detectorName}, TAGM counter # {tagmNumber}") tagmDetectorData.id[tagmIndex] = tagmNumber tagmDetectorData.rate[tagmIndex] = detectorValue tagmIndex += 1 else : branchName = detectorName # branchVecDict[branchName].clear() branchVecDict[branchName].push_back(detectorValue) rootTree.Fill() rootFile.Write() return rootTree def simulateScanData( scanMap ) : upstreamMotor = scanMap["uMotor"] tagmSettings = scanMap["tagm"] detectorArray = scanMap["detectors"] numberofRowsTAGM = len(tagmSettings.data[0] ) numberOfColumnsTAGM = len(detectorArray) nMotorPositions = len(upstreamMotor.data) minMotor = min( upstreamMotor.data ) - ( max(upstreamMotor.data) - min(upstreamMotor.data) ) / (2. * (nMotorPositions-1)) maxMotor = max( upstreamMotor.data ) + ( max(upstreamMotor.data) - min(upstreamMotor.data) ) / (2. * (nMotorPositions-1)) NumberOfSimulatedEvents = 5000 ROOT.gRandom = ROOT.TRandom3() meanHeight = 10.0 sigmaHeight = 2.0 simHisto = [None] * numberofRowsTAGM for tagmRowIndex in range( 0, numberofRowsTAGM ) : # simHisto[tagmRowIndex] = [None] * numberOfColumnsTAGM histName = "Sim_{0}".format(tagmRowIndex) histTitle = "Simulation TAGM Row {0}".format(tagmRowIndex) simHisto[tagmRowIndex] = ROOT.TH2D( histName, histTitle, numberOfColumnsTAGM, 0.5, numberOfColumnsTAGM+0.5, nMotorPositions, minMotor, maxMotor ) for tagmColumn in range( 1, numberOfColumnsTAGM+1 ): meanHeight = .5 + (0.6/60.) * (tagmColumn-1) for iEvt in range(NumberOfSimulatedEvents): simData = ROOT.gRandom.Gaus( meanHeight, sigmaHeight ) simHisto[tagmRowIndex].Fill( tagmColumn, simData ) return simHisto def fillScanHistos( scanMap, simFlag = False ) : if simFlag: simHisto = simulateScanData(scanMap) upstreamMotor = scanMap["uMotor"] tagmSettings = scanMap["tagm"] detectorArray = scanMap["detectors"] numberofRowsTAGM = len(tagmSettings.data[0] ) numberOfColumnsTAGM = len(detectorArray) nMotorPositions = len(upstreamMotor.data) minMotor = min( upstreamMotor.data ) - ( max(upstreamMotor.data) - min(upstreamMotor.data) ) / (2. * (nMotorPositions-1)) maxMotor = max( upstreamMotor.data ) + ( max(upstreamMotor.data) - min(upstreamMotor.data) ) / (2. * (nMotorPositions-1)) rowHisto = [None] * numberofRowsTAGM for tagmColumn in range( 0, numberOfColumnsTAGM ): 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() if( simFlag ) : measuredData = simHisto[tagmRowIndex].GetBinContent(tagmColumn+1, motorIndex+1) else : measuredData = detectorArray[tagmColumn].data[motorIndex][tagmRowIndex] # print( "Data in [[{0}][{1}] is {2}".format(tagmColumn, motorPosition, measuredData ) ) rowHisto[tagmRowIndex].Fill(tagmColumn, motorPosition, measuredData ) return rowHisto def drawHistos( rowHisto ): numberofRowsTAGM = len(rowHisto) 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() profileHisto = [None] * numberofRowsTAGM 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(), 5) fitFun.SetParNames("Background", "Amplitude", "Sigma_Y", "P0", "P1" ) canvas[rowIndex] = ROOT.TCanvas() # canvas[rowIndex].Divide(2, 1) # rightPad.cd( rowIndex+1 ) rowHisto[rowIndex].SetContour(50) # canvas[rowIndex].cd(1) # canvas[rowIndex].GetPad(1).SetLogz(1) # fitGauss2D( rowHisto[rowIndex], fitFun ) # rowHisto[rowIndex].DrawCopy("LEGO2") # canvas[rowIndex].cd(2) # profileHisto[rowIndex] = rowHisto[rowIndex].ProfileX() # profileHisto[rowIndex].Draw() # profileHisto[rowIndex].Fit( "pol1", "Q") rowHisto[rowIndex].DrawCopy("COLZ") canvas[rowIndex].Update() return canvas # This gets executed if used as main if __name__ == "__main__": parser = OptionParser(usage = "usage: %prog [options] ") parser.add_option( "-d", "--dir", action="store", dest="dir", type="string", metavar="InputDir", default="/gluex/data/MicroscopeScans/sscanData/", help="Define directory for MDA file" ) parser.add_option( "-f", "--file", action="store", dest="file", type="string", metavar="InputFile", default="tagm_scan_0010.mda", help="Define file name for MDA file" ) parser.add_option( "-s", "--show", action="store_true", dest="show", default=False, help="Plot the graphs on ROOT canvas ?" ) parser.add_option( "-t", "--tree", action="store_true", dest="tree", default=False, help="Make the ROOT tree ?" ) (opts, arguments) = parser.parse_args( sys.argv[1:] ) mdaFullFileName = opts.dir + "/" + opts.file rootFileName = mdaFullFileName + ".root" scanMap = getScanData( mdaFullFileName ) if opts.tree: rootFile = openRootFile( rootFileName ) rootFile.cd() #rowHistograms = fillScanHistos( scanMap, simFlag=True ) rowHistograms = fillScanHistos( scanMap, simFlag=False ) if opts.show : canvasList = drawHistos( rowHistograms ) if opts.tree: rootFile.cd() rootTree = createNtuple( scanMap ) rootTree.Write() closeRootFile( rootFile ) input("Press Enter to close the windows...")