{ Промежуточный объект - виток масс-спектрометра

Содержит ДОПОЛНИТЕЛЬНЫЕ функции доступа, хранения данных и управления

}
{
Реперные массы:
  238,050784
  257,049187
  273,04073
  276,04759
  292,039133
  295,045993
  311,037536
  314,044396
  330,035939
  333,042799
  349,034342
  352,041202
}

unit MCAD_MI1201_Thread1;

interface
  uses
    Windows, Classes, SysUtils, SyncObjs, Registry,
    xSystem, MMTimer, MiscFunc,
    MCAD_MI1201_FileOutput, MCAD_MI1201_Range,
    c_ISSB, MassClbr, c_Mi1201, MITypes, MCAD_MI1201_Registry,
    MCAD_Event, McadUserTypes, MCAD_MI1201_TChartSeries,
    MCAD_MI1201_Thread_Types, MCAD_MI1201_Thread0;

Const
  RightScan = +1;
  LeftScan  = -1;

  DefaultPeakLevel = 0.5;
  DefaultPeakDelta = 2;


type
  // Сравнивающая функция
  tCompareFunc = Function (a,b:double):boolean;

  // Флаги и команды
  tPCFlags = (pcError,                 // Ошибка при работе процедуры
              pcBreak,                 // Прервать процедуру
              pcStarted,               // Процедура запущена
              pcFinished);             // Процедура завершена

  // Информация о пике
  tPeakData = record
    LeftLimit      :boolean;
    RightLimit     :boolean;  // Достигнуты левый или правый пределы

    IndexOfMaximum :longint;  // Индекс максимума пика
    IndexOfLeft    :longint;
    IndexOfRight   :longint;  // Индексы левого и правого предела на полувысоте

    LeftBound      :tMass;
    RightBound     :tMass;
    PeakCenter     :tMass;    // Найденая масса центра пика
    Height         :double;

    HalfWidth      :double;   // Полуширина пика
    end;


  tMsThread = class(MCAD_MI1201_Thread0.tMsThread)
    private
      prPeakData :tPeakData;
      prPeakDelta:double;  // Максимальная ширина пика
      prPeakLevel:double; // Часть от высоты на которой некоторый максимум

      function  GetLimitState:boolean;
      function  IndexCorrect(Var i:longint):boolean;

      function  GetIndexForY(Y:double; StartIndex,delta:longint; CompareTrue:tCompareFunc):longint;
      function  GetIndexForYAboveEqual(Y:double; StartIndex,delta:longint):longint;
      function  GetIndexForYBelowEqual(Y:double; StartIndex,delta:longint):longint;

      function  RangeTotal(LeftIndex,RightIndex:longint):double;

      function  DirectSignalGet(AMass:tMass; AChannel:tSeries):tChannelSignal;
      function  GetMassMinStep:tMass;

      procedure SetPeakDelta(PD:double);
      procedure SetPeakLevel(PL:double);

    protected

    public
      PCFlag    :tPCFlags;
      property  PeakData     :tPeakData read prPeakData;
      property  LimitOfRange :boolean   read GetLimitState;
      property  PeakDelta    :double    read prPeakDelta write SetPeakDelta;
      property  PeakLevel    :double    read prPeakLevel write SetPeakLevel;

      constructor Create(const Name:string);

      function  CenterOfPeak(StartMass:tMass; DoRemeasure:boolean):tMass;
  end;  // tMSThread 1


function Less(a,b:double):boolean;
function Great(a,b:double):boolean;


implementation
  Uses MCAD_MiscFuncs;


constructor tMSThread.Create;
begin
  Inherited Create(Name);
  FillChar(prPeakData, SizeOf(prPeakData), -1);
  prPeakDelta := DefaultPeakDelta;
  prPeakLevel := DefaultPeakLevel;
end;


procedure tMSThread.SetPeakDelta(PD:double);
begin

end;


procedure tMSThread.SetPeakLevel(PL:double);
begin
  if (PL > 0)and(PL < 1) then
    prPEakLevel := PL;
end;


// По значению Y ищет индекс, начиная со StartIndex
function Less(a,b:double):boolean;far;
begin
  Result := a < b;
end;


function Great(a,b:double):boolean;far;
begin
  Result := a > b;
end;


