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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.8 - (hide 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 pam-fi 1.6 ************************************************************************
2 mocchiut 1.1 *
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 mocchiut 1.7 view = view
69 mocchiut 1.1 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 bongi 1.5 c------------------------------------------------------------------------
137 mocchiut 1.1
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 mocchiut 1.3 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 mocchiut 1.1
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 mocchiut 1.3 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 mocchiut 1.1
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 bongi 1.5 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 mocchiut 1.1 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 pam-fi 1.4 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 mocchiut 1.1 sx=ladder
456     sy=sen
457     sz=npl(view)
458    
459 pam-fi 1.4
460 mocchiut 1.1 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 pam-fi 1.2 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 mocchiut 1.8 real incut
491 pam-fi 1.2 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 mocchiut 1.1
543 mocchiut 1.8 real incut
544 pam-fi 1.2 integer nbad
545     nbad = 0
546 mocchiut 1.1
547 pam-fi 1.2 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 mocchiut 1.3 c if(ncog.ge.5)then
609     c print*,'function CLBAD(NCOG,IC) ==> WARNING!! NCOG=',NCOG
610     c $ ,' not implemented'
611     c endif
612 pam-fi 1.2
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 mocchiut 1.3 c print*,'function CLBAD(NCOG,IC) ==> WARNING!! NCOG=',NCOG
654     c $ ,' not implemented'
655 pam-fi 1.2
656    
657     endif
658    
659     nbadstrips = nbad
660    
661     return
662     end

  ViewVC Help
Powered by ViewVC 1.1.23