Skip to content
Snippets Groups Projects
Commit e6c83679 authored by Hans Rabus's avatar Hans Rabus
Browse files

Upload New File

parent 1ce75f36
No related branches found
No related tags found
1 merge request!1Add new directory
; Read one data file from the PTra simulations
function get_one,file,xpos,zpos,eventf
line=strarr(file_lines(file)) & openr,lun,/get,file & readf,lun,line & free_lun,lun
nhline=where(strpos(line,'GEANT4') gt 0,nw)
if nw gt 0 then begin
eventf=strsplit(line[nhline],'= ',/extract)
; eventf=float(eventf[10+3*indgen(7)])/float(eventf[4]) ; 23.07.2024
eventf=float(eventf[10+3*indgen(7)])/1e8
endif
nhline=max(where(strmid(line,0,1) eq '#')) & line=line[nhline+1:*]
dd=float([[strmid(line,2,7)],[strmid(line,9,7)],[strmid(line,16,12)]])
; 23.07.2024 >>>>>>>>>>>>>>>>>>>>>>>
if strpos(file,'other') ne -1 then begin
eventfile=find_file(strmid(file,0,strpos(file,'PTra'))+'Event_Summary_other.txt',/full,count=nef)
if nef gt 0 then begin
lines=strarr(9)
openr,lun,/get,eventfile & readf,lun,lines & free_lun,lun
line=strsplit(lines[2],' ',/extra) & evtot=float(line[3])
line=strsplit(lines[4],' ',/extra) & evtva=float(line[3])
eventf[4] = evtva/1e8 ; use events through 16.5 nmm aperture
dd[*,2] *= evtot/evtva ; correct results accordingly
endif
endif
; 23.04.2024 <<<<<<<<<<<<<<<<<<<<<<<
uni=uniq(dd[*,0]) & xpos=dd[uni,0] & zpos=dd[0:uni[0],1]
return,reform(dd[*,2],n_elements(zpos),n_elements(xpos))
end ;function
; construct a header line for output files
function headline,tabstart=tabstart,param,what
if keyword_set(tabstart) then txt=string(9b) else txt=''
if vtype(param) eq 'FLOAT' then thick=string(param,f='(f11.7)') else thick=param
thick=strtrim(thick,2)
if n_params() eq 1 or what eq '' then begin
bra='' & ket=''
endif else begin
bra='[' & ket=']'
endelse
txt += strtrim(what+bra+thick[0]+ket,2)
for j=1,n_elements(thick)-1 do txt += string(9b)+strtrim(what+bra+thick[j]+ket,2)
return,txt
end ; function headline
; Write synopses of M1(d) and M1(0) to output files
pro output_M1_vs_d,subdir,which,thick,events,xpos,dat,unc,yplot=yplot,the_case_minus_1
nthick=n_elements(thick)
; write results for M1(d) -------------------------------------------
openw,lun,/get,subdir+'PTra4HIT_M1(d)'+which+'.dat'
printf,lun,headline(thick,'ev',/tab)
printf,lun,headline(events,'',/tab)
printf,lun,string('x / mm',f='(a11)')+headline(thick,'M1',/tab)+headline(thick,'uM1',/tab)
txt=string(xpos,f='(f11.5)')
for j=0,nthick-1 do txt=txt+string(9b)+string(dat[*,j],f='(f11.5)')
for j=0,nthick-1 do txt=txt+string(9b)+string(unc[*,j],f='(f11.5)')
printf,lun,transpose(txt)
free_lun,lun
print,f='(a)','# Output written to '+subdir+'PTra4HIT_M1(d)'+which+'.dat '
; write results for M1(d) in 1 mm steps (for comparison with experiment)
openw,lun,/get,subdir+'PTra4HIT_M1(d)'+which+'_1mm.dat'
printf,lun,headline(thick,'ev',/tab)
printf,lun,headline(events,'',/tab)
printf,lun,string('x / mm',f='(a11)')+headline(thick,'M1',/tab)+headline(thick,'uM1',/tab)
printf,lun,transpose(txt[4*indgen(8)])
free_lun,lun
; write results for M1(d) in 1 mm steps (for comparison with experiment)
openw,lun,/get,subdir+'PTra4HIT_M1(d)'+which+'_ave.dat'
printf,lun,headline(thick,'ev',/tab)
printf,lun,headline(events,'',/tab)
printf,lun,string('d / mm',f='(a11)')+headline(thick,'M1',/tab)+headline(thick,'uM1',/tab)
dmean=dblarr(1+n_elements(xpos),nthick) & umean=dmean
nxpos=n_elements(xpos)
for i=0,nxpos-1 do begin
case i of
0: begin & range=indgen(3) & weights=[1.,4.,1.]/6. & end
1: begin & range=indgen(4)-1 & weights=[4.,3.,4.,1.]/12. & end
nxpos-2: begin & range=1-indgen(4) & weights=[1.,3.,3.,1.]/8. & end
nxpos-1: begin & range=-indgen(3) & weights=[1.,4.,1.]/6. & end
else: begin & range=indgen(5)-2 & weights=[1.,4.,2.,4.,1.]/12. & end
endcase
dmean[i,*]=weights#dat[i+range,*]
umean[i,*]=sqrt(weights#unc[i+range,*]^2)
txt=string(xpos[i],f='(f11.5)')
for j=0,nthick-1 do txt=txt+string(9b)+string(dmean[i,j],f='(f11.5)')
for j=0,nthick-1 do txt=txt+string(9b)+string(umean[i,j],f='(f11.5)')
printf,lun,txt
endfor
free_lun,lun
print,f='(20X,a)',subdir+'PTra4HIT_M1(d)'+which+'_ave.dat '
; write results for M1(0) -------------------------------------------
file=subdir+'PTra4HIT_'
if the_case_minus_1 le 1 then begin
file += 'M1(0)'
nnull=0
endif else begin
file += 'M1bgd'
nnull=3
endelse
file += which+'_ave.dat'
openw,lun,/get,file
; printf,lun,headline(thick,'ev')
; printf,lun,headline(events,'')
; printf,lun,headline(thick,'M1')+headline(thick,'uM1',/tab)
printf,lun,'d_PMMA/mm'+string(9b)+'f_events'+headline('0','M1',/tab)+headline('0','uM1',/tab)
for i=0,nthick-1 do begin
txt=thick[i]
txt += string(9b)
txt += string(events[i],f='(f9.7)')
txt += string(9b)
txt += string(dmean[nnull,i],f='(f11.5)')
txt += string(9b)
txt += string(umean[nnull,i],f='(f11.5)')
printf,lun,txt
endfor
free_lun,lun
print,f='(20X,a)',file
; plot M1(d) -------------------------------------------------------
if keyword_set(yplot) then begin
plot,xpos,dat[*,0],/ylog,ps=2
for i=0,nthick-1 do oplot,xpos,dat[*,i],ps=2,col=i+1
endif
end ; ------------------------------------------------------------------
; Write synopses of M1(x) and M1(3) to output files
pro output_M1_vs_x,subdir,which,thick,events,xpos,dat,unc,yplot=yplot,the_case_minus_1
nfiles=n_elements(thick)
if the_case_minus_1 gt 1 then begin
; write results for M1(x) -------------------------------------------
openw,lun,/get,subdir+'PTra4HIT_M1(x)'+which+'.dat'
printf,lun,headline(thick,'ev',/tab)
printf,lun,headline(events,'',/tab)
printf,lun,string('x / mm',f='(a11)')+headline(thick,'M1',/tab)+headline(thick,'uM1',/tab)
txt=string(xpos,f='(f11.5)')
for j=0,nfiles-1 do txt=txt+string(9b)+string(dat[*,j],f='(f11.5)')
for j=0,nfiles-1 do txt=txt+string(9b)+string(unc[*,j],f='(f11.5)')
printf,lun,transpose(txt)
free_lun,lun
print,f='(a)','# Output written to '+subdir+'PTra4HIT_M1(x)'+which+'.dat '
endif
; write results for M1bgd -------------------------------------------
weights=[1.,4.,2.,4.,1.]/12.
if the_case_minus_1 eq 1 then begin
weights /= 10. ; normalize to width of PSD
range=indgen(5)
; impact parameters up to 2 are twice as frequent -> sum them fist and then again
dmean=weights#dat[range,*]+weights#dat[range+4,*]
umean2=weights#unc[range,*]^2+weights#unc[range+4,*]^2
for i=0,6 do begin
dmean += weights#dat[range+4*i,*]
umean2=weights#unc[range+4*i,*]^2
endfor
umean=sqrt(umean2)
endif else begin
range=min(where(xpos eq 3))+indgen(5)-2
dmean=weights#dat[range,*]
umean=sqrt(weights#unc[range,*]^2)
endelse
dmean=reform(dmean)
umean=reform(umean)
events=reform(events)
file=subdir+'PTra4HIT_M1bgd'+which+'_ave.dat'
openw,lun,/get,file
printf,lun,'d_PMMA/mm'+string(9b)+'f_events'+string(9b)+'M1bgd'+string(9b)+'uM1bgd'
txt = thick+string(9b)+string(events,f='(f9.7)')
txt += string(9b)+string(dmean,f='(f11.5)')
txt += string(9b)+string(umean,f='(f11.5)')
printf,lun,txt,f='(a)'
free_lun,lun
print,f='(20X,a)',file
; plot M1(d) -------------------------------------------------------
if keyword_set(yplot) then begin
plot,xpos,dat[*,0],/ylog,ps=2
for i=0,nfiles-1 do oplot,xpos,dat[*,i],ps=2,col=i+1
endif
end ; ------------------------------------------------------------------
; Write synopses of M1(z,d=x) to output files
pro output_M1_vs_z_at_x,subdir,which,thick,events,zpos,M1z,zplot=zplot
nfiles=n_elements(thick)
; write results for M1(0,z) ----------------------------------------'
; openw,lun,/get,subdir+'PTra4HIT_M1(z)'+which+'.dat'
openw,lun,/get,subdir+'PTra4HIT_'+which+'.dat'
printf,lun,headline(thick,'ev',/tab)
printf,lun,headline(events,'',/tab)
printf,lun,string('z / mm',f='(a11)')+headline(thick,'M1',/tab)
txt=string(zpos,f='(f8.3)')
for i=0,nfiles-1 do txt=txt+string(9b)+string(M1z[*,i],f='(f8.5)')
printf,lun,transpose(txt)
free_lun,lun
; print,f='(20X,a)',subdir+'PTra4HIT_M1(z)'+which+'.dat'
print,f='(20X,a)',subdir+'PTra4HIT_'+which+'.dat'
; plot M1(0,z) ---------------------------------------------------------
if keyword_set(zplot) then begin
plot,zpos,M1z[*,0],/ylog,ps=2
for i=0,nfiles-1 do oplot,zpos,M1z[*,i],ps=2,col=i+1
endif
end ; ------------------------------------------------------------------
; Master subroutine to handle the different cases of C or B ion triggering,
; C or B ion missing and events without C and B ions
pro Ptra4MScaseFinal,subdir,date,bor=bor,out=out,other=other,yplot=yplot,zplot=zplot
if not keyword_set(date) then date='*'
; identify which files to read --------------------------------------
the_case_minus_1=keyword_set(bor)+2*keyword_set(out)+4*keyword_set(other)
case the_case_minus_1 of
0: files=file_search(subdir+'*mm','*_HIT_both_ICSD_*'+date+'*.dat',count=nf)
1: files=file_search(subdir+'*mm','*_HIT_bothB_ICSD_*'+date+'*.dat',count=nf)
2: files=file_search(subdir+'*mm','*_HIT_out_ICSD_*'+date+'*.dat',count=nf)
3: files=file_search(subdir+'*mm','*_HIT_outB_ICSD_*'+date+'*.dat',count=nf)
else: files=file_search(subdir+'*mm','*_HIT_other_ICSD_*'+date+'*.dat',count=nf)
endcase
if nf eq 0 then stop
; find the most recent ones -------------------------------------
linux=max(strpos(files,'\')) eq -1
if linux then files=files[uniq(strmid(files,0,max(strpos(files,'/P'))))] $
else files=files[uniq(strmid(files,0,max(strpos(files,'\P'))))]
nfiles=n_elements(files)
;
if nfiles eq 0 then begin & print,'ERROR: not enough files' & stop & end
; sort file names ascending in absorber thickness ---------------'
thick=files
the_date=''
for i=0,n_elements(files)-1 do begin
if linux then pos=strpos(thick[i],'/') else pos=strpos(thick[i],'\')
thick[i]=strmid(thick[i],pos+1)
the_date=max([the_date,strmid(files[i],strlen(files[i])-12,8)])
if pos gt 0 then begin
thick[i]=strmid(thick[i],0,strpos(thick[i],'mm'))
endif else begin
case strmid(thick[i],0,1) of
'b': thick[i]='1'
'g': thick[i]='2'
'p': thick[i]='3'
'r': thick[i]='0'
endcase
endelse
endfor
so=sort(float(thick))
files=files[so] & thick=thick[so]
if strpos(files[0],'2/') gt 0 then begin
for i=0,n_elements(files)-1 do begin
case thick[i] of
'0': thick[i]='red'
'1': thick[i]='blue'
'2': thick[i]='green'
'3': thick[i]='pink'
endcase
endfor
endif
; define the suffices of the output files
; which=strmid(files[0],0,strlen(files[0])-17)
; which=strmid(which,strpos(which,'HIT')+3)
; if which eq '_both_' then which='_bothC_'
; if which eq '_out_' then which='_outC_'
case the_case_minus_1 of
0: which='_C-triggered_'
1: which='_bothB-bgd_'
2: which='_C-outsider_'
3: which='_outB-bgd_'
else: which='_other-bgd_'
endcase
which=which+the_date
; read the data into array "what" ------------------------------
for i=0,n_elements(files)-1 do begin
print,files[i]
if i eq 0 then what=get_one(files[i],xpos,zpos,eventf) $
else begin
what=[[what],[get_one(files[i],xpos,zpos,events)]]
eventf=[[eventf],[events]]
endelse
endfor
what=reform(what,n_elements(zpos),n_elements(xpos),nfiles)
; calculate M1(d) and uncertainties -----------------------------'
dat=dblarr(n_elements(xpos),nfiles)
unc=dat
qz=where(zpos ge 92.0 and zpos le 227.0,nz)
; if the_case_minus_1 le 1 then begin
if the_case_minus_1 le 1 then begin
q0=where(zpos ge 92.0 and zpos le 227.0,nq)
if the_case_minus_1 eq 0 then qx=0 else qx=3
endif else begin
q0=where(zpos ge 128.8 and zpos le 138.8,nq)
dumm=min(abs(xpos-3.),qx)
endelse
if the_case_minus_1 eq 0 then label='M1(z,0)' else label='M1(z,3)'
for i=0,nq-1 do begin
dat += reform(what[q0[i],*,*],n_elements(xpos),nfiles)
unc += reform(what[q0[i],*,*],n_elements(xpos),nfiles)^2
endfor
;
dat=dat/nq
if total(unc/nq-dat^2 lt 0) gt 1 or nq eq 1 then stop
unc=sqrt((unc/nq-dat^2)*nq/(nq-1))
; write results for M1(d) -------------------------------------------
if the_case_minus_1 eq 0 then begin
output_M1_vs_d,subdir,which,thick,eventf[the_case_minus_1,*], $
xpos,dat,unc,yplot=yplot,the_case_minus_1
endif else begin ; 28.07.2024
output_M1_vs_x,subdir,which,thick,eventf[the_case_minus_1,*], $
xpos,dat,unc,yplot=yplot,the_case_minus_1
endelse ; 28.07.2024
; Save the data for M1(z,x) ---------------------------------------
M1z=reform(what[qz,qx,*],nz,nfiles)
if the_case_minus_1 eq 1 then begin
weights=replicate(1.,n_elements(xpos))
for i=1,n_elements(xpos)-1,2 do weights[i]=4.
for i=2,n_elements(xpos)-1,2 do weights[i]=2.
weights[0:7] *=2. & weights[8]=3.
weights /= 120.
for i=0,nfiles-1 do M1z[*,i]=what[qz,*,i]#weights
endif
output_M1_vs_z_at_x,subdir,label+which,thick,eventf[the_case_minus_1,*],zpos,M1z,zplot=zplot
end ; ------------------------------------------------------------------
; Extract from the M1(d) for triggered events the values at integer
; values of d (for direct comparison with experiment)
pro write_trigger_1mm,subdir,date
;
files1=file_search(subdir,'P*M1(d)_C-trigged_*'+date+'*.dat')
files1=files1[where(strpos(strmid(files1,strlen(subdir)),'PTra4H') eq 0)]
stop
;
for i=0,n_elements(files1)-1 do begin
pos=strpos(files1[i],'.dat')
file=strmid(files1[i],0,pos)+'_1mm.dat'
print,f='(20X,a)',file
lines1=strarr(file_lines(files1[i]))
openr,lun,/get,files1[i] & readf,lun,lines1 & free_lun,lun
openw,lun,/get,file
for j=0,2 do printf,lun,lines1[j]
for j=3,n_elements(lines1)-1,4 do printf,lun,lines1[j]
free_lun,lun
endfor
end
; Extract from the M1(d) averaged over 1 mm of impact parameter those
; values at integer values of d (for direct comparison with experiment
pro write_ave_1mm,subdir,date
;
files1=file_search(subdir,'P*M1(d*_'+date+'_ave.dat')
files1=files1[where(strpos(strmid(files1,strlen(subdir)),'PTra4H') eq 0)]
;
for i=0,n_elements(files1)-1 do begin
pos=strpos(files1[i],'.dat')
file=strmid(files1[i],0,pos)+'_1mm.dat'
print,f='(20X,a)',file
lines1=strarr(file_lines(files1[i]))
openr,lun,/get,files1[i] & readf,lun,lines1 & free_lun,lun
openw,lun,/get,file
for j=0,2 do printf,lun,lines1[j]
for j=3,n_elements(lines1)-1,4 do printf,lun,lines1[j]
free_lun,lun
endfor
end
; Merge the M1(d) and M1(0) data of C and B ions hitting or missing
; the trigger detector.
pro write_C_B_unite,subdir,date,what
;
files1=file_search(subdir,'P*M1*'+what+'C*'+date+'*')
q1=where(strpos(files1,'_ave') ne -1,n1)
q1=where(strpos(files1,'_ave') eq -1,n2)
if n2 gt 0 then begin
files1=files1[q1]
pos=strpos(files1[0],what)+strlen(what)+2
files1=files1[uniq(strmid(files1,0,pos))]
endif
if n1 gt 0 then begin
ufiles=file_search(subdir,'P*M1*'+what+'C*'+date+'*_ave.dat')
pos=strpos(ufiles[0],what)+strlen(what)+2
files1=[files1,ufiles[uniq(strmid(ufiles,0,pos))]]
endif
;
files2=file_search(subdir,'P*M1*'+what+'B*'+date+'*')
q1=where(strpos(files2,'_ave') ne -1,n1)
q1=where(strpos(files2,'_ave') eq -1,n2)
if n2 gt 0 then begin
files2=files2[q1]
pos=strpos(files2[0],what)+strlen(what)+2
files2=files2[uniq(strmid(files2,0,pos))]
endif
if n1 gt 0 then begin
ufiles=file_search(subdir,'P*M1*'+what+'B*'+date+'*_ave.dat')
pos=strpos(ufiles[0],what)+strlen(what)+2
files2=[files2,ufiles[uniq(strmid(ufiles,0,pos))]]
endif
;
for i=0,n_elements(files1)-1 do begin
pos=strpos(files1[i],what)+strlen(what)
file=strmid(files1[i],0,pos)+strmid(files1[i],pos+1)
print,f='(20X,a)',file
lines1=strarr(file_lines(files1[i])) & lines2=lines1
openr,lun,/get,files1[i] & readf,lun,lines1 & free_lun,lun
openr,lun,/get,files2[i] & readf,lun,lines2 & free_lun,lun
openw,lun,/get,file
printf,lun,lines1[0]
if strpos(files1[i],'(0)') gt 0 or strpos(files1[i],'bgd') gt 0 then begin
lines1=lines1[1:*] & lines2=lines2[1:*]
x1=1. & u1=0. & x2=1. & u2=0. & ev1=0. & ev2=0. & thick=0.
for j=0,n_elements(lines1)-1 do begin
reads,lines1[j],thick,ev1,x1,u1
reads,lines2[j],thick,ev2,x2,u2
xx=(ev1*x1+ev2*x2)/(ev1+ev2)
uu=sqrt((ev1*u1^2+ev2*u2^2)/(ev1+ev2))
txt=strtrim(thick,2)+string(9b)+string(ev1+ev2,f='(f9.7)')
txt += string(9b)+strtrim(string(xx,f='(F12.6)'),2)
txt += string(9b)+strtrim(string(uu,f='(F12.6)'),2)
printf,lun,txt
endfor
endif else begin $
ev1=float(strsplit(lines1[1],string(9b),/extra))
ev2=float(strsplit(lines2[1],string(9b),/extra))
nfiles=n_elements(ev1)
txt=string(9b)+string(ev1[0]+ev2[0],f='(f11.7)')
for j=1,nfiles-1 do txt += string(9b)+string(ev1[j]+ev2[j],f='(f11.7)')
printf,lun,txt
printf,lun,lines1[2]
lines1=lines1[3:*] & lines2=lines2[3:*]
if strpos(files1[i],'(z)') gt 0 then begin
x1=dblarr(nfiles) & x2=x1
pos=1.
for j=0,n_elements(lines1)-1 do begin
reads,lines1[j],pos,x1
reads,lines2[j],pos,x2
xx=(ev1*x1+ev2*x2)/(ev1+ev2)
txt=string(pos,f='(f11.5)')
for k=0,nfiles-1 do txt += string(9b)+strtrim(string(xx[k],f='(F12.6)'),2)
printf,lun,txt
endfor
endif else begin
x1=dblarr(nfiles) & u1=x1 & x2=x1 & u2=x1
pos=1.
for j=0,n_elements(lines1)-1 do begin
reads,lines1[j],pos,x1,u1
reads,lines2[j],pos,x2,u2
xx=(ev1*x1+ev2*x2)/(ev1+ev2)
uu=sqrt((ev1*u1^2+ev2*u2^2)/(ev1+ev2))
txt=string(pos,f='(f11.5)')
for k=0,nfiles-1 do txt += string(9b)+strtrim(string(xx[k],f='(F12.6)'),2)
for k=0,nfiles-1 do txt += string(9b)+strtrim(string(uu[k],f='(F12.6)'),2)
printf,lun,txt
endfor
endelse
endelse
free_lun,lun
endfor
end
pro write_corrback_get,subdir,trunc,header=header,thick=thick,ev, M1,uM1,file=file
header=''
file=file_search(subdir,trunc,count=nf)
file=file[nf-1]
print,'GET',file
lines=fltarr(4,file_lines(file)-1)
openr,lun,/get,file & readf,lun,header & readf,lun,lines & free_lun,lun
lines=transpose(lines)
thick=lines[*,0] & ev=lines[*,1] & M1=lines[*,2] & uM1=lines[*,3]
end
pro write_corrback,subdir,pcorr
if not keyword_set(pcorr) then pcorr=0.056
if not keyword_set(subdir) then subdir='largeBeamExp1/'
write_corrback_get,subdir,'P*M1(0)_C-triggered_2024*',header=header,thick=thick,ev1,M10CB,uM10
write_corrback_get,subdir,'P*M1bgd_C-outsider*',ev2,M1bCB,uM1b
write_corrback_get,subdir,'P*M1bgd_other*',ev3,M1oth,uM1oth,file=file
write_corrback_get,subdir,'P*M1bgd_bothB*',ev4,M1bB,uM1bB
write_corrback_get,subdir,'P*M1bgd_outB*',ev5,M1oB,uM1oB
nfiles=n_elements(ev1)
;
M1bgd=(M1oth*ev3+M1bB*ev4+M1oB*ev5)/(ev3+ev4+ev5)
uM1bgd=sqrt((uM1oth^2*ev3+uM1bB^2*ev4+uM1oB^2*ev5)/(ev3+ev4+ev5))
; create file for output of M1bgd_no_carbon
file=strmid(file,0,strpos(file,'other'))+'no-carbon_'+strmid(file,strpos(file,'2024'))
print,f='(20X,a)',file
openw,lun,/get,file
printf,lun,'d_PMMA/mm'+string(9b)+'f_events'+string(9b)+'M1bgd'+string(9b)+'uM1bgd'
txt=strtrim(thick,2)+string(9b)+strtrim(string(ev3+ev4+ev5,f='(F12.6)'),2)
txt=txt+string(9b)+strtrim(string(M1bgd,f='(F12.6)'),2)
txt=txt+string(9b)+strtrim(string(sqrt(uM1bgd^2),f='(F12.6)'),2)
printf,lun,f='(a)',txt
free_lun,lun
;
factor_alignment_unc_200um = 1.09091
ufactor_alignment_unc_200um = 0.03728
urel_pcorr=sqrt(1./12.)
;
coincid=(ev2+ev3+ev4+ev5)/ev1*pcorr*factor_alignment_unc_200um
ucoinci=coincid*sqrt(ufactor_alignment_unc_200um^2+urel_pcorr^2)
M1bgd=(M1bCB*ev2+M1oth*ev3+M1bB*ev4+M1oB*ev5)/(ev2+ev3+ev4+ev5)
uM1bgd=sqrt((uM1b^2*ev2+uM1oth^2*ev3+uM1bB^2*ev4+uM1oB^2*ev5)/(ev2+ev3+ev4+ev5))
urel_M1bgd=sqrt((uM1bgd/M1bgd)^2+(ucoinci/coincid)^2)
M1bgd *= coincid
uM1bgd = M1bgd*urel_M1bgd
; create file for output of corrected M1(0)
file=file_search(subdir,'P*M1(0)_C-triggered_2024*',count=nf)
file=file[nf-1]
pos=strpos(file,'C-triggered_')+strlen('C-triggered_')
file=strmid(file,0,pos)+'corr_'+strmid(file,pos)
print,f='(20X,a)',file
openw,lun,/get,file
printf,lun,'d_PMMA/mm'+string(9b)+'M1(0)'+string(9b)+'uM1(0)'
txt=strtrim(thick,2)+string(9b)+strtrim(string(M10CB+M1bgd,f='(F12.6)'),2)
txt=txt+string(9b)+strtrim(string(sqrt(uM10^2+uM1bgd^2),f='(F12.6)'),2)
printf,lun,f='(a)',txt
free_lun,lun
; ============================= M1(d) ==============================
files4=file_search(subdir,'P*M1(d)_C-triggered_2024*ave*')
files4=files4[where(strpos(strmid(files4,strlen(subdir)),'PTra4H') eq 0)]
files4=files4[sort(strmid(files4,strpos(files4[0],'.')))]
files4=files4[uniq(strmid(files4,strpos(files4[0],'.')))]
print,files4
n4=n_elements(files4)
; create files for output of background versus d full range (for plotting)
lines4=strarr(file_lines(files4[0]))
openr,lun,/get,files4[0] & readf,lun,lines4 & free_lun,lun
header=lines4[0:2] & lines4=lines4[3:*]
pos=strpos(files4[0],'M1(d)_C-triggered')
file=strmid(files4[0],0,pos)+'M1(d)_bgd_'
pos=strpos(files4[0],'_2024')
file += strmid(files4[0],pos+1)
print,f='(20X,a)',file
openw,lun,/get,file
for j=2,2 do printf,lun,header[j]
for i=0,n_elements(lines4)-1 do begin
reads,lines4[i],xpos
txt = strtrim(string(xpos,f='(F12.3)'),2)
for j=0,nfiles-1 do txt += string(9b)+strtrim(string(M1bgd[j],f='(F12.6)'),2)
for j=0,nfiles-1 do txt += string(9b)+strtrim(string(uM1bgd[j],f='(F12.6)'),2)
printf,lun,txt
endfor
free_lun,lun
; create file for output of corrected M1(d)
pos=strpos(files4[0],'C-triggered_')+strlen('C-triggered_')
file=strmid(files4[0],0,pos)+'corr_'+strmid(files4[0],pos)
print,f='(20X,a)',file
openw,lun,/get,file
for j=0,2 do printf,lun,header[j]
for i=0,n_elements(lines4)-1 do begin
temp=float(strsplit(lines4[i],string(9b),/extract))
temp[1:nfiles] += M1bgd
txt = strtrim(string(temp[0],f='(F12.3)'),2)
for j=1,nfiles do txt += string(9b)+strtrim(string(temp[j],f='(F12.6)'),2)
for j=0,nfiles-1 do txt += string(9b)+strtrim(string(sqrt(uM10[j]^2+uM1bgd[j]^2),f='(F12.6)'),2)
printf,lun,txt
endfor
free_lun,lun
end
pro Ptra4MSfinal,date,yplot=yplot,zplot=zplot
if not keyword_set(date) then date=''
subdirs=['largeBeamExp1','largeBeamExp2','smallBeamExp1','smallBeamExp2']+'/'
for i=0,0 do begin
subdir=subdirs[i]
f=file_search('largeBeamExp1/','PTra4HIT_M1*')
;if strpos(f[0],'/') eq 0 then sp
; Get triggered events and unite them
Ptra4MScaseFinal,subdir,date,bor=0,out=0,yplot=yplot,zplot=zplot
; Produce output files with 1 mm step of impact parameter
Ptra4MScaseFinal,subdir,date,bor=1,out=0,yplot=yplot,zplot=zplot
;write_C_B_unite,subdir,date,'both'
; Get outlier events and unite them
Ptra4MScaseFinal,subdir,date,bor=0,out=1,yplot=yplot,zplot=zplot
Ptra4MScaseFinal,subdir,date,bor=1,out=1,yplot=yplot,zplot=zplot
;write_C_B_unite,subdir,date,'out'
; Get other events
Ptra4MScaseFinal,subdir,date,other=1,yplot=yplot,zplot=zplot
; Write the background signal and the background-corrected results
write_corrback,subdir,pcorr
write_ave_1mm,subdir,date
endfor
end
\ No newline at end of file
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment