combinedMassScaleAndExtractSumSpecAPi.jl

push!(LOAD_PATH, pwd())


import HDF5
using PyCall
pygui(:tk) # :tk, :gtk3, :gtk, :qt5, :qt4, :qt, or :wx
import PyPlot
import APiTOFFunctions
import InterpolationFunctions

function correctMassScaleAndExtractSumSpecAPi(
  filepath,
  masslistMasses,
  masslistElements,
  masslistElementsMasses,
  masslistCompositions,
  referenceFile,
  calibRegions; # Regions which are pattern matched for mass shifts
  filefilterRegexp=r"\.h5$",
  outputfilename="results/_result.hdf5",
  firstNFiles=0, # only analyze the first N files, set to 0 for all files
  lastNFiles=0, # only analyze the last N files, set to 0 for all files
  debuglevel=3,
  searchWidth = 0.7, # Width of Mass Scale Search Regions in AMU
  dynamicMassScaleCorrection = true,
  recalibInterval = 60,
  peakshapeRegions = 10,
  createTotalAvg = true, # Needed for manualPeakFitter, not needed for timetraces
  onlyUseAverages = false, # Fast mode using only total file averages instead of individual spectra in file
  storeSubSpectra = true,
  plotControlMass = false, # plot the control mass peak for each file after correction
  testRangeStart = 137.0, # the mass shift of this region will be shown if plot control mass is set true. Should not be part of calibRegions
  testRangeEnd = 137.5,
  massBorderCalculation = 1, # How to calculate borders? 0 = Cernter -0.1 to Center + 0.4,  1 = based on resolution and peak distance, 2 = constant bin width
  binWidth = 6,
  resolution = 7500)



  tic()
  if (isfile(joinpath(filepath, outputfilename)))
    mv(joinpath(filepath, outputfilename), joinpath(filepath, "$outputfilename.bak"), remove_destination=true)
    if debuglevel > 0   println("Found existing result file, moving to '$outputfilename.bak'") end
  end

  #files = filter(r"2016-07-04.*\.h5", readdir(filepath)) #one day only
  files = filter(filefilterRegexp, readdir(filepath))




  nFiles = size(files,1)

  if firstNFiles != 0
    if length(files) > firstNFiles
      files = files[1:firstNFiles]
    end
  end

  if lastNFiles != 0
    if length(files) > lastNFiles
      files = files[length(files)-lastNFiles:end]
    end
  end

  nFiles = size(files,1)

  validFiles, timeSortIndices = APiTOFFunctions.validateHDF5Files(filepath, files)
  if !issorted(timeSortIndices)
      println("Reordering Files according to start acquisition time!")
  else
      println("File order seems fine, continuing...")
  end
  files = validFiles[timeSortIndices]
  nFiles = size(files,1)


  if (referenceFile == "")
    referenceFile = joinpath(filepath, files[1])
    if (debuglevel > 1)
      println("Reference File Autodetected (using first file in folder): $referenceFile")
    end
  end


  referenceSpectrum = APiTOFFunctions.getAvgSpectrumFromFile(referenceFile)
  referenceMassScaleMode, referenceMassScaleParameters = APiTOFFunctions.getMassCalibParametersFromFile(referenceFile)
  APiTOFFunctions.setMassScaleReferenceSpectrum(referenceSpectrum, calibRegions, searchWidth, referenceMassScaleMode, referenceMassScaleParameters, plotControlMass=plotControlMass, testRangeStart=testRangeStart, testRangeEnd=testRangeEnd)
  referenceMassAxis = []
  fh = HDF5.h5open(referenceFile,"r")
  if HDF5.exists(HDF5.attrs(fh), "InstrumentType")
    lenMassAxis = length(HDF5.h5read(referenceFile, "SPECdata/AverageSpec"))
    timeBins = vcat(linspace(1,lenMassAxis,lenMassAxis))
    pV = HDF5.h5read(referenceFile,"/CALdata/Spectrum")
    p1 = pV[1,1]
    p2 = pV[2,1]
    referenceMassAxis = ((timeBins[:]-p2[1])/p1[1]).^2
  else
    referenceMassAxis = HDF5.h5read(referenceFile, "FullSpectra/MassAxis")
  end
  HDF5.close(fh)
  # Tofwerk does not update the mass axis after manual recalibration in TofDAQ Viewer.
  # ==> recalculate based on calib parameters.
  println("Recalculating mass axis")
  for i=1:length(referenceMassAxis)
      referenceMassAxis[i] = APiTOFFunctions.timebin2mass(i, referenceMassScaleMode, referenceMassScaleParameters)
  end
  #totalAvgSpectrum = SharedArray{Float64}(length(referenceSpectrum))
  #totalAvgSubSpectrum = SharedArray{Float64}(length(referenceSpectrum))
  #totalMinSpectrum = SharedArray{Float64}(length(referenceSpectrum))
  #totalMaxSpectrum = SharedArray{Float64}(length(referenceSpectrum))
  totalAvgSpectrum = Array{Float64}(length(referenceSpectrum))
  totalAvgSubSpectrum = Array{Float64}(length(referenceSpectrum))
  totalMinSpectrum = Array{Float64}(length(referenceSpectrum))
  totalMaxSpectrum = Array{Float64}(length(referenceSpectrum))

  #@sync @parallel
  for i = 1 : length(referenceSpectrum)
    totalAvgSpectrum[i] = 0
    totalAvgSubSpectrum[i] = 0
    totalMinSpectrum[i] = 1e99
    totalMaxSpectrum[i] = 0
  end



  nMasses=length(masslistMasses)
  println("Calculating Stick CPS for $nMasses masses.")
  timebinshifts = zeros(Float64,(length(calibRegions),nFiles))
  intensities = zeros(Float64,(length(calibRegions),nFiles))
  monitorTimetrace = zeros(nFiles)

  #time = SharedArray{DateTime}(nFiles)
  #stickcps=SharedArray{Float64}(nMasses,nFiles)
  time = Array{DateTime}(nFiles)
  stickcps=Array{Float64}(nMasses,nFiles)

  # Calculate Integration Borders
  if (massBorderCalculation == 0)
  mlow = masslistMasses - 0.1
  mhigh = masslistMasses + 0.4
  elseif  (massBorderCalculation == 1)
  #mlow = masslistMasses .* resolution./(resolution+0.5)
  #mhigh = masslistMasses .* resolution./(resolution-0.5)
  massborders = MasslistFunctions.IntegrationBorders(masslistMasses, resolution)
  end


  if (plotControlMass == true)
    PyPlot.figure()
  end

  if (createTotalAvg == true)
      #interpolatedSpectrum = SharedArray{Float64}(length(referenceSpectrum))
      interpolatedSpectrum = Array{Float64}(length(referenceSpectrum))
  end

  ############## Check output path and remove existing files #####################
  if (!isdir("$filepath/results"))
    mkdir("$filepath/results")
  end
  outfilepath = joinpath(filepath, outputfilename)
  if isfile(outfilepath)
    rm(outfilepath)
  end

  fid = HDF5.h5open(outfilepath, "w")

  ################################################################################

  ############## open outputfile and create extendable SumSpecs Array ############
  if (createTotalAvg == true)

    dsAvgSumWidth = length(referenceSpectrum)
    dspaceAvgSumSpecs = HDF5.dataspace((dsAvgSumWidth,1)::Dims, max_dims=(dsAvgSumWidth,typemax(Int64)))
    dtypeAvgSumSpecs = HDF5.datatype(Float32)
    dsetAvgSumSpecs = HDF5.d_create(fid, "SumSpecs", dtypeAvgSumSpecs, dspaceAvgSumSpecs, "chunk", (dsAvgSumWidth,1), "compress", 3)
  end

  if (!onlyUseAverages)
    specCount = 0
    dsStickCpsWidth = nMasses
    dspaceStickCps = HDF5.dataspace((1,dsStickCpsWidth)::Dims, max_dims=(typemax(Int64),dsStickCpsWidth))
    dtypeStickCps = HDF5.datatype(Float32)
    dsetStickCps = HDF5.d_create(fid, "StickCps", dtypeStickCps, dspaceStickCps, "chunk", (1,dsStickCpsWidth), "compress", 3)
    dsetStickCpsErr = HDF5.d_create(fid, "StickCpsErr", dtypeStickCps, dspaceStickCps, "chunk", (1,dsStickCpsWidth), "compress", 3)

    dspaceTimes = HDF5.dataspace((1,)::Dims, max_dims=(typemax(Int64),))
    dtypeTimes = HDF5.datatype(Float64)
    dsetTimes = HDF5.d_create(fid, "Times", dtypeTimes, dspaceTimes, "chunk", (1,))

  end
  ################################################################################
  badFiles = Array{String}(0)

  for j=1:nFiles
    totalPath = joinpath(filepath, files[j])
    if j < nFiles
        totalPrecachePath = joinpath(filepath, files[j+1])
    else
        totalPrecachePath = ""
    end
    if debuglevel > 0   println("Processing File $j/$nFiles :  $totalPath") end

    fileIsBad = false;
    massAxis = []
    fh = HDF5.h5open(totalPath,"r")
    #println("fh= ",fh)
    if HDF5.exists(HDF5.attrs(fh), "InstrumentType")
      try obj = fh["PROCESSED"]
          massAxis = HDF5.h5read(filename, "PROCESSED/AverageSpectrum/MassAxis")
      catch
          lenMassAxis = length(HDF5.h5read(totalPath, "SPECdata/AverageSpec"))
          timeBins = vcat(linspace(1,lenMassAxis,lenMassAxis))
          pV = HDF5.h5read(totalPath,"/CALdata/Spectrum")
          p1 = pV[1,1]
          p2 = pV[2,1]
          massAxis = ((timeBins[:]-p2[1])/p1[1]).^2
      end
    else
      massAxis = HDF5.h5read(totalPath, "FullSpectra/MassAxis")
    end
    HDF5.close(fh)

    time[j] = APiTOFFunctions.getTimeFromFile(totalPath)

    avgSpectrum = APiTOFFunctions.getAvgSpectrumFromFile(totalPath)
    newParams, success, tbs, ins = APiTOFFunctions.recalibrateMassScale(avgSpectrum, referenceSpectrum, calibRegions, searchWidth, referenceMassScaleMode, referenceMassScaleParameters)

    #if debuglevel > 0   println("Processing File $totalPath") end
    if (createTotalAvg == true && !fileIsBad)
      print("Spectrum interpolation for total average... ")
      indexesExact = APiTOFFunctions.mass2timebin(referenceMassAxis, referenceMassScaleMode, newParams)
      interpolatedSpectrum = InterpolationFunctions.interpolate(indexesExact,avgSpectrum)
      totalAvgSpectrum += interpolatedSpectrum
      # calculate min and max spectra
      #@sync @parallel
      for bin=1:length(interpolatedSpectrum)
        if interpolatedSpectrum[bin] > totalMaxSpectrum[bin]
          totalMaxSpectrum[bin] = interpolatedSpectrum[bin]
        end
        if interpolatedSpectrum[bin] < totalMinSpectrum[bin]
          totalMinSpectrum[bin] = interpolatedSpectrum[bin]
        end
      end
      # Write to hdf5, line by line, so there is no limit to number of files that can fit in RAM
      if (storeSubSpectra == true)
        currDims = size(dsetAvgSumSpecs)[2]
        dsetAvgSumSpecs[:,currDims] = interpolatedSpectrum
        if (j != nFiles)
          HDF5.set_dims!(dsetAvgSumSpecs, (dsAvgSumWidth,currDims+1)::Dims)
        end
      end
      println("DONE")
    end
    if debuglevel > 1   println() end
    print("Average Stick integration... ")
    for i=(1:nMasses)
      #if debuglevel > 0   println("Processing region mass($(mlow[i]):$(mhigh[i])) --> timebin($(round(APiTOFFunctions.mass2timebin(mlow[i],massCalibMode,newParams))):$(round(APiTOFFunctions.mass2timebin(mhigh[i],massCalibMode,newParams))))") end
      if (massBorderCalculation == 2)
        centerIndex = Int64(floor(APiTOFFunctions.mass2timebin(masslistMasses[i], referenceMassScaleMode, newParams)))
        if ((centerIndex-binWidth) > 0) && ((centerIndex+binWidth)<length(avgSpectrum))
          raw = sum(view(avgSpectrum, (centerIndex-binWidth : centerIndex + binWidth)))
        else
          raw = 0
        end
      else
        subIdxStartExact=APiTOFFunctions.mass2timebin(massborders.lowMass[i],referenceMassScaleMode,newParams)
        subIdxEndExact = APiTOFFunctions.mass2timebin(massborders.highMass[i],referenceMassScaleMode,newParams)
        raw = InterpolationFunctions.interpolatedSum(subIdxStartExact,subIdxEndExact,avgSpectrum)
      end
      stickcps[i,j] = raw
    end
    println("DONE")
    if (!onlyUseAverages)
      spectrumMultFactor = APiTOFFunctions.getSpecMultiplicator(totalPath)
      println("   spectrumMultFactor = $spectrumMultFactor")
      #subSpecStickCps = SharedArray{Float32}(nMasses)
      subSpecStickCps = Array{Float32}(nMasses)
      totalSubSpectra = APiTOFFunctions.getSubSpectraCount(totalPath)

      ################## prepare arrays for averaging #########################
      fileInternalLocalAvg = APiTOFFunctions.getSubSpectrumFromFile(totalPath,1, preloadFile=totalPrecachePath)
      #println("referenceSpectrum is a $(summary(referenceSpectrum))")
      #println("SubSpec is a $(summary(fileInternalLocalAvg))")
      #println("spectrumMultFactor is a $(summary(spectrumMultFactor))")
      print("Sub spectrum Stick integration... ")

      fill!(fileInternalLocalAvg,0)
      fileInternalLocalAvgCount = 0
      print("Precalculating average for mass scale correction... ")
      ################## Get First Set of Spectra for first mass scale calib
      if totalSubSpectra >= recalibInterval
        # Prepare first calib beforehead
        fileInternalLocalAvg = APiTOFFunctions.getSubSpectrumFromFile(totalPath,1, preloadFile=totalPrecachePath)
        println(" minimum recalInterval: $(minimum([recalibInterval totalSubSpectra]))")
        for avgIdx=2:minimum([recalibInterval totalSubSpectra])
          #println("####-----------#######1 and $totalSubSpectra ## avgIdx = $avgIdx ## recalibInterval = $recalibInterval #####")
          fileInternalLocalAvg += APiTOFFunctions.getSubSpectrumFromFile(totalPath,avgIdx, preloadFile=totalPrecachePath)
        end
        newParams, success, tbs, ins = APiTOFFunctions.recalibrateMassScale(fileInternalLocalAvg, referenceSpectrum, calibRegions, searchWidth, referenceMassScaleMode, referenceMassScaleParameters)
        ######
        fill!(fileInternalLocalAvg,0)
        fileInternalLocalAvgCount = 0
      end
      #######################################################################

      ################### Loop over Subspectra Start ########################
      ################### Averaging and Mass Scale Calib ####################
      print("Recalibrating at: ")
      specCount += totalSubSpectra
      for subSpecIdx=1:totalSubSpectra
        subSpectrum = APiTOFFunctions.getSubSpectrumFromFile(totalPath,subSpecIdx, preloadFile=totalPrecachePath)
        fileInternalLocalAvg += subSpectrum
        fileInternalLocalAvgCount += 1
        #println("####-----------#######2 and $subSpecIdx/$totalSubSpectra #####-----------------#############")

        if (fileInternalLocalAvgCount > recalibInterval) || fileInternalLocalAvgCount == totalSubSpectra
          print("$subSpecIdx...")
          newParams, success, tbs, ins = APiTOFFunctions.recalibrateMassScale(fileInternalLocalAvg, referenceSpectrum, calibRegions, searchWidth, referenceMassScaleMode, referenceMassScaleParameters)

          indexesExact = APiTOFFunctions.mass2timebin(referenceMassAxis, referenceMassScaleMode, newParams)
          interpolatedSpectrum = InterpolationFunctions.interpolate(indexesExact,fileInternalLocalAvg)

          totalAvgSubSpectrum += interpolatedSpectrum*spectrumMultFactor

          fill!(fileInternalLocalAvg,0)
          fileInternalLocalAvgCount = 0
        end
        ################## Peak Integration ################################
        #@sync @parallel
        mbIndicesLowMass = APiTOFFunctions.mass2timebin(massborders.lowMass,referenceMassScaleMode,newParams)
        mbIndicesHighMass = APiTOFFunctions.mass2timebin(massborders.highMass,referenceMassScaleMode,newParams)
        Threads.@threads for i=(1:nMasses)
          if (mod(i,100) == 0)
            #println("Processing mass $(masslistMasses[i])")
          end
          if (massBorderCalculation == 2)
            centerIndex = Int64(floor(APiTOFFunctions.mass2timebin(masslistMasses[i], referenceMassScaleMode, newParams)))
            if ((centerIndex-binWidth) > 0) && ((centerIndex+binWidth)<length(avgSpectrum))
              raw = sum(view(subSpectrum, (centerIndex-binWidth : centerIndex + binWidth)))
            else
              raw = 0
            end
            subSpecStickCps[i]=raw
          else
            #if debuglevel > 0   println("Processing region mass($(mlow[i]):$(mhigh[i])) --> timebin($(round(APiTOFFunctions.mass2timebin(mlow[i],massCalibMode,newParams))):$(round(APiTOFFunctions.mass2timebin(mhigh[i],massCalibMode,newParams))))") end
            subSpecStickCps[i]= InterpolationFunctions.interpolatedSum(mbIndicesLowMass[i],mbIndicesHighMass[i],subSpectrum)
          end
        end
        rawTime = APiTOFFunctions.getSubSpectrumTimeFromFile(totalPath,subSpecIdx)
        #println("####-----------rawTime = $rawTime-----------------#############")
        currDimsTime = size(dsetTimes)[1]
        absTimeOfCurrSample = Dates.datetime2unix(time[j]) + rawTime
        dsetTimes[currDimsTime] = absTimeOfCurrSample
        dsetStickCps[currDimsTime,:] = subSpecStickCps*spectrumMultFactor
        #println("####-----------#######6+$subSpecIdx#####-----------------#############")
        if (currDimsTime == 1)
          deltaT = 5 #rawTime # Not correct???
        else
          deltaT = absTimeOfCurrSample - dsetTimes[currDimsTime-1]
        end
        #println("DeltaT: $deltaT")
        dsetStickCpsErr[currDimsTime,:] = sqrt.(abs.(subSpecStickCps./deltaT))
        if !((j==nFiles) && (subSpecIdx== APiTOFFunctions.getSubSpectraCount(totalPath)))
          HDF5.set_dims!(dsetTimes, (currDimsTime+1,)::Dims)
          HDF5.set_dims!(dsetStickCps, (currDimsTime+1,dsStickCpsWidth)::Dims)
          HDF5.set_dims!(dsetStickCpsErr, (currDimsTime+1,dsStickCpsWidth)::Dims)
        end
      end
      println("DONE")
    end
    println("###################################################################\n\n")
  end
  close(fid)

  stickcpsFlat = zeros(nMasses,nFiles)
  stickcpsErrFlat = zeros(nMasses,nFiles)

  timeFlat = zeros(nFiles)
  timeFlat = Dates.datetime2unix.(time)
  # We need more than one timestamp to calculate a difference
  if length(timeFlat) > 1
      deltaT = median(diff(timeFlat))
  else
      deltaT = 1
  end

  totalAvgSpectrumFlat = zeros(length(totalAvgSpectrum))
  totalAvgSubSpectrumFlat = zeros(length(totalAvgSubSpectrum))


  for i=1:nMasses, j=1:nFiles
    stickcpsFlat[i,j] = stickcps[i,j]
    stickcpsErrFlat[i,j] = sqrt(abs(stickcps[i,j])/deltaT)
  end

  for i=1:length(totalAvgSpectrum)
    totalAvgSpectrumFlat[i] = totalAvgSpectrum[i]/nFiles
  end
  if !onlyUseAverages
    for i=1:length(totalAvgSpectrum)
      totalAvgSubSpectrumFlat[i] = totalAvgSubSpectrum[i]/specCount
    end
  end

  HDF5.h5write(outfilepath, "AvgStickCps", stickcpsFlat')
  HDF5.h5write(outfilepath, "AvgStickCpsErr", stickcpsErrFlat')
  HDF5.h5write(outfilepath, "MassList", masslistMasses)
  println("Wrote $(length(masslistMasses)) masses to file.")
  HDF5.h5write(outfilepath, "MassListIntegrationBordersLow", massborders.lowMass)
  HDF5.h5write(outfilepath, "MassListIntegrationBordersHigh", massborders.highMass)

  HDF5.h5write(outfilepath, "MassListIdx", APiTOFFunctions.mass2timebin(masslistMasses, referenceMassScaleMode, referenceMassScaleParameters))
  HDF5.h5write(outfilepath, "MassListIntegrationBordersIdxLow", APiTOFFunctions.mass2timebin(massborders.lowMass, referenceMassScaleMode, referenceMassScaleParameters))
  HDF5.h5write(outfilepath, "MassListIntegrationBordersIdxHigh", APiTOFFunctions.mass2timebin(massborders.highMass, referenceMassScaleMode, referenceMassScaleParameters))
  HDF5.h5write(outfilepath, "MassCalibMode", referenceMassScaleMode)
  HDF5.h5write(outfilepath, "MassCalibParameters", referenceMassScaleParameters)

  HDF5.h5write(outfilepath, "AvgStickCpsTimes", timeFlat)
  HDF5.h5write(outfilepath, "MassAxis", referenceMassAxis)
  if onlyUseAverages
    HDF5.h5write(outfilepath, "AvgSpectrum", totalAvgSpectrumFlat)
  else
    HDF5.h5write(outfilepath, "AvgSpectrum", totalAvgSubSpectrumFlat)
  end

  if (createTotalAvg == true)
    #HDF5.h5write(outfilepath, "SumSpecs", sdata(interpolatedSpectra))
    HDF5.h5write(outfilepath, "SumSpecMax", sdata(totalMaxSpectrum))
    HDF5.h5write(outfilepath, "SumSpecMin", sdata(totalMinSpectrum))
  end
  HDF5.h5write(outfilepath, "ElementNames", masslistElements)
  HDF5.h5write(outfilepath, "ElementMasses", masslistElementsMasses)
  compFlat = Array{Int64}(length(masslistCompositions[1]),length(masslistCompositions))
  [compFlat[:,i]=masslistCompositions[i] for i=1:length(masslistCompositions)]
  HDF5.h5write(outfilepath, "ElementalCompositions", compFlat)
  print("Processing took: ")
  toc()
end