RPC_sim_7.f
Go to the documentation of this file.
1  program rpc_sim
2 
3 * Author J.Repond. Version of April 21, 2010
4 * added transfer of interaction layer information to atuple
5 
6  implicit none
7 
8  integer i1,i2,i3,i4,i5,i6
9  integer i_event,n_event,n_point
10  integer i_pad_x,i_pad_y
11  integer hit(0:350,0:350,0:150)
12  integer z_store(0:10000),z_tmp(10000),f_tmp(10000),n_tmp
13  integer i_int,n_int,inter_layer,inter_layer_new
14 
15  real r1,r2,r3,r4,r5
16  real e,x,y,z,q
17  real t,q0
18  real par(3)
19  real x_store(0:10000),y_store(0:10000)
20  real x_tmp(10000),y_tmp(10000),dist_cut,e_cut
21  real x_int,y_int,radius_max
22  real q_int,q_tot,sigma1,sigma2,ratio
23  real q_integral1,q_integral2,q_integral,q_check
24  real q_pad(0:600,0:600,0:600)
25  real x_deposit,y_deposit
26  real x_crack(2),x_scaler,y_crack(4),y_scaler
27  real dist(10000)
28 
29  character*12 file0
30  character*3 file1
31  character*7 file2
32  character*4 file3
33  character*25 file4
34 
35 * data par /0.00057,1.735,0.0012/ ! 6.2 kV
36 * data par /0.00007,1.951,0.0010/ ! 6.3 kV data par /0.00001,2.167,0.0008/ ! 6.4 kV
37  data par /0.00001,2.167,0.0008/ ! 6.4 kV
38 
39 
40 * units are cm
41 
42  data n_int /1000/ ! Number of points in distributed charge
43  data radius_max /4.0/ ! Maximum extent of charge, default = 4.0
44  data y_crack /0.005,32.525,65.065,97.595/ ! Cracks in y
45 
46 * the famous 6 parameters
47 
48  data sigma1 /0.078/ ! Width of narrow Gaussian
49  data sigma2 /0.28/ ! Width of wide Gaussian
50  data ratio /0.310/ ! Ratio of contribution from slope1 and slope2; default = 0.345
51  data t /0.27/ ! Threshold (in pC); default = 0.3645
52  data dist_cut /0.092/ ! Minimum distance between hits in x,y plane; default = 0.092
53  data e_cut /0.15/ ! Distance at which charge not attenuated anymore, value needs to be > dist_cut
54  data q0 /+0.125/ ! offset in charge; default = +0.201
55  data file2 /'.atuple'/
56 
57 * MC output files
58 
59 * data file0 /'../mc1/Muon.'/
60  data file0 /'../mc1/Posi.'/
61 * data file0 /'../mc1/Pion.'/
62  data file3 /'.txt'/
63 * data file4 /'/HEPdata/ILC/repond/Posi.'/
64  data file4 /'/HEPdata/ILC/repond/Muon.'/
65 * data file4 /'/HEPdata/ILC/repond/Test.'/
66 
67 * check distance cuts
68 
69  if(dist_cut.gt.e_cut) then
70  write(6,*) 'Problem with distance cut values'
71  stop
72  endif
73 
74 * setup
75 
76  i_event=0
77  n_event=0
78 
79  write(6,*) 'Type run number'
80  read(5,*) file1
81 * open(unit=61,status='OLD',file=file0//file1//file3)
82  open(unit=61,status='OLD',file=file4//file1//file3) ! for Posi.007.txt and beyond
83  open(unit=62,status='UNKNOWN',file=file4//file1//file2) ! the output atuple
84 
85  read(61,*) z_store(0),x_store(0),y_store(0)
86  inter_layer=x_store(0)
87 
88 * calculate integral over charge distribution
89 
90  q_integral1=sigma1**2
91  q_integral2=sigma2**2
92 
93 ** get x,y,z of energy deposit
94 
95  1 continue
96  i_event=i_event+1
97  n_event=n_event+1
98  if(n_event.ge.100) then
99  write(6,*) i_event
100  n_event=0
101  endif
102  if(i_event.gt.60000) stop
103 
104 * read in from file (generated by MC)
105 
106  n_point=1
107  5 continue
108 
109  read(61,*,err=4,end=4) z_store(n_point),
110  > x_store(n_point),y_store(n_point)
111 
112  if(z_store(n_point).ge.0) then
113  x_store(n_point)=x_store(n_point)/10.
114  y_store(n_point)=y_store(n_point)/10.
115  n_point=n_point+1
116  goto 5
117  else
118  x_store(0)=x_store(n_point)
119  y_store(0)=y_store(n_point)
120  z_store(0)=z_store(n_point)
121  n_point=n_point-1
122  endif
123  inter_layer_new=x_store(0)
124 
125 * clean up hits (eliminate close by hits)
126 
127  do i1=1,n_point
128  f_tmp(i1)=1 ! 0 = discard point, 1 = keep point
129  enddo
130 
131  do i1=1,n_point-1
132  if(f_tmp(i1).eq.1) then
133  do i2=i1+1,n_point
134  if(z_store(i1).eq.z_store(i2).and.
135  > f_tmp(i2).eq.1) then
136  r1=sqrt((x_store(i1)-x_store(i2))**2+
137  > (y_store(i1)-y_store(i2))**2)
138  if(r1.lt.dist_cut) f_tmp(i2)=0
139  endif
140  enddo
141  endif
142  enddo
143 
144  n_tmp=0
145  do i1=1,n_point
146  if(f_tmp(i1).eq.1) then
147  n_tmp=n_tmp+1
148  x_tmp(n_tmp)=x_store(i1)
149  y_tmp(n_tmp)=y_store(i1)
150  z_tmp(n_tmp)=z_store(i1)
151  endif
152  enddo
153 
154  do i1=1,n_tmp
155  x_store(i1)=x_tmp(i1)
156  y_store(i1)=y_tmp(i1)
157  z_store(i1)=z_tmp(i1)
158  enddo
159  n_point=n_tmp
160 
161 * calculate distance of each point to closest point
162 
163  do i1=1,n_point
164  dist(i1)=999.
165  do i2=1,n_point
166  if(i2.ne.i1) then
167  if(z_store(i1).eq.z_store(i2)) then
168  r1=sqrt((x_store(i1)-x_store(i2))**2+(y_store(i1)-y_store(i2))**2)
169  if(r1.lt.dist(i1)) then
170  dist(i1)=r1
171  endif
172  endif
173  endif
174  enddo
175  enddo
176 
177 * reset pad array and pad charges
178 
179  do i1=0,150
180  do i2=0,150
181  do i3=0,150
182  hit(i3,i2,i1)=0
183  q_pad(i3,i2,i1)=0.
184  enddo
185  enddo
186  enddo
187 
188 * loop over all points in event
189 
190  do i1=1,n_point
191 
192 * get x,y
193 
194  x=x_store(i1)
195  y=y_store(i1)
196 
197 * determine id of pad hit
198 
199  i_pad_x=int(x)
200  i_pad_y=int(y)
201 
202 * determine minimum distance to boundaries in x and y
203 
204  x_crack(1)=x+48.0
205  x_crack(2)=96.0-(x+48.0)
206  x_scaler=amin1(x_crack(1),x_crack(2))
207 
208  y_scaler=999.
209  do i2=1,4
210  r1=abs(y-y_crack(i2)+48.8)
211  if(r1.le.y_scaler) y_scaler=r1
212  enddo
213 
214 * generate charge: Distributed like RABBIT data
215 
216  2 r1=rand(0)*12000.
217  r2=par(1)*r1**par(2)*exp(-r1*par(3))
218  r3=rand(0)*40.
219  if(r2.gt.r3) then
220  q=r1/1100.
221  else
222  goto 2
223  endif
224 
225 * Introduce offset
226 
227  q=q-q0
228 
229 * Adjust charge if close to other points
230 
231  if(dist(i1).gt.dist_cut.and.dist(i1).le.e_cut) then
232  r1=0.5*(1.+(dist(i1)-dist_cut)/(e_cut-dist_cut))
233  q=r1*q
234  endif
235 
236 * Adjust charge close to boundary
237 
238  if(x_scaler.lt.10.0) q=q*((x_scaler/10.0)**0.7)
239  if(y_scaler.lt.10.0) q=q*((y_scaler/10.0)**0.25)
240  if(y_scaler.lt.1.0) q=q/6.
241 
242 * Distribute charge according to STAR measurements
243 
244 
245 * integrate with n_int points
246 
247  i_int=0
248  q_check=0.
249  3 continue
250  x_int=2.*radius_max*(rand(0)-.5)
251  y_int=2.*radius_max*(rand(0)-.5)
252  r1=sqrt(x_int**2+y_int**2)
253  if(r1.gt.radius_max) goto 3
254  i_int=i_int+1
255 
256 * calculate charge
257 
258  q_int=q/real(n_int)*8.*
259  > ((1.-ratio)/q_integral1*exp(-0.5*((r1/sigma1)**2))+
260  > ratio/q_integral2*exp(-0.5*((r1/sigma2)**2)))
261 
262  q_check=q_check+q_int
263 * write(63,*) Q,q_int,r1
264 
265 * deposit charge on corresponding pad
266 
267  x_deposit=x+x_int+48.0 ! x of deposited charge in cm, centered onto plane
268  y_deposit=y+y_int+48.8 ! y of deposited charge in cm, centered onto plane
269 
270  if(x_deposit.lt.0.0.or.x_deposit.gt.96.0) q_int=0. ! set charge to zero if off the boards in x
271 
272  i2=int(x_deposit) ! pad index in x
273 
274  if(y_deposit.lt.0.27) then
275  q_int=0. ! set charge to zero if off the boards in y
276  elseif(y_deposit.gt.32.27.and.y_deposit.lt.32.8) then
277  q_int=0.
278  elseif(y_deposit.gt.64.8.and.y_deposit.lt.65.33) then
279  q_int=0.
280  elseif(y_deposit.gt.97.33) then
281  q_int=0.
282  endif
283 
284  if(y_deposit.le.32.27) then ! pad index in y
285  i3=int(y_deposit-0.27)
286  elseif(y_deposit.le.64.80) then
287  i3=int(y_deposit-0.8)
288  elseif(y_deposit.le.97.33) then
289  i3=int(y_deposit-1.33)
290  else
291  i3=0
292  endif
293 
294 * write(70,*) 'x,y,x_int,y_int',x,y,x_int,y_int
295 * write(70,*) 'x_deposit,y_deposit',x_deposit,y_deposit
296 * write(70,*) 'padx,pady',i2,i3
297 
298  i4=z_store(i1) ! pad index in z
299 
300  q_pad(i2,i3,i4)=q_pad(i2,i3,i4)+q_int
301  if(i_int.lt.n_int) goto 3
302 
303  enddo
304 
305 * determine pads hit
306 
307  do i1=0,150
308  do i2=0,150
309  do i3=0,150
310  if(q_pad(i1,i2,i3).gt.t) hit(i1,i2,i3)=1
311  enddo
312  enddo
313  enddo
314 
315 * write out hits
316 
317  write(62,*) i_event,inter_layer,-1,-1
318  do i3=0,150
319  do i2=0,200
320  do i1=0,200
321  if(hit(i1,i2,i3).eq.1) then
322  write(62,*) i_event,i1,i2,i3
323  endif
324  enddo
325  enddo
326  enddo
327 
328 * check total charge
329 
330  write(60,*) q,q_check
331 
332 * get next event
333 
334  inter_layer=inter_layer_new
335  goto 1
336 
337 * finish up
338 
339  4 continue
340 
341  stop
342  end
rpc_sim
program rpc_sim
Definition: RPC_sim_7.f:1
Generated by   doxygen 1.8.18