ForwardTracking/macros/ePHENIXFieldRICH.m

0001 close all
0002 clear all
0003
0004 %% RICH input
0005
0006 p_RICH_Ref = 10;
0007 % n_index = 1.00137; % C4F10
0008 % n_index = 1.00054; % CF4
0009 n_index = 1.00062; % CF4
0010 Theta_max = acos(1./n_index)*1000; % mrad
0011
0012
0013 MaxLength = 350;
0014 etas = [1.02:.05:1.3,1.4:.2:3, 3:.5:4]';
0015 N_step = 1000;
0016
0017
0018 %% open field map
0019
0020 file1 = 'DrawFieldMap_hfield_sPHENIX.2d.root.data';
0021 % file2 = 'DrawFieldMap_hfield_fsPHENIX.2d.root.data';
0022
0023 fileID = fopen(file1);
0024 data1 = fscanf(fileID,'%f%f%f%f\n', [4 Inf])';
0025 fclose(fileID);
0026
0027 %% resampling
0028
0029 bc1 = griddata(data1(:,2),data1(:,1),data1(:,4),0,0,'natural');
0030 data1_scale = 1.4/1.5;
0031
0032 [Z,R]= meshgrid(-400:2:400,0:2:300);
0033 % [Z,R]= meshgrid(-200:2:200,0:2:200);
0034
0035 br1 = griddata(data1(:,2),data1(:,1),data1(:,3),Z,R,'natural') .* data1_scale;
0036 bz1 = griddata(data1(:,2),data1(:,1),data1(:,4),Z,R,'natural') .* data1_scale;
0037 b1 = sqrt(br1.^2 + bz1.^2);
0038
0039 %% RICH trajectory sampling
0040
0041 thetas = (2. .* atan( exp(-etas))) * ones(1, N_step);
0042
0043 trajectory_step = ones(size(etas)) * linspace(0,MaxLength,N_step);
0044 step = MaxLength ./ (N_step-1);
0045
0046 trajectory_r = sin(thetas) .* trajectory_step;
0047 trajectory_z = cos(thetas) .* trajectory_step;
0048
0049
0050 %% Field mapping
0051
0052 trajectory_br = griddata(data1(:,2),data1(:,1),data1(:,3),trajectory_z,trajectory_r,'natural') .* data1_scale;
0053 trajectory_bz = griddata(data1(:,2),data1(:,1),data1(:,4),trajectory_z,trajectory_r,'natural') .* data1_scale;
0054
0055 % bending B
0056 trajectory_bt = trajectory_bz.*cos(thetas + pi ./ 2) + trajectory_br.*sin(thetas + pi ./ 2);
0057 % longitudinal B
0058 trajectory_bl = trajectory_bz.*cos(thetas) + trajectory_br.*sin(thetas);
0059
0060
0061
0062 %% angular kicks
0063
0064 %BdL in T*m
0065 trajectory_bdL = 0.01 .* step.* cumsum(trajectory_bt,2);
0066 trajectory_angle_kick = 0.3 .* trajectory_bdL ./ p_RICH_Ref;
0067
0068 %% ring spread
0069
0070 RICH_volumne_cut = IsInCurvedRich(trajectory_z, etas * ones(1, N_step));
0071
0072 std_RICH = zeros(size(etas));
0073 for trajectory = 1:size(etas)
0074     std_RICH(trajectory) = std(trajectory_angle_kick(trajectory,RICH_volumne_cut(trajectory,:)))  ;
0075 end
0076
0077
0078 %% ring plot
0079
0080 ylim =  8.01e-2*Theta_max; % max mrad
0081 font_size = 17;
0082
0083 m2014 = matfile('BABAR_V11_StationCuts_Z5.0cm_300.0cm.mat');
0084
0085 figure('name',sprintf('Curved_RICH_Dispersion'),'PaperPositionMode','auto', ...
0086     'position',[100,0,800,600]) ;
0087
0088
0089
0090 std_RICH_draw = std_RICH* 1e3; % rad -> mrad
0091
0092 hold on;
0093
0094
0095
0096 plot(m2014.etas_fine, interp1(m2014.etas, m2014.std_RICH_draw, m2014.etas_fine,'spline'),'-',...
0097     'LineWidth',2, 'Color',[0,0,.5]);
0098
0099
0100 plot(etas, std_RICH_draw,'-',...
0101     'LineWidth',2, 'Color',[0,.5,.5]);
0102 RICH_fit = polyfit(etas, std_RICH_draw,3);
0103 std_RICH_fit = polyval(RICH_fit,etas);
0104 % plot(etas, std_RICH_fit,'--');
0105 disp(RICH_fit./sqrt(2));
0106
0107 set(gca,'XLim',[.5,4.5],'FontSize',font_size);
0108 set(gca,'YLim',[0, ylim],'FontSize',font_size);
0109 % title(sprintf('RICH Ring RMS Dispersion for track of p = %.1f GeV/c due to Field Bending', p_RICH_Ref ),'FontSize',15);
0110 ylabel('RICH Ring Dispersion, \Delta\phi (mrad)','FontSize',font_size);
0111 xlabel(['\eta'],'FontSize',font_size);
0112
0113 box on;
0114
0115 legend('2014 Concept: arXiv:1402.1209','2018 Update: sPH-cQCD-2018-001',...
0116     'Location','East');
0117
0118 % grid on;
0119
0120 % text(3,2,...
0121 %     'PID error \delta R = \Delta\phi / \sqrt{2N_{\gamma}}(10 GeV/c)/p',...
0122 %     'Interpreter','Latex'...
0123 %     ,'FontSize',font_size...
0124 % );
0125 text(2.8,2.3,...
0126     'RICH ring error $\delta R = \Delta\phi / \sqrt{2N_{\gamma}} \times (10 GeV/c)/p$',...
0127     'Interpreter','Latex'...
0128     ,'FontSize',font_size...
0129     ,'HorizontalAlignment','center'...
0130     );
0131
0132
0133 ax1 = gca;
0134
0135 ax2 = axes('Position',get(ax1,'Position'),...
0136     'XAxisLocation','bottom',...
0137     'YAxisLocation','right',...
0138     'Color','none',...
0139     'LineWidth',1, ...
0140     'YColor','b');
0141 set(ax2,'YLim',[0,ylim/Theta_max*100],'FontSize',font_size);
0142 set(get(ax2,'YLabel'),...
0143     'String','RICH Ring Dispersion, \Delta\phi (percentage of \theta_{MAX})',...
0144     'FontSize',font_size);
0145 set(ax2,'XTick',[]);
0146
0147 SaveCavas('ePHENIXFieldRICH');
0148
0149
0150
0151 %% Geometry Check
0152
0153 downsample = 4;
0154
0155
0156
0157 figure('name',['DrawField'],'PaperPositionMode','auto', ...
0158     'position',[100,0,2000,800]) ;
0159
0160 % subplot(2,1,1);
0161
0162
0163 contourf(Z,R,b1,0:.1:2.5,'ShowText','on')
0164 daspect([1 1 1])
0165 % colormap cool
0166 hold on
0167 quiver(Z(1:downsample:end,1:downsample:end),R(1:downsample:end,1:downsample:end),...
0168     bz1(1:downsample:end,1:downsample:end),br1(1:downsample:end,1:downsample:end)...
0169     )
0170
0171 caxis([0, 2.5])
0172 box on;
0173 grid off;
0174
0175
0176 h = colorbar;
0177 set(get(h,'title'),'string','Field Strength [T]','FontSize',18);
0178 xlabel('Z [cm]','FontSize',20);
0179 ylabel('R [cm]','FontSize',20);
0180 title('Magnetic field strength and vector in sPHENIX inner detector region','FontSize',20);
0181 set(gca,'FontSize',18)
0182 % text(0,80,'Central Tracking Volume','VerticalAlignment','bottom','HorizontalAlignment','center','FontSize',20,'Color',[.3 .3 .3]*.1)
0183 % text(-110,80,'e-going Tracker','VerticalAlignment','bottom','HorizontalAlignment','right','FontSize',20,'Color',[.3 .3 .3]*.1)
0184 % text(160,80,'h-going Tracker','VerticalAlignment','bottom','HorizontalAlignment','left','FontSize',20,'Color',[.3 .3 .3]*.1)
0185
0186
0187
0188 h = rectangle('Position',[-102 0  +2*102  78-0]);
0189 set(h,'EdgeColor',[.3 .3 .3]);
0190 set(h,'LineWidth',4);
0191
0192 % E-GEM
0193 h = line([-1.1,-1.1].*100, [.1, 0.8].*100 );
0194 set(h,'Color',[.3 .3 .3]);
0195 set(h,'LineWidth',4,'LineStyle','-');
0196
0197 h = line([-1.35,-1.35].*100, [.1, 0.8].*100 );
0198 set(h,'Color',[.3 .3 .3]);
0199 set(h,'LineWidth',4,'LineStyle','-');
0200
0201
0202 % H-GEM
0203 h = line([1.2,1.2].*100, [.1, 0.8].*100 );
0204 set(h,'Color',[.3 .3 .3]);
0205 set(h,'LineWidth',4,'LineStyle','-');
0206
0207 h = line([1.5,1.5].*100, [.1, 1].*100 );
0208 set(h,'Color',[.3 .3 .3]);
0209 set(h,'LineWidth',4,'LineStyle','-');
0210
0211 % h = line([2.7,2.7,2.5].*100, [.1, .8, 1.3].*100 );
0212 % set(h,'Color',[.3 .3 .3]);
0213 % set(h,'LineWidth',4,'LineStyle','-');
0214
0215 for trajectory = 1:size(etas)
0216     h = plot(trajectory_z(trajectory,:), trajectory_r(trajectory,:),'-','Color',[1,1,1]*.5 );
0217     h = plot(trajectory_z(trajectory,RICH_volumne_cut(trajectory,:)), trajectory_r(trajectory,RICH_volumne_cut(trajectory,:)),...
0218         '-','Color',[1,1,1]*.2,'LineWidth',2);
0219 end
0220
0221 SaveCavas('ePHENIXFieldRICH');
0222
0223
0224