Index: applgrid/tags/applgrid-01-06-25/depcomp
===================================================================
--- applgrid/tags/applgrid-01-06-25/depcomp (revision 0)
+++ applgrid/tags/applgrid-01-06-25/depcomp (revision 1980)
@@ -0,0 +1,791 @@
+#! /bin/sh
+# depcomp - compile a program generating dependencies as side-effects
+
+scriptversion=2013-05-30.07; # UTC
+
+# Copyright (C) 1999-2014 Free Software Foundation, Inc.
+
+# This program 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.
+
+# This program 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, see .
+
+# As a special exception to the GNU General Public License, if you
+# distribute this file as part of a program that contains a
+# configuration script generated by Autoconf, you may include it under
+# the same distribution terms that you use for the rest of that program.
+
+# Originally written by Alexandre Oliva .
+
+case $1 in
+ '')
+ echo "$0: No command. Try '$0 --help' for more information." 1>&2
+ exit 1;
+ ;;
+ -h | --h*)
+ cat <<\EOF
+Usage: depcomp [--help] [--version] PROGRAM [ARGS]
+
+Run PROGRAMS ARGS to compile a file, generating dependencies
+as side-effects.
+
+Environment variables:
+ depmode Dependency tracking mode.
+ source Source file read by 'PROGRAMS ARGS'.
+ object Object file output by 'PROGRAMS ARGS'.
+ DEPDIR directory where to store dependencies.
+ depfile Dependency file to output.
+ tmpdepfile Temporary file to use when outputting dependencies.
+ libtool Whether libtool is used (yes/no).
+
+Report bugs to .
+EOF
+ exit $?
+ ;;
+ -v | --v*)
+ echo "depcomp $scriptversion"
+ exit $?
+ ;;
+esac
+
+# Get the directory component of the given path, and save it in the
+# global variables '$dir'. Note that this directory component will
+# be either empty or ending with a '/' character. This is deliberate.
+set_dir_from ()
+{
+ case $1 in
+ */*) dir=`echo "$1" | sed -e 's|/[^/]*$|/|'`;;
+ *) dir=;;
+ esac
+}
+
+# Get the suffix-stripped basename of the given path, and save it the
+# global variable '$base'.
+set_base_from ()
+{
+ base=`echo "$1" | sed -e 's|^.*/||' -e 's/\.[^.]*$//'`
+}
+
+# If no dependency file was actually created by the compiler invocation,
+# we still have to create a dummy depfile, to avoid errors with the
+# Makefile "include basename.Plo" scheme.
+make_dummy_depfile ()
+{
+ echo "#dummy" > "$depfile"
+}
+
+# Factor out some common post-processing of the generated depfile.
+# Requires the auxiliary global variable '$tmpdepfile' to be set.
+aix_post_process_depfile ()
+{
+ # If the compiler actually managed to produce a dependency file,
+ # post-process it.
+ if test -f "$tmpdepfile"; then
+ # Each line is of the form 'foo.o: dependency.h'.
+ # Do two passes, one to just change these to
+ # $object: dependency.h
+ # and one to simply output
+ # dependency.h:
+ # which is needed to avoid the deleted-header problem.
+ { sed -e "s,^.*\.[$lower]*:,$object:," < "$tmpdepfile"
+ sed -e "s,^.*\.[$lower]*:[$tab ]*,," -e 's,$,:,' < "$tmpdepfile"
+ } > "$depfile"
+ rm -f "$tmpdepfile"
+ else
+ make_dummy_depfile
+ fi
+}
+
+# A tabulation character.
+tab=' '
+# A newline character.
+nl='
+'
+# Character ranges might be problematic outside the C locale.
+# These definitions help.
+upper=ABCDEFGHIJKLMNOPQRSTUVWXYZ
+lower=abcdefghijklmnopqrstuvwxyz
+digits=0123456789
+alpha=${upper}${lower}
+
+if test -z "$depmode" || test -z "$source" || test -z "$object"; then
+ echo "depcomp: Variables source, object and depmode must be set" 1>&2
+ exit 1
+fi
+
+# Dependencies for sub/bar.o or sub/bar.obj go into sub/.deps/bar.Po.
+depfile=${depfile-`echo "$object" |
+ sed 's|[^\\/]*$|'${DEPDIR-.deps}'/&|;s|\.\([^.]*\)$|.P\1|;s|Pobj$|Po|'`}
+tmpdepfile=${tmpdepfile-`echo "$depfile" | sed 's/\.\([^.]*\)$/.T\1/'`}
+
+rm -f "$tmpdepfile"
+
+# Avoid interferences from the environment.
+gccflag= dashmflag=
+
+# Some modes work just like other modes, but use different flags. We
+# parameterize here, but still list the modes in the big case below,
+# to make depend.m4 easier to write. Note that we *cannot* use a case
+# here, because this file can only contain one case statement.
+if test "$depmode" = hp; then
+ # HP compiler uses -M and no extra arg.
+ gccflag=-M
+ depmode=gcc
+fi
+
+if test "$depmode" = dashXmstdout; then
+ # This is just like dashmstdout with a different argument.
+ dashmflag=-xM
+ depmode=dashmstdout
+fi
+
+cygpath_u="cygpath -u -f -"
+if test "$depmode" = msvcmsys; then
+ # This is just like msvisualcpp but w/o cygpath translation.
+ # Just convert the backslash-escaped backslashes to single forward
+ # slashes to satisfy depend.m4
+ cygpath_u='sed s,\\\\,/,g'
+ depmode=msvisualcpp
+fi
+
+if test "$depmode" = msvc7msys; then
+ # This is just like msvc7 but w/o cygpath translation.
+ # Just convert the backslash-escaped backslashes to single forward
+ # slashes to satisfy depend.m4
+ cygpath_u='sed s,\\\\,/,g'
+ depmode=msvc7
+fi
+
+if test "$depmode" = xlc; then
+ # IBM C/C++ Compilers xlc/xlC can output gcc-like dependency information.
+ gccflag=-qmakedep=gcc,-MF
+ depmode=gcc
+fi
+
+case "$depmode" in
+gcc3)
+## gcc 3 implements dependency tracking that does exactly what
+## we want. Yay! Note: for some reason libtool 1.4 doesn't like
+## it if -MD -MP comes after the -MF stuff. Hmm.
+## Unfortunately, FreeBSD c89 acceptance of flags depends upon
+## the command line argument order; so add the flags where they
+## appear in depend2.am. Note that the slowdown incurred here
+## affects only configure: in makefiles, %FASTDEP% shortcuts this.
+ for arg
+ do
+ case $arg in
+ -c) set fnord "$@" -MT "$object" -MD -MP -MF "$tmpdepfile" "$arg" ;;
+ *) set fnord "$@" "$arg" ;;
+ esac
+ shift # fnord
+ shift # $arg
+ done
+ "$@"
+ stat=$?
+ if test $stat -ne 0; then
+ rm -f "$tmpdepfile"
+ exit $stat
+ fi
+ mv "$tmpdepfile" "$depfile"
+ ;;
+
+gcc)
+## Note that this doesn't just cater to obsosete pre-3.x GCC compilers.
+## but also to in-use compilers like IMB xlc/xlC and the HP C compiler.
+## (see the conditional assignment to $gccflag above).
+## There are various ways to get dependency output from gcc. Here's
+## why we pick this rather obscure method:
+## - Don't want to use -MD because we'd like the dependencies to end
+## up in a subdir. Having to rename by hand is ugly.
+## (We might end up doing this anyway to support other compilers.)
+## - The DEPENDENCIES_OUTPUT environment variable makes gcc act like
+## -MM, not -M (despite what the docs say). Also, it might not be
+## supported by the other compilers which use the 'gcc' depmode.
+## - Using -M directly means running the compiler twice (even worse
+## than renaming).
+ if test -z "$gccflag"; then
+ gccflag=-MD,
+ fi
+ "$@" -Wp,"$gccflag$tmpdepfile"
+ stat=$?
+ if test $stat -ne 0; then
+ rm -f "$tmpdepfile"
+ exit $stat
+ fi
+ rm -f "$depfile"
+ echo "$object : \\" > "$depfile"
+ # The second -e expression handles DOS-style file names with drive
+ # letters.
+ sed -e 's/^[^:]*: / /' \
+ -e 's/^['$alpha']:\/[^:]*: / /' < "$tmpdepfile" >> "$depfile"
+## This next piece of magic avoids the "deleted header file" problem.
+## The problem is that when a header file which appears in a .P file
+## is deleted, the dependency causes make to die (because there is
+## typically no way to rebuild the header). We avoid this by adding
+## dummy dependencies for each header file. Too bad gcc doesn't do
+## this for us directly.
+## Some versions of gcc put a space before the ':'. On the theory
+## that the space means something, we add a space to the output as
+## well. hp depmode also adds that space, but also prefixes the VPATH
+## to the object. Take care to not repeat it in the output.
+## Some versions of the HPUX 10.20 sed can't process this invocation
+## correctly. Breaking it into two sed invocations is a workaround.
+ tr ' ' "$nl" < "$tmpdepfile" \
+ | sed -e 's/^\\$//' -e '/^$/d' -e "s|.*$object$||" -e '/:$/d' \
+ | sed -e 's/$/ :/' >> "$depfile"
+ rm -f "$tmpdepfile"
+ ;;
+
+hp)
+ # This case exists only to let depend.m4 do its work. It works by
+ # looking at the text of this script. This case will never be run,
+ # since it is checked for above.
+ exit 1
+ ;;
+
+sgi)
+ if test "$libtool" = yes; then
+ "$@" "-Wp,-MDupdate,$tmpdepfile"
+ else
+ "$@" -MDupdate "$tmpdepfile"
+ fi
+ stat=$?
+ if test $stat -ne 0; then
+ rm -f "$tmpdepfile"
+ exit $stat
+ fi
+ rm -f "$depfile"
+
+ if test -f "$tmpdepfile"; then # yes, the sourcefile depend on other files
+ echo "$object : \\" > "$depfile"
+ # Clip off the initial element (the dependent). Don't try to be
+ # clever and replace this with sed code, as IRIX sed won't handle
+ # lines with more than a fixed number of characters (4096 in
+ # IRIX 6.2 sed, 8192 in IRIX 6.5). We also remove comment lines;
+ # the IRIX cc adds comments like '#:fec' to the end of the
+ # dependency line.
+ tr ' ' "$nl" < "$tmpdepfile" \
+ | sed -e 's/^.*\.o://' -e 's/#.*$//' -e '/^$/ d' \
+ | tr "$nl" ' ' >> "$depfile"
+ echo >> "$depfile"
+ # The second pass generates a dummy entry for each header file.
+ tr ' ' "$nl" < "$tmpdepfile" \
+ | sed -e 's/^.*\.o://' -e 's/#.*$//' -e '/^$/ d' -e 's/$/:/' \
+ >> "$depfile"
+ else
+ make_dummy_depfile
+ fi
+ rm -f "$tmpdepfile"
+ ;;
+
+xlc)
+ # This case exists only to let depend.m4 do its work. It works by
+ # looking at the text of this script. This case will never be run,
+ # since it is checked for above.
+ exit 1
+ ;;
+
+aix)
+ # The C for AIX Compiler uses -M and outputs the dependencies
+ # in a .u file. In older versions, this file always lives in the
+ # current directory. Also, the AIX compiler puts '$object:' at the
+ # start of each line; $object doesn't have directory information.
+ # Version 6 uses the directory in both cases.
+ set_dir_from "$object"
+ set_base_from "$object"
+ if test "$libtool" = yes; then
+ tmpdepfile1=$dir$base.u
+ tmpdepfile2=$base.u
+ tmpdepfile3=$dir.libs/$base.u
+ "$@" -Wc,-M
+ else
+ tmpdepfile1=$dir$base.u
+ tmpdepfile2=$dir$base.u
+ tmpdepfile3=$dir$base.u
+ "$@" -M
+ fi
+ stat=$?
+ if test $stat -ne 0; then
+ rm -f "$tmpdepfile1" "$tmpdepfile2" "$tmpdepfile3"
+ exit $stat
+ fi
+
+ for tmpdepfile in "$tmpdepfile1" "$tmpdepfile2" "$tmpdepfile3"
+ do
+ test -f "$tmpdepfile" && break
+ done
+ aix_post_process_depfile
+ ;;
+
+tcc)
+ # tcc (Tiny C Compiler) understand '-MD -MF file' since version 0.9.26
+ # FIXME: That version still under development at the moment of writing.
+ # Make that this statement remains true also for stable, released
+ # versions.
+ # It will wrap lines (doesn't matter whether long or short) with a
+ # trailing '\', as in:
+ #
+ # foo.o : \
+ # foo.c \
+ # foo.h \
+ #
+ # It will put a trailing '\' even on the last line, and will use leading
+ # spaces rather than leading tabs (at least since its commit 0394caf7
+ # "Emit spaces for -MD").
+ "$@" -MD -MF "$tmpdepfile"
+ stat=$?
+ if test $stat -ne 0; then
+ rm -f "$tmpdepfile"
+ exit $stat
+ fi
+ rm -f "$depfile"
+ # Each non-empty line is of the form 'foo.o : \' or ' dep.h \'.
+ # We have to change lines of the first kind to '$object: \'.
+ sed -e "s|.*:|$object :|" < "$tmpdepfile" > "$depfile"
+ # And for each line of the second kind, we have to emit a 'dep.h:'
+ # dummy dependency, to avoid the deleted-header problem.
+ sed -n -e 's|^ *\(.*\) *\\$|\1:|p' < "$tmpdepfile" >> "$depfile"
+ rm -f "$tmpdepfile"
+ ;;
+
+## The order of this option in the case statement is important, since the
+## shell code in configure will try each of these formats in the order
+## listed in this file. A plain '-MD' option would be understood by many
+## compilers, so we must ensure this comes after the gcc and icc options.
+pgcc)
+ # Portland's C compiler understands '-MD'.
+ # Will always output deps to 'file.d' where file is the root name of the
+ # source file under compilation, even if file resides in a subdirectory.
+ # The object file name does not affect the name of the '.d' file.
+ # pgcc 10.2 will output
+ # foo.o: sub/foo.c sub/foo.h
+ # and will wrap long lines using '\' :
+ # foo.o: sub/foo.c ... \
+ # sub/foo.h ... \
+ # ...
+ set_dir_from "$object"
+ # Use the source, not the object, to determine the base name, since
+ # that's sadly what pgcc will do too.
+ set_base_from "$source"
+ tmpdepfile=$base.d
+
+ # For projects that build the same source file twice into different object
+ # files, the pgcc approach of using the *source* file root name can cause
+ # problems in parallel builds. Use a locking strategy to avoid stomping on
+ # the same $tmpdepfile.
+ lockdir=$base.d-lock
+ trap "
+ echo '$0: caught signal, cleaning up...' >&2
+ rmdir '$lockdir'
+ exit 1
+ " 1 2 13 15
+ numtries=100
+ i=$numtries
+ while test $i -gt 0; do
+ # mkdir is a portable test-and-set.
+ if mkdir "$lockdir" 2>/dev/null; then
+ # This process acquired the lock.
+ "$@" -MD
+ stat=$?
+ # Release the lock.
+ rmdir "$lockdir"
+ break
+ else
+ # If the lock is being held by a different process, wait
+ # until the winning process is done or we timeout.
+ while test -d "$lockdir" && test $i -gt 0; do
+ sleep 1
+ i=`expr $i - 1`
+ done
+ fi
+ i=`expr $i - 1`
+ done
+ trap - 1 2 13 15
+ if test $i -le 0; then
+ echo "$0: failed to acquire lock after $numtries attempts" >&2
+ echo "$0: check lockdir '$lockdir'" >&2
+ exit 1
+ fi
+
+ if test $stat -ne 0; then
+ rm -f "$tmpdepfile"
+ exit $stat
+ fi
+ rm -f "$depfile"
+ # Each line is of the form `foo.o: dependent.h',
+ # or `foo.o: dep1.h dep2.h \', or ` dep3.h dep4.h \'.
+ # Do two passes, one to just change these to
+ # `$object: dependent.h' and one to simply `dependent.h:'.
+ sed "s,^[^:]*:,$object :," < "$tmpdepfile" > "$depfile"
+ # Some versions of the HPUX 10.20 sed can't process this invocation
+ # correctly. Breaking it into two sed invocations is a workaround.
+ sed 's,^[^:]*: \(.*\)$,\1,;s/^\\$//;/^$/d;/:$/d' < "$tmpdepfile" \
+ | sed -e 's/$/ :/' >> "$depfile"
+ rm -f "$tmpdepfile"
+ ;;
+
+hp2)
+ # The "hp" stanza above does not work with aCC (C++) and HP's ia64
+ # compilers, which have integrated preprocessors. The correct option
+ # to use with these is +Maked; it writes dependencies to a file named
+ # 'foo.d', which lands next to the object file, wherever that
+ # happens to be.
+ # Much of this is similar to the tru64 case; see comments there.
+ set_dir_from "$object"
+ set_base_from "$object"
+ if test "$libtool" = yes; then
+ tmpdepfile1=$dir$base.d
+ tmpdepfile2=$dir.libs/$base.d
+ "$@" -Wc,+Maked
+ else
+ tmpdepfile1=$dir$base.d
+ tmpdepfile2=$dir$base.d
+ "$@" +Maked
+ fi
+ stat=$?
+ if test $stat -ne 0; then
+ rm -f "$tmpdepfile1" "$tmpdepfile2"
+ exit $stat
+ fi
+
+ for tmpdepfile in "$tmpdepfile1" "$tmpdepfile2"
+ do
+ test -f "$tmpdepfile" && break
+ done
+ if test -f "$tmpdepfile"; then
+ sed -e "s,^.*\.[$lower]*:,$object:," "$tmpdepfile" > "$depfile"
+ # Add 'dependent.h:' lines.
+ sed -ne '2,${
+ s/^ *//
+ s/ \\*$//
+ s/$/:/
+ p
+ }' "$tmpdepfile" >> "$depfile"
+ else
+ make_dummy_depfile
+ fi
+ rm -f "$tmpdepfile" "$tmpdepfile2"
+ ;;
+
+tru64)
+ # The Tru64 compiler uses -MD to generate dependencies as a side
+ # effect. 'cc -MD -o foo.o ...' puts the dependencies into 'foo.o.d'.
+ # At least on Alpha/Redhat 6.1, Compaq CCC V6.2-504 seems to put
+ # dependencies in 'foo.d' instead, so we check for that too.
+ # Subdirectories are respected.
+ set_dir_from "$object"
+ set_base_from "$object"
+
+ if test "$libtool" = yes; then
+ # Libtool generates 2 separate objects for the 2 libraries. These
+ # two compilations output dependencies in $dir.libs/$base.o.d and
+ # in $dir$base.o.d. We have to check for both files, because
+ # one of the two compilations can be disabled. We should prefer
+ # $dir$base.o.d over $dir.libs/$base.o.d because the latter is
+ # automatically cleaned when .libs/ is deleted, while ignoring
+ # the former would cause a distcleancheck panic.
+ tmpdepfile1=$dir$base.o.d # libtool 1.5
+ tmpdepfile2=$dir.libs/$base.o.d # Likewise.
+ tmpdepfile3=$dir.libs/$base.d # Compaq CCC V6.2-504
+ "$@" -Wc,-MD
+ else
+ tmpdepfile1=$dir$base.d
+ tmpdepfile2=$dir$base.d
+ tmpdepfile3=$dir$base.d
+ "$@" -MD
+ fi
+
+ stat=$?
+ if test $stat -ne 0; then
+ rm -f "$tmpdepfile1" "$tmpdepfile2" "$tmpdepfile3"
+ exit $stat
+ fi
+
+ for tmpdepfile in "$tmpdepfile1" "$tmpdepfile2" "$tmpdepfile3"
+ do
+ test -f "$tmpdepfile" && break
+ done
+ # Same post-processing that is required for AIX mode.
+ aix_post_process_depfile
+ ;;
+
+msvc7)
+ if test "$libtool" = yes; then
+ showIncludes=-Wc,-showIncludes
+ else
+ showIncludes=-showIncludes
+ fi
+ "$@" $showIncludes > "$tmpdepfile"
+ stat=$?
+ grep -v '^Note: including file: ' "$tmpdepfile"
+ if test $stat -ne 0; then
+ rm -f "$tmpdepfile"
+ exit $stat
+ fi
+ rm -f "$depfile"
+ echo "$object : \\" > "$depfile"
+ # The first sed program below extracts the file names and escapes
+ # backslashes for cygpath. The second sed program outputs the file
+ # name when reading, but also accumulates all include files in the
+ # hold buffer in order to output them again at the end. This only
+ # works with sed implementations that can handle large buffers.
+ sed < "$tmpdepfile" -n '
+/^Note: including file: *\(.*\)/ {
+ s//\1/
+ s/\\/\\\\/g
+ p
+}' | $cygpath_u | sort -u | sed -n '
+s/ /\\ /g
+s/\(.*\)/'"$tab"'\1 \\/p
+s/.\(.*\) \\/\1:/
+H
+$ {
+ s/.*/'"$tab"'/
+ G
+ p
+}' >> "$depfile"
+ echo >> "$depfile" # make sure the fragment doesn't end with a backslash
+ rm -f "$tmpdepfile"
+ ;;
+
+msvc7msys)
+ # This case exists only to let depend.m4 do its work. It works by
+ # looking at the text of this script. This case will never be run,
+ # since it is checked for above.
+ exit 1
+ ;;
+
+#nosideeffect)
+ # This comment above is used by automake to tell side-effect
+ # dependency tracking mechanisms from slower ones.
+
+dashmstdout)
+ # Important note: in order to support this mode, a compiler *must*
+ # always write the preprocessed file to stdout, regardless of -o.
+ "$@" || exit $?
+
+ # Remove the call to Libtool.
+ if test "$libtool" = yes; then
+ while test "X$1" != 'X--mode=compile'; do
+ shift
+ done
+ shift
+ fi
+
+ # Remove '-o $object'.
+ IFS=" "
+ for arg
+ do
+ case $arg in
+ -o)
+ shift
+ ;;
+ $object)
+ shift
+ ;;
+ *)
+ set fnord "$@" "$arg"
+ shift # fnord
+ shift # $arg
+ ;;
+ esac
+ done
+
+ test -z "$dashmflag" && dashmflag=-M
+ # Require at least two characters before searching for ':'
+ # in the target name. This is to cope with DOS-style filenames:
+ # a dependency such as 'c:/foo/bar' could be seen as target 'c' otherwise.
+ "$@" $dashmflag |
+ sed "s|^[$tab ]*[^:$tab ][^:][^:]*:[$tab ]*|$object: |" > "$tmpdepfile"
+ rm -f "$depfile"
+ cat < "$tmpdepfile" > "$depfile"
+ # Some versions of the HPUX 10.20 sed can't process this sed invocation
+ # correctly. Breaking it into two sed invocations is a workaround.
+ tr ' ' "$nl" < "$tmpdepfile" \
+ | sed -e 's/^\\$//' -e '/^$/d' -e '/:$/d' \
+ | sed -e 's/$/ :/' >> "$depfile"
+ rm -f "$tmpdepfile"
+ ;;
+
+dashXmstdout)
+ # This case only exists to satisfy depend.m4. It is never actually
+ # run, as this mode is specially recognized in the preamble.
+ exit 1
+ ;;
+
+makedepend)
+ "$@" || exit $?
+ # Remove any Libtool call
+ if test "$libtool" = yes; then
+ while test "X$1" != 'X--mode=compile'; do
+ shift
+ done
+ shift
+ fi
+ # X makedepend
+ shift
+ cleared=no eat=no
+ for arg
+ do
+ case $cleared in
+ no)
+ set ""; shift
+ cleared=yes ;;
+ esac
+ if test $eat = yes; then
+ eat=no
+ continue
+ fi
+ case "$arg" in
+ -D*|-I*)
+ set fnord "$@" "$arg"; shift ;;
+ # Strip any option that makedepend may not understand. Remove
+ # the object too, otherwise makedepend will parse it as a source file.
+ -arch)
+ eat=yes ;;
+ -*|$object)
+ ;;
+ *)
+ set fnord "$@" "$arg"; shift ;;
+ esac
+ done
+ obj_suffix=`echo "$object" | sed 's/^.*\././'`
+ touch "$tmpdepfile"
+ ${MAKEDEPEND-makedepend} -o"$obj_suffix" -f"$tmpdepfile" "$@"
+ rm -f "$depfile"
+ # makedepend may prepend the VPATH from the source file name to the object.
+ # No need to regex-escape $object, excess matching of '.' is harmless.
+ sed "s|^.*\($object *:\)|\1|" "$tmpdepfile" > "$depfile"
+ # Some versions of the HPUX 10.20 sed can't process the last invocation
+ # correctly. Breaking it into two sed invocations is a workaround.
+ sed '1,2d' "$tmpdepfile" \
+ | tr ' ' "$nl" \
+ | sed -e 's/^\\$//' -e '/^$/d' -e '/:$/d' \
+ | sed -e 's/$/ :/' >> "$depfile"
+ rm -f "$tmpdepfile" "$tmpdepfile".bak
+ ;;
+
+cpp)
+ # Important note: in order to support this mode, a compiler *must*
+ # always write the preprocessed file to stdout.
+ "$@" || exit $?
+
+ # Remove the call to Libtool.
+ if test "$libtool" = yes; then
+ while test "X$1" != 'X--mode=compile'; do
+ shift
+ done
+ shift
+ fi
+
+ # Remove '-o $object'.
+ IFS=" "
+ for arg
+ do
+ case $arg in
+ -o)
+ shift
+ ;;
+ $object)
+ shift
+ ;;
+ *)
+ set fnord "$@" "$arg"
+ shift # fnord
+ shift # $arg
+ ;;
+ esac
+ done
+
+ "$@" -E \
+ | sed -n -e '/^# [0-9][0-9]* "\([^"]*\)".*/ s:: \1 \\:p' \
+ -e '/^#line [0-9][0-9]* "\([^"]*\)".*/ s:: \1 \\:p' \
+ | sed '$ s: \\$::' > "$tmpdepfile"
+ rm -f "$depfile"
+ echo "$object : \\" > "$depfile"
+ cat < "$tmpdepfile" >> "$depfile"
+ sed < "$tmpdepfile" '/^$/d;s/^ //;s/ \\$//;s/$/ :/' >> "$depfile"
+ rm -f "$tmpdepfile"
+ ;;
+
+msvisualcpp)
+ # Important note: in order to support this mode, a compiler *must*
+ # always write the preprocessed file to stdout.
+ "$@" || exit $?
+
+ # Remove the call to Libtool.
+ if test "$libtool" = yes; then
+ while test "X$1" != 'X--mode=compile'; do
+ shift
+ done
+ shift
+ fi
+
+ IFS=" "
+ for arg
+ do
+ case "$arg" in
+ -o)
+ shift
+ ;;
+ $object)
+ shift
+ ;;
+ "-Gm"|"/Gm"|"-Gi"|"/Gi"|"-ZI"|"/ZI")
+ set fnord "$@"
+ shift
+ shift
+ ;;
+ *)
+ set fnord "$@" "$arg"
+ shift
+ shift
+ ;;
+ esac
+ done
+ "$@" -E 2>/dev/null |
+ sed -n '/^#line [0-9][0-9]* "\([^"]*\)"/ s::\1:p' | $cygpath_u | sort -u > "$tmpdepfile"
+ rm -f "$depfile"
+ echo "$object : \\" > "$depfile"
+ sed < "$tmpdepfile" -n -e 's% %\\ %g' -e '/^\(.*\)$/ s::'"$tab"'\1 \\:p' >> "$depfile"
+ echo "$tab" >> "$depfile"
+ sed < "$tmpdepfile" -n -e 's% %\\ %g' -e '/^\(.*\)$/ s::\1\::p' >> "$depfile"
+ rm -f "$tmpdepfile"
+ ;;
+
+msvcmsys)
+ # This case exists only to let depend.m4 do its work. It works by
+ # looking at the text of this script. This case will never be run,
+ # since it is checked for above.
+ exit 1
+ ;;
+
+none)
+ exec "$@"
+ ;;
+
+*)
+ echo "Unknown depmode $depmode" 1>&2
+ exit 1
+ ;;
+esac
+
+exit 0
+
+# Local Variables:
+# mode: shell-script
+# sh-indentation: 2
+# eval: (add-hook 'write-file-hooks 'time-stamp)
+# time-stamp-start: "scriptversion="
+# time-stamp-format: "%:y-%02m-%02d.%02H"
+# time-stamp-time-zone: "UTC"
+# time-stamp-end: "; # UTC"
+# End:
Property changes on: applgrid/tags/applgrid-01-06-25/depcomp
___________________________________________________________________
Added: svn:executable
## -0,0 +1 ##
+*
\ No newline at end of property
Index: applgrid/tags/applgrid-01-06-25/src/appl_igrid.cxx
===================================================================
--- applgrid/tags/applgrid-01-06-25/src/appl_igrid.cxx (revision 0)
+++ applgrid/tags/applgrid-01-06-25/src/appl_igrid.cxx (revision 1980)
@@ -0,0 +1,2778 @@
+/**
+ ** @file appl_igrid.cxx
+ **
+ ** @brief grid class - all the functions needed to create and
+ ** fill the grid from an NLO calculation program.
+ **
+ ** @author mark sutton
+ ** @date 2007/10/16 17:01:39
+ **
+ ** @copyright (C) 2002-2019 mark sutton (sutt @ cern.ch)
+ **
+ ** $Id: appl_igrid.cxx, v1.00 2007/10/16 17:01:39 sutt $
+ **
+ **/
+
+#include
+#include
+#include
+#include
+
+#include "amconfig.h"
+
+#include "appl_igrid.h"
+#include "appl_grid/appl_grid.h"
+#include "appl_grid/appl_pdf.h"
+#include "appl_grid/appl_timer.h"
+#include "appl_grid/vector_stream.h"
+#include "appl_grid/stream_grid.h"
+
+#include "hoppet_init.h"
+#include "threadManager.h"
+
+#ifdef USEROOT
+#include "TFile.h"
+#include "TObject.h"
+#include "TObjString.h"
+#include "TH3D.h"
+#include "TVectorT.h"
+#include "TFileString.h"
+#endif
+
+#include "Splitting.h"
+
+
+#include "appl_grid/lumi_pdf.h"
+
+#include "appl_grid/stream_vector.h"
+
+// splitting function code
+
+
+// pdf reweighting
+// bool igrid::m_reweight = false;
+// bool igrid::m_symmetrise = false;
+
+// variable tranformation parameters
+double appl::igrid::transvar = 5;
+bool appl::igrid::threads_disabled = true;
+
+double appl::igrid::lambda = 0.0625;
+
+
+#if __cplusplus <= 199711L
+// #warning Not using C++11 compilation - using substitute functions
+std::string std::to_string( int i );
+std::string std::to_string( double f );
+#endif
+
+
+static int ithread = 0;
+
+std::string label( int i ) {
+ char lab[64];
+ std::sprintf( lab, "thread-%d", i );
+ return lab;
+}
+
+
+double (*fuserscale)( double ) = 0;
+
+
+// static bool dbg = true;
+
+appl::igrid::igrid() :
+ threadManager( label(ithread++) ),
+ mfy(0), mfx(0),
+ m_parent(0),
+ m_Ny1(0), m_y1min(0), m_y1max(0), m_deltay1(0),
+ m_Ny2(0), m_y2min(0), m_y2max(0), m_deltay2(0),
+ m_yorder(0),
+ m_Ntau(0), m_taumin(0), m_taumax(0), m_deltatau(0), m_tauorder(0),
+ m_Nproc(0),
+ m_transform(""),
+ m_qtransform(""),
+ m_transvar(transvar),
+ m_lambda(lambda),
+ m_reweight(false),
+ m_symmetrise(false),
+ m_optimised(false),
+ m_weight(0),
+ m_fg1(0), m_fg2(0),
+ m_fsplit1(0), m_fsplit2(0),
+ m_fsplit12(0), m_fsplit22(0),
+ m_alphas(0),
+ m_taufilledmin(-1), m_taufilledmax(-1),
+ m_y1filledmin(-1), m_y1filledmax(-1),
+ m_y2filledmin(-1), m_y2filledmax(-1),
+ m_partons(13)
+{
+ init_fmap();
+ // std::string qtransform = "h0";
+ // if ( m_lambda==0.25 ) qtransform = "h1";
+ set_transforms( m_transform, mfx, mfy );
+ set_transforms( m_qtransform, mfQ2, mftau );
+
+ // std::cout << "igrid() (default) Ntau=" << m_Ntau << "\t" << fQ2(m_taumin) << " - " << fQ2(m_taumax) << std::endl;
+
+}
+
+
+
+
+
+// standard constructor
+appl::igrid::igrid(int NQ2, double Q2min, double Q2max, int Q2order,
+ int Nx, double xmin, double xmax, int xorder,
+ std::string transform, std::string qtransform,
+ int Nproc, bool disflag ):
+ threadManager( label(ithread++) ),
+ mfy(0), mfx(0),
+ m_parent(0),
+ m_Ny1(Nx), m_Ny2( disflag ? 1 : Nx ), m_yorder(xorder),
+ m_Ntau(NQ2), m_tauorder(Q2order),
+ m_Nproc(Nproc),
+ m_transform(transform),
+ m_qtransform(qtransform),
+ m_transvar(transvar),
+ m_lambda(lambda),
+ m_reweight(false),
+ m_symmetrise(false),
+ m_optimised(false),
+ m_weight(0),
+ m_fg1(0), m_fg2(0),
+ m_fsplit1(0), m_fsplit2(0),
+ m_fsplit12(0), m_fsplit22(0),
+ m_alphas(0),
+ m_DISgrid(disflag),
+ m_taufilledmin(-1), m_taufilledmax(-1),
+ m_y1filledmin(-1), m_y1filledmax(-1),
+ m_y2filledmin(-1), m_y2filledmax(-1),
+ m_partons(13)
+{
+ // std::cout << "igrid::igrid() transform=" << m_transform << std::endl;
+ init_fmap();
+
+ // std::string qtransform = "h0";
+ // if ( m_lambda==0.25 ) qtransform = "h1";
+ set_transforms( m_transform, mfx, mfy );
+ set_transforms( m_qtransform, mfQ2, mftau );
+
+
+ double ymin1 = fy(xmax);
+ double ymax1 = fy(xmin);
+
+ m_y1min = ymin1 ;
+ m_y1max = ymax1 ;
+ m_y2min = ymin1 ;
+ m_y2max = ymax1 ;
+
+ if ( m_DISgrid ) m_y2min = m_y2max = 1;
+
+ if ( m_Ny1>1 ) m_deltay1 = (m_y1max-m_y1min)/(m_Ny1-1);
+ else m_deltay1 = 0;
+
+ if ( m_Ny2>1 ) m_deltay2 = (m_y2max-m_y2min)/(m_Ny2-1);
+ else m_deltay2 = 0;
+
+ double taumin=ftau(Q2min);
+ double taumax=ftau(Q2max);
+ m_taumin = taumin;
+ m_taumax = taumax;
+ if ( m_Ntau>1 ) m_deltatau = (taumax-taumin)/(m_Ntau-1);
+ else m_deltatau = 0;
+
+ if ( m_Ny1-1start_thread();
+}
+
+
+
+
+
+// copy constructor
+appl::igrid::igrid(const appl::igrid& g) :
+ threadManager( label(ithread++) ),
+ mfy(0), mfx(0),
+ m_parent(0),
+ m_Ny1(g.m_Ny1),
+ m_y1min(g.m_y1min), m_y1max(g.m_y1max), m_deltay1(g.m_deltay1),
+ m_Ny2(g.m_Ny2),
+ m_y2min(g.m_y2min), m_y2max(g.m_y2max), m_deltay2(g.m_deltay2),
+ m_yorder(g.m_yorder),
+ m_Ntau(g.m_Ntau),
+ m_taumin(g.m_taumin), m_taumax(g.m_taumax), m_deltatau(g.m_deltatau), m_tauorder(g.m_tauorder),
+ m_Nproc(g.m_Nproc),
+ m_transform(g.m_transform),
+ m_qtransform(g.m_qtransform),
+ m_transvar(g.m_transvar),
+ m_lambda(g.m_lambda),
+ m_reweight(g.m_reweight),
+ m_symmetrise(g.m_symmetrise),
+ m_optimised(g.m_optimised),
+ m_weight(0),
+ m_fg1(0), m_fg2(0),
+ m_fsplit1(0), m_fsplit2(0),
+ m_fsplit12(0), m_fsplit22(0),
+ m_alphas(0),
+ m_taufilledmin(g.m_taufilledmin), m_taufilledmax(g.m_taufilledmax),
+ m_y1filledmin(g.m_y1filledmin), m_y1filledmax(g.m_y1filledmax),
+ m_y2filledmin(g.m_y2filledmin), m_y2filledmax(g.m_y2filledmax),
+ m_partons(g.m_partons)
+{
+ init_fmap();
+
+ // std::string qtransform = "h0";
+ // if ( m_lambda==0.25 ) qtransform = "h1";
+ set_transforms( m_transform, mfx, mfy );
+ set_transforms( m_qtransform, mfQ2, mftau );
+
+
+ m_weight = new SparseMatrix3d*[m_Nproc];
+
+
+ for( int ip=0 ; ipstart_thread();
+}
+
+
+
+
+void _setlimits( int& _min, int& _max, const int _mint, const int _maxt ) {
+ if ( _mint<=_maxt ) {
+ if ( _min==-1 || _min>_mint ) _min = _mint;
+ if ( _max==-1 || _max<_maxt ) _max = _maxt;
+ }
+}
+
+
+
+
+#ifdef USEROOT
+// read from a file
+appl::igrid::igrid(TFile& f, const std::string& s) :
+ threadManager( label(ithread++) ),
+ mfy(0), mfx(0),
+ m_parent(0),
+ m_Ny1(0), m_y1min(0), m_y1max(0), m_deltay1(0),
+ m_Ny2(0), m_y2min(0), m_y2max(0), m_deltay2(0),
+ m_yorder(0),
+ m_Ntau(0), m_taumin(0), m_taumax(0), m_deltatau(0), m_tauorder(0),
+ m_Nproc(0),
+ m_transform(""),
+ m_qtransform(""),
+ m_transvar(transvar),
+ m_lambda(lambda),
+ m_reweight(false),
+ m_symmetrise(false),
+ m_optimised(false),
+ m_weight(0),
+ m_fg1(0), m_fg2(0),
+ m_fsplit1(0), m_fsplit2(0),
+ m_fsplit12(0), m_fsplit22(0),
+ m_alphas(0),
+ m_taufilledmin(-1), m_taufilledmax(-1),
+ m_y1filledmin(-1), m_y1filledmax(-1),
+ m_y2filledmin(-1), m_y2filledmax(-1),
+ m_partons(13)
+{
+ // std::cout << "igrid::igrid()" << std::endl;
+
+ // get the name of the transform pair
+ TFileString _tag = *(TFileString*)f.Get((s+"/Transform").c_str());
+ m_transform = _tag[0];
+
+ TFileString* _qtag = (TFileString*)f.Get((s+"/QTransform").c_str());
+
+ if ( _qtag != 0 ) m_qtransform = _qtag->at(0);
+ else m_qtransform = "h0";
+
+
+ init_fmap();
+
+ // std::string qtransform = "h0";
+ // TFileString* _qtag = (TFileString*)f.Get((s+"/TransformQ").c_str());
+ // if ( _qtag ) qtransform = (*_qtag)[0];
+
+ // if ( m_lambda==0.25 ) qtransform = "h1";
+ set_transforms( m_transform, mfx, mfy );
+ set_transforms( m_qtransform, mfQ2, mftau );
+
+
+ // delete _transform;
+
+ // retrieve setup parameters
+ TVectorT* setup=(TVectorT*)f.Get((s+"/Parameters").c_str());
+ // f.GetObject((s+"/Parameters").c_str(), setup);
+
+ // NB: round integer variables to nearest integer
+ // in case (unlikely) truncation error during
+ // conversion to double when they were stored
+ m_Ny1 = int((*setup)(0)+0.5);
+ m_y1min = (*setup)(1);
+ m_y1max = (*setup)(2);
+
+ m_Ny2 = int((*setup)(3)+0.5);
+ m_y2min = (*setup)(4);
+ m_y2max = (*setup)(5);
+
+ m_yorder = int((*setup)(6)+0.5);
+
+ m_Ntau = int((*setup)(7)+0.5);
+ m_taumin = (*setup)(8);
+ m_taumax = (*setup)(9);
+ m_tauorder = int((*setup)(10)+0.5);
+
+ m_transvar = (*setup)(11);
+
+ m_Nproc = int((*setup)(12)+0.5);
+
+ m_reweight = ((*setup)(13)!=0 ? true : false );
+ m_symmetrise = ((*setup)(14)!=0 ? true : false );
+ m_optimised = ((*setup)(15)!=0 ? true : false );
+ m_DISgrid = ( setup->GetNoElements()>16 && (*setup)[16]!=0 ) ;
+
+ m_lambda = (*setup)(17);
+
+ if ( m_lambda==0 ) m_lambda = lambda;
+
+ delete setup;
+
+ // std::cout << "igrid::igrid() read setup" << std::endl;
+
+ // create grids
+ m_deltay1 = (m_y1max-m_y1min)/(m_Ny1-1);
+ m_deltay2 = (m_y2max-m_y2min)/(m_Ny2-1);
+
+ m_deltatau = (m_taumax-m_taumin)/(m_Ntau-1);
+
+ // int rawsize=0;
+ // int trimsize=0;
+
+ m_weight = new SparseMatrix3d*[m_Nproc];
+
+ for( int ip=0 ; iptrim();
+
+ // delete storage histogram
+ delete htmp;
+
+ // DON'T trim on reading unless the user wants it!!
+ // he can trim himself if need be!!
+ // rawsize += m_weight[ip]->size();
+ // m_weight[ip]->trim(); // trim the grid and do some book keeping
+ // trimsize += m_weight[ip]->size();
+ }
+
+ /// now calculate the actual limits of this grid
+
+ // static double _ctime = 0;
+
+ setlimits();
+
+ if ( !threads_disabled ) this->start_thread();
+
+}
+#endif
+
+
+
+
+
+
+// read from a file
+appl::igrid::igrid( appl::file& file, const std::string& s) :
+ threadManager( label(ithread++) ),
+ mfy(0), mfx(0),
+ m_parent(0),
+ m_Ny1(0), m_y1min(0), m_y1max(0), m_deltay1(0),
+ m_Ny2(0), m_y2min(0), m_y2max(0), m_deltay2(0),
+ m_yorder(0),
+ m_Ntau(0), m_taumin(0), m_taumax(0), m_deltatau(0), m_tauorder(0),
+ m_Nproc(0),
+ m_transform(""),
+ m_qtransform(""),
+ m_transvar(transvar),
+ m_lambda(lambda),
+ m_reweight(false),
+ m_symmetrise(false),
+ m_optimised(false),
+ m_weight(0),
+ m_fg1(0), m_fg2(0),
+ m_fsplit1(0), m_fsplit2(0),
+ m_fsplit12(0), m_fsplit22(0),
+ m_alphas(0),
+ m_taufilledmin(-1), m_taufilledmax(-1),
+ m_y1filledmin(-1), m_y1filledmax(-1),
+ m_y2filledmin(-1), m_y2filledmax(-1),
+ m_partons(13)
+{
+ // get the name of the transform pair
+
+ stream_vector tag = file.Read >( s+"/Transform" );
+ // stream_vector tag; file.Read(tag);
+
+ m_transform = tag[0];
+
+ stream_vector qtag = file.Read >( s+"/QTransform" );
+ // stream_vector qtag; file.Read( qtag );
+
+ if ( qtag.size()>0 ) m_qtransform = qtag.at(0);
+ else m_qtransform = "h0";
+
+ init_fmap();
+
+ // if ( m_lambda==0.25 ) qtransform = "h1";
+ set_transforms( m_transform, mfx, mfy );
+ set_transforms( m_qtransform, mfQ2, mftau );
+
+
+ // delete _transform;
+
+ // retrieve setup parameters
+
+ stream_vector setup = file.Read >( s+"/Parameters" );
+ // stream_vector setup; file.Read( setup );
+
+ // NB: round integer variables to nearest integer
+ // in case (unlikely) truncation error during
+ // conversion to double when they were stored
+ m_Ny1 = int(setup[0]+0.5);
+ m_y1min = setup[1];
+ m_y1max = setup[2];
+
+ m_Ny2 = int(setup[3]+0.5);
+ m_y2min = setup[4];
+ m_y2max = setup[5];
+
+ m_yorder = int(setup[6]+0.5);
+
+ m_Ntau = int(setup[7]+0.5);
+ m_taumin = setup[8];
+ m_taumax = setup[9];
+ m_tauorder = int(setup[10]+0.5);
+
+ m_transvar = setup[11];
+
+ m_Nproc = int(setup[12]+0.5);
+
+ m_reweight = (setup[13]!=0 ? true : false );
+ m_symmetrise = (setup[14]!=0 ? true : false );
+ m_optimised = (setup[15]!=0 ? true : false );
+ m_DISgrid = ( setup.size()>16 && setup[16]!=0 ) ;
+
+ m_lambda = setup[17];
+
+ if ( m_lambda==0 ) m_lambda = lambda;
+
+ // std::cout << "igrid::igrid() read setup" << std::endl;
+
+ // create grids
+ m_deltay1 = (m_y1max-m_y1min)/(m_Ny1-1);
+ m_deltay2 = (m_y2max-m_y2min)/(m_Ny2-1);
+
+ m_deltatau = (m_taumax-m_taumin)/(m_Ntau-1);
+
+ // int rawsize=0;
+ // int trimsize=0;
+
+ m_weight = new SparseMatrix3d*[m_Nproc];
+
+ for( int ip=0 ; ip( _name );
+ // stream_grid htmp; file.Read( htmp );
+
+ // create grid
+ m_weight[ip]=new SparseMatrix3d(htmp);
+
+ // save some space
+ m_weight[ip]->trim();
+
+ // delete storage histogram
+
+ // DON'T trim on reading unless the user wants it!!
+ // he can trim himself if need be!!
+ // rawsize += m_weight[ip]->size();
+ // m_weight[ip]->trim(); // trim the grid and do some book keeping
+ // trimsize += m_weight[ip]->size();
+ }
+
+ /// now calculate the actual limits of this grid
+
+ // static double _ctime = 0;
+
+ setlimits();
+
+ if ( !threads_disabled ) this->start_thread();
+
+
+}
+
+
+
+
+
+
+
+void appl::igrid::setlimits() {
+ if ( m_weight ) {
+ // struct timeval _tstart = appl_timer_start();
+ for ( int i=0 ; iempty() ) continue;
+ _setlimits( m_taufilledmin, m_taufilledmax, _weight->xmin(), _weight->xmax() );
+ _setlimits( m_y1filledmin, m_y1filledmax, _weight->ymin(), _weight->ymax() );
+ _setlimits( m_y2filledmin, m_y2filledmax, _weight->zmin(), _weight->zmax() );
+ }
+ }
+ }
+}
+
+
+
+// constructor common internals
+void appl::igrid::construct()
+{
+ // Initialize histograms representing the weight grid
+ for( int ip=0 ; ipWrite();
+
+ // write the name of the transform pair
+ TFileString("Transform", m_transform).Write();
+ TFileString("QTransform", m_qtransform).Write();
+
+
+ TVectorT* setup=new TVectorT(20); // a few spare
+
+ (*setup)(0) = m_Ny1;
+ (*setup)(1) = m_y1min;
+ (*setup)(2) = m_y1max;
+
+ (*setup)(3) = m_Ny2;
+ (*setup)(4) = m_y2min;
+ (*setup)(5) = m_y2max;
+
+ (*setup)(6) = m_yorder;
+
+ (*setup)(7) = m_Ntau;
+ (*setup)(8) = m_taumin;
+ (*setup)(9) = m_taumax;
+ (*setup)(10) = m_tauorder;
+
+ (*setup)(11) = m_transvar;
+
+ (*setup)(12) = m_Nproc;
+
+ (*setup)(13) = ( m_reweight ? 1 : 0 );
+ (*setup)(14) = ( m_symmetrise ? 1 : 0 );
+ (*setup)(15) = ( m_optimised ? 1 : 0 );
+ (*setup)(16) = ( m_DISgrid ? 1 : 0 );
+
+ (*setup)(17) = m_lambda;
+
+ setup->Write("Parameters");
+
+ int igridsize = 0;
+ int igridtrimsize = 0;
+
+ for ( int ip=0 ; ipsize();
+
+ igridsize += m_weight[ip]->size();
+
+ // trim it so that it's quicker to copy into the TH3D
+ m_weight[ip]->trim();
+
+ igridtrimsize += m_weight[ip]->size();
+
+ TH3D* h=m_weight[ip]->getTH3D(hname);
+ h->SetDirectory(0);
+ h->Write();
+ delete h; // is this dengerous??? will root try to delete it later? I guess not if we SetDirectory(0)
+ }
+
+#if 0
+ // std::cout << _name << " trimmed" << std::endl;
+ std::cout << _name << "\tsize=" << igridsize << "\t-> " << igridtrimsize;
+ if ( igridsize ) std::cout << "\t( " << igridtrimsize*100/igridsize << "% )";
+ std::cout << std::endl;
+#endif
+
+ d.pop();
+}
+
+#endif
+
+
+
+
+
+// write to file
+void appl::igrid::write( appl::file& _file, const std::string& _name) {
+
+ // write the name of the transform pair
+
+ _file.Write( stream_vector( _name+"/Transform", std::vector(1,m_transform) ) );
+ _file.Write( stream_vector( _name+"/QTransform", std::vector(1,m_qtransform) ) );
+
+ std::vector setup(20); // a few spare
+
+ setup[0] = m_Ny1;
+ setup[1] = m_y1min;
+ setup[2] = m_y1max;
+
+ setup[3] = m_Ny2;
+ setup[4] = m_y2min;
+ setup[5] = m_y2max;
+
+ setup[6] = m_yorder;
+
+ setup[7] = m_Ntau;
+ setup[8] = m_taumin;
+ setup[9] = m_taumax;
+ setup[10] = m_tauorder;
+
+ setup[11] = m_transvar;
+
+ setup[12] = m_Nproc;
+
+ setup[13] = ( m_reweight ? 1 : 0 );
+ setup[14] = ( m_symmetrise ? 1 : 0 );
+ setup[15] = ( m_optimised ? 1 : 0 );
+ setup[16] = ( m_DISgrid ? 1 : 0 );
+
+ setup[17] = m_lambda;
+
+ _file.Write( stream_vector( _name+"/Parameters", setup ) );
+
+ int igridsize = 0;
+ int igridtrimsize = 0;
+
+ for ( int ip=0 ; ipsize();
+
+ igridsize += m_weight[ip]->size();
+
+ // trim it so that it's quicker to copy into the TH3D
+ m_weight[ip]->trim();
+
+ igridtrimsize += m_weight[ip]->size();
+
+ stream_grid* h = m_weight[ip]->get(hname);
+
+ // std::cout << *h << std::endl;
+
+
+ _file.Write( *h );
+
+ delete h;
+ }
+
+#if 0
+ // std::cout << _name << " trimmed" << std::endl;
+ std::cout << _name << "\tsize=" << igridsize << "\t-> " << igridtrimsize;
+ if ( igridsize ) std::cout << "\t( " << igridtrimsize*100/igridsize << "% )";
+ std::cout << std::endl;
+#endif
+
+ // d.pop();
+
+
+}
+
+
+
+
+
+void appl::igrid::fill_DIS(const double x1, const double Q2, const double* weight)
+{
+
+ // find preferred node for low end of interpolation range
+ int k1=fk1(x1);
+ int k3=fkappa(Q2);
+
+ double u_y1 = ( fy(x1)-gety1(k1) )/deltay1();
+ double u_tau = ( ftau(Q2)-gettau(k3))/deltatau();
+
+ double fillweight, fI_factor;
+
+ // cache the interpolation coefficients so only
+ // have to calculate each one once
+ // NB: a (very naughty) hardcoded maximum interpolation order
+ // of 16 for code simplicity.
+ double _fI1[16];
+ double _fI3[16];
+
+ for( int i1=0 ; i1<=m_yorder ; i1++ ) _fI1[i1] = fI(i1, m_yorder, u_y1);
+ for( int i3=0 ; i3<=m_tauorder ; i3++ ) _fI3[i3] = fI(i3, m_tauorder, u_tau);
+
+ double invwfun = 1/(weightfun(x1));
+
+ // double fI_tmp0 = 0;
+ // double fI_tmp1 = 0;
+
+ for( int i3=0 ; i3<=m_tauorder ; i3++ ) { // "interpolation" loop in Q2
+
+ // (*m_weight[ip])(k3+i3, k1+i1, k2+i2) += fillweight;
+
+ // if ( m_reweight ) fI_tmp0 = _fI3[i3]*invwfun;
+ // else fI_tmp0 = _fI3[i3];
+
+ for( int i1=0 ; i1<=m_yorder ; i1++ ) { // "interpolation" loop in x1
+
+ // fI_tmp1 = _fI1[i1] * fI_tmp0;
+
+ fI_factor = _fI1[i1] * _fI3[i3];
+
+ if ( m_reweight ) fI_factor *= invwfun;
+
+ for( int ip=0 ; ip ";
+
+ (*m_weight[ip])(k3+i3, k1+i1, 0) += fillweight;
+
+ // std::cout << (*m_weight[ip])(k3+i3, k1+i1, k2+i2) << ")";
+
+ // m_weight[ip]->print();
+ // m_weight[ip]->fill_fast(k3+i3, k1+i1, k2+i2) += fillweight/wfun;
+ }
+
+ }
+ }
+}
+
+
+
+
+
+void appl::igrid::fill(const double x1, const double x2, const double Q2, const double* weight)
+{
+ if ( isDISgrid() ) return fill_DIS( x1, Q2, weight );
+
+ // find preferred node for low end of interpolation range
+ int k1=fk1(x1);
+ int k2=fk2(x2);
+ int k3=fkappa(Q2);
+
+ double u_y1 = ( fy(x1)-gety1(k1) )/deltay1();
+ double u_y2 = ( fy(x2)-gety2(k2) )/deltay2();
+ double u_tau = ( ftau(Q2)-gettau(k3))/deltatau();
+
+ double fillweight, fI_factor;
+
+ // cache the interpolation coefficients so only
+ // have to calculate each one once
+ // NB: a (very naughty) hardcoded maximum interpolation order
+ // of 16 for code simplicity.
+ double _fI1[16];
+ double _fI2[16];
+ double _fI3[16];
+
+ for( int i1=0 ; i1<=m_yorder ; i1++ ) _fI1[i1] = fI(i1, m_yorder, u_y1);
+ for( int i2=0 ; i2<=m_yorder ; i2++ ) _fI2[i2] = fI(i2, m_yorder, u_y2);
+ for( int i3=0 ; i3<=m_tauorder ; i3++ ) _fI3[i3] = fI(i3, m_tauorder, u_tau);
+
+ double invwfun = 1/(weightfun(x1)*weightfun(x2));
+
+ // double fI_tmp0 = 0;
+ // double fI_tmp1 = 0;
+
+ for( int i3=0 ; i3<=m_tauorder ; i3++ ) { // "interpolation" loop in Q2
+
+ // (*m_weight[ip])(k3+i3, k1+i1, k2+i2) += fillweight;
+
+ // if ( m_reweight ) fI_tmp0 = _fI3[i3]*invwfun;
+ // else fI_tmp0 = _fI3[i3];
+
+ for( int i1=0 ; i1<=m_yorder ; i1++ ) { // "interpolation" loop in x1
+
+ // fI_tmp1 = _fI1[i1] * fI_tmp0;
+
+ for( int i2=0 ; i2<=m_yorder ; i2++ ) { // "interpolation" loop in x2
+
+ // fI_factor = fI_tmp1 * _fI2[i2];
+
+ fI_factor = _fI1[i1] * _fI2[i2] * _fI3[i3];
+
+ if ( m_reweight ) fI_factor *= invwfun;
+
+ for( int ip=0 ; ip ";
+
+ (*m_weight[ip])(k3+i3, k1+i1, k2+i2) += fillweight;
+
+ // std::cout << (*m_weight[ip])(k3+i3, k1+i1, k2+i2) << ")";
+
+ // m_weight[ip]->print();
+ // m_weight[ip]->fill_fast(k3+i3, k1+i1, k2+i2) += fillweight/wfun;
+ }
+
+ // std::cout << std::endl;
+ }
+ }
+ }
+}
+
+
+void appl::igrid::fill_phasespace(const double x1, const double x2, const double Q2, const double* weight) {
+
+ if ( isDISgrid() ) return fill_DIS_phasespace( x1, Q2, weight );
+
+ int k1=fk1(x1);
+ int k2=fk2(x2);
+ int k3=fkappa(Q2);
+
+ // std::cout << "appl::igrid::fill_phasespace() k: " << k1 << " " << k2 << " " << k3 << std::endl;
+ // std::cout << "\tweights: ";
+
+ // for ( int ip=0 ; ipsize() << std::endl;
+
+ // for ( int ip=0 ; ipsize() << std::endl;
+
+ // for ( int ip=0 ; ip( Ntau(), 0 );
+
+ m_fg1 = std::vector< std::vector< std::vector > >( Ntau(), std::vector >(Ny1(), std::vector(14,0) ) );
+ if ( !isSymmetric() && !isDISgrid() ) {
+ m_fg2 = std::vector< std::vector< std::vector > >( Ntau(), std::vector >(Ny2(), std::vector(14,0) ) );
+ }
+
+ if ( nloop>=1 && use_split ) {
+ m_fsplit1 = std::vector< std::vector< std::vector > >( Ntau(), std::vector >(Ny1(), std::vector(14,0) ) );
+ if ( !isSymmetric() && !isDISgrid() ) {
+ m_fsplit2 = std::vector< std::vector< std::vector > >( Ntau(), std::vector >(Ny2(), std::vector(14,0) ) );
+ }
+ if ( nloop==2 ) {
+ m_fsplit12 = std::vector< std::vector< std::vector > >( Ntau(), std::vector >(Ny1(), std::vector(14,0) ) );
+ if ( !isSymmetric() && !isDISgrid() ) {
+ m_fsplit22 = std::vector< std::vector< std::vector > >( Ntau(), std::vector >(Ny2(), std::vector(14,0) ) );
+ }
+ }
+ }
+
+
+ const int n_y1 = Ny1();
+ const int n_y2 = Ny2();
+
+
+ bool scale_beams = false;
+ if ( beam_scale0!=1 || beam_scale1!=1 ) scale_beams = true;
+
+ if ( initialise_hoppet ) hoppet_init::assign( pdf0->pdf() );
+
+ // set up pdf grid, splitting function grid
+ // and alpha_s grid
+ // for ( int itau=m_Ntau ; itau-- ; ) {
+ for ( int itau=0 ; itaum_taufilledmax ) continue;
+
+ double tau = gettau(itau);
+ double Q2 = fQ2(tau);
+ double fQ = std::sqrt(Q2);
+
+ double Q = fQ;
+
+#if 1
+ /// not iplemented yet ...
+ if ( fuserscale ) {
+ Q = fuserscale(fQ);
+ double scale_factor = Q/fQ;
+ logf = logf_save + 2*std::log(scale_factor);
+ }
+#endif
+
+ // alpha_s table
+ m_alphas[itau] = alphas(rscale_factor*Q)*invtwopi;
+
+#if 0
+
+ std::cout << itau << "\ttau " << tau
+ << "\tQ2 " << Q2 << "\tQ " << Q
+ << "\talphas " << m_alphas[itau] << std::endl;
+
+
+ int iymin1 = m_weight[0]->ymin();
+ int iymax1 = m_weight[0]->ymax();
+
+ if ( isSymmetric() ) {
+ int iymin2 = iymin1;
+ int iymax2 = iymax1;
+ }
+ else {
+ int iymin2 = m_weight[0]->zmin();
+ int iymax2 = m_weight[0]->zmax();
+ }
+#endif
+
+ // grid not filled for iyiymax so no need to
+ // consider outside this range
+
+ // y1 tables
+ for ( int iy=n_y1 ; iy-- ; ) {
+
+ if ( iym_y1filledmax ) continue;
+
+ double y = gety1(iy);
+ double x = fx(y);
+ double fun = 1;
+
+#if 0
+ static bool ff = true;
+
+ if ( ff ) { std::cout << "ff: " << m_reweight << " " << fun << std::endl; ff = false; }
+
+ if ( itau==0 ) {
+ std::cout << iy << " y: " << y << " x: " << x << std::endl;
+ }
+#endif
+
+ if ( scale_beams ) x *= beam_scale0;
+
+ if ( x>=1 ) {
+ for ( int ip=0 ; ip=1 && use_split ) {
+ for ( int ip=0 ; ipevaluate(x, fscale_factor*Q, m_fg1[itau][iy]);
+
+ // std::cout << "PDF1: itau " << itau << " Q: " << Q << " ( " << fscale_factor*Q << " ) iy1: " << iy << " x: " << x << " g-pdf: " << m_fg1[itau][iy][6] << " - " << (m_fg1[itau][iy][6]/x) << std::endl;
+
+ double invx = 1/x;
+ for ( int ip=0 ; ip=1 && use_split ) {
+
+ double PLO[14]; for ( int ip=14 ; ip-- ; ) PLO[ip] = 0;
+
+ splitting( x, fscale_factor*Q, PLO, 1 ); /// PLO x pdf
+
+ for ( int ip=m_partons ; ip-- ; ) m_fsplit1[itau][iy][ip] = PLO[ip]*invx;
+
+ if ( m_reweight ) for ( int ip=m_partons ; ip-- ; ) m_fsplit1[itau][iy][ip] *= fun;
+
+ if (nloop>1) {
+ // splitting( x, fscale_factor*Q, m_fsplit12[itau][iy], 1 ); placeholder for hoppet call
+
+ double P2LO[14]; for ( int ip=14 ; ip-- ; ) P2LO[ip] = 0;
+
+ double PNLO[14]; for ( int ip=14 ; ip-- ; ) PNLO[ip] = 0;
+
+ splitting( x, fscale_factor*Q, P2LO, 11 ); /// PLO x PLO x pdf
+ splitting( x, fscale_factor*Q, PNLO, 2 ); /// PNLO x pdf
+
+ // for ( int ip=0 ; ip<13 ; ip++ ) m_fsplit12[itau][iy][ip] = logf*(beta0*(0.5*logf - logr)*PLO[ip] + 0.5*logf*P2LO[ip] - PNLO[ip] )*invx ;
+ // for ( int ip=0 ; ip<13 ; ip++ ) m_fsplit12[itau][iy][ip] = logf*(beta0*(0.5*logf - logr)*PLO[ip] + 0.5*logf*P2LO[ip] - PNLO[ip] )*invx ;
+
+ for ( int ip=0 ; ippdf() );
+
+ if ( ( !isSymmetric() && !isDISgrid() ) || ( pdf1 && pdf1!=pdf0 ) ) {
+
+ for ( int itau=0 ; itaum_taufilledmax ) continue;
+
+ double tau = gettau(itau);
+ double Q2 = fQ2(tau);
+ double fQ = std::sqrt(Q2);
+
+ double Q = fQ;
+
+#if 1
+ /// not iplemented yet ...
+ if ( fuserscale ) {
+ Q = fuserscale(fQ);
+ double scale_factor = Q/fQ;
+ logf = logf_save + 2*std::log(scale_factor);
+ }
+#endif
+
+ /// alpha_s table has already been filled
+ /// m_alphas[itau] = alphas(rscale_factor*Q)*invtwopi;
+
+ // y2 tables
+ // for ( int iy=iymin2 ; iy<=iymax2 ; iy++ ) {
+ for ( int iy=n_y2 ; iy-- ; ) {
+
+ if ( iym_y2filledmax ) continue;
+
+ double y = gety2(iy);
+ double x = fx(y);
+ double fun = 1;
+
+ if ( scale_beams ) x *= beam_scale1;
+
+ if ( x>=1 ) {
+ for ( int ip=0 ; ip=1 && use_split ) {
+ for ( int ip=0 ; ipevaluate(x, fscale_factor*Q, m_fg2[itau][iy]);
+
+ // std::cout << "PDF2: itau " << itau << " Q: " << Q << " ( " << fscale_factor*Q << " ) iy1: " << iy << " x: " << x << " g-pdf: " << m_fg2[itau][iy][6] << " - " << (m_fg2[itau][iy][6]/x) << std::endl;
+
+ double invx = 1/x;
+ for ( int ip=0 ; ip=1 && use_split ) {
+
+ double PLO[14]; for ( int ip=14 ; ip-- ; ) PLO[ip] = 0;
+ splitting( x, fscale_factor*Q, PLO, 1 ); /// PLO x pdf
+ for ( int ip=m_partons ; ip-- ; ) m_fsplit2[itau][iy][ip] = PLO[ip]*invx;
+ if ( m_reweight ) for ( int ip=m_partons ; ip-- ; ) m_fsplit2[itau][iy][ip] *= fun;
+
+ if (nloop > 1) {
+ //splitting( x, fscale_factor*Q, m_fsplit22[itau][iy], 1 );
+
+ double P2LO[14]; for ( int ip=14 ; ip-- ; ) P2LO[ip] = 0;
+ double PNLO[14]; for ( int ip=14 ; ip-- ; ) PNLO[ip] = 0;
+
+ splitting( x, fscale_factor*Q, P2LO, 11 ); /// PLO x PLO x pdf
+ splitting( x, fscale_factor*Q, PNLO, 2 ); /// PNLO x pdf
+
+ // for ( int ip=0 ; ip<13 ; ip++ ) m_fsplit22[itau][iy][ip] = 0.5*logf*(beta0*(logf - 2*logr)*PLO[ip] + logf*P2LO[ip] - 2*PNLO[ip] )*invx ;
+
+ for ( int ip=0 ; ip(genpdf) << std::endl;
+
+#ifndef PDFTHREAD
+
+ // is the grid empty
+ int size=0;
+ for ( int ip=0 ; ipempty() << std::endl;
+ if ( m_weight[ip]->empty() ) continue;
+ if ( !m_weight[ip]->trimmed() ) m_weight[ip]->trim();
+ // size += m_weight[ip]->xmax() - m_weight[ip]->xmin() + 1;
+ size++;
+ }
+
+ // std::cout << "\t" << name() << " :: size: " << size << std::endl;
+
+ // grid is empty
+ if ( size==0 ) return 0;
+
+ int nloop = std::fabs(_nloop);
+
+ setuppdf( alphas, pdf0, pdf1, nloop, rscale_factor, fscale_factor, Escale0, Escale1 );
+
+#endif
+
+ if ( threads_disabled ) convolute_internal();
+ else process();
+
+ dsigma = m_conv_param.dsigma;
+ // dsigmaNLO = m_conv_param.dsigmaNLO;
+ // dsigmaNNLO = m_conv_param.dsigmaNNLO;
+
+ return dsigma;
+}
+
+
+
+
+
+
+
+void appl::igrid::convolute_internal() {
+
+ // struct timeval mytimer = appl_timer_start();
+
+ int lo_order = m_conv_param.lo_order;
+ int _nloop = m_conv_param._nloop;
+
+ int nloop = std::fabs(_nloop);
+
+ double rscale_factor = m_conv_param.rscale_factor;
+ double fscale_factor = m_conv_param.fscale_factor;
+
+ appl_pdf* genpdf = m_conv_param.genpdf;
+
+ // static const double twopi = 2*M_PI;
+ // static const int nc = 3;
+ // //TC const int nf = 6;
+ // static const int nf = 5;
+ // static double beta0=(11.*nc-2.*nf)/(6.*twopi);
+ // const bool debug=false;
+
+
+ double dsigma = 0.; //, xsigma = 0.;
+ double dsigmaNLO = 0.; //, xsigma = 0.;
+ double dsigmaNNLO = 0.; //, xsigma = 0.;
+
+ int size=0;
+ for ( int ip=0 ; ipempty() ) continue;
+ if ( !m_weight[ip]->trimmed() ) m_weight[ip]->trim();
+ // size += m_weight[ip]->xmax() - m_weight[ip]->xmin() + 1;
+ size++;
+ }
+
+ // grid is empty
+ if ( size==0 ) return;
+
+#ifdef PDFTHREAD
+
+ double (*alphas)(const double& ) = m_conv_param.alphas;
+
+ NodeCache* pdf0 = m_conv_param.pdf0;
+ NodeCache* pdf1 = m_conv_param.pdf1;
+
+ double Escale0 = m_conv_param.Escale0;
+ double Escale1 = m_conv_param.Escale1;
+
+ // double _alphasn = 1.;
+ // double alphanplus1 = 0.;
+ // do the convolution
+ // if (debug) std::cout<empty() ) continue;
+
+ if ( !m_weight[ip]->trimmed() ) {
+ /// std::cerr << "igrid::convolute() naughty, naughty!" << std::endl;
+ m_weight[ip]->trim();
+ }
+ size += m_weight[ip]->xmax() - m_weight[ip]->xmin() + 1;
+ }
+
+ // std::cout << "\t" << name() << " " << ip << " :: size " << size << std::endl;
+
+ // grid is empty
+ if ( size==0 ) return;
+
+
+ setuppdf( alphas, pdf0, pdf1, nloop, rscale_factor, fscale_factor, Escale0, Escale1 );
+
+
+#endif
+
+ double* sig = new double[m_Nproc]; // weights from grid
+
+ double* H = new double[m_Nproc]; // generalised pdf
+
+ double* HA = 0; // generalised splitting functions
+ double* HB = 0; // generalised splitting functions
+
+ double* HA2 = 0; // generalised splitting functions at NNLO
+ double* HB2 = 0; // generalised splitting functions at NNLO
+ double* HAB = 0; // generalised splitting functions at NNLO
+
+ if ( nloop>=1 && fscale_factor!=1 ) {
+ HA = new double[m_Nproc]; // generalised splitting functions
+ HB = new double[m_Nproc]; // generalised splitting functions
+ if (nloop > 1 ) {
+ HA2 = new double[m_Nproc]; // generalised splitting functions at NNLO
+ HB2 = new double[m_Nproc]; // generalised splitting functions at NNLO
+ HAB = new double[m_Nproc]; // generalised splitting functions at NNLO
+ }
+ }
+
+ // cross section for this igrid
+
+ //char name[]="appl_grid:igrid::convolute(): ";
+ // static const double twopi = 2*M_PI;
+ static const int nc = 3;
+ //TC const int nf = 6;
+ static const int nf = 5;
+ // static double beta0=(11.*nc - 2.*nf)/(6.*twopi);
+ // static double beta0=(11.*nc - 2.*nf)/(6.*twopi);
+ static double beta0_twopi=(11.*nc - 2.*nf)/6;
+
+ // static double beta1=(34.*nc*nc - 3.*(nc-1./nc)*nf - 10.*nc*nf)/(12.*twopi*twopi);
+ static double beta1_twopi_twopi=(34.*nc*nc - 3.*(nc-1./nc)*nf - 10.*nc*nf)/12.;
+
+
+ double alphas_tmp = 0.;
+
+ double _alphasn = 1.;
+ double alphanplus1 = 0.;
+ // double alphanplus2 = 0.;
+
+
+ m_conv_param.dsigma = 0;
+ m_conv_param.dsigmaNLO = 0;
+ m_conv_param.dsigmaNNLO = 0;
+
+
+ // loop over the grid
+ //
+
+ double logr2 = std::log(rscale_factor*rscale_factor);
+ double logf2 = std::log(fscale_factor*fscale_factor);
+
+ // std::vector ccks( m_Nproc, 0 );
+
+
+ double rlocal = rscale_factor;
+ double flocal = fscale_factor;
+
+ if ( genpdf->size() != m_Nproc ) {
+ throw exception( "igrid::convolute() parton luminosity mismatch: "+
+ std::to_string(genpdf->size())+" "+
+ std::to_string(m_Nproc) );
+ }
+
+ for ( int itau=0 ; itautrimmed(itau,iy1,iy2) ) continue;
+ bool nonzero = false;
+ // basic convolution order component for either the born level
+ // or the convolution of the nlo grid with the pdf
+
+ for ( int ip=0 ; ipevaluate( &m_fg1[itau][iy1][0], &m_fg2[itau][iy2][0], H );
+ else genpdf->evaluate( &m_fg1[itau][iy1][0], 0, H );
+
+ // std::cout << "genpdf: " << *genpdf << "\tH: " << H[0] << " " << H[1] << " " << H[2] << std::endl;
+
+ // std::cout << iy1 << " " << iy2 << " " << itau << " " << std::endl;
+
+ // for ( int ip=0 ; ip(genpdf) << std::endl;
+
+ // do the convolution
+ double xsigma=0.;
+
+ if ( m_parent && m_parent->subproc()!=-1 ) {
+ int ip=m_parent->subproc();
+ xsigma+= sig[ip]*H[ip];
+ }
+ else {
+ for ( int ip=0 ; ip= 0 ) dsigma += _alphasn*xsigma;
+
+ // std::cout << "\t\tdsigma: " << dsigma << " ( alphasn: " << _alphasn << " )" << std::endl;
+
+
+ // now do the convolution for the variation of factorisation and
+ // renormalisation scales, proportional to the leading order weights
+ if ( nloop>=1 ) {
+
+ // renormalisation scale dependent bit
+
+ if ( rlocal!=1 ) {
+ // nlo relative ln mu_R^2 term
+
+ double fac = alphanplus1*lo_order*logr2*xsigma;
+
+ dsigmaNLO += fac*beta0_twopi;
+
+ if (nloop>1) {
+ dsigmaNNLO += fac*alphas_tmp*( beta1_twopi_twopi + 0.5*beta0_twopi*beta0_twopi*(1+lo_order)*logr2 );
+ }
+
+ }
+
+
+ // factorisation scale dependent bit
+
+ // nlo relative ln mu_F^2 term
+ if ( flocal!=1 ) {
+
+ if ( !isDISgrid() ) {
+ genpdf->evaluate( &m_fg1[itau][iy1][0], &m_fsplit2[itau][iy2][0], HA);
+ genpdf->evaluate( &m_fsplit1[itau][iy1][0], &m_fg2[itau][iy2][0], HB);
+ }
+ else {
+ for ( int ip=0 ; ipevaluate( &m_fsplit1[itau][iy1][0], 0, HB);
+ }
+
+ xsigma=0.; /// setting this to 0 ??? maybe this is ok because we have already used it for the rscale part
+
+ if ( m_parent && m_parent->subproc()!=-1 ) {
+ int ip=m_parent->subproc();
+ xsigma += sig[ip]*(HA[ip]+HB[ip]);
+ }
+ else {
+ for ( int ip=0 ; ip1 ) {
+
+ /// this is most strange - NNLOJET needs logr2*(lo_order+1)
+ /// whereas an additional calculation suggests it should be
+ /// logr2*lo_order
+ /// hmmm
+ dsigmaNNLO += fac*alphas_tmp*beta0_twopi*( 0.5*logf2 - logr2*(lo_order+1) );
+
+ /// plus the double splitting function components ...
+ /// PDFA2 PDFB0 alphas2 +PDFA1 PDFB1 alphas2 + PDFA0 PDFB2 alphas2
+
+ if ( !isDISgrid() ) {
+ genpdf->evaluate( &m_fsplit12[itau][iy1][0], &m_fg2[itau][iy2][0], HA2 );
+ genpdf->evaluate( &m_fg1[itau][iy1][0], &m_fsplit22[itau][iy2][0], HB2 );
+ genpdf->evaluate( &m_fsplit1[itau][iy1][0], &m_fsplit2[itau][iy2][0], HAB );
+ }
+ else {
+ for ( int ip=0 ; ipevaluate( &m_fsplit12[itau][iy1][0], 0, HA2 );
+ }
+
+ double fac2 = 0;
+ double fac3 = 0;
+
+ if ( m_parent && m_parent->subproc()!=-1 ) {
+ int ip=m_parent->subproc();
+ fac2 += sig[ip]*( HA2[ip] + HB2[ip] );
+ fac3 += sig[ip]*( HAB[ip] );
+ }
+ else {
+ for ( int ip=0 ; iptrimmed() ) {
+ // std::cout << "igrid::convolute() naughty, naughty!" << std::endl;
+ m_weight[ip]->trim();
+ }
+ size += m_weight[ip]->xmax() - m_weight[ip]->xmin() + 1;
+ }
+
+ // grid is empty
+ if ( size==0 ) return 0;
+
+ setuppdf( alphas, pdf0, pdf1, nloop, rscale_factor, fscale_factor, Escale0, Escale1 );
+
+ if ( threads_disabled ) amc_convolute_internal();
+ else process();
+
+ return m_conv_param.dsigma;
+}
+
+
+
+
+
+
+
+void appl::igrid::amc_convolute_internal() {
+
+ int lo_order = m_conv_param.lo_order;
+
+ // int _nloop = m_conv_param._nloop;
+
+ // int nloop = std::fabs(_nloop);
+
+ // double rscale_factor = m_conv_param.rscale_factor;
+ // double fscale_factor = m_conv_param.fscale_factor;
+
+ appl_pdf* genpdf = m_conv_param.genpdf;
+
+ static const double eightpisquared = 8*M_PI*M_PI;
+
+ double alphas_tmp = 0.;
+ double dsigma = 0.; //, xsigma = 0.;
+ double _alphasn = 1.;
+
+ double* sig = new double[m_Nproc]; // weights from grid
+ double* H = new double[m_Nproc]; // generalised pdf
+ double* HA = 0; // generalised splitting functions
+ double* HB = 0; // generalised splitting functions
+ // if ( nloop==1 && fscale_factor!=1 ) {
+ // HA = new double[m_Nproc]; // generalised splitting functions
+ // HB = new double[m_Nproc]; // generalised splitting functions
+ // }
+
+ // cross section for this igrid
+
+ // loop over the grid
+ //
+
+ for ( int itau=0 ; itautrimmed(itau,iy1,iy2) ) continue;
+ bool nonzero = false;
+ // basic convolution order component for either the born level
+ // or the convolution of the nlo grid with the pdf
+ for ( int ip=0 ; ipevaluate( &m_fg1[itau][iy1][0], &m_fg2[itau][iy2][0], H );
+
+ // do the convolution
+
+ double xsigma=0.;
+ for ( int ip=0 ; ip& keep ) {
+
+ /// save the old grids
+ int Nproc = m_Nproc;
+ SparseMatrix3d** w = m_weight;
+
+ /// resize
+ m_Nproc = keep.size();
+ m_weight = new SparseMatrix3d*[m_Nproc];
+
+ /// copy across the grids we want to save
+ for ( unsigned ip=0 ; ip " << ip << std::endl;
+
+ m_weight[ip] = w[keep[ip]]; /// move across the processes we want to keep
+ w[keep[ip]] = 0; /// now flag as 0
+ }
+
+ for ( int ip=0 ; ip=m_Nproc ) return false;
+
+ SparseMatrix3d** mw0 = m_weight;
+ SparseMatrix3d** mw = new SparseMatrix3d*[m_Nproc-1];
+
+ int ir=0;
+ int ip0=0;
+ for ( int ip=0 ; ipsize();
+
+ // std::cout << "ymin=" << gety(0) << "\tymax=" << gety(m_Ny-1)
+ // << "\txmin=" << fx(gety(m_Ny-1)) << "\txmax=" << fx(gety(0)) << std::endl;
+
+ trim();
+
+ // std::cout << "\tsize(trimmed)=" << m_weight[0]->size() << std::endl;
+
+#if 1
+
+ // overall igrid optimisation limits
+
+ int _y1setmin = Ny1();
+ int _y1setmax = -1;
+
+ int _y2setmin = Ny2();
+ int _y2setmax = -1;
+
+ int _tausetmin = Ntau();
+ int _tausetmax = -1;
+
+ // double oldy1min = m_y1min;
+ // double oldy1max = m_y1max;
+
+ // double oldy2min = m_y2min;
+ // double oldy2max = m_y2max;
+
+ // go through all the subprocess to get the limits
+ // FIXME: this is actually redundant, at the moment, it assumes all the subprocesses
+ // occupy the same phase space, so only the first need be tested, and if
+ // they don;t all occupy the same phase space, then they should each be
+ // optimised seperately
+
+ // find limits for this grid
+ // initialise to original grid limits
+ int y1setmin = _y1setmin;
+ int y1setmax = _y1setmax;
+
+ int y2setmin = _y2setmin;
+ int y2setmax = _y2setmax;
+
+ int tausetmin = _tausetmin;
+ int tausetmax = _tausetmax;
+
+ bool firstdebug = true;
+ bool firstempty = true;
+
+ for ( int ip=0 ; ipxmax() << " - " << m_weight[ip]->xmin()
+ << "\tx1 " << m_weight[ip]->ymax() << " - " << m_weight[ip]->ymin()
+ << "\tx2 " << m_weight[ip]->zmax() << " - " << m_weight[ip]->zmin();
+
+ // is it empty?
+ if ( m_weight[ip]->xmax()-m_weight[ip]->xmin()+1 == 0 ) {
+ if ( firstdebug && firstempty ) std::cout << "\tempty" << std::endl;
+ firstempty = false;
+ firstdebug = false;
+ continue;
+ }
+
+
+ firstdebug = false;
+
+ // m_weight[ip]->print();
+
+ // m_weight[ip]->yaxis().print(fx);
+
+ // std::cout << "initial ip=" << ip
+ // << "\tm_y2min=" << m_y2min << "\tm_y2max=" << m_y2max
+ // << "\tm_y1min=" << m_y1min << "\tm_y1max=" << m_y1max
+ // << std::endl;
+
+ // m_weight[ip]->print();
+
+ // find actual limits
+
+ // y1 optimisation
+ int ymin1 = m_weight[ip]->ymin();
+ int ymax1 = m_weight[ip]->ymax();
+
+ if ( ymin1zmin();
+ int ymax2 = m_weight[ip]->zmax();
+
+ if ( ymin2xmin();
+ int taumax = m_weight[ip]->xmax();
+
+ if ( taumintausetmax ) tausetmax = taumax;
+
+ }
+
+ // if grid is empty, do "nothing" ie create the grid with the same
+ // limits as before but with the new required number of bins
+ if ( y1setmax==-1 || y2setmax==-1 || tausetmax==-1 ) {
+ m_Ny1 = Nx1;
+ m_Ny2 = Nx2;
+ m_Ntau = NQ2;
+ m_deltay1 = (m_y1max-m_y1min)/(m_Ny1-1);
+ m_deltay2 = (m_y2max-m_y2min)/(m_Ny2-1);
+ m_deltatau = (m_taumax-m_taumin)/(m_Ntau-1);
+ }
+ else {
+
+
+
+ // y1 optimisation
+ // double oldy1min = m_y1min;
+ // double oldy1max = m_y1max;
+
+
+ if ( isOptimised() ) {
+ // add a bit on each side
+ y1setmin -= extrabins;
+ y1setmax += extrabins;
+ }
+ else {
+ // not optimised yet so filled with phase space only so add
+ // the order to the max filled grid element position also
+ y1setmin -= extrabins;
+ y1setmax += m_yorder+extrabins;
+ }
+
+ if ( y1setmin<0 ) y1setmin = 0;
+ if ( y1setmax>m_Ny1-1 ) y1setmax = m_Ny1-1;
+
+
+ double _min = gety1(y1setmin);
+ double _max = gety1(y1setmax);
+
+ m_Ny1 = Nx1;
+ m_y1min = _min;
+ m_y1max = _max;
+ m_deltay1 = (m_y1max-m_y1min)/(m_Ny1-1);
+
+
+
+ // y2 optimisation
+ // double oldy2min = m_y2min;
+ // double oldy2max = m_y2max;
+
+ if ( isOptimised() ) {
+ // add a bit on each side
+ y2setmin -= extrabins;
+ y2setmax += extrabins;
+ }
+ else {
+ // not optimised yet so add the order to the
+ // max filled grid element position also
+ y2setmin -= extrabins;
+ y2setmax += m_yorder+extrabins;
+ }
+
+ if ( y2setmin<0 ) y2setmin = 0;
+ if ( y2setmax>=m_Ny2 ) y2setmax = m_Ny2-1;
+
+
+ // m_Ny2 = y2setmax-y2setmin+1;
+ _min = gety2(y2setmin);
+ _max = gety2(y2setmax);
+
+ m_Ny2 = Nx2;
+ m_y2min = _min;
+ m_y2max = _max;
+ m_deltay2 = (m_y2max-m_y2min)/(m_Ny2-1);
+
+ std::cout << "\t-> "
+ << y1setmin << " - " << y1setmax << " : "
+ << y2setmin << " - " << y2setmax
+ << std::endl;
+
+ // tau optimisation
+ // double oldtaumin = m_taumin;
+ // double oldtaumax = m_taumax;
+
+ // add a bit on each side
+ if ( isOptimised() ) {
+ // add a bit on each side
+ tausetmin--;
+ tausetmax++;
+ }
+ else {
+ // not optimised yet so add the order to the
+ // max filled grid element position also
+ tausetmax += m_tauorder+1;
+ tausetmin -= 1;
+ }
+
+ if ( tausetmin<0 ) tausetmin = 0;
+ if ( tausetmax>=m_Ntau-1 ) tausetmax = m_Ntau-1;
+
+
+ _min = gettau(tausetmin);
+ _max = gettau(tausetmax);
+
+ m_Ntau = NQ2;
+ m_taumin = _min;
+ m_taumax = _max;
+ m_deltatau = (m_taumax-m_taumin)/(m_Ntau-1);
+
+ // std::cout << "done ip=" << ip << std::endl;
+
+ }
+
+#endif
+
+ // now create the new subprocess grids with optimised limits
+ for ( int ip=0 ; ipempty() ) {
+ if ( m_weight[ip]->empty() ) {
+ delete m_weight[ip];
+ m_weight[ip] = new SparseMatrix3d( *g.m_weight[ip] );
+ }
+ else if ( m_weight[ip]->compare_axes( *g.m_weight[ip] ) ) (*m_weight[ip]) += (*g.m_weight[ip]);
+ else {
+ throw exception("igrid::operator+=() grids do not match");
+ }
+ }
+ }
+ }
+ return *this;
+}
+
+
+
+void appl::igrid::add( igrid* g ) {
+
+ int nx1 = g->Ny1();
+ int nx2 = g->Ny2();
+ int nQ2 = g->Ntau();
+
+ /// map nx, ny and ny to x1, x2 and Q2, then just
+ /// call fill for each !!!
+
+ igrid* ig[2] = { this, g };
+
+ std::cout << "appl::grid::add() combining bins: " << std::endl;
+
+ for ( int i=0 ; i<2 ; i++ ) std::cout << i << "\t" << *ig[i] << std::endl;
+
+ // for ( int i=0 ; im_weight[i]->untrim();
+
+ std::vector w(m_Nproc);
+
+ for ( int i=0 ; im_weight[i];
+
+ for ( int i1=0 ; i1gety1(i1);
+ double x1 = g->fx(y1);
+
+ for ( int i2=0 ; i2gety2(i2);
+ double x2 = g->fx(y2);
+
+ for ( int it=0 ; it weights(m_Nproc,0);
+
+ double nuffin = 0;
+
+ for ( int i=0 ; igettau(it);
+ double Q2 = g->fQ2(tau);
+
+ /// now fill this grid !!!
+ fill( x1, x2, Q2, &weights[0] );
+
+ }
+
+ } /// tau
+ } /// y2
+ } /// y1
+
+}
+
+#include
+
+void nodes( const std::string& s, double a, double b, int n ) {
+ double d = (b-a)/n;
+ std::printf( "%s\t", s.c_str() );
+ for ( int i=0 ; i sig(m_Nproc);
+
+ for ( int itau=0 ; itaufill( x1, x2, Q2, &sig[0] );
+ else g->fill_DIS( getx1(iy1), getQ2(itau), &sig[0] );
+ }
+ }
+ }
+ }
+
+ std::cout << "remap() out" << std::endl;
+
+}
+
+
+
+
+void appl::igrid::merge( igrid* g ) {
+
+
+ for ( int ip=0 ; ipNy1();
+ nx2[i] = g->Ny2();
+ nQ2[i] = g->Ntau();
+
+ y1[i] = g->gety1(0);
+ y1m[i] = g->gety1(nx1[i]-1);
+ dx1[i] = (y1m[i]-y1[i])/nx1[i];
+
+ y2[i] = g->gety2(0);
+ y2m[i] = g->gety2(nx2[i]-1);
+ dx2[i] = (y2m[i]-y2[i])/nx2[i];
+
+ tau[i] = g->gettau(0);
+ taum[i] = g->gettau(nQ2[i]-1);
+ dt[i] = (taum[i]-tau[i])/nQ2[i];
+
+ std::cout << "grid:" << i;
+ std::cout << "\t" << nx1[i] << "\tx1: " << "\t" << y1[i] << "\t " << y1m[i];
+ std::cout << "\t" << nx2[i] << "\tx2: " << "\t" << y2[i] << "\t " << y2m[i];
+ std::cout << "\t" << nQ2[i] << "\tQ2: " << "\t" << tau[i] << "\t " << taum[i];
+ std::cout << std::endl;
+
+ nodes( "y1: ", y1[i], y1m[i], nx1[i] );
+ // nodes( "y2: ", y2[i], y2m[i], nx2[i] );
+ // nodes( "tau: ", tau[i], taum[i], nQ2[i] );
+
+ }
+
+ bool rebuild1 = false;
+ bool rebuild2 = false;
+ bool rebuildt = false;
+
+ if ( y1[1]y1m[0] ) rebuild1 = true;
+
+ if ( y2[1]y2m[0] ) rebuild2 = true;
+
+ if ( tau[1]taum[0] ) rebuildt = true;
+
+ bool rebuild = rebuild1 || rebuild2 || rebuildt;
+
+ /// gt limits, now create the new igrid
+
+ std::cout << "rebuild: " << rebuild << std::endl;
+
+ m_yorder = 5;
+ m_tauorder = 5;
+
+ if ( rebuild ) {
+
+ m_y1min = std::min( y1[0], y1[1] );
+ m_y1max = std::max( y1m[0], y1m[1] );
+
+ double d1 = std::min( dx1[0], dx1[1] );
+
+ m_y2min = std::min( y2[0], y2[1] );
+ m_y2max = std::max( y2m[0], y2m[1] );
+
+ double d2 = std::min( dx2[0], dx2[1] );
+
+ m_taumin = std::min( tau[0], tau[1] );
+ m_taumax = std::max( taum[0], taum[1] );
+
+ double ndt = std::min( dt[0], dt[1] );
+
+ std::cout << "d: " << ndt << "\t" << d1 << "\t" << d2 << std::endl;
+ std::cout << "N: " << m_Ntau << "\t" << m_Ny1 << "\t" << m_Ny2 << std::endl;
+
+ if ( rebuildt ) m_Ntau = int((m_taumax-m_taumin)/ndt + 1);
+ if ( rebuild1 ) m_Ny1 = int((m_y1max-m_y1min)/d1 + 1);
+ if ( rebuild2 ) m_Ny2 = int((m_y2max-m_y2min)/d2 + 1);
+
+
+ nodes( "y1 ", m_y1min, m_y1max, m_Ny1 );
+ // nodes( "y2 ", m_y2min, m_y2max, m_Ny2 );
+ // nodes( "tau", m_taumin, m_taumax, m_Ntau );
+
+ std::cout << "N: " << m_Ntau << "\t" << m_Ny1 << "\t" << m_Ny2 << std::endl;
+
+ std::cout << "xorder: " << m_yorder << "\tQorder: " << m_tauorder << std::endl;
+
+ std::cout << "copy grid" << std::endl;
+
+ igrid* gsave = new igrid( *g0 );
+
+ std::cout << "new grids" << std::endl;
+
+ // now create the new subprocess grids with optimised limits
+ for ( int ip=0 ; ipm_weight[i]->untrim();
+
+ std::vector w(m_Nproc);
+
+ for ( int i=0 ; im_weight[i];
+
+ for ( int i1=0 ; i1gety1(i1);
+ double x1 = g->fx(y1);
+
+ for ( int i2=0 ; i2gety2(i2);
+ double x2 = g->fx(y2);
+
+ for ( int it=0 ; it weights(m_Nproc,0);
+
+ double nuffin = 0;
+
+ for ( int i=0 ; igettau(it);
+ double Q2 = g->fQ2(tau);
+
+ /// now fill this grid !!!
+ fill( x1, x2, Q2, &weights[0] );
+
+ }
+
+ } /// tau
+ } /// y2
+ } /// y1
+
+#endif
+
+}
+
+
+
+
+
+
+// debug print
+std::ostream& appl::igrid::debug(std::ostream& s) const {
+ header(std::cout);
+ for ( int i=0 ; ixaxis().min();
+ double xmax = m_weight[i]->xaxis().max();
+
+ double ymin = m_weight[i]->yaxis().min();
+ double ymax = m_weight[i]->yaxis().max();
+
+ double zmin = m_weight[i]->zaxis().min();
+ double zmax = m_weight[i]->zaxis().max();
+
+ std::cout << "ranges: Q2 : " << fQ2(xmin) << " " << fQ2(xmax) << " : " << xmin << " " << xmax << std::endl;
+ std::cout << "ranges: x1 : " << fx(ymin) << " " << fx(ymax) << " : " << ymin << " " << ymax << std::endl;
+ std::cout << "ranges: x2 : " << fx(zmin) << " " << fx(zmax) << " : " << zmin << " " << zmax << std::endl;
+
+ }
+ return s;
+}
+
+
+std::ostream& appl::igrid::header(std::ostream& s) const {
+ // s << "interpolation orders: x=" << g.yorder() << "\tQ2=" << g.tauorder() << std::endl;
+ s << "\t x: [ " << std::setw(2)
+ // << Ny1() << " ;\t" << std::setw(7) << fx(y1max()) << " - " << std::setw(7) << std::setprecision(6) << fx(y1min()) << "\t : "
+ // << Ny2() << " ;\t" << std::setw(7) << fx(y2max()) << " - " << std::setw(7) << std::setprecision(6) << fx(y2min())
+ << Ny1() << " :\t " << std::setw(7) << std::setprecision(6) << fx(y1max()) << " - " << std::setw(7) << std::setprecision(6) << fx(y1min()) << "\t : "
+ << Ny2() << " :\t " << std::setw(7) << std::setprecision(6) << fx(y2max()) << " - " << std::setw(7) << std::setprecision(6) << fx(y2min())
+ << "\t : " << "\t( order=" << yorder() << " ) ]";
+ s << "\t Q2: [ "
+ << Ntau() << " :\t " << std::setw(7) << std::setprecision(6) << fQ2(taumin()) << " - " << std::setw(7) << std::setprecision(6) << fQ2(taumax())
+ << "\t( order=" << tauorder() << " - reweight " << ( m_reweight ? "on " : "off" ) << ") ]";
+ return s;
+}
+
+
+std::ostream& operator<<(std::ostream& s, const appl::igrid& g) {
+ g.header(s);
+ // g.print(s);
+ return s;
+}
+
+
+
+
+void appl::igrid::run_thread() {
+
+ lock_proc();
+ mprocessing = false;
+ unlock_proc();
+
+
+ while( true ) {
+
+ // std::cout << "thread " << this << " " << mname << " waiting ... " << std::endl;
+
+ /// supend immediately
+ suspend();
+
+ if ( mterminate ) break;
+
+ // std::cout << "thread " << this << " " << mname << " running ... " << std::endl;
+
+ // struct timeval mytimer = appl_timer_start();
+
+ /// put the actual convoluting steps in here
+ /// once round the colbolution step, it puts itself to
+ /// sleep again
+
+ if ( parent()->calculation()==grid::AMCATNLO ) amc_convolute_internal();
+ else convolute_internal();
+
+ // std::printf("thread done: param: %lf internal\n", m_conv_param.dsigma );
+
+ /// after starting the processing step, the calling thread must
+ /// wait for all the threads to finish
+
+ // double mytime = appl_timer_stop(mytimer);
+
+ // std::printf("thread done: param: %lf internal time %lf ms\n", m_conv_param.dsigma, mytime );
+
+ /// signal computation has finished
+
+ // std::cout << "thread " << this << " " << mname << " done [" << mytime << " ms" << "]" << std::endl;
+ }
+
+ suspend(false);
+
+}
+
+
Index: applgrid/tags/applgrid-01-06-25/src/TFileVector.h
===================================================================
--- applgrid/tags/applgrid-01-06-25/src/TFileVector.h (revision 0)
+++ applgrid/tags/applgrid-01-06-25/src/TFileVector.h (revision 1980)
@@ -0,0 +1,82 @@
+/** emacs: this is -*- c++ -*- **/
+/**
+ ** @file TFileVector.h
+ **
+ ** @brief root TObject std::string std::vector class for writing std::string std::vectors
+ ** to root files
+ **
+ ** @author mark sutton
+ ** @date Sat Mar 15 19:49:16 GMT 2008
+ **
+ ** @copyright (C) 2002-2019 mark sutton (sutt @ cern.ch)
+ **
+ ** $Id: TFileVector.h, v0.0 Sat Mar 15 19:49:16 GMT 2008 sutt $
+ **
+ **/
+
+
+#ifndef TFILESTRING_H
+#define TFILESTRING_H
+
+#include "appl_grid/appl_root.h"
+
+#ifdef USEROOT
+
+#include
+#include
+#include
+
+#include "TObjString.h"
+#include "TObject.h"
+
+
+// template
+class TFileVector : public TObjString {
+
+public:
+
+ TFileVector(const std::string& name="") : TObjString(name.c_str()) { }
+
+ std::vector >& histos() { return mv; }
+ const std::vector >& histos() const { return mv; }
+
+ // get the name
+ std::string name() const { return GetName(); }
+
+ // get a value
+ std::vector& operator[](int i) { return mv[i]; }
+ const std::vector& operator[](int i) const { return mv[i]; }
+
+ // get the size
+ unsigned size() const { return mv.size(); }
+
+ // add an element
+ void add(const std::vector& s) {
+ mv.push_back(s);
+ }
+
+private:
+
+ std::vector > mv;
+
+ ClassDef(TFileVector, 1)
+
+};
+
+
+std::ostream& operator<<(std::ostream& s, const TFileVector& fv);
+
+
+#endif
+
+#endif // TFILEVECTOR_H
+
+
+
+
+
+
+
+
+
+
Index: applgrid/tags/applgrid-01-06-25/src/fappl_fitter.cxx
===================================================================
--- applgrid/tags/applgrid-01-06-25/src/fappl_fitter.cxx (revision 0)
+++ applgrid/tags/applgrid-01-06-25/src/fappl_fitter.cxx (revision 1980)
@@ -0,0 +1,152 @@
+/**
+ ** @file fappl_fitter.cxx
+ **
+ ** @author mark sutton
+ ** @date Fri 25 Oct 2013 00:49:13 CEST
+ **
+ ** @copyright (C) 2002-2019 mark sutton (sutt @ cern.ch)
+ **
+ ** $Id: fappl_fitter.cxx, v0.0 Fri 25 Oct 2013 00:49:13 CEST sutt $
+ **
+ **/
+
+
+
+#include "fappl_grid.h"
+
+
+// /// externally defined alpha_s and pdf routines for fortran
+// /// callable convolution wrapper
+// extern "C" double appl_fnalphas_(const double& Q);
+// extern "C" void appl_fnpdf_(const double& x, const double& Q, double* xf);
+
+// extern "C" double fnalphas_(const double& Q) { return appl_fnalphas_( Q); }
+// extern "C" void fnpdf_(const double& x, const double& Q, double* xf) { return appl_fnpdf_( x, Q, xf); }
+
+
+/// create a grid
+// extern "C" void appl_bookgrid_(int& id, const int& Nobs, const double* binlims) { bookgrid_( id, Nobs, binlims); }
+
+/// delete a grid
+extern "C" void appl_releasegrid_(int& id) { return releasegrid_( id ); }
+
+/// delete all grids
+extern "C" void appl_releasegrids_() { return releasegrids_(); }
+
+/// read a grid from a file
+extern "C" void appl_readgrid_(int& id, const char* s, int _len ) { return readgrid_( id, s, _len ); }
+
+/// write to a file
+extern "C" void appl_writegrid_(int& id, const char* s, int _len) { return writegrid_( id, s, _len ); }
+
+
+// /// add an entry
+// extern "C" void appl_fillgrid_(int& id,
+// const int& ix1, const int& ix2, const int& iQ,
+// const int& iobs,
+// const double* w,
+// const int& iorder ) {
+// return fillgrid_( id, ix1, ix2, iQ, iobs, w, iorder );
+// }
+
+// /// redefine the grid dimensions
+// extern "C" void appl_redefine_(int& id,
+// const int& iobs, const int& iorder,
+// const int& NQ2, const double& Q2min, const double& Q2max,
+// const int& Nx, const double& xmin, const double& xmax) {
+
+// return redefine_( id, iobs, iorder, NQ2, Q2min, Q2max, Nx, xmin, xmax );
+// }
+
+
+/// get the ckm matrix for a grid
+extern "C" void appl_getckm_( int& id, double* ckm ) { return getckm_( id, ckm ); }
+
+/// get the ckm matrix for a grid
+extern "C" void appl_setckm_( int& id, const double* ckm ) { return setckm_( id, ckm ); }
+
+/// get number of observable bins for a grid
+extern "C" int appl_getnbins_(int& id) { return getnbins_(id); }
+
+/// fins a bin number based on an ordinate
+extern "C" int appl_getbinnumber_(int& id, double& x) { return getbinnumber_(id,x); }
+
+/// get low edge for a grid/bin number
+extern "C" double appl_getbinlowedge_(int& id, int& bin) { return getbinlowedge_(id,bin); }
+
+/// get bin width for a grid/bin number
+extern "C" double appl_getbinwidth_(int& id, int& bin) { return getbinwidth_(id,bin); }
+
+
+/// do the convolution!! hooray!!
+
+extern "C" void appl_convolute_(int& id, double* data) { return convolute_( id, data ); }
+
+
+extern "C" void appl_convoluteorder_(int& id, int& nloops, double* data) { return convoluteorder_( id, nloops, data); }
+
+
+extern "C" void appl_convolutewrap_(int& id, double* data,
+ void (*pdf)(const double& , const double&, double* ),
+ double (*alphas)(const double& ) ) {
+ return convolutewrap_( id, data, pdf, alphas );
+}
+
+
+
+
+
+extern "C" void appl_fullconvolutewrap_(int& id, double* data,
+ void (*pdf)(const double& , const double&, double* ),
+ double (*alphas)(const double& ),
+ int& nloops,
+ double& rscale, double& fscale ) {
+ return fullconvolutewrap_( id, data, pdf, alphas, nloops, rscale, fscale );
+}
+
+
+
+extern "C" void appl_escaleconvolute_(int& id, double* data, double& Escale) {
+ return escaleconvolute_( id, data, Escale);
+}
+
+
+
+extern "C" void appl_escaleconvolutewrap_(int& id, double* data,
+ void (*pdf)(const double& , const double&, double* ),
+ double (*alphas)(const double& ), double& Escale ) {
+ return escaleconvolutewrap_( id, data, pdf, alphas, Escale );
+}
+
+
+
+extern "C" void appl_escalefullconvolutewrap_(int& id, double* data,
+ void (*pdf)(const double& , const double&, double* ),
+ double (*alphas)(const double& ),
+ int& nloops,
+ double& rscale, double& fscale,
+ double& Escale ) {
+ return escalefullconvolutewrap_( id, data, pdf, alphas, nloops, rscale, fscale, Escale );
+}
+
+
+
+/// print a grid
+extern "C" void appl_printgrid_(int& id) { return printgrid_( id ); }
+
+/// print all grids
+extern "C" void appl_printgrids_() { return printgrids_(); }
+
+/// print the grid documentation
+extern "C" void appl_printgriddoc_(int& id) { return printgriddoc_(id); }
+
+/// create grids from fastnlo
+extern "C" void appl_readfastnlogrids_( int* ids, const char* s, int _len ) { return readfastnlogrids_( ids, s, _len ); }
+
+
+/// grid std::map management
+
+extern "C" void appl_ngrids_(int& n) { return ngrids_( n ); }
+
+extern "C" void appl_gridids_(int* ids) { return gridids_( ids ); }
+
Index: applgrid/tags/applgrid-01-06-25/src/equals.h
===================================================================
--- applgrid/tags/applgrid-01-06-25/src/equals.h (revision 0)
+++ applgrid/tags/applgrid-01-06-25/src/equals.h (revision 1980)
@@ -0,0 +1,37 @@
+/* emacs: this is -*- c++ -*- */
+/**
+ ** @file equals.h
+ **
+ ** @author sutt
+ ** @date Wed 30 Dec 2020 11:26:54 GMT
+ **
+ ** $Id: equals.h, v0.0 Wed 30 Dec 2020 11:26:54 GMT sutt $
+ **
+ ** Copyright (C) 2020 sutt (sutt@cern.ch)
+ **
+ **/
+
+
+#ifndef EQUALS_H
+#define EQUALS_H
+
+
+inline bool equals( double a, double b, double limit=1 ) {
+ if ( a==b ) return true;
+ if ( std::fabs(a-b)<(std::fabs(std::min(a,b))*std::numeric_limits::epsilon()*limit ) ) return true;
+ if ( ( a==0 || b==0 ) && ( std::fabs(a-b)::epsilon()*limit ) ) return true;
+ return false;
+}
+
+
+#endif // EQUALS_H
+
+
+
+
+
+
+
+
+
+
Index: applgrid/tags/applgrid-01-06-25/src/convert.cxx
===================================================================
--- applgrid/tags/applgrid-01-06-25/src/convert.cxx (revision 0)
+++ applgrid/tags/applgrid-01-06-25/src/convert.cxx (revision 1980)
@@ -0,0 +1,66 @@
+
+
+#include "appl_grid/appl_grid.h"
+#include "amconfig.h"
+
+int Usage( std::ostream& s=std::cout, int i=0 ) {
+ s << "Usage: applgrid-convert [options] [grid1] [grid2] ...\n";
+ s << "\noption:\n";
+ s << " -p | --pdf value \t record specifed pdf to the grid\n";
+ s << " -i | --ipdf value \t record specifed ipdf from set to the grid\n\n";
+ s << " -h | --help \t this help\n";
+ s << "\nSee " << PACKAGE_URL << " for more details\n";
+ s << "\nReport bugs to <" << PACKAGE_BUGREPORT << ">";
+ s << std::endl;
+ return i;
+}
+
+
+int main( int argc, char** argv ) {
+
+ std::vector grids;
+
+ std::string pdf = "";
+ int ipdf = 0;
+
+ for ( int i=1 ; i " << newname << std::endl;
+ g.Write( newname );
+ }
+
+ return 0;
+}
Index: applgrid/tags/applgrid-01-06-25/src/nlojetpp_pdf.cxx
===================================================================
--- applgrid/tags/applgrid-01-06-25/src/nlojetpp_pdf.cxx (revision 0)
+++ applgrid/tags/applgrid-01-06-25/src/nlojetpp_pdf.cxx (revision 1980)
@@ -0,0 +1,29 @@
+/**
+ ** @file nlojetpp_pdf.cxx
+ **
+ ** @brief pdf transform function for nlojet
+ **
+ ** @author mark sutton
+ ** @date Mon Dec 10 01:36:04 GMT 2007
+ **
+ ** @copyright (C) 2002-2019 mark sutton (sutt @ cern.ch)
+ **
+ ** $Id: nlojetpp_pdf.cxx, v1.0 Mon Dec 10 01:36:04 GMT 2007 sutt $
+ **
+ **/
+
+
+
+#include "appl_grid/appl_pdf.h"
+
+#include "nlojetpp_pdf.h"
+
+
+extern "C" void fnlojetpp_pdf__(const double* fA, const double* fB, double* H) {
+ static nlojetpp_pdf pdf;
+ pdf.evaluate(fA, fB, H);
+}
+
+
+
+
Index: applgrid/tags/applgrid-01-06-25/src/mcfmw_pdf.h
===================================================================
--- applgrid/tags/applgrid-01-06-25/src/mcfmw_pdf.h (revision 0)
+++ applgrid/tags/applgrid-01-06-25/src/mcfmw_pdf.h (revision 1980)
@@ -0,0 +1,169 @@
+/** emacs: this is -*- c++ -*- **/
+/**
+ ** @file mcfmw_pdf.h
+ **
+ ** @brief pdf transform functions for W production
+ **
+ ** @author mark sutton
+ ** @date Mon Dec 10 01:36:04 GMT 2007
+ **
+ ** @copyright (C) 2002-2019 mark sutton (sutt @ cern.ch)
+ **
+ ** $Id: mcfmw_pdf.h, v1.0 Mon Dec 10 01:36:04 GMT 2007 sutt $
+ **
+ **/
+
+
+
+#ifndef __MCFMW_PDF_H
+#define __MCFMW_PDF_H
+
+
+#include "appl_grid/appl_pdf.h"
+
+//
+// MCFM W production
+//
+
+
+//
+// W+
+//
+class mcfmwp_pdf : public appl::appl_pdf {
+
+public:
+
+ mcfmwp_pdf(const std::string& s="mcfm-wp") : appl_pdf(s) {
+ m_Nproc = 6;
+ make_ckm( true );
+ }
+
+ ~mcfmwp_pdf() { }
+
+ virtual void evaluate(const double* fA, const double* fB, double* H) const;
+
+};
+
+
+
+//
+// W-
+//
+class mcfmwm_pdf : public appl::appl_pdf {
+
+public:
+
+ mcfmwm_pdf(const std::string& s="mcfm-wm") : appl_pdf(s) {
+ m_Nproc = 6;
+ make_ckm( false );
+ }
+
+ ~mcfmwm_pdf() { }
+
+ virtual void evaluate(const double* fA, const double* fB, double* H) const;
+
+};
+
+
+// actual funtion to evaluate the pdf combinations
+// for the W+
+
+inline void mcfmwp_pdf::evaluate(const double* fA, const double* fB, double* H) const {
+
+ const int nQuark = 6;
+ const int iQuark = 5;
+
+ // offset so we can use fA[-6]
+ // fA += 6;
+ // fB += 6;
+
+ double GA=fA[6];
+ double GB=fB[6];
+ double QA=0; double QB=0; double QbA=0; double QbB=0;
+
+ for(int i = 1; i <= iQuark; i++)
+ {
+ QA += fA[nQuark + i]*m_ckmsum[nQuark + i];
+ QB += fB[nQuark + i]*m_ckmsum[nQuark + i];
+ }
+ for(int i = -iQuark; i < 0; i++)
+ {
+ QbA += fA[nQuark + i]*m_ckmsum[nQuark + i];
+ QbB += fB[nQuark + i]*m_ckmsum[nQuark + i];
+ }
+
+ H[2]=QbA * GB;
+ H[3]= QA * GB;
+ H[4]= GA * QbB;
+ H[5]= GA * QB;
+
+ // have to zero H[0..1] first
+ H[0]=H[1]=0;
+
+ for (int i1 = 3; i1 <= 5; i1 += 2)
+ {
+ for(int i2 = 8; i2 <= 10; i2 += 2)
+ {
+ H[0] += fA[i1]*fB[i2]*m_ckm2[i1][i2];
+ H[1] += fA[i2]*fB[i1]*m_ckm2[i2][i1];
+ }
+ }
+}
+
+
+
+
+//
+// actual funtion to evaluate the pdf combinations
+// for the W-
+//
+inline void mcfmwm_pdf::evaluate(const double* fA, const double* fB, double* H) const {
+
+ const int nQuark = 6;
+ const int iQuark = 5;
+
+ // offset so we can use fA[-6]
+ // fA += 6;
+ // fB += 6;
+
+ double GA=fA[6];
+ double GB=fB[6];
+ double QA=0; double QB=0; double QbA=0; double QbB=0;
+
+ for(int i = 1; i <= iQuark; i++)
+ {
+ QA += fA[nQuark + i]*m_ckmsum[nQuark + i];
+ QB += fB[nQuark + i]*m_ckmsum[nQuark + i];
+ }
+ for(int i = -iQuark; i < 0; i++)
+ {
+ QbA += fA[nQuark + i]*m_ckmsum[nQuark + i];
+ QbB += fB[nQuark + i]*m_ckmsum[nQuark + i];
+ }
+
+ H[2]= QA * GB;
+ H[3]= QbA * GB;
+ H[4]= GA * QB;
+ H[5]= GA * QbB;
+
+ // have to zero H[0..1] first
+ H[0]=H[1]=0;
+
+ for (int i1 = 7; i1 <= 9; i1 += 2)
+ {
+ for(int i2 = 2; i2 <= 4; i2 += 2)
+ {
+ H[0] += fA[i1]*fB[i2]*m_ckm2[i1][i2];
+ H[1] += fA[i2]*fB[i1]*m_ckm2[i2][i1];
+ }
+ }
+}
+
+
+// fortran callable wrappers
+extern "C" void fmcfmwp_pdf__(const double* fA, const double* fB, double* H);
+extern "C" void fmcfmwm_pdf__(const double* fA, const double* fB, double* H);
+
+
+
+#endif // __MCFMW_PDF_H
Index: applgrid/tags/applgrid-01-06-25/src/jetrad_pdf.h
===================================================================
--- applgrid/tags/applgrid-01-06-25/src/jetrad_pdf.h (revision 0)
+++ applgrid/tags/applgrid-01-06-25/src/jetrad_pdf.h (revision 1980)
@@ -0,0 +1,85 @@
+/** emacs: this is -*- c++ -*- **/
+/**
+ ** @file jetrad_pdf.h
+ **
+ ** @brief pdf transform function for jetrad
+ **
+ ** @author mark sutton
+ ** @date Mon Dec 10 01:36:04 GMT 2007
+ **
+ ** @copyright (C) 2002-2019 mark sutton (sutt @ cern.ch)
+ **
+ ** $Id: jetrad_pdf.h, v1.0 Mon Dec 10 01:36:04 GMT 2007 sutt $
+ **
+ **/
+
+
+
+#ifndef __JETRAD_PDF_H
+#define __JETRAD_PDF_H
+
+
+#include "appl_grid/appl_pdf.h"
+
+
+
+// jetrad pdf combination class
+class jetrad_pdf : public appl::appl_pdf {
+
+public:
+
+ jetrad_pdf() : appl_pdf("jetrad") { m_Nproc=7; }
+
+ void evaluate(const double* _fA, const double* _fB, double* H) const;
+
+};
+
+
+
+inline void jetrad_pdf::evaluate(const double* _fA, const double* _fB, double* H) const {
+
+ // remapping from pdg -6..6 convention to u..t ubar..tbar g internal
+ // jetrad convention
+ static int index_map[13] = { 2, 1, 3, 4, 5, 6, -2, -1, -3, -4, -5, -6, 0 };
+
+ double f1[13];
+ double f2[13];
+
+ const double* _f1tmp = _fA+6;
+ const double* _f2tmp = _fB+6;
+
+ for ( int i=0 ; i<13 ; i++ ) f1[i] = _f1tmp[index_map[i]];
+ for ( int i=0 ; i<13 ; i++ ) f2[i] = _f2tmp[index_map[i]];
+
+
+ // now calculate the combinations
+
+ double Q1=0, Q1bar=0, Q2=0, Q2bar=0, D=0, Dbar=0;
+
+ for ( int i=0 ; i<6 ; i++ ) { Q1 += f1[i]; Q2 += f2[i]; }
+ for ( int i=6 ; i<12 ; i++ ) { Q1bar += f1[i]; Q2bar += f2[i]; }
+
+ for ( int i=0 ; i<12 ; i++ ) D += f1[i]*f2[i];
+
+ for ( int i=0 ; i<6 ; i++ ) {
+ Dbar += f1[i]*f2[i+6];
+ Dbar += f1[i+6]*f2[i];
+ }
+
+ H[0] = D;
+ H[1] = Q1*Q2+Q1bar*Q2bar-D ;
+ H[2] = Dbar ;
+ H[3] = Q1*Q2bar+Q1bar*Q2-Dbar ;
+ H[4] = (Q1+Q1bar)*f2[12] + f1[12]*(Q2+Q2bar) ;
+ H[5] = f1[12]*f2[12];
+
+}
+
+
+
+// fortran callable wrapper
+extern "C" void fjetrad_pdf__(const double* fA, const double* fB, double* H);
+
+
+#endif // __JETRAD_PDF_H
+
Index: applgrid/tags/applgrid-01-06-25/src/TFileString.cxx
===================================================================
--- applgrid/tags/applgrid-01-06-25/src/TFileString.cxx (revision 0)
+++ applgrid/tags/applgrid-01-06-25/src/TFileString.cxx (revision 1980)
@@ -0,0 +1,33 @@
+/**
+ ** @file TFileString.cxx
+ **
+ ** @brief root TObject std::string std::vector class for writing std::string std::vectors
+ ** to root files
+ **
+ ** @author mark sutton
+ ** @date Sat Mar 15 19:50:15 GMT 2008
+ **
+ ** @copyright (C) 2002-2019 mark sutton (sutt @ cern.ch)
+ **
+ ** $Id: TFileString.cxx, v0.0 Sat Mar 15 19:50:15 GMT 2008 sutt $
+ **
+ **/
+
+
+#include "appl_grid/appl_root.h"
+
+#ifdef USEROOT
+
+#include "TFileString.h"
+
+ClassImp(TFileString)
+
+
+std::ostream& operator<<(std::ostream& s, const TFileString& fs) {
+ s << fs.name() << ":";
+ for ( unsigned i=0 ; i
+
+#include "appl_grid/appl_pdf.h"
+
+
+//
+// MCFM Z production
+//
+class mcfmz_pdf : public appl::appl_pdf {
+
+public:
+
+ mcfmz_pdf() : appl_pdf("mcfm-z") { m_Nproc = 12; }
+
+ ~mcfmz_pdf() { }
+
+ virtual void evaluate(const double* fA, const double* fB, double* H) const;
+
+};
+
+
+// actual funtion to evaluate the pdf combinations
+
+inline void mcfmz_pdf::evaluate(const double* fA, const double* fB, double* H) const
+{
+
+ // const int nQuark = 6;
+ const int iQuark = 5;
+
+ // offset psd ptrs so can use [-6..6] indexing rather than [nQuark+..] indexing
+
+ fA += 6;
+ fB += 6;
+
+ // double GA=fA[6];
+ // double GB=fB[6];
+ double GA=fA[0];
+ double GB=fB[0];
+
+ double UpA=0; double UpB=0; double DnA=0; double DnB=0;
+ double UpbarA=0; double UpbarB=0; double DnbarA=0; double DnbarB=0;
+
+
+ for(int i = 1; i <= iQuark; i++)
+ {
+ if ((i % 2) == 0)
+ {
+ // UpA += fA[nQuark + i];
+ // UpB += fB[nQuark + i];
+ UpA += fA[i];
+ UpB += fB[i];
+ }
+ else
+ {
+ // DnA += fA[nQuark + i];
+ // DnB += fB[nQuark + i];
+ DnA += fA[i];
+ DnB += fB[i];
+ }
+ }
+
+ for(int i = -iQuark; i < 0; i++)
+ {
+ if (((int)(std::abs((double)i)) % 2) == 0)
+ {
+ // UpbarA += fA[nQuark + i];
+ // UpbarB += fB[nQuark + i];
+ UpbarA += fA[i];
+ UpbarB += fB[i];
+ }
+ else
+ {
+ // DnbarA += fA[nQuark + i];
+ // DnbarB += fB[nQuark + i];
+ DnbarA += fA[i];
+ DnbarB += fB[i];
+ }
+ }
+
+ // zero H first
+ H[0] = H[1] = H[2] = H[3] = 0;
+
+ static int _choice[13] = { 2, 3, 2, 3, 2, 3, 0, 1, 0, 1, 0, 1, 0 } ;
+ static int* choice = _choice+6;
+
+ for(int i = -iQuark; i <= iQuark; i++)
+ {
+ if (i == 0) continue;
+
+ // int Choice = ( i>0 ? i%2 : abs(i%2)+2 );
+ int Choice = choice[i];
+
+ H[Choice] += fA[i] * fB[-i];
+ }
+
+ H[4] = GA * UpB;
+ H[5] = GA * UpbarB;
+ H[6] = GA * DnB;
+ H[7] = GA * DnbarB;
+ H[8] = UpA *GB;
+ H[9] = UpbarA *GB;
+ H[10]= DnA *GB;
+ H[11]= DnbarA *GB;
+
+ return;
+}
+
+
+
+
+#endif // __MCFMZ_PDF_H
+
Index: applgrid/tags/applgrid-01-06-25/src/appl_file.cxx
===================================================================
--- applgrid/tags/applgrid-01-06-25/src/appl_file.cxx (revision 0)
+++ applgrid/tags/applgrid-01-06-25/src/appl_file.cxx (revision 1980)
@@ -0,0 +1,213 @@
+/**
+ ** @file appl_file.cxx
+ **
+ ** @author sutt
+ ** @date Wed 30 Dec 2020 15:25:42 GMT
+ **
+ ** $Id: appl_file.cxx, v0.0 Wed 30 Dec 2020 15:25:42 GMT sutt $
+ **
+ ** Copyright (C) 2020 sutt (sutt@cern.ch)
+ **
+ **/
+
+
+#include "appl_grid/appl_file.h"
+
+#include
+
+
+appl::file::file( const std::string& name, const std::string& option ) :
+ mfilename(name), mfile(0), mopen(false), msize(0) {
+
+ mopt = option;
+
+ /// make sure it is always binary mode
+ if ( mopt.find("b")==std::string::npos ) mopt += "b";
+
+ // std::cout << "appl::file::opening file: " << filename() << "\toptions: " << mopt << std::endl;
+
+ /// file sizes - actual filesize and uncopmpressed size
+ int filesize = 0;
+ int usize = 0;
+
+ if ( mopt.find("r")!=std::string::npos ) {
+
+ /// first get file size ...
+
+ struct stat stat_buf;
+ int rc = stat(filename().c_str(), &stat_buf);
+ if ( rc == 0 ) filesize = stat_buf.st_size;
+ else {
+ std::cerr << "appl::file: can not determine file size for file: " << filename() << std::endl;
+ return;
+ }
+
+ /// next check if it is in zipped format ...
+
+ FILE* tmp_file = fopen( filename().c_str(), "rb" );
+ int zip_signature = 0;
+ fread( &zip_signature, sizeof(int), 1, tmp_file );
+ if ( (zip_signature&0xffffff)==0x88b1f ) {
+ /// if it is read the file size, and the uncompressed filesize from the file ...
+ /// uncompressed size
+ int offset = filesize - 4;
+ fseek( tmp_file, offset, SEEK_SET );
+ fread( &usize, sizeof(int), 1, tmp_file );
+ }
+ else usize = filesize;
+
+ // std::cout << "\tfile size: " << filesize << std::endl;
+ // std::cout << "\tuncompressed size: " << usize << std::endl;
+
+ fclose( tmp_file );
+
+ }
+
+ /// now open the file ...
+
+ mfile = gzopen( mfilename.c_str(), mopt.c_str() );
+
+ mopen = true;
+
+ /// write header
+ if ( mopt.find("w")!=std::string::npos ) {
+
+ double v = SB::WRITEGUARD;
+ gzwrite( mfile, (void*)&v, sizeof(SB::TYPE) );
+ mindex.add( "header", sizeof(SB::TYPE) );
+
+ msize += sizeof(SB::TYPE);
+
+ /// the index is added at the end before the file is closed
+ }
+ if ( mopt.find("r")!=std::string::npos ) {
+
+ /// read the header ...
+
+ double v = 0;
+ gzread( mfile, (void*)&v, sizeof(SB::TYPE) );
+
+ /// if this one of our files ?
+ if ( v!=SB::WRITEGUARD ) {
+ std::cerr << "appl::file: incorrect file format file: " << filename() << std::endl;
+ Close();
+ mopen = false;
+ return;
+ }
+
+ /// yes it is so ...
+
+ /// first get the (uncompressed) file size from the last 4 bytes
+ /// of the compressed file, then use that get the offsets to the
+ /// index stored at the end of the file - *why* does this have to
+ /// be soooo difficult
+
+ // std::cout << "opened file: " << filename() << "\tsize: " << filesize << std::endl;
+
+ /// read the uncompressed file size from the comprfessed file
+ /// so that we can navigate through the commpressed file,
+ /// to get the index in case we want unsequential access
+
+ // std::cout << "uncompressed size: " << usize << std::endl;
+
+ /// get the total file size and offset of the index offset
+
+ SB::TYPE vtrailer[3];
+
+ gzseek( mfile, usize-sizeof(SB::TYPE)*3, SEEK_SET );
+ gzread( mfile, (void*)vtrailer, sizeof(double)*3 );
+
+ int index_size = (vtrailer[1]-vtrailer[0]-sizeof(SB::TYPE)*3); /// in bytes
+
+ std::vector vindex(index_size/sizeof(SB::TYPE));
+
+ int indexptr = vtrailer[0];
+
+ /// std::cout << "vtrailer[0]: " << vtrailer[0] << std::endl;
+ /// std::cout << "vtrailer[1]: " << vtrailer[1] << std::endl;
+ /// std::cout << "vtrailer[2]: " << vtrailer[2] << std::endl;
+
+ // gzrewind( mfile );
+ gzseek( mfile, indexptr, SEEK_SET );
+ gzread( mfile, (void*)&vindex[0], index_size );
+
+ mindex.deserialise( vindex );
+
+ // std::cout << mindex << std::endl;
+
+ /// rewind to open the file for proper reading ...
+ gzrewind( mfile );
+ /// read the header again, so that we are aligned with the
+ /// first object in the file
+ gzread( mfile, (void*)&v, sizeof(SB::TYPE) );
+
+ // std::cout << "write guard " << v << std::endl;
+
+ }
+}
+
+
+appl::file::~file() { Close(); }
+
+
+void appl::file::Close() {
+
+ if ( isopen() ) {
+
+ /// if writing the file, update the global size
+
+ if ( mopt.find("w")!=std::string::npos ) {
+ double content_trailer = SB::WRITEGUARD;
+ int bytes = gzwrite( mfile, (void*)&content_trailer, sizeof(SB::TYPE) );
+ mindex.add( "trailer", bytes );
+
+ msize += bytes;
+
+ double index_offset = msize;
+
+ // std::cout << "writing: index offset: " << index_offset << std::endl;
+
+ /// this is a thorny one - should we add the index to the index ?
+ /// until the index has been written, it should not be in the index
+ /// but once it has been written, then it *should* be in the index
+ /// so we would need to add the index to itself, *before* writing
+ /// to the file
+ /// basically, I think it is more sensible *not* to include the
+ /// index itself into the index, but of course, that could could
+ /// be changed if necessary
+ /// what we *should* do, is to add some pointer to the index at
+ /// the start of the file, so that we can navigate to it directly
+ /// on openning the file
+ Write( mindex );
+
+ /// this doesn;t work - need to work out why
+ // gzrewind( mfile );
+ // gzseek( mfile, sizeof(double), SEEK_SET );
+ // double tsize = msize;
+
+ SB::TYPE trailer[3] = { 0, 0, 0 };
+ trailer[0] = index_offset;
+ trailer[1] = msize+3*sizeof(SB::TYPE);
+ trailer[2] = SB::WRITEGUARD;
+
+ bytes = gzwrite( mfile, (void*)&trailer, sizeof(SB::TYPE)*3 );
+ // std::cout << "writing trailer: " << std::endl;
+ // std::cout << "\tindex offset: " << trailer[0] << std::endl;
+ // std::cout << "\ttotal size: " << trailer[1] << std::endl;
+ // std::cout << "\twrite guard: " << trailer[2] << std::endl;
+
+ mindex.add("file trailer", bytes );
+ }
+
+ gzclose(mfile);
+ mopen = false;
+
+ mindex.clear();
+
+ }
+
+}
+
+
+
+
Index: applgrid/tags/applgrid-01-06-25/src/tsparse1d.h
===================================================================
--- applgrid/tags/applgrid-01-06-25/src/tsparse1d.h (revision 0)
+++ applgrid/tags/applgrid-01-06-25/src/tsparse1d.h (revision 1980)
@@ -0,0 +1,197 @@
+/** emacs: this is -*- c++ -*- **/
+/**
+ ** @file tsparse1d.h
+ **
+ ** @brief very basic 2d sparse matrix class, needs to be improved,
+ ** and extended to 3d
+ **
+ ** @author mark sutton
+ ** @date Fri Nov 9 00:17:31 CET 2007
+ **
+ ** @copyright (C) 2002-2019 mark sutton (sutt @ cern.ch)
+ **
+ ** $Id$
+ **
+ **/
+
+
+#ifndef __TSPARSE1D_H
+#define __TSPARSE1D_H
+
+#include
+
+#include
+
+#include
+
+#include
+
+#include "tsparse_base.h"
+
+
+template
+class tsparse1d : public tsparse_base {
+
+public:
+
+ tsparse1d(int nx) : tsparse_base(nx), m_v(NULL) {
+ m_v = new T[m_Nx];
+ for ( int i=0 ; im_ux ) return 0;
+ return m_v[i-m_lx];
+ }
+
+
+#if 1
+ void print() const {
+ std::cout << " \tsize=" << size() << "\t";
+ for ( int i=0 ; i=m_lx && i<=m_ux ) return;
+
+ // is it empty?
+ if ( m_lx>m_ux ) {
+ m_v = new T[1];
+ (*m_v) = 0;
+ m_lx = m_ux = i;
+ return;
+ }
+
+ T* newv = new T[ ( ii ; m_lx-- ) (*newp++) = 0;
+ for ( ; lx<=ux ; lx++ ) (*newp++) = (*mvp++);
+ for ( ; m_uxstd::ostream& operator<<(std::ostream& s, const tsparse1d& sp) {
+ for ( int i=0 ; i
+#include
+
+#include "appl_grid/appl_timer.h"
+
+
+pthread_mutex_t time_lock = PTHREAD_MUTEX_INITIALIZER;
+
+
+typedef struct {
+ struct timeval start_time;
+ struct timeval stop_time;
+ struct timeval since_time;
+ struct timeval diff_time;
+} __appl_timer;
+
+
+
+int appl_timersub(struct timeval* stop_time, struct timeval* start_time, struct timeval* diff_time)
+{
+ diff_time->tv_sec = stop_time->tv_sec - start_time->tv_sec;
+ diff_time->tv_usec = stop_time->tv_usec - start_time->tv_usec;
+ return 0;
+}
+
+
+
+/** get the time, (probably needlessly) mutex protected
+ ** call to gettimeofday
+ **/
+
+void __appl__gettime(struct timeval* t)
+{
+ pthread_mutex_lock(&time_lock);
+ gettimeofday (t, NULL);
+ pthread_mutex_unlock(&time_lock);
+}
+
+
+
+/** significantly (0.02ms) faster and simpler timer start/stop
+ ** functions
+ **/
+
+struct timeval appl_timer_start(void)
+{
+ struct timeval start_time;
+ __appl__gettime (&start_time);
+ return start_time;
+}
+
+
+// double appl_d(struct timeval time) { return (time.tv_sec*1000.0) + (time.tv_usec/1000.0); }
+
+double appl_timer_stop(struct timeval start_time)
+{
+ double time = 0;
+ struct timeval stop_time;
+ struct timeval diff_time;
+ // double diff;
+
+ __appl__gettime (&stop_time);
+ pthread_mutex_lock(&time_lock);
+ appl_timersub( &stop_time, &start_time, &diff_time );
+ // printf("diff: %d\n", _timersub( &stop_time, &start_time, &diff_time ) );
+ // diff = appl_d(stop_time)-appl_d(start_time);
+ // printf("timer: %12.5f %12.5f %12.5f %12.5f\n", appl_d(stop_time), appl_d(start_time), appl_d(diff_time), diff );
+ pthread_mutex_unlock(&time_lock);
+ time = (diff_time.tv_sec*1000.0) + (diff_time.tv_usec/1000.0);
+ return time;
+}
+
+
+
+
+
+
+
+
+
+
+
Index: applgrid/tags/applgrid-01-06-25/src/tsparse3d.h
===================================================================
--- applgrid/tags/applgrid-01-06-25/src/tsparse3d.h (revision 0)
+++ applgrid/tags/applgrid-01-06-25/src/tsparse3d.h (revision 1980)
@@ -0,0 +1,376 @@
+/** emacs: this is -*- c++ -*- **/
+/**
+ ** @file tsparse3d.h
+ **
+ ** @brief very basic 3d sparse matrix class, needs to be improved
+ **
+ **
+ ** @author mark sutton
+ ** @date Fri Nov 9 00:17:31 CET 2007
+ **
+ ** @copyright (C) 2002-2019 mark sutton (sutt @ cern.ch)
+ **
+ ** $Id: tsparse3d.h, v Fri Nov 9 00:17:31 CET 2007 sutt $
+ **
+ **/
+
+
+#ifndef __TSPARSE3D_H
+#define __TSPARSE3D_H
+
+#include
+
+#include "tsparse2d.h"
+
+
+template
+class tsparse3d : public tsparse_base {
+
+public:
+
+ tsparse3d(int nx, int ny, int nz)
+ : tsparse_base(nx), m_Ny(ny), m_Nz(nz), m_v(NULL), m_trimmed(false) {
+ m_v = new tsparse2d*[m_Nx];
+ for ( int i=0 ; i(m_Ny, m_Nz);
+ // setup_fast();
+ }
+
+ // Fixme: need to rewrite this constructor properly, so that
+ // a trimmed matrix is copied as trimmed and not copied in full
+ // and then trimmed
+ tsparse3d(const tsparse3d& t)
+ : tsparse_base(t.m_Nx), m_Ny(t.m_Ny), m_Nz(t.m_Nz), m_v(NULL), m_trimmed(t.m_trimmed) {
+ m_v = new tsparse2d*[m_Nx];
+ for ( int i=0 ; i(m_Ny, m_Nz);
+ // setup_fast();
+ // deep copy of all elements
+ for ( int i=0 ; itrim();
+
+ int xmin = m_lx;
+ int xmax = m_ux;
+
+ // now find new limits
+ for ( ; xmin** m_vnew = new tsparse2d*[xmax-xmin+1];
+ tsparse2d** mnp = m_vnew;
+ tsparse2d** mvp = m_v+xmin-m_lx;
+ for ( int j=xmin ; j<=xmax ; j++ ) (*mnp++) = (*mvp++);
+ delete[] m_v;
+
+ // update the sparse data and corresponding limits
+ m_lx = xmin;
+ m_ux = xmax;
+ m_v = m_vnew;
+
+ m_empty = false;
+ if ( empty() ) m_empty = true;
+ }
+
+
+ void untrim() {
+ m_trimmed = false;
+ grow(0);
+ grow(Nx()-1);
+ tsparse2d** mvp = m_v;
+ for ( int i=m_lx ; i<=m_ux ; i++ ) (*mvp++)->untrim();
+ }
+
+ bool trimmed() const { return m_trimmed; }
+
+ bool trimmed(int i, int j, int k) const {
+ if ( im_ux ) return false;
+ else return m_v[i-m_lx]->trimmed(j,k);
+ }
+
+
+ T operator()(int i, int j, int k) const {
+ // range_check(i);
+ if ( im_ux ) return 0;
+ return (*((const tsparse2d** const)m_v)[i-m_lx])(j,k);
+ }
+
+
+ const tsparse2d* operator[](int i) const {
+ // range_check(i);
+ if ( im_ux ) return NULL;
+ return m_v[i-m_lx];
+ }
+
+
+ T& operator()(int i, int j, int k) {
+ // range_check(i);
+ grow(i);
+ return (*m_v[i-m_lx])(j,k);
+ }
+
+ int size() const {
+ int N=0;
+ for ( int i=m_ux-m_lx+1 ; i-- ; ) N += m_v[i]->size();
+ return N; // +3*sizeof(int);
+ }
+
+#if 0
+ void setup_fast() {
+ // set up fast lookup table into the (untrimmed) 3d array.
+ m_fastindex = new T*[Nx()*Ny()*Nz()];
+
+ for ( int i=0 ; iv()[j])->v()[k];
+ }
+ }
+ }
+ }
+
+
+ void empty_fast() {
+ if ( m_fastindex ) delete[] m_fastindex;
+ m_fastindex=NULL;
+ }
+
+ T& fill_fast(int i, int j, int k) {
+ return *m_fastindex[(i*Ny()+j)*Nz()+k];
+ }
+
+ T fill_fast(int i, int j, int k) const {
+ return *m_fastindex[(i*Ny()+j)*Nz()+k];
+ }
+#endif
+
+
+ void print() const {
+ if ( m_ux-m_lx+1==0 ) std::cout << "-" << "\n";
+ else {
+ for ( int i=0 ; iprint();
+ }
+ else {
+ std::cout << "- \t";
+ }
+ std::cout << "\n";
+ }
+ }
+ }
+
+ int ymin() const {
+ int minx = m_Ny;
+ tsparse2d** mvp = m_v;
+ for ( int i=m_lx ; i<=m_ux ; i++ ) {
+ int _minx=(*mvp++)->xmin();
+ if ( _minx** mvp = m_v;
+ for ( int i=m_lx ; i<=m_ux ; i++ ) {
+ int _minx=(*mvp)->xmin();
+ int _maxx=(*mvp++)->xmax();
+ if ( _minx<=_maxx && _maxx>maxx ) maxx=_maxx;
+ }
+ if ( maxx==-1 ) maxx=m_Ny-1;
+ return maxx;
+ }
+
+
+ int zmin() const {
+ int minx = m_Nz;
+ tsparse2d** mvp = m_v;
+ for ( int i=m_lx ; i<=m_ux ; i++ ) {
+ int _minx=(*mvp++)->ymin();
+ if ( _minx** mvp = m_v;
+ for ( int i=m_lx ; i<=m_ux ; i++ ) {
+ int _minx=(*mvp)->ymin();
+ int _maxx=(*mvp++)->ymax();
+ if ( _minx<=_maxx && _maxx>maxx ) maxx=_maxx;
+ }
+ if ( maxx==-1 ) maxx=m_Nz-1;
+ return maxx;
+ }
+
+
+ // algebraic operators
+
+ tsparse3d& operator=(const tsparse3d& t) {
+ // clearout what is already in the matrix
+ if ( m_v ) {
+ for ( int i=m_ux-m_lx+1 ; i-- ; ) if ( m_v[i] ) delete m_v[i];
+ delete[] m_v;
+ m_v = NULL;
+ }
+
+ m_trimmed = t.m_trimmed;
+
+ // now copy everything else
+ m_Nx = t.m_Nx;
+ m_Ny = t.m_Ny;
+ m_Nz = t.m_Nz;
+
+ m_lx = 0;
+ m_ux = m_Nx-1;
+ m_v = new tsparse2d*[m_Nx];
+ for ( int i=0 ; i(m_Ny, m_Nz);
+ // setup_fast();
+ // deep copy of all elements
+ for ( int i=0 ; i=m_lx && i<=m_ux ) return;
+
+ // is it empty? add a single row
+ if ( m_lx>m_ux ) {
+ m_v = new tsparse2d*[1];
+ (*m_v) = new tsparse2d(m_Ny, m_Nz, m_Ny, m_Nz, 0, 0);
+ m_lx = m_ux = i;
+ return;
+ }
+
+ tsparse2d** newv = new tsparse2d*[ ( i** newp = newv;
+ tsparse2d** mvp = m_v;
+
+ int lx = m_lx;
+ int ux = m_ux;
+
+ // grow at front
+ for ( ; m_lx>i ; m_lx-- ) (*newp++) = new tsparse2d(m_Ny, m_Nz, m_Ny, m_Nz, 0, 0);
+ // copy existing
+ for ( ; lx<=ux ; lx++ ) (*newp++) = (*mvp++);
+ // grow at back
+ for ( ; m_ux(m_Ny, m_Nz, m_Ny, m_Nz, 0, 0);
+
+ delete[] m_v;
+ m_v = newv;
+ }
+
+
+protected:
+
+ int m_Ny;
+ int m_Nz;
+
+ tsparse2d** m_v;
+
+ // T** m_fastindex;
+
+ bool m_trimmed;
+
+};
+
+
+// stream IO template
+template std::ostream& operator<<(std::ostream& s, const tsparse3d& sp) {
+ for ( int i=0 ; i0 ) create( nbins, limits );
+ else throw exception("histogram: not enough bins creating histogram: "+name);
+}
+
+
+histogram::histogram( const std::string& name, const std::vector& limits ) : mname(name) {
+ if ( limits.size()>1 ) create( limits.size()-1, &limits[0] );
+ else throw exception("histogram: not enough bin limits creating histogram: "+name);
+}
+
+
+histogram::histogram( const std::string& name, size_t nbins, double lo, double hi ) : mname(name) {
+ if ( nbins==0 ) throw exception("histogram: not enough bins creating histogram: "+name);
+ std::vector limits(nbins+1,0);
+ for ( size_t i=nbins ; i-- ; ) limits[i] = (lo*(nbins-i) + hi*i)/nbins;
+ /// make sure the limits at least are exact
+ limits[0] = lo;
+ limits[nbins] = hi;
+ create( nbins, &limits[0] );
+}
+
+
+/// constructor from stream
+
+histogram::histogram( const std::vector stream ) {
+ deserialise( stream );
+ mx.resize(my.size());
+ for ( size_t i=my.size() ; i-- ; ) mx[i] = 0.5*(mxlimits[i]+mxlimits[i+1]);
+}
+
+
+histogram::histogram( const histogram& h ) :
+ mname(h.mname),
+ mxlimits(h.mxlimits),
+ mx(h.mx),
+ my(h.my),
+ mye(h.mye),
+ myelo(h.myelo)
+{ }
+
+
+/// merge bins i and bin i+1
+
+void histogram::merge_bins( size_t i, bool norm ) {
+
+ // std::cout << "before:\n" << *this << std::endl;
+
+ size_t n = my.size();
+
+ if ( n<2 || i>=n-1 ) {
+ std::cerr << "app::grid::merge_bins() cannot merge" << std::endl;
+ return;
+ }
+
+ int j=i+1;
+
+ double w0 = mxlimits[i+1]-mxlimits[i];
+ double w1 = mxlimits[j+1]-mxlimits[j];
+ double w = mxlimits[j+1]-mxlimits[i];
+
+ if ( !norm ) w0 = w1 = w = 1;
+
+ mx[i] = 0.5*(mxlimits[j+1] + mxlimits[i]);
+ my[i] = (my[i]*w0 + my[j]*w1)/w;
+ mye[i] = std::sqrt( w0*w0*mye[i]*mye[i] + w1*w1*mye[j]*mye[j] )/w;
+
+ mxlimits.erase( mxlimits.begin()+j );
+ mx.erase( mx.begin()+j );
+ my.erase( my.begin()+j );
+ mye.erase( mye.begin()+j );
+
+ if ( myelo.size()>0 ) {
+ myelo[i] = std::sqrt( w0*w0*myelo[i]*myelo[i] + w1*w1*myelo[j]*myelo[j] )/w;
+ myelo.erase( myelo.begin()+j );
+ }
+
+ // std::cout << "after:\n" << *this << std::endl;
+
+}
+
+
+#if 0
+void histogram::merge_bin( int i ) {
+
+ std::cout << "before:\n" << *this << std::endl;
+
+ size_t n = my.size();
+ size_t nn = my.size()-1;
+
+
+ double w0 = mxlimits[n-1]-mxlimits[n-2];
+ double w1 = mxlimits[n]-mxlimits[n-1];
+
+ double w = mxlimits[n]-mxlimits[n-2];
+
+ if ( n<2 ) return;
+
+ mxlimits[nn] = mxlimits[n];
+ mxlimits.resize(n);
+
+ mx[nn-1] = 0.5*(mxlimits[nn] - mxlimits[nn-1]);
+
+ mx.resize(nn);
+
+ my[nn-1] = (my[nn-1]*w0 + my[n-1]*w1)/w;
+ my.resize(nn);
+
+ mye[nn-1] = std::sqrt( w0*w0*mye[nn-1]*mye[nn-1] + w1*w1*mye[n-1]*mye[n-1] )/w;
+ mye.resize(nn);
+
+ if ( myelo.size()>0 ) {
+ myelo[nn-1] = std::sqrt( w0*w0*myelo[nn-1]*myelo[nn-1] + w1*w1*myelo[n-1]*myelo[n-1] )/w;
+ myelo.resize(nn);
+ }
+
+ std::cout << "after:\n" << *this << std::endl;
+
+}
+#endif
+
+
+histogram& histogram::operator=( histogram&& h ) {
+ std::swap( mname, h.mname );
+
+ std::swap( mxlimits, h.mxlimits );
+ std::swap( mx, h.mx );
+
+ std::swap( my, h.my );
+ std::swap( mye, h.mye );
+ std::swap( myelo, h.myelo );
+
+ return *this;
+}
+
+
+histogram& histogram::operator=( const histogram& h ) {
+ mname = h.mname;
+ mxlimits = h.mxlimits;
+ mx = h.mx;
+ my = h.my;
+ mye = h.mye;
+ myelo = h.myelo;
+ return *this;
+}
+
+
+void histogram::fill( double x, double w ) {
+ int ibin = index( x );
+ if ( ibin != -1 ) {
+ my[ibin] += w;
+ mye[ibin] = std::sqrt( mye[ibin]*mye[ibin] + w );
+ }
+ if ( myelo.size()>0 ) myelo[ibin] = std::sqrt( myelo[ibin]*mye[ibin] + w );
+}
+
+
+/// serialisation and deserialisation ...
+
+void histogram::serialise_internal( std::vector& s ) const {
+ SB::serialise( s, name() );
+ SB::serialise( s, mxlimits );
+ SB::serialise( s, my );
+ std::vector ye(mye);
+ if ( myelo.size()>0 ) ye.insert( ye.end(), myelo.begin(), myelo.end() );
+ SB::serialise( s, ye );
+}
+
+
+void histogram::deserialise_internal( std::vector::const_iterator& itr ) {
+ SB::deserialise( itr, mname );
+ SB::deserialise( itr, mxlimits );
+ SB::deserialise( itr, my );
+ std::vector ye;
+ SB::deserialise( itr, ye );
+ if ( ye.size()==my.size() ) mye = ye;
+ else {
+ mye.clear();
+ mye.insert( mye.begin(), ye.begin(), ye.begin()+(ye.size()/2) );
+ myelo.insert( myelo.begin(), ye.begin()+(ye.size()/2), ye.end() );
+ }
+}
+
+
+
+void histogram::create( size_t nbins, const double* limits ) {
+
+ mxlimits.resize(nbins+1);
+
+ for ( size_t i=nbins+1 ; i-- ; ) mxlimits[i] = limits[i];
+
+ mx.resize(nbins);
+
+ for ( size_t i=nbins ; i-- ; ) mx[i] = 0.5*(limits[i]+limits[i+1]);
+
+ my = std::vector(nbins,0);
+ mye = std::vector(nbins,0);
+
+}
+
+
+
+void histogram::set( const std::vector& v,
+ const std::vector& ve,
+ const std::vector& velo )
+{
+ if ( v.size()!=size() ) throw exception( "histogram: number of histogram and value bins don't match" );
+ my = v;
+
+ if ( ve.size()!=0 ) {
+ if ( ve.size()!=size() ) throw exception( "histogram: number of histogram and value bins don't match" );
+ mye = ve;
+ }
+ else mye = std::vector(v.size(),0);
+
+ if ( velo.size()!=0 ) {
+ if ( velo.size()!=size() ) throw exception( "histogram: number of histogram and value bins don't match" );
+ myelo = velo;
+ }
+ else myelo.clear();
+}
+
+
+void histogram::set_errors( const std::vector& ve,
+ const std::vector& velo )
+{
+ if ( ve.size()!=size() ) throw exception( "histogram: number of histogram and value bins don't match" );
+ mye = ve;
+
+ if ( velo.size()!=0 ) {
+ if ( velo.size()!=size() ) throw exception( "histogram: number of histogram and value bins don't match" );
+ myelo = velo;
+ }
+ else myelo.clear();
+
+}
+
+
+
+
+
+
+
+
Index: applgrid/tags/applgrid-01-06-25/src/dis_pdf.h
===================================================================
--- applgrid/tags/applgrid-01-06-25/src/dis_pdf.h (revision 0)
+++ applgrid/tags/applgrid-01-06-25/src/dis_pdf.h (revision 1980)
@@ -0,0 +1,59 @@
+/** emacs: this is -*- c++ -*- **/
+/**
+ ** @file dis_pdf.h
+ **
+ ** @brief pdf transform functions for nlojet
+ **
+ ** @author mark sutton
+ ** @date Mon Dec 10 01:36:04 GMT 2007
+ **
+ ** @copyright (C) 2002-2019 mark sutton (sutt @ cern.ch)
+ **
+ ** $Id: dis_pdf.h, v1.0 Mon Dec 10 01:36:04 GMT 2007 sutt $
+ **
+ **/
+
+#ifndef __DIS_PDF_H
+#define __DIS_PDF_H
+
+#include "appl_grid/appl_pdf.h"
+
+
+
+// nlojet pdf combination class
+
+class dis_pdf : public appl::appl_pdf {
+
+public:
+
+ dis_pdf() : appl_pdf("dis") { m_Nproc=3; }
+
+ void evaluate(const double* fA, const double* fB, double* H) const;
+
+};
+
+
+inline void dis_pdf::evaluate(const double* fA, const double* , double* H) const {
+
+ fA += 6; // offset internal pointers so can use fA[-6]..fA[6] for simplicity
+
+ double U=0, D=0, G=fA[0];
+
+ // don't include the gluon (fA[0], fB[0])
+
+ for ( int i=1 ; i<=5 ; i+=2 ) D += fA[i]+fA[-i] ;
+ for ( int i=2 ; i<=6 ; i+=2 ) U += fA[i]+fA[-i] ;
+
+ H[0] = G ;
+ H[1] = D+U ;
+ H[2] = (D+4*U)/9 ;
+}
+
+
+extern "C" void fdis_pdf__(const double* fA, const double* fB, double* H);
+
+
+#endif // __DIS_PDF_H
+
+
+
Index: applgrid/tags/applgrid-01-06-25/src/mcfmzjet_pdf.h
===================================================================
--- applgrid/tags/applgrid-01-06-25/src/mcfmzjet_pdf.h (revision 0)
+++ applgrid/tags/applgrid-01-06-25/src/mcfmzjet_pdf.h (revision 1980)
@@ -0,0 +1,124 @@
+/** emacs: this is -*- c++ -*- **/
+/**
+ ** @file mcfmzjet_pdf.h
+ **
+ ** @brief pdf transform functions
+ **
+ ** @author mark sutton
+ ** @date Mon Dec 10 01:36:04 GMT 2007
+ **
+ ** @copyright (C) 2002-2019 mark sutton (sutt @ cern.ch)
+ **
+ ** $Id: mcfmz_pdf.h, v1.0 Mon Dec 10 01:36:04 GMT 2007 sutt $
+ **
+ **/
+
+
+#ifndef __MCFMZJET_PDF_H
+#define __MCFMZJET_PDF_H
+
+#include
+
+#include "appl_grid/appl_pdf.h"
+
+
+//
+// MCFM Z-jet production
+//
+class mcfmzjet_pdf : public appl::appl_pdf {
+
+public:
+ mcfmzjet_pdf() : appl_pdf("mcfm-zjet") { m_Nproc = 33; }
+
+ ~mcfmzjet_pdf() { }
+
+ virtual void evaluate(const double* fA, const double* fB, double* H) const;
+
+};
+
+
+// actual funtion to evaluate the pdf combinations
+
+inline void mcfmzjet_pdf::evaluate(const double* fA, const double* fB, double* H) const {
+ // const int nQuark = 6;
+ const int iQuark = 5;
+
+ // offset psd ptrs so can use [-6..6] indexing rather than [nQuark+..] indexing
+ fA += 6;
+ fB += 6;
+
+ double GA=fA[0];
+ double GB=fB[0];
+
+ double UpA=0; double UpB=0; double DnA=0; double DnB=0;
+ double UpbarA=0; double UpbarB=0; double DnbarA=0; double DnbarB=0;
+
+ for (int i = 1; i <= iQuark; i++) {
+ if ((i % 2) == 0) {
+ UpA += fA[i];
+ UpB += fB[i];
+ }else{
+ DnA += fA[i];
+ DnB += fB[i];
+ }
+ }
+
+ for(int i = -iQuark; i < 0; i++) {
+ if (((int)(std::abs((double)i)) % 2) == 0) {
+ UpbarA += fA[i];
+ UpbarB += fB[i];
+ } else {
+ DnbarA += fA[i];
+ DnbarB += fB[i];
+ }
+ }
+
+ // zero H first
+ for(int i = 0; i0 ? i%2 : abs(i%2)+2 );
+ int Choice = choice[i];
+ H[Choice] += fA[i] * fB[-i];
+ }
+
+ H[4] = GA * UpB;
+ H[5] = GA * UpbarB;
+ H[6] = GA * DnB;
+ H[7] = GA * DnbarB;
+ H[8] = UpA *GB;
+ H[9] = UpbarA *GB;
+ H[10]= DnA *GB;
+ H[11]= DnbarA *GB;
+
+ H[12] = GA*GB;
+
+
+ for(int i = -iQuark; i <= iQuark; i++) {
+ if(i==0) continue;
+ for(int j = -iQuark; j <= iQuark; j++) {
+ if(j==0 || i==j || i==-j) continue;
+ H[13+choice[i]+4*choice[j]] += fA[i] * fB[j];
+ }
+ }
+
+ for(int i = -iQuark; i <= iQuark; i++) {
+ if (i == 0) continue;
+ // int Choice = ( i>0 ? i%2 : abs(i%2)+2 );
+ int Choice = choice[i];
+ H[29+Choice] += fA[i] * fB[i];
+ }
+
+
+ return;
+}
+
+
+extern "C" void fmcfmzjet_pdf__(const double* fA, const double* fB, double* H);
+
+#endif // __MCFMZJET_PDF_H
+
Index: applgrid/tags/applgrid-01-06-25/src/appl_igrid.h
===================================================================
--- applgrid/tags/applgrid-01-06-25/src/appl_igrid.h (revision 0)
+++ applgrid/tags/applgrid-01-06-25/src/appl_igrid.h (revision 1980)
@@ -0,0 +1,818 @@
+/** emacs: this is -*- c++ -*- **/
+/**
+ ** @file appl_igrid.h
+ **
+ ** @brief grid class header - all the functions needed to create and
+ ** fill the grid from an NLO calculation program, based on the
+ ** class from D.Clements and T.Carli.
+ **
+ ** See contribution from T.Carli et al from the HERA-LHC
+ ** workshop - working group 3
+ **
+ ** @author mark sutton
+ ** @date Fri Nov 2 05:31:03 GMT 2007
+ **
+ ** @copyright (C) 2002-2019 mark sutton (sutt @ cern.ch)
+ **
+ ** $Id: appl_igrid.h, v2.00 Fri Nov 2 05:31:03 GMT 2007 sutt $
+ **
+ **/
+
+// Fixme: this needs to be tidied up. eg there are too many different,
+// and too many version of, accessors for x/y, Q2/tau etc there
+// should be only one set, for x and Q2 *or* y and tau, but
+// not both. In fact maybe they should be for x and Q2, since y
+// and tau should perhaps be purely an internal grid issue of no
+// concern for the user.
+
+
+#ifndef __APPL_IGRID_H
+#define __APPL_IGRID_H
+
+#include
+#include
+#include
+#include