1 |
************************************************************************* |
2 |
* |
3 |
* functions.f |
4 |
* |
5 |
* - !??? |
6 |
* |
7 |
* needs: |
8 |
* - !??? |
9 |
* |
10 |
* output variables: |
11 |
* - !??? |
12 |
* |
13 |
* to be called inside !??? |
14 |
* |
15 |
* |
16 |
* MODIFIED in order to have in input a |
17 |
* REAL-defined strip number instead of INTEGER |
18 |
* |
19 |
************************************************************************* |
20 |
|
21 |
|
22 |
function pitch(view) !it gives the strip pitch, knowing the view number |
23 |
|
24 |
real pitch |
25 |
integer view |
26 |
|
27 |
include './commontracker.f' |
28 |
|
29 |
if(mod(view,2).eq.0) then !X |
30 |
pitch=pitchX |
31 |
else !Y |
32 |
pitch=pitchY |
33 |
endif |
34 |
|
35 |
end |
36 |
|
37 |
|
38 |
|
39 |
c------------------------------------------------------------------------ |
40 |
|
41 |
|
42 |
|
43 |
function npl(view) !it gives the plane number, knowing the view number. |
44 |
! plane 1 = views 11+12, calorimeter side |
45 |
! ... |
46 |
! plane 6 = views 1+2, TRD side |
47 |
integer npl,view |
48 |
|
49 |
npl=7-(INT((view-1)/2)+1) |
50 |
|
51 |
end |
52 |
|
53 |
|
54 |
|
55 |
c------------------------------------------------------------------------ |
56 |
|
57 |
|
58 |
|
59 |
function nld(istrip,view) |
60 |
* it gives the number of the ladder, knowing the |
61 |
* strip number (1..3072) and the view number. |
62 |
* the first strip belongs to ladder 1 |
63 |
|
64 |
integer istrip,view,nld |
65 |
|
66 |
include './commontracker.f' |
67 |
|
68 |
|
69 |
nld=INT((istrip-1)/nstrips_ladder)+1 |
70 |
|
71 |
end |
72 |
|
73 |
|
74 |
c------------------------------------------------------------------------ |
75 |
|
76 |
|
77 |
function nviewx(iplane) !it gives the view number of a X plane |
78 |
|
79 |
integer nviewx,iplane |
80 |
|
81 |
nviewx=2*(7-iplane) |
82 |
|
83 |
end |
84 |
|
85 |
|
86 |
c------------------------------------------------------------------------ |
87 |
|
88 |
function nviewy(iplane) !it gives the view number of a Y plane |
89 |
|
90 |
integer nviewy,iplane |
91 |
|
92 |
nviewy=2*(7-iplane)-1 |
93 |
|
94 |
end |
95 |
|
96 |
c------------------------------------------------------------------------ |
97 |
|
98 |
|
99 |
|
100 |
|
101 |
function nvk(istrip) |
102 |
|
103 |
* it gives the number of the VA1, knowing the strip |
104 |
* number (1..3072). |
105 |
* the first strip belongs to VA1 1 |
106 |
integer istrip,nvk |
107 |
|
108 |
include './commontracker.f' |
109 |
|
110 |
nvk=INT((istrip-1)/nstrips_va1)+1 |
111 |
|
112 |
end |
113 |
|
114 |
|
115 |
|
116 |
c------------------------------------------------------------------------ |
117 |
|
118 |
|
119 |
|
120 |
function nst(istrip) |
121 |
|
122 |
* it gives the VA1 strip, knowing the strip number |
123 |
* (1..3072). |
124 |
* the first strip belongs to VA1 1 |
125 |
|
126 |
integer istrip,nst |
127 |
|
128 |
include './commontracker.f' |
129 |
|
130 |
nst=INT(mod((istrip-1),nstrips_va1))+1 |
131 |
|
132 |
|
133 |
end |
134 |
|
135 |
|
136 |
c$$$c------------------------------------------------------------------------ |
137 |
c$$$ |
138 |
c$$$c NB: le coordinate in mech_pos.dat sono calcolate a partire da alcuni dati |
139 |
c$$$c contenuti in commontracker.f. forse si puo' evitare mech_pos.dat e mettere |
140 |
c$$$c tutto in commontracker.f |
141 |
c$$$ |
142 |
c$$$ |
143 |
c$$$ subroutine mech_sensor !it reads sensors coordinates (in PAMELA reference |
144 |
c$$$ ! frame) from a text file and it uses them to fill |
145 |
c$$$ ! x/y/z_mech_sensor variables, taking into account |
146 |
c$$$ ! last plane inversion |
147 |
c$$$ |
148 |
c$$$ include './commontracker.f' |
149 |
c$$$ include './common_tracks.f' |
150 |
c$$$ |
151 |
c$$$ real xvec(nladders_view),yvec(2),zvec(nplanes) |
152 |
c$$$ |
153 |
c$$$ integer id !file identifier |
154 |
c$$$ logical od !.true. if the specified unit is connected to a file |
155 |
c$$$ |
156 |
c$$$ do id=20,100,1 !opens the file using a free file id |
157 |
c$$$ inquire (id, opened=od) |
158 |
c$$$ if(.not.od) goto 666 |
159 |
c$$$ enddo |
160 |
c$$$ 666 continue |
161 |
c$$$ |
162 |
c$$$ open(id,FILE='../common/mech_pos.dat') !sensors centres coordinates in mm in |
163 |
c$$$ ! PAMELA reference frame: |
164 |
c$$$ ! the first plane is the one with lowest Z (the one |
165 |
c$$$ ! nearest the calorimeter) |
166 |
c$$$ ! the first ladder is the one with lowest X (the |
167 |
c$$$ ! one on which the first X strip is) |
168 |
c$$$ ! the first sensor is the one with lowest Y (the |
169 |
c$$$ ! one on which the first Y strip is) for planes |
170 |
c$$$ ! 2..6. for plane 1 the first sensor has higher Y |
171 |
c$$$ |
172 |
c$$$ read(20,*) xvec |
173 |
c$$$ read(20,*) yvec |
174 |
c$$$ read(20,*) zvec |
175 |
c$$$ |
176 |
c$$$ do i=1,nplanes |
177 |
c$$$ do j=1,nladders_view |
178 |
c$$$ do k=1,2 |
179 |
c$$$ x_mech_sensor(i,j,k)=xvec(j) |
180 |
c$$$ y_mech_sensor(i,j,k)=yvec(k) |
181 |
c$$$ z_mech_sensor(i,j,k)=zvec(i) |
182 |
c$$$ if(i.eq.1) then !y coordinates of first plane (11th view) are |
183 |
c$$$ y_mech_sensor(i,j,k)=-yvec(k) ! exchanged due to last plane inversion |
184 |
c$$$ endif |
185 |
c$$$ enddo |
186 |
c$$$ enddo |
187 |
c$$$ enddo |
188 |
c$$$ |
189 |
c$$$ close(id) |
190 |
c$$$ |
191 |
c$$$ |
192 |
c$$$c$$$ ! *** INIZIO DEBUG *** |
193 |
c$$$c$$$ do i=1,6 |
194 |
c$$$c$$$ do j=1,3 |
195 |
c$$$c$$$ do k=1,2 |
196 |
c$$$c$$$ c print*,x_mech_sensor(1,j,k) |
197 |
c$$$c$$$ print*,y_mech_sensor(i,j,k) |
198 |
c$$$c$$$ c print*,z_mech_sensor(i,j,k) |
199 |
c$$$c$$$ enddo |
200 |
c$$$c$$$ enddo |
201 |
c$$$c$$$ print*,' ' |
202 |
c$$$c$$$ enddo |
203 |
c$$$c$$$ ! *** FINE DEBUG *** |
204 |
c$$$ |
205 |
c$$$ |
206 |
c$$$ return |
207 |
c$$$ end |
208 |
|
209 |
|
210 |
c$$$c------------------------------------------------------------------------ |
211 |
c$$$ |
212 |
c$$$c NB: le coordinate in mech_pos.dat sono calcolate a partire da alcuni dati |
213 |
c$$$c contenuti in commontracker.f. forse si puo' evitare mech_pos.dat e mettere |
214 |
c$$$c tutto in commontracker.f |
215 |
c$$$ |
216 |
c$$$ |
217 |
c$$$ subroutine mech_sensor |
218 |
c$$$c !it reads sensors coordinates (in PAMELA reference |
219 |
c$$$c ! frame) from a text file and it uses them to fill |
220 |
c$$$c ! x/y/z_mech_sensor variables, taking into account |
221 |
c$$$c ! last plane inversion |
222 |
c$$$ |
223 |
c$$$ include './commontracker.f' |
224 |
c$$$ include './common_tracks.f' |
225 |
c$$$ |
226 |
c$$$ real xvec(nladders_view),yvec(2),zvec(nplanes) |
227 |
c$$$ |
228 |
c$$$ integer id !file identifier |
229 |
c$$$ logical od !.true. if the specified unit is connected to a file |
230 |
c$$$ |
231 |
c$$$ do id=20,100,1 !opens the file using a free file id |
232 |
c$$$ inquire (id, opened=od) |
233 |
c$$$ if(.not.od) goto 666 |
234 |
c$$$ enddo |
235 |
c$$$ 666 continue |
236 |
c$$$ |
237 |
c$$$c open(id,FILE='../common/mech_pos.dat') !sensors centres coordinates in mm in |
238 |
c$$$c open(id,FILE='source/common/mech_pos.dat') |
239 |
c$$$c call system('cp $TRK_GRND/source/common/mech_pos.dat .') |
240 |
c$$$ print *,'Opening file: mech_pos.dat' |
241 |
c$$$ open(id,FILE='./bin-aux/mech_pos.dat',IOSTAT=iostat) |
242 |
c$$$c !sensors centres coordinates in mm in |
243 |
c$$$c ! PAMELA reference frame: |
244 |
c$$$c ! the first plane is the one with lowest Z (the one |
245 |
c$$$c ! nearest the calorimeter) |
246 |
c$$$c ! the first ladder is the one with lowest X (the |
247 |
c$$$c ! one on which the first X strip is) |
248 |
c$$$c ! the first sensor is the one with lowest Y (the |
249 |
c$$$c ! one on which the first Y strip is) for planes |
250 |
c$$$c ! 2..6. for plane 1 the first sensor has higher Y |
251 |
c$$$ |
252 |
c$$$ if(iostat.ne.0)then |
253 |
c$$$ print*,'MECH_SENSOR: *** Error in opening file ***' |
254 |
c$$$ return |
255 |
c$$$ endif |
256 |
c$$$ |
257 |
c$$$ read(id,*) xvec |
258 |
c$$$ read(id,*) yvec |
259 |
c$$$ read(id,*) zvec |
260 |
c$$$ |
261 |
c$$$ do i=1,nplanes |
262 |
c$$$ do j=1,nladders_view |
263 |
c$$$ do k=1,2 |
264 |
c$$$ x_mech_sensor(i,j,k)=xvec(j) |
265 |
c$$$ y_mech_sensor(i,j,k)=yvec(k) |
266 |
c$$$ z_mech_sensor(i,j,k)=zvec(i) |
267 |
c$$$ if(i.eq.1) then !y coordinates of first plane (11th view) are |
268 |
c$$$ y_mech_sensor(i,j,k)=-yvec(k) ! exchanged due to last plane inversion |
269 |
c$$$ endif |
270 |
c$$$ enddo |
271 |
c$$$ enddo |
272 |
c$$$ enddo |
273 |
c$$$ |
274 |
c$$$ close(id) |
275 |
c$$$c call system('rm -f mech_pos.dat') |
276 |
c$$$ |
277 |
c$$$c$$$ ! *** INIZIO DEBUG *** |
278 |
c$$$c$$$ do i=1,6 |
279 |
c$$$c$$$c do j=1,3 |
280 |
c$$$c$$$ do k=1,2 |
281 |
c$$$c$$$ j=1 |
282 |
c$$$c$$$c print*,x_mech_sensor(1,j,k) |
283 |
c$$$c$$$ print*,y_mech_sensor(i,j,k) |
284 |
c$$$c$$$c print*,z_mech_sensor(i,j,k) |
285 |
c$$$c$$$ enddo |
286 |
c$$$c$$$c enddo |
287 |
c$$$c$$$ print*,' ' |
288 |
c$$$c$$$ enddo |
289 |
c$$$c$$$ ! *** FINE DEBUG *** |
290 |
c$$$ |
291 |
c$$$ |
292 |
c$$$ return |
293 |
c$$$ end |
294 |
c$$$ |
295 |
|
296 |
c------------------------------------------------------------------------ |
297 |
|
298 |
|
299 |
function coordsi(istrip,view) |
300 |
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * |
301 |
* it gives the strip coordinate in micrometers, |
302 |
* knowing the strip number (1..3072) and the view |
303 |
* number. the origin of the coordinate is on the |
304 |
* centre of the sensor the strip belongs to. |
305 |
* the axes directions are the same as in the PAMELA |
306 |
* reference frame (i.e.: the 11th view coordinate |
307 |
* direction has to be inverted here) |
308 |
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * |
309 |
|
310 |
c integer is,view,istrip |
311 |
|
312 |
integer view,is,istrip |
313 |
real coordsi |
314 |
|
315 |
include './commontracker.f' |
316 |
|
317 |
c NB mettere il 1024 nel commontracker...!??? |
318 |
|
319 |
|
320 |
|
321 |
is=istrip !it stores istrip number |
322 |
is=mod(is-1,1024)+1 !it puts all clusters on a single ladder |
323 |
|
324 |
coordsi=0. |
325 |
|
326 |
if(mod(view,2).eq.0) then !X view |
327 |
|
328 |
if((is.le.3).or.(is.ge.1022)) then !X has 1018 strips... |
329 |
print*,'functions: WARNING: false X strip: strip ',is |
330 |
endif |
331 |
|
332 |
is=is-3 !4 =< is =< 1021 --> 1 =< is =< 1018 |
333 |
|
334 |
edge=edgeX |
335 |
dim=SiDimX |
336 |
|
337 |
elseif(mod(view,2).eq.1) then !Y view |
338 |
|
339 |
edge=edgeY |
340 |
dim=SiDimY |
341 |
|
342 |
c$$$ if(view.eq.11) then !INVERSIONE!??? |
343 |
c$$$ is=1025-is |
344 |
c$$$ endif |
345 |
|
346 |
endif |
347 |
|
348 |
p=pitch(view) |
349 |
|
350 |
coord1=(is-1)*p !referred to 1st sensor strip |
351 |
coord1=coord1+edge !referred to sensor edge |
352 |
|
353 |
coordsi=coord1-dim/2 !referred to the centre of the sensor |
354 |
|
355 |
if(view.eq.11) then !INVERSION: it puts y axis in the same direction for all views |
356 |
coordsi=-coordsi |
357 |
endif |
358 |
|
359 |
end |
360 |
|
361 |
|
362 |
c------------------------------------------------------------------------ |
363 |
|
364 |
|
365 |
function acoordsi(strip,view) |
366 |
* |
367 |
* same as COORDSI, but accept a real value of strip!!! |
368 |
* |
369 |
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * |
370 |
* it gives the strip coordinate in micrometers, |
371 |
* knowing the strip number (1..3072) and the view |
372 |
* number. the origin of the coordinate is on the |
373 |
* centre of the sensor the strip belongs to. |
374 |
* the axes directions are the same as in the PAMELA |
375 |
* reference frame (i.e.: the 11th view coordinate |
376 |
* direction has to be inverted here) |
377 |
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * |
378 |
|
379 |
c integer is,view,istrip |
380 |
|
381 |
integer view,is,istrip |
382 |
real coordsi,acoordsi |
383 |
real strip,stripladder |
384 |
|
385 |
|
386 |
include './commontracker.f' |
387 |
|
388 |
c NB mettere il 1024 nel commontracker...!??? |
389 |
|
390 |
istrip = int(strip+0.5) !istrip stores the closest integer to strip |
391 |
|
392 |
is=istrip !it stores istrip number |
393 |
is=mod(is-1,1024)+1 !it puts all clusters on a single ladder |
394 |
|
395 |
coordsi=0. |
396 |
|
397 |
if(mod(view,2).eq.0) then !X view |
398 |
|
399 |
if((is.le.3).or.(is.ge.1022)) then !X has 1018 strips... |
400 |
print*,'functions: WARNING: false X strip: strip ',is |
401 |
endif |
402 |
|
403 |
is=is-3 !4 =< is =< 1021 --> 1 =< is =< 1018 |
404 |
|
405 |
edge=edgeX |
406 |
dim=SiDimX |
407 |
|
408 |
elseif(mod(view,2).eq.1) then !Y view |
409 |
|
410 |
edge=edgeY |
411 |
dim=SiDimY |
412 |
|
413 |
c$$$ if(view.eq.11) then !INVERSIONE!??? |
414 |
c$$$ is=1025-is |
415 |
c$$$ endif |
416 |
|
417 |
endif |
418 |
|
419 |
|
420 |
stripladder = float(is)+(strip-float(istrip))!cluster position relative to ladder |
421 |
p=pitch(view) |
422 |
|
423 |
ccccc coord1=(is-1)*p !referred to 1st sensor strip |
424 |
coord1=(stripladder-1)*p !referred to 1st sensor strip |
425 |
coord1=coord1+edge !referred to sensor edge |
426 |
acoordsi=coord1-dim/2 !referred to the centre of the sensor |
427 |
|
428 |
if(view.eq.11) then !INVERSION: it puts y axis in the same direction for all views |
429 |
acoordsi=-acoordsi |
430 |
endif |
431 |
|
432 |
end |
433 |
|
434 |
|
435 |
|
436 |
c------------------------------------------------------------------------ |
437 |
|
438 |
|
439 |
function coord(coordsi,view,ladder,sen) |
440 |
* it gives the coordinate in |
441 |
* micrometers, knowing the coordinate in the sensor |
442 |
* frame, the view, the ladder and the sensor numbers. |
443 |
* the origin is in the centre of the magnet (PAMELA |
444 |
* reference frame) |
445 |
|
446 |
include './commontracker.f' |
447 |
include './common_tracks.f' |
448 |
|
449 |
integer view,ladder,sen |
450 |
integer sx,sy,sz |
451 |
|
452 |
real coord,coordsi,trasl |
453 |
|
454 |
c$$$c parameter (offset=4365.) !??? ! in um !CONTROLLARE SE HA SENSO: |
455 |
c$$$ ! dalle misure sul piano dovrebbe essere 4970, |
456 |
c$$$ ! dallo shift dei residui viene 4365 |
457 |
c$$$ ! va messo .ne.0. se in mech_sensor assegno ai |
458 |
c$$$ ! sensori del sesto piano coordinate Y uguali |
459 |
c$$$ ! a quelle degli altri sensori |
460 |
c$$$ parameter (offset=0.) !??? altrimenti se il sesto piano ha coordinate |
461 |
c$$$ ! Y diverse offset dovrebbe essere .eq.0. |
462 |
c$$$ ! CONTROLLARE CON I GRAFICI DEI RESIDUI!!! |
463 |
|
464 |
|
465 |
coord=0. |
466 |
|
467 |
sx=ladder |
468 |
sy=sen |
469 |
sz=npl(view) |
470 |
|
471 |
if(mod(view,2).eq.0) then !X view |
472 |
|
473 |
trasl=x_mech_sensor(sz,sx,sy) !in mm |
474 |
|
475 |
elseif(mod(view,2).eq.1) then !Y view |
476 |
|
477 |
trasl=y_mech_sensor(sz,sx,sy) !in mm |
478 |
|
479 |
c$$$ if(view.eq.11) then !INVERSIONE!???INUTILE, ne e' gia' tenuto conto |
480 |
c$$$ coordsi=coordsi+offset ! in y_mech_pos... |
481 |
c$$$ endif |
482 |
|
483 |
endif |
484 |
|
485 |
coord=coordsi+trasl*1000. |
486 |
|
487 |
end |
488 |
|
489 |
|
490 |
c------------------------------------------------------------------------ |
491 |
c------------------------------------------------------------------------ |
492 |
|
493 |
c double precision version of the above subroutine |
494 |
|
495 |
double precision function dcoord(coordsi,view,ladder,sen) !it gives the coordinate in |
496 |
! micrometers, knowing the coordinate in the sensor |
497 |
! frame, the view, the ladder and the sensor numbers. |
498 |
! the origin is in the centre of the magnet (PAMELA |
499 |
! reference frame) |
500 |
|
501 |
include './commontracker.f' |
502 |
include './common_tracks.f' |
503 |
|
504 |
integer view,ladder,sen |
505 |
integer sx,sy,sz |
506 |
|
507 |
c double precision dcoord |
508 |
double precision coordsi,trasl |
509 |
|
510 |
c$$$c parameter (offset=4365.) !??? ! in um !CONTROLLARE SE HA SENSO: |
511 |
c$$$ ! dalle misure sul piano dovrebbe essere 4970, |
512 |
c$$$ ! dallo shift dei residui viene 4365 |
513 |
c$$$ ! va messo .ne.0. se in mech_sensor assegno ai |
514 |
c$$$ ! sensori del sesto piano coordinate Y uguali |
515 |
c$$$ ! a quelle degli altri sensori |
516 |
c$$$ parameter (offset=0.) !??? altrimenti se il sesto piano ha coordinate |
517 |
c$$$ ! Y diverse offset dovrebbe essere .eq.0. |
518 |
c$$$ ! CONTROLLARE CON I GRAFICI DEI RESIDUI!!! |
519 |
|
520 |
|
521 |
dcoord=0. |
522 |
|
523 |
sx=ladder |
524 |
sy=sen |
525 |
sz=npl(view) |
526 |
|
527 |
if(mod(view,2).eq.0) then !X view |
528 |
|
529 |
trasl=x_mech_sensor(sz,sx,sy) !in mm |
530 |
|
531 |
elseif(mod(view,2).eq.1) then !Y view |
532 |
|
533 |
trasl=y_mech_sensor(sz,sx,sy) !in mm |
534 |
|
535 |
c$$$ if(view.eq.11) then !INVERSIONE!???INUTILE, ne e' gia' tenuto conto |
536 |
c$$$ dcoordsi=dcoordsi+offset ! in y_mech_pos... |
537 |
c$$$ endif |
538 |
|
539 |
endif |
540 |
|
541 |
dcoord=coordsi+trasl*1000. |
542 |
|
543 |
end |
544 |
|
545 |
|
546 |
c------------------------------------------------------------------------ |
547 |
|
548 |
|
549 |
|