source: MML/trunk/mml/doc_html/mml/quaderrors.html @ 4

Last change on this file since 4 was 4, checked in by zhangj, 10 years ago

Initial import--MML version from SOLEIL@2013
File size: 13.2 KB
Line 
1<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"
2                "http://www.w3.org/TR/REC-html40/loose.dtd">
3<html>
4<head>
5  <title>Description of quaderrors</title>
6  <meta name="keywords" content="quaderrors">
7  <meta name="description" content="QUADERRORS - Computes the error bars for quadrupole center measurement">
8  <meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
9  <meta name="generator" content="m2html &copy; 2003 Guillaume Flandin">
10  <meta name="robots" content="index, follow">
11  <link type="text/css" rel="stylesheet" href="../m2html.css">
12</head>
13<body>
14<a name="_top"></a>
15<div><a href="../index.html">Home</a> &gt;  <a href="index.html">mml</a> &gt; quaderrors.m</div>
16
17<!--<table width="100%"><tr><td align="left"><a href="../index.html"><img alt="<" border="0" src="../left.png">&nbsp;Master index</a></td>
18<td align="right"><a href="index.html">Index for mml&nbsp;<img alt=">" border="0" src="../right.png"></a></td></tr></table>-->
19
20<h1>quaderrors
21</h1>
22
23<h2><a name="_name"></a>PURPOSE <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
24<div class="box"><strong>QUADERRORS - Computes the error bars for quadrupole center measurement</strong></div>
25
26<h2><a name="_synopsis"></a>SYNOPSIS <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
27<div class="box"><strong>function QMS = quaderrors(Input1, FigureHandle, MontiCarloFlag) </strong></div>
28
29<h2><a name="_description"></a>DESCRIPTION <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
30<div class="fragment"><pre class="comment">QUADERRORS - Computes the error bars for quadrupole center measurement
31  QMS = quaderrors(FileName, FigureHandle, MontiCarloFlag)
32
33 Written by Greg Portmann</pre></div>
34
35<!-- crossreference -->
36<h2><a name="_cross"></a>CROSS-REFERENCE INFORMATION <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
37This function calls:
38<ul style="list-style-image:url(../matlabicon.gif)">
39</ul>
40This function is called by:
41<ul style="list-style-image:url(../matlabicon.gif)">
42<li><a href="quadplot.html" class="code" title="function [QMS, WarningString] = quadplot(Input1, FigureHandle, sigmaBPM)">quadplot</a>  QUADPLOT - Plots quadrupole centering data</li><li><a href="quadplotphysics.html" class="code" title="function [QMS, WarningString] = quadplotphysics(Input1, FigureHandle, sigmaBPM)">quadplotphysics</a>      QUADPLOTPHYSICS - Plots quadrupole centering data in physics units</li></ul>
43<!-- crossreference -->
44
45
46<h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
47<div class="fragment"><pre>0001 <a name="_sub0" href="#_subfunctions" class="code">function QMS = quaderrors(Input1, FigureHandle, MontiCarloFlag)</a>
480002 <span class="comment">%QUADERRORS - Computes the error bars for quadrupole center measurement</span>
490003 <span class="comment">%  QMS = quaderrors(FileName, FigureHandle, MontiCarloFlag)</span>
500004 <span class="comment">%</span>
510005 <span class="comment">% Written by Greg Portmann</span>
520006
530007
540008 Buffer = .03;
550009 HeightBuffer = .08;
560010
570011     
580012 <span class="comment">% Don't do the monticarlo plot if being called from quadplot</span>
590013 MontiCarloFlag = 1;
600014 [ST,I] = dbstack;
610015 <span class="keyword">for</span> i = 1:length(ST)
620016     <span class="keyword">if</span> strcmpi(ST(i).file,<span class="string">'quadplot.m'</span>)
630017         MontiCarloFlag = 0;
640018     <span class="keyword">end</span>
650019 <span class="keyword">end</span>
660020
670021
680022 <span class="comment">% Inputs</span>
690023 <span class="keyword">try</span>
700024     WarningString = [];
710025     <span class="keyword">if</span> nargin == 0
720026         [FileName, PathName] = uigetfile(<span class="string">'*.mat'</span>, <span class="string">'Input Quadrupole Center file.'</span>);
730027         <span class="keyword">if</span> ~isstr(FileName)
740028             <span class="keyword">return</span>
750029         <span class="keyword">else</span>
760030             eval([<span class="string">'load '</span>,PathName,FileName]);
770031         <span class="keyword">end</span>
780032     <span class="keyword">else</span>
790033         <span class="keyword">if</span> isstr(Input1)
800034             FileName = Input1;
810035             eval([<span class="string">'load '</span>, FileName]);
820036         <span class="keyword">elseif</span> isempty(Input1)
830037             [FileName, PathName] = uigetfile(<span class="string">'*.mat'</span>, <span class="string">'Input Quadrupole Center file.'</span>);
840038             <span class="keyword">if</span> ~isstr(FileName)
850039                 <span class="keyword">return</span>
860040             <span class="keyword">else</span>
870041                 eval([<span class="string">'load '</span>,PathName,FileName]);
880042             <span class="keyword">end</span>
890043         <span class="keyword">else</span>
900044             QMS = Input1;
910045             FileName = [];
920046         <span class="keyword">end</span>
930047     <span class="keyword">end</span>
940048 <span class="keyword">catch</span>
950049     error(<span class="string">'Problem getting input data'</span>);
960050 <span class="keyword">end</span>
970051
980052 <span class="comment">%QMS = quadplot(QMS);</span>
990053
1000054 <span class="keyword">if</span> nargin &lt; 2
1010055     FigureHandle = gcf;
1020056 <span class="keyword">else</span>
1030057     <span class="keyword">if</span> all(FigureHandle ~= 0)
1040058         figure(FigureHandle);
1050059     <span class="keyword">end</span>
1060060 <span class="keyword">end</span>
1070061
1080062 <span class="keyword">if</span> all(FigureHandle ~= 0)
1090063     clf reset
1100064     <span class="keyword">if</span> MontiCarloFlag
1110065         AxesHandles(1) = subplot(3,1,1);
1120066         AxesHandles(2) = subplot(3,1,2);
1130067         AxesHandles(3) = subplot(3,1,3);
1140068     <span class="keyword">else</span>
1150069         AxesHandles(1) = subplot(2,1,1);
1160070         AxesHandles(2) = subplot(2,1,2);
1170071     <span class="keyword">end</span>
1180072
1190073     <span class="keyword">if</span> QMS.QuadPlane ~= 1
1200074         set(FigureHandle,<span class="string">'units'</span>,<span class="string">'normal'</span>,<span class="string">'position'</span>,[Buffer .25+Buffer .5-Buffer-.003 .75-1.2*Buffer-HeightBuffer]);
1210075         <span class="comment">%set(FigureHandle,'units','normal','position',[0+Buffer .09+Buffer .5-2*Buffer .9-2*Buffer-HeightBuffer]);</span>
1220076     <span class="keyword">else</span>
1230077         set(FigureHandle,<span class="string">'units'</span>,<span class="string">'normal'</span>,<span class="string">'position'</span>,[.5+.003 .25+Buffer .5-Buffer-.003 .75-1.2*Buffer-HeightBuffer]);
1240078         <span class="comment">%set(FigureHandle,'units','normal','position',[.5+Buffer .09+Buffer .5-2*Buffer .9-2*Buffer-HeightBuffer]);</span>
1250079     <span class="keyword">end</span>
1260080
1270081     <span class="comment">%subplot(3,2,1);</span>
1280082     axes(AxesHandles(1));
1290083     errorbar(1:length(QMS.FitParameters(:,1)), QMS.FitParameters(:,1), QMS.FitParametersStd(:,1),<span class="string">'b'</span>);
1300084     hold on;
1310085     errorbar(1:length(QMS.FitParameters(:,1)), QMS.FitParameters(:,1), QMS.FitParametersStd(:,1),<span class="string">'.r'</span>);
1320086     errorbar(QMS.GoodIndex, QMS.FitParameters(QMS.GoodIndex,1), QMS.FitParametersStd(QMS.GoodIndex,1),<span class="string">'.b'</span>);
1330087     hold off
1340088     title(<span class="string">'Linear Fit'</span>);
1350089     ylabel(<span class="string">'Y-Intercept'</span>);
1360090     <span class="comment">%xlabel('BPM Number');</span>
1370091     grid
1380092     xaxis([0 length(QMS.FitParameters(:,1))+1]);
1390093
1400094
1410095     <span class="comment">%subplot(3,1,2);</span>
1420096     axes(AxesHandles(2));
1430097     errorbar(1:length(QMS.FitParameters(:,1)), QMS.FitParameters(:,2), QMS.FitParametersStd(:,2),<span class="string">'b'</span>);
1440098     hold on
1450099     errorbar(1:length(QMS.FitParameters(:,1)), QMS.FitParameters(:,2), QMS.FitParametersStd(:,2),<span class="string">'.r'</span>);
1460100     errorbar(QMS.GoodIndex, QMS.FitParameters(QMS.GoodIndex,2), QMS.FitParametersStd(QMS.GoodIndex,2),<span class="string">'.b'</span>);
1470101     hold off
1480102     xlabel(<span class="string">'BPM Number'</span>);
1490103     ylabel(<span class="string">'Slope'</span>);
1500104     <span class="comment">%xlabel('BPM Number');</span>
1510105     grid
1520106     xaxis([0 length(QMS.FitParameters(:,1))+1]);
1530107 <span class="keyword">end</span>
1540108
1550109
1560110 <span class="keyword">if</span> MontiCarloFlag
1570111     <span class="comment">% Monte Carlo the mean and standard deviation of all the good offsets</span>
1580112     N = 20000;
1590113     QMS.OffsetMean = NaN*ones(length(QMS.FitParameters(:,2)),1);
1600114     QMS.OffsetStd  = NaN*ones(length(QMS.FitParameters(:,2)),1);
1610115
1620116     OffsetMean1 = QMS.OffsetMean;
1630117     OffsetStd1  = QMS.OffsetStd;
1640118     NormalRV1 = randn(N,1);
1650119     NormalRV2 = randn(N,1);
1660120     <span class="keyword">if</span> QMS.QuadraticFit == 1
1670121         NormalRV3 = randn(N,1);
1680122     <span class="keyword">end</span>
1690123     <span class="keyword">for</span> i = QMS.GoodIndex(:)'   <span class="comment">% 1:length(QMS.Offset)</span>
1700124         rvinter = QMS.FitParameters(i,1) + QMS.FitParametersStd(i,1)*NormalRV1;
1710125         rvslope = QMS.FitParameters(i,2) + QMS.FitParametersStd(i,2)*NormalRV2;
1720126         <span class="keyword">if</span> QMS.QuadraticFit == 1
1730127             rvquadradic = QMS.FitParameters(i,3) + QMS.FitParametersStd(i,3)*NormalRV3;
1740128             x = (-rvslope + sqrt(rvslope.^2 - 4*rvquadradic.*rvinter)) ./ (2*rvquadradic);
1750129         <span class="keyword">else</span>
1760130             x = -rvinter ./ rvslope;
1770131         <span class="keyword">end</span>
1780132         OffsetMean1(i) = mean(x);
1790133         OffsetStd1(i)  = std(x);
1800134     <span class="keyword">end</span>
1810135
1820136     OffsetMean2 = QMS.OffsetMean;
1830137     OffsetStd2  = QMS.OffsetStd;
1840138     NormalRV1 = randn(N,1);
1850139     NormalRV2 = randn(N,1);
1860140     <span class="keyword">if</span> QMS.QuadraticFit == 1
1870141         NormalRV3 = randn(N,1);
1880142     <span class="keyword">end</span>
1890143     <span class="keyword">for</span> i = QMS.GoodIndex(:)'   <span class="comment">% 1:length(QMS.Offset)</span>
1900144         rvinter = QMS.FitParameters(i,1) + QMS.FitParametersStd(i,1)*NormalRV1;
1910145         rvslope = QMS.FitParameters(i,2) + QMS.FitParametersStd(i,2)*NormalRV2;
1920146         <span class="keyword">if</span> QMS.QuadraticFit == 1
1930147             rvquadradic = QMS.FitParameters(i,3) + QMS.FitParametersStd(i,3)*NormalRV3;
1940148             x = (-rvslope + sqrt(rvslope.^2 - 4*rvquadradic.*rvinter)) ./ (2*rvquadradic);
1950149         <span class="keyword">else</span>
1960150             x = -rvinter ./ rvslope;
1970151         <span class="keyword">end</span>
1980152         OffsetMean2(i) = mean(x);
1990153         OffsetStd2(i)  = std(x);
2000154     <span class="keyword">end</span>
2010155
2020156     QMS.OffsetMean(:,1) = [OffsetMean1 + OffsetMean2] / 2;
2030157     QMS.OffsetStd(:,1)  = [OffsetStd1  + OffsetStd2]  / 2;
2040158
2050159
2060160     <span class="keyword">if</span> all(FigureHandle ~= 0)
2070161
2080162         <span class="comment">% plot([OffsetStd1  OffsetStd2  OffsetStd3],'x')</span>
2090163
2100164         <span class="comment">%subplot(3,1,3);</span>
2110165         axes(AxesHandles(3));
2120166         <span class="comment">%errorbar(1:length(QMS.Offset), QMS.Offset, QMS.OffsetStd,'.r');</span>
2130167         <span class="comment">%hold on</span>
2140168         errorbar(QMS.GoodIndex, QMS.Offset(QMS.GoodIndex), QMS.OffsetStd(QMS.GoodIndex),<span class="string">'.b'</span>);
2150169
2160170
2170171         <span class="comment">% errorbar(1:length(QMS.Offset), QMS.Offset, QMS.OffsetStd,'.b');</span>
2180172         <span class="comment">% hold on</span>
2190173         <span class="comment">% errorbar(1:length(QMS.Offset), QMS.Offset, OffsetStd1,'.r');</span>
2200174         <span class="comment">% errorbar(1:length(QMS.Offset), QMS.Offset, OffsetStd2,'.g');</span>
2210175         <span class="comment">% errorbar(1:length(QMS.Offset), QMS.Offset, OffsetStd3,'.k');</span>
2220176
2230177         <span class="comment">% errorbar(QMS.GoodIndex, QMS.Offset(QMS.GoodIndex), QMS.OffsetStd(QMS.GoodIndex),'.b');</span>
2240178         <span class="comment">% hold on</span>
2250179         <span class="comment">% errorbar(QMS.GoodIndex, QMS.Offset(QMS.GoodIndex), OffsetStd1(QMS.GoodIndex),'.r');</span>
2260180         <span class="comment">% errorbar(QMS.GoodIndex, QMS.Offset(QMS.GoodIndex), OffsetStd2(QMS.GoodIndex),'.g');</span>
2270181         <span class="comment">% errorbar(QMS.GoodIndex, QMS.Offset(QMS.GoodIndex), OffsetStd3(QMS.GoodIndex),'.k');</span>
2280182
2290183         hold off
2300184         ylabel(<span class="string">'X-Intercept = - (Y-Intercept)/Slope'</span>);
2310185         ylabel(<span class="string">'Monte Carlo'</span>);
2320186         xlabel(<span class="string">'BPM Number'</span>);
2330187         grid
2340188         xaxis([0 length(QMS.Offset)+1]);
2350189         <span class="comment">%yaxis([-.025 .025]+QMS.Center);</span>
2360190     <span class="keyword">end</span>
2370191
2380192     QMS.OffsetSTDMontiCarlo = sqrt((sum(QMS.OffsetStd(QMS.GoodIndex).^2))/length(QMS.GoodIndex));
2390193 <span class="keyword">end</span></pre></div>
240<hr><address>Generated on Mon 21-May-2007 15:29:18 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/">m2html</a></strong> &copy; 2003</address>
241</body>
242</html>
Note: See TracBrowser for help on using the repository browser.