Function tMsThread.GetIndexForY(Y:double; StartIndex,delta:longint;
                                          CompareTrue:tCompareFunc):longint;
begin
  if not Assigned(DataSeries[MainChannel]) then
    begin
    Result := -1;
    exit;
    end;

  With DataSeries[MainChannel] do
    While IndexCorrect(StartIndex)and CompareTrue(Y,YValues[StartIndex]) do
      StartIndex := StartIndex + delta;

  Result := StartIndex;
end;


//
Function tMsThread.GetIndexForYAboveEqual(Y:double; StartIndex,delta:longint):longint;
begin
  Result := GetIndexForY(Y, StartIndex, Delta, Less);
end;


Function tMsThread.GetIndexForYBelowEqual(Y:double; StartIndex,delta:longint):longint;
begin
  Result := GetIndexForY(Y, StartIndex, Delta, Great);
end;


// Сумма значений в интервале индексов
function  tMSThread.RangeTotal(LeftIndex,RightIndex:longint):double;
Var
  i:longint;
begin
  with DataSeries[MainChannel] do
    begin
    Result := 0;
    if not IndexCorrect(LeftIndex) or not IndexCorrect(RightIndex) then
      Exit;

    if LeftIndex > RightIndex then
      begin
      i := LeftIndex;
      LeftIndex  := RightIndex;
      RightIndex := i;
      end;

    For i := LeftIndex to RightIndex do
      Result := Result + YValues[i];

    end; //with
end;


Function  tMSThread.GetLimitState;
begin
  with prPeakData do
    Result := LeftLimit or RightLimit;
end;


//Прверка правильности индекса и его корректировка к ближайшему
Function tMsThread.IndexCorrect;
begin
  with prPeakData, DataSeries[MainChannel] do
    begin
    LeftLimit    := (i < 0);
    RightLimit   := (i >= Count);

    if LeftLimit then
      i := 0;

    If RightLimit then
      i := Count - 1;

    Result := not LeftLimit and not RightLimit;
    end; //with
end;


function tMsThread.DirectSignalGet(AMass:tMass; AChannel:tSeries):tChannelSignal;
begin
  if WaitReady then
    begin
    try
      DataGetAndStoreEx(AMass,False);
      if IsNoError then
        begin
        Result.Time:=Signals.Time;
        Result.Mass:=Signals.Mass;
        if AChannel<>sIntegrationTime then
          begin
          Result.Signal:=Signals.Signals[tChannel(AChannel)];
          end
        else
          Result.Signal:=Signals.Time;
        end
      else
        begin
        Result.Mass:=0;
        end;
    except
      end;
    SetReady;
    end
  else
    begin
    Result.Mass:=0;
    end;
end;

function  tMsThread.GetMassMinStep:tMass;
  var c:integer; dm:tMass;
begin
  Result:=MassStep;
  if Result=0 then
    Exit;
  c:=prMsSp.Counter;
  dm:=prMsSp.Mass;
  if Result<0 then
    begin
    dm:=(prMsSp.Counter2Mass(Pred(c))-dm);
    if Result>dm then
      Result:=dm;
    end
  else
    begin
    dm:=(prMsSp.Counter2Mass(Succ(c))-dm);
    if Result<dm then
      Result:=dm;
    end;
end;

function tMSThread.CenterOfPeak(StartMass:tMass; DoRemeasure:boolean):tMass;

  procedure DeleteRange(m1,m2:tMass);
    // Удаляем  точки из диапазона m1..m2, не включая концы
    Var
      t:tMass;
      Cur{,Last}:longint;
    begin
      with DataSeries[MainChannel] do
        begin
        if Count = 0 then
          Exit;

        if m1 > m2 then
          begin
//          t  := m1;
          m1 := m2;
//          m2 := t;
          end;

        // Следующий индекс для минимальной массы
        Cur := FindNearestX(m1) + 1;
        m2  := XValues[FindNearestX(m1)];
        DataBeginWrite;

        while (Cur < Count)and(XValues[Cur] < m2) do
          DataSeries[MainChannel].Delete(Cur);

        DataEndWrite;
        end; // with
    end;


