/[PAMELA software]/tracker/ground/utilities/template/template.f
ViewVC logotype

Annotation of /tracker/ground/utilities/template/template.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1 - (hide annotations) (download)
Wed Mar 8 15:00:49 2006 UTC (18 years, 9 months ago) by pam-fi
Branch point for: MAIN, trk-ground
Initial revision

1 pam-fi 1.1
2     *************************************************************************
3     * Program template.f
4     *
5     * - template to access LEVEL2 tracker data
6     * and use the track parameters to evaluate the
7     * the particle trajectory in the apparatus
8     *
9     *************************************************************************
10    
11     program template
12    
13    
14     include 'trk_level2.f'
15    
16     c----- HBOOK
17     INTEGER HMEM
18     parameter (NWPAWC=8500000)
19     common/PAWC/HMEM(NWPAWC)
20    
21     parameter (ntp_level2=22)
22     c------------------------------------------------------------------------
23     c
24     c local variables
25     c
26     c------------------------------------------------------------------------
27     character*74 data_file !data file name
28     character*74 data_dir !data directory
29     character*74 data_file_level2
30    
31     parameter (lun_data_level2=72) !data file id number
32    
33    
34     logical general,single,column,trackok
35     integer l1,s1,l6,s6
36    
37    
38     COMMON/QUEST/IQUEST(100)
39    
40     c =======================================
41     c variables for tracking routine
42     c =======================================
43     parameter(NPOINT_MAX=100)
44     DOUBLE PRECISION ZIN(NPOINT_MAX)
45     DOUBLE PRECISION XOUT(NPOINT_MAX),YOUT(NPOINT_MAX)
46     DOUBLE PRECISION AL_P(5)
47     c =======================================
48    
49     c define TOF Z-coordinates
50     parameter (NPTOF=3)
51     REAL ZTOF(NPTOF)
52     DATA ZTOF/55.,25.,-23/ !NB!!! from TOP to BOTTOM
53    
54     c------------------------------------------------------------------------
55     c
56     c HBOOK initialization
57     c
58     c------------------------------------------------------------------------
59    
60     call HLIMIT(NWPAWC)
61     c------------------------------------------------------------------------
62     c
63     c reads input informations
64     c
65     c------------------------------------------------------------------------
66     111 format(a)
67     print*,'Data dir:'
68     read(*,111)data_dir
69     print*,data_dir
70     print*,'Data file'
71     print*,'(without estention: output_YYMMDD_XXX )'
72     read(*,*)data_file
73     print*,data_file
74     minevent=1
75     print*,'Maximum number of events to be analysed:'
76     read(*,*) ntotev
77     print*,ntotev
78     print*,'-----------------------------------------'
79     print*,''
80     c------------------------------------------------------------------------
81     c
82     c opens level2 file
83     c
84     c------------------------------------------------------------------------
85     504 format(a,a,'_level2.rz')
86     write(data_file_level2,504)
87     $ data_dir(1:LNBLNK(data_dir))
88     $ ,data_file(1:LNBLNK(data_file))
89     print*,'__________ opening LEVEL2 rz file __________'
90     print*,data_file_level2
91     IQUEST(10)=65000
92     call HROPEN(lun_data_level2,
93     $ 'LEVEL2',data_file_level2,'QP',4096,istat) !opens rz
94     if(istat.ne.0) goto 19
95     print*,'reading LEVEL2 n-tuple...'
96     call HRIN(ntp_level2,9999,0)
97    
98     * -----------------------------------------------
99     CALL HBNAME(ntp_level2,' ',0,'$CLEAR')
100     CALL HBNAME(ntp_level2,'GENERAL',GOOD2,'$SET')
101     CALL HBNAME(ntp_level2,'CPU',pkt_type,'$SET')
102     CALL HBNAME(ntp_level2,'TRACKS',ntrk,'$SET')
103     CALL HBNAME(ntp_level2,'SINGLETX',nclsx,'$SET')
104     CALL HBNAME(ntp_level2,'SINGLETY',nclsy,'$SET')
105     * -----------------------------------------------
106    
107     call HPRNTU(ntp_level2)
108     call HNOENT(ntp_level2,iemax0)
109     print*,'ok'
110     print*,' '
111    
112    
113     c------------------------------------------------------------------------
114     c
115     c =======================
116     c read magnetic field map
117     c =======================
118     c
119     c------------------------------------------------------------------------
120     print*,'- read magnetic field map'
121     print*,' '
122     call read_B
123     print*,' '
124    
125     call HCDIR('//LEVEL2',' ')
126    
127     c------------------------------------------------------------------------
128     c ////////////////////
129     c
130     c start loop on events
131     c
132     c \\\\\\\\\\\\\\\\\\\\
133     c------------------------------------------------------------------------
134     maxevent=minevent+ntotev
135     do iev = minevent,MIN(iemax0,maxevent) !loop on events
136    
137     * ----------------------------------------------
138     call HGNT(ntp_level2,iev,ierr) !reads an event
139     if(ierr.ne.0) goto 21
140     * ----------------------------------------------
141    
142    
143     *=======> INSTERT HERE YOUR CODE
144    
145    
146     * initializations
147     do itof=1,NPTOF
148     XOUT(itof) = 0
149     YOUT(itof) = 0
150     enddo
151    
152     * ########################################
153     * GENERAL CUTS
154     general = .true.
155     if (
156     + .not.GOOD2.or.
157     + ntrk.eq.0.or.
158     + .false.) general = .false.
159     if(.not.general)goto 100
160    
161     * ########################################
162     * SINGLE TRACK
163     single = .false.
164     itrk = 0
165     if(ntrk.eq.1)then
166     single = .true.
167     itrk = 1
168     elseif(ntrk.eq.2..and.image(1).ne.0)then
169     single = .true.
170     * -----------------
171     * CHOSE THE IMAGE!!
172     * -----------------
173     itrk = 1
174     if(CHI2(image(1)).lt.CHI2(1))itrk = 2
175     * -----------------
176     endif
177     if(.not.single)goto 100
178    
179     * ########################################
180     * GEOMETRY
181    
182     l1 = int( p1_ladder(XV(1,itrk)))
183     l6 = int( p1_ladder(XV(6,itrk)))
184     s1 = int( p1_sensor(YV(1,itrk)))
185     s6 = int( p1_sensor(YV(6,itrk)))
186     column=.false.
187     if(
188     + l1.eq.l6.and.
189     + l1.ne.0.and.
190     + s1.eq.s6.and.
191     + s1.ne.0.and.
192     + .true.)column=.true.
193     c$$$ if(.not.column)goto 100
194    
195    
196     * ########################################
197     * TRACK SELECTION
198     trackok = .false.
199     npx = 0
200     npy = 0
201     do i=1,6
202     npx = npx + XGOOD(i,itrk)
203     npy = npy + YGOOD(i,itrk)
204     enddo
205     * POINTS ----------------------------------
206     if(
207     + npy.ge.4.and.
208     + npx.ge.5.and.
209     c + npy.ge.6.and.
210     c + npx.ge.6.and.
211     + XGOOD(1,itrk).eq.1.and.
212     + XGOOD(6,itrk).eq.1.and.
213     + CHI2(itrk).gt.0.and.
214     + .true.) trackok = .true.
215    
216     rig=abs(1./AL(5,itrk))
217    
218     * CHI^2 ------------------------------------
219     chicut = 10. !
220     chicut = chicut *(
221     $ 3.7032*exp(-0.41911E-01*rig)
222     $ +0.85355
223     $ +10.905*exp(-0.61893*rig) )
224     if(
225     + CHI2(itrk).gt.chicut.and.
226     + .true.)trackok = .false.
227    
228     * CLUSTERS OUT OF THE TRACK --------------------------
229     ncls=nclsx+nclsy
230     if(ncls.gt.0)trackok = .false.
231    
232     * RIGIDITY ----------------------------------
233     * *******************************************
234     if(
235     + rig.lt.5..and.
236     + .true.)trackok = .false.
237     * *******************************************
238    
239     if(.not.trackok)goto 100
240    
241     *
242    
243    
244     print*,'----------------------------------'
245     print*,'Event ',iev,' selected'
246     * assigned input parameters for track routine
247     * 1) Z-coordinates where the trajectory is evaluated
248     do itof=1,NPTOF
249     ZIN(itof) = ZTOF(itof)
250     enddo
251     * 2) track status vector
252     do i=1,5
253     AL_P(i) = AL(i,itrk)
254     enddo
255     * -------- *** tracking routine *** --------
256     call track(NPTOF,ZIN,XOUT,YOUT,AL_P,IFAIL)
257     * ------------------------------------------
258     do itof=1,NPTOF
259     print*,'S',itof,' --- (z,x,y) = '
260     $ ,ZIN(itof),XOUT(itof),YOUT(itof)
261     enddo
262    
263    
264    
265    
266    
267     100 continue
268     enddo !end loop on events
269     GOTO 9000 !got to end
270    
271     c------------------------------------------------------------------------
272     c
273     c data file opening error
274     c
275     c------------------------------------------------------------------------
276     19 continue
277    
278     print*,' '
279     print*,'ERROR OPENING DATA FILE: ',data_file
280     print*,' '
281     print*,' '
282    
283     goto 9000 !the end
284    
285     c------------------------------------------------------------------------
286     c
287     c level2 ntuple event reading error
288     c
289     c------------------------------------------------------------------------
290    
291     21 continue
292    
293     print*,' '
294     print*,'ERROR WHILE READING LEVEL2 NTUPLE, AT EVENT
295     $ : ',iev
296     print*,' '
297     print*,' '
298    
299     goto 9000 !the end
300    
301     c------------------------------------------------------------------------
302     c
303     c closes files and exits
304     c
305     c------------------------------------------------------------------------
306    
307     9000 continue
308    
309     call HCDIR('//LEVEL2',' ')
310     call HREND('level2')
311     close(lun_data_level2)
312    
313     return
314     end
315    
316    
317     c include 'access_level2.f'
318     c include 'track.f'
319    
320     function p1_ladder(x)
321    
322     real xlimit(4)
323     data xlimit/-7.995,-2.665,2.665,7.995/
324     parameter (xmargin=0.5)
325     c parameter (xmargin=1.)
326     real xlo,xhi
327    
328    
329    
330     nladder=0
331     do i=2,4
332     xlo=xlimit(i-1)+xmargin
333     xhi=xlimit(i)-xmargin
334     if(
335     + x.lt.xhi.and.
336     + x.ge.xlo)then
337     nladder=i-1
338     goto 10
339     endif
340     enddo
341     10 continue
342     p1_ladder=nladder
343    
344     return
345     end
346    
347     function p1_sensor(y)
348    
349     real ylimit(3)
350     data ylimit/-7.,0.,7./
351     parameter (ymargin=0.5)
352     c parameter (ymargin=2.)
353     real ylo,yhi
354    
355     nsensor=0
356     do i=2,3
357     ylo=ylimit(i-1)+ymargin
358     yhi=ylimit(i)-ymargin
359     if(
360     + y.lt.yhi.and.
361     + y.ge.ylo)then
362     nsensor=i-1
363     goto 10
364     endif
365     enddo
366     10 continue
367     p1_sensor=nsensor
368    
369     return
370     end
371    

  ViewVC Help
Powered by ViewVC 1.1.23