/[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.8 - (show annotations) (download)
Tue May 15 14:31:45 2012 UTC (12 years, 6 months ago) by mocchiut
Branch: MAIN
CVS Tags: v10RED, v10REDr01, HEAD
Changes since 1.7: +2 -1 lines
Reprocessing bugs fixed, ToF bugs fixed, new software versions, new quaternions, IGRF bug fixed, code cleanup

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

  ViewVC Help
Powered by ViewVC 1.1.23