Var
  Sig:tChannelSignal;
  LeftScanned :boolean;
  RightScanned:boolean;     // Флаги просканированности сторон

  Direct      :longint;
  CurIndex,PredIndex:longint;
  CurMass     :tMass;
  PredMass    :tMass;
  MaxMass     :tMass;

  Height      :double;       // Высота, до которой нужно мерить


begin
  // Измерение запущено
  PCFlag := pcStarted;

  // Для начала заполняем всё нулями
  with prPeakData do
    begin
    FillChar(prPeakData, SizeOf(prPeakData), -1);
    Result := 0;
    end;

  if not assigned(DataSeries[MainChannel])or not WaitReady then
    begin
    PCFlag := pcError;
    Notify(evPeakCenterEnd);
    Exit;
    end;

  with prPeakData, DataSeries[MainChannel] do
    begin
    // Останавливаем измерение
    Stop;

    if DoRemeasure then
      begin
      Sig := DirectSignalGet(StartMass,MainChannel);  // Перемериваем стартовую массу
      StartMass := Sig.Mass;
      Height    := Sig.Signal;
      end
    else
      Height := YValues[FindNearestX(StartMass)];

    LeftScanned  := False;
    RightScanned := False;

    Direct         := LeftScan;              // Начинаем влево (например)
    IndexOfMaximum := FindNearestX(StartMass);      // Пока считаем это максимумом
    CurMass        := StartMass;

    Repeat
      // Проверка на прерывание
      if PCFlag = pcBreak then
        begin
        PCFlag := pcError;
        SetReady;
        Notify(evPeakCenterEnd);
        Exit;
        end;

      // Формируем очередную массу
      PredMass := MassAjust(CurMass);
      CurMass  := MassAjust(PredMass + Direct * GetMassMinStep);
      MaxMass  := XValues[IndexOfMaximum];

      if DoRemeasure then
        begin
        // Убираем все промежуточные точки...

        // Удаляем старые массы в пределах массового шага
        DeleteRange(PredMass, CurMass);

        // ...и перемериваем для нового индекса
        Sig      := DirectSignalGet(CurMass,MainChannel);
        CurMass  := Sig.Mass;
        IndexOfMaximum := FindNearestX(MaxMass);
      //  Mass     := PMass;
        end;

//      CurMass := MassAjust(CurMass);
      CurIndex  := FindNearestX(CurMass);
      PredIndex := FindNearestX(PredMass);
      if not DoRemeasure then
        Sig.Signal := YValues[CurIndex];

      // Проверяем на максимум
      if Sig.Signal > Height then // Пик здесь
        begin
        IndexOfMaximum := CurIndex;
        MaxMass        := CurMass;
        Height         := Sig.Signal;
        end;

      // Проверяем: достигнута ли граница?
      // Здесь LeftXIndex и RightXIndex либо будут заполнены,
      // либо возникнет ошибка границ пика
      if Sig.Signal <= Height * prPeakLevel then // Направление сосканировано
        begin
        if Direct * LeftScan  >= 0 then
          LeftScanned := True
        else
          RightScanned := True;

        Direct  := -Direct;  // Противоположное направление
        Sig     := DirectSignalGet(MaxMass,MainChannel); // Другое направление сканируем от максимума
        CurMass := Sig.Mass;
        end;

      // Проверка: не слишком ли мы увлеклись
      if abs(CurMass - MaxMass) > PeakDelta then
        begin
        Break;
        end;

    Until LeftScanned and RightScanned;

    // пересчитываем границы
    DirectSignalGet(PredMass,MainChannel); // сброс данных

    IndexOfLeft  := GetIndexForYAboveEqual(Height * prPeakLevel,
                                           IndexOfMaximum, LeftScan);
    IndexOfRight := GetIndexForYAboveEqual(Height * prPeakLevel,
                                           IndexOfMaximum, RightScan);

    LeftBound    := XValues[IndexOfLeft];
    RightBound   := XValues[IndexOfRight];

    HalfWidth  := abs(RightBound - LeftBound);

    PeakCenter := (RightBound + LeftBound) / 2;
    Result     := PeakCenter;
    if HalfWidth <> 0 then
      HalfWidth := PeakCenter / HalfWidth
    else
      HalfWidth := 0;
    MassX      := PeakCenter;
    Pause;
    end; //with
  SetReady;
  PCFlag := pcFinished;
  Notify(evPeakCenterEnd);
end;  // CenterOfPeak

END.