/[PAMELA software]/DarthVader/TrackerLevel2/src/F77/functions.f
ViewVC logotype

Contents of /DarthVader/TrackerLevel2/src/F77/functions.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.4 - (show annotations) (download)
Thu May 24 16:45:48 2007 UTC (17 years, 6 months ago) by pam-fi
Branch: MAIN
Changes since 1.3: +15 -0 lines
several upgrades

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
137 function coordsi(istrip,view)
138 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
139 * it gives the strip coordinate in micrometers,
140 * knowing the strip number (1..3072) and the view
141 * number. the origin of the coordinate is on the
142 * centre of the sensor the strip belongs to.
143 * the axes directions are the same as in the PAMELA
144 * reference frame (i.e.: the 11th view coordinate
145 * direction has to be inverted here)
146 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
147
148 c integer is,view,istrip
149
150 integer view,is,istrip
151 real coordsi
152
153 include 'commontracker.f'
154
155 c NB mettere il 1024 nel commontracker...!???
156
157
158
159 is=istrip !it stores istrip number
160 is=mod(is-1,1024)+1 !it puts all clusters on a single ladder
161
162 coordsi=0.
163
164 if(mod(view,2).eq.0) then !X view
165
166 c if((is.le.3).or.(is.ge.1022)) then !X has 1018 strips...
167 c print*,'functions: WARNING: false X strip: strip ',is
168 c endif
169
170 is=is-3 !4 =< is =< 1021 --> 1 =< is =< 1018
171
172 edge=edgeX
173 dim=SiDimX
174
175 elseif(mod(view,2).eq.1) then !Y view
176
177 edge=edgeY
178 dim=SiDimY
179
180 c$$$ if(view.eq.11) then !INVERSIONE!???
181 c$$$ is=1025-is
182 c$$$ endif
183
184 endif
185
186 p=pitch(view)
187
188 coord1=(is-1)*p !referred to 1st sensor strip
189 coord1=coord1+edge !referred to sensor edge
190
191 coordsi=coord1-dim/2 !referred to the centre of the sensor
192
193 if(view.eq.11) then !INVERSION: it puts y axis in the same direction for all views
194 coordsi=-coordsi
195 endif
196
197 end
198
199
200 c------------------------------------------------------------------------
201
202
203 function acoordsi(strip,view)
204 *
205 * same as COORDSI, but accept a real value of strip!!!
206 *
207 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
208 * it gives the strip coordinate in micrometers,
209 * knowing the strip number (1..3072) and the view
210 * number. the origin of the coordinate is on the
211 * centre of the sensor the strip belongs to.
212 * the axes directions are the same as in the PAMELA
213 * reference frame (i.e.: the 11th view coordinate
214 * direction has to be inverted here)
215 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
216
217 c integer is,view,istrip
218
219 integer view,is,istrip
220 real coordsi,acoordsi
221 real strip,stripladder
222
223
224 include 'commontracker.f'
225
226 c NB mettere il 1024 nel commontracker...!???
227
228 istrip = int(strip+0.5) !istrip stores the closest integer to strip
229
230 is=istrip !it stores istrip number
231 is=mod(is-1,1024)+1 !it puts all clusters on a single ladder
232
233 coordsi=0.
234
235 if(mod(view,2).eq.0) then !X view
236
237 c if((is.le.3).or.(is.ge.1022)) then !X has 1018 strips...
238 c print*,'functions: WARNING: false X strip: strip ',is
239 c endif
240
241 is=is-3 !4 =< is =< 1021 --> 1 =< is =< 1018
242
243 edge=edgeX
244 dim=SiDimX
245
246 elseif(mod(view,2).eq.1) then !Y view
247
248 edge=edgeY
249 dim=SiDimY
250
251 c$$$ if(view.eq.11) then !INVERSIONE!???
252 c$$$ is=1025-is
253 c$$$ endif
254
255 endif
256
257
258 stripladder = float(is)+(strip-float(istrip))!cluster position relative to ladder
259 p=pitch(view)
260
261 ccccc coord1=(is-1)*p !referred to 1st sensor strip
262 coord1=(stripladder-1)*p !referred to 1st sensor strip
263 coord1=coord1+edge !referred to sensor edge
264 acoordsi=coord1-dim/2 !referred to the centre of the sensor
265
266 if(view.eq.11) then !INVERSION: it puts y axis in the same direction for all views
267 acoordsi=-acoordsi
268 endif
269
270 end
271
272
273
274 c------------------------------------------------------------------------
275
276
277 function coord(coordsi,view,ladder,sen)
278 * it gives the coordinate in
279 * micrometers, knowing the coordinate in the sensor
280 * frame, the view, the ladder and the sensor numbers.
281 * the origin is in the centre of the magnet (PAMELA
282 * reference frame)
283
284 include 'commontracker.f'
285 include 'common_tracks.f'
286
287 integer view,ladder,sen
288 integer sx,sy,sz
289
290 real coord,coordsi,trasl
291
292 c$$$c parameter (offset=4365.) !??? ! in um !CONTROLLARE SE HA SENSO:
293 c$$$ ! dalle misure sul piano dovrebbe essere 4970,
294 c$$$ ! dallo shift dei residui viene 4365
295 c$$$ ! va messo .ne.0. se in mech_sensor assegno ai
296 c$$$ ! sensori del sesto piano coordinate Y uguali
297 c$$$ ! a quelle degli altri sensori
298 c$$$ parameter (offset=0.) !??? altrimenti se il sesto piano ha coordinate
299 c$$$ ! Y diverse offset dovrebbe essere .eq.0.
300 c$$$ ! CONTROLLARE CON I GRAFICI DEI RESIDUI!!!
301
302
303 coord=0.
304
305 sx=ladder
306 sy=sen
307 sz=npl(view)
308
309 if(mod(view,2).eq.0) then !X view
310
311 trasl=x_mech_sensor(sz,sx,sy) !in mm
312
313 elseif(mod(view,2).eq.1) then !Y view
314
315 trasl=y_mech_sensor(sz,sx,sy) !in mm
316
317 c$$$ if(view.eq.11) then !INVERSIONE!???INUTILE, ne e' gia' tenuto conto
318 c$$$ coordsi=coordsi+offset ! in y_mech_pos...
319 c$$$ endif
320
321 endif
322
323 coord=coordsi+trasl*1000.
324
325 end
326
327
328 c------------------------------------------------------------------------
329 c------------------------------------------------------------------------
330
331 c double precision version of the above subroutine
332
333 double precision function dcoord(coordsi,view,ladder,sen) !it gives the coordinate in
334 ! micrometers, knowing the coordinate in the sensor
335 ! frame, the view, the ladder and the sensor numbers.
336 ! the origin is in the centre of the magnet (PAMELA
337 ! reference frame)
338
339 include 'commontracker.f'
340 include 'common_tracks.f'
341
342 integer view,ladder,sen
343 integer sx,sy,sz
344
345 c double precision dcoord
346 double precision coordsi,trasl
347
348 c$$$c parameter (offset=4365.) !??? ! in um !CONTROLLARE SE HA SENSO:
349 c$$$ ! dalle misure sul piano dovrebbe essere 4970,
350 c$$$ ! dallo shift dei residui viene 4365
351 c$$$ ! va messo .ne.0. se in mech_sensor assegno ai
352 c$$$ ! sensori del sesto piano coordinate Y uguali
353 c$$$ ! a quelle degli altri sensori
354 c$$$ parameter (offset=0.) !??? altrimenti se il sesto piano ha coordinate
355 c$$$ ! Y diverse offset dovrebbe essere .eq.0.
356 c$$$ ! CONTROLLARE CON I GRAFICI DEI RESIDUI!!!
357
358
359 dcoord=0.
360
361 if(
362 $ ladder.lt.1.or.
363 $ ladder.gt.nladders_view.or.
364 $ sen.lt.1.or.
365 $ sen.gt.2.or.
366 $ view.lt.1.or.
367 $ view.gt.nviews.or.
368 $ .false.)then
369 print*,'dcoord ---> wrong input: ladder ',ladder
370 $ ,' sensor ',sen
371 $ ,' view ',view
372 return
373 endif
374
375 sx=ladder
376 sy=sen
377 sz=npl(view)
378
379
380 if(mod(view,2).eq.0) then !X view
381
382 trasl=x_mech_sensor(sz,sx,sy) !in mm
383
384 elseif(mod(view,2).eq.1) then !Y view
385
386 trasl=y_mech_sensor(sz,sx,sy) !in mm
387
388 c$$$ if(view.eq.11) then !INVERSIONE!???INUTILE, ne e' gia' tenuto conto
389 c$$$ dcoordsi=dcoordsi+offset ! in y_mech_pos...
390 c$$$ endif
391
392 endif
393
394 dcoord=coordsi+trasl*1000.
395
396 end
397
398
399 c------------------------------------------------------------------------
400 integer function nsatstrips(ic)
401 *--------------------------------------------------------------
402 * this function returns the number of saturated strips
403 * inside a cluster
404 *--------------------------------------------------------------
405 include 'commontracker.f'
406 include 'level1.f'
407 include 'calib.f'
408
409
410
411 integer nsat
412 nsat = 0
413 iv=VIEW(ic)
414 if(mod(iv,2).eq.1)incut=incuty
415 if(mod(iv,2).eq.0)incut=incutx
416 istart = INDSTART(IC)
417 istop = TOTCLLENGTH
418 if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
419 do i = INDMAX(IC),istart,-1
420 c cut = incut*CLSIGMA(i)
421 c if(CLSIGNAL(i).ge.cut)then
422 if( (mod(iv,2).eq.1.and.CLADC(i).lt.ADCsatx)
423 $ .or.
424 $ (mod(iv,2).eq.0.and.CLADC(i).gt.ADCsaty) )then
425 nsat = nsat +1
426 else
427 goto 10
428 endif
429 enddo
430 10 continue
431 do i = INDMAX(IC)+1,istop
432 c cut = incut*CLSIGMA(i)
433 c if(CLSIGNAL(i).ge.cut)then
434 if( (mod(iv,2).eq.1.and.CLADC(i).lt.ADCsatx)
435 $ .or.
436 $ (mod(iv,2).eq.0.and.CLADC(i).gt.ADCsaty) )then
437 nsat = nsat +1
438 else
439 goto 20
440 endif
441 enddo
442 20 continue
443
444 nsatstrips = nsat
445 return
446 end
447
448 c------------------------------------------------------------------------
449 integer function nbadstrips(ncog,ic)
450 *--------------------------------------------------------------
451 * this function returns the number of BAD strips
452 * inside a cluster:
453 * - if NCOG=0, the number BAD strips inside the whole cluster
454 * are given, according to the cluster multiplicity
455 *
456 * - if NCOG>0, the number BAD strips is evaluated using NCOG
457 * strips, even if they have a negative signal (according to Landi)
458 *--------------------------------------------------------------
459 include 'commontracker.f'
460 include 'level1.f'
461 include 'calib.f'
462
463 integer nbad
464 nbad = 0
465
466 if (ncog.gt.0) then
467
468 * --> signal of the central strip
469 sc = CLSIGNAL(INDMAX(ic)) !center
470 * signal of adjacent strips
471 sl1 = -100000 !left 1
472 if(
473 $ (INDMAX(ic)-1).ge.INDSTART(ic)
474 $ )
475 $ sl1 = CLSIGNAL(INDMAX(ic)-1)
476
477 sl2 = -100000 !left 2
478 if(
479 $ (INDMAX(ic)-2).ge.INDSTART(ic)
480 $ )
481 $ sl2 = CLSIGNAL(INDMAX(ic)-2)
482
483 sr1 = -100000 !right 1
484 if(
485 $ (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))
486 $ .or.
487 $ (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)
488 $ )
489 $ sr1 = CLSIGNAL(INDMAX(ic)+1)
490
491 sr2 = -100000 !right 2
492 if(
493 $ (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))
494 $ .or.
495 $ (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)
496 $ )
497 $ sr2 = CLSIGNAL(INDMAX(ic)+2)
498
499 if(ncog.ge.1)nbad = nbad+1-CLBAD(INDMAX(ic))
500
501 if(ncog.ge.2)then
502 if(sl1.gt.sr1.and.(INDMAX(ic)-1).ge.1)then
503 nbad=nbad+1-CLBAD(INDMAX(ic)-1)
504 elseif(sl1.le.sr1.and.(INDMAX(ic)+1).le.NCLSTR1)then
505 nbad=nbad+1-CLBAD(INDMAX(ic)+1)
506 endif
507 endif
508
509 if(ncog.ge.3)then
510 if(sl1.gt.sr1.and.(INDMAX(ic)+1).le.NCLSTR1)then
511 nbad=nbad+1-CLBAD(INDMAX(ic)+1)
512 elseif(sl1.le.sr1.and.(INDMAX(ic)-1).ge.1)then
513 c if(INDMAX(ic)-1.eq.0)
514 c $ print*,' ======= ',sl2,sl1,sc,sr1,sr2
515 nbad=nbad+1-CLBAD(INDMAX(ic)-1)
516 endif
517 endif
518
519 if(ncog.ge.4)then
520 if(sl2.gt.sr2.and.(INDMAX(ic)-2).ge.1)then
521 nbad=nbad+1-CLBAD(INDMAX(ic)-2)
522 elseif(sl2.le.sr2.and.(INDMAX(ic)+2).le.NCLSTR1)then
523 nbad=nbad+1-CLBAD(INDMAX(ic)+2)
524 endif
525 endif
526
527 c if(ncog.ge.5)then
528 c print*,'function CLBAD(NCOG,IC) ==> WARNING!! NCOG=',NCOG
529 c $ ,' not implemented'
530 c endif
531
532 elseif(ncog.eq.0)then
533 * =========================
534 * COG computation
535 * =========================
536
537
538
539 iv=VIEW(ic)
540 if(mod(iv,2).eq.1)incut=incuty
541 if(mod(iv,2).eq.0)incut=incutx
542
543 istart = INDSTART(IC)
544 istop = TOTCLLENGTH
545 if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
546 nbad = 0
547 c$$$ do i=istart,istop
548 c$$$ cut = incut*CLSIGMA(i)
549 c$$$ if(CLSIGNAL(i).ge.cut)nbad = nbad +1 -CLBAD(i)
550 c$$$ enddo
551 do i = INDMAX(IC),istart,-1
552 cut = incut*CLSIGMA(i)
553 if(CLSIGNAL(i).ge.cut)then
554 nbad = nbad +1 -CLBAD(i)
555 else
556 goto 10
557 endif
558 enddo
559 10 continue
560 do i = INDMAX(IC)+1,istop
561 cut = incut*CLSIGMA(i)
562 if(CLSIGNAL(i).ge.cut)then
563 nbad = nbad +1 -CLBAD(i)
564 else
565 goto 20
566 endif
567 enddo
568 20 continue
569
570 else
571
572 c print*,'function CLBAD(NCOG,IC) ==> WARNING!! NCOG=',NCOG
573 c $ ,' not implemented'
574
575
576 endif
577
578 nbadstrips = nbad
579
580 return
581 end

  ViewVC Help
Powered by ViewVC 1.1.23