function  [t,d,s,eventT]=PrepareData(iIFO,...
            ifoName,FrameType,ChannelName,filterState,...
            FilterPar,event,DataType,CalibrationFlag,...
	    sizeJobs,simulPar,MC,logsFid);
%=====================================================================
% This is a script to perform pre-process and rstat test over 
% end-of-pipeline burst events
% It is invoked by CrossCorrTriple.m, which defines the inputs
%
% function  [time,data,sample_rate]=PrepareData(iIFO,...
%            ifoName,FrameType,ChannelName,filterState,...
%            FilterPar,event,DataType,CalibrationFlag,...
%	     sizeJobs,simulPar,MC,logsFid);
% Required
%   ifoName
%   FrameType
%   ChannelName
%   filterState
%   event  = structure 
%            event.start 
%            event.duration 
%            event.del(iIFO)
%   sizeJobs  = number of events in the job - needed by the framecache 
%   simulPar  = structure
%             simulPar.doit
%             simulPar.wf
%             simulPar.id
%             simulPar.h0
%             simulPar.injTime (start of 1 sec of data with simul wf)
%   logsFid = id of log file
%
%  if DataType=0,10: Uses FRAMEPATH environment variable if defined
%     Otherwise uses hard-coded frame_dir
%
% INPUTS TO CHECK: hard-coded preprocess parameters
%
% L. Cadonati,Mar 22, 2004 
%      adapted from MainCorrC.m
% Keith Thorne, Feb 12, 2004
%      uses chanstruct/chanvector for LDR access
% $Id: PrepareData.m,v 1.1 2007/09/20 20:25:53 kathorne Exp $
%=====================================================================

Narg=13;
argOk = nargchk(Narg,Narg,nargin);
if(argOk)
    fprintf(1,' PrepareData: Bad Inputs\n - was %d should be %d\n',nargin,Narg);
    return
end

if(DataType==0) 
  DEFAULT_FRAME_DIR='/data2/S3frames/ASQ/<obs>/<epoch>';
  persistent frame_dir;
  if isempty(frame_dir)
    frame_dir = getenv('FRAMEPATH');
    if isempty(frame_dir)
        frame_dir = DEFAULT_FRAME_DIR;
        fprintf(' FRAMEPATH not defined - use %s\n',frame_dir);
    end    
  end
  
elseif(DataType==2 | DataType==4) 

DEFAULT_DATA_CACHE='/archive/home/bhughey/CorrHFS/cache/S5final.cache' ;
DEFAULT_MDC_CACHE='/archive/home/bhughey/CorrHFS/cache/S5final.cache' ;


  persistent data_cache;
  persistent mdc_cache;
  persistent datacache;
  persistent mdccache;

  if isempty(data_cache)
    data_cache = getenv('DATA_CACHE_FILE');
    if isempty(data_cache)
      data_cache = DEFAULT_DATA_CACHE;
      fprintf(logsFid,' DATA_CACHE_FILE not defined - use %s\n',data_cache);
    end    
  end
  if isempty(mdc_cache)
    mdc_cache = getenv('MDC_CACHE_FILE');
    if isempty(mdc_cache)
      mdc_cache = DEFAULT_MDC_CACHE;
      fprintf(logsFid,' MDC_CACHE_FILE not defined - use %s\n',mdc_cache);
    end    
  end

%=====================================================================
% The following lines of code were modified to never re-define
% data_cache and mdc_cache, since the way it was written appears to
% always load the default data_cache_file and mdc_cache_file despite
% having the corresponding environment variables properly defined.
%=====================================================================

  if isempty(data_cache)
    fprintf(1,'using /archive/home/bhughey/CorrHFS/cache/S5.txt as cachefile ');
    data_cache='/archive/home/bhughey/CorrHFS/cache/S5.txt';
  end
  if isempty(mdc_cache)
    mdc_cache='/archive/home/bhughey/CorrHFS/cache/S5.txt';
  end
  datacache=loadframecache(data_cache);
  mdccache=loadframecache(mdc_cache);

end


%{
  if isempty(datacache)
    fprintf(1,'using /archive/home/bhughey/CorrHFS/cache/S5.txt as cachefile ');
    data_cache='/archive/home/bhughey/CorrHFS/cache/S5.txt';
    datacache=loadframecache(data_cache);
  end
  mdc_cache;
  if isempty(mdccache)
    mdc_cache='/archive/home/bhughey/CorrHFS/cache/S5.txt';
    mdccache=loadframecache(mdc_cache);
  end

end
%}

