Index: branches/attic/lumilink/kinds.f90 =================================================================== --- branches/attic/lumilink/kinds.f90 (revision 8598) +++ branches/attic/lumilink/kinds.f90 (revision 8599) @@ -1,23 +0,0 @@ -!******************************************** -!* CompHEP version 3.2 * -!******************************************** -!* F90 format by Wolfgang Kilian Jul 19 1999 -!******************************************** -module kinds - implicit none - private - -! Three types of precision. double is the default, usually. - public :: single, double, quadruple - public :: default - - integer, parameter :: single = & - & selected_real_kind (precision(1.), range(1.)) - integer, parameter :: double = & - & selected_real_kind (precision(1._single) + 1, range(1._single) + 1) - integer, parameter :: quadruple = & - & selected_real_kind (precision (1._double) + 1, range (1._double)) - - integer, parameter :: default = double - -end module kinds Index: branches/attic/lumilink/energy_spread_states.f90 =================================================================== --- branches/attic/lumilink/energy_spread_states.f90 (revision 8598) +++ branches/attic/lumilink/energy_spread_states.f90 (revision 8599) @@ -1,91 +0,0 @@ -module energy_spread_states - - use kinds, only: default - use file_utils, only: free_unit - - implicit none - private - - public :: energy_spread_state_create - - - integer, parameter, public :: n_mcs_pcs_max=36 - - type, public :: energy_spread_state - character(len=200), dimension(36) :: single_file_name - integer :: seed - integer :: i_sqrts - integer :: i_photon - integer :: n_single_file - end type energy_spread_state - - - type(energy_spread_state), public :: mcs - - - - -contains - - - subroutine energy_spread_state_create(energy_spread_in_name) - character(len=*), intent(in) :: energy_spread_in_name - integer :: iu - - mcs%single_file_name = " " - call system_clock(mcs%seed) - mcs%i_sqrts = 500 - mcs%i_photon = 1 - mcs%n_single_file = 0 - - - - - iu = free_unit() - open(unit=iu, file=energy_spread_in_name, action='read', status='old') - print *, " just opened iu,energy_spread_in_name=", iu,energy_spread_in_name - call energy_spread_state_read_unit(iu) - close(unit=iu) - print "(' seed=',i12)", mcs%seed - print "(' i_sqrts=',i12)", mcs%i_sqrts - print "(' i_photon=',i12)", mcs%i_photon - print "(' n_single_file=',i12)", mcs%n_single_file - print "(' single_file_name='/(1x,a))", mcs%single_file_name(1:mcs%n_single_file) - end subroutine energy_spread_state_create - - subroutine energy_spread_state_read_unit(u) - integer, intent(in) :: u - character(len=200), dimension(36) :: single_file_name - integer :: seed - integer :: i_sqrts - integer :: i_photon - namelist /energy_spread_input/ single_file_name - namelist /energy_spread_input/ seed - namelist /energy_spread_input/ i_sqrts - namelist /energy_spread_input/ i_photon - - single_file_name=mcs%single_file_name - seed=mcs%seed - i_sqrts=mcs%i_sqrts - i_photon=mcs%i_photon - - print *, " about to do namelist read u=", u - - read(u, nml=energy_spread_input) - - mcs%single_file_name=single_file_name - mcs%seed=seed - mcs%i_sqrts=i_sqrts - mcs%i_photon=i_photon - - mcs%n_single_file=count(mcs%single_file_name.ne." ") - - return - - end subroutine energy_spread_state_read_unit - - - - - -end module energy_spread_states Index: branches/attic/lumilink/pydatr_common.f90 =================================================================== --- branches/attic/lumilink/pydatr_common.f90 (revision 8598) +++ branches/attic/lumilink/pydatr_common.f90 (revision 8599) @@ -1,16 +0,0 @@ -module pydatr_common - - use kinds, only: double - - implicit none - - integer, dimension(6) :: mrpy - real(kind=double), dimension(100) :: rrpy - common/pydatr/mrpy,rrpy - - -end module pydatr_common - - - - Index: branches/attic/lumilink/pyr_inter.f90 =================================================================== --- branches/attic/lumilink/pyr_inter.f90 (revision 8598) +++ branches/attic/lumilink/pyr_inter.f90 (revision 8599) @@ -1,13 +0,0 @@ - module pyr_inter -! - interface - function pyr(i_dum) - use kinds, only: double - implicit none - real(kind=double) :: pyr - integer, intent(in) :: i_dum - end function pyr - end interface -! - end module pyr_inter - Index: branches/attic/lumilink/Makefile.am =================================================================== --- branches/attic/lumilink/Makefile.am (revision 8598) +++ branches/attic/lumilink/Makefile.am (revision 8599) @@ -1,60 +0,0 @@ -## Makefile.am -- Makefile for WHIZARD -## -## Process this file with automake to produce Makefile.in -# -# Copyright (C) 1999-2014 by -# Wolfgang Kilian -# Thorsten Ohl -# Juergen Reuter -# with contributions from -# Christian Speckner -# -# WHIZARD is free software; you can redistribute it and/or modify it -# under the terms of the GNU General Public License as published by -# the Free Software Foundation; either version 2, or (at your option) -# any later version. -# -# WHIZARD is distributed in the hope that it will be useful, but -# WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with this program; if not, write to the Free Software -# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. -# -######################################################################## - -## If (parts of) LHAPDF are not available, we create a dummy to stub them. -noinst_LTLIBRARIES = liblumilink.la -liblumilink_la_SOURCES = \ - a6f_driver.f90 \ - energy_spread_states.f90 \ - espntp_common.f90 \ - ini_evt_ptr.f90 \ - main.f90 \ - pydatr_common.f90 \ - pur_inter.f90 \ - user.f90 - -AM_FCFLAGS = - -######################################################################## -## Default Fortran compiler options - -## Profiling -if FC_USE_PROFILING -AM_FCFLAGS += $(FCFLAGS_PROFILING) -endif - -## OpenMP -if FC_USE_OPENMP -AM_FCFLAGS += $(FCFLAGS_OPENMP) -endif - -######################################################################## -## Non-standard cleanup tasks - -## Remove backup files -maintainer-clean-local: - -rm -f *~ Index: branches/attic/lumilink/user.f90 =================================================================== --- branches/attic/lumilink/user.f90 (revision 8598) +++ branches/attic/lumilink/user.f90 (revision 8599) @@ -1,913 +0,0 @@ -!======================================================================== -! File: user.f90.default -! -! This is a template for user-defined functions to be called by WHIZARD. -! Copy this code to 'user.f90' and modify it according to your needs. -! The code will automatically be compiled and included in the WHIZARD -! executable. -! -! There are four places where user code can be inserted: -! 1) User-defined beam spectra -! 2) User-defined structure functions -! 3) User-defined cut algorithms -! 4) User-defined reweighting functions for the squared matrix element -! 5) User-defined fragmentation -!======================================================================== - -module user - - use kinds, only: default, double, single !NODEP! - use kinds, only: i32, i64 !NODEP! - use file_utils, only: free_unit !NODEP! - use limits, only: PYTHIA_PARAMETERS_LEN - use lorentz - use particles - use diagnostics, only: msg_level - implicit none - private - - public :: spectrum_single, spectrum_double - public :: spectrum_prt_out, spectrum_beam_remnant - - public :: strfun_single, strfun_double - public :: strfun_prt_out, strfun_beam_remnant - - public :: pdf_init, pdf_strfun_array - - public :: cut - public :: weight - - public :: fragment_init, fragment_call, fragment_end - - integer :: ncall_spectrum=0 - integer :: ncall_spectrum_overflow=0 - - !character(len=100) :: lumi_ee_file='/tmp/lumi_linker_directory/lumi_ee_linker_000' - !character(len=100) :: lumi_eg_file='/tmp/lumi_linker_directory/lumi_eg_linker_000' - !character(len=100) :: lumi_ge_file='/tmp/lumi_linker_directory/lumi_ge_linker_000' - !character(len=100) :: lumi_gg_file='/tmp/lumi_linker_directory/lumi_gg_linker_000' - !character(len=100) :: lumi_file='/tmp/lumi_linker_directory/lumi_linker_000' - !character(len=100) :: photons_beam1_file='/tmp/lumi_linker_directory/photons_beam1_linker_000' - !character(len=100) :: photons_beam2_file='/tmp/lumi_linker_directory/photons_beam2_linker_000' - !character(len=100) :: ebeam_in_file='/tmp/lumi_linker_directory/ebeam_in_linker_000' - !character(len=100) :: pbeam_in_file='/tmp/lumi_linker_directory/pbeam_in_linker_000' - - character(len=200) :: lumi_ee_file - character(len=200) :: lumi_eg_file - character(len=200) :: lumi_ge_file - character(len=200) :: lumi_gg_file - character(len=200) :: lumi_file - character(len=200) :: photons_beam1_file - character(len=200) :: photons_beam2_file - character(len=200) :: ebeam_in_file - character(len=200) :: pbeam_in_file - - real(kind=double) :: yy_electron_peak - real(kind=double) :: photons_beam1_factor - real(kind=double) :: photons_beam2_factor - - integer :: ndiv_eg_lumi - real(kind=double), dimension(:,:), pointer :: yl_eg_lumi - real(kind=double), dimension(:,:), pointer :: yu_eg_lumi - real(kind=double), dimension(:,:), pointer :: rho_single_eg_lumi - real(kind=double), dimension(:,:), pointer :: rho_corr_eg_lumi - real(kind=double), dimension(2) :: avg_energy_eg_lumi - - integer :: ndiv_ge_lumi - real(kind=double), dimension(:,:), pointer :: yl_ge_lumi - real(kind=double), dimension(:,:), pointer :: yu_ge_lumi - real(kind=double), dimension(:,:), pointer :: rho_single_ge_lumi - real(kind=double), dimension(:,:), pointer :: rho_corr_ge_lumi - real(kind=double), dimension(2) :: avg_energy_ge_lumi - - integer :: ndiv_gg_lumi - real(kind=double), dimension(:,:), pointer :: yl_gg_lumi - real(kind=double), dimension(:,:), pointer :: yu_gg_lumi - real(kind=double), dimension(:,:), pointer :: rho_single_gg_lumi - real(kind=double), dimension(:,:), pointer :: rho_corr_gg_lumi - real(kind=double), dimension(2) :: avg_energy_gg_lumi - - integer :: ndiv_lumi - real(kind=double), dimension(:,:), pointer :: yl_lumi - real(kind=double), dimension(:,:), pointer :: yu_lumi - real(kind=double), dimension(:,:), pointer :: rho_single_lumi - real(kind=double), dimension(:,:), pointer :: rho_corr_lumi - real(kind=double), dimension(2) :: avg_energy_lumi - - integer :: ndiv_photons_beam1 - real(kind=double), dimension(:), pointer :: yl_photons_beam1 - real(kind=double), dimension(:), pointer :: yu_photons_beam1 - real(kind=double), dimension(:), pointer :: rho_single_photons_beam1 - real(kind=double) :: avg_energy_photons_beam1 - - integer :: ndiv_photons_beam2 - real(kind=double), dimension(:), pointer :: yl_photons_beam2 - real(kind=double), dimension(:), pointer :: yu_photons_beam2 - real(kind=double), dimension(:), pointer :: rho_single_photons_beam2 - real(kind=double) :: avg_energy_photons_beam2 - - integer :: ndiv_ebeam_in - real(kind=double), dimension(:), pointer :: yl_ebeam_in - real(kind=double), dimension(:), pointer :: yu_ebeam_in - real(kind=double), dimension(:), pointer :: rho_single_ebeam_in - real(kind=double) :: avg_energy_ebeam_in - - integer :: ndiv_pbeam_in - real(kind=double), dimension(:), pointer :: yl_pbeam_in - real(kind=double), dimension(:), pointer :: yu_pbeam_in - real(kind=double), dimension(:), pointer :: rho_single_pbeam_in - real(kind=double) :: avg_energy_pbeam_in - - ! module variables used by the structure function part - integer :: pdf_code_in, pdf_ngroup, pdf_nset, pdf_nfl, pdf_lo - real(kind=double) :: pdf_QCDL4, pdf_top_mass - - logical :: old_style_lumi_linker - -contains - - !======================================================================== - - ! While the physical meaning of a 'spectrum' and a 'structure function' - ! is different, the implementation is exactly analogous. - ! The difference is that a user spectrum will be applied BEFORE all built-in - ! spectra and structure functions (if any), while a user structure functions - ! will be applied AFTER all built-ins. - ! - ! Four functions have to be coded by the user: - ! spectrum_single, spectrum_double, spectrum_prt_out, spectrum_beam_remnant - ! Except for 'spectrum_single', the templates below can be left unchanged, - ! if the default behavior is appropriate. - ! - ! The user spectrum functions take two parameters which can be set in - ! the WHIZARD input file: - ! sqrts : default real (if not set, the total c.m. energy of the process) - ! mode : integer (to distinguish among different user spectra) - ! In the input file, these are called USER_spectrum_sqrts and - ! USER_spectrum_mode, resp. USER_strfun_sqrts and USER_strfun_mode - - !------------------------------------------------------------------------ - ! Return the outgoing particle PDG code, given the incoming particle code - ! and the mode. Note that named constants for PDG codes are defined in the - ! 'particles.f90' module and are available here. - ! For a beam spectrum, this is usually trivial. - - function spectrum_prt_out (prt_in, mode) result (prt_out) - integer, intent(in) :: prt_in, mode - integer :: prt_out - select case (mode) - case default - prt_out = prt_in - end select - end function spectrum_prt_out - - !------------------------------------------------------------------------ - ! Return the PDG code of the beam remnant, given the incoming particle code - ! and the mode. - ! For a beam spectrum, normally there is no beam remnant, and we return - ! UNDEFINED (=0) - - function spectrum_beam_remnant (prt_in, mode) result (prt_out) - integer, intent(in) :: prt_in, mode - integer :: prt_out - check_prt_in: if(prt_in.eq.22) then - check_mode: if(mode.gt.0) then - prt_out=11 - else check_mode - prt_out=-11 - end if check_mode - else check_prt_in - prt_out = UNDEFINED - end if check_prt_in - if (msg_level > 2 ) print *, " spectrum_beam_remnant prt_in,mode,prt_out=", prt_in,mode,prt_out - end function spectrum_beam_remnant - - !----------------------------------------------------------------------- - ! The correlated spectrum for a beam pair. - ! Given the two input parameters x(1), x(2), which are uniformly distributed - ! between 0 and 1, the spectrum function returns output x values together - ! with 'factor', which is the spectrum times the Jacobian of the x mapping, - ! if any. - ! - ! The numbers returned may also depend on the helicities of the incoming - ! and outgoing particles. These are specified in the form of complex density - ! matrices rho_in(:,:,i) and rho_out(:,:,i), where i is the beam index. - ! - ! Each density matrix is a two-dimensional array, where the diagonal elements - ! should be real, and the trace should be normalized to unity. For states - ! with definite helicity, only the diagonal elements are set. Off-diagonal - ! elements refer to superpositions of helicity states, e.g., transversal - ! polarization of fermions. - ! The norm of the spectrum, i.e., the spectrum averaged over initial - ! helicities and summed over outgoing helicities, should be returned - ! as the overall 'factor'. - ! The non-zero entries in the density matrices are defined as +-1 for fermions, - ! +-1 for massless vector bosons, +-1 and 0 for massive vector bosons, and 0 - ! for scalars. Note that the correct assignment is not checked. - ! - ! If the spectrum does not depend on polarization, one should set - ! rho_out = rho_in if the in- and out-particles are identical, resp. - ! rho_out = unpolarized_density_matrix (prt_out) if they are not. - ! The function 'unpolarized_density_matrix' is defined in the module - ! 'particles', as are the PDG codes of the particles. - ! - ! For factorized spectra, the code below is appropriate. Modify it if - ! spectra are correlated, or if you want to improve the performance by - ! mapping the unit square onto itself. - - subroutine spectrum_double (factor, x, prt_in, sqrts, rho_in, rho_out, mode) - real(kind=default), intent(out) :: factor - real(kind=default), dimension(2), intent(inout) :: x - integer, dimension(2), intent(in) :: prt_in - real(kind=default), dimension(2), intent(in) :: sqrts - complex(kind=default), dimension(-2:2,-2:2,2), intent(in) :: rho_in - complex(kind=default), dimension(-2:2,-2:2,2), intent(out) :: rho_out - integer, dimension(2), intent(in) :: mode - real(kind=default), dimension(2) :: x_in - integer, dimension(2) :: ii - real(kind=default), dimension(2) :: factor_single - integer :: i - check_ncall_overflow: if(ncall_spectrum.lt.huge(0)) then - ncall_spectrum=ncall_spectrum+1 - elseif(ncall_spectrum_overflow.eq.0) then check_ncall_overflow - ncall_spectrum_overflow=1 - if (msg_level > 0 ) print *, " encountered ncall_spectrum_overflow; ncall_spectrum,huge(0)=", ncall_spectrum,huge(0) - end if check_ncall_overflow - if(ncall_spectrum.eq.1) call spectrum_ini(mode(1)) - x_in=x - check_prt_in: if(prt_in(1).eq.11.and.prt_in(2).eq.-11) then - rho_out=rho_in - ii=min(ndiv_lumi,int(x_in*dble(ndiv_lumi))+1) - x(1)=yl_lumi(1,ii(1))+(x_in(1)*dble(ndiv_lumi)-dble(ii(1)-1))*(yu_lumi(1,ii(1))-yl_lumi(1,ii(1))) - x(2)=yl_lumi(2,ii(2))+(x_in(2)*dble(ndiv_lumi)-dble(ii(2)-1))*(yu_lumi(2,ii(2))-yl_lumi(2,ii(2))) - x=x/avg_energy_ebeam_in - factor=rho_corr_lumi(ii(1),ii(2))/rho_single_lumi(1,ii(1))/rho_single_lumi(2,ii(2)) - check_ncall_print_ee: if(ncall_spectrum.le.100 .and. msg_level > 3 ) then - print *, " " - print *, " ncall,ii,prt_in,x_in,x,factor=", ncall_spectrum,ii,prt_in,x_in,x,factor - print *, " ncall,yl_lumi(1,ii(1)),yu_lumi(1,ii(1))-yl_lumi(1,ii(1)),x_in(1)*dble(ndiv_lumi)-dble(ii(1)-1)=", ncall_spectrum,yl_lumi(1,ii(1)),yu_lumi(1,ii(1))-yl_lumi(1,ii(1)),x_in(1)*dble(ndiv_lumi)-dble(ii(1)-1) - print *, " ncall,yl_lumi(2,ii(2)),yu_lumi(2,ii(2))-yl_lumi(2,ii(2)),x_in(2)*dble(ndiv_lumi)-dble(ii(2)-1)=", ncall_spectrum,yl_lumi(2,ii(2)),yu_lumi(2,ii(2))-yl_lumi(2,ii(2)),x_in(2)*dble(ndiv_lumi)-dble(ii(2)-1) - end if check_ncall_print_ee - else if((.not.old_style_lumi_linker).and.prt_in(1).eq.11.and.prt_in(2).eq.22) then check_prt_in - rho_out=rho_in - ii=min(ndiv_eg_lumi,int(x_in*dble(ndiv_eg_lumi))+1) - x(1)=yl_eg_lumi(1,ii(1))+(x_in(1)*dble(ndiv_eg_lumi)-dble(ii(1)-1))*(yu_eg_lumi(1,ii(1))-yl_eg_lumi(1,ii(1))) - x(2)=yl_eg_lumi(2,ii(2))+(x_in(2)*dble(ndiv_eg_lumi)-dble(ii(2)-1))*(yu_eg_lumi(2,ii(2))-yl_eg_lumi(2,ii(2))) - x=x/avg_energy_ebeam_in - factor=rho_corr_eg_lumi(ii(1),ii(2))/rho_single_eg_lumi(1,ii(1))/rho_single_eg_lumi(2,ii(2)) - else if((.not.old_style_lumi_linker).and.prt_in(1).eq.22.and.prt_in(2).eq.-11) then check_prt_in - rho_out=rho_in - ii=min(ndiv_ge_lumi,int(x_in*dble(ndiv_ge_lumi))+1) - x(1)=yl_ge_lumi(1,ii(1))+(x_in(1)*dble(ndiv_ge_lumi)-dble(ii(1)-1))*(yu_ge_lumi(1,ii(1))-yl_ge_lumi(1,ii(1))) - x(2)=yl_ge_lumi(2,ii(2))+(x_in(2)*dble(ndiv_ge_lumi)-dble(ii(2)-1))*(yu_ge_lumi(2,ii(2))-yl_ge_lumi(2,ii(2))) - x=x/avg_energy_ebeam_in - factor=rho_corr_ge_lumi(ii(1),ii(2))/rho_single_ge_lumi(1,ii(1))/rho_single_ge_lumi(2,ii(2)) - else if((.not.old_style_lumi_linker).and.prt_in(1).eq.22.and.prt_in(2).eq.22) then check_prt_in - rho_out=rho_in - ii=min(ndiv_gg_lumi,int(x_in*dble(ndiv_gg_lumi))+1) - x(1)=yl_gg_lumi(1,ii(1))+(x_in(1)*dble(ndiv_gg_lumi)-dble(ii(1)-1))*(yu_gg_lumi(1,ii(1))-yl_gg_lumi(1,ii(1))) - x(2)=yl_gg_lumi(2,ii(2))+(x_in(2)*dble(ndiv_gg_lumi)-dble(ii(2)-1))*(yu_gg_lumi(2,ii(2))-yl_gg_lumi(2,ii(2))) - x=x/avg_energy_ebeam_in - factor=rho_corr_gg_lumi(ii(1),ii(2))/rho_single_gg_lumi(1,ii(1))/rho_single_gg_lumi(2,ii(2)) - else check_prt_in - factor=1.d0 - loop_singles: do i=1,2 - call spectrum_single(factor_single(i),x(i),prt_in(i),sqrts(i),rho_in(:,:,i),rho_out(:,:,i),mode(i)) - factor=factor*factor_single(i) - end do loop_singles - check_ncall_print_other: if(ncall_spectrum.le.100 .and. msg_level > 3 ) then - print *, " " - print *, " ncall,ii,prt_in(1:2),x_in(1:2),x(1:2),factor_single(1:2),factor=", ncall_spectrum,ii,prt_in,x_in,x,factor_single,factor - end if check_ncall_print_other - end if check_prt_in - if(ncall_spectrum.le.100) print *, " exit from spectrum_double ncall_spectrum= ",ncall_spectrum - - end subroutine spectrum_double - - !----------------------------------------------------------------------- - ! The spectrum for a single beam - ! Given the input parameter x, which is initially uniformly distributed - ! between 0 and 1, the spectrum function returns an output x value together - ! with 'factor', which is the unpolarized spectrum times the Jacobian of the - ! x mapping, if any. - ! The outgoing density matrix must also be set. See the comment for - ! spectrum_double. - ! - ! The code below describes a trivial spectrum which is unity over the - ! whole kinematical range. Modify this according to your needs. - ! For correlated spectra, leave the code below as is and modify the subroutine - ! spectrum_double instead. - - subroutine spectrum_single (factor, x, prt_in, sqrts, rho_in, rho_out, mode) - real(kind=default), intent(out) :: factor - real(kind=default), intent(inout) :: x - integer, intent(in) :: prt_in - real(kind=default), intent(in) :: sqrts - complex(kind=default), dimension(-2:2,-2:2), intent(in) :: rho_in - complex(kind=default), dimension(-2:2,-2:2), intent(out) :: rho_out - integer, intent(in) :: mode - real(kind=default) :: x_in - integer :: ii - check_ncall_overflow: if(ncall_spectrum.lt.huge(0)) then - ncall_spectrum=ncall_spectrum+1 - elseif(ncall_spectrum_overflow.eq.0) then check_ncall_overflow - ncall_spectrum_overflow=1 - if (msg_level > 0 ) print *, " encountered ncall_spectrum_overflow; ncall_spectrum,huge(0)=", ncall_spectrum,huge(0) - end if check_ncall_overflow - if(ncall_spectrum.eq.1) call spectrum_ini(mode) - rho_out=rho_in - x_in=x - check_prt_in: if(abs(prt_in).eq.11) then - ii=min(ndiv_lumi,int(x_in*dble(ndiv_lumi))+1) - check_mode_sign_lumi: if(mode.gt.0) then - x=yl_lumi(1,ii)+(x_in*dble(ndiv_lumi)-dble(ii-1))*(yu_lumi(1,ii)-yl_lumi(1,ii)) - x=x/avg_energy_ebeam_in - else check_mode_sign_lumi - x=yl_lumi(2,ii)+(x_in*dble(ndiv_lumi)-dble(ii-1))*(yu_lumi(2,ii)-yl_lumi(2,ii)) - x=x/avg_energy_ebeam_in - end if check_mode_sign_lumi - factor=1.d0 - check_ncall_print_lumi: if(ncall_spectrum.le.100 .and. msg_level > 3 ) then - print *, " " - print *, " ncall,prt_in,mode,ii,x_in,x,factor=", ncall_spectrum,prt_in,mode,ii,x_in,x,factor - check_mode_print: if(mode.gt.0) then - print *, " ncall,yl_lumi(1,ii),yu_lumi(1,ii)-yl_lumi(1,ii),x_in*dble(ndiv_lumi)-dble(ii-1)=", ncall_spectrum,yl_lumi(1,ii),yu_lumi(1,ii)-yl_lumi(1,ii),x_in*dble(ndiv_lumi)-dble(ii-1) - else check_mode_print - print *, " ncall,yl_lumi(2,ii),yu_lumi(2,ii)-yl_lumi(2,ii),x_in*dble(ndiv_lumi)-dble(ii-1)=", ncall_spectrum,yl_lumi(2,ii),yu_lumi(2,ii)-yl_lumi(2,ii),x_in*dble(ndiv_lumi)-dble(ii-1) - end if check_mode_print - end if check_ncall_print_lumi - else check_prt_in - check_mode_sign_photons: if(mode.gt.0) then - ii=min(ndiv_photons_beam1,int(x_in*dble(ndiv_photons_beam1))+1) - x=yl_photons_beam1(ii)+(x_in*dble(ndiv_photons_beam1)-dble(ii-1))*(yu_photons_beam1(ii)-yl_photons_beam1(ii)) - x=x/avg_energy_ebeam_in - factor=photons_beam1_factor - else check_mode_sign_photons - ii=min(ndiv_photons_beam2,int(x_in*dble(ndiv_photons_beam2))+1) - x=yl_photons_beam2(ii)+(x_in*dble(ndiv_photons_beam2)-dble(ii-1))*(yu_photons_beam2(ii)-yl_photons_beam2(ii)) - x=x/avg_energy_ebeam_in - factor=photons_beam2_factor - end if check_mode_sign_photons - check_ncall_print_photons: if(ncall_spectrum.le.100 .and. msg_level > 3) then - print *, " " - print *, " ncall,prt_in,mode,ii,x_in,x,factor=", ncall_spectrum,prt_in,mode,ii,x_in,x,factor - check_mode_print_photons: if(mode.gt.0) then - print *, " ncall,yl_photons_beam1(ii),yu_photons_beam1(ii)-yl_photons_beam1(ii),x_in*dble(ndiv_photons_beam1)-dble(ii-1)=", & - & ncall_spectrum,yl_photons_beam1(ii),yu_photons_beam1(ii)-yl_photons_beam1(ii),x_in*dble(ndiv_photons_beam1)-dble(ii-1) - else check_mode_print_photons - print *, " ncall,yl_photons_beam2(ii),yu_photons_beam2(ii)-yl_photons_beam2(ii),x_in*dble(ndiv_photons_beam2)-dble(ii-1)=", & - & ncall_spectrum,yl_photons_beam2(ii),yu_photons_beam2(ii)-yl_photons_beam2(ii),x_in*dble(ndiv_photons_beam2)-dble(ii-1) - end if check_mode_print_photons - end if check_ncall_print_photons - end if check_prt_in - end subroutine spectrum_single - - !======================================================================== - - ! The structure function definitions are analogous to spectrum definitions - ! - ! Four functions have to be coded by the user: - ! strfun_single, strfun_double, strfun_prt_out, strfun_beam_remnant - ! Except for 'strfun_single', the templates below can be left unchanged, - ! if the default behavior is appropriate. - ! - ! The user functions take two parameters which can be set in - ! the WHIZARD input file: - ! sqrts : default real (if not set, the total c.m. energy of the process) - ! mode : integer (to distinguish among different user spectra) - ! In the input file, these are called USER_strfun_sqrts and - ! USER_strfun_mode, resp. USER_strfun_sqrts and USER_strfun_mode - - !------------------------------------------------------------------------ - ! Return the outgoing particle PDG code, given the incoming particle code - ! and the mode. Note that named constants for PDG codes are defined in the - ! 'particles.f90' module and are available here. - ! Default: incoming=outgoing (like, e.g., ISR) - - function strfun_prt_out (prt_in, mode) result (prt_out) - integer, intent(in) :: prt_in, mode - integer :: prt_out - select case (mode) - case default - prt_out = prt_in - end select - end function strfun_prt_out - - !------------------------------------------------------------------------ - ! Return the PDG code of the beam remnant, given the incoming particle code - ! and the mode. - ! Insert something sensible here if you want to have it in the event record. - - function strfun_beam_remnant (prt_in, mode) result (prt_out) - integer, intent(in) :: prt_in, mode - integer :: prt_out - select case (mode) - case default - prt_out = UNDEFINED - end select - end function strfun_beam_remnant - - !----------------------------------------------------------------------- - ! The correlated structure functions for a beam pair. - ! Given the two input parameters x(1), x(2), which are uniformly distributed - ! between 0 and 1, the strfun function returns output x values together - ! with 'factor', which is the strfun times the Jacobian of the x mapping, - ! if any. - ! For the helicity treatment, see the comments for spectrum_double above. - ! - ! Given the fact that structure functions are normally factorized, this may - ! be left as is. However, you can improve the performance by - ! mapping the unit square. Mapping should take place only if the last - ! argument 'map' is .true. If 'map' is .false., this function should return - ! the structure-function f(x) without modifying the x values - - subroutine strfun_double & - & (factor, x, prt_in, sqrts, rho_in, rho_out, mode, map) - real(kind=default), intent(out) :: factor - real(kind=default), dimension(2), intent(inout) :: x - integer, dimension(2), intent(in) :: prt_in - real(kind=default), dimension(2), intent(in) :: sqrts - complex(kind=default), dimension(-2:2,-2:2,2), intent(in) :: rho_in - complex(kind=default), dimension(-2:2,-2:2,2), intent(out) :: rho_out - integer, dimension(2), intent(in) :: mode - logical, intent(in) :: map - real(kind=default), dimension(2) :: f - integer :: i - factor = 1 - do i = 1, 2 - call strfun_single (f(i), x(i), prt_in(i), sqrts(i), & - & rho_in(:,:,i), rho_out(:,:,i), mode(i), map) - factor = factor * f(i) - end do - end subroutine strfun_double - - !----------------------------------------------------------------------- - ! The structure function for a single beam/parton - ! Given the input parameter x, which is initially uniformly distributed - ! between 0 and 1, the strfun function returns an output x value together - ! with 'factor', which is the strfun times the Jacobian of the x mapping, - ! if any. - ! - ! The code below describes a trivial structure function which is unity over the - ! whole kinematical range. Modify this according to your needs. - ! - ! If the structure function is strongly peaked (normally at x=0 or x=1) or - ! even infinite there (it still has to be integrable), you may apply a mapping, - ! so the input and output values of x differ. Mapping should take place only - ! if the last argument 'map' is .true. If 'map' is .false., this function - ! should return the structure-function f(x) without modifying the x value - - subroutine strfun_single & - & (factor, x, prt_in, sqrts, rho_in, rho_out, mode, map) - real(kind=default), intent(out) :: factor - real(kind=default), intent(inout) :: x - integer, intent(in) :: prt_in - real(kind=default), intent(in) :: sqrts - complex(kind=default), dimension(-2:2,-2:2), intent(in) :: rho_in - complex(kind=default), dimension(-2:2,-2:2), intent(out) :: rho_out - integer, intent(in) :: mode - logical, intent(in) :: map - integer :: prt_out - prt_out = spectrum_prt_out (prt_in, mode) - select case (mode) - case default ! Trivial function without polarization dependence - factor = 1 - if (prt_in == prt_out) then - rho_out = rho_in - else - rho_out = unpolarized_density_matrix (prt_out) - end if - end select - end subroutine strfun_single - - !======================================================================= - ! This interface allow for inserting user-defined parton distribution - ! functions or PDFs which are not (yet) contained in the standard PDFlib - - ! This interface will be selected if PDF_ngroup is negative in the - ! WHIZARD input file, while a positive value will select the corresponding - ! PDFLIB structure functions. - - ! Available parameters and module variables: - ! code_in : PDG code of the beam particle (e.g., PROTON) - ! ngroup : Absolute value of PDF_ngroup - ! nset : Value of PDF_nset - ! nfl : Number of active flavors (PDF_nfl) - ! lo : QCD order (PDF_lo) - ! QCDL4 : Value of Lambda(QCD)_4 (PDF_QCDL4) - ! top_mass : Top mass value to be used - - ! Initialization call: - subroutine pdf_init (pdg_code_in, ngroup, nset, nfl, lo, QCDL4, top_mass) - integer, intent(in) :: pdg_code_in, ngroup, nset, nfl, lo - real(kind=double), intent(in) :: QCDL4, top_mass - pdf_code_in = pdg_code_in - pdf_ngroup = ngroup - pdf_nset = nset - pdf_nfl = nfl - pdf_lo = lo - pdf_QCDL4 = QCDL4 - pdf_top_mass = top_mass - ! More initialization here - end subroutine pdf_init - - ! Structure function call: - ! Input parameters: - ! x : Bjorken x value - ! scale : Q scale (set by PDF_running_scale or PDF_scale) - ! Output: - ! dxpdf : array of PDF values [x*f(x)] - ! : (tbar,bbar,cbar,sbar,ubar,dbar,g,d,u,s,c,b,t) - ! : -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 - subroutine pdf_strfun_array (x, scale, dxpdf) - real(kind=double), intent(in) :: x, scale - real(kind=double), dimension(-6:6), intent(out) :: dxpdf - dxpdf = 0 ! Modify this as needed - end subroutine pdf_strfun_array - - !======================================================================= - - ! This minimal cut function accepts everything - - function cut (p, code, mode) result (accept) - type(four_momentum), dimension(:), intent(in) :: p - integer, dimension(:), intent(in) :: code - integer, intent(in) :: mode - logical :: accept - select case (mode) - case default - accept = .true. - end select - end function cut - - ! ! This is an example of a cut function which simulates a 'trigger' at - ! ! a hadron collider detector. It may be adapted to your special needs, - ! ! placed instead of the no-op routine above, and activated by setting - ! ! user_cut_mode nonzero in the input file. The value of user_cut_mode - ! ! is transferred as the mode argument of the cut function. - ! ! The routine is written such that it scans the particle codes, therefore - ! ! it applies to any process. This is recommended practice, although it - ! ! is possible to write process-dependent code by explicitly adressing - ! ! particle indices. - ! ! In this sample cut function, the event is accepted if - ! ! (1) all partons are 0.4 units in eta-phi space apart from each other - ! ! (2) There is at least - ! ! - a lepton with pT > 20 or - ! ! - missing pT > 50 - ! ! - a b quark with pT > 80 or - ! ! - a jet (quark/gluon) with pT > 100 GeV - ! ! Note that kinematical functions are available in the module lorentz.f90 - ! ! while the PDG codes E_LEPTON etc. and functions such as [[visible]] - ! ! are defined in particles.f90 - - ! function cut (p, code, mode) result (accept) - ! ! function arguments and result - ! type(four_momentum), dimension(:), intent(in) :: p - ! integer, dimension(:), intent(in) :: code - ! integer, intent(in) :: mode - ! logical :: accept - ! ! local variables - ! integer :: n, i, j - ! logical, dimension (size(code)) :: is_visible - ! ! number of particles - ! n = size (code) - ! ! flag the visible particles in a logical array - ! is_visible = visible (code) - ! ! There is one nontrivial set of cuts for mode=1 - ! select case (mode) - ! case (1) - ! ! check if all particles (except invisible ones) are separated - ! accept = .true. - ! do i=1, n - ! if (is_visible(i)) then - ! do j=i+1, n - ! if (is_visible(j) .and. eta_phi_distance (p(i), p(j)) < 0.4) & - ! & accept = .false. - ! end do - ! end if - ! end do - ! if (.not.accept) return ! No need to check further - ! ! Check for particles with minimum pT - ! accept = .false. - ! do i=1, n - ! ! Don't distinguish particle/antiparticle - ! select case (abs(code(i))) - ! case (E_LEPTON,MU_LEPTON,TAU_LEPTON) - ! if (transverse_momentum (p(i)) > 20) accept = .true. - ! case (B_QUARK) - ! if (transverse_momentum (p(i)) > 80) accept = .true. - ! case (GLUON, GLUON2, U_QUARK, D_QUARK, S_QUARK, C_QUARK) - ! if (transverse_momentum (p(i)) > 100) accept = .true. - ! end select - ! end do - ! if (accept) return ! No need to check further - ! ! Check the total visible and invisible tranverse momentum - ! if (transverse_momentum (sum (p, mask=is_visible)) > 50) then - ! accept = .true. - ! ! else if (transverse_momentum (sum (p, mask=.not.is_visible)) > 80) then - ! ! accept = .true. - ! end if - ! case default - ! ! No operation for other modes - ! accept = .true. - ! end select - ! end function cut - - - !======================================================================== - - ! This weight function returns unity, i.e. no reweighting. - ! Change this code if you need to reweight the matrix element according - ! to a known (positive-semidefinite) function of the particle momenta - ! and activate it by setting user_weight_mode nonzero (the mode) in the - ! input file. - ! You may use the kinematical functions, operators and particle codes from the - ! lorentz.f90 and particles.f90 modules for that purpose. - - function weight (p, code, mode) - type(four_momentum), dimension(:), intent(in) :: p - integer, dimension(:), intent(in) :: code - integer, intent(in) :: mode - real(kind=default) :: weight - select case (mode) - case default - weight = 1 - end select - end function weight - - !======================================================================== - - ! A standard interface for fragmentation has been defined in - ! E. Boos et al., Proc. Les Houches 2001, hep-ph/0109068 - ! The interface consists of two common blocks, HEPRUP and HEPEUP - ! which are used to transfer all relevant information from the - ! Monte Carlo generator to the fragmentation program. - - ! This subroutine can be modified for initializing user fragmentation. - ! The fragmentation method is determined by the input variable - ! user_fragmentation_mode. It can be used to select among different versions - ! of user fragmentation. - ! For convenience, the process energy and the random generator seed are given. - ! The PYTHIA parameters string is also made accessible here and can thus be - ! abused for transferring any further information if PYTHIA is not invoked - ! otherwise. - ! You can assume that the common block HEPRUP has been filled. - ! Use HEPRUP to access all relevant information. - - subroutine fragment_init (mode, sqrts, seed, pythia_parameters) - integer, intent(in) :: mode - real(kind=default), intent(in) :: sqrts - integer, intent(in) :: seed - character(*), intent(in) :: pythia_parameters - select case (mode) - case default - ! (insert code here) - call ilc_fragment_init(seed,pythia_parameters,sqrts) - end select - end subroutine fragment_init - - ! This subroutine can be modified for executing user fragmentation of - ! a generated event. - ! You may access the event counter as a 64-bit integer. - ! You can assume that the common block HEPEUP has been filled. - ! Use HEPEUP to access all relevant information. - ! The result of fragmentation should be stored in the standard HEPEVT common - ! block. - - subroutine fragment_call (mode, event_count) - integer, intent(in) :: mode - integer(i64), intent(inout) :: event_count - select case (mode) - case default - ! (insert code here) - call ilc_fragment_call(event_count) - end select - end subroutine fragment_call - - ! This subroutine can be modified for finalizing user fragmentation. - ! You may access the event counter (64-bit integer) and the integral. - - subroutine fragment_end (mode, event_count, integral) - integer, intent(in) :: mode - integer(i64), intent(in) :: event_count - real(default), intent(in) :: integral - select case (mode) - case default - ! (insert code here) - call ilc_fragment_print(event_count,integral) - end select - end subroutine fragment_end - - subroutine spectrum_ini (mode) - integer, intent(in) :: mode - integer :: iu - integer, dimension(1) :: i1 - integer :: i - integer :: ios - - real(kind=double), dimension(:,:), pointer :: yl_lumi_loc - real(kind=double), dimension(:,:), pointer :: yu_lumi_loc - real(kind=double), dimension(:,:), pointer :: rho_single_lumi_loc - real(kind=double), dimension(:,:), pointer :: rho_corr_lumi_loc - real(kind=double), dimension(2) :: avg_energy_lumi_loc - - - call getenv("LUMI_EE_LINKER",lumi_ee_file) - call getenv("LUMI_EG_LINKER",lumi_eg_file) - call getenv("LUMI_GE_LINKER",lumi_ge_file) - call getenv("LUMI_GG_LINKER",lumi_gg_file) - call getenv("LUMI_LINKER",lumi_file) - call getenv("PHOTONS_B1",photons_beam1_file) - call getenv("PHOTONS_B2",photons_beam2_file) - call getenv("EBEAM",ebeam_in_file) - call getenv("PBEAM",pbeam_in_file) - lumi_file=adjustr(lumi_file) - write(lumi_file(198:200),FMT='(i3.3)') abs(mode) - lumi_file=adjustl(lumi_file) - if ( msg_level > 0 ) print *, " lumi_file=", lumi_file - - lumi_ee_file=adjustr(lumi_ee_file) - write(lumi_ee_file(198:200),FMT='(i3.3)') abs(mode) - lumi_ee_file=adjustl(lumi_ee_file) - print *, " lumi_ee_file=", lumi_ee_file - - lumi_eg_file=adjustr(lumi_eg_file) - write(lumi_eg_file(198:200),FMT='(i3.3)') abs(mode) - lumi_eg_file=adjustl(lumi_eg_file) - print *, " lumi_eg_file=", lumi_eg_file - - lumi_ge_file=adjustr(lumi_ge_file) - write(lumi_ge_file(198:200),FMT='(i3.3)') abs(mode) - lumi_ge_file=adjustl(lumi_ge_file) - print *, " lumi_ge_file=", lumi_ge_file - - lumi_gg_file=adjustr(lumi_gg_file) - write(lumi_gg_file(198:200),FMT='(i3.3)') abs(mode) - lumi_gg_file=adjustl(lumi_gg_file) - print *, " lumi_gg_file=", lumi_gg_file - - - iu=free_unit() - open(unit=iu,file=lumi_file,action='read',status='old',form='unformatted',iostat=ios) - check_old_or_new_style: if(ios.eq.0) then - old_style_lumi_linker=.true. - else check_old_or_new_style - old_style_lumi_linker=.false. - close(iu) - open(unit=iu,file=lumi_ee_file,action='read',status='old',form='unformatted') - end if check_old_or_new_style - read(iu) ndiv_lumi - allocate (rho_corr_lumi(ndiv_lumi,ndiv_lumi)) - allocate (rho_single_lumi(2,ndiv_lumi)) - allocate (yl_lumi(2,ndiv_lumi)) - allocate (yu_lumi(2,ndiv_lumi)) - allocate (rho_corr_lumi_loc(ndiv_lumi,ndiv_lumi)) - allocate (rho_single_lumi_loc(2,ndiv_lumi)) - allocate (yl_lumi_loc(2,ndiv_lumi)) - allocate (yu_lumi_loc(2,ndiv_lumi)) - read(iu) yl_lumi_loc,yu_lumi_loc,rho_single_lumi_loc,rho_corr_lumi_loc,avg_energy_lumi_loc - close(iu) - check_mode_5: if(abs(mode).ne.5) then - yl_lumi=yl_lumi_loc - yu_lumi=yu_lumi_loc - rho_single_lumi=rho_single_lumi_loc - rho_corr_lumi=rho_corr_lumi_loc - avg_energy_lumi=avg_energy_lumi_loc - else check_mode_5 - - yl_lumi(1,:)=yl_lumi_loc(2,:) - yu_lumi(1,:)=yu_lumi_loc(2,:) - rho_single_lumi(1,:)=rho_single_lumi_loc(2,:) - avg_energy_lumi(1)=avg_energy_lumi_loc(2) - - yl_lumi(2,:)=yl_lumi_loc(1,:) - yu_lumi(2,:)=yu_lumi_loc(1,:) - rho_single_lumi(2,:)=rho_single_lumi_loc(1,:) - avg_energy_lumi(2)=avg_energy_lumi_loc(1) - - rho_corr_lumi=transpose(rho_corr_lumi_loc) - - end if check_mode_5 - i1=maxloc(rho_single_lumi(1,:)) - yy_electron_peak=0.5d0*(yl_lumi(1,i1(1))+yu_lumi(1,i1(1))) - if ( msg_level > 3 ) then - print *, " ndiv_lumi,avg_energy_lumi=", ndiv_lumi,avg_energy_lumi - loop_i_lumi: do i=1,ndiv_lumi - print *, " i,yl_lumi(:,i),yu_lumi(:,i),rho_single_lumi(:,i),rho_corr_lumi(i,i)=", i,yl_lumi(:,i),yu_lumi(:,i),rho_single_lumi(:,i),rho_corr_lumi(i,i) - end do loop_i_lumi - print *, " i1,yy_electron_peak=", i1,yy_electron_peak - endif - - ebeam_in_file=adjustr(ebeam_in_file) - write(ebeam_in_file(198:200),FMT='(i3.3)') abs(mode) - ebeam_in_file=adjustl(ebeam_in_file) - if ( msg_level > 0 ) print *, " ebeam_in_file=", ebeam_in_file - open(unit=iu,file=ebeam_in_file,action='read',status='old',form='unformatted') - read(iu) ndiv_ebeam_in - allocate (rho_single_ebeam_in(ndiv_ebeam_in)) - allocate (yl_ebeam_in(ndiv_ebeam_in)) - allocate (yu_ebeam_in(ndiv_ebeam_in)) - read(iu) yl_ebeam_in,yu_ebeam_in,rho_single_ebeam_in,avg_energy_ebeam_in - close(iu) - print *, " ndiv_ebeam_in,avg_energy_ebeam_in=", ndiv_ebeam_in,avg_energy_ebeam_in - loop_i_ebeam_in: do i=1,ndiv_ebeam_in - print *, " i,yl_ebeam_in(i),yu_ebeam_in(i),rho_single_ebeam_in(i)=", i,yl_ebeam_in(i),yu_ebeam_in(i),rho_single_ebeam_in(i) - end do loop_i_ebeam_in - - check_old_or_new_style_photons: if(old_style_lumi_linker) then - - pbeam_in_file=adjustr(pbeam_in_file) - write(pbeam_in_file(198:200),FMT='(i3.3)') abs(mode) - pbeam_in_file=adjustl(pbeam_in_file) - print *, " pbeam_in_file=", pbeam_in_file - open(unit=iu,file=pbeam_in_file,action='read',status='old',form='unformatted') - read(iu) ndiv_pbeam_in - allocate (rho_single_pbeam_in(ndiv_pbeam_in)) - allocate (yl_pbeam_in(ndiv_pbeam_in)) - allocate (yu_pbeam_in(ndiv_pbeam_in)) - read(iu) yl_pbeam_in,yu_pbeam_in,rho_single_pbeam_in,avg_energy_pbeam_in - close(iu) - print *, " ndiv_pbeam_in,avg_energy_pbeam_in=", ndiv_pbeam_in,avg_energy_pbeam_in - loop_i_pbeam_in: do i=1,ndiv_pbeam_in - print *, " i,yl_pbeam_in(i),yu_pbeam_in(i),rho_single_pbeam_in(i)=", i,yl_pbeam_in(i),yu_pbeam_in(i),rho_single_pbeam_in(i) - end do loop_i_pbeam_in - - photons_beam1_file=adjustr(photons_beam1_file) - write(photons_beam1_file(198:200),FMT='(i3.3)') abs(mode) - photons_beam1_file=adjustl(photons_beam1_file) - print *, " photons_beam1_file=", photons_beam1_file - open(unit=iu,file=photons_beam1_file,action='read',status='old',form='unformatted') - read(iu) ndiv_photons_beam1 - allocate (rho_single_photons_beam1(ndiv_photons_beam1)) - allocate (yl_photons_beam1(ndiv_photons_beam1)) - allocate (yu_photons_beam1(ndiv_photons_beam1)) - read(iu) yl_photons_beam1,yu_photons_beam1,rho_single_photons_beam1,avg_energy_photons_beam1 - close(iu) - if ( msg_level > 3 ) then - print *, " ndiv_photons_beam1,avg_energy_photons_beam1=", ndiv_photons_beam1,avg_energy_photons_beam1 - loop_i_photons_beam1: do i=1,ndiv_photons_beam1 - print *, " i,yl_photons_beam1(i),yu_photons_beam1(i),rho_single_photons_beam1(i)=", i,yl_photons_beam1(i),yu_photons_beam1(i),rho_single_photons_beam1(i) - end do loop_i_photons_beam1 - endif - - photons_beam2_file=adjustr(photons_beam2_file) - write(photons_beam2_file(198:200),FMT='(i3.3)') abs(mode) - photons_beam2_file=adjustl(photons_beam2_file) - if ( msg_level > 0 ) print *, " photons_beam2_file=", photons_beam2_file - open(unit=iu,file=photons_beam2_file,action='read',status='old',form='unformatted') - read(iu) ndiv_photons_beam2 - allocate (rho_single_photons_beam2(ndiv_photons_beam2)) - allocate (yl_photons_beam2(ndiv_photons_beam2)) - allocate (yu_photons_beam2(ndiv_photons_beam2)) - read(iu) yl_photons_beam2,yu_photons_beam2,rho_single_photons_beam2,avg_energy_photons_beam2 - close(iu) - if ( msg_level > 3 ) then - print *, " ndiv_photons_beam2,avg_energy_photons_beam2=", ndiv_photons_beam2,avg_energy_photons_beam2 - loop_i_photons_beam2: do i=1,ndiv_photons_beam2 - print *, " i,yl_photons_beam2(i),yu_photons_beam2(i),rho_single_photons_beam2(i)=", i,yl_photons_beam2(i),yu_photons_beam2(i),rho_single_photons_beam2(i) - end do loop_i_photons_beam2 - endif - - - photons_beam1_factor=(avg_energy_ebeam_in-avg_energy_lumi(1))/avg_energy_photons_beam1 - photons_beam2_factor=(avg_energy_pbeam_in-avg_energy_lumi(2))/avg_energy_photons_beam2 - - print *, " photons_beam1_factor,photons_beam2_factor=", photons_beam1_factor,photons_beam2_factor - - else check_old_or_new_style_photons - - open(unit=iu,file=lumi_eg_file,action='read',status='old',form='unformatted') - read(iu) ndiv_eg_lumi - allocate (rho_corr_eg_lumi(ndiv_eg_lumi,ndiv_eg_lumi)) - allocate (rho_single_eg_lumi(2,ndiv_eg_lumi)) - allocate (yl_eg_lumi(2,ndiv_eg_lumi)) - allocate (yu_eg_lumi(2,ndiv_eg_lumi)) - read(iu) yl_eg_lumi,yu_eg_lumi,rho_single_eg_lumi,rho_corr_eg_lumi,avg_energy_eg_lumi - close(iu) - - open(unit=iu,file=lumi_ge_file,action='read',status='old',form='unformatted') - read(iu) ndiv_ge_lumi - allocate (rho_corr_ge_lumi(ndiv_ge_lumi,ndiv_ge_lumi)) - allocate (rho_single_ge_lumi(2,ndiv_ge_lumi)) - allocate (yl_ge_lumi(2,ndiv_ge_lumi)) - allocate (yu_ge_lumi(2,ndiv_ge_lumi)) - read(iu) yl_ge_lumi,yu_ge_lumi,rho_single_ge_lumi,rho_corr_ge_lumi,avg_energy_ge_lumi - close(iu) - - open(unit=iu,file=lumi_gg_file,action='read',status='old',form='unformatted') - read(iu) ndiv_gg_lumi - allocate (rho_corr_gg_lumi(ndiv_gg_lumi,ndiv_gg_lumi)) - allocate (rho_single_gg_lumi(2,ndiv_gg_lumi)) - allocate (yl_gg_lumi(2,ndiv_gg_lumi)) - allocate (yu_gg_lumi(2,ndiv_gg_lumi)) - read(iu) yl_gg_lumi,yu_gg_lumi,rho_single_gg_lumi,rho_corr_gg_lumi,avg_energy_gg_lumi - close(iu) - - end if check_old_or_new_style_photons - - - end subroutine spectrum_ini - - !======================================================================= - -end module user Index: branches/attic/lumilink/file_utils.f90 =================================================================== --- branches/attic/lumilink/file_utils.f90 (revision 8598) +++ branches/attic/lumilink/file_utils.f90 (revision 8599) @@ -1,32 +0,0 @@ -! WHIZARD file utilities -! -! Automatically generated file, do not edit - -module file_utils - implicit none - private - -! Function to allow for generic I/O -! Returns the next free unit, or -1 if it fails - public :: free_unit - - integer, parameter :: MIN_UNIT = 11, MAX_UNIT = 99 - -contains - - function free_unit () result (unit) - integer :: unit - logical :: exists, is_open - integer :: i, status - do i = MIN_UNIT, MAX_UNIT - inquire (unit=i, exist=exists, opened=is_open, iostat=status) - if (status == 0) then - if (exists .and. .not. is_open) then - unit = i; return - end if - end if - end do - unit = -1 - end function free_unit - -end module file_utils Index: branches/attic/lumilink/espntp_common.f90 =================================================================== --- branches/attic/lumilink/espntp_common.f90 (revision 8598) +++ branches/attic/lumilink/espntp_common.f90 (revision 8599) @@ -1,28 +0,0 @@ -module espntp_common - - use kinds, only: single - - implicit none - - real(kind=single) :: x_nm - real(kind=single) :: y_nm - real(kind=single) :: xp_rad - real(kind=single) :: yp_rad - - real(kind=single) :: e1_gev - real(kind=single) :: e2_gev - real(kind=single) :: x_u - real(kind=single) :: y_u - real(kind=single) :: z_u - real(kind=single) :: xp_urad - real(kind=single) :: yp_urad - common/espntp/e1_gev,e2_gev,x_u,y_u,z_u,xp_urad,yp_urad - - character(len=60), parameter :: ch_espntp='e1_gev:R,e2_gev:R,x_u:R,y_u:R,z_u:R,xp_urad:R,yp_urad:R' -! 1234567890123456789012345678901234567890123456789012345 - -end module espntp_common - - - - Index: branches/attic/lumilink/ini_evt_prt.f90 =================================================================== --- branches/attic/lumilink/ini_evt_prt.f90 (revision 8598) +++ branches/attic/lumilink/ini_evt_prt.f90 (revision 8599) @@ -1,519 +0,0 @@ -module ini_evt_prt - - use kinds, only: single, double - use pyr_inter - use pydatr_common - use energy_spread_states - use file_utils, only: free_unit - - implicit none - private - - public :: a6f_ini - - Common/PAWC/H(10000000) - real(kind=single) :: h - - real(kind=double), parameter :: rn_unit=0.001d0 - real(kind=single), parameter :: e_rn_min=1.0 - real(kind=single) :: sqrts - real(kind=single) :: e1_gev - real(kind=single) :: e2_gev - real(kind=single) :: x_mm - real(kind=single) :: y_mm - real(kind=single) :: x_nm - real(kind=single) :: y_nm - real(kind=single) :: x_u - real(kind=single) :: y_u - real(kind=single) :: z_u - real(kind=single) :: betx1 - real(kind=single) :: bety1 - real(kind=single) :: betx2 - real(kind=single) :: bety2 - real(kind=single) :: photons_xp_rad - real(kind=single) :: photons_yp_rad - real(kind=single) :: unknown - real(kind=single) :: elec_xp_urad - real(kind=single) :: elec_yp_urad - real(kind=single) :: posi_xp_urad - real(kind=single) :: posi_yp_urad - real(kind=double), dimension(2) :: avg_energy_sum=(/0.d0,0.d0/) - real(kind=double), dimension(2) :: avg_energy_test=(/0.d0,0.d0/) - real(kind=double), dimension(2) :: avg_energy=(/0.d0,0.d0/) - - real(kind=double) :: rcount - real(kind=double) :: rcount_1 - real(kind=double) :: rcount_2 - real(kind=double) :: rcount_sum - - real(kind=double), dimension(:,:), pointer :: egev - real(kind=double), dimension(:,:), pointer :: rho_corr - real(kind=double), dimension(:,:), pointer :: rho_single - real(kind=double), dimension(:,:), pointer :: yl - real(kind=double), dimension(:,:), pointer :: yu - - logical, dimension(:,:), pointer :: fg_ord - logical :: f_photons - logical :: f_ebeam_in - logical :: f_pbeam_in - logical :: f_lumi - logical, parameter :: do_fits=.true. - ! logical, parameter :: do_fits=.false. - - integer :: ndiv - integer, dimension(:,:), pointer :: ig_ord - ! integer, parameter :: ndiv_max=20 - ! integer, parameter :: ndiv_max=200 - integer, parameter :: ndiv_max=300 - integer, parameter :: nscale=10 - integer :: a1 - integer :: a2 - integer :: a3 - integer :: a4 - integer :: a5 - integer :: a6 - integer :: a7 - integer :: t_int - - -contains - - subroutine a6f_ini - - real(kind=double) :: sqrt_shat - integer :: icycle - integer :: istat - integer :: ios - integer :: ifile - integer :: iu - integer :: iu_hbook - integer :: irec - integer :: nrec - integer, dimension(:), pointer :: nrec_file - integer :: iu_ntup - character(len=200) :: grid_file - character(len=200) :: grid_file_1 - character(len=200) :: grid_file_n - character(len=200) :: hbook_file - character(len=1) :: c_photon - logical :: f_ip - logical :: f_lumi - integer :: ndel - integer :: i - integer :: i1 - integer :: i2 - integer :: j - integer :: j1 - integer :: j2 - ! - call energy_spread_state_create('energy_spread.in') - - mrpy(1)=mcs%seed - mrpy(2)=0 - - iu_hbook=free_unit() - iu=free_unit() - call hlimit(10000000) - sqrts=1.024*real(mcs%i_sqrts) - - write(unit=c_photon,fmt='(i1.1)') mcs%i_photon - - grid_file_1=adjustr(mcs%single_file_name(1)) - grid_file_n=adjustr(mcs%single_file_name(mcs%n_single_file)) - - f_photons=index(grid_file_1,'photons').ne.0 - f_ebeam_in=index(grid_file_1,'beam1').ne.0 - f_pbeam_in=index(grid_file_1,'beam2').ne.0 - f_lumi=index(grid_file_1,'lumi').ne.0 - - - - check_f_photons_grid_file: if(f_photons) then -! grid_file=grid_file_1(11:190)//'_beam'//c_photon//grid_file_1(192:196)//'-'//grid_file_n(193:196)//'.grd' - grid_file=grid_file_1(12:191)//'_beam'//c_photon//grid_file_1(192:196)//'-'//grid_file_n(193:196)//'.grd' - else check_f_photons_grid_file - ! grid_file=grid_file_1(5:196)//'-'//grid_file_n(194:196)//'.grd' - grid_file=grid_file_1(6:196)//'-'//grid_file_n(193:196)//'.grd' - end if check_f_photons_grid_file - print *, " index(grid_file,'guinea-pig')= ", index(grid_file,'guinea-pig') - - grid_file='/nfs/slac/g/lcd/mc/prj/sw/dist/whizard/'//grid_file(index(grid_file,'guinea-pig'):200) - -! grid_file='/nfs/slac/g/lcd/mc/prj/sw/dist/whizard/guinea-pig/ilc_0250_jun08/lumi/ilc250n_lumi_0xxx-0xxx.grd' - grid_file=adjustr(grid_file) - - hbook_file=grid_file(1:196)//'.hbk' - grid_file=adjustl(grid_file) - hbook_file=adjustl(hbook_file) - hbook_file='energy_spread.hbk' - print *, " " - print *, " grid_file=", grid_file - print *, " hbook_file=", hbook_file - print *, " " - - ! if(1.ne.2) stop - - call hropen(iu_hbook,'ANAHBK',hbook_file,'NX',1024,istat) - check_f_hbook: if(f_photons) then - call hbook1(211,' photon beam '//c_photon//' energy ', 100,0.,0.02*sqrts,0.) - call hbook1(212,' photon beam '//c_photon//' energy ', 128,0.,0.500*sqrts,0.) - call hbook1(213,' photon beam '//c_photon//' energy ', 1280,0.,0.500*sqrts,0.) - elseif(f_ebeam_in) then check_f_hbook - call hbook1(211,' e- energy ', 100,0.480*sqrts,0.500*sqrts,0.) - call hbook1(212,' e- energy ', 128,0.,0.500*sqrts,0.) - call hbook1(213,' e- energy ', 1280,0.,0.500*sqrts,0.) - call hbook1(711,' e- xpos (nm) ', 200,-2000.,2000.,0.) - call hbook1(712,' e- ypos (nm) ', 200,-20.,20.,0.) - call hbook1(713,' e- zpos (um) ', 200,-2000.,2000.,0.) - elseif(f_pbeam_in) then check_f_hbook - call hbook1(211,' e+ energy ', 100,0.480*sqrts,0.500*sqrts,0.) - call hbook1(212,' e+ energy ', 128,0.,0.500*sqrts,0.) - call hbook1(213,' e+ energy ', 1280,0.,0.500*sqrts,0.) - call hbook1(711,' e+ xpos (nm) ', 200,-2000.,2000.,0.) - call hbook1(712,' e+ ypos (nm) ', 200,-20.,20.,0.) - call hbook1(713,' e+ zpos (um) ', 200,-2000.,2000.,0.) - elseif(f_lumi) then check_f_hbook - call hbook1(211,' e- energy ', 100,0.480*sqrts,0.500*sqrts,0.) - call hbook1(212,' e- energy ', 128,0.,0.500*sqrts,0.) - call hbook1(213,' e- energy ', 1280,0.,0.500*sqrts,0.) - call hbook1(214,' e- energy ', 1500,1489.,1490.5,0.) - call hbook1(311,' e+ energy ', 100,0.480*sqrts,0.500*sqrts,0.) - call hbook1(312,' e+ energy ', 128,0.,0.500*sqrts,0.) - call hbook1(314,' e+ energy ', 1200,840.,852.,0.) - call hbook1(411,' sqrt(s) ', 300,0.883*sqrts,1.000*sqrts,0.) - call hbook1(412,' sqrt(s) ', 128,0.,1.000*sqrts,0.) - call hbook2(611,' e+ vs e- energy ', 60,0.480*sqrts,0.500*sqrts,60,0.480*sqrts,0.500*sqrts,0.) - call hbook2(612,' e+ vs e- energy ', 60,0.2*sqrts,0.500*sqrts,60,0.2*sqrts,0.500*sqrts,0.) - call hbook1(711,' xpos (nm) ', 200,-2000.,2000.,0.) - call hbook1(712,' ypos (nm) ', 200,-50.,50.,0.) - call hbook1(713,' zpos (um) ', 200,-2000.,2000.,0.) - end if check_f_hbook - - call hidopt(0,'STAT') - call hidopt(0,'INTE') - - allocate (nrec_file(mcs%n_single_file)) - - nrec=0 - nrec_file=0 - - loop_ifile: do ifile=1,mcs%n_single_file - - print *, " ifile=", ifile, " about to open ", mcs%single_file_name(ifile) - open(unit=iu,file=mcs%single_file_name(ifile),action='read',status='old',form='formatted') -! print *, " immediately after opening ", mcs%single_file_name(ifile) - loop_read: do - check_f_read: if(f_lumi) then - read(unit=iu,fmt=*,iostat=ios) e1_gev,e2_gev,x_nm,y_nm,z_u - check_ios_lumi: if(ios.eq.0) then - nrec_file(ifile)=nrec_file(ifile)+1 - nrec=nrec+1 - else check_ios_lumi - exit loop_read - end if check_ios_lumi - elseif(f_ebeam_in.or.f_pbeam_in) then check_f_read - read(unit=iu,fmt=*,iostat=ios) e1_gev,x_mm,y_mm,z_u,unknown -! print *, " ios,ifile,nrec_file(ifile)= ", ios,ifile,nrec_file(ifile) - check_ios_ebeam_in: if(ios.eq.0) then - nrec_file(ifile)=nrec_file(ifile)+1 - nrec=nrec+1 - else check_ios_ebeam_in - exit loop_read - end if check_ios_ebeam_in - elseif(f_photons) then check_f_read - read(unit=iu,fmt=*,iostat=ios) e1_gev,photons_xp_rad,photons_yp_rad - check_ios_photons: if(ios.eq.0) then - nrec_file(ifile)=nrec_file(ifile)+1 - if((e1_gev.gt.0..and.mcs%i_photon.eq.1).or.(e1_gev.lt.0..and.mcs%i_photon.eq.2)) nrec=nrec+1 - else check_ios_photons - exit loop_read - end if check_ios_photons - end if check_f_read - end do loop_read - close(unit=iu) - end do loop_ifile - - ndiv=min(ndiv_max,nrec/nscale) - print *, "ndiv,nrec,nrec_file=", ndiv,nrec,nrec_file - - allocate (egev(2,nrec)) - allocate (ig_ord(2,nrec)) - allocate (fg_ord(2,nrec)) - allocate (rho_corr(ndiv,ndiv)) - allocate (rho_single(2,ndiv)) - allocate (yl(2,ndiv)) - allocate (yu(2,ndiv)) - - - egev=0.d0 - ig_ord=0 - fg_ord=.true. - rho_corr=0.d0 - rho_single=0.d0 - yl=0.d0 - yu=1.d0 - - - irec=0 - - loop_ifile_again: do ifile=1,mcs%n_single_file - - print *, " ifile=", ifile, " about to again open ", mcs%single_file_name(ifile) - open(unit=iu,file=mcs%single_file_name(ifile),action='read',status='old',form='formatted') - loop_read_again: do i=1,nrec_file(ifile) - check_f_read_again: if(f_lumi) then - read(unit=iu,fmt=*,iostat=ios) e1_gev,e2_gev,x_nm,y_nm,z_u - irec=irec+1 - if(e1_gev.gt.e_rn_min) egev(1,irec)=dble(e1_gev)-rn_unit+2.*rn_unit*pyr(0) - if(e2_gev.gt.e_rn_min) egev(2,irec)=dble(e2_gev)-rn_unit+2.*rn_unit*pyr(0) - print *, " irec,egev(1,irec),egev(2,irec)= ", irec,egev(1,irec),egev(2,irec) - sqrt_shat=sqrt((egev(1,irec)+egev(2,irec))**2-(egev(1,irec)-egev(2,irec))**2) - avg_energy_test=avg_energy_test+egev(:,irec) - call hfill(211,real(egev(1,irec)),0.,1.) - call hfill(212,real(egev(1,irec)),0.,1.) - call hfill(213,real(egev(1,irec)),0.,1.) - call hfill(214,real(egev(1,irec)),0.,1.) - call hfill(311,real(egev(2,irec)),0.,1.) - call hfill(312,real(egev(2,irec)),0.,1.) - call hfill(314,real(egev(2,irec)),0.,1.) - call hfill(411,real(sqrt_shat),0.,1.) - call hfill(412,real(sqrt_shat),0.,1.) - call hfill(611,real(egev(1,irec)),real(egev(2,irec)),1.) - call hfill(612,real(egev(1,irec)),real(egev(2,irec)),1.) - call hfill(711,x_nm,0.,1.) - call hfill(712,y_nm,0.,1.) - call hfill(713,z_u,0.,1.) -! if(mod(irec,100).eq.1) print *, "irec,egev(:,irec)=", irec,egev(:,irec) - elseif(f_ebeam_in.or.f_pbeam_in) then check_f_read_again - read(unit=iu,fmt=*,iostat=ios) e1_gev,x_mm,y_mm,z_u,unknown - irec=irec+1 - egev(1,irec)=dble(e1_gev)-rn_unit+2.*rn_unit*pyr(0) - avg_energy_test=avg_energy_test+egev(:,irec) - x_nm=1.e6*x_mm - y_nm=1.e6*y_mm - call hfill(211,real(egev(1,irec)),0.,1.) - call hfill(212,real(egev(1,irec)),0.,1.) - call hfill(213,real(egev(1,irec)),0.,1.) - call hfill(711,x_nm,0.,1.) - call hfill(712,y_nm,0.,1.) - call hfill(713,z_u,0.,1.) - if(mod(irec,100).eq.1) print *, "irec,egev(1,irec)=", irec,egev(1,irec) - elseif(f_photons) then check_f_read_again - read(unit=iu,fmt=*,iostat=ios) e1_gev,photons_xp_rad,photons_yp_rad - check_i_photon: if((e1_gev.gt.0..and.mcs%i_photon.eq.1).or.(e1_gev.lt.0..and.mcs%i_photon.eq.2)) then - irec=irec+1 - egev(1,irec)=abs(dble(e1_gev)) - avg_energy_test=avg_energy_test+egev(:,irec) - call hfill(211,real(egev(1,irec)),0.,1.) - call hfill(212,real(egev(1,irec)),0.,1.) - call hfill(213,real(egev(1,irec)),0.,1.) - if(mod(irec,100).eq.1) print *, "irec,egev(1,irec)=", irec,egev(1,irec) - end if check_i_photon - end if check_f_read_again - end do loop_read_again - - close(unit=iu) - end do loop_ifile_again - - avg_energy_sum=sum(egev,dim=2) - - - - print *, " irec,nrec,avg_energy_sum,avg_energy_test=", irec,nrec,avg_energy_sum,avg_energy_test - - avg_energy=avg_energy_sum/dble(nrec) - - - check_do_fits: if(do_fits) then - loop_order: do i=1,nrec - ig_ord(1,i:i)=minloc(egev(1,:),mask=fg_ord(1,:)) - fg_ord(1,ig_ord(1,i))=.false. - check_f_lumi_order: if(f_lumi) then - ig_ord(2,i:i)=minloc(egev(2,:),mask=fg_ord(2,:)) - fg_ord(2,ig_ord(2,i))=.false. - end if check_f_lumi_order - if(mod(i,100).eq.0) print *, "i,ig_ord(:,i:i)=", i,ig_ord(:,i:i) - end do loop_order - - - - ndel=nrec/ndiv - - print *, " ndel=", ndel - - loop_div: do i=1,ndiv - i1=1+ndel*(i-1) - i2=ndel*i - check_1: if(i.eq.1) then - yl(1,i)=max(0.d0,egev(1,ig_ord(1,i1))-epsilon(1.d0)) - if(f_lumi) yl(2,i)=max(0.d0,egev(2,ig_ord(2,i1))-epsilon(1.d0)) - else check_1 - yl(:,i)=yu(:,i-1) - end if check_1 - - - check_div: if(i.eq.ndiv) then - yu(1,i)=egev(1,ig_ord(1,i2))+epsilon(1.d0) - print *, " i,i1,i2,yl(1,i),egev(1,ig_ord(1,i1)),egev(1,ig_ord(1,i2)),yu(1,i)=", i,i1,i2,yl(1,i),egev(1,ig_ord(1,i1)),egev(1,ig_ord(1,i2)),yu(1,i) - check_f_lumi_div_1: if(f_lumi) then - yu(2,i)=egev(2,ig_ord(2,i2))+epsilon(1.d0) - print *, " i,i1,i2,yl(2,i),egev(2,ig_ord(2,i1)),egev(2,ig_ord(2,i2)),yu(2,i)=", i,i1,i2,yl(2,i),egev(2,ig_ord(2,i1)),egev(2,ig_ord(2,i2)),yu(2,i) - end if check_f_lumi_div_1 - else check_div - yu(1,i)=0.5d0*(egev(1,ig_ord(1,i2))+egev(1,ig_ord(1,i2+1))) - print *, " i,i1,i2,yl(1,i),egev(1,ig_ord(1,i1)),egev(1,ig_ord(1,i2)),egev(1,ig_ord(1,i2+1)),yu(1,i)=", i,i1,i2,yl(1,i),egev(1,ig_ord(1,i1)),egev(1,ig_ord(1,i2)),egev(1,ig_ord(1,i2+1)),yu(1,i) - check_f_lumi_div_2: if(f_lumi) then - yu(2,i)=0.5d0*(egev(2,ig_ord(2,i2))+egev(2,ig_ord(2,i2+1))) - print *, " i,i1,i2,yl(2,i),egev(2,ig_ord(2,i1)),egev(2,ig_ord(2,i2)),egev(2,ig_ord(2,i2+1)),yu(2,i)=", i,i1,i2,yl(2,i),egev(2,ig_ord(2,i1)),egev(2,ig_ord(2,i2)),egev(2,ig_ord(2,i2+1)),yu(2,i) - end if check_f_lumi_div_2 - end if check_div - - rho_single(:,i)=1.d0/(yu(:,i)-yl(:,i)) - - end do loop_div - - rho_single=rho_single/dble(ndiv) - - loop_print_single: do i=1,ndiv - print *, " i,yl(1,i),yu(1,i),rho_single(1,i)=", i,yl(1,i),yu(1,i),rho_single(1,i) - if(f_lumi) print *, " i,yl(2,i),yu(2,i),rho_single(2,i)=", i,yl(2,i),yu(2,i),rho_single(2,i) - end do loop_print_single - - check_f_lumi_rho: if(f_lumi) then - rcount_sum=0.d0 - - loop_div_i: do i=1,ndiv - loop_div_j: do j=1,ndiv - rcount_1=dble(count(egev(1,:).ge.yl(1,i).and.egev(1,:).lt.yu(1,i))) - rcount_2=dble(count(egev(2,:).ge.yl(2,j).and.egev(2,:).lt.yu(2,j))) - rcount=dble(count(egev(1,:).ge.yl(1,i).and.egev(1,:).lt.yu(1,i) & - & .and.egev(2,:).ge.yl(2,j).and.egev(2,:).lt.yu(2,j))) - rcount_sum=rcount_sum+rcount - rho_corr(i,j)=rcount/(yu(1,i)-yl(1,i))/(yu(2,j)-yl(2,j)) - print *, " i,j,rcount_1,rcount_2,rcount,rcount_sum=", i,j,rcount_1,rcount_2,rcount,rcount_sum - end do loop_div_j - end do loop_div_i - - rho_corr=rho_corr/rcount_sum - loop_print_corr_i: do i=1,ndiv - loop_print_corr_j: do j=1,ndiv - print *, " i,j,yl(1,i),yu(1,i),yl(2,j),yu(2,j),rho_corr(i,j),rho_single(1,i)*rho_single(2,j)=", & - & i,j,yl(1,i),yu(1,i),yl(2,j),yu(2,j),rho_corr(i,j),rho_single(1,i)*rho_single(2,j) - end do loop_print_corr_j - end do loop_print_corr_i - - end if check_f_lumi_rho - - - - - open(unit=iu,file=grid_file,action='write',status='replace',form='unformatted') - - write(iu) ndiv - check_f_lumi_write: if(f_lumi) then - write(iu) yl,yu,rho_single,rho_corr,avg_energy - print *, " avg_energy 1,2=", avg_energy - else check_f_lumi_write - write(iu) yl(1,:),yu(1,:),rho_single(1,:),avg_energy(1) - print *, " avg_energy 1=", avg_energy(1) - end if check_f_lumi_write - - close(unit=iu) - - - check_f_hbfun: if(f_photons) then - call hbfun1(1211,' func photon beam '//c_photon//' energy ',100,0,0.02*sqrts,func_eminus) - call hbfun1(1213,' func photon beam '//c_photon//' energy ',1280,0.,0.500*sqrts,func_eminus) - elseif(f_ebeam_in) then check_f_hbfun - call hbfun1(1211,' func e- energy ',100,0.480*sqrts,0.500*sqrts,func_eminus) - call hbfun1(1213,' func e- energy ',1280,0.,0.500*sqrts,func_eminus) - elseif(f_pbeam_in) then check_f_hbfun - call hbfun1(1211,' func e+ energy ',100,0.480*sqrts,0.500*sqrts,func_eminus) - call hbfun1(1213,' func e+ energy ',1280,0.,0.500*sqrts,func_eminus) - elseif(f_lumi) then check_f_hbfun - call hbfun1(1211,' func e- energy ',100,0.480*sqrts,0.500*sqrts,func_eminus) - call hbfun1(1213,' func e- energy ',1280,0.,0.500*sqrts,func_eminus) - call hbfun1(1311,' func e+ energy ',100,0.480*sqrts,0.500*sqrts,func_eplus) - call hbfun2(1611,' func e+ vs e- energy ',60,0.480*sqrts,0.500*sqrts,60,0.480*sqrts,0.500*sqrts,rho_corr_y) - end if check_f_hbfun - - end if check_do_fits - - deallocate (yu) - deallocate (yl) - deallocate (rho_single) - deallocate (rho_corr) - deallocate (ig_ord) - deallocate (fg_ord) - deallocate (egev) - - call hcdir('//ANAHBK',' ') - call hprint(0) - call hrout(0,icycle,' ') - call hrend('ANAHBK') - - call hdelet(0) - - - return - - end subroutine a6f_ini - - function rho_corr_y(yy) result(rho) - real(kind=single), dimension(2), intent(in) :: yy - real(kind=single) :: rho - real(kind=double), dimension(2) :: yd - integer, dimension(1) :: i1 - integer, dimension(1) :: i2 - yd=dble(yy) - check_range: if(yd(1).ge.yl(1,1).and.yd(1).lt.yu(1,ndiv).and.yd(2).ge.yl(2,1).and.yd(2).lt.yu(2,ndiv)) then - i1=minloc(abs(yd(1)-yl(1,:))) - i2=minloc(abs(yd(2)-yl(2,:))) - if(yd(1).lt.yl(1,i1(1))) i1=i1-1 - if(yd(2).lt.yl(2,i2(1))) i2=i2-1 - rho=real(rho_corr(i1(1),i2(1))) - else check_range - rho=0. - end if check_range - return - end function rho_corr_y - - function rho_single_y(i,yy) result(rho) - integer, intent(in) :: i - real(kind=single), intent(in) :: yy - real(kind=single) :: rho - real(kind=double) :: yd - integer, dimension(1) :: i1 - yd=dble(yy) - check_range: if(yd.ge.yl(i,1).and.yd.lt.yu(i,ndiv)) then - i1=minloc(abs(yd-yl(i,:))) - if(yd.lt.yl(i,i1(1))) i1=i1-1 - rho=real(rho_single(i,i1(1))) - else check_range - rho=0. - end if check_range - return - end function rho_single_y - - - - function func_eminus(yy) result(rho) - real(kind=single), intent(in) :: yy - real(kind=single) :: rho - rho=rho_single_y(1,yy) - ! print *, " yy,func_eminus=", yy,rho - return - end function func_eminus - - function func_eplus(yy) result(rho) - real(kind=single), intent(in) :: yy - real(kind=single) :: rho - rho=rho_single_y(2,yy) - ! print *, " yy,func_eplus=", yy,rho - return - end function func_eplus - - - - -end module ini_evt_prt Index: branches/attic/lumilink/main.f90 =================================================================== --- branches/attic/lumilink/main.f90 (revision 8598) +++ branches/attic/lumilink/main.f90 (revision 8599) @@ -1,10 +0,0 @@ -program main - - call a6f_driver - - - stop - -end program main - - Index: branches/attic/lumilink/a6f_driver.f90 =================================================================== --- branches/attic/lumilink/a6f_driver.f90 (revision 8598) +++ branches/attic/lumilink/a6f_driver.f90 (revision 8599) @@ -1,17 +0,0 @@ - subroutine a6f_driver - - use ini_evt_prt - - implicit none - - call a6f_ini - -! call a6f_evt - -! call a6f_prt - - return - - end subroutine a6f_driver - -