--- DarthVader/ToFLevel2/src/tofl2com.for 2008/03/03 09:51:03 1.7 +++ DarthVader/ToFLevel2/src/tofl2com.for 2008/03/31 19:24:04 1.8 @@ -26,6 +26,9 @@ C then in a second step we check the residuals of the single C measurements, reject if > 10 sigma, calculate chi2 and "quality" C beta is taken as good if chi2<20 and quality>10 +C mar-08 WM: Call to "newbeta" changed, now a flag tells the function if the +C call comes from "tofl2com" or form "toftrack" +C mar-08 WM: Bug found in dEdx if check_charge>1 C****************************************************************************** INTEGER FUNCTION TOFL2COM() @@ -1054,7 +1057,7 @@ tof11(left,i,iadc) = tof11(left,i,iadc)*cos(theta13) xkorr = atten(left,11,i,yhelp) c write(40+i,*) yhelp,xkorr - if (iz.le.1) xkorr=xkorr/hepratio + xkorr=xkorr/hepratio adctof_c(ch11a(i),hb11a(i))=tof11(left,i,iadc)/xkorr endif @@ -1063,7 +1066,7 @@ tof11(right,i,iadc) = tof11(right,i,iadc)*cos(theta13) xkorr = atten(right,11,i,yhelp) c write(40+i,*) yhelp,xkorr - if (iz.le.1) xkorr=xkorr/hepratio + xkorr=xkorr/hepratio adctof_c(ch11b(i),hb11b(i))=tof11(right,i,iadc)/xkorr endif ENDIF @@ -1080,7 +1083,7 @@ tof12(left,i,iadc) = tof12(left,i,iadc)*cos(theta13) xkorr = atten(left,12,i,xhelp) c write(50+i,*) xhelp,xkorr - if (iz.le.1) xkorr=xkorr/hepratio + xkorr=xkorr/hepratio adctof_c(ch12a(i),hb12a(i))=tof12(left,i,iadc)/xkorr endif @@ -1089,7 +1092,7 @@ tof12(right,i,iadc) = tof12(right,i,iadc)*cos(theta13) xkorr = atten(right,12,i,xhelp) c write(50+i,*) xhelp,xkorr - if (iz.le.1) xkorr=xkorr/hepratio + xkorr=xkorr/hepratio adctof_c(ch12b(i),hb12b(i))=tof12(right,i,iadc)/xkorr endif ENDIF @@ -1108,7 +1111,7 @@ tof21(left,i,iadc) = tof21(left,i,iadc)*cos(theta13) xkorr = atten(left,21,i,xhelp) c write(60+i,*) xhelp,xkorr - if (iz.le.1) xkorr=xkorr/hepratio + xkorr=xkorr/hepratio adctof_c(ch21a(i),hb21a(i))=tof21(left,i,iadc)/xkorr endif @@ -1118,7 +1121,7 @@ xkorr=adcx21(right,i,1)*exp(xhelp/adcx21(right,i,2)) xkorr = atten(right,21,i,xhelp) c write(60+i,*) xhelp,xkorr - if (iz.le.1) xkorr=xkorr/hepratio + xkorr=xkorr/hepratio adctof_c(ch21b(i),hb21b(i))=tof21(right,i,iadc)/xkorr endif ENDIF @@ -1136,7 +1139,7 @@ tof22(left,i,iadc) = tof22(left,i,iadc)*cos(theta13) xkorr = atten(left,22,i,yhelp) c write(70+i,*) yhelp,xkorr - if (iz.le.1) xkorr=xkorr/hepratio + xkorr=xkorr/hepratio adctof_c(ch22a(i),hb22a(i))=tof22(left,i,iadc)/xkorr endif @@ -1145,7 +1148,7 @@ tof22(right,i,iadc) = tof22(right,i,iadc)*cos(theta13) xkorr = atten(right,22,i,yhelp) c write(70+i,*) yhelp,xkorr - if (iz.le.1) xkorr=xkorr/hepratio + xkorr=xkorr/hepratio adctof_c(ch22b(i),hb22b(i))=tof22(right,i,iadc)/xkorr endif ENDIF @@ -1164,7 +1167,7 @@ tof31(left,i,iadc) = tof31(left,i,iadc)*cos(theta13) xkorr = atten(left,31,i,yhelp) c write(80+i,*) yhelp,xkorr - if (iz.le.1) xkorr=xkorr/hepratio + xkorr=xkorr/hepratio adctof_c(ch31a(i),hb31a(i))=tof31(left,i,iadc)/xkorr endif @@ -1173,7 +1176,7 @@ tof31(right,i,iadc) = tof31(right,i,iadc)*cos(theta13) xkorr = atten(right,31,i,yhelp) c write(80+i,*) yhelp,xkorr - if (iz.le.1) xkorr=xkorr/hepratio + xkorr=xkorr/hepratio adctof_c(ch31b(i),hb31b(i))=tof31(right,i,iadc)/xkorr endif ENDIF @@ -1190,7 +1193,7 @@ tof32(left,i,iadc) = tof32(left,i,iadc)*cos(theta13) xkorr = atten(left,32,i,xhelp) c write(90+i,*) xhelp,xkorr - if (iz.le.1) xkorr=xkorr/hepratio + xkorr=xkorr/hepratio adctof_c(ch32a(i),hb32a(i))=tof32(left,i,iadc)/xkorr endif @@ -1199,7 +1202,7 @@ tof32(right,i,iadc) = tof32(right,i,iadc)*cos(theta13) xkorr = atten(right,32,i,xhelp) c write(90+i,*) xhelp,xkorr - if (iz.le.1) xkorr=xkorr/hepratio + xkorr=xkorr/hepratio adctof_c(ch32b(i),hb32b(i))=tof32(right,i,iadc)/xkorr endif ENDIF @@ -1584,13 +1587,14 @@ C if (icount.gt.0) beta_mean=sxw/sw C betatof_a(13) = beta_mean C + C-------- New mean beta calculation ----------------------- do i=1,12 btemp(i) = betatof_a(i) enddo - betatof_a(13)=newbeta(btemp,hitvec,10.,10.,20.) + betatof_a(13)=newbeta(1,btemp,hitvec,10.,10.,20.) C-------------------------------------------------------------- C write(*,*) betatof_a @@ -2113,7 +2117,7 @@ C**************************************************************************** C**************************************************************************** - function newbeta(b,hitvec,resmax,qualitycut,chi2cut) + function newbeta(iflag,b,hitvec,resmax,qualitycut,chi2cut) include 'input_tof.txt' include 'output_tof.txt' @@ -2123,8 +2127,9 @@ REAL resmax,qualitycut,chi2cut REAL w_i(12),w_il(6),quality,res,betachi,beta_mean_inv REAL sw,sxw,b(12),beta_mean,chi2,xhelp + REAL tdcfl(4,12) - INTEGER icount,hitvec(6) + INTEGER iflag,icount,hitvec(6) INTEGER itop(12),ibot(12) DATA itop /1,1,2,2,3,3,4,4,1,1,2,2/ @@ -2139,9 +2144,22 @@ tof31_i = hitvec(5) tof32_i = hitvec(6) -c write(*,*) '------------ In NEWBETA ----------------' -c write(*,*) hitvec -c write(*,*) b + if (iflag.eq.1) then ! call from tofl2com + do i=1,4 + do j=1,12 + tdcfl(i,j) = tdcflagtof(i,j) + enddo + enddo + endif + + if (iflag.eq.2) then ! call from toftrk + do i=1,4 + do j=1,12 + tdcfl(i,j) = tdcflag(i,j) + enddo + enddo + endif + C--- Find out ToF layers with artificial TDC values ------------- @@ -2149,15 +2167,13 @@ w_il(jj) = 1000. enddo -C write(*,*) tdcflagtof if (tof11_i.gt.0) then if ((tofmask(ch11a(tof11_i),hb11a(tof11_i)).gt.0).or. & (tofmask(ch11b(tof11_i),hb11b(tof11_i)).gt.0)) then w_il(1)=0 - i1=tdcflagtof(ch11a(tof11_i),hb11a(tof11_i)) - i2=tdcflagtof(ch11b(tof11_i),hb11b(tof11_i)) -c write(*,*) '11 ',i1,i2 + i1=tdcfl(ch11a(tof11_i),hb11a(tof11_i)) + i2=tdcfl(ch11b(tof11_i),hb11b(tof11_i)) if ((i1.eq.1).or.(i2.eq.1)) w_il(1) = 1 ! tdcflag endif endif @@ -2166,9 +2182,8 @@ if ((tofmask(ch12a(tof12_i),hb12a(tof12_i)).gt.0).or. & (tofmask(ch12b(tof12_i),hb12b(tof12_i)).gt.0)) then w_il(2)=0 - i1=tdcflagtof(ch12a(tof12_i),hb12a(tof12_i)) - i2=tdcflagtof(ch12b(tof12_i),hb12b(tof12_i)) -c write(*,*) '12 ',i1,i2 + i1=tdcfl(ch12a(tof12_i),hb12a(tof12_i)) + i2=tdcfl(ch12b(tof12_i),hb12b(tof12_i)) if ((i1.eq.1).or.(i2.eq.1)) w_il(2) = 1 ! tdcflag endif endif @@ -2177,9 +2192,8 @@ if ((tofmask(ch21a(tof21_i),hb21a(tof21_i)).gt.0).or. & (tofmask(ch21b(tof21_i),hb21b(tof21_i)).gt.0)) then w_il(3)=0 - i1=tdcflagtof(ch21a(tof21_i),hb21a(tof21_i)) - i2=tdcflagtof(ch21b(tof21_i),hb21b(tof21_i)) -c write(*,*) '21 ',i1,i2 + i1=tdcfl(ch21a(tof21_i),hb21a(tof21_i)) + i2=tdcfl(ch21b(tof21_i),hb21b(tof21_i)) if ((i1.eq.1).or.(i2.eq.1)) w_il(3) = 1 ! tdcflag endif endif @@ -2188,9 +2202,8 @@ if ((tofmask(ch22a(tof22_i),hb22a(tof22_i)).gt.0).or. & (tofmask(ch22b(tof22_i),hb22b(tof22_i)).gt.0)) then w_il(4)=0 - i1=tdcflagtof(ch22a(tof22_i),hb22a(tof22_i)) - i2=tdcflagtof(ch22b(tof22_i),hb22b(tof22_i)) -c write(*,*) '22 ',i1,i2 + i1=tdcfl(ch22a(tof22_i),hb22a(tof22_i)) + i2=tdcfl(ch22b(tof22_i),hb22b(tof22_i)) if ((i1.eq.1).or.(i2.eq.1)) w_il(4) = 1 ! tdcflag endif endif @@ -2199,9 +2212,8 @@ if ((tofmask(ch31a(tof31_i),hb11a(tof31_i)).gt.0).or. & (tofmask(ch31b(tof31_i),hb31b(tof31_i)).gt.0)) then w_il(5)=0 - i1=tdcflagtof(ch31a(tof31_i),hb31a(tof31_i)) - i2=tdcflagtof(ch31b(tof31_i),hb31b(tof31_i)) -c write(*,*) '31 ',i1,i2 + i1=tdcfl(ch31a(tof31_i),hb31a(tof31_i)) + i2=tdcfl(ch31b(tof31_i),hb31b(tof31_i)) if ((i1.eq.1).or.(i2.eq.1)) w_il(5) = 1 ! tdcflag endif endif @@ -2210,14 +2222,12 @@ if ((tofmask(ch32a(tof32_i),hb32a(tof32_i)).gt.0).or. & (tofmask(ch32b(tof32_i),hb32b(tof32_i)).gt.0)) then w_il(6)=0 - i1=tdcflagtof(ch32a(tof32_i),hb32a(tof32_i)) - i2=tdcflagtof(ch32b(tof32_i),hb32b(tof32_i)) -c write(*,*) '32 ',i1,i2 + i1=tdcfl(ch32a(tof32_i),hb32a(tof32_i)) + i2=tdcfl(ch32b(tof32_i),hb32b(tof32_i)) if ((i1.eq.1).or.(i2.eq.1)) w_il(6) = 1 ! tdcflag endif endif -c write(*,*) w_il C------------------------------------------------------------------------ C--- Set weights for the 12 measurements using information for top and bottom: C--- if no measurements: weight = set to very high value=> not used @@ -2237,8 +2247,6 @@ w_i(jj) = 1./xhelp ENDDO -c write(*,*) w_i - C======================================================================== C--- Calculate mean beta for the first time ----------------------------- C--- We are using "1/beta" since its error is gaussian ------------------ @@ -2259,7 +2267,6 @@ if (icount.gt.0) beta_mean=1./(sxw/sw) beta_mean_inv = 1./beta_mean -c write(*,*) icount,beta_mean C--- Calculate beta for the second time, use residuals of the single C--- measurements to get a chi2 value @@ -2275,7 +2282,6 @@ IF ((abs(1./b(jj)).gt.0.1).and.(abs(1./b(jj)).lt.15.) & .and.(w_i(jj).GT.0.01)) THEN res = beta_mean_inv - (1./b(jj)) ; -C write(*,*) jj,abs(res*w_i(jj)) if (abs(res*w_i(jj)).lt.resmax) THEN chi2 = chi2 + (res*w_i(jj))**2. icount = icount+1 @@ -2296,7 +2302,6 @@ beta_mean=100. if ((chi2.lt.chi2cut).and.(quality.gt.qualitycut)) & beta_mean = betachi -c write(*,*) icount,chi2,quality,beta_mean newbeta = beta_mean END @@ -2304,4 +2309,3 @@ C**************************************************************************** -