--- DarthVader/ToFLevel2/src/tofl2com.for 2006/08/10 06:32:04 1.3 +++ DarthVader/ToFLevel2/src/tofl2com.for 2010/02/17 11:50:53 1.14 @@ -1,4 +1,40 @@ -***************************************************************************** + +C****************************************************************************** +C +C 08-12-06 WM: adc_c-bug : The raw ADc value was multiplied with cos(theta) +C and AFTER that there was an if statement "if tof32(right,i,iadc) < 4095" +C +C jan-07 GF: ADCflags(4,12) inserted to flag artificial ADC values +C jan-07 WM: artificial ADC values created using attenuation calibration +C jan-07 WM: modified xtofpos flag "101". xtofpos must be inside physical +C dimension of the paddle +/- 10 cm +C jan-07 WM: if xtofpos=101 then this paddle is not used for beta +C calculation +C jan-07 WM: the definition for a "hit" is changed: Now we must have a +C valid TDC signal on both sides +C jan-07 WM: flag for PMTs #10 and #35 added, TDC=819 due to bit-shift +C jan-07 WM: bug removed: in some cases tdc_tw was calculated due to a +C leftover "xhelp" value +C apr-07 WM: attenuation fit curve is now a double exponential fit +C conversion from raw ADC to pC using calibration function +C variables xtr_tof and ytr_tof inserted (filled with default) +C jan-08 WM: Major Update: Time Walk correction introduced +C Additionalyl we use the information from the "check_charge" +C function to fill artificial ADC values and make small corrections +C to the k1-parameter (for Z>2) +C feb-08 WM: Calculation of beta(13) changed: First a mean beta is calculated, +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 oct-08 WM: Calculation of zenith angle debugged, sometimes strange values +C were possible +C nov-09 WM: the dEdx part ("adctof_c") moved to the new dEdx routine from Napoli +C feb-10 WM: k1 values now for Z=1, Z=2, Z>2, k2 values are fix +C****************************************************************************** + INTEGER FUNCTION TOFL2COM() c IMPLICIT NONE @@ -13,20 +49,20 @@ LOGICAL check REAL secure - INTEGER j - REAL xhelp_a,xhelp_t + INTEGER j,hitvec(6) REAL dx,dy,dr,ds - REAL yhelp,xdummy,xkorr0,xhelp,xhelp1,xhelp2 - REAL c1,c2,sw,sxw,w_i - INTEGER icount + REAL yhelp,yhelp1,yhelp2,xhelp,xhelp1,xhelp2 + REAL c1,c2 + REAL dist + +C REAL sw,sxw,w_i +C INTEGER icount +C REAL beta_mean INTEGER tof11_j,tof21_j,tof31_j INTEGER tof12_j,tof22_j,tof32_j - REAL beta_mean - - c value for status of each PM-data c first index : 1 = left, 2 = right c second index : 1... number of paddle @@ -34,19 +70,66 @@ INTEGER tof21_event(2,2),tof22_event(2,2) INTEGER tof31_event(2,3),tof32_event(2,3) + + REAL y_coor_lin11c(8,2),x_coor_lin12c(6,2) + REAL x_coor_lin21c(2,2),y_coor_lin22c(2,2) + REAL y_coor_lin31c(3,2),x_coor_lin32c(3,2) + + DATA y_coor_lin11c(1,1),y_coor_lin11c(1,2) /-20.66,-2.497/ + DATA y_coor_lin11c(2,1),y_coor_lin11c(2,2) /-9.10, -2.52/ + DATA y_coor_lin11c(3,1),y_coor_lin11c(3,2) /-24.07,-2.12/ + DATA y_coor_lin11c(4,1),y_coor_lin11c(4,2) /-13.40,-2.47/ + DATA y_coor_lin11c(5,1),y_coor_lin11c(5,2) /-31.07,-2.32/ + DATA y_coor_lin11c(6,1),y_coor_lin11c(6,2) /-21.69,-2.63/ + DATA y_coor_lin11c(7,1),y_coor_lin11c(7,2) /-12.37,-2.65/ + DATA y_coor_lin11c(8,1),y_coor_lin11c(8,2) /-10.81,-3.15/ + + DATA x_coor_lin12c(1,1),x_coor_lin12c(1,2) /12.96, -2.65/ + DATA x_coor_lin12c(2,1),x_coor_lin12c(2,2) /17.12,-2.44/ + DATA x_coor_lin12c(3,1),x_coor_lin12c(3,2) /7.26, -1.98/ + DATA x_coor_lin12c(4,1),x_coor_lin12c(4,2) /-22.52,-2.27/ + DATA x_coor_lin12c(5,1),x_coor_lin12c(5,2) /-18.54,-2.28/ + DATA x_coor_lin12c(6,1),x_coor_lin12c(6,2) /-7.67,-2.15/ + + DATA x_coor_lin21c(1,1),x_coor_lin21c(1,2) /22.56,-1.56/ + DATA x_coor_lin21c(2,1),x_coor_lin21c(2,2) /13.94,-1.56/ + + DATA y_coor_lin22c(1,1),y_coor_lin22c(1,2) /-24.24,-2.23/ + DATA y_coor_lin22c(2,1),y_coor_lin22c(2,2) /-45.99,-1.68/ + + DATA y_coor_lin31c(1,1),y_coor_lin31c(1,2) /-22.99,-3.54/ + DATA y_coor_lin31c(2,1),y_coor_lin31c(2,2) /-42.28,-4.10/ + DATA y_coor_lin31c(3,1),y_coor_lin31c(3,2) /-41.29,-3.69/ + + DATA x_coor_lin32c(1,1),x_coor_lin32c(1,2) /0.961, -3.22/ + DATA x_coor_lin32c(2,1),x_coor_lin32c(2,2) /4.98,-3.48/ + DATA x_coor_lin32c(3,1),x_coor_lin32c(3,2) /-22.08,-3.37/ + - REAL theta12,theta13,theta23 -C-- DATA ZTOF/53.74,53.04,23.94,23.44,-23.49,-24.34/ !Sergio 9.05.2006 + REAL theta13 + + DOUBLE PRECISION ZTOF(6) + DATA ZTOF/53.74,53.04,23.94,23.44,-23.49,-24.34/ !Sergio 9.05.2006 + REAL tofarm12 PARAMETER (tofarm12 = 29.70) ! from 53.39 to 23.69 - REAL tofarm23 + REAL tofarm23 PARAMETER (tofarm23 = 47.61) ! from 23.69 to -23.92 REAL tofarm13 PARAMETER (tofarm13 = 77.31) ! from 53.39 to -23.92 - + + REAL hepratio INTEGER ihelp - REAL xkorr + REAL xkorr,btemp(12) + + REAL atten,pc_adc,check_charge,newbeta + + INTEGER IZ + + + INTEGER ifst + DATA ifst /0/ C--------------------------------------- C @@ -56,28 +139,42 @@ C C CALCULATE COMMON VARIABLES C +C------------------------------------------------------------------- -******************************************************************* - icounter = icounter + 1 + if (ifst.eq.0) then + + ifst=1 + +C amplitude has to be 'secure' higher than pedestal for an adc event + secure = 2. -* amplitude has to be 'secure' higher than pedestal for an adc event - secure = 2. +C ratio between helium and proton ca. 4 + hepratio = 4. ! + offset = 1 + slope = 2 + left = 1 + right = 2 + none_ev = 0 + none_find = 0 + tdc_ev = 1 + adc_ev = 1 + itdc = 1 + iadc = 2 + + ENDIF +C--------------------------------------------------------------------- + + icounter = icounter + 1 - offset = 1 - slope = 2 - left = 1 - right = 2 - none_ev = 0 - none_find = 0 - tdc_ev = 1 - adc_ev = 1 - itdc = 1 - iadc = 2 do i=1,13 betatof_a(i) = 100. ! As in "troftrk.for" enddo + do i=1,6 + hitvec(i) = -1 + enddo + do i=1,4 do j=1,12 adctof_c(i,j) = 1000. @@ -99,55 +196,76 @@ enddo -c the calibration files are read in the main program from xxx_tofcalib.rz +c gf adc falg: + do i=1,4 + do j=1,12 + adcflagtof(i,j) = 0 + enddo + enddo + +c gf tdc falg: + do i=1,4 + do j=1,12 + tdcflagtof(i,j) = 0 + enddo + enddo + + +C--- Fill xtr_tof and ytr_tof: positions from tracker at ToF layers +C--- since this is standalone ToF fill with default values + do j=1,6 + xtr_tof(j) = 101. + ytr_tof(j) = 101. + enddo +c the calibration files are read in the main program from xxx_tofcalib.rz c-------------------------get ToF data -------------------------------- c put the adc and tdc values from ntuple into tofxx(i,j,k) variables - +c adc valueas are then pC do j=1,8 - tof11(1,j,2) = adc(ch11a(j),hb11a(j)) - tof11(2,j,2) = adc(ch11b(j),hb11b(j)) - tof11(1,j,1) = tdc(ch11a(j),hb11a(j)) - tof11(2,j,1) = tdc(ch11b(j),hb11b(j)) + tof11(1,j,2) = pc_adc(adc(ch11a(j),hb11a(j))) + tof11(2,j,2) = pc_adc(adc(ch11b(j),hb11b(j))) + tof11(1,j,1) = (tdc(ch11a(j),hb11a(j))) + tof11(2,j,1) = (tdc(ch11b(j),hb11b(j))) enddo do j=1,6 - tof12(1,j,2) = adc(ch12a(j),hb12a(j)) - tof12(2,j,2) = adc(ch12b(j),hb12b(j)) - tof12(1,j,1) = tdc(ch12a(j),hb12a(j)) - tof12(2,j,1) = tdc(ch12b(j),hb12b(j)) + tof12(1,j,2) = pc_adc(adc(ch12a(j),hb12a(j))) + tof12(2,j,2) = pc_adc(adc(ch12b(j),hb12b(j))) + tof12(1,j,1) = (tdc(ch12a(j),hb12a(j))) + tof12(2,j,1) = (tdc(ch12b(j),hb12b(j))) enddo do j=1,2 - tof21(1,j,2) = adc(ch21a(j),hb21a(j)) - tof21(2,j,2) = adc(ch21b(j),hb21b(j)) - tof21(1,j,1) = tdc(ch21a(j),hb21a(j)) - tof21(2,j,1) = tdc(ch21b(j),hb21b(j)) + tof21(1,j,2) = pc_adc(adc(ch21a(j),hb21a(j))) + tof21(2,j,2) = pc_adc(adc(ch21b(j),hb21b(j))) + tof21(1,j,1) = (tdc(ch21a(j),hb21a(j))) + tof21(2,j,1) = (tdc(ch21b(j),hb21b(j))) enddo do j=1,2 - tof22(1,j,2) = adc(ch22a(j),hb22a(j)) - tof22(2,j,2) = adc(ch22b(j),hb22b(j)) - tof22(1,j,1) = tdc(ch22a(j),hb22a(j)) - tof22(2,j,1) = tdc(ch22b(j),hb22b(j)) + tof22(1,j,2) = pc_adc(adc(ch22a(j),hb22a(j))) + tof22(2,j,2) = pc_adc(adc(ch22b(j),hb22b(j))) + tof22(1,j,1) = (tdc(ch22a(j),hb22a(j))) + tof22(2,j,1) = (tdc(ch22b(j),hb22b(j))) enddo do j=1,3 - tof31(1,j,2) = adc(ch31a(j),hb31a(j)) - tof31(2,j,2) = adc(ch31b(j),hb31b(j)) - tof31(1,j,1) = tdc(ch31a(j),hb31a(j)) - tof31(2,j,1) = tdc(ch31b(j),hb31b(j)) + tof31(1,j,2) = pc_adc(adc(ch31a(j),hb31a(j))) + tof31(2,j,2) = pc_adc(adc(ch31b(j),hb31b(j))) + tof31(1,j,1) = (tdc(ch31a(j),hb31a(j))) + tof31(2,j,1) = (tdc(ch31b(j),hb31b(j))) enddo do j=1,3 - tof32(1,j,2) = adc(ch32a(j),hb32a(j)) - tof32(2,j,2) = adc(ch32b(j),hb32b(j)) - tof32(1,j,1) = tdc(ch32a(j),hb32a(j)) - tof32(2,j,1) = tdc(ch32b(j),hb32b(j)) + tof32(1,j,2) = pc_adc(adc(ch32a(j),hb32a(j))) + tof32(2,j,2) = pc_adc(adc(ch32b(j),hb32b(j))) + tof32(1,j,1) = (tdc(ch32a(j),hb32a(j))) + tof32(2,j,1) = (tdc(ch32b(j),hb32b(j))) enddo C---------------------------------------------------------------------- @@ -195,8 +313,66 @@ if (abs(tof32(2,i,iadc)).gt.10000.) tof32(2,i,iadc)= 10000. ENDDO +C---------------------------------------------------------------------- +C------------------ set ADC & TDC flag = 0 ------------------------ +C---------------------------------------------------------------------- + + do j=1,8 + if (adc(ch11a(j),hb11a(j)).LT.4096)adcflagtof(ch11a(j),hb11a(j))=0 + if (adc(ch11b(j),hb11b(j)).LT.4096)adcflagtof(ch11b(j),hb11b(j))=0 + if (tdc(ch11a(j),hb11a(j)).LT.4096)tdcflagtof(ch11a(j),hb11a(j))=0 + if (tdc(ch11b(j),hb11b(j)).LT.4096)tdcflagtof(ch11b(j),hb11b(j))=0 + enddo + do j=1,6 + if (adc(ch12a(j),hb12a(j)).LT.4096)adcflagtof(ch12a(j),hb12a(j))=0 + if (adc(ch12b(j),hb12b(j)).LT.4096)adcflagtof(ch12b(j),hb12b(j))=0 + if (tdc(ch12a(j),hb12a(j)).LT.4096)tdcflagtof(ch12a(j),hb12a(j))=0 + if (tdc(ch12b(j),hb12b(j)).LT.4096)tdcflagtof(ch12b(j),hb12b(j))=0 + enddo + do j=1,2 + if (adc(ch21a(j),hb21a(j)).LT.4096)adcflagtof(ch21a(j),hb21a(j))=0 + if (adc(ch21b(j),hb21b(j)).LT.4096)adcflagtof(ch21b(j),hb21b(j))=0 + if (tdc(ch21a(j),hb21a(j)).LT.4096)tdcflagtof(ch21a(j),hb21a(j))=0 + if (tdc(ch21b(j),hb21b(j)).LT.4096)tdcflagtof(ch21b(j),hb21b(j))=0 + enddo + do j=1,2 + if (adc(ch22a(j),hb22a(j)).LT.4096)adcflagtof(ch22a(j),hb22a(j))=0 + if (adc(ch22b(j),hb22b(j)).LT.4096)adcflagtof(ch22b(j),hb22b(j))=0 + if (tdc(ch22a(j),hb22a(j)).LT.4096)tdcflagtof(ch22a(j),hb22a(j))=0 + if (tdc(ch22b(j),hb22b(j)).LT.4096)tdcflagtof(ch22b(j),hb22b(j))=0 + enddo + do j=1,3 + if (adc(ch31a(j),hb31a(j)).LT.4096)adcflagtof(ch31a(j),hb31a(j))=0 + if (adc(ch31b(j),hb31b(j)).LT.4096)adcflagtof(ch31b(j),hb31b(j))=0 + if (tdc(ch31a(j),hb31a(j)).LT.4096)tdcflagtof(ch31a(j),hb31a(j))=0 + if (tdc(ch31b(j),hb31b(j)).LT.4096)tdcflagtof(ch31b(j),hb31b(j))=0 + enddo + do j=1,3 + if (adc(ch32a(j),hb32a(j)).LT.4096)adcflagtof(ch32a(j),hb32a(j))=0 + if (adc(ch32b(j),hb32b(j)).LT.4096)adcflagtof(ch32b(j),hb32b(j))=0 + if (tdc(ch32a(j),hb32a(j)).LT.4096)tdcflagtof(ch32a(j),hb32a(j))=0 + if (tdc(ch32b(j),hb32b(j)).LT.4096)tdcflagtof(ch32b(j),hb32b(j))=0 + enddo + C---------------------------------------------------------------- -C------------Check Paddles for hits ----------------------- +C---------- Check PMTs 10 and 35 for strange TDC values---------- +C---------------------------------------------------------------- + +C---- S116A TDC=819 + if (tof11(1,6,1).EQ.819) then + tof11(1,6,1) = 4095 + tdcflagtof(ch11a(6),hb11a(6))=2 + endif + +C---- S222B TDC=819 + if (tof22(2,2,1).EQ.819) then + tof22(2,2,1) = 4095 + tdcflagtof(ch22b(2),hb22b(2))=2 + endif + +C---------------------------------------------------------------- +C------------ Check Paddles for hits ----------------------- +C------ a "hit" means TDC values<4095 on both sides ------------ C---------------------------------------------------------------- C upper tof S11 @@ -206,9 +382,6 @@ tof11_event(j,i) = none_ev IF ((tof11(j,i,itdc).LT.2000).AND.(tof11(j,i,itdc).GT.100)) + tof11_event(j,i) = tof11_event(j,i) + tdc_ev - IF ((tof11(j,i,iadc).GT.secure).AND. - + (tof11(j,i,iadc).LT.4095)) - + tof11_event(j,i) = tof11_event(j,i) + adc_ev ENDDO ENDDO @@ -239,9 +412,6 @@ tof12_event(j,i) = none_ev IF ((tof12(j,i,itdc).LT.2000).AND.(tof12(j,i,itdc).GT.100)) + tof12_event(j,i) = tof12_event(j,i) + tdc_ev - IF ((tof12(j,i,iadc).GT.secure).AND. - + (tof12(j,i,iadc).LT.4095)) - + tof12_event(j,i) = tof12_event(j,i) + adc_ev ENDDO ENDDO @@ -272,9 +442,6 @@ tof21_event(j,i) = none_ev IF ((tof21(j,i,itdc).LT.2000).AND.(tof21(j,i,itdc).GT.100)) + tof21_event(j,i) = tof21_event(j,i) + tdc_ev - IF ((tof21(j,i,iadc).GT.secure).AND. - + (tof21(j,i,iadc).LT.4095)) - + tof21_event(j,i) = tof21_event(j,i) + adc_ev ENDDO ENDDO @@ -304,9 +471,6 @@ tof22_event(j,i) = none_ev IF ((tof22(j,i,itdc).LT.2000).AND.(tof22(j,i,itdc).GT.100)) + tof22_event(j,i) = tof22_event(j,i) + tdc_ev - IF ((tof22(j,i,iadc).GT.secure).AND. - + (tof22(j,i,iadc).LT.4095)) - + tof22_event(j,i) = tof22_event(j,i) + adc_ev ENDDO ENDDO @@ -337,9 +501,6 @@ tof31_event(j,i) = none_ev IF ((tof31(j,i,itdc).LT.2000).AND.(tof31(j,i,itdc).GT.100)) + tof31_event(j,i) = tof31_event(j,i) + tdc_ev - IF ((tof31(j,i,iadc).GT.secure).AND. - + (tof31(j,i,iadc).LT.4095)) - + tof31_event(j,i) = tof31_event(j,i) + adc_ev ENDDO ENDDO @@ -369,9 +530,6 @@ tof32_event(j,i) = none_ev IF ((tof32(j,i,itdc).LT.2000).AND.(tof32(j,i,itdc).GT.100)) + tof32_event(j,i) = tof32_event(j,i) + tdc_ev - IF ((tof32(j,i,iadc).GT.secure).AND. - + (tof32(j,i,iadc).LT.4095)) - + tof32_event(j,i) = tof32_event(j,i) + adc_ev ENDDO ENDDO @@ -393,130 +551,58 @@ ENDIF ENDIF ENDIF - ENDDO + ENDDO - do i=1,6 + do i=1,6 tof_i_flag(i)=0 tof_j_flag(i)=0 - enddo - - tof_i_flag(1)=tof11_i - tof_i_flag(2)=tof12_i - tof_i_flag(3)=tof21_i - tof_i_flag(4)=tof22_i - tof_i_flag(5)=tof31_i - tof_i_flag(6)=tof32_i - - tof_j_flag(1)=tof11_j - tof_j_flag(2)=tof12_j - tof_j_flag(3)=tof21_j - tof_j_flag(4)=tof22_j - tof_j_flag(5)=tof31_j - tof_j_flag(6)=tof32_j + enddo + + tof_i_flag(1)=tof11_i + tof_i_flag(2)=tof12_i + tof_i_flag(3)=tof21_i + tof_i_flag(4)=tof22_i + tof_i_flag(5)=tof31_i + tof_i_flag(6)=tof32_i + + tof_j_flag(1)=tof11_j + tof_j_flag(2)=tof12_j + tof_j_flag(3)=tof21_j + tof_j_flag(4)=tof22_j + tof_j_flag(5)=tof31_j + tof_j_flag(6)=tof32_j + + hitvec(1)=tof11_i + hitvec(2)=tof12_i + hitvec(3)=tof21_i + hitvec(4)=tof22_i + hitvec(5)=tof31_i + hitvec(6)=tof32_i -C-------------------------------------------------------------------- -C--------------------Time walk correction ------------------------- -C-------------------------------------------------------------------- - - DO i=1,8 - xhelp_a = tof11(left,i,iadc) - xhelp_t = tof11(left,i,itdc) - if(xhelp_a>0) xhelp = tw11(left,i)/sqrt(xhelp_a) - tof11(left,i,itdc) = xhelp_t + xhelp - tdc_c(ch11a(i),hb11a(i))=tof11(left,i,itdc) - xhelp_a = tof11(right,i,iadc) - xhelp_t = tof11(right,i,itdc) - if(xhelp_a>0) xhelp = tw11(right,i)/sqrt(xhelp_a) - tof11(right,i,itdc) = xhelp_t + xhelp - tdc_c(ch11b(i),hb11b(i))=tof11(right,i,itdc) - ENDDO - - DO i=1,6 - xhelp_a = tof12(left,i,iadc) - xhelp_t = tof12(left,i,itdc) - if(xhelp_a>0) xhelp = tw12(left,i)/sqrt(xhelp_a) - tof12(left,i,itdc) = xhelp_t + xhelp - tdc_c(ch12a(i),hb12a(i))=tof12(left,i,itdc) - xhelp_a = tof12(right,i,iadc) - xhelp_t = tof12(right,i,itdc) - if(xhelp_a>0) xhelp = tw12(right,i)/sqrt(xhelp_a) - tof12(right,i,itdc) = xhelp_t + xhelp - tdc_c(ch12b(i),hb12b(i))=tof12(right,i,itdc) - ENDDO -C---- - DO i=1,2 - xhelp_a = tof21(left,i,iadc) - xhelp_t = tof21(left,i,itdc) - if(xhelp_a>0) xhelp = tw21(left,i)/sqrt(xhelp_a) - tof21(left,i,itdc) = xhelp_t + xhelp - tdc_c(ch21a(i),hb21a(i))=tof21(left,i,itdc) - xhelp_a = tof21(right,i,iadc) - xhelp_t = tof21(right,i,itdc) - if(xhelp_a>0) xhelp = tw21(right,i)/sqrt(xhelp_a) - tof21(right,i,itdc) = xhelp_t + xhelp - tdc_c(ch21b(i),hb21b(i))=tof21(right,i,itdc) - ENDDO - - DO i=1,2 - xhelp_a = tof22(left,i,iadc) - xhelp_t = tof22(left,i,itdc) - if(xhelp_a>0) xhelp = tw22(left,i)/sqrt(xhelp_a) - tof22(left,i,itdc) = xhelp_t + xhelp - tdc_c(ch22a(i),hb22a(i))=tof22(left,i,itdc) - xhelp_a = tof22(right,i,iadc) - xhelp_t = tof22(right,i,itdc) - if(xhelp_a>0) xhelp = tw22(right,i)/sqrt(xhelp_a) - tof22(right,i,itdc) = xhelp_t + xhelp - tdc_c(ch22b(i),hb22b(i))=tof22(right,i,itdc) - ENDDO -C---- - - DO i=1,3 - xhelp_a = tof31(left,i,iadc) - xhelp_t = tof31(left,i,itdc) - if(xhelp_a>0) xhelp = tw31(left,i)/sqrt(xhelp_a) - tof31(left,i,itdc) = xhelp_t + xhelp - tdc_c(ch31a(i),hb31a(i))=tof31(left,i,itdc) - xhelp_a = tof31(right,i,iadc) - xhelp_t = tof31(right,i,itdc) - if(xhelp_a>0) xhelp = tw31(right,i)/sqrt(xhelp_a) - tof31(right,i,itdc) = xhelp_t + xhelp - tdc_c(ch31b(i),hb31b(i))=tof31(right,i,itdc) - ENDDO - - DO i=1,3 - xhelp_a = tof32(left,i,iadc) - xhelp_t = tof32(left,i,itdc) - if(xhelp_a>0) xhelp = tw32(left,i)/sqrt(xhelp_a) - tof32(left,i,itdc) = xhelp_t + xhelp - tdc_c(ch32a(i),hb32a(i))=tof32(left,i,itdc) - xhelp_a = tof32(right,i,iadc) - xhelp_t = tof32(right,i,itdc) - if(xhelp_a>0) xhelp = tw32(right,i)/sqrt(xhelp_a) - tof32(right,i,itdc) = xhelp_t + xhelp - tdc_c(ch32b(i),hb32b(i))=tof32(right,i,itdc) - ENDDO -C---- - C------------------------------------------------------------------ -C--- calculate track position in paddle using timing difference +C-- calculate track position in paddle using timing difference +C-- this calculation is preliminary and uses some standard +C-- calibration values, but we need to find a rough position to +C-- be able to calculate artificial ADC values (needed for the +C-- timewalk... C------------------------------------------------------------------ - do i=1,3 + do i=1,3 xtofpos(i)=100. ytofpos(i)=100. - enddo + enddo + C-----------------------------S1 -------------------------------- IF (tof11_i.GT.none_find) THEN ytofpos(1) = ((tof11(1,tof11_i,itdc)-tof11(2,tof11_i,itdc))/2. - + -y_coor_lin11(tof11_i,offset))/y_coor_lin11(tof11_i,slope) + + -y_coor_lin11c(tof11_i,offset))/y_coor_lin11c(tof11_i,slope) endif IF (tof12_i.GT.none_find) THEN xtofpos(1) = ((tof12(1,tof12_i,itdc)-tof12(2,tof12_i,itdc))/2. - + -x_coor_lin12(tof12_i,offset))/x_coor_lin12(tof12_i,slope) + + -x_coor_lin12c(tof12_i,offset))/x_coor_lin12c(tof12_i,slope) endif @@ -524,12 +610,12 @@ IF (tof21_i.GT.none_find) THEN xtofpos(2) = ((tof21(1,tof21_i,itdc)-tof21(2,tof21_i,itdc))/2. - + -x_coor_lin21(tof21_i,offset))/x_coor_lin21(tof21_i,slope) + + -x_coor_lin21c(tof21_i,offset))/x_coor_lin21c(tof21_i,slope) endif IF (tof22_i.GT.none_find) THEN ytofpos(2) = ((tof22(1,tof22_i,itdc)-tof22(2,tof22_i,itdc))/2. - + -y_coor_lin22(tof22_i,offset))/y_coor_lin22(tof22_i,slope) + + -y_coor_lin22c(tof22_i,offset))/y_coor_lin22c(tof22_i,slope) endif @@ -537,220 +623,446 @@ IF (tof31_i.GT.none_find) THEN ytofpos(3) = ((tof31(1,tof31_i,itdc)-tof31(2,tof31_i,itdc))/2. - + -y_coor_lin31(tof31_i,offset))/y_coor_lin31(tof31_i,slope) + + -y_coor_lin31c(tof31_i,offset))/y_coor_lin31c(tof31_i,slope) endif IF (tof32_i.GT.none_find) THEN xtofpos(3) = ((tof32(1,tof32_i,itdc)-tof32(2,tof32_i,itdc))/2. - + -x_coor_lin32(tof32_i,offset))/x_coor_lin32(tof32_i,slope) + + -x_coor_lin32c(tof32_i,offset))/x_coor_lin32c(tof32_i,slope) endif - do i=1,3 - if (abs(xtofpos(i)).gt.100.) then - xtofpos(i)=101. - endif - if (abs(ytofpos(i)).gt.100.) then - ytofpos(i)=101. - endif - enddo - C---------------------------------------------------------------------- -C--------------------Corrections on ADC-data ------------------------- -C---------------------zenith angle theta --------------------------- +C--------------------- zenith angle theta --------------------------- C---------------------------------------------------------------------- - dx=0. - dy=0. - dr=0. - theta13 = 0. - - IF ((tof12_i.GT.none_find).AND.(tof32_i.GT.none_find)) - & dx = xtofpos(1) - xtofpos(3) - IF ((tof11_i.GT.none_find).AND.(tof31_i.GT.none_find)) - & dy = ytofpos(1) - ytofpos(3) - dr = sqrt(dx*dx+dy*dy) - theta13 = atan(dr/tofarm13) - - dx=0. - dy=0. - dr=0. - theta12 = 0. - - IF ((tof12_i.GT.none_find).AND.(tof21_i.GT.none_find)) - & dx = xtofpos(1) - xtofpos(2) - IF ((tof11_i.GT.none_find).AND.(tof22_i.GT.none_find)) - & dy = ytofpos(1) - ytofpos(2) - dr = sqrt(dx*dx+dy*dy) - theta12 = atan(dr/tofarm12) - - dx=0. - dy=0. - dr=0. - theta23 = 0. - - IF ((tof21_i.GT.none_find).AND.(tof32_i.GT.none_find)) - & dx = xtofpos(2) - xtofpos(3) - IF ((tof22_i.GT.none_find).AND.(tof31_i.GT.none_find)) - & dy = ytofpos(2) - ytofpos(3) - dr = sqrt(dx*dx+dy*dy) - theta23 = atan(dr/tofarm23) - - + xhelp1=0. + if (tof11_i.GT.none_find) xhelp1=tof11_x(tof11_i) + if (xtofpos(1).lt.100) xhelp1=xtofpos(1) + + yhelp1=0. + if (tof12_i.GT.none_find) yhelp1=tof12_y(tof12_i) + if (ytofpos(1).lt.100) yhelp1=ytofpos(1) + + + yhelp2=0. + if (tof32_i.GT.none_find) yhelp2=tof32_y(tof32_i) + if (ytofpos(3).lt.100) yhelp2=ytofpos(3) + + xhelp2=0. + if (tof31_i.GT.none_find) xhelp2=tof31_x(tof31_i) + if (xtofpos(3).lt.100) xhelp2=xtofpos(3) + + + dx=0. + dy=0. + dr=0. + theta13 = 0. + + dx = xhelp1 - xhelp2 + dy = yhelp1 - yhelp2 + dr = sqrt(dx*dx+dy*dy) + theta13 = atan(dr/tofarm13) + + C---------------------------------------------------------------------- -C------------------angle and ADC(x) correction +C--- check charge: +C--- if Z=2 we should use the attenuation curve for helium to +C--- fill the artificail ADC values and NOT divide by "hepratio" +C--- if Z>2 we should do a correction to +C--- the k1 constants in the beta calculation C---------------------------------------------------------------------- -C-----------------------------S1 -------------------------------- + + iz = int(check_charge(theta13,hitvec)) +C write(*,*) 'charge in tofl2com',iz + +C-------------------------------------------------------------------- +C---- if TDCleft.and.TDCright and NO ADC insert artificial ADC +C---- values +C-------------------------------------------------------------------- c middle y (or x) position of the upper and middle ToF-Paddle c DATA tof11_x/ -17.85,-12.75,-7.65,-2.55,2.55,7.65,12.75,17.85/ c DATA tof12_y/ -13.75,-8.25,-2.75,2.75,8.25,13.75/ -c DATA tof21_y/ -3.75,3.75/ +c DATA tof21_y/ 3.75,-3.75/ ! paddles in different order c DATA tof22_x/ -4.5,4.5/ c DATA tof31_x/ -6.0,0.,6.0/ c DATA tof32_y/ -5.0,0.0,5.0/ - yhelp=0. - if (tof12_i.GT.none_find) yhelp=tof12_y(tof12_i) - if (ytofpos(1).lt.100) yhelp=ytofpos(1) - IF (tof11_i.GT.none_find.AND.abs(yhelp).lt.100) THEN +C---------------------------- S1 ------------------------------------- +c yhelp=0. + yhelp=100. ! WM + if (tof12_i.GT.none_find) yhelp=tof12_y(tof12_i) + if (ytofpos(1).lt.100) yhelp=ytofpos(1) + + IF (tof11_i.GT.none_find.AND.abs(yhelp).lt.100) THEN i = tof11_i - xdummy=tof11(left,i,iadc) - tof11(left,i,iadc) = tof11(left,i,iadc)*cos(theta13) - if (tof11(left,i,iadc).lt.4095) then - xkorr=adcx11(left,i,1)*exp(-yhelp/adcx11(left,i,2)) - xkorr0=adcx11(left,i,1) - adctof_c(ch11a(i),hb11a(i))=tof11(left,i,iadc)/xkorr + if (adc(ch11a(i),hb11a(i)).eq.4095) then + xkorr = atten(left,11,i,yhelp) + if (iz.le.1) xkorr=xkorr/hepratio + tof11(left,i,iadc)=xkorr/cos(theta13) + adcflagtof(ch11a(i),hb11a(i)) = 1 endif - - tof11(right,i,iadc) = tof11(right,i,iadc)*cos(theta13) - if (tof11(right,i,iadc).lt.4095) then - xkorr=adcx11(right,i,1)*exp(yhelp/adcx11(right,i,2)) - xkorr0=adcx11(right,i,1) - adctof_c(ch11b(i),hb11b(i))=tof11(right,i,iadc)/xkorr + if (adc(ch11b(i),hb11b(i)).eq.4095) then + xkorr = atten(right,11,i,yhelp) + if (iz.le.1) xkorr=xkorr/hepratio + tof11(right,i,iadc)=xkorr/cos(theta13) + adcflagtof(ch11b(i),hb11b(i)) = 1 endif - ENDIF - - xhelp=0. - if (tof11_i.GT.none_find) xhelp=tof11_x(tof11_i) - if (xtofpos(1).lt.100) xhelp=xtofpos(1) + ENDIF - IF (tof12_i.GT.none_find.AND.abs(xhelp).lt.100) THEN +c xhelp=0. + xhelp=100. ! WM + if (tof11_i.GT.none_find) xhelp=tof11_x(tof11_i) + if (xtofpos(1).lt.100) xhelp=xtofpos(1) + IF (tof12_i.GT.none_find.AND.abs(xhelp).lt.100) THEN i = tof12_i - tof12(left,i,iadc) = tof12(left,i,iadc)*cos(theta13) - if (tof12(left,i,iadc).lt.4095) then - xkorr=adcx12(left,i,1)*exp(-xhelp/adcx12(left,i,2)) - xkorr0=adcx12(left,i,1) - adctof_c(ch12a(i),hb12a(i))=tof12(left,i,iadc)/xkorr + if (adc(ch12a(i),hb12a(i)).eq.4095) then + xkorr = atten(left,12,i,xhelp) + if (iz.le.1) xkorr=xkorr/hepratio + tof12(left,i,iadc) = xkorr/cos(theta13) + adcflagtof(ch12a(i),hb12a(i)) = 1 endif - - tof12(right,i,iadc) = tof12(right,i,iadc)*cos(theta13) - if (tof12(right,i,iadc).lt.4095) then - xkorr=adcx12(right,i,1)*exp(xhelp/adcx12(right,i,2)) - xkorr0=adcx12(right,i,1) - adctof_c(ch12b(i),hb12b(i))=tof12(right,i,iadc)/xkorr + if (adc(ch12b(i),hb12b(i)).eq.4095) then + xkorr = atten(right,12,i,xhelp) + if (iz.le.1) xkorr=xkorr/hepratio + tof12(right,i,iadc) = xkorr/cos(theta13) + adcflagtof(ch12b(i),hb12b(i)) = 1 endif - ENDIF + ENDIF C-----------------------------S2 -------------------------------- - xhelp=0. - if (tof22_i.GT.none_find) xhelp=tof22_x(tof22_i) - if (xtofpos(2).lt.100) xhelp=xtofpos(2) - - IF (tof21_i.GT.none_find.AND.abs(xhelp).lt.100) THEN +c xhelp=0. + xhelp=100. ! WM + if (tof22_i.GT.none_find) xhelp=tof22_x(tof22_i) + if (xtofpos(2).lt.100) xhelp=xtofpos(2) + IF (tof21_i.GT.none_find.AND.abs(xhelp).lt.100) THEN i = tof21_i - tof21(left,i,iadc) = tof21(left,i,iadc)*cos(theta13) - if (tof21(left,i,iadc).lt.4095) then - xkorr=adcx21(left,i,1)*exp(-xhelp/adcx21(left,i,2)) - xkorr0=adcx21(left,i,1) - adctof_c(ch21a(i),hb21a(i))=tof21(left,i,iadc)/xkorr + if (adc(ch21a(i),hb21a(i)).eq.4095) then + xkorr = atten(left,21,i,xhelp) + if (iz.le.1) xkorr=xkorr/hepratio + tof21(left,i,iadc) = xkorr/cos(theta13) + adcflagtof(ch21a(i),hb21a(i)) = 1 endif - - tof21(right,i,iadc) = tof21(right,i,iadc)*cos(theta13) - if (tof21(right,i,iadc).lt.4095) then - xkorr=adcx21(right,i,1)*exp(xhelp/adcx21(right,i,2)) - xkorr0=adcx21(right,i,1) - adctof_c(ch21b(i),hb21b(i))=tof21(right,i,iadc)/xkorr + if (adc(ch21b(i),hb21b(i)).eq.4095) then + xkorr = atten(right,21,i,xhelp) + if (iz.le.1) xkorr=xkorr/hepratio + tof21(right,i,iadc) = xkorr/cos(theta13) + adcflagtof(ch21b(i),hb21b(i)) = 1 endif - ENDIF + ENDIF - yhelp=0. - if (tof21_i.GT.none_find) yhelp=tof21_y(tof21_i) - if (ytofpos(2).lt.100) yhelp=ytofpos(2) - - IF (tof22_i.GT.none_find.AND.abs(yhelp).lt.100) THEN +c yhelp=0. + yhelp=100. ! WM + if (tof21_i.GT.none_find) yhelp=tof21_y(tof21_i) + if (ytofpos(2).lt.100) yhelp=ytofpos(2) + IF (tof22_i.GT.none_find.AND.abs(yhelp).lt.100) THEN i = tof22_i - tof22(left,i,iadc) = tof22(left,i,iadc)*cos(theta13) - if (tof22(left,i,iadc).lt.4095) then - xkorr=adcx22(left,i,1)*exp(-yhelp/adcx22(left,i,2)) - xkorr0=adcx22(left,i,1) - adctof_c(ch22a(i),hb22a(i))=tof22(left,i,iadc)/xkorr + if (adc(ch22a(i),hb22a(i)).eq.4095) then + xkorr = atten(left,22,i,yhelp) + if (iz.le.1) xkorr=xkorr/hepratio + tof22(left,i,iadc) = xkorr/cos(theta13) + adcflagtof(ch22a(i),hb22a(i)) = 1 endif - - tof22(right,i,iadc) = tof22(right,i,iadc)*cos(theta13) - if (tof22(right,i,iadc).lt.4095) then - xkorr=adcx22(right,i,1)*exp(yhelp/adcx22(right,i,2)) - xkorr0=adcx22(right,i,1) - adctof_c(ch22b(i),hb22b(i))=tof22(right,i,iadc)/xkorr + if (adc(ch22b(i),hb22b(i)).eq.4095) then + xkorr = atten(right,22,i,yhelp) + if (iz.le.1) xkorr=xkorr/hepratio + tof22(right,i,iadc) = xkorr/cos(theta13) + adcflagtof(ch22b(i),hb22b(i)) = 1 endif - ENDIF + ENDIF C-----------------------------S3 -------------------------------- - yhelp=0. - if (tof32_i.GT.none_find) yhelp=tof32_y(tof32_i) - if (ytofpos(3).lt.100) yhelp=ytofpos(3) - - IF (tof31_i.GT.none_find.AND.abs(yhelp).lt.100) THEN +c yhelp=0. + yhelp=100. ! WM + if (tof32_i.GT.none_find) yhelp=tof32_y(tof32_i) + if (ytofpos(3).lt.100) yhelp=ytofpos(3) + IF (tof31_i.GT.none_find.AND.abs(yhelp).lt.100) THEN i = tof31_i - tof31(left,i,iadc) = tof31(left,i,iadc)*cos(theta13) - if (tof31(left,i,iadc).lt.4095) then - xkorr=adcx31(left,i,1)*exp(-yhelp/adcx31(left,i,2)) - xkorr0=adcx31(left,i,1) - adctof_c(ch31a(i),hb31a(i))=tof31(left,i,iadc)/xkorr + if (adc(ch31a(i),hb31a(i)).eq.4095) then + xkorr = atten(left,31,i,yhelp) + if (iz.le.1) xkorr=xkorr/hepratio + tof31(left,i,iadc) = xkorr/cos(theta13) + adcflagtof(ch31a(i),hb31a(i)) = 1 endif - - tof31(right,i,iadc) = tof31(right,i,iadc)*cos(theta13) - if (tof31(right,i,iadc).lt.4095) then - xkorr=adcx31(right,i,1)*exp(yhelp/adcx31(right,i,2)) - xkorr0=adcx31(right,i,1) - adctof_c(ch31b(i),hb31b(i))=tof31(right,i,iadc)/xkorr + if (adc(ch31b(i),hb31b(i)).eq.4095) then + xkorr = atten(right,31,i,yhelp) + if (iz.le.1) xkorr=xkorr/hepratio + tof31(right,i,iadc) = xkorr/cos(theta13) + adcflagtof(ch31b(i),hb31b(i)) = 1 endif - ENDIF - - xhelp=0. - if (tof31_i.GT.none_find) xhelp=tof31_x(tof31_i) - if (xtofpos(3).lt.100) xhelp=xtofpos(3) + ENDIF - IF (tof32_i.GT.none_find.AND.abs(xhelp).lt.100) THEN +c xhelp=0. + xhelp=100. ! WM + if (tof31_i.GT.none_find) xhelp=tof31_x(tof31_i) + if (xtofpos(3).lt.100) xhelp=xtofpos(3) + IF (tof32_i.GT.none_find.AND.abs(xhelp).lt.100) THEN i = tof32_i - tof32(left,i,iadc) = tof32(left,i,iadc)*cos(theta13) - if (tof32(left,i,iadc).lt.4095) then - xkorr=adcx32(left,i,1)*exp(-xhelp/adcx32(left,i,2)) - xkorr0=adcx32(left,i,1) - adctof_c(ch32a(i),hb32a(i))=tof32(left,i,iadc)/xkorr + if (adc(ch32a(i),hb32a(i)).eq.4095) then + xkorr = atten(left,32,i,xhelp) + if (iz.le.1) xkorr=xkorr/hepratio + tof32(left,i,iadc) = xkorr/cos(theta13) + adcflagtof(ch32a(i),hb32a(i)) = 1 endif - - tof32(right,i,iadc) = tof32(right,i,iadc)*cos(theta13) - if (tof32(right,i,iadc).lt.4095) then - xkorr=adcx32(right,i,1)*exp(xhelp/adcx32(right,i,2)) - xkorr0=adcx32(right,i,1) - adctof_c(ch32b(i),hb32b(i))=tof32(right,i,iadc)/xkorr + if (adc(ch32b(i),hb32b(i)).eq.4095) then + xkorr = atten(right,32,i,xhelp) + if (iz.le.1) xkorr=xkorr/hepratio + tof32(right,i,iadc) = xkorr/cos(theta13) + adcflagtof(ch32b(i),hb32b(i)) = 1 endif - ENDIF + ENDIF + + +C------------------------------------------------------------------- +C--------------------Time walk correction ------------------------- +C------------------------------------------------------------------- +C------------------------------------------------------------------- +C Now there is for each hitted paddle a TDC and ADC value, if the +C TDC was < 4095. +C There might be also TDC-ADC pairs in paddles not hitted + +C------------------------------------------------------------------- +C If we have multiple paddles hit, so that no artificial ADC value +C is created, we set the raw TDC value as "tdc_c" +C------------------------------------------------------------------- +c +c do i=1,4 +c do j=1,12 +c tdc_c(i,j) = tdc(i,j) +c enddo +c enddo +c +C---- Let's correct the raw TDC value with the time walk --------- + + DO i=1,8 + if ((tdc(ch11a(i),hb11a(i)).lt.4095).and. + & (tof11(left,i,iadc).lt.3786)) THEN + xhelp = tw11(left,i)/(tof11(left,i,iadc)**0.5) + tof11(left,i,itdc) = tof11(left,i,itdc) + xhelp + tdc_c(ch11a(i),hb11a(i))=tof11(left,i,itdc) + ENDIF + + if ((tdc(ch11b(i),hb11b(i)).lt.4095).and. + & (tof11(right,i,iadc).lt.3786)) THEN + xhelp = tw11(right,i)/(tof11(right,i,iadc)**0.5) + tof11(right,i,itdc) = tof11(right,i,itdc) + xhelp + tdc_c(ch11b(i),hb11b(i))=tof11(right,i,itdc) + ENDIF + ENDDO + + + DO i=1,6 + if ((tdc(ch12a(i),hb12a(i)).lt.4095).and. + & (tof12(left,i,iadc).lt.3786)) THEN + xhelp = tw12(left,i)/(tof12(left,i,iadc)**0.5) + tof12(left,i,itdc) = tof12(left,i,itdc) + xhelp + tdc_c(ch12a(i),hb12a(i))=tof12(left,i,itdc) + ENDIF + + if ((tdc(ch12b(i),hb12b(i)).lt.4095).and. + & (tof12(right,i,iadc).lt.3786)) THEN + xhelp = tw12(right,i)/(tof12(right,i,iadc)**0.5) + tof12(right,i,itdc) = tof12(right,i,itdc) + xhelp + tdc_c(ch12b(i),hb12b(i))=tof12(right,i,itdc) + ENDIF + ENDDO + +C---- + DO I=1,2 + if ((tdc(ch21a(i),hb21a(i)).lt.4095).and. + & (tof21(left,i,iadc).lt.3786)) THEN + xhelp = tw21(left,i)/(tof21(left,i,iadc)**0.5) + tof21(left,i,itdc) = tof21(left,i,itdc) + xhelp + tdc_c(ch21a(i),hb21a(i))=tof21(left,i,itdc) + ENDIF + + if ((tdc(ch21b(i),hb21b(i)).lt.4095).and. + & (tof21(right,i,iadc).lt.3786)) THEN + xhelp = tw21(right,i)/(tof21(right,i,iadc)**0.5) + tof21(right,i,itdc) = tof21(right,i,itdc) + xhelp + tdc_c(ch21b(i),hb21b(i))=tof21(right,i,itdc) + ENDIF + ENDDO + + DO I=1,2 + if ((tdc(ch22a(i),hb22a(i)).lt.4095).and. + & (tof22(left,i,iadc).lt.3786)) THEN + xhelp = tw22(left,i)/(tof22(left,i,iadc)**0.5) + tof22(left,i,itdc) = tof22(left,i,itdc) + xhelp + tdc_c(ch22a(i),hb22a(i))=tof22(left,i,itdc) + ENDIF + + if ((tdc(ch22b(i),hb22b(i)).lt.4095).and. + & (tof22(right,i,iadc).lt.3786)) THEN + xhelp = tw22(right,i)/(tof22(right,i,iadc)**0.5) + tof22(right,i,itdc) = tof22(right,i,itdc) + xhelp + tdc_c(ch22b(i),hb22b(i))=tof22(right,i,itdc) + ENDIF + ENDDO + +C---- + DO I=1,3 + if ((tdc(ch31a(i),hb31a(i)).lt.4095).and. + & (tof31(left,i,iadc).lt.3786)) THEN + xhelp = tw31(left,i)/(tof31(left,i,iadc)**0.5) + tof31(left,i,itdc) = tof31(left,i,itdc) + xhelp + tdc_c(ch31a(i),hb31a(i))=tof31(left,i,itdc) + ENDIF + + if ((tdc(ch31b(i),hb31b(i)).lt.4095).and. + & (tof31(right,i,iadc).lt.3786)) THEN + xhelp = tw31(right,i)/(tof31(right,i,iadc)**0.5) + tof31(right,i,itdc) = tof31(right,i,itdc) + xhelp + tdc_c(ch31b(i),hb31b(i))=tof31(right,i,itdc) + ENDIF + ENDDO + + DO I=1,3 + if ((tdc(ch32a(i),hb32a(i)).lt.4095).and. + & (tof32(left,i,iadc).lt.3786)) THEN + xhelp = tw32(left,i)/(tof32(left,i,iadc)**0.5) + tof32(left,i,itdc) = tof32(left,i,itdc) + xhelp + tdc_c(ch32a(i),hb32a(i))=tof32(left,i,itdc) + ENDIF + + if ((tdc(ch32b(i),hb32b(i)).lt.4095).and. + & (tof32(right,i,iadc).lt.3786)) THEN + xhelp = tw32(right,i)/(tof32(right,i,iadc)**0.5) + tof32(right,i,itdc) = tof32(right,i,itdc) + xhelp + tdc_c(ch32b(i),hb32b(i))=tof32(right,i,itdc) + ENDIF + ENDDO + + +C--------------------------------------------------------------- +C--- calculate track position in paddle using timing difference +C--- now using the time-walk corrected TDC values +C--------------------------------------------------------------- + + do i=1,3 + xtofpos(i)=100. + ytofpos(i)=100. + enddo + +C-----------------------------S1 -------------------------------- -C----------------------------------------------------------------------- + IF (tof11_i.GT.none_find) THEN + ytofpos(1) = ((tof11(1,tof11_i,itdc)-tof11(2,tof11_i,itdc))/2. + + -y_coor_lin11(tof11_i,offset))/y_coor_lin11(tof11_i,slope) + i=tof11_i + endif + + IF (tof12_i.GT.none_find) THEN + xtofpos(1) = ((tof12(1,tof12_i,itdc)-tof12(2,tof12_i,itdc))/2. + + -x_coor_lin12(tof12_i,offset))/x_coor_lin12(tof12_i,slope) + i=tof12_i + endif + + +C-----------------------------S2 -------------------------------- + + IF (tof21_i.GT.none_find) THEN + xtofpos(2) = ((tof21(1,tof21_i,itdc)-tof21(2,tof21_i,itdc))/2. + + -x_coor_lin21(tof21_i,offset))/x_coor_lin21(tof21_i,slope) + i=tof21_i + endif + + IF (tof22_i.GT.none_find) THEN + ytofpos(2) = ((tof22(1,tof22_i,itdc)-tof22(2,tof22_i,itdc))/2. + + -y_coor_lin22(tof22_i,offset))/y_coor_lin22(tof22_i,slope) + i=tof22_i + endif + + +C-----------------------------S3 -------------------------------- + + IF (tof31_i.GT.none_find) THEN + ytofpos(3) = ((tof31(1,tof31_i,itdc)-tof31(2,tof31_i,itdc))/2. + + -y_coor_lin31(tof31_i,offset))/y_coor_lin31(tof31_i,slope) + i=tof31_i + endif + + IF (tof32_i.GT.none_find) THEN + xtofpos(3) = ((tof32(1,tof32_i,itdc)-tof32(2,tof32_i,itdc))/2. + + -x_coor_lin32(tof32_i,offset))/x_coor_lin32(tof32_i,slope) + i=tof32_i + endif + + +c do i=1,3 +c if (abs(xtofpos(i)).gt.100.) then +c xtofpos(i)=101. +c endif +c if (abs(ytofpos(i)).gt.100.) then +c ytofpos(i)=101. +c endif +c enddo + + +C-- restrict TDC measurements to physical paddle dimensions +/- 10 cm +C-- this cut is now stronger than in the old versions + + if (abs(xtofpos(1)).gt.31.) xtofpos(1)=101. + if (abs(xtofpos(2)).gt.19.) xtofpos(2)=101. + if (abs(xtofpos(3)).gt.19.) xtofpos(3)=101. + + if (abs(ytofpos(1)).gt.26.) ytofpos(1)=101. + if (abs(ytofpos(2)).gt.18.) ytofpos(2)=101. + if (abs(ytofpos(3)).gt.18.) ytofpos(3)=101. + +C---------------------------------------------------------------------- +C--------------------- zenith angle theta --------------------------- +C---------------------------------------------------------------------- +C------------------- improved calculation --------------------------- + + xhelp1=0. + if (tof11_i.GT.none_find) xhelp1=tof11_x(tof11_i) + if (xtofpos(1).lt.100) xhelp1=xtofpos(1) + + yhelp1=0. + if (tof12_i.GT.none_find) yhelp1=tof12_y(tof12_i) + if (ytofpos(1).lt.100) yhelp1=ytofpos(1) + + yhelp2=0. + if (tof32_i.GT.none_find) yhelp2=tof32_y(tof32_i) + if (ytofpos(3).lt.100) yhelp2=ytofpos(3) + + xhelp2=0. + if (tof31_i.GT.none_find) xhelp2=tof31_x(tof31_i) + if (xtofpos(3).lt.100) xhelp2=xtofpos(3) + + + dx=0. + dy=0. + dr=0. + theta13 = 0. + + dx = xhelp1 - xhelp2 + dy = yhelp1 - yhelp2 + dr = sqrt(dx*dx+dy*dy) + theta13 = atan(dr/tofarm13) + + +C------------------------------------------------------------------ +C------------------------------------------------------------------ +C-------angle and ADC(x) correction: moved to new dEdx routine +C------------------------------------------------------------------ +C------------------------------------------------------------------ + +C-------------------------------------------------------------------- C----------------------calculate Beta ------------------------------ -C----------------------------------------------------------------------- -C-------------------difference of sums --------------------------- +C-------------------------------------------------------------------- +C-------------------difference of sums ----------------------------- C C DS = (t1+t2) - t3+t4) C DS = c1 + c2/beta*cos(theta) @@ -759,13 +1071,19 @@ C since TDC resolution varies slightly c2 has to be calibrated C S11 - S31 - IF (tof11_i.GT.none_find.AND.tof31_i.GT.none_find) THEN + + dist = ZTOF(1) - ZTOF(5) + c2 = (2.*0.01*dist)/(3.E08*50.E-12 ) + + IF ((tof11_i.GT.none_find).AND.(tof31_i.GT.none_find).AND. + & (ytofpos(1).NE.101.).AND.(ytofpos(3).NE.101.)) THEN xhelp1 = tof11(1,tof11_i,itdc)+tof11(2,tof11_i,itdc) xhelp2 = tof31(1,tof31_i,itdc)+tof31(2,tof31_i,itdc) ds = xhelp1-xhelp2 ihelp=(tof11_i-1)*3+tof31_i - c1 = k_S11S31(1,ihelp) - c2 = k_S11S31(2,ihelp) + if (iz.le.1) c1 = k_S11S31(1,ihelp) + if (iz.eq.2) c1 = k_S11S31(2,ihelp) + if (iz.gt.2) c1 = k_S11S31(3,ihelp) betatof_a(1) = c2/(cos(theta13)*(ds-c1)) C------- ToF Mask - S11 - S31 @@ -785,13 +1103,19 @@ ENDIF C S11 - S32 - IF (tof11_i.GT.none_find.AND.tof32_i.GT.none_find) THEN + + dist = ZTOF(1) - ZTOF(6) + c2 = (2.*0.01*dist)/(3.E08*50.E-12 ) + + IF ((tof11_i.GT.none_find).AND.(tof32_i.GT.none_find).AND. + & (ytofpos(1).NE.101.).AND.(xtofpos(3).NE.101.)) THEN xhelp1 = tof11(1,tof11_i,itdc)+tof11(2,tof11_i,itdc) xhelp2 = tof32(1,tof32_i,itdc)+tof32(2,tof32_i,itdc) ds = xhelp1-xhelp2 ihelp=(tof11_i-1)*3+tof32_i - c1 = k_S11S32(1,ihelp) - c2 = k_S11S32(2,ihelp) + if (iz.le.1) c1 = k_S11S32(1,ihelp) + if (iz.eq.2) c1 = k_S11S32(2,ihelp) + if (iz.gt.2) c1 = k_S11S32(3,ihelp) betatof_a(2) = c2/(cos(theta13)*(ds-c1)) C------- ToF Mask - S11 - S32 @@ -811,13 +1135,19 @@ ENDIF C S12 - S31 - IF (tof12_i.GT.none_find.AND.tof31_i.GT.none_find) THEN + + dist = ZTOF(2) - ZTOF(5) + c2 = (2.*0.01*dist)/(3.E08*50.E-12 ) + + IF ((tof12_i.GT.none_find).AND.(tof31_i.GT.none_find).AND. + & (xtofpos(1).NE.101.).AND.(ytofpos(3).NE.101.)) THEN xhelp1 = tof12(1,tof12_i,itdc)+tof12(2,tof12_i,itdc) xhelp2 = tof31(1,tof31_i,itdc)+tof31(2,tof31_i,itdc) ds = xhelp1-xhelp2 ihelp=(tof12_i-1)*3+tof31_i - c1 = k_S12S31(1,ihelp) - c2 = k_S12S31(2,ihelp) + if (iz.le.1) c1 = k_S12S31(1,ihelp) + if (iz.eq.2) c1 = k_S12S31(2,ihelp) + if (iz.gt.2) c1 = k_S12S31(3,ihelp) betatof_a(3) = c2/(cos(theta13)*(ds-c1)) C------- ToF Mask - S12 - S31 @@ -837,13 +1167,19 @@ ENDIF C S12 - S32 - IF (tof12_i.GT.none_find.AND.tof32_i.GT.none_find) THEN + + dist = ZTOF(2) - ZTOF(6) + c2 = (2.*0.01*dist)/(3.E08*50.E-12 ) + + IF ((tof12_i.GT.none_find).AND.(tof32_i.GT.none_find).AND. + & (xtofpos(1).NE.101.).AND.(xtofpos(3).NE.101.)) THEN xhelp1 = tof12(1,tof12_i,itdc)+tof12(2,tof12_i,itdc) xhelp2 = tof32(1,tof32_i,itdc)+tof32(2,tof32_i,itdc) ds = xhelp1-xhelp2 ihelp=(tof12_i-1)*3+tof32_i - c1 = k_S12S32(1,ihelp) - c2 = k_S12S32(2,ihelp) + if (iz.le.1) c1 = k_S12S32(1,ihelp) + if (iz.eq.2) c1 = k_S12S32(2,ihelp) + if (iz.gt.2) c1 = k_S12S32(3,ihelp) betatof_a(4) = c2/(cos(theta13)*(ds-c1)) C------- ToF Mask - S12 - S32 @@ -863,14 +1199,20 @@ ENDIF C S21 - S31 - IF (tof21_i.GT.none_find.AND.tof31_i.GT.none_find) THEN + + dist = ZTOF(3) - ZTOF(5) + c2 = (2.*0.01*dist)/(3.E08*50.E-12 ) + + IF ((tof21_i.GT.none_find).AND.(tof31_i.GT.none_find).AND. + & (xtofpos(2).NE.101.).AND.(ytofpos(3).NE.101.)) THEN xhelp1 = tof21(1,tof21_i,itdc)+tof21(2,tof21_i,itdc) xhelp2 = tof31(1,tof31_i,itdc)+tof31(2,tof31_i,itdc) ds = xhelp1-xhelp2 ihelp=(tof21_i-1)*3+tof31_i - c1 = k_S21S31(1,ihelp) - c2 = k_S21S31(2,ihelp) - betatof_a(5) = c2/(cos(theta23)*(ds-c1)) + if (iz.le.1) c1 = k_S21S31(1,ihelp) + if (iz.eq.2) c1 = k_S21S31(2,ihelp) + if (iz.gt.2) c1 = k_S21S31(3,ihelp) + betatof_a(5) = c2/(cos(theta13)*(ds-c1)) C------- ToF Mask - S21 - S31 @@ -889,14 +1231,21 @@ ENDIF C S21 - S32 - IF (tof21_i.GT.none_find.AND.tof32_i.GT.none_find) THEN + + dist = ZTOF(3) - ZTOF(6) + c2 = (2.*0.01*dist)/(3.E08*50.E-12 ) + + + IF ((tof21_i.GT.none_find).AND.(tof32_i.GT.none_find).AND. + & (xtofpos(2).NE.101.).AND.(xtofpos(3).NE.101.)) THEN xhelp1 = tof21(1,tof21_i,itdc)+tof21(2,tof21_i,itdc) xhelp2 = tof32(1,tof32_i,itdc)+tof32(2,tof32_i,itdc) ds = xhelp1-xhelp2 ihelp=(tof21_i-1)*3+tof32_i - c1 = k_S21S32(1,ihelp) - c2 = k_S21S32(2,ihelp) - betatof_a(6) = c2/(cos(theta23)*(ds-c1)) + if (iz.le.1) c1 = k_S21S32(1,ihelp) + if (iz.eq.2) c1 = k_S21S32(2,ihelp) + if (iz.gt.2) c1 = k_S21S32(3,ihelp) + betatof_a(6) = c2/(cos(theta13)*(ds-c1)) C------- ToF Mask - S21 - S32 @@ -915,13 +1264,19 @@ ENDIF C S22 - S31 - IF (tof22_i.GT.none_find.AND.tof31_i.GT.none_find) THEN + + dist = ZTOF(4) - ZTOF(5) + c2 = (2.*0.01*dist)/(3.E08*50.E-12 ) + + IF ((tof22_i.GT.none_find).AND.(tof31_i.GT.none_find).AND. + & (ytofpos(2).NE.101.).AND.(ytofpos(3).NE.101.)) THEN xhelp1 = tof22(1,tof22_i,itdc)+tof22(2,tof22_i,itdc) xhelp2 = tof31(1,tof31_i,itdc)+tof31(2,tof31_i,itdc) ds = xhelp1-xhelp2 ihelp=(tof22_i-1)*3+tof31_i - c1 = k_S22S31(1,ihelp) - c2 = k_S22S31(2,ihelp) + if (iz.le.1) c1 = k_S22S31(1,ihelp) + if (iz.eq.2) c1 = k_S22S31(2,ihelp) + if (iz.gt.2) c1 = k_S22S31(3,ihelp) betatof_a(7) = c2/(cos(theta13)*(ds-c1)) C------- ToF Mask - S22 - S31 @@ -941,13 +1296,19 @@ ENDIF C S22 - S32 - IF (tof22_i.GT.none_find.AND.tof32_i.GT.none_find) THEN + + dist = ZTOF(4) - ZTOF(6) + c2 = (2.*0.01*dist)/(3.E08*50.E-12 ) + + IF ((tof22_i.GT.none_find).AND.(tof32_i.GT.none_find).AND. + & (ytofpos(2).NE.101.).AND.(xtofpos(3).NE.101.)) THEN xhelp1 = tof22(1,tof22_i,itdc)+tof22(2,tof22_i,itdc) xhelp2 = tof32(1,tof32_i,itdc)+tof32(2,tof32_i,itdc) ds = xhelp1-xhelp2 ihelp=(tof22_i-1)*3+tof32_i - c1 = k_S22S32(1,ihelp) - c2 = k_S22S32(2,ihelp) + if (iz.le.1) c1 = k_S22S32(1,ihelp) + if (iz.eq.2) c1 = k_S22S32(2,ihelp) + if (iz.gt.2) c1 = k_S22S32(3,ihelp) betatof_a(8) = c2/(cos(theta13)*(ds-c1)) C------- ToF Mask - S22 - S32 @@ -967,13 +1328,19 @@ ENDIF C S11 - S21 - IF (tof11_i.GT.none_find.AND.tof21_i.GT.none_find) THEN + + dist = ZTOF(1) - ZTOF(3) + c2 = (2.*0.01*dist)/(3.E08*50.E-12 ) + + IF ((tof11_i.GT.none_find).AND.(tof21_i.GT.none_find).AND. + & (ytofpos(1).NE.101.).AND.(xtofpos(2).NE.101.)) THEN xhelp1 = tof11(1,tof11_i,itdc)+tof11(2,tof11_i,itdc) xhelp2 = tof21(1,tof21_i,itdc)+tof21(2,tof21_i,itdc) ds = xhelp1-xhelp2 ihelp=(tof11_i-1)*2+tof21_i - c1 = k_S11S21(1,ihelp) - c2 = k_S11S21(2,ihelp) + if (iz.le.1) c1 = k_S11S21(1,ihelp) + if (iz.eq.2) c1 = k_S11S21(2,ihelp) + if (iz.gt.2) c1 = k_S11S21(3,ihelp) betatof_a(9) = c2/(cos(theta13)*(ds-c1)) C------- ToF Mask - S11 - S21 @@ -993,13 +1360,19 @@ ENDIF C S11 - S22 - IF (tof11_i.GT.none_find.AND.tof22_i.GT.none_find) THEN + + dist = ZTOF(1) - ZTOF(4) + c2 = (2.*0.01*dist)/(3.E08*50.E-12 ) + + IF ((tof11_i.GT.none_find).AND.(tof22_i.GT.none_find).AND. + & (ytofpos(1).NE.101.).AND.(ytofpos(2).NE.101.)) THEN xhelp1 = tof11(1,tof11_i,itdc)+tof11(2,tof11_i,itdc) xhelp2 = tof22(1,tof22_i,itdc)+tof22(2,tof22_i,itdc) ds = xhelp1-xhelp2 ihelp=(tof11_i-1)*2+tof22_i - c1 = k_S11S22(1,ihelp) - c2 = k_S11S22(2,ihelp) + if (iz.le.1) c1 = k_S11S22(1,ihelp) + if (iz.eq.2) c1 = k_S11S22(2,ihelp) + if (iz.gt.2) c1 = k_S11S22(3,ihelp) betatof_a(10) = c2/(cos(theta13)*(ds-c1)) C------- ToF Mask - S11 - S22 @@ -1019,13 +1392,19 @@ ENDIF C S12 - S21 - IF (tof12_i.GT.none_find.AND.tof21_i.GT.none_find) THEN + + dist = ZTOF(2) - ZTOF(3) + c2 = (2.*0.01*dist)/(3.E08*50.E-12 ) + + IF ((tof12_i.GT.none_find).AND.(tof21_i.GT.none_find).AND. + & (xtofpos(1).NE.101.).AND.(xtofpos(2).NE.101.)) THEN xhelp1 = tof12(1,tof12_i,itdc)+tof12(2,tof12_i,itdc) xhelp2 = tof21(1,tof21_i,itdc)+tof21(2,tof21_i,itdc) ds = xhelp1-xhelp2 ihelp=(tof12_i-1)*2+tof21_i - c1 = k_S12S21(1,ihelp) - c2 = k_S12S21(2,ihelp) + if (iz.le.1) c1 = k_S12S21(1,ihelp) + if (iz.eq.2) c1 = k_S12S21(2,ihelp) + if (iz.gt.2) c1 = k_S12S21(3,ihelp) betatof_a(11) = c2/(cos(theta13)*(ds-c1)) C------- ToF Mask - S12 - S21 @@ -1045,13 +1424,19 @@ ENDIF C S12 - S22 - IF (tof12_i.GT.none_find.AND.tof22_i.GT.none_find) THEN + + dist = ZTOF(2) - ZTOF(4) + c2 = (2.*0.01*dist)/(3.E08*50.E-12 ) + + IF ((tof12_i.GT.none_find).AND.(tof22_i.GT.none_find).AND. + & (xtofpos(1).NE.101.).AND.(ytofpos(2).NE.101.)) THEN xhelp1 = tof12(1,tof12_i,itdc)+tof12(2,tof12_i,itdc) xhelp2 = tof22(1,tof22_i,itdc)+tof22(2,tof22_i,itdc) ds = xhelp1-xhelp2 ihelp=(tof12_i-1)*2+tof22_i - c1 = k_S12S22(1,ihelp) - c2 = k_S12S22(2,ihelp) + if (iz.le.1) c1 = k_S12S22(1,ihelp) + if (iz.eq.2) c1 = k_S12S22(2,ihelp) + if (iz.gt.2) c1 = k_S12S22(3,ihelp) betatof_a(12) = c2/(cos(theta13)*(ds-c1)) C------- ToF Mask - S12 - S22 @@ -1071,28 +1456,786 @@ ENDIF C--------------------------------------------------------- +C +C icount=0 +C sw=0. +C sxw=0. +C beta_mean=100. +C +C do i=1,12 +C if ((betatof_a(i).gt.-1.5).and.(betatof_a(i).lt.1.5)) then +C icount= icount+1 +C if (i.le.4) w_i=1./(0.13**2.) +C if ((i.ge.5).and.(i.le.8)) w_i=1./(0.16**2.) +C if (i.ge.9) w_i=1./(0.25**2.) ! to be checked +C sxw=sxw + betatof_a(i)*w_i +C sw =sw + w_i +C endif +C enddo +C +C if (icount.gt.0) beta_mean=sxw/sw +C betatof_a(13) = beta_mean +C - icount=0 - sw=0. - sxw=0. - beta_mean=100. +C-------- New mean beta calculation ----------------------- - do i=1,12 - if ((betatof_a(i).gt.-1.5).and.(betatof_a(i).lt.1.5)) then - icount= icount+1 - if (i.le.4) w_i=1./(0.13**2.) - if ((i.ge.5).and.(i.le.8)) w_i=1./(0.16**2.) - if (i.ge.9) w_i=1./(0.25**2.) ! to be checked - sxw=sxw + betatof_a(i)*w_i - sw =sw + w_i - endif - enddo + do i=1,12 + btemp(i) = betatof_a(i) + enddo + + betatof_a(13)=newbeta(1,btemp,hitvec,10.,10.,20.) + +C-------------------------------------------------------------- +C write(*,*) betatof_a +c write(*,*) xtofpos +c write(*,*) ytofpos +c write(*,*)'tofl2com beta', betatof_a +C write(*,*) adcflagtof +c write(*,*) 'tofl2com' +c write(*,*) xtofpos +c write(*,*) ytofpos +c write(*,*) xtr_tof +c write(*,*) ytr_tof - if (icount.gt.0) beta_mean=sxw/sw - betatof_a(13) = beta_mean - - 100 continue +c 100 continue + continue C RETURN END + + +C------------------------------------------------------------------ +C------------------------------------------------------------------ + + function atten(is,ilay,ipad,x) + include 'input_tof.txt' + real atten + real x + real xmin,xmax + integer ilay,ipad + +* S11 8 paddles 33.0 x 5.1 cm +* S12 6 paddles 40.8 x 5.5 cm +* S21 2 paddles 18.0 x 7.5 cm +* S22 2 paddles 15.0 x 9.0 cm +* S31 3 paddles 15.0 x 6.0 cm +* S32 3 paddles 18.0 x 5.0 cm + + +c if (ilay.eq.11) write(*,*) 'start ',ipad,is,adcx11(is,ipad,1), +c & adcx11(is,ipad,2),adcx11(is,ipad,3),adcx11(is,ipad,4) +c if (ilay.eq.12) write(*,*) 'start ',ipad,is,adcx12(is,ipad,1), +c & adcx12(is,ipad,2),adcx12(is,ipad,3),adcx12(is,ipad,4) + + + if (ilay.eq.11) xmin=-33.0/2. + if (ilay.eq.11) xmax= 33.0/2. + if (ilay.eq.12) xmin=-40.8/2. + if (ilay.eq.12) xmax= 40.8/2. + + if (ilay.eq.21) xmin=-18.0/2. + if (ilay.eq.21) xmax= 18.0/2. + if (ilay.eq.22) xmin=-15.0/2. + if (ilay.eq.22) xmax= 15.0/2. + + if (ilay.eq.31) xmin=-15.0/2. + if (ilay.eq.31) xmax= 15.0/2. + if (ilay.eq.32) xmin=-18.0/2. + if (ilay.eq.32) xmax= 18.0/2. + + if (x .lt. xmin) x=xmin + if (x .gt. xmax) x=xmax + + + if (ilay.eq.11) atten= + & adcx11(is,ipad,1)*exp(x*adcx11(is,ipad,2)) + & + adcx11(is,ipad,3)*exp(x*adcx11(is,ipad,4)) + + if (ilay.eq.12) atten= + & adcx12(is,ipad,1)*exp(x*adcx12(is,ipad,2)) + & + adcx12(is,ipad,3)*exp(x*adcx12(is,ipad,4)) + + if (ilay.eq.21) atten= + & adcx21(is,ipad,1)*exp(x*adcx21(is,ipad,2)) + & + adcx21(is,ipad,3)*exp(x*adcx21(is,ipad,4)) + + if (ilay.eq.22) atten= + & adcx22(is,ipad,1)*exp(x*adcx22(is,ipad,2)) + & + adcx22(is,ipad,3)*exp(x*adcx22(is,ipad,4)) + + if (ilay.eq.31) atten= + & adcx31(is,ipad,1)*exp(x*adcx31(is,ipad,2)) + & + adcx31(is,ipad,3)*exp(x*adcx31(is,ipad,4)) + + if (ilay.eq.32) atten= + & adcx32(is,ipad,1)*exp(x*adcx32(is,ipad,2)) + & + adcx32(is,ipad,3)*exp(x*adcx32(is,ipad,4)) + + if (atten.gt.10000) atten=10000. + + end + +C------------------------------------------------------------------ +C------------------------------------------------------------------ + + function pc_adc(ix) + include 'input_tof.txt' + real pc_adc + integer ix + + pc_adc=28.0407 + 0.628929*ix + & - 5.80901e-05*ix*ix + 3.14092e-08*ix*ix*ix +c write(*,*) ix,pc_adc + end + +C------------------------------------------------------------------ +C------------------------------------------------------------------ + + function check_charge(theta,hitvec) + + include 'input_tof.txt' + include 'tofcomm.txt' + + real check_charge + integer hitvec(6) + REAL CHARGE, theta + +C upper and lower limits for the helium selection + REAL A_l(24),A_h(24) + DATA A_l /200,190,300,210,220,200,210,60,60,120,220, + & 120,160,50,300,200,120,250,350,300,350,250,280,300/ + DATA A_h /550,490,800,600,650,600,600,260,200,380, + & 620,380,550,200,850,560,400,750,900,800,880,800,750,800/ + +C The k1 constants for the beta calculation, only for S1-S3 +C k2 constant is taken to be the standard 2D/c + REAL k1(84) + DATA k1 /50,59.3296,28.4328,-26.0818,5.91253,-19.588, + & -9.26316,24.7544,2.32465,-50.5058,-15.3195,-39.1443, + & -91.2546,-58.6243,-84.5641,-63.1516,-32.2091,-58.3358, + & 13.8084,45.5322,33.2416,-11.5313,51.3271,75,-14.1141, + & 42.8466,15.1794,-63.6672,-6.07739,-32.164,-41.771,10.5274, + & -9.46096,-81.7404,-28.783,-52.7167,-127.394,-69.6166, + & -93.4655,-98.9543,-42.863,-67.8244,-19.3238,31.1221,8.7319, + & -43.1627,5.55573,-14.4078,-83.4466,-47.4647,-77.8379, + & -108.222,-75.986,-101.297,-96.0205,-63.1881,-90.1372, + & -22.7347,8.31409,-19.6912,-7.49008,23.6979,-1.66677, + & 1.81556,34.4668,6.23693,-100,-59.5861,-90.9159,-141.639, + & -89.2521,-112.881,-130.199,-77.0357,-98.4632,-60.2086, + & -4.82097,-29.3705,-43.6469,10.5884,-9.31304,-35.3329, + & 25.2514,25.6/ + + + + REAL zin(6) + DATA zin /53.74, 53.04, 23.94, 23.44, -23.49, -24.34/ + + REAL c1,c2,xhelp,xhelp1,xhelp2,ds,dist,F + REAL sw,sxw,beta_mean_tof,w_i + INTEGER ihelp + INTEGER ipmt(4) + REAL time(4),beta1(4) + + REAL adca(48), tdca(48) + + REAL a1,a2 + INTEGER jj + +c get rid of warnings EMILIANO + i = 0 + slope = 0 + offset = 0 + none_find = 0 + none_ev = 0 + adc_ev = 0 + tdc_ev = 0 + iadc = 0 + itdc = 0 + right = 0 + left = 0 + tof12_y(1) = tof12_y(1) + tof11_x(1) = tof11_x(1) + tof21_y(1) = tof21_y(1) + tof22_x(1) = tof22_x(1) + tof32_y(1) = tof32_y(1) + tof31_x(1) = tof31_x(1) +c get rid of warnings + +C----------------------------------------------------------- +C--- get data +C----------------------------------------------------------- + + do j=1,8 + ih = 1 + 0 +((j-1)*2) + adca(ih) = adc(ch11a(j),hb11a(j)) + adca(ih+1) = adc(ch11b(j),hb11b(j)) + tdca(ih) = tdc(ch11a(j),hb11a(j)) + tdca(ih+1) = tdc(ch11b(j),hb11b(j)) + enddo + + do j=1,6 + ih = 1 + 16+((j-1)*2) + adca(ih) = adc(ch12a(j),hb12a(j)) + adca(ih+1) = adc(ch12b(j),hb12b(j)) + tdca(ih) = tdc(ch12a(j),hb12a(j)) + tdca(ih+1) = tdc(ch12b(j),hb12b(j)) + enddo + + do j=1,2 + ih = 1 + 28+((j-1)*2) + adca(ih) = adc(ch21a(j),hb21a(j)) + adca(ih+1) = adc(ch21b(j),hb21b(j)) + tdca(ih) = tdc(ch21a(j),hb21a(j)) + tdca(ih+1) = tdc(ch21b(j),hb21b(j)) + enddo + + do j=1,2 + ih = 1 + 32+((j-1)*2) + adca(ih) = adc(ch22a(j),hb22a(j)) + adca(ih+1) = adc(ch22b(j),hb22b(j)) + tdca(ih) = tdc(ch22a(j),hb22a(j)) + tdca(ih+1) = tdc(ch22b(j),hb22b(j)) + enddo + + do j=1,3 + ih = 1 + 36+((j-1)*2) + adca(ih) = adc(ch31a(j),hb31a(j)) + adca(ih+1) = adc(ch31b(j),hb31b(j)) + tdca(ih) = tdc(ch31a(j),hb31a(j)) + tdca(ih+1) = tdc(ch31b(j),hb31b(j)) + enddo + + do j=1,3 + ih = 1 + 42+((j-1)*2) + adca(ih) = adc(ch32a(j),hb32a(j)) + adca(ih+1) = adc(ch32b(j),hb32b(j)) + tdca(ih) = tdc(ch32a(j),hb32a(j)) + tdca(ih+1) = tdc(ch32b(j),hb32b(j)) + enddo + + +c write(*,*) adca +c write(*,*) tdca + + +C============ calculate beta and select charge > Z=1 =============== + + ICHARGE=1 + +C find hitted paddle by looking for ADC values on both sides +C since we looking for Z>1 this gives decent results + + tof11_i = hitvec(1)-1 + tof12_i = hitvec(2)-1 + tof21_i = hitvec(3)-1 + tof22_i = hitvec(4)-1 + tof31_i = hitvec(5)-1 + tof32_i = hitvec(6)-1 + +c write(*,*) ' in charge check' +c write(*,*) theta,tof11_i,tof12_i,tof21_i,tof22_i,tof31_i,tof32_i + +C---------------------------------------------------------------- + + beta_help=100. + beta_mean_tof=100. + + do jj=1,4 + beta1(jj) = 100. + enddo + +C---------------------------------------------------------------- +C--------- S1 - S3 --------------------------------------------- +C---------------------------------------------------------------- + +C--------- S11 - S31 ------------------------------------------- + + if ((tof11_i.gt.-1).and.(tof31_i.gt.-1)) then + + dist = zin(1) - zin(5) + c2 = (2.*0.01*dist)/(3.E08*50.E-12) + F = 1./cos(theta) + + ipmt(1) = (tof11_i)*2+1 + ipmt(2) = (tof11_i)*2+2 + ipmt(3) = 36+(tof31_i)*2+1 + ipmt(4) = 36+(tof31_i)*2+2 + +c write(*,*) ipmt + + do jj=1,4 + time(jj) = tdca(ipmt(jj)) + enddo + +c write(*,*) time + + if ((time(1).lt.4095).and.(time(2).lt.4095).and. + & (time(3).lt.4095).and.(time(4).lt.4095)) then + xhelp1 = time(1) + time(2) + xhelp2 = time(3) + time(4) + ds = xhelp1-xhelp2 + ihelp=0+(tof11_i)*3+tof31_i + c1 = k1(ihelp+1) + beta1(1) = c2*F/(ds-c1); + endif +c write(*,*) beta1(1) + endif ! tof_.... + + +C--------- S11 - S32 ------------------------------------------- + + if ((tof11_i.gt.-1).and.(tof32_i.gt.-1)) then + + dist = zin(1) - zin(6) + c2 = (2.*0.01*dist)/(3.E08*50.E-12) + F = 1./cos(theta) + + ipmt(1) = (tof11_i)*2+1 + ipmt(2) = (tof11_i)*2+2 + ipmt(3) = 42+(tof32_i)*2+1 + ipmt(4) = 42+(tof32_i)*2+2 + + do jj=1,4 + time(jj) = tdca(ipmt(jj)) + enddo + + if ((time(1).lt.4095).and.(time(2).lt.4095).and. + & (time(3).lt.4095).and.(time(4).lt.4095)) then + xhelp1 = time(1) + time(2) + xhelp2 = time(3) + time(4) + ds = xhelp1-xhelp2 + ihelp=24+(tof11_i)*3+tof32_i + c1 = k1(ihelp+1) + beta1(2) = c2*F/(ds-c1); + endif + endif ! tof_.... + + +C--------- S12 - S31 ------------------------------------------- + + if ((tof12_i.gt.-1).and.(tof31_i.gt.-1)) then + + dist = zin(2) - zin(5) + c2 = (2.*0.01*dist)/(3.E08*50.E-12) + F = 1./cos(theta) + + ipmt(1) = 16+(tof12_i)*2+1 + ipmt(2) = 16+(tof12_i)*2+2 + ipmt(3) = 36+(tof31_i)*2+1 + ipmt(4) = 36+(tof31_i)*2+2 + + do jj=1,4 + time(jj) = tdca(ipmt(jj)) + enddo + + if ((time(1).lt.4095).and.(time(2).lt.4095).and. + & (time(3).lt.4095).and.(time(4).lt.4095)) then + xhelp1 = time(1) + time(2) + xhelp2 = time(3) + time(4) + ds = xhelp1-xhelp2 + ihelp=48+(tof12_i)*3+tof31_i + c1 = k1(ihelp+1) + beta1(3) = c2*F/(ds-c1); + endif + endif ! tof_.... + + +C--------- S12 - S32 ------------------------------------------- + + if ((tof12_i.gt.-1).and.(tof32_i.gt.-1)) then + + dist = zin(2) - zin(6) + c2 = (2.*0.01*dist)/(3.E08*50.E-12) + F = 1./cos(theta) + + ipmt(1) = 16+(tof12_i)*2+1 + ipmt(2) = 16+(tof12_i)*2+2 + ipmt(3) = 42+(tof32_i)*2+1 + ipmt(4) = 42+(tof32_i)*2+2 + + do jj=1,4 + time(jj) = tdca(ipmt(jj)) + enddo + + if ((time(1).lt.4095).and.(time(2).lt.4095).and. + & (time(3).lt.4095).and.(time(4).lt.4095)) then + xhelp1 = time(1) + time(2) + xhelp2 = time(3) + time(4) + ds = xhelp1-xhelp2 + ihelp=56+(tof12_i)*3+tof32_i + c1 = k1(ihelp+1) + beta1(4) = c2*F/(ds-c1); + endif + + endif ! tof_.... + +c write(*,*) beta1 + +C---- calculate beta mean, only downward going particles are interesting ---- + + sw=0. + sxw=0. + beta_mean_tof=100. + + do jj=1,4 + if ((beta1(jj).gt.0.1).and.(beta1(jj).lt.2.0)) then + w_i=1./(0.13*0.13) + sxw=sxw + beta1(jj)*w_i + sw =sw + w_i ; + endif + enddo + + if (sw.gt.0) beta_mean_tof=sxw/sw; + +C write(*,*) 'beta_mean_tof ',beta_mean_tof + + beta_help = beta_mean_tof ! pow(beta_mean_tof,1.0) gave best results + +CCCCC endif ! if tof11_i > -1 && ...... beta calculation + +C----------------------- Select charge -------------------------- + + charge=0 + + if ((beta_mean_tof.gt.0.2).and.(beta_mean_tof.lt.2.0)) then + + icount1=0 + icount2=0 + icount3=0 + + do jj=0,23 + a1 = adca(2*jj+1) + a2 = adca(2*jj+2) + if ((a1.lt.4095).and.(a2.lt.4095)) then + a1 = adca(2*jj+1)*cos(theta) + a2 = adca(2*jj+2)*cos(theta) + xhelp = 100000. + xhelp1 = 100000. + xhelp = sqrt(a1*a2) ! geometric mean + xhelp1 = beta_help*xhelp +C if geometric mean multiplied by beta_help is inside/outside helium +C limits, increase counter + if (xhelp1.lt.A_l(jj+1)) icount1=icount1+1 + if ((xhelp1.gt.A_l(jj+1)).and.(xhelp1.lt.A_h(jj+1))) + & icount2=icount2+1 + if (xhelp1.gt.A_h(jj+1)) icount3=icount3+1 + endif + enddo + + +C if more than three paddles see the same... + + if (icount1 .gt. 3) charge=1 + if (icount2 .gt. 3) charge=2 + if (icount3 .gt. 3) charge=3 + + endif ! 0.2 not used +C--- top or bottom artificial: weight*sqrt(2) +C--- top and bottom artificial: weight*sqrt(2)*sqrt(2) + + DO jj=1,12 + if (jj.le.4) xhelp = 0.11 ! S1-S3 + if ((jj.gt.4).and.(jj.le.8)) xhelp = 0.18 ! S2-S3 + if (jj.gt.8) xhelp = 0.28 ! S1-S2 + if ((w_il(itop(jj)).eq.1000.).and.(w_il(ibot(jj)).eq.1000.)) + & xhelp = 1.E09 + if ((w_il(itop(jj)).eq.1).or.(w_il(ibot(jj)).eq.1.)) + & xhelp = xhelp*1.414 + if ((w_il(itop(jj)).eq.1).and.(w_il(ibot(jj)).eq.1.)) + & xhelp = xhelp*2. + w_i(jj) = 1./xhelp + ENDDO + +C======================================================================== +C--- Calculate mean beta for the first time ----------------------------- +C--- We are using "1/beta" since its error is gaussian ------------------ + + icount=0 + sw=0. + sxw=0. + beta_mean=100. + + DO jj=1,12 + IF ((abs(1./b(jj)).gt.0.1).and.(abs(1./b(jj)).lt.15.)) THEN + icount = icount+1 + sxw = sxw + (1./b(jj))*w_i(jj)*w_i(jj) + sw = sw + w_i(jj)*w_i(jj) + ENDIF + ENDDO + + if (icount.gt.0) beta_mean=1./(sxw/sw) + beta_mean_inv = 1./beta_mean + + +C--- Calculate beta for the second time, use residuals of the single +C--- measurements to get a chi2 value + + icount = 0 + sw = 0. + sxw = 0. + betachi = 100. + chi2 = 0. + quality = 0. + + DO jj=1,12 + 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)) ; + if (abs(res*w_i(jj)).lt.resmax) THEN + chi2 = chi2 + (res*w_i(jj))**2. + icount = icount+1 + sxw = sxw + (1./b(jj))*w_i(jj)*w_i(jj) + sw = sw + w_i(jj)*w_i(jj) + ENDIF + ENDIF + ENDDO + +c quality = sw + quality = sqrt(sw) + + if (icount.eq.0) chi2 = 1000. + if (icount.gt.0) chi2 = chi2/(icount) + + if (icount.gt.0) betachi=1./(sxw/sw); + + beta_mean=100. + if ((chi2.lt.chi2cut).and.(quality.gt.qualitycut)) + & beta_mean = betachi + newbeta = beta_mean + + END + +C**************************************************************************** +