% Clear up before we start.
close all; clear all
% Load hit data from Mark's file of (x, y) coordinate pairs
load('hitlist.txt');
% Check number of pulses
whos hitlist

%% Make plot of x and y hit predictions (separately) vs. pulse number
% (Remember there may be missing pulses, not necessarily consecutive in time).

figure;
pulseid=[1:120];
subplot(2,2,1,'align')
plot(pulseid, hitlist(:,1),'.')
xlabel('pulse #'); ylabel('x hits/mm');

subplot(2,2,2,'align')
plot(pulseid, hitlist(:,2),'.')
xlabel('pulse #'); ylabel('y hits/mm');
title('Megarun, sample # 1, 21 Feb 2008');

%% Identify ranges of hits per region in y - NB: all by hand based on the
% above plot of y vs. pulse number.

subplot(2,2,3,'align')
% Plot all pulses
plot(pulseid(:), hitlist(:,1),'.b')

% Next, superimpose the groups of pulses corresponding to nominally
% different y values in our scan, each coloured separately.
%hold on
% First for x positions
ihitgrp(1).pulses = 1:14;    ihitgrp(1).col='r';
ihitgrp(2).pulses = 17:32;   ihitgrp(2).col='g';
ihitgrp(3).pulses = 35:50;   ihitgrp(3).col='y';
ihitgrp(4).pulses = 54:68;   ihitgrp(4).col='m';
ihitgrp(5).pulses = 72:86;   ihitgrp(5).col='c';
ihitgrp(6).pulses = 90:102;  ihitgrp(6).col='k';
ihitgrp(7).pulses = 106:120; ihitgrp(7).col='g';
hold on
%%
% Not the most elegant way, but...
for ii=1:size(ihitgrp,2)
  plot(pulseid(ihitgrp(ii).pulses), hitlist(ihitgrp(ii).pulses,1),['.' ihitgrp(ii).col]);
end
xlabel('pulse #'); ylabel('x hits/mm');

%%
% Next, for y positions.
subplot(2,2,4,'align')
plot(pulseid(:), hitlist(:,2),'.b')
hold on

for ii=1:size(ihitgrp,2)
  plot(pulseid(ihitgrp(ii).pulses), hitlist(ihitgrp(ii).pulses,2),['.' ihitgrp(ii).col]);
end
xlabel('pulse #'); ylabel('y hits/mm');

% Save output.
set(gcf,'PaperPositionMode','auto')
print( gcf, '-djpeg', '-r96', 'megarun_fig1.jpeg')
print( gcf, '-dpdf', 'megarun_fig1.pdf')

%% Plot of x vs. y of sample
figure;
% All in one large figure to start with...

% We make use of the ellipse.m (from David Long, see
% http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=
% 289&objectType=file)
% Let's define ellipse semi-major axes (for x, y, assuming not rotated), in
% mm.
ra=0.02; rb=0.01; 
plot(hitlist(:,1), hitlist(:,2),'.b'); xlabel('x hits/mm'); ylabel('y hits/mm');

title('Megarun, sample # 1, 21 Feb 2008');
axis equal % Ensure same units on both axes
hold on;
for ii=1:size(ihitgrp,2)
  plot(hitlist(ihitgrp(ii).pulses,1), hitlist(ihitgrp(ii).pulses,2),['.' ihitgrp(ii).col]);
  % Useful on this magnification...?
%  ellipse(ra,rb,0.,hitlist(ihitgrp(ii).pulses,1),hitlist(ihitgrp(ii).pulses,2),'k');
  grid on
end
hold off;
% Save output.
set(gcf,'PaperPositionMode','auto')
print( gcf, '-djpeg', '-r 96','megarun_fig2.jpeg')
print( gcf, '-dpdf', 'megarun_fig2.pdf')
%%


% And then separated into RoIs
% Plot centroids and a "typical" (yeah, right) beam ellipse to guide the
% eye.
grid on
for ii=1:size(ihitgrp,2)
  %  subplot(3,3,ii,'align')
  figure
  plot(hitlist(ihitgrp(ii).pulses,1), hitlist(ihitgrp(ii).pulses,2),['.' ihitgrp(ii).col]);
  ellipse(ra,rb,0.,hitlist(ihitgrp(ii).pulses,1),hitlist(ihitgrp(ii).pulses,2),'k');
  xlabel('x hits/mm'); ylabel('y hits/mm'); axis equal;
  title(['Megarun, sample 1, 21 Feb 2008, Region of Interest ' int2str(ii)])
  set(gcf,'PaperPositionMode','auto')
  print( gcf, '-djpeg', '-r 96', ['megarun_fig3_RoI_' int2str(ii) '.jpg'])
  print( gcf, '-dpdf', ['megarun_fig3_RoI_' int2str(ii) '.pdf'])
end


% Conclusion: we need field of view of 300um in x, 100um in y for each
% RoI, i.e. x400 on microscope would be sufficient, see:
% https://twiki.cern.ch/twiki/bin/view/ILCBDSColl/DamageDataAnalysis
%


%% C