Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ePHENIXFieldRICH.m
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file ePHENIXFieldRICH.m
1 close all
2 clear all
3 
4 %% RICH input
5 
6 p_RICH_Ref = 10;
7 % n_index = 1.00137; % C4F10
8 % n_index = 1.00054; % CF4
9 n_index = 1.00062; % CF4
10 Theta_max = acos(1./n_index)*1000; % mrad
11 
12 
13 MaxLength = 350;
14 etas = [1.02:.05:1.3,1.4:.2:3, 3:.5:4]';
15 N_step = 1000;
16 
17 
18 %% open field map
19 
20 file1 = 'DrawFieldMap_hfield_sPHENIX.2d.root.data';
21 % file2 = 'DrawFieldMap_hfield_fsPHENIX.2d.root.data';
22 
23 fileID = fopen(file1);
24 data1 = fscanf(fileID,'%f%f%f%f\n', [4 Inf])';
25 fclose(fileID);
26 
27 %% resampling
28 
29 bc1 = griddata(data1(:,2),data1(:,1),data1(:,4),0,0,'natural');
30 data1_scale = 1.4/1.5;
31 
32 [Z,R]= meshgrid(-400:2:400,0:2:300);
33 % [Z,R]= meshgrid(-200:2:200,0:2:200);
34 
35 br1 = griddata(data1(:,2),data1(:,1),data1(:,3),Z,R,'natural') .* data1_scale;
36 bz1 = griddata(data1(:,2),data1(:,1),data1(:,4),Z,R,'natural') .* data1_scale;
37 b1 = sqrt(br1.^2 + bz1.^2);
38 
39 %% RICH trajectory sampling
40 
41 thetas = (2. .* atan( exp(-etas))) * ones(1, N_step);
42 
43 trajectory_step = ones(size(etas)) * linspace(0,MaxLength,N_step);
44 step = MaxLength ./ (N_step-1);
45 
46 trajectory_r = sin(thetas) .* trajectory_step;
47 trajectory_z = cos(thetas) .* trajectory_step;
48 
49 
50 %% Field mapping
51 
52 trajectory_br = griddata(data1(:,2),data1(:,1),data1(:,3),trajectory_z,trajectory_r,'natural') .* data1_scale;
53 trajectory_bz = griddata(data1(:,2),data1(:,1),data1(:,4),trajectory_z,trajectory_r,'natural') .* data1_scale;
54 
55 % bending B
56 trajectory_bt = trajectory_bz.*cos(thetas + pi ./ 2) + trajectory_br.*sin(thetas + pi ./ 2);
57 % longitudinal B
58 trajectory_bl = trajectory_bz.*cos(thetas) + trajectory_br.*sin(thetas);
59 
60 
61 
62 %% angular kicks
63 
64 %BdL in T*m
65 trajectory_bdL = 0.01 .* step.* cumsum(trajectory_bt,2);
66 trajectory_angle_kick = 0.3 .* trajectory_bdL ./ p_RICH_Ref;
67 
68 %% ring spread
69 
70 RICH_volumne_cut = IsInCurvedRich(trajectory_z, etas * ones(1, N_step));
71 
72 std_RICH = zeros(size(etas));
73 for trajectory = 1:size(etas)
74  std_RICH(trajectory) = std(trajectory_angle_kick(trajectory,RICH_volumne_cut(trajectory,:))) ;
75 end
76 
77 
78 %% ring plot
79 
80 ylim = 8.01e-2*Theta_max; % max mrad
81 font_size = 17;
82 
83 m2014 = matfile('BABAR_V11_StationCuts_Z5.0cm_300.0cm.mat');
84 
85 figure('name',sprintf('Curved_RICH_Dispersion'),'PaperPositionMode','auto', ...
86  'position',[100,0,800,600]) ;
87 
88 
89 
90 std_RICH_draw = std_RICH* 1e3; % rad -> mrad
91 
92 hold on;
93 
94 
95 
96 plot(m2014.etas_fine, interp1(m2014.etas, m2014.std_RICH_draw, m2014.etas_fine,'spline'),'-',...
97  'LineWidth',2, 'Color',[0,0,.5]);
98 
99 
100 plot(etas, std_RICH_draw,'-',...
101  'LineWidth',2, 'Color',[0,.5,.5]);
102 RICH_fit = polyfit(etas, std_RICH_draw,3);
103 std_RICH_fit = polyval(RICH_fit,etas);
104 % plot(etas, std_RICH_fit,'--');
105 disp(RICH_fit./sqrt(2));
106 
107 set(gca,'XLim',[.5,4.5],'FontSize',font_size);
108 set(gca,'YLim',[0, ylim],'FontSize',font_size);
109 % title(sprintf('RICH Ring RMS Dispersion for track of p = %.1f GeV/c due to Field Bending', p_RICH_Ref ),'FontSize',15);
110 ylabel('RICH Ring Dispersion, \Delta\phi (mrad)','FontSize',font_size);
111 xlabel(['\eta'],'FontSize',font_size);
112 
113 box on;
114 
115 legend('2014 Concept: arXiv:1402.1209','2018 Update: sPH-cQCD-2018-001',...
116  'Location','East');
117 
118 % grid on;
119 
120 % text(3,2,...
121 % 'PID error \delta R = \Delta\phi / \sqrt{2N_{\gamma}}(10 GeV/c)/p',...
122 % 'Interpreter','Latex'...
123 % ,'FontSize',font_size...
124 % );
125 text(2.8,2.3,...
126  'RICH ring error $\delta R = \Delta\phi / \sqrt{2N_{\gamma}} \times (10 GeV/c)/p$',...
127  'Interpreter','Latex'...
128  ,'FontSize',font_size...
129  ,'HorizontalAlignment','center'...
130  );
131 
132 
133 ax1 = gca;
134 
135 ax2 = axes('Position',get(ax1,'Position'),...
136  'XAxisLocation','bottom',...
137  'YAxisLocation','right',...
138  'Color','none',...
139  'LineWidth',1, ...
140  'YColor','b');
141 set(ax2,'YLim',[0,ylim/Theta_max*100],'FontSize',font_size);
142 set(get(ax2,'YLabel'),...
143  'String','RICH Ring Dispersion, \Delta\phi (percentage of \theta_{MAX})',...
144  'FontSize',font_size);
145 set(ax2,'XTick',[]);
146 
147 SaveCavas('ePHENIXFieldRICH');
148 
149 
150 
151 %% Geometry Check
152 
153 downsample = 4;
154 
155 
156 
157 figure('name',['DrawField'],'PaperPositionMode','auto', ...
158  'position',[100,0,2000,800]) ;
159 
160 % subplot(2,1,1);
161 
162 
163 contourf(Z,R,b1,0:.1:2.5,'ShowText','on')
164 daspect([1 1 1])
165 % colormap cool
166 hold on
167 quiver(Z(1:downsample:end,1:downsample:end),R(1:downsample:end,1:downsample:end),...
168  bz1(1:downsample:end,1:downsample:end),br1(1:downsample:end,1:downsample:end)...
169  )
170 
171 caxis([0, 2.5])
172 box on;
173 grid off;
174 
175 
176 h = colorbar;
177 set(get(h,'title'),'string','Field Strength [T]','FontSize',18);
178 xlabel('Z [cm]','FontSize',20);
179 ylabel('R [cm]','FontSize',20);
180 title('Magnetic field strength and vector in sPHENIX inner detector region','FontSize',20);
181 set(gca,'FontSize',18)
182 % text(0,80,'Central Tracking Volume','VerticalAlignment','bottom','HorizontalAlignment','center','FontSize',20,'Color',[.3 .3 .3]*.1)
183 % text(-110,80,'e-going Tracker','VerticalAlignment','bottom','HorizontalAlignment','right','FontSize',20,'Color',[.3 .3 .3]*.1)
184 % text(160,80,'h-going Tracker','VerticalAlignment','bottom','HorizontalAlignment','left','FontSize',20,'Color',[.3 .3 .3]*.1)
185 
186 
187 
188 h = rectangle('Position',[-102 0 +2*102 78-0]);
189 set(h,'EdgeColor',[.3 .3 .3]);
190 set(h,'LineWidth',4);
191 
192 % E-GEM
193 h = line([-1.1,-1.1].*100, [.1, 0.8].*100 );
194 set(h,'Color',[.3 .3 .3]);
195 set(h,'LineWidth',4,'LineStyle','-');
196 
197 h = line([-1.35,-1.35].*100, [.1, 0.8].*100 );
198 set(h,'Color',[.3 .3 .3]);
199 set(h,'LineWidth',4,'LineStyle','-');
200 
201 
202 % H-GEM
203 h = line([1.2,1.2].*100, [.1, 0.8].*100 );
204 set(h,'Color',[.3 .3 .3]);
205 set(h,'LineWidth',4,'LineStyle','-');
206 
207 h = line([1.5,1.5].*100, [.1, 1].*100 );
208 set(h,'Color',[.3 .3 .3]);
209 set(h,'LineWidth',4,'LineStyle','-');
210 
211 % h = line([2.7,2.7,2.5].*100, [.1, .8, 1.3].*100 );
212 % set(h,'Color',[.3 .3 .3]);
213 % set(h,'LineWidth',4,'LineStyle','-');
214 
215 for trajectory = 1:size(etas)
216  h = plot(trajectory_z(trajectory,:), trajectory_r(trajectory,:),'-','Color',[1,1,1]*.5 );
217  h = plot(trajectory_z(trajectory,RICH_volumne_cut(trajectory,:)), trajectory_r(trajectory,RICH_volumne_cut(trajectory,:)),...
218  '-','Color',[1,1,1]*.2,'LineWidth',2);
219 end
220 
221 SaveCavas('ePHENIXFieldRICH');
222 
223 
224