/[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.3 - (show annotations) (download)
Tue Jun 27 13:57:41 2006 UTC (18 years, 5 months ago) by mocchiut
Branch: MAIN
CVS Tags: v1r01beta, v1r01
Changes since 1.2: +1 -1 lines
Workaround in CaloCore and ToFCore programs + new ToF default calibration file to process flight data

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

  ViewVC Help
Powered by ViewVC 1.1.23