1/******************************************************************************* 2* 3* McStas, neutron ray-tracing package 4* Copyright (C) 1997-2011, All rights reserved 5* Risoe National Laboratory, Roskilde, Denmark 6* Institut Laue Langevin, Grenoble, France 7* 8* %I 9* Written by: Peter Willendrup, derived from TOF_lambda_monitor.comp 10* Date: May 23, 2012 11* Origin: DTU Physics 12* 13* Special "Brilliance" monitor. 14* 15* %D 16* If used in the right setting, will output "instantaneous" and "mean" brilliances in units of Neutrons/cm^2/ster/AA/s. Conditions for proper units: 17* <ul> 18* <li>Use a with a source of area 1x1cm 19* <li>The source must illuminate/focus to an area of 1x1cm a 1m distance 20* <li>Parametrise the Brilliance_monitor with the frequency of the source 21* <li>To not change the source TOF distribution, place the Brilliance monitor close to the source! 22* </ul> 23* 24* with a source of area 1x1cm illuminating/focusing to an area of 1x1cm a 1m distance, this monitor will output "instantaneous" and "mean" brilliances in units of Neutrons/cm^2/ster/AA/s 25* 26* Here is an example of the use of the component. Note how the mentioned Unit conditions are implemented in instrument code. 27* 28*COMPONENT Source = ESS_moderator_long( 29* l_low = lambdamin, l_high = lambdamax, dist = 1, xw = 0.01, yh = 0.01, 30* freq = 14, T=50, tau=287e-6, tau1=0, tau2=20e-6, 31* n=20, n2=5, d=0.00286, chi2=0.9, I0=6.9e11, I2=27.6e10, 32* branch1=0, branch2=0.5, twopulses=0, size=0.01) 33* AT (0, 0, 0) RELATIVE Origin 34* 35*COMPONENT BRIL = Brilliance_monitor(nlam=196,nt=401,filename="bril.sim", 36* t_0=0,t_1=4000,lambda_0=lambdamin, 37* lambda_1=lambdamax, Freq=14) 38*AT (0,0,0.000001) RELATIVE Source 39* 40* %P 41* INPUT PARAMETERS: 42* 43* nlam: [1] Number of bins in wavelength 44* nt: [1] Number of bins in TOF 45* t_0: [us] Minimum time 46* t_1: [us] Maximum time 47* lambda_0: [AA] Minimum wavelength detected 48* lambda_1: [AA] Maximum wavelength detected 49* filename: [string] Defines filenames for the detector images. Stored as:<br>Peak_<filename> and Mean_<filename> 50* restore_neutron: [1] If set, the monitor does not influence the neutron state 51* Freq: [Hz] Source frequency. Use freq=1 for reactor source 52* srcarea: Source area [cm^2] 53* tofcuts: [1] Flag to generate TOF-distributions as function of wavelength 54* toflambda: [1] Flag to generate TOF-lambda distribution output 55* source_dist: [m] Distance from source. Beware when transporting through neutron optics! 56* xwidth: [m] width of monitor 57* yheight: [m] height of monitor 58* nowritefile: [1] If set, monitor will skip writing to disk 59* 60* OUTPUT PARAMETERS: 61* 62* Div_N: [] Array of neutron counts 63* Div_p: [] Array of neutron weight counts 64* Div_p2: [] Array of second moments 65* 66* %E 67*******************************************************************************/ 68DEFINE COMPONENT Brilliance_monitor 69DEFINITION PARAMETERS (nlam=101, nt=1001, srcarea=1) 70SETTING PARAMETERS (lambda_0=0, lambda_1=20, restore_neutron=0, Freq, int tofcuts=0, int toflambda=0, xwidth = 0.01, yheight=0.01, source_dist=1, string filename, t_0=0, t_1=20000, int nowritefile=0) 71OUTPUT PARAMETERS (tt_0, tt_1, BRIL_N, BRIL_p, BRIL_p2, BRIL_mean, BRIL_meanN, BRIL_meanE, BRIL_peak, BRIL_peakN, BRIL_peakE, BRIL_shape, BRIL_shapeN, BRIL_shapeE, xmin, xmax, ymin, ymax, ster, prsec, dlam, dt) 72//STATE PARAMETERS (x,y,z,vx,vy,vz,t,s1,s2,p) 73//POLARISATION PARAMETERS (sx,sy,sz) 74 75DECLARE 76%{ 77double BRIL_N[nt][nlam]; 78double BRIL_p[nt][nlam]; 79double BRIL_p2[nt][nlam]; 80double BRIL_mean[nlam]; 81double BRIL_meanN[nlam]; 82double BRIL_meanE[nlam]; 83double BRIL_peak[nlam]; 84double BRIL_peakN[nlam]; 85double BRIL_peakE[nlam]; 86 87double BRIL_shape[nt]; 88double BRIL_shapeN[nt]; 89double BRIL_shapeE[nt]; 90 91double tt_0, tt_1; 92double xmin, xmax, ymin, ymax; 93double ster; 94double prsec; 95double dlam; 96double dt; 97%} 98 99INITIALIZE 100%{ 101prsec=1e-6; 102 103int i,j; 104 105tt_0 = t_0*prsec; 106tt_1 = t_1*prsec; 107dt=(t_1-t_0)*prsec/nt; 108 dlam=(lambda_1-lambda_0)/(nlam-1); 109 for (i=0; i<nlam; i++) 110 { 111 BRIL_mean[i] = 0; 112 BRIL_peak[i] = 0; 113 BRIL_meanN[i] = 0; 114 BRIL_peakN[i] = 0; 115 BRIL_meanE[i] = 0; 116 BRIL_peakE[i] = 0; 117 for (j=0; j<nt; j++) 118 { 119 BRIL_N[j][i] = 0; 120 BRIL_p[j][i] = 0; 121 BRIL_p2[j][i] = 0; 122 if (i==0) { 123 BRIL_shape[j] = 0; 124 BRIL_shapeN[j] = 0; 125 BRIL_shapeE[j] = 0; 126 } 127 } 128 } 129 130 xmax = xwidth/2.0; 131 ymax = yheight/2.0; 132 xmin = -xmax; 133 ymin = -ymax; 134 135 ster = xwidth * yheight/(source_dist*source_dist); 136 137 %} 138 139TRACE 140 %{ 141 int i,j; 142 double div; 143 double lambda; 144 double Pnorm; 145 146 PROP_Z0; 147 lambda = (2*PI/V2K)/sqrt(vx*vx + vy*vy + vz*vz); 148 if (x>xmin && x<xmax && y>ymin && y<ymax && 149 lambda > lambda_0 && lambda < lambda_1) 150 { 151 if (t < tt_1 && t > tt_0) 152 { 153 i = floor((lambda - lambda_0)*nlam/(lambda_1 - lambda_0)); 154 j = floor((t-tt_0)*nt/(tt_1-tt_0)); 155 Pnorm=p/dlam/ster/srcarea; 156 BRIL_meanN[i]++; 157 BRIL_mean[i] += Pnorm; 158 BRIL_meanE[i] += Pnorm*Pnorm; 159 Pnorm=Pnorm/Freq/dt; 160 BRIL_N[j][i]++; 161 BRIL_p[j][i] += Pnorm; 162 BRIL_p2[j][i] += Pnorm*Pnorm; 163 164 } 165 } 166 if (restore_neutron) { 167 RESTORE_NEUTRON(INDEX_CURRENT_COMP, x, y, z, vx, vy, vz, t, sx, sy, sz, p); 168 } 169 %} 170 171SAVE 172%{ 173 if (!nowritefile) { 174 /* First, dump the 2D monitor */ 175 176 /* For each Wavelength channel, find peak brilliance */ 177 int i,j,jmax; 178 double Pnorm; 179 char ff[256]; 180 char tt[256]; 181 182 for (i=0; i<nlam; i++) 183 { 184 Pnorm = -1; 185 jmax = -1; 186 for (j=0; j<nt; j++) 187 { 188 if (BRIL_p[j][i]>=Pnorm) 189 { 190 Pnorm = BRIL_p[j][i]; 191 jmax=j; 192 } 193 BRIL_shape[j] = BRIL_p[j][i]; 194 BRIL_shapeN[j] = BRIL_N[j][i]; 195 BRIL_shapeE[j] = BRIL_p2[j][i]; 196 197 } 198 if (tofcuts == 1) { 199 sprintf(ff, "Shape_%s_%g",filename,lambda_0+i*dlam); 200 sprintf(tt, "Peak shape at %g AA",lambda_0+i*dlam); 201 DETECTOR_OUT_1D( 202 tt, 203 "TOF [us]", 204 "Peak Brilliance", 205 "Shape", t_0, t_1, nt, 206 &BRIL_shapeN[0],&BRIL_shape[0],&BRIL_shapeE[0], 207 ff); 208 } 209 BRIL_peakN[i] = BRIL_N[jmax][i]; 210 BRIL_peak[i] = BRIL_p[jmax][i]; 211 BRIL_peakE[i] = BRIL_p2[jmax][i]; 212 } 213 sprintf(ff, "Mean_%s",filename); 214 DETECTOR_OUT_1D( 215 "Mean brilliance", 216 "Wavelength [AA]", 217 "Mean Brilliance", 218 "Mean", lambda_0, lambda_1, nlam, 219 &BRIL_meanN[0],&BRIL_mean[0],&BRIL_meanE[0], 220 ff); 221 sprintf(ff, "Peak_%s",filename); 222 DETECTOR_OUT_1D( 223 "Peak brilliance", 224 "Wavelength [AA]", 225 "Peak Brilliance", 226 "Peak", lambda_0, lambda_1, nlam, 227 &BRIL_peakN[0],&BRIL_peak[0],&BRIL_peakE[0], 228 ff); 229 230 /* MPI related NOTE: Order is important here! The 2D-data used to generate wavelength-slices and calculate 231 the peak brilliance should be done LAST, otherwise we will get a factor of MPI_node_count too much as 232 scatter/gather has been performed on the arrays... */ 233 if (toflambda == 1) { 234 sprintf(ff, "TOFL_%s",filename); 235 DETECTOR_OUT_2D( 236 "TOF-wavelength brilliance", 237 "Time-of-flight [\\gms]", "Wavelength [AA]", 238 t_0, t_1, lambda_0, lambda_1, 239 nt, nlam, 240 &BRIL_N[0][0],&BRIL_p[0][0],&BRIL_p2[0][0], 241 filename); 242 } 243 } 244%} 245 246MCDISPLAY 247%{ 248 249 multiline(5, (double)xmin, (double)ymin, 0.0, 250 (double)xmax, (double)ymin, 0.0, 251 (double)xmax, (double)ymax, 0.0, 252 (double)xmin, (double)ymax, 0.0, 253 (double)xmin, (double)ymin, 0.0); 254%} 255 256END 257