%=====================================================================
% event start time and duration: 
% Start Offset seconds before the event start time, allowing for
% twice the duration of the event plus a additional Offset seconds 
% afterwards (probably not needed, but since filters are applied in 
% zerophase I would rather load more data than less than needed) 
%=====================================================================
Offset=ceil(FilterPar.train_lpef+FilterPar.offset_lpef+1);

gps_start=floor(event.start);
event_start_offset=event.start-gps_start;
Duration=2*Offset+ceil(event.duration)+1;
    
gps_start_load = gps_start-Offset;
gps_stop_load  = gps_start_load+Duration;
    
% apply time lag as specified in input file
gps_start_load1 = gps_start_load + event.del(iIFO);
gps_stop_load1  = gps_start_load1 + Duration;
eventT = Offset + event_start_offset;


fprintf(logsFid,...
   '-------------------------------------------------------------------------------\n');
fprintf(logsFid,'Loading %s with lag = %10.4f : %14.4f - %14.4f\n',...
    ChannelName,event.del(iIFO),gps_start_load1,gps_start_load1+Duration);

%=====================================================================
% load data, using the method associated to DataType:
%=====================================================================

if (DataType==0)
  % read data from local folder
  [d, t, s]=getframedata(frame_dir,ChannelName, ...
        floor(gps_start_load1-1),ceil(gps_stop_load1+1),ifoName(1),FrameType);  
  if (s==0)
    sample_rate=0; time=0; data=0;
    fprintf(logsFid,' Skipped Event: %14.4f %10.4f\n',event.start,event.duration);
    return;
  else
    indices=find(t>=gps_start_load1 & t<gps_stop_load1);
    t=t(indices);
    d=d(indices);  
  end                                                                                                                                 

elseif (DataType==1)
  % use LDR and Channel class
  ifoChan = chanstruct(ChannelName);
  ifoChan.type = FrameType;
    
  % get list of frame files from cache to cover event
  [gpsTimes,frameFiles,durList]=eventframecache(ifoChan,...
                                gps_start_load1,Duration,sizeJobs);
  [d,s,seriesOK] = chanvector(ifoChan,gps_start_load1,Duration,...
                   gpsTimes,frameFiles,durList);

  if (seriesOK ~= 0)
    sample_rate=0; time=0; data=0;
    s=0; t=0; d=0;
    fprintf(logsFid,' Skipped Event: %14.4f %10.4f\n',event.start,event.duration);
    return;
  end
  
  lenSeries = length(d);
  t = gps_start_load1 + ((0:(lenSeries-1))/s);
  
  if (simulPar.MDC==1)        
    mdcChan = chanstruct(char(simulPar.MDCchannelName(iIFO)));
    mdcChan.instrument = char(simulPar.MDCinstrument(iIFO));
    mdcChan.type = char(simulPar.MDCframeType(iIFO));

    % specify a default sample rate to ensure we get a 
    % valid MDC data vector even if some frames are missing
    % (Keith Thorne, April 29, 2004)
    mdcChan.rate = s;
       
    if (event.duration > 60) 
       % only load MDC for the event duration, w/out padding
       % (this is important when searching over the playground: frames
       % do not exist outside the interval being searched and the 
       % padding makes the job bomb).
       
       [gpsTimesMDC,frameFilesMDC,durListMDC]=eventframecache(mdcChan,...
                     gps_start_load1+Offset,ceil(event.duration),sizeJobs);

       [dMDC,sMDC,mdcOK] = chanvector(mdcChan,gps_start_load1+Offset, ...
             ceil(event.duration),gpsTimesMDC,frameFilesMDC,durListMDC);

      pad=zeros(1,Offset*sMDC);
      dMDC=[pad dMDC];
      pad2=zeros(1,length(d)-length(dMDC));
      dMDC=[dMDC pad2];
    else
       [gpsTimesMDC,frameFilesMDC,durListMDC]=eventframecache(mdcChan,...
                                gps_start_load1,Duration,sizeJobs);
       [dMDC,sMDC,mdcOK] = chanvector(mdcChan,gps_start_load1,Duration,...
                                 gpsTimesMDC,frameFilesMDC,durListMDC);
    end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%save /archive/home/mzanolin/CorrHFSrebuild_devel/debugging/dMDC_H1time_t0_1.mat dMDC;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%
%  Sept, 2007 Keith Thorne
%  - To handle MDC data that is h(t) but to use DARM_ERR data, we 
%   need to apply a h(t)-->DARM_ERR conversion
%
    if(simulPar.MDCcvrtFlag == true)
        dMdc = h2darm_cvrt(dMdc,ifoChan,gps_start_load1,sMdc);    
    end
%
    mdcDurSamples = length(dMDC);
    if ((length(dMDC) ~= length(d)) | (sMDC ~= s))
        fprintf(logsFid, ...
       '# samples %d rate %f of MDC does not match raw data - abort\n', ...
          mdcDurSamples,sMDC);
        sample_rate=0; time=0; data=0;
        s=0; t=0; d=0;
        return
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  fprintf(resultFid, ' MDCscaleFactor=%3.6f ', MDCscaleFactor);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    d = d + dMDC*simulPar.MDCscaleFactor;
    clear dMDC;
  end
  indices=find(t>=gps_start_load1 & t<gps_stop_load1);
  t=t(indices);
  d=d(indices);

elseif (DataType==2 | DataType==4)
  if(DataType==2)
    % use caches to find frames
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%savefile = '/archive/home/mzanolin/CorrHFSrebuild_devel/debugging/framedata.data';
%save(savefile, 'datacache', 'ChannelName', 'FrameType', 'gps_start_load1', 'gps_stop_load1');
fprintf(logsFid,'\nPrepareData.m making function call [d,s,t]=readframedata(%s,%s,%s,%d,%d)\n',data_cache,ChannelName,FrameType,floor(gps_start_load1-1),ceil(gps_stop_load1+1));
%data_cache
%iIFO
%datacache
%ChannelName
%gps_start_load1
%gps_stop_load1
%ifoName
%FrameType
%duration=gps_stop_load1 - gps_start_load1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

	[d, s, t] = readframedata(datacache, ChannelName, FrameType, ...
                               floor(gps_start_load1-1), ceil(gps_stop_load1+1) );

startLoadT = gps_start_load1;
gps_start_load;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%if (iIFO == 1)
%	preparedataH1d = d;
%	save /archive/home/mzanolin/CorrHFSrebuild_devel/debugging/noiseH1.mat preparedataH1d;
%elseif (iIFO == 2)
%	preparedataL1d = d;
%	save /archive/home/mzanolin/CorrHFSrebuild_devel/debugging/noiseL1.mat preparedataL1d;
%elseif (iIFO == 3)
%	preparedataV1d = d;
%	save /archive/home/mzanolin/CorrHFSrebuild_devel/debugging/noiseV1.mat preparedataV1d;
%end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


    if (s==0)
      sample_rate=0; time=0; data=0;
      fprintf(logsFid,' Skipped Event: %14.4f %10.4f\n',event.start,event.duration);
      return;
    end                                                                                                                                 
  else
    state=randn('state');
    s=filterState.rate;    

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%MC.rms(iIFO)
%Duration
%s
%state
%simulPar
%logsFid
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    [t,d,ds,state]=GaussNoiseMC(MC.rms(iIFO),Duration,...
                     s,state,simulPar,logsFid);
    t=t+gps_start_load1;
  end

  if (simulPar.MDC==1)        
    mdc_type = char(simulPar.MDCframeType(iIFO));
    mdc_channel = char(simulPar.MDCchannelName(iIFO));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%mdccache
%mdc_channel
%mdc_type
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if (event.duration > 60) 
      % only load MDC for the event duration, w/out padding
      % (this is important when searching over the playground: frames
      % do not exist outside the interval being searched and the 
      % padding makes the job bomb).

      [dMDC, sMDC, tMDC] = readframedata(mdccache, mdc_channel, mdc_type, ...
                               floor(gps_start_load1+Offset), ceil(event.duration));    
      % zero pad with the Offset
      Npad=find(t==tMDC(1))-1;
      if (tMDC(1) ~=t(Npad+1))
        ERROR(' MISMATCH IN MDC ADDITION TO DATA ');
      end
      pad=zeros(1,Npad);
      dMDC=[pad dMDC];
      pad2=zeros(1,length(d)-length(dMDC));
      dMDC=[dMDC pad2];
    else
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%if (iIFO==1)
%	save /archive/home/mzanolin/CorrHFSrebuild_devel/debugging/readframedata_input_1_mdccache.mat mdccache;
%	save /archive/home/mzanolin/CorrHFSrebuild_devel/debugging/readframedata_input_1_mdc_channel.mat mdc_channel;
%	save /archive/home/mzanolin/CorrHFSrebuild_devel/debugging/readframedata_input_1_mdc_type.mat mdc_type;
%	save /archive/home/mzanolin/CorrHFSrebuild_devel/debugging/readframedata_input_1_gpsstartload1.mat gps_start_load1;
%	save /archive/home/mzanolin/CorrHFSrebuild_devel/debugging/readframedata_input_1_gpsstopload1.mat gps_stop_load1;
%end
%if (iIFO==2)
%        save /archive/home/mzanolin/CorrHFSrebuild_devel/debugging/readframedata_input_2_mdccache.mat mdccache;
%        save /archive/home/mzanolin/CorrHFSrebuild_devel/debugging/readframedata_input_2_mdc_channel.mat mdc_channel;
%        save /archive/home/mzanolin/CorrHFSrebuild_devel/debugging/readframedata_input_2_mdc_type.mat mdc_type;
%        save /archive/home/mzanolin/CorrHFSrebuild_devel/debugging/readframedata_input_2_gpsstartload1.mat gps_start_load1;
%        save /archive/home/mzanolin/CorrHFSrebuild_devel/debugging/readframedata_input_2_gpsstopload1.mat gps_stop_load1;
%end
%if (iIFO==3)
%       save /archive/home/mzanolin/CorrHFSrebuild_devel/debugging/readframedata_input_3_mdccache.mat mdccache;
%        save /archive/home/mzanolin/CorrHFSrebuild_devel/debugging/readframedata_input_3_mdc_channel.mat mdc_channel;
%        save /archive/home/mzanolin/CorrHFSrebuild_devel/debugging/readframedata_input_3_mdc_type.mat mdc_type;
%        save /archive/home/mzanolin/CorrHFSrebuild_devel/debugging/readframedata_input_3_gpsstartload1.mat gps_start_load1;
%        save /archive/home/mzanolin/CorrHFSrebuild_devel/debugging/readframedata_input_3_gpsstopload1.mat gps_stop_load1;
%end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      [dMDC, sMDC, tMDC] = readframedata(mdccache, mdc_channel, mdc_type, ...
                               floor(gps_start_load1-1), ceil(gps_stop_load1+1));

startTime = int2str(floor(event.start));

%      if(length(mdc_channel) == 7)
%         if(mdc_channel == 'H1:GW-H')
%	   eval(['save MDC-H_' , startTime , '.mat dMDC;']);
%         elseif(mdc_channel == 'L1:GW-H')
%	   eval(['save MDC-L_' , startTime , '.mat dMDC;']);
%	 end
%      else  
%	eval(['save MDC-V_' , startTime , '.mat dMDC;']);
%      end
    end

    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~          
    mdcDurSamples = length(dMDC);
    if ((length(dMDC) ~= length(d)) | (sMDC ~= s))
      fprintf(logsFid, ...
       '# samples %d rate %f of MDC does not match raw data (%d %f)-abort\n',...
          mdcDurSamples,sMDC,length(d),s);
        sample_rate=0; time=[]; data=[];
        s=0; t=0; d=0;
        return
    end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%	if (iIFO == 1)
%		dMDC_H1=dMDC;
%		save /archive/home/mzanolin/CorrHFSrebuild_devel/debugging/dMDC_H1time_1.mat dMDC_H1;
%		clear dMDC_H1;
%	elseif (iIFO == 2)
%		dMDC_H2=dMDC;
%		save /archive/home/mzanolin/CorrHFSrebuild_devel/debugging/dMDC_H1time_2.mat dMDC_H2;
%		clear dMDC_H2;
%	elseif (iIFO == 3)
%		dMDC_H3=dMDC;
%		save /archive/home/mzanolin/CorrHFSrebuild_devel/debugging/dMDC_H1time_3.mat dMDC_H3;
%		clear dMDC_H3;
%	end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%if ( strcmp(mdc_type,'JW1_WNB'))
%	if (iIFO == 1)
%	noisefile = ['/archive/home/mzanolin/CorrHFSrebuild_devel/debugging/noisedata_' num2str(iIFO) '.mat' ];
%	injectionfile = ['/archive/home/mzanolin/CorrHFSrebuild_devel/debugging/wavedata_' num2str(iIFO) '.mat' ];
%		noiseH1 = d; 
%		waveH1 = dMDC;
%		save /archive/home/mzanolin/CorrHFSrebuild_devel/debugging/noisedata_H1.mat noiseH1;
%		save /archive/home/mzanolin/CorrHFSrebuild_devel/debugging/wavedata_H1.mat waveH1;
%	elseif (iIFO == 2)
%	noisefile = ['/archive/home/mzanolin/CorrHFSrebuild_devel/debugging/noisedata_' num2str(iIFO) '.mat' ];
%	injectionfile = ['/archive/home/mzanolin/CorrHFSrebuild_devel/debugging/wavedata_' num2str(iIFO) '.mat' ];
%	        noiseH2 = d;
%	        waveH2 = dMDC;
%		save /archive/home/mzanolin/CorrHFSrebuild_devel/debugging/noisedata_L1.mat noiseH2;
%		save /archive/home/mzanolin/CorrHFSrebuild_devel/debugging/waveedata_L1.mat waveH2;
%	elseif (iIFO == 3)
%	noisefile = ['/archive/home/mzanolin/CorrHFSrebuild_devel/debugging/noisedata_' num2str(iIFO) '.mat' ];
%	injectionfile = ['/archive/home/mzanolin/CorrHFSrebuild_devel/debugging/wavedata_' num2str(iIFO) '.mat' ];
%	        noiseH3 = d;
%	        waveH3 = dMDC;
%		save /archive/home/mzanolin/CorrHFSrebuild_devel/debugging/noisedata_V1.mat noiseH3;
%		save /archive/home/mzanolin/CorrHFSrebuild_devel/debugging/wavedata_V1.mat waveH3;
%	end
%end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    d = d + dMDC*simulPar.MDCscaleFactor;
 
%    if(length(mdc_channel) == 7)
%       if(mdc_channel == 'H1:GW-H')
%	 eval(['save data-H_' , startTime , '.mat d;']);
%       elseif(mdc_channel == 'L1:GW-H')
%	 eval(['save data-L_' , startTime , '.mat d;']);
%       end
%    else
%      eval(['save data-V_' , startTime , '.mat d;']);
%    end
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%if (strcmp(mdc_type,'JW1_WNB') )
%	if (iIFO == 1)
%		merge1 = d;
%		save /archive/home/mzanolin/CorrHFSrebuild_devel/debugging/totaldata_H1.mat merge1;
%		clear merge1;
%	elseif (iIFO == 2)
%		merge2 = d;
%		save /archive/home/mzanolin/CorrHFSrebuild_devel/debugging/totaldata_L1.mat merge2;
%		clear merge2;
%	elseif (iIFO == 3)
%		merge3 = d;
%	        save /archive/home/mzanolin/CorrHFSrebuild_devel/debugging/totaldata_V1.mat merge3;
%	        clear merge3;
%	end
%end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    clear dMDC;
  end
indices=find(t>=gps_start_load1 & t<gps_stop_load1);
  t=t(indices);
  d=d(indices);
  
else
  fprintf(logsFid,...
   'PREPAREDATA: I do not know what to do with DataType = %d \n',DataType);
  return
end
  fprintf(1,'%14.4f %10.4f\n',event.start,event.duration);  
  
%======================================================  
% set t=0 at event.start (was t(1)=gps_start_load)
%======================================================  
t = (t - t(1)) - (event.start - gps_start_load );

%======================================================  
% add simulated signal, if requested
%======================================================  
if (simulPar.doit==1),
  [tsig,sigh,sigADC]=WFsimul(gps_start, 1, ifoName, simulPar, ...
                     s, logsFid, CalibrationFlag, 0);
  ind =find(t<(simulPar.injTime-event.start)+1/s/2 & ...
            t>=(simulPar.injTime-event.start)-1/s/2);

  if (CalibrationFlag==0) 
    %%%%% NOTE: SIGNAL SHOULD BE BAND-PASSED!!!!!!
    %%%%% this should be taken care of later on, but I need to think about it.
    d(ind:ind+s-1)=d(ind:ind+s-1)+sigh; 
  else
    d(ind:ind+s-1)=d(ind:ind+s-1)+sigADC; 
  end 

end  
%======================================================  
% CALIBRATE DATA (AS_Q->h)
%======================================================  
if (DataType ~=5 & (CalibrationFlag==2 | CalibrationFlag==12) ) 
  [dcal, scal]=CalibF(gps_start,ifoName,d,t,s,0,CalibrationFlag,0);
  d=dcal;
  if (scal~=s),
    fprintf(1,' Skipped Event: %14.4f %10.4f\n',event.start,event.duration);
   return;
  end  
end    

%======================================================  

if (s==0),
    fprintf(1,' Skipped Event: %14.4f %10.4f\n',event.start,event.duration);
   return;
end                                                                                                                                 

return