1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
|
/*******************************************************************************
*
* McStas, neutron ray-tracing package
* Copyright (C) 1997-2009, All rights reserved
* Risoe National Laboratory, Roskilde, Denmark
* Institut Laue Langevin, Grenoble, France
*
* Component: StatisticalChopper_Monitor
*
* %Identification
* Written by: C. Monzat/E. Farhi
* Date: 2009
* Origin: ILL
*
* Monitor designed to compute the autocorrelation signal for the Statistical Chopper
*
* %Description
* This component is a time sensitive monitor which calculates the cross
* correlation between the pseudo random sequence of a statistical chopper
* and the signal received. It mainly uses fonctions of component Monitor_nD.
* It is possible to use the various options of the Monitor_nD but the user MUST NOT
* specify "time" in the options. Auto detection of the time limits is possible if the
* user chooses tmin=>tmax.
*
* StatisticalChopper_Monitor(options ="banana bins=500, abs theta limits=[5,105],bins=1000")
*
* %Parameters
* INPUT PARAMETERS:
* comp: [StatisticalChopper] quoted name of the component monitored
* tmin: [s] minimal time of detection
* tmax: [s] maximal time of detection
*
* CALCULATED PARAMETERS:
* delta_t: [s] interval of time of the detection
* T_p: [s] period of the statistical chopper comp
*
* %L
* R. Von Jan and R. Scherm. The statistical chopper for neutron time-of-flight spectroscopy. Nuclear Instruments and Methods, 80 (1970) 69-76.
*
* %End
*******************************************************************************/
DEFINE COMPONENT StatisticalChopper_Monitor COPY Monitor_nD
SETTING PARAMETERS (string comp, tmin=0.0, tmax=0.0)
/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */
DECLARE COPY Monitor_nD
INITIALIZE COPY Monitor_nD EXTEND
%{
char op[CHAR_BUF_LENGTH];
/* check if time is already in the variables to monitor */
if (!strstr(Vars.option, "time")) {
sprintf(op, "time %s", Vars.option);
strcpy(Vars.option,op);
} else {
/* we must make sure time is the first variable */
if (Vars.Coord_Number < 1 || strcmp(Vars.Coord_Var[1], "t"))
exit(fprintf(stderr, "StatisticalChopper_Monitor: %s: First variable must be 'time', currently '%s (%s)'\n.", NAME_CURRENT_COMP, Vars.Coord_Label[1], Vars.Coord_Var[1]));
if (strcmp(Vars.Coord_Var[0], "I") && strcmp(Vars.Coord_Var[0], "p"))
exit(fprintf(stderr, "StatisticalChopper_Monitor: %s: Must record signal intensity 'Intensity', currently '%s (%s)'\n.", NAME_CURRENT_COMP, Vars.Coord_Label[0], Vars.Coord_Var[0]));
}
/* set up Monitor_nD time limits */
if (tmin<tmax){
sprintf(op,"limits=[%g %g] %s",tmin, tmax, Vars.option);
}else{
sprintf(op,"auto %s",Vars.option);
}
strcpy(Vars.option,op);
/* re-initialize Monitor_nD */
Monitor_nD_Init(&DEFS, &Vars, xwidth, yheight,zdepth,xmin,xmax,ymin,ymax,zmin,zmax,0);
%}
TRACE COPY Monitor_nD
SAVE
%{
// avoid name clash with Detector structure members and component parameters.
#undef min
#undef xmin
#undef xmax
#undef ymin
#undef ymax
/* the detector file written by the Monitor_nD is stored in a 'MCDETECTOR detector' structure */
/* get back information from the StatisticalChopper */
int m = *COMP_GETPAR3(StatisticalChopper, comp, m); /* number of appertures in the sequence */
int nslit = *COMP_GETPAR3(StatisticalChopper, comp, nslit); /* length of the sequence (number of possible slits around disk) */
double nu = *COMP_GETPAR3(StatisticalChopper, comp, nu);
int *Sequence= *COMP_GETPAR3(StatisticalChopper, comp, Sequence);
double c = (double)(m-1)/(double)(nslit-1); /* duty cycle */
int f; /* number of periods in the detected time range */
/* save results, but do not free pointers */
detector = Monitor_nD_Save(&DEFS, &Vars);
f = (int)floor((detector.xmax-detector.xmin)*nu);
if (f < 1) f=1;
if (detector.intensity && m >= 1 && f >= 1 && nslit > 1 && detector.m > 1 && c < 1 && nu > 0 && Sequence) {
double S_sum = 0;
int i,j,k; /* indices for loops */
/* copy raw detector as correlation base */
long correlation_m = detector.m/f; /* new time binning for autocorrelation */
double *p0=malloc(correlation_m*detector.n*detector.p*sizeof(double)); /* Arrays to store correlation monitor */
double *p1=malloc(correlation_m*detector.n*detector.p*sizeof(double));
double *p2=malloc(correlation_m*detector.n*detector.p*sizeof(double));
/* initialize arrays to zero */
for (i=0;i<correlation_m*detector.n*detector.p;i++) p0[i]=p1[i]=p2[i]=0;
/* indices in loops for detector */
/* p1[i][j] = p1[index] with index= (!detector.istransposed ? i*n*p + j : i+j*m); */
/* compute scattering function S: Von Jan and Scherm, Eq (18) */
for (i=0;i<correlation_m;i++){
for (k=0;k<detector.n*detector.p;k++){
long index_new= !detector.istransposed ? i*detector.n*detector.p + k : i+k*correlation_m; /* [i][k] index in correlation */
for (j=0;j<f*nslit;j++){
long i2 = ( j*(detector.m-1)/(f*nslit-1) + i*(detector.m-1)/(correlation_m-1) ) % detector.m; /* time index in raw detector */
long index_old= !detector.istransposed ? i2*detector.n*detector.p + k : i2+k*detector.m; /* full [i][k] index in raw detector */
p1[index_new]+=(Sequence[j % nslit]-c) * detector.p1[ index_old ];
}
p1[index_new] /= m*f*(1-c);
p1[index_new] += -detector.min/m;
p0[index_new] = detector.events/correlation_m;
S_sum += p1[index_new];
}
}
/* Normalization of p1 compared to detector.p1: must have same integrated intensity */
/* normalizations:
coeff_ZS=(S_sum/detector.intensity)*f;
coeff_ZS=(S_sum/(m*detector.intensity-detector.min/nu))*(m*f*(1-m/nslit)+m);
*/
double coeff_ZS=(S_sum/(m*detector.intensity-detector.min/nu))*(m*f*(1.0-(double)m/nslit)+m);
for (i=0;i<correlation_m*detector.n*detector.p;i++) {
p1[i] /= coeff_ZS;
}
/* Calculation of delta(S)^2: Von Jan and Scherm, Eq (19) */
for (i=0;i<correlation_m;i++){
for (k=0;k<detector.n*detector.p;k++){
long index_new= !detector.istransposed ? i*detector.n*detector.p + k : i+k*correlation_m; /* [i][k] index in correlation */
for (j=0;j<f*nslit;j++){
long i2 = ( j*(detector.m-1)/(f*nslit-1) + i*(detector.m-1)/(correlation_m-1) ) % detector.m; /* time index in raw detector */
long index_old= !detector.istransposed ? i2*detector.n*detector.p + k : i2+k*detector.m; /* full [i][k] index in raw detector */
p2[index_new]+= (Sequence[j % nslit]-c)*(Sequence[j % nslit]-c)
* detector.p1[index_old]*detector.p1[index_old];
}
}
}
/* Transformation from sigma to p2 before output */
for (i=0;i<correlation_m*detector.n*detector.p;i++){
p2[i] /= (double)(m*m*f*f*(1-c)*(1-c));
if (p2[i] > 0) p2[i] = (p0[i] > 1 ? ((p0[i]-1)*p2[i]*p2[i] + p1[i]*p1[i]/p0[i])/p0[i]
: p1[i]);
}
char file[CHAR_BUF_LENGTH];
sprintf(file, "%s_corr", filename ? filename : NAME_CURRENT_COMP);
if (detector.rank == 1)
DETECTOR_OUT_1D("Correlation monitor",
detector.xlabel,detector.ylabel,
detector.xvar,detector.xmin,detector.xmax,correlation_m,
p0,p1,p2,file);
else if (detector.rank == 2)
DETECTOR_OUT_2D("Correlation monitor",
detector.xlabel,detector.ylabel,
detector.xmin,detector.xmax,detector.ymin,detector.ymax,
correlation_m,detector.n,
p0,p1,p2,file);
free(p0);free(p1);free(p2);
}
%}
FINALLY COPY Monitor_nD
MCDISPLAY COPY Monitor_nD
END
|
