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