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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.4 - (show annotations) (download)
Fri Aug 4 08:18:06 2006 UTC (18 years, 5 months ago) by pam-fi
Branch: MAIN
Changes since 1.3: +27 -15 lines
some memory leak bugs fixed + CN computation modified

1 *************************************************************************
2 *
3 * Program reductionflight.f
4 *
5 * - reads readraw.f output files: LEVEL0 ntuple, and ped, sig and bad histograms
6 * - decodes raw data (DATATRACKER) using DSP ped, sig and bad values
7 * - looks for clusters information using ped, sig and bad values from
8 * DSP histograms
9 * - fills LEVEL1 ntuple
10 *
11 *************************************************************************
12
13 subroutine reductionflight(ierror)
14
15 include 'commontracker.f'
16 include 'level0.f'
17 include 'level1.f'
18 include 'common_reduction.f'
19 include 'calib.f'
20
21 integer ierror
22 ierror = 0
23
24 * -------------------------------------------------------
25 * STRIP MASK
26 * -------------------------------------------------------
27
28 c call stripmask !called later, after CN computation
29 call init_level1
30
31 good1 = good0
32 c--------------------------------------------------
33 c read the variable DATATRACKER from LEVEL0
34 c and fill the variable ADC (inverting view 11)
35 c--------------------------------------------------
36 call filladc(iflag)
37 if(iflag.ne.0)then
38 good1=0
39 c if(DEBUG)print*,'event ',eventn(1),' >>>>> decode ERROR'
40 ierror = 220
41 goto 200
42 endif
43
44 c--------------------------------------------------
45 c computes common noise for each VA1
46 c (excluding strips affected by signal,
47 c tagged with the flag CLSTR)
48 c--------------------------------------------------
49 do iv=1,nviews
50 do ik=1,nva1_view
51 cn(iv,ik)=0
52 mask_vk_ev(iv,ik)=1
53 iflag=0
54 if(mask_vk(iv,ik).eq.1)call cncomp(iv,ik,iflag)
55 c if(iflag.ne.0)good1=0
56 if(iflag.ne.0)then
57 mask_vk_ev(iv,ik)=0
58 ierror = 220
59 endif
60 enddo
61 enddo
62 c if(good1.eq.0)then
63 c ierror = 220
64 c endif
65
66 call stripmask !compute mask(i,j,k)
67 c---------------------------------------------
68 c loops on views, VA1 and strips,
69 c and computes strips signals using
70 c badstrip, pedestals, and
71 c sigma informations from histograms
72 c---------------------------------------------
73 flag_shower = .false.
74 ind=1 !clsignal array index
75 do iv=1,nviews !loop on views
76 do is=1,nstrips_view !loop on strips (1)
77 if(mod(iv,2).eq.1) then
78 C=== > Y view
79 value(is)= -(DBLE(adc(iv,nvk(is),nst(is)))
80 $ -cn(iv,nvk(is))-pedestal(iv,nvk(is),nst(is)))
81 $ *mask(iv,nvk(is),nst(is))
82 clseedcut(is)=clcuty*sigma(iv,nvk(is),nst(is))
83 $ *mask(iv,nvk(is),nst(is))
84 clinclcut(is)=incuty*sigma(iv,nvk(is),nst(is))
85 $ *mask(iv,nvk(is),nst(is))
86 ccc print*,"value(",is,")(reduction)= ",value(is)
87 else
88 C=== > X view
89 value(is)= (DBLE(adc(iv,nvk(is),nst(is)))
90 $ -cn(iv,nvk(is))-pedestal(iv,nvk(is),nst(is)))
91 $ *mask(iv,nvk(is),nst(is))
92 clseedcut(is)=clcutx*sigma(iv,nvk(is),nst(is))
93 $ *mask(iv,nvk(is),nst(is))
94 clinclcut(is)=incutx*sigma(iv,nvk(is),nst(is))
95 $ *mask(iv,nvk(is),nst(is))
96 endif
97 c$$$ print*,iv,is,' --- ',adc(iv,nvk(is),nst(is)),cn(iv,nvk(is))
98 c$$$ $ ,pedestal(iv,nvk(is),nst(is)),value(is)
99 c$$$ $ ,sigma(iv,nvk(is),nst(is))
100 c if(value(is).gt.clseedcut(is))
101 c $ print*,iv,is,' --- (ADC_PED_CN) ',value(is),clseedcut(is)
102 enddo !end loop on strips (1)
103 call search_cluster(iv)
104 if(flag_shower.eqv..true.)then
105 call init_level1
106 good1=0
107 goto 200 !jump to next event
108 endif
109 enddo ! end loop on views
110 do iv=1,nviews
111 do ik=1,nva1_view
112 cnev(iv,ik)=cn(iv,ik) !assigns computed CN to ntuple variables
113 cnevflag(iv,ik)=cnflag(iv,ik) !assigns computed CN to ntuple variables
114 ccc print*,"cnev(",iv,",",ik,")(reduction)= ",cnev(iv,ik)
115 enddo
116 enddo
117 C---------------------------------------------
118 C come here if GOOD1=0
119 C or the event has too many clusters
120 C---------------------------------------------
121 200 continue
122 c------------------------------------------------------------------------
123 c
124 c closes files and exits
125 c
126 c------------------------------------------------------------------------
127 RETURN
128 END
129
130 ***...***...***...***...***...***...***...***...***...***...***...***...***...***...***...***
131 *
132 *
133 *
134 *
135 *
136 *
137 *
138 *
139 *
140 ***...***...***...***...***...***...***...***...***...***...***...***...***...***...***...***
141
142
143 subroutine init_level1
144
145 include 'commontracker.f'
146 include 'level1.f'
147 include 'level0.f'
148
149 good1=0
150 nclstr1=0
151 totCLlength=0
152 do ic=1,nclstrmax
153 view(ic)=0
154 ladder(ic)=0
155 indstart(ic)=0
156 indmax(ic)=0
157 maxs(ic)=0
158 mult(ic)=0
159 dedx(ic)=0
160 enddo
161 do id=1,maxlength !???
162 clsignal(id)=0.
163 enddo
164 do iv=1,nviews
165 c crc1(iv)=0
166 do ik=1,nva1_view
167 cnev(iv,ik)=0
168 enddo
169 enddo
170
171 return
172 end
173 *---***---***---***---***---***---***---***---***
174 *
175 *
176 *
177 *
178 *
179 *---***---***---***---***---***---***---***---***
180
181 subroutine search_cluster(iv)
182
183 include 'commontracker.f'
184 include 'common_reduction.f'
185 include 'level0.f'
186 include 'level1.f'
187 include 'calib.f'
188
189
190
191 c local variables
192 integer rmax,lmax !estremi del cluster
193 integer rstop,lstop !per decidere quali strip includere nel cluster
194 ! oltre il seed
195 integer first,last,diff !per includere le strip giuste... !???
196
197 integer multtemp !temporary multiplicity variable
198
199 integer CLlength !lunghezza in strip del cluster
200
201 external nst
202
203 c------------------------------------------------------------------------
204 c looks for clusters on each view
205 C : CERCO STRIP SOPRA CLSEEDCUT, POI SCORRO A DX FINCHE'
206 c NON TROVO
207 C STRIP PIU' BASSA (in segnale/rumore)
208 C => L'ULTIMA DELLA SERIE CRESCENTE
209 C (LA PIU' ALTA) E' IL
210 C CLUSTER SEED. POI SCORRO A SX E DX INCLUDENDO TUTTE
211 C LE STRIP (FINO A 17 AL
212 C MAX) CHE SUPERANO CLINCLCUT.
213 C QUANDO CERCO IL CLUSTER SEED SUCCESSIVO SALTO LA STRIP
214 C ADIACENTE A DESTRA
215 C DELL'ULTIMO CLUSTER SEED (CHE SARA' NECESSARIAMENTE
216 C PIU' BASSA) E PRENDO
217 C COME SEED UNA STRIP SOLO SE IL SUO SEGNALE E'
218 C MAGGIORE DI QUELLO DELLA STRIP
219 C PRECEDENTE (PRATICAMENTE PER EVITARE CHE L'ULTIMA
220 C STRIP DI UN GRUPPO DI STRIP
221 C TUTTE SOPRA IL CLSEEDCUT VENGA AUTOMATICAMENTE PRESA
222 C COME SEED... DEVE ESSERE
223 C PRESA SOLO SE IL CLUSTER E' DOUBLE PEAKED...)
224 c------------------------------------------------------------------------
225 c 6 ottobre 2003
226 c Elena: CLSEEDCUT = 7 (old value 10)
227 c Elena: CLINCLCUT = 4 (old value 5)
228
229 iseed=-999 !cluster seed index initialization
230
231 do jl=1,nladders_view !1..3 !loops on ladders
232 first=1+nstrips_ladder*(jl-1) !1,1025,2049
233 last=nstrips_ladder*jl !1024,2048,3072
234 c X views have 1018 strips instead of 1024
235 if(mod(iv,2).eq.0) then
236 first=first+3
237 last=last-3
238 endif
239 do is=first,last !loop on strips in each ladder
240 if(is.le.iseed+1) goto 220
241 c-----------------------------------------
242 c after a cluster seed as been found,
243 c look for next one skipping one strip on the right
244 c (i.e. look for double peak cluster)
245 c-----------------------------------------
246 if(is.ne.first) then
247 if(value(is).le.value(is-1)) goto 220
248 endif
249 c-----------------------------------------
250 c skips cluster seed
251 c finding if strips values are descreasing (a strip
252 c can be a cluster seed only if previous strip value
253 c is lower)
254 c-----------------------------------------
255 if(value(is).gt.clseedcut(is)) then
256 ccc print*,"value(",is,")=",value(is),
257 ccc $ " .gt.clseedcut(",is,")=",clseedcut(is)
258 c-----------------------------------------
259 c possible SEED...
260 c-----------------------------------------
261 itemp=is
262 if(itemp.eq.last) goto 230 !estremo...
263 do while(value(itemp)
264 $ /sigma(iv,nvk(itemp),nst(itemp))
265 $ .le.value(itemp+1)
266 $ /sigma(iv,nvk(itemp+1),nst(itemp+1))) !BIAS: aggiustare il caso uguale!???
267 itemp=itemp+1
268 if(itemp.eq.last) goto 230 !stops if reaches last strip
269 enddo ! of the ladder
270 230 continue
271 c-----------------------------------------
272 c fownd SEED!!!
273 c-----------------------------------------
274 iseed=itemp
275 c----------------------------------------------------------
276 c after finding a cluster seed, checks also adjacent strips,
277 C and marks the ones exceeding clinclcut
278 c----------------------------------------------------------
279 ir=iseed !indici destro
280 il=iseed ! e sinistro
281
282 rmax=ir !estremo destro del cluster
283 lmax=il ! e sinistro
284
285 rstop=0 !initialize flags used to exit from
286 lstop=0 ! inclusion loop
287
288 do while(lstop.eq.0.or.rstop.eq.0) !shifts left and right from
289 ir=ir+1 !position index for strips on right side of
290 ! cluster seed
291 il=il-1 !and for left side
292 c------------------------------------------------------------------------
293 c checks for last or first strip of the ladder
294 c------------------------------------------------------------------------
295 if(ir.gt.last) then !when index goes beyond last strip
296 rstop=1 ! of the ladder, change rstop flag in order
297 ! to "help" exiting from loop
298 endif
299
300 if(il.lt.first) then !idem when index goes beyond
301 lstop=1 ! first strip of the ladder
302 endif
303
304 c------------------------------------------------------------------------
305 c check for clusters including more than nclstrp strips
306 c------------------------------------------------------------------------
307 if((rmax-lmax+1).ge.nclstrp) then
308 goto 210 !exits inclusion loop:
309 ! lmax and rmax maintain last value
310 ! NB .ge.!???
311 endif
312 c------------------------------------------------------------------------
313 c marks strips exceeding inclusion cut
314 c------------------------------------------------------------------------
315 if(rstop.eq.0) then !if last strip of the ladder or last
316 ! over-cut strip has not been reached
317 if(value(ir).gt.clinclcut(ir)) then !puts in rmax the
318 rmax=ir ! last right over-cut strip
319 else
320 rstop=1 !otherwise cluster ends on right and rstop
321 endif ! flag=1 signals it
322 endif
323 if(lstop.eq.0) then
324 if(value(il).gt.clinclcut(il)) then
325 lmax=il
326 else
327 lstop=1
328 endif
329 endif
330
331 enddo !ends strip inclusion loop
332 210 continue !jumps here if more than nclstrp have been included
333
334 multtemp=rmax-lmax+1 !stores multiplicity in temp
335 ! variable. NB rmax and lmax can change later in
336 ! order to include enough strips to calculate eta3
337 ! and eta4. so mult is not always equal to cllength
338 c------------------------------------------------------------------------
339 c NB per essere sicuro di poter calcolare eta3 e eta4 devo includere
340 c sempre e comunque le 2 strip adiacenti al cluster seed e quella
341 c adiacente ulteriore dalla parte della piu' alta fra queste due
342 c (vedi oltre...)!???
343 c------------------------------------------------------------------------
344
345 c nel caso di estremi del ladder...!???
346
347 c ho meno di 4 strip nel cluster --> se sono sui bordi o quasi del ladder
348 c costruisco il cluster ad hoc e poi esco, se non sono sui bordi o quasi
349 c vado oltre (aggiungero' quindi strip a sx e dx in modo da poter calcolare
350 c eta3e4)
351 if((rmax-lmax+1).lt.4) then
352
353 if(iseed.eq.first) then !estremi...
354 rmax=iseed+2 !NB in questo modo puo' anche capitare di
355 lmax=iseed ! includere strip sotto taglio di inclusione
356 goto 250 ! che non serviranno per eta3e4!???
357 endif
358
359 if(iseed.eq.last) then !estremi...
360 rmax=iseed
361 lmax=iseed-2 !NB 2 e non 3, perche' altrimenti sarei in
362 goto 250 ! ((rmax-lmax+1).lt.4).eq.false. !???
363 endif !NMB questo e' l'unico caso di cllength=3!???
364
365 if(iseed.eq.first+1) then !quasi estremi...
366 rmax=iseed+2
367 lmax=iseed-1
368 goto 250
369 endif
370 if(iseed.eq.last-1) then
371 rmax=iseed+1
372 lmax=iseed-2
373 goto 250
374 endif
375 c se ho 4 o piu' strip --> se sono sui bordi esco, se sono sui quasi bordi
376 c includo la strip del bordo
377 else
378
379 if(iseed.eq.first) goto 250 !estremi... non includo altro
380 if(iseed.eq.last) goto 250
381 if(iseed.eq.first+1) then !quasi estremi... mi assicuro di
382 lmax=first ! avere le strip adiacenti al seed
383 if((rmax-lmax+1).gt.nclstrp) rmax=rmax-1 !NB effetto
384 goto 250 ! coperta: se la lunghezza del cluster era gia'
385 endif ! al limite (nclstrp), per poter aggiungere questa
386 ! strip a sinistra devo toglierne una a destra...!???
387 if(iseed.eq.last-1) then
388 rmax=last
389 if((rmax-lmax+1).gt.nclstrp) lmax=lmax+1
390 goto 250
391 endif
392 endif
393 c------------------------------------------------------------------------
394 c be sure to include in the cluster the cluster seed with its 2 adjacent
395 c strips, and the one adjacent to the greatest between this two strip, as the
396 c fourth one. if the strips have the same value (!) the fourth one is chosen
397 c as the one having the greatest value between the second neighbors
398 c------------------------------------------------------------------------
399 if(value(iseed+1).eq.value(iseed-1)) then
400 if(value(iseed+2).ge.value(iseed-2)) then !??? qui cmq c'e'
401 diff=(iseed+2)-rmax
402 if(diff.gt.0) then
403 rmax=rmax+diff
404 if((rmax-lmax+1).gt.nclstrp) then
405 lmax=rmax-nclstrp+1
406 endif
407 endif
408 diff=(iseed-1)-lmax
409 if(diff.lt.0) then
410 lmax=lmax+diff
411 if((rmax-lmax+1).gt.nclstrp) then
412 rmax=lmax+nclstrp-1
413 endif
414 endif
415 else
416 diff=(iseed-2)-lmax
417 if(diff.lt.0) then
418 lmax=lmax+diff
419 if((rmax-lmax+1).gt.nclstrp) then
420 rmax=lmax+nclstrp-1
421 endif
422 endif
423 diff=(iseed+1)-rmax
424 if(diff.gt.0) then
425 rmax=rmax+diff
426 if((rmax-lmax+1).gt.nclstrp) then
427 lmax=rmax-nclstrp+1
428 endif
429 endif
430 endif
431 elseif(value(iseed+1).gt.value(iseed-1)) then
432 c !??? sposto il limite del cluster a destra per includere sempre le strip
433 c necessarie al calcolo di eta-i
434 c se il cluster diventa troppo lungo lo accorcio a sinistra per avere non piu'
435 c di nclstrp (in questo caso sono sicuro di aver gia' incluso le strip
436 c necessarie al calcolo di eta-i a sinistra, quindi se voglio posso uscire)
437 diff=(iseed+2)-rmax
438 if(diff.gt.0) then
439 rmax=rmax+diff
440 if((rmax-lmax+1).gt.nclstrp) then
441 lmax=rmax-nclstrp+1
442 c goto 250
443 endif
444 endif
445 diff=(iseed-1)-lmax
446 if(diff.lt.0) then
447 lmax=lmax+diff
448 if((rmax-lmax+1).gt.nclstrp) then
449 rmax=lmax+nclstrp-1
450 c goto 250 !inutile!???
451 endif
452 endif
453 else
454 diff=(iseed-2)-lmax
455 if(diff.lt.0) then
456 lmax=lmax+diff
457 if((rmax-lmax+1).gt.nclstrp) then
458 rmax=lmax+nclstrp-1
459 c goto 250
460 endif
461 endif
462 diff=(iseed+1)-rmax
463 if(diff.gt.0) then
464 rmax=rmax+diff
465 if((rmax-lmax+1).gt.nclstrp) then
466 lmax=rmax-nclstrp+1
467 c goto 250 !inutile!???
468 endif
469 endif
470 endif
471 250 continue
472
473 c--------------------------------------------------------
474 c fills cluster variables
475 c--------------------------------------------------------
476 nclstr1=nclstr1+1 !cluster number
477 ccc print*,nclstr1,multtemp
478 if(nclstr1.gt.nclstrmax) then !too many clusters for the event:
479 if(verbose)print*,'Event ',eventn(1),
480 $ ': more than ',nclstrmax,' clusters'
481 good1=0 ! event
482 nclstr1=0
483 totCLlength=0
484 flag_shower = .true.
485 goto 2000
486 endif
487 view(nclstr1)=iv !vista del cluster
488 ladder(nclstr1)=nld(iseed,iv) !ladder a cui appartiene il cluster seed
489 maxs(nclstr1)=iseed !strip del cluster seed
490 mult(nclstr1)=multtemp !molteplicita'
491
492 indstart(nclstr1)=ind !posizione dell'inizio del cluster nell'
493 ! array clsignal
494 indmax(nclstr1)=indstart(nclstr1)+(iseed-lmax) !posizione del
495 ! cluster seed nell'array clsignal
496
497 CLlength=rmax-lmax+1 !numero di strip del cluster
498 totCLlength=totCLlength+CLlength
499 dedx(nclstr1)=0
500 do j=lmax,rmax !stores sequentially cluter strip values in
501 clsignal(ind)=value(j) ! clsignal array
502 ind=ind+1
503 c if(value(j).gt.0)
504 if(value(j).gt.clinclcut(j))
505 $ dedx(nclstr1)=dedx(nclstr1)+value(j) !cluster charge
506 enddo
507 c--------------------------------------------------------
508 c
509 c--------------------------------------------------------
510 endif !end possible seed conditio
511 220 continue !jumps here to skip strips left of last seed
512
513 enddo ! end loop on strips
514 enddo !end loop on ladders
515 2000 continue
516 return
517 end
518
519
520 *---***---***---***---***---***---***---***---***
521 *
522 *
523 *
524 *
525 *
526 *---***---***---***---***---***---***---***---***
527
528
529 subroutine stripmask
530
531 * this routine set va1 and single-strip masks,
532 * on the basis of the VA1 mask saved in the DB
533 *
534 * mask(nviews,nva1_view,nstrips_va1) !strip mask
535 * mask_vk(nviews,nva1_view) !VA1 mask
536 *
537 include 'commontracker.f'
538 c include 'level1.f'
539 include 'common_reduction.f'
540 include 'calib.f'
541
542 * init mask
543 do iv=1,nviews
544 do ivk=1,nva1_view
545 do is=1,nstrips_va1
546 c mask(iv,ivk,is) = mask_vk(iv,ivk)
547 mask(iv,ivk,is) = mask_vk(iv,ivk) * mask_vk_ev(iv,ivk)
548 enddo
549 enddo
550 enddo
551
552
553 return
554 end
555

  ViewVC Help
Powered by ViewVC 1.1.23