Index: trunk/HiggsSignals-2/configure =================================================================== --- trunk/HiggsSignals-2/configure (revision 599) +++ trunk/HiggsSignals-2/configure (revision 600) @@ -1,136 +1,136 @@ #! /bin/sh #configure script for HiggsSignals #last modified 20.12.18 TS # this tells the fortran program where to look to read in tables of experimental data chmod u+x ./create_store_pathname_HS.bat ./create_store_pathname_HS.bat > store_pathname_HS.f90 echo ' ******************************************************** ' echo ' Configuring HiggsSignals...' 1>&2 echo ' ******************************************************** ' echo ' You can specify the correct paths to HiggsBounds and ' echo ' (if needed) FeynHiggs via the options ' echo ' --hbpath=/PATH/TO/HIGGSBOUNDS' echo ' --fhpath=/PATH/TO/FEYNHIGGS' echo ' or manually in the configure script.' echo '' echo ' The Fortran compiler has to be set manually in the ' echo ' configure script. The default is gfortran.' echo ' ******************************************************** ' #---------------------------- # Note that you have to build HiggsBounds successfully BEFORE you build HiggsSignals. -CONF_HBPATH="/Users/timstefaniak/Work/HiggsBounds/trunk/HiggsBounds-5" +CONF_HBPATH="/Users/timstefaniak/Work/higgsbounds_SVN/trunk/HiggsBounds-5" CONF_FHPATH="/Users/timstefaniak/Work/Codes/FeynHiggs-2.14.3" # For FeynHiggs: CONF_OS=`uname -s` CONF_MACH=`uname -m` CONF_DEFPREFIX="$CONF_MACH-$CONF_OS" for arg in "$@" ; do case "$arg" in --hbpath=*) CONF_HBPATH=`expr "$arg" : ".*--hbpath=\(.*\)"` ;; --fhpath=*) CONF_FHPATH=`expr "$arg" : ".*--fhpath=\(.*\)"` ;; -*) echo "Warning: $arg is not a valid option." ;; *=*) eval `echo $arg\" | sed 's/=/="/'` ;; *) echo "Warning: $arg is not a valid argument." ;; esac done # Remove slash in entered paths if [ "${CONF_HBPATH: -1}" = "/" ]; then CONF_HBPATH="${CONF_HBPATH%?}" fi if [ "${CONF_FHPATH: -1}" = "/" ]; then CONF_FHPATH="${CONF_FHPATH%?}" fi echo " hbpath=${CONF_HBPATH}" echo " fhpath=${CONF_FHPATH}" #---------------------------- cat - makefile.in > makefile << _EOF_ # --- variables defined by configure --- # If you want to edit these variables, change ./configure, not ./makefile #---------------------------- # Put you favourite compiler and compiler options in here: F90C = gfortran F77C = gfortran F90FLAGS = -fbounds-check -ffixed-line-length-none #F90FLAGS = -fbounds-check -ffixed-line-length-none -fPIC #F90FLAGS = -fbounds-check -Wall -W #other useful gfortran flags: # -pg can then use gprof ./programname to look at how efficient program is # -fopenmp switches on parallel processing (see HiggsBounds.F90) # -O3 sets a high level of optimisation #F90C = g95 #F77C = g95 #F90FLAGS = -fbounds-check -ffree-line-length-huge #F90C = pgf90 #F77C = pgf90 #F90FLAGS = -C -Ktrap=fp -Mbounds #F90C = f95 #F77C = f95 #F90FLAGS = -C #F90C = ifort #F77C = ifort #F90FLAGS = -C #if using the NAG compiler, you need the compiler flag -DNAGf90Fortran #because the subroutines flush,iargc,getarg need to use modules #caution: the use of the NAG compiler is unsupported HBLIBS =-L${CONF_HBPATH} HBINCLUDE =-I${CONF_HBPATH} #------------------------------ # If you wish to use HiggsSignals in conjunction with FeynHiggs, # make sure these paths indicate where you've stored these packages on your system, FHINCLUDE = -I${CONF_FHPATH}/${CONF_DEFPREFIX}/include FHLIBS = -L${CONF_FHPATH}/${CONF_DEFPREFIX}/lib -lFH # Note that care must be taken to ensure that consistent compilers are used # If you get error messages similar to # ... undefined reference to initialize_higgsbounds__... # and the paths the libraries are correct, it's a good idea to check # that the same compilers are used for each package. # (you may be able to work around this # e.g. if FeynHiggs was compiled with g77 it may be possible to compile HiggsBounds # with gfortran and add -lg2c to the end of the variable FHLIBS) #----------------------------- HSLIBS = -L./ -lHS EXE = HiggsSignals #----------------------------- _EOF_ echo ' ******************************************************** ' echo ' ...finished configure script.' 1>&2 echo ' ******************************************************** ' Index: trunk/HiggsSignals-2/XScov_13TeV.in =================================================================== --- trunk/HiggsSignals-2/XScov_13TeV.in (revision 0) +++ trunk/HiggsSignals-2/XScov_13TeV.in (revision 600) @@ -0,0 +1,11 @@ + 0.00254666 4.39899e-07 4.27442e-07 1.63212e-06 0.00113514 0.00254666 2.98794e-06 2.21462e-06 -1.35343e-06 4.22512e-07 9.46349e-06 + 4.39899e-07 0.00139354 0.000390055 0.000329841 -2.4179e-06 4.39899e-07 3.8443e-06 -3.10503e-06 -2.45439e-07 0.000390201 1.69209e-08 + 4.27442e-07 0.000390055 0.000396957 0.000507496 -1.38692e-06 4.27442e-07 -3.38238e-06 -1.93738e-07 2.21157e-07 0.00039392 0.00130827 + 1.63212e-06 0.000329841 0.000507496 0.00142065 -6.65353e-07 1.63212e-06 -1.12692e-05 1.70796e-06 8.99267e-07 0.000491329 0.00752954 + 0.00113514 -2.4179e-06 -1.38692e-06 -6.65353e-07 0.00698249 0.00113514 -1.78693e-05 8.44596e-06 1.75528e-06 -1.40258e-06 4.74238e-06 + 0.00254666 4.39899e-07 4.27442e-07 1.63212e-06 0.00113514 0.00254666 2.98794e-06 2.21462e-06 -1.35343e-06 4.22512e-07 9.46349e-06 + 2.98794e-06 3.8443e-06 -3.38238e-06 -1.12692e-05 -1.78693e-05 2.98794e-06 0.0485626 8.74098e-06 3.00546e-06 -3.21038e-06 -6.82845e-05 + 2.21462e-06 -3.10503e-06 -1.93738e-07 1.70796e-06 8.44596e-06 2.21462e-06 8.74098e-06 0.0131502 -6.68913e-07 -2.12361e-07 1.60028e-05 + -1.35343e-06 -2.45439e-07 2.21157e-07 8.99267e-07 1.75528e-06 -1.35343e-06 3.00546e-06 -6.68913e-07 0.000931484 2.06182e-07 6.40571e-06 + 4.22512e-07 0.000390201 0.00039392 0.000491329 -1.40258e-06 4.22512e-07 -3.21038e-06 -2.12361e-07 2.06182e-07 0.000391226 0.00120381 + 9.46349e-06 1.69209e-08 0.00130827 0.00752954 4.74238e-06 9.46349e-06 -6.82845e-05 1.60028e-05 6.40571e-06 0.00120381 0.0494146 \ No newline at end of file Index: trunk/HiggsSignals-2/XScovSM_13TeV.in =================================================================== --- trunk/HiggsSignals-2/XScovSM_13TeV.in (revision 0) +++ trunk/HiggsSignals-2/XScovSM_13TeV.in (revision 600) @@ -0,0 +1,11 @@ + 0.00254666 4.39899e-07 4.27442e-07 1.63212e-06 0.00113514 0.00254666 2.98794e-06 2.21462e-06 -1.35343e-06 4.22512e-07 9.46349e-06 + 4.39899e-07 0.00139354 0.000390055 0.000329841 -2.4179e-06 4.39899e-07 3.8443e-06 -3.10503e-06 -2.45439e-07 0.000390201 1.69209e-08 + 4.27442e-07 0.000390055 0.000396957 0.000507496 -1.38692e-06 4.27442e-07 -3.38238e-06 -1.93738e-07 2.21157e-07 0.00039392 0.00130827 + 1.63212e-06 0.000329841 0.000507496 0.00142065 -6.65353e-07 1.63212e-06 -1.12692e-05 1.70796e-06 8.99267e-07 0.000491329 0.00752954 + 0.00113514 -2.4179e-06 -1.38692e-06 -6.65353e-07 0.00698249 0.00113514 -1.78693e-05 8.44596e-06 1.75528e-06 -1.40258e-06 4.74238e-06 + 0.00254666 4.39899e-07 4.27442e-07 1.63212e-06 0.00113514 0.00254666 2.98794e-06 2.21462e-06 -1.35343e-06 4.22512e-07 9.46349e-06 + 2.98794e-06 3.8443e-06 -3.38238e-06 -1.12692e-05 -1.78693e-05 2.98794e-06 0.0485626 8.74098e-06 3.00546e-06 -3.21038e-06 -6.82845e-05 + 2.21462e-06 -3.10503e-06 -1.93738e-07 1.70796e-06 8.44596e-06 2.21462e-06 8.74098e-06 0.0131502 -6.68913e-07 -2.12361e-07 1.60028e-05 + -1.35343e-06 -2.45439e-07 2.21157e-07 8.99267e-07 1.75528e-06 -1.35343e-06 3.00546e-06 -6.68913e-07 0.000931484 2.06182e-07 6.40571e-06 + 4.22512e-07 0.000390201 0.00039392 0.000491329 -1.40258e-06 4.22512e-07 -3.21038e-06 -2.12361e-07 2.06182e-07 0.000391226 0.00120381 + 9.46349e-06 1.69209e-08 0.00130827 0.00752954 4.74238e-06 9.46349e-06 -6.82845e-05 1.60028e-05 6.40571e-06 0.00120381 0.0494146 \ No newline at end of file Index: trunk/HiggsSignals-2/example_data/SLHA/MSSMexample.fh.1 =================================================================== --- trunk/HiggsSignals-2/example_data/SLHA/MSSMexample.fh.1 (revision 599) +++ trunk/HiggsSignals-2/example_data/SLHA/MSSMexample.fh.1 (revision 600) @@ -1,481 +1,481 @@ BLOCK SPINFO 1 FeynHiggs 2 2.14.3 2 built on Aug 01, 2018 BLOCK MODSEL 1 0 # Model 2 1 # GridPts 3 0 # Content 4 0 # RPV 5 0 # CPV 6 0 # FV BLOCK SMINPUTS 1 1.28952896E+02 # invAlfaMZ 2 1.16637002E-05 # GF 3 1.18000000E-01 # AlfasMZ 4 9.11876000E+01 # MZ 5 4.17999983E+00 # Mb 6 1.72500000E+02 # Mt 7 1.77703000E+00 # Mtau 11 5.10998902E-04 # Me 13 1.05658357E-01 # Mmu 21 6.00000000E-03 # Md 22 3.00000000E-03 # Mu 23 9.50000000E-02 # Ms 24 1.28600000E+00 # Mc BLOCK MINPAR 3 2.00000000E+01 # TB BLOCK EXTPAR 0 0.00000000E+00 # Q 1 1.80000000E+02 # M1 2 3.00000000E+02 # M2 3 2.50000000E+03 # M3 11 2.85000000E+03 # At 12 2.85000000E+03 # Ab 13 8.00000000E+02 # Atau 23 1.00000000E+03 # MUE 25 2.00000000E+01 # TB 26 1.00000000E+03 # MA0 27 1.00322567E+03 # MHp 31 2.00000000E+03 # MSL(1) 32 2.00000000E+03 # MSL(2) 33 3.50000000E+02 # MSL(3) 34 2.00000000E+03 # MSE(1) 35 2.00000000E+03 # MSE(2) 36 3.50000000E+02 # MSE(3) 41 2.00000000E+03 # MSQ(1) 42 2.00000000E+03 # MSQ(2) 43 1.50000000E+03 # MSQ(3) 44 2.00000000E+03 # MSU(1) 45 2.00000000E+03 # MSU(2) 46 1.50000000E+03 # MSU(3) 47 2.00000000E+03 # MSD(1) 48 2.00000000E+03 # MSD(2) 49 1.50000000E+03 # MSD(3) BLOCK MASS 1000012 1.99896552E+03 # MSf(1,1,1) 1000011 2.00057313E+03 # MSf(1,2,1) 2000011 2.00046095E+03 # MSf(2,2,1) 1000002 1.99927301E+03 # MSf(1,3,1) 2000002 1.99969262E+03 # MSf(2,3,1) 1000001 2.00088111E+03 # MSf(1,4,1) 2000001 2.00015290E+03 # MSf(2,4,1) 1000014 1.99896552E+03 # MSf(1,1,2) 1000013 2.00097334E+03 # MSf(1,2,2) 2000013 2.00006064E+03 # MSf(2,2,2) 1000004 1.99855848E+03 # MSf(1,3,2) 2000004 2.00040758E+03 # MSf(2,3,2) 1000003 2.00103927E+03 # MSf(1,4,2) 2000003 1.99999468E+03 # MSf(2,4,2) 1000016 3.44039455E+02 # MSf(1,1,3) 1000015 3.00752176E+02 # MSf(1,2,3) 2000015 3.98360763E+02 # MSf(2,2,3) 1000006 1.33965932E+03 # MSf(1,3,3) 2000006 1.66153199E+03 # MSf(2,3,3) 1000005 1.48076973E+03 # MSf(1,4,3) 2000005 1.52035582E+03 # MSf(2,4,3) 25 1.25446331E+02 # Mh0 35 9.99859026E+02 # MHH 36 1.00000000E+03 # MA0 37 1.00352099E+03 # MHp 1000022 1.79456889E+02 # MNeu(1) 1000023 2.97205684E+02 # MNeu(2) 1000025 -1.00294285E+03 # MNeu(3) 1000035 1.00628028E+03 # MNeu(4) 1000024 2.97193887E+02 # MCha(1) 1000037 1.00727320E+03 # MCha(2) 1000021 2.50000000E+03 # MGl BLOCK DMASS 0 1.72500000E+02 # Q 25 1.98374186E+00 # Delta Mh0 35 3.20482021E-02 # Delta MHH 36 0.00000000E+00 # Delta MA0 37 5.13172661E-02 # Delta MHp BLOCK NMIX 1 1 9.98876349E-01 # ZNeu(1,1) 1 2 -8.40320482E-03 # ZNeu(1,2) 1 3 4.54802158E-02 # ZNeu(1,3) 1 4 -1.03429432E-02 # ZNeu(1,4) 2 1 1.27116592E-02 # ZNeu(2,1) 2 2 9.95548058E-01 # ZNeu(2,2) 2 3 -8.83657886E-02 # ZNeu(2,3) 2 4 3.02318502E-02 # ZNeu(2,4) 3 1 2.44324084E-02 # ZNeu(3,1) 3 2 -4.14182340E-02 # ZNeu(3,2) 3 3 -7.05140858E-01 # ZNeu(3,3) 3 4 -7.07434773E-01 # ZNeu(3,4) 4 1 -3.85682707E-02 # ZNeu(4,1) 4 2 8.42495085E-02 # ZNeu(4,2) 4 3 7.02067951E-01 # ZNeu(4,3) 4 4 -7.06056019E-01 # ZNeu(4,4) BLOCK UMIX 1 1 9.92105591E-01 # UCha(1,1) 1 2 -1.25405328E-01 # UCha(1,2) 2 1 1.25405328E-01 # UCha(2,1) 2 2 9.92105591E-01 # UCha(2,2) BLOCK VMIX 1 1 9.99077591E-01 # VCha(1,1) 1 2 -4.29414463E-02 # VCha(1,2) 2 1 4.29414463E-02 # VCha(2,1) 2 2 9.99077591E-01 # VCha(2,2) BLOCK STAUMIX 1 1 7.04779339E-01 # USf(1,1) 1 2 7.09426587E-01 # USf(1,2) 2 1 7.09426587E-01 # USf(2,1) 2 2 -7.04779339E-01 # USf(2,2) BLOCK STOPMIX 1 1 7.07720635E-01 # USf(1,1) 1 2 -7.06492394E-01 # USf(1,2) 2 1 7.06492394E-01 # USf(2,1) 2 2 7.07720635E-01 # USf(2,2) BLOCK SBOTMIX 1 1 6.98400871E-01 # USf(1,1) 1 2 7.15706800E-01 # USf(1,2) 2 1 7.15706800E-01 # USf(2,1) 2 2 -6.98400871E-01 # USf(2,2) BLOCK ALPHA -5.11239780E-02 # Alpha BLOCK DALPHA 7.94658769E-05 # Delta Alpha BLOCK HMIX Q= -0.99900000E+03 1 1.00000000E+03 # MUE 2 2.00000000E+01 # TB BLOCK MSOFT Q= 0.00000000E+00 1 1.80000000E+02 # M1 2 3.00000000E+02 # M2 3 2.50000000E+03 # M3 31 2.00000000E+03 # MSL(1) 32 2.00000000E+03 # MSL(2) 33 3.50000000E+02 # MSL(3) 34 2.00000000E+03 # MSE(1) 35 2.00000000E+03 # MSE(2) 36 3.50000000E+02 # MSE(3) 41 2.00000000E+03 # MSQ(1) 42 2.00000000E+03 # MSQ(2) 43 1.50000000E+03 # MSQ(3) 44 2.00000000E+03 # MSU(1) 45 2.00000000E+03 # MSU(2) 46 1.50000000E+03 # MSU(3) 47 2.00000000E+03 # MSD(1) 48 2.00000000E+03 # MSD(2) 49 1.50000000E+03 # MSD(3) BLOCK AE Q= 0.00000000E+00 1 1 2.85000000E+03 # Af(1,1) 2 2 2.85000000E+03 # Af(2,2) 3 3 8.00000000E+02 # Af(3,3) BLOCK AU Q= 0.00000000E+00 1 1 2.85000000E+03 # Af(1,1) 2 2 2.85000000E+03 # Af(2,2) 3 3 2.85000000E+03 # Af(3,3) BLOCK AD Q= 0.00000000E+00 1 1 2.85000000E+03 # Af(1,1) 2 2 2.85000000E+03 # Af(2,2) 3 3 2.85000000E+03 # Af(3,3) BLOCK YE Q= 0.00000000E+00 1 1 5.87736719E-05 # Yf(1,1) 2 2 1.21525302E-02 # Yf(2,2) 3 3 2.04389046E-01 # Yf(3,3) BLOCK YU Q= 0.00000000E+00 1 1 1.72525826E-05 # Yf(1,1) 2 2 7.39560709E-03 # Yf(2,2) 3 3 9.92023502E-01 # Yf(3,3) BLOCK YD Q= 0.00000000E+00 1 1 6.35951590E-04 # Yf(1,1) 2 2 1.00679551E-02 # Yf(2,2) 3 3 3.98258811E-01 # Yf(3,3) BLOCK VCKMIN 1 2.25300000E-01 # lambda 2 8.08000000E-01 # A 3 1.32000000E-01 # rhobar 4 3.41000000E-01 # etabar BLOCK MSL2 Q= 0.00000000E+00 1 1 4.00000000E+06 # MSL2(1,1) 2 2 4.00000000E+06 # MSL2(2,2) 3 3 1.22500000E+05 # MSL2(3,3) BLOCK MSE2 Q= 0.00000000E+00 1 1 4.00000000E+06 # MSE2(1,1) 2 2 4.00000000E+06 # MSE2(2,2) 3 3 1.22500000E+05 # MSE2(3,3) BLOCK MSQ2 Q= 0.00000000E+00 1 1 4.00000000E+06 # MSQ2(1,1) 2 2 4.00000000E+06 # MSQ2(2,2) 3 3 2.25000000E+06 # MSQ2(3,3) BLOCK MSU2 Q= 0.00000000E+00 1 1 4.00000000E+06 # MSU2(1,1) 2 2 4.00000000E+06 # MSU2(2,2) 3 3 2.25000000E+06 # MSU2(3,3) BLOCK MSD2 Q= 0.00000000E+00 1 1 4.00000000E+06 # MSD2(1,1) 2 2 4.00000000E+06 # MSD2(2,2) 3 3 2.25000000E+06 # MSD2(3,3) BLOCK TE Q= 0.00000000E+00 1 1 1.67504965E-01 # Tf(1,1) 2 2 3.46347112E+01 # Tf(2,2) 3 3 1.63511237E+02 # Tf(3,3) BLOCK TU Q= 0.00000000E+00 1 1 4.91698605E-02 # Tf(1,1) 2 2 2.10774802E+01 # Tf(2,2) 3 3 2.82726698E+03 # Tf(3,3) BLOCK TD Q= 0.00000000E+00 1 1 1.81246203E+00 # Tf(1,1) 2 2 2.86936721E+01 # Tf(2,2) 3 3 1.13503761E+03 # Tf(3,3) BLOCK SELMIX 1 1 9.99809281E-01 # UASf(1,1) 1 4 -1.95295227E-02 # UASf(1,4) 2 2 7.49270774E-01 # UASf(2,2) 2 5 -6.62263775E-01 # UASf(2,5) 3 3 7.04779339E-01 # UASf(3,3) 3 6 7.09426587E-01 # UASf(3,6) 4 1 1.95295227E-02 # UASf(4,1) 4 4 9.99809281E-01 # UASf(4,4) 5 2 6.62263775E-01 # UASf(5,2) 5 5 7.49270774E-01 # UASf(5,5) 6 3 7.09426587E-01 # UASf(6,3) 6 6 -7.04779339E-01 # UASf(6,6) BLOCK USQMIX 1 1 9.99987470E-01 # UASf(1,1) 1 4 -5.00595551E-03 # UASf(1,4) 2 2 7.83236011E-01 # UASf(2,2) 2 5 -6.21724498E-01 # UASf(2,5) 3 3 7.07720635E-01 # UASf(3,3) 3 6 -7.06492394E-01 # UASf(3,6) 4 1 5.00595551E-03 # UASf(4,1) 4 4 9.99987470E-01 # UASf(4,4) 5 2 6.21724498E-01 # UASf(5,2) 5 5 7.83236011E-01 # UASf(5,5) 6 3 7.06492394E-01 # UASf(6,3) 6 6 7.07720635E-01 # UASf(6,6) BLOCK DSQMIX 1 1 9.99469672E-01 # UASf(1,1) 1 4 -3.25633860E-02 # UASf(1,4) 2 2 9.20773138E-01 # UASf(2,2) 2 5 -3.90098486E-01 # UASf(2,5) 3 3 6.98400871E-01 # UASf(3,3) 3 6 7.15706800E-01 # UASf(3,6) 4 1 3.25633860E-02 # UASf(4,1) 4 4 9.99469672E-01 # UASf(4,4) 5 2 3.90098486E-01 # UASf(5,2) 5 5 9.20773138E-01 # UASf(5,5) 6 3 7.15706800E-01 # UASf(6,3) 6 6 -6.98400871E-01 # UASf(6,6) BLOCK CVHMIX 1 1 9.99999944E-01 # UH(1,1) 1 2 3.33494451E-04 # UH(1,2) 1 3 0.00000000E+00 # UH(1,3) 2 1 -3.33494451E-04 # UH(2,1) 2 2 9.99999944E-01 # UH(2,2) 2 3 0.00000000E+00 # UH(2,3) 3 1 0.00000000E+00 # UH(3,1) 3 2 0.00000000E+00 # UH(3,2) 3 3 1.00000000E+00 # UH(3,3) BLOCK PRECOBS 1 3.70801438E-02 # DeltaR 2 5.73493744E-05 # DeltaRho 3 8.03580769E+01 # MWMSSM 4 8.03532994E+01 # MWSM 5 2.31519688E-01 # SW2effMSSM 6 2.31537613E-01 # SW2effSM 11 1.55628199E-10 # gminus2mu 21 0.00000000E+00 # EDMeTh 22 0.00000000E+00 # EDMn 23 0.00000000E+00 # EDMHg 31 4.42059407E-04 # bsgammaMSSM 32 3.92583816E-04 # bsgammaSM 33 2.09482130E+01 # DeltaMsMSSM 34 2.08767455E+01 # DeltaMsSM 35 6.06876785E-08 # BsmumuMSSM 36 3.35592960E-09 # BsmumuSM DECAY 25 4.03250828E-03 # Gamma(h0) 2.31128319E-03 2 22 22 # BR(h0 -> photon photon) 1.50311516E-03 2 22 23 # BR(h0 -> photon Z) 2.81285945E-02 2 23 23 # BR(h0 -> Z Z) 2.30327715E-01 2 -24 24 # BR(h0 -> W W) 6.91739636E-02 2 21 21 # BR(h0 -> gluon gluon) 5.16700736E-09 2 -11 11 # BR(h0 -> Electron electron) 2.29837940E-04 2 -13 13 # BR(h0 -> Muon muon) 6.63595508E-02 2 -15 15 # BR(h0 -> Tau tau) 2.04401275E-07 2 -2 2 # BR(h0 -> Up up) 2.86468550E-02 2 -4 4 # BR(h0 -> Charm charm) 8.63906818E-07 2 -1 1 # BR(h0 -> Down down) 2.16956394E-04 2 -3 3 # BR(h0 -> Strange strange) 5.73101055E-01 2 -5 5 # BR(h0 -> Bottom bottom) DECAY 35 4.84215254E+00 # Gamma(HH) 7.01532225E-07 2 22 22 # BR(HH -> photon photon) 1.48943143E-06 2 22 23 # BR(HH -> photon Z) 5.88610434E-05 2 23 23 # BR(HH -> Z Z) 7.90218210E-05 2 -24 24 # BR(HH -> W W) 2.90176556E-05 2 21 21 # BR(HH -> gluon gluon) 1.26713589E-08 2 -11 11 # BR(HH -> Electron electron) 5.63997808E-04 2 -13 13 # BR(HH -> Muon muon) 1.69514757E-01 2 -15 15 # BR(HH -> Tau tau) 1.60333008E-12 2 -2 2 # BR(HH -> Up up) 2.24577074E-07 2 -4 4 # BR(HH -> Charm charm) 1.20781980E-02 2 -6 6 # BR(HH -> Top top) 1.28529222E-06 2 -1 1 # BR(HH -> Down down) 3.22736531E-04 2 -3 3 # BR(HH -> Strange strange) 7.07323382E-01 2 -5 5 # BR(HH -> Bottom bottom) 1.32833438E-02 2 -1000024 1000024 # BR(HH -> Chargino1 chargino1) 8.17648278E-04 2 1000022 1000022 # BR(HH -> neutralino1 neutralino1) 4.77219537E-03 2 1000022 1000023 # BR(HH -> neutralino1 neutralino2) 6.45783335E-03 2 1000023 1000023 # BR(HH -> neutralino2 neutralino2) 9.13887154E-04 2 25 25 # BR(HH -> h0 h0) 4.68826918E-02 2 -1000015 1000015 # BR(HH -> Stau1 stau1) 1.09393449E-06 2 -1000015 2000015 # BR(HH -> Stau1 stau2) 1.09393449E-06 2 -2000015 1000015 # BR(HH -> Stau2 stau1) 3.68631158E-02 2 -2000015 2000015 # BR(HH -> Stau2 stau2) DECAY 36 4.86015881E+00 # Gamma(A0) 9.57405804E-07 2 22 22 # BR(A0 -> photon photon) 2.82671580E-06 2 22 23 # BR(A0 -> photon Z) 1.28712044E-04 2 21 21 # BR(A0 -> gluon gluon) 1.24418936E-08 2 -11 11 # BR(A0 -> Electron electron) 5.53785038E-04 2 -13 13 # BR(A0 -> Muon muon) 1.66442344E-01 2 -15 15 # BR(A0 -> Tau tau) 1.53459338E-12 2 -2 2 # BR(A0 -> Up up) 2.14918378E-07 2 -4 4 # BR(A0 -> Charm charm) 1.32833979E-02 2 -6 6 # BR(A0 -> Top top) 1.26199784E-06 2 -1 1 # BR(A0 -> Down down) 3.16887390E-04 2 -3 3 # BR(A0 -> Strange strange) 6.94548917E-01 2 -5 5 # BR(A0 -> Bottom bottom) 2.19326967E-02 2 -1000024 1000024 # BR(A0 -> Chargino1 chargino1) 9.79429746E-04 2 1000022 1000022 # BR(A0 -> neutralino1 neutralino1) 6.43193791E-03 2 1000022 1000023 # BR(A0 -> neutralino1 neutralino2) 1.06677808E-02 2 1000023 1000023 # BR(A0 -> neutralino2 neutralino2) 8.51998309E-05 2 23 25 # BR(A0 -> Z h0) 6.37923825E-37 2 25 25 # BR(A0 -> h0 h0) 4.23118188E-02 2 -1000015 2000015 # BR(A0 -> Stau1 stau2) 4.23118188E-02 2 -2000015 1000015 # BR(A0 -> Stau2 stau1) DECAY 37 4.64425412E+00 # Gamma(Hp) 1.42932534E-08 2 -11 12 # BR(Hp -> Electron nu_e) 6.11081106E-04 2 -13 14 # BR(Hp -> Muon nu_mu) 1.72853411E-01 2 -15 16 # BR(Hp -> Tau nu_tau) 1.29523838E-06 2 -1 2 # BR(Hp -> Down up) 1.49763113E-05 2 -3 2 # BR(Hp -> Strange up) 7.73575675E-06 2 -5 2 # BR(Hp -> Bottom up) 7.48090776E-08 2 -1 4 # BR(Hp -> Down charm) 3.24402846E-04 2 -3 4 # BR(Hp -> Strange charm) 1.08328003E-03 2 -5 4 # BR(Hp -> Bottom charm) 1.48273253E-06 2 -1 6 # BR(Hp -> Down top) 3.27765366E-05 2 -3 6 # BR(Hp -> Strange top) 7.24519543E-01 2 -5 6 # BR(Hp -> Bottom top) 1.14825567E-02 2 1000022 1000024 # BR(Hp -> neutralino1 chargino1) 1.73930211E-06 2 1000023 1000024 # BR(Hp -> neutralino2 chargino1) 9.06810606E-05 2 24 25 # BR(Hp -> W h0) 1.92294506E-10 2 24 35 # BR(Hp -> W HH) 1.58042563E-10 2 24 36 # BR(Hp -> W A0) 4.59711000E-02 2 -1000015 1000016 # BR(Hp -> Stau1 snu_tau1) 4.30038481E-02 2 -2000015 1000016 # BR(Hp -> Stau2 snu_tau1) DECAY 6 1.35305515E+00 # Gamma(top) 1.00000000E+00 2 5 24 # BR(top -> bottom W) # Block HiggsCouplingsBosons # For exact definitions of NormEffCoup see HiggsBounds manual 0.999999 3 25 24 24 # higgs-W-W effective coupling, normalised to SM 0.116558E-02 3 35 24 24 # higgs-W-W effective coupling, normalised to SM 0.00000 3 36 24 24 # higgs-W-W effective coupling, normalised to SM 0.999999 3 25 23 23 # higgs-Z-Z effective coupling, normalised to SM 0.116558E-02 3 35 23 23 # higgs-Z-Z effective coupling, normalised to SM 0.00000 3 36 23 23 # higgs-Z-Z effective coupling, normalised to SM 0.952977 3 25 21 21 # higgs-gluon-gluon effective coupling, normalised to SM 0.501981E-01 3 35 21 21 # higgs-gluon-gluon effective coupling, normalised to SM 0.113261 3 36 21 21 # higgs-gluon-gluon effective coupling, normalised to SM 0.00000 3 25 25 23 # higgs-higgs-Z effective coupling, normalised 0.00000 3 35 25 23 # higgs-higgs-Z effective coupling, normalised 0.00000 3 35 35 23 # higgs-higgs-Z effective coupling, normalised 0.145434E-03 3 36 25 23 # higgs-higgs-Z effective coupling, normalised 0.124773 3 36 35 23 # higgs-higgs-Z effective coupling, normalised 0.00000 3 36 36 23 # higgs-higgs-Z effective coupling, normalised # Block HiggsCouplingsFermions # For exact definitions of NormEffCoup see HiggsBounds manual # ScalarNormEffCoup PseudoSNormEffCoup NP IP1 IP2 IP3 # Scalar, Pseudoscalar Normalised Effective Coupling 1.0193000674322237 1.7565995677947101E-029 3 25 5 5 # higgs-b-b eff. coupling, normalised to SM 16.557714685947133 1.5070568169751837E-026 3 35 5 5 # higgs-b-b eff. coupling, normalised to SM 1.5070578407059939E-026 16.558891516303152 3 36 5 5 # higgs-b-b eff. coupling, normalised to SM 0.99994104160649944 0.0000000000000000 3 25 6 6 # higgs-top-top eff. coupling, normalised to SM 5.1165548085686831E-002 0.0000000000000000 3 35 6 6 # higgs-top-top eff. coupling, normalised to SM 0.0000000000000000 5.0000000000000003E-002 3 36 6 6 # higgs-top-top eff. coupling, normalised to SM 1.0233109617137364 0.0000000000000000 3 25 15 15 # higgs-tau-tau eff. coupling, normalised to SM 19.998820832129994 0.0000000000000000 3 35 15 15 # higgs-tau-tau eff. coupling, normalised to SM 0.0000000000000000 20.000000000000000 3 36 15 15 # higgs-tau-tau eff. coupling, normalised to SM # Block ChargedHiggsLHC8 5 6 37 0.239848E-03 # t-b-Hp production (in pb) 4 5 37 0.00000 # c-b-Hp production (in pb) 2 5 37 0.00000 # u-b-Hp production (in pb) 3 4 37 0.00000 # c-s-Hp production (in pb) 1 4 37 0.00000 # c-d-Hp production (in pb) 1 2 37 0.00000 # u-d-Hp production (in pb) 2 3 37 0.00000 # u-s-Hp production (in pb) 0 24 37 0.00000 # W-Hp production (in pb) 0 23 37 0.00000 # Z-Hp production (in pb) 1 1 37 0.00000 # Hp VBF production (in pb) 0 -37 37 0.00000 # Hp-Hm production (in pb) 0 25 37 0.00000 # h0-Hp production (in pb) 0 35 37 0.00000 # h0-Hp production (in pb) 0 36 37 0.00000 # h0-Hp production (in pb) # Block ChargedHiggsLHC13 5 6 37 0.147262E-02 # t-b-Hp production (in pb) 4 5 37 0.00000 # c-b-Hp production (in pb) 2 5 37 0.00000 # u-b-Hp production (in pb) 3 4 37 0.00000 # c-s-Hp production (in pb) 1 4 37 0.00000 # c-d-Hp production (in pb) 1 2 37 0.00000 # u-d-Hp production (in pb) 2 3 37 0.00000 # u-s-Hp production (in pb) 0 24 37 0.00000 # W-Hp production (in pb) 0 23 37 0.00000 # Z-Hp production (in pb) 1 1 37 0.00000 # Hp VBF production (in pb) 0 -37 37 0.00000 # Hp-Hm production (in pb) 0 25 37 0.00000 # h0-Hp production (in pb) 0 35 37 0.00000 # h0-Hp production (in pb) 0 36 37 0.00000 # h0-Hp production (in pb) Block HiggsBoundsResults # results from HiggsBounds http://projects.hepforge.org/higgsbounds # HBresult : scenario allowed flag (1: allowed, 0: excluded, -1: unphysical) # chan id number: most sensitive channel (see below). chan=0 if no channel applies # obsratio : ratio [sig x BR]_model/[sig x BR]_limit (<1: allowed, >1: excluded) # ncomb : number of Higgs bosons combined in most sensitive channel # Note that the HB channel id number varies depending on the HB version and setting "whichanalyses" # 0 5.3.0beta ||LandH|| # version of HB used to produce these results,the HB setting "whichanalyses" # #CHANNEL info: ranked from highest statistical sensitivity 1 1 296 # channel id number 1 2 1 # HBresult 1 3 0.51704807517546159 # obsratio 1 4 1 # ncombined 1 5 ||(p p)->h1 ->Z Z-> l l l l (low mass) where h1 is SM-like (CMS-PAS-HIG-13-002)|| # text description of channel 2 1 278 # channel id number 2 2 1 # HBresult 2 3 0.43492453977958845 # obsratio 2 4 1 # ncombined 2 5 ||(p p)->h1 ->Z Z-> l l l l where h1 is SM-like (ATLAS-CONF-2013-013)|| # text description of channel 3 1 230 # channel id number 3 2 1 # HBresult 3 3 0.73995763062265141 # obsratio 3 4 1 # ncombined 3 5 ||(p p)->h1 +...->W W +... where h1 is SM-like (CMS-PAS-HIG-13-003)|| # text description of channel # BLOCK HiggsSignalsResults 0 ||2.2.3beta|| # HiggsSignals version - 1 ||LHC7+8|| # experimental data set + 1 ||LHC13|| # experimental data set 2 1 # Chi-squared method ("peak"(1) or "mass"(2)-centered or "both"(3)) 3 2 # Parametrization of Higgs mass uncertainty (1:box, 2:gaussian, 3:box+gaussian) - 4 85 # Number of signal strength peak observables - 5 0 # Number of simplified template cross section (STXS) signal rate observables + 4 56 # Number of signal strength peak observables + 5 24 # Number of simplified template cross section (STXS) signal rate observables 6 20 # Number of LHC Run-1 signal rate observables - 7 5 # Number of Higgs mass observables - 8 110 # Number of observables (total) - 9 67.48060416 # chi^2 (signal strength) from peak observables - 10 0.00000000 # chi^2 (signal strength) from STXS observables + 7 2 # Number of Higgs mass observables + 8 102 # Number of observables (total) + 9 44.57180674 # chi^2 (signal strength) from peak observables + 10 28.96698911 # chi^2 (signal strength) from STXS observables 11 21.77184515 # chi^2 (signal strength) from LHC Run-1 observables - 12 7.17088971 # chi^2 (Higgs mass) from peak observables + 12 0.00865033 # chi^2 (Higgs mass) from peak observables 13 0.00000000 # chi^2 (Higgs mass) from STXS observables 14 0.03179993 # chi^2 (Higgs mass) from LHC Run-1 observables - 15 89.25244931 # chi^2 (signal strength) (total) - 16 7.20268964 # chi^2 (Higgs mass) (total) - 17 96.45513895 # chi^2 (total) - 18 0.86185478 # Probability for peak observables + 15 95.31064100 # chi^2 (signal strength) (total) + 16 0.04045026 # chi^2 (Higgs mass) (total) + 17 95.35109126 # chi^2 (total) + 18 0.88428140 # Probability for peak observables 19 0.41089986 # Probability for LHC-Run1 observables - 20 1.00000000 # Probability for STXS observables - 21 0.81810798 # Probability (total chi^2, total number observables) + 20 0.22137570 # Probability for STXS observables + 21 0.66606708 # Probability (total chi^2, total number observables) Index: trunk/HiggsSignals-2/STXS.f90 =================================================================== --- trunk/HiggsSignals-2/STXS.f90 (revision 599) +++ trunk/HiggsSignals-2/STXS.f90 (revision 600) @@ -1,984 +1,991 @@ module STXS ! Still to do: ! 1: Read in correlation matrix ! 2: Write chi^2 test ! use numerics ! use combinatorics use usefulbits_hs implicit none ! integer :: i,j,k ! double precision,parameter :: pi=3.14159265358979323846264338328D0 ! integer, allocatable :: peakindices_best(:,:) type STXS_observable integer :: id character(LEN=100) :: label ! Reference character(LEN=100) :: desc ! Description character(LEN=3) :: expt ! Experiment character(LEN=10) :: collider character(LEN=10) :: collaboration double precision :: lumi,dlumi,energy character(LEN=100) :: assignmentgroup integer :: rate_SM_normalized integer :: mhchisq double precision :: massobs, dmassobs ! This one enters the chi^2 for the mass! double precision :: mass, dmass ! This one is the mass position for the measurement and the "experimentally allowed assignment range" double precision :: eff_ref_mass ! This is the mass for which the signal efficiencies are given. double precision, allocatable :: model_rate_per_Higgs(:,:) double precision, allocatable :: inclusive_SM_rate(:) integer :: Nc double precision :: model_total_rate double precision :: rate, rate_up, rate_low, drate_up, drate_low double precision :: SMrate, SMrate_up, SMrate_low, dSMrate_up, dSMrate_low ! SM rate used/quoted by the experiment ! At the moment, interpret STXS observables as "pure" channels (production, or decay rate) character(LEN=5), allocatable :: channel_id_str(:) ! Channels array as string, dim(Nc) ! integer, allocatable :: channel_id(:) integer, allocatable :: channel_p_id(:) ! Production channels array, dim(Nc) integer, allocatable :: channel_d_id(:) ! Decay channels array, dim(Nc) double precision, allocatable :: channel_efficiency(:) ! SM signal efficiency of inclusive rates (analysis-specific) double precision, allocatable :: relative_efficiency(:,:) ! Model signal efficiency relative to SM per Higgs double precision :: chisq ! character(LEN=10),allocatable :: channel_description(:,:) ! TODO: How do we deal with ratio of BRs? end type type(STXS_observable), allocatable :: STXSlist(:) type(correlation_info), allocatable :: STXScorrlist(:) contains !------------------------------------------------------------------------------------ subroutine load_STXS(dataset) !------------------------------------------------------------------------------------ use store_pathname_HS use usefulbits, only: file_id_common2, file_id_common3, np, Hneut use datatables, only : read_in_mass_resolution_and_assignment_group ! implicit none character(LEN=*), intent(in) :: dataset character(LEN=100) :: datafile(500) character(LEN=pathname_length+150) :: fullfilename integer, allocatable :: skip(:) integer :: i, n, n_datafiles,n_correlations, n_correlations_tmp, ios, k, m, int1, int2 double precision :: db1 character(LEN=200) :: comment character(LEN=1) :: firstchar character(LEN=100) :: line integer :: id, posperiod - call system('basename -a `ls -1 -p '//trim(adjustl(pathname_HS))// & -& 'Expt_tables/'//trim(adjustl(dataset))//'/*.stxs 2>/dev/null` > STXS_analyses.txt 2>/dev/null') +! call system('basename -a `ls -1 -p '//trim(adjustl(pathname_HS))// & +! & 'Expt_tables/'//trim(adjustl(dataset))//'/*.stxs 2>/dev/null` > STXS_analyses.txt 2>/dev/null') + + call system('ls -1 -p '//trim(adjustl(pathname_HS))// & +& 'Expt_tables/'//trim(adjustl(dataset))//'/*.stxs | xargs -L 1 basename > STXS_analyses.txt') open(file_id_common3, file="STXS_analyses.txt",form='formatted') print *, "Reading in STXS measurements from analysis-set "//& trim(adjustl(dataset))//":" n = 0 n_datafiles = 0 do n = n+1 read(file_id_common3,'(A)', iostat=ios) datafile(n) if(ios.ne.0) exit write(*,'(I4,2X,A)') n, datafile(n) enddo n_datafiles = n - 1 close(file_id_common3) allocate(STXSlist(n_datafiles),skip(n_datafiles)) do n=1,n_datafiles skip(n)=11 open(file_id_common3, file=trim(adjustl(pathname_HS)) //'Expt_tables/'// & & trim(adjustl(dataset))//'/' // datafile(n)) do read(file_id_common3,'(A)') comment comment = trim(adjustl(comment)) write(firstchar,'(A1)') comment if(firstchar.ne.'#') then exit else skip(n)=skip(n)+1 endif enddo backspace(file_id_common3) read(file_id_common3,*) STXSlist(n)%id read(file_id_common3,'(A)') STXSlist(n)%label read(file_id_common3,*) STXSlist(n)%collider,STXSlist(n)%collaboration, & & STXSlist(n)%expt read(file_id_common3,'(A)') STXSlist(n)%desc read(file_id_common3,*) STXSlist(n)%energy, STXSlist(n)%lumi, STXSlist(n)%dlumi read(file_id_common3,*) STXSlist(n)%mhchisq, STXSlist(n)%rate_SM_normalized if(STXSlist(n)%mhchisq == 1) then read(file_id_common3,*) STXSlist(n)%massobs, STXSlist(n)%dmassobs else read(file_id_common3,*) STXSlist(n)%massobs = 0.0D0 STXSlist(n)%dmassobs = 0.0D0 endif !--CHECK FOR ASSIGNMENT GROUP AS SECOND COLUMN: read(file_id_common3,*) STXSlist(n)%mass read(file_id_common3,'(A)') line call read_in_mass_resolution_and_assignment_group(line, STXSlist(n)%dmass,& & STXSlist(n)%assignmentgroup) read(file_id_common3,*) STXSlist(n)%Nc, STXSlist(n)%eff_ref_mass allocate(STXSlist(n)%channel_id_str(STXSlist(n)%Nc)) allocate(STXSlist(n)%channel_p_id(STXSlist(n)%Nc)) allocate(STXSlist(n)%channel_d_id(STXSlist(n)%Nc)) read(file_id_common3,*) (STXSlist(n)%channel_id_str(i),i=1,STXSlist(n)%Nc) do i=1,STXSlist(n)%Nc posperiod = index(STXSlist(n)%channel_id_str(i),'.') if(posperiod.eq.0) then if(len(trim(adjustl(STXSlist(n)%channel_id_str(i)))).eq.2) then read(STXSlist(n)%channel_id_str(i),*) id STXSlist(n)%channel_p_id(i) = int((id-modulo(id,10))/dble(10)) STXSlist(n)%channel_d_id(i) = modulo(id,10) else write(*,*) " For observable ID = ",STXSlist(n)%id stop " Error: Cannot handle channel IDs!" endif else read(STXSlist(n)%channel_id_str(i)(:posperiod-1),*) STXSlist(n)%channel_p_id(i) read(STXSlist(n)%channel_id_str(i)(posperiod+1:),*) STXSlist(n)%channel_d_id(i) endif enddo ! write(*,*) "Production channels = ",STXSlist(n)%channel_p_id ! write(*,*) "Decay channels = ",STXSlist(n)%channel_d_id ! allocate(STXSlist(n)%channel_id(STXSlist(n)%Nc)) ! read(file_id_common3,*) (STXSlist(n)%channel_id(i),i=1,STXSlist(n)%Nc) allocate(STXSlist(n)%channel_efficiency(STXSlist(n)%Nc)) if(STXSlist(n)%eff_ref_mass.ge.0D0) then read(file_id_common3,*) (STXSlist(n)%channel_efficiency(i),i=1,STXSlist(n)%Nc) else do i=1,STXSlist(n)%Nc STXSlist(n)%channel_efficiency(i)=1.0D0 enddo read(file_id_common3,*) endif ! read(file_id_common3,*) STXSlist(n)%channel_id ! read(file_id_common3,*) STXSlist(n)%relative_efficiency read(file_id_common3,*) STXSlist(n)%rate_low, STXSlist(n)%rate, STXSlist(n)%rate_up if(STXSlist(n)%rate_SM_normalized.eq.1) then read(file_id_common3,*) comment STXSlist(n)%SMrate_low = 0.0D0 STXSlist(n)%SMrate = 0.0D0 STXSlist(n)%SMrate_up = 0.0D0 else read(file_id_common3,*) STXSlist(n)%SMrate_low, STXSlist(n)%SMrate, STXSlist(n)%SMrate_up endif STXSlist(n)%drate_low = STXSlist(n)%rate - STXSlist(n)%rate_low STXSlist(n)%drate_up = STXSlist(n)%rate_up - STXSlist(n)%rate STXSlist(n)%dSMrate_low = STXSlist(n)%SMrate - STXSlist(n)%SMrate_low STXSlist(n)%dSMrate_up = STXSlist(n)%SMrate_up - STXSlist(n)%SMrate close(file_id_common3) allocate(STXSlist(n)%relative_efficiency(np(Hneut),STXSlist(n)%Nc)) do k=1, np(Hneut) STXSlist(n)%relative_efficiency(k,:)=1.0D0 enddo enddo close(file_id_common3) !NEW: - call system('basename -a `ls -1 -p '//trim(adjustl(pathname_HS))// & -& 'Expt_tables/'//trim(adjustl(dataset))//'/*.stxscorr 2>/dev/null` > STXS_correlations.txt 2>/dev/null') + + call system('ls -1 -p '//trim(adjustl(pathname_HS))// & +& 'Expt_tables/'//trim(adjustl(dataset))//'/*.stxscorr | xargs -L 1 basename > STXS_correlations.txt') + +! call system('basename -a `ls -1 -p '//trim(adjustl(pathname_HS))// & +! & 'Expt_tables/'//trim(adjustl(dataset))//'/*.stxscorr 2>/dev/null` > STXS_correlations.txt 2>/dev/null') call system('rm -rf STXS_ncorrelations.txt') open(file_id_common3, file="STXS_correlations.txt",form='formatted') print *, "Reading in correlations from the following datafiles in analysis-set "// & trim(adjustl(Exptdir))//":" n = 0 n_datafiles = 0 n_correlations = 0 do n = n+1 read(file_id_common3,'(A)', iostat=ios) datafile(n) if(ios.ne.0) exit fullfilename=trim(adjustl(pathname_HS))//'Expt_tables/'//trim(adjustl(dataset))//'/'& & //trim(datafile(n)) call system('cat '//trim(adjustl(fullfilename))//' | wc -l > STXS_ncorrelations.txt') open(file_id_common2,file="STXS_ncorrelations.txt",form='formatted') read(file_id_common2,'(I10)') n_correlations_tmp close(file_id_common2) write(*,'(2I4,2X,A)') n, n_correlations_tmp, datafile(n) n_correlations = n_correlations + n_correlations_tmp enddo n_datafiles = n - 1 close(file_id_common3) allocate(STXScorrlist(n_correlations)) m=0 do n=1,n_datafiles fullfilename=trim(adjustl(pathname_HS))//'Expt_tables/'//trim(adjustl(dataset))//'/'& & //trim(datafile(n)) open(file_id_common3,file=fullfilename) do m= m+1 read(file_id_common3,*,iostat=ios) int1, int2, db1 ! write(*,*) m, int1, int2, db1, ios if(ios.ne.0) exit STXScorrlist(m)%obsID1 = int1 STXScorrlist(m)%obsID2 = int2 STXScorrlist(m)%corr = db1 enddo m=m-1 close(file_id_common3) enddo end subroutine load_STXS !------------------------------------------------------------------------------------ subroutine assign_modelefficiencies_to_STXS(obsID, Nc, relative_efficiency) !------------------------------------------------------------------------------------ use usefulbits, only : np, Hneut implicit none integer, intent(in) :: obsID integer, intent(in) :: Nc double precision, dimension(np(Hneut),Nc), intent(in) :: relative_efficiency integer :: i logical :: foundid = .False. do i=lbound(STXSlist, dim=1),ubound(STXSlist, dim=1) if(STXSlist(i)%id.eq.obsID) then if(Nc.ne.STXSlist(i)%Nc) then stop 'Error: Number of channels does not match!' else STXSlist(i)%relative_efficiency = relative_efficiency foundid = .True. endif endif enddo if(.not.foundid) write(*,*) "WARNING in assign_modelefficiencies_to_STXS: ",& & "observable ID ",obsID," not known!" end subroutine assign_modelefficiencies_to_STXS !------------------------------------------------------------------------------------ subroutine get_chisq_from_STXS(chisq_tot, pval) !------------------------------------------------------------------------------------ use usefulbits, only : vsmall use usefulbits_hs, only : Nparam use numerics, only : invmatrix, matmult, gammp implicit none double precision, intent(out) :: chisq_tot, pval integer :: i,j,m,N double precision :: cov logical :: correlationfound, somecorrelationsmissing double precision, allocatable :: covmat(:,:),vmat(:,:),invcovmat(:,:) double precision, allocatable :: v(:), v2(:) N=size(STXSlist) allocate(covmat(N,N),invcovmat(N,N)) allocate(v(N),v2(N)) allocate(vmat(N,1)) somecorrelationsmissing = .False. do i=1,N do j=1,N correlationfound=.False. do m=lbound(STXScorrlist,dim=1), ubound(STXScorrlist,dim=1) if((STXScorrlist(m)%obsID1.eq.STXSlist(i)%id.and.STXScorrlist(m)%obsID2.eq.STXSlist(j)%id)& &.or.(STXScorrlist(m)%obsID2.eq.STXSlist(i)%id.and.STXScorrlist(m)%obsID1.eq.STXSlist(j)%id)) then covmat(i,j) = STXScorrlist(m)%corr*get_drate(i)*get_drate(j) correlationfound=.True. endif enddo if(.not.correlationfound) then ! Use a unit-matrix for the correlations here. ! if(.not.somecorrelationsmissing) then ! write(*,*) "Warning: Correlation matrix element not found for observable ids: ",STXSlist(i)%id, STXSlist(j)%id ! write(*,*) " Suppressing future warnings about missing correlation matrix elements..." ! endif covmat(i,j) = 0.0D0 somecorrelationsmissing = .True. if(STXSlist(i)%id.eq.STXSlist(j)%id) then covmat(i,j) = get_drate(i)*get_drate(j) endif endif enddo v(i) = STXSlist(i)%rate - STXSlist(i)%model_total_rate vmat(i,1) = v(i) enddo ! if(somecorrelationsmissing) then ! write(*,*) "Warning: Some correlation matrix elements were not found." ! endif call invmatrix(covmat,invcovmat) call matmult(invcovmat,vmat,v2,N,1) chisq_tot= 0.0D0 do i=1,N STXSlist(i)%chisq = v(i)*v2(i) chisq_tot = chisq_tot + STXSlist(i)%chisq enddo pval = 1.0D0 if(chisq_tot.gt.vsmall.and.(N-Nparam).gt.0) then pval = 1 - gammp(dble(N-Nparam)/2,chisq_tot/2) endif deallocate(covmat,invcovmat,v,v2,vmat) end subroutine get_chisq_from_STXS !------------------------------------------------------------------------------------ subroutine get_number_of_STXS_observables(Nobs_rates, Nobs_mh) integer, intent(out) :: Nobs_rates, Nobs_mh Nobs_rates=size(STXSlist) Nobs_mh = 0 end subroutine get_number_of_STXS_observables !------------------------------------------------------------------------------------ function get_drate(i) !------------------------------------------------------------------------------------ implicit none integer :: i double precision get_drate if(STXSlist(i)%model_total_rate.le.STXSlist(i)%rate) then get_drate = STXSlist(i)%drate_low else get_drate = STXSlist(i)%drate_up endif end function get_drate !------------------------------------------------------------------------------------ subroutine calculate_model_predictions_for_STXS() !------------------------------------------------------------------------------------ use usefulbits, only : theo use theo_manip, only : HB5_complete_theo integer :: i call HB5_complete_theo do i=lbound(STXSlist,dim=1), ubound(STXSlist,dim=1) call evaluate_model_for_STXS(STXSlist(i),theo(1)) enddo end subroutine calculate_model_predictions_for_STXS !------------------------------------------------------------------------------------ subroutine evaluate_model_for_STXS(STXSobs, t) !------------------------------------------------------------------------------------ use usefulbits, only : theo, div, small, np, Hneut, dataset, vsmall use usefulbits_HS, only : normalize_rates_to_reference_position, & & normalize_rates_to_reference_position_outside_dmtheo, & & assignmentrange_STXS ! use_SMrate_at_reference_position_for_STXS, use theory_XS_SM_functions use theory_BRfunctions use theo_manip, only : HB5_complete_theo implicit none type(STXS_observable), intent(inout) :: STXSobs type(dataset), intent(in) :: t double precision :: norm_rate, SMrate, SMrate_refmass, refmass, BR_SMref integer :: i, j, id, p, d STXSobs%model_total_rate = 0.0D0 if(.not.allocated(theo))then stop 'subroutine HiggsSignals_initialize must be called first' endif if(.not.allocated(STXSobs%model_rate_per_Higgs)) then allocate(STXSobs%model_rate_per_Higgs(np(Hneut),STXSobs%Nc)) endif if(.not.allocated(STXSobs%inclusive_SM_rate)) then ! allocate(STXSobs%inclusive_SM_rate(np(Hneut),STXSobs%Nc)) allocate(STXSobs%inclusive_SM_rate(STXSobs%Nc)) endif ! write(*,*) 'DEBUG HS: id = ', STXSobs%id ! write(*,*) 'DEBUG HS, channel = ',STXSobs%channel_id refmass = STXSobs%mass do i=1,STXSobs%Nc ! id = STXSobs%channel_id(i) ! p = int((id-modulo(id,10))/dble(10)) ! d = modulo(id,10) p = STXSobs%channel_p_id(i) d = STXSobs%channel_d_id(i) do j=1, np(Hneut) ! write(*,*) 'DEBUG HS, m = ', t%particle(Hneut)%M(j) !--Do the production rate for the relevant experiment and cms-energy if(STXSobs%collider.eq.'LHC') then if(abs(STXSobs%energy-7.0D0).le.small) then if(p.eq.1) then norm_rate=t%lhc7%XS_hj_ratio(j) SMrate=t%lhc7%XS_H_SM(j) SMrate_refmass=XS_lhc7_gg_H_SM(refmass)+XS_lhc7_bb_H_SM(refmass) ! STXSobs%channel_description(i,1)='singleH' else if(p.eq.2) then norm_rate=t%lhc7%XS_vbf_ratio(j) SMrate=t%lhc7%XS_vbf_SM(j) SMrate_refmass=XS_lhc7_vbf_SM(refmass) ! STXSobs%channel_description(i,1)='VBF' else if(p.eq.3) then norm_rate=t%lhc7%XS_hjW_ratio(j) SMrate=t%lhc7%XS_HW_SM(j) SMrate_refmass=XS_lhc7_HW_SM(refmass) ! STXSobs%channel_description(i,1)='HW' else if(p.eq.4) then norm_rate=t%lhc7%XS_hjZ_ratio(j) SMrate=t%lhc7%XS_HZ_SM(j) SMrate_refmass=ZH_cpmix_nnlo_ggqqbb(refmass,'LHC7 ',1.0D0,1.0D0,1.0D0,0.0D0,0.0D0,.True.) ! STXSobs%channel_description(i,1)='HZ' else if(p.eq.5) then norm_rate=t%lhc7%XS_tthj_ratio(j) SMrate=t%lhc7%XS_ttH_SM(j) SMrate_refmass=XS_lhc7_ttH_SM(refmass) ! STXSobs%channel_description(i,1)='ttH' else if(p.eq.6) then norm_rate=t%lhc7%XS_gg_hj_ratio(j) SMrate=t%lhc7%XS_gg_H_SM(j) SMrate_refmass=XS_lhc7_gg_H_SM(refmass) ! mutab%channel_description(i,1)='ggH' else if(p.eq.7) then norm_rate=t%lhc7%XS_bb_hj_ratio(j) SMrate=t%lhc7%XS_bb_H_SM(j) SMrate_refmass=XS_lhc7_bb_H_SM(refmass) ! mutab%channel_description(i,1)='bbH' else if(p.eq.8) then norm_rate=t%lhc7%XS_thj_tchan_ratio(j) SMrate=t%lhc7%XS_tH_tchan_SM(j) SMrate_refmass=XS_lhc7_tH_tchan_SM(refmass) ! mutab%channel_description(i,1)='tH (t-channel)' else if(p.eq.9) then norm_rate=t%lhc7%XS_thj_schan_ratio(j) SMrate=t%lhc7%XS_tH_schan_SM(j) SMrate_refmass=XS_lhc7_tH_schan_SM(refmass) ! mutab%channel_description(i,1)='tH (s-channel)' else if(p.eq.10) then norm_rate=t%lhc7%XS_qq_hjZ_ratio(j) SMrate=t%lhc7%XS_qq_HZ_SM(j) SMrate_refmass=ZH_cpmix_nnlo_qqbb(refmass,'LHC7 ',1.0D0,1.0D0,1.0D0,0.0D0,0.0D0,.True.) ! mutab%channel_description(i,1)='qq-HZ' else if(p.eq.11) then norm_rate=t%lhc7%XS_gg_hjZ_ratio(j) SMrate=t%lhc7%XS_gg_HZ_SM(j) SMrate_refmass=ZH_cpmix_nnlo_gg(refmass,'LHC7 ',1.0D0,1.0D0,1.0D0,0.0D0,0.0D0,.True.) ! mutab%channel_description(i,1)='gg-HZ' else if(p.eq.0) then norm_rate=1.0D0 SMrate=1.0D0 SMrate_refmass=1.0D0 ! STXSobs%channel_description(i,1)='none' else write(*,*) "WARNING: Unknown production mode id p=",p," for STXS observable id = ",STXSobs%id endif else if(abs(STXSobs%energy-8.0D0).le.small) then if(p.eq.1) then norm_rate=t%lhc8%XS_hj_ratio(j) SMrate=t%lhc8%XS_H_SM(j) SMrate_refmass=XS_lhc8_gg_H_SM(refmass)+XS_lhc8_bb_H_SM(refmass) ! STXSobs%channel_description(i,1)='singleH' else if(p.eq.2) then norm_rate=t%lhc8%XS_vbf_ratio(j) SMrate=t%lhc8%XS_vbf_SM(j) SMrate_refmass=XS_lhc8_vbf_SM(refmass) ! STXSobs%channel_description(i,1)='VBF' else if(p.eq.3) then norm_rate=t%lhc8%XS_hjW_ratio(j) SMrate=t%lhc8%XS_HW_SM(j) SMrate_refmass=XS_lhc8_HW_SM(refmass) ! STXSobs%channel_description(i,1)='HW' else if(p.eq.4) then norm_rate=t%lhc8%XS_hjZ_ratio(j) SMrate=t%lhc8%XS_HZ_SM(j) SMrate_refmass=ZH_cpmix_nnlo_ggqqbb(refmass,'LHC8 ',1.0D0,1.0D0,1.0D0,0.0D0,0.0D0,.True.) ! STXSobs%channel_description(i,1)='HZ' else if(p.eq.5) then norm_rate=t%lhc8%XS_tthj_ratio(j) SMrate=t%lhc8%XS_ttH_SM(j) SMrate_refmass=XS_lhc8_ttH_SM(refmass) ! STXSobs%channel_description(i,1)='ttH' else if(p.eq.6) then norm_rate=t%lhc8%XS_gg_hj_ratio(j) SMrate=t%lhc8%XS_gg_H_SM(j) SMrate_refmass=XS_lhc8_gg_H_SM(refmass) ! mutab%channel_description(i,1)='ggH' else if(p.eq.7) then norm_rate=t%lhc8%XS_bb_hj_ratio(j) SMrate=t%lhc8%XS_bb_H_SM(j) SMrate_refmass=XS_lhc8_bb_H_SM(refmass) ! mutab%channel_description(i,1)='bbH' else if(p.eq.8) then norm_rate=t%lhc8%XS_thj_tchan_ratio(j) SMrate=t%lhc8%XS_tH_tchan_SM(j) SMrate_refmass=XS_lhc8_tH_tchan_SM(refmass) ! mutab%channel_description(i,1)='tH (t-channel)' else if(p.eq.9) then norm_rate=t%lhc8%XS_thj_schan_ratio(j) SMrate=t%lhc8%XS_tH_schan_SM(j) SMrate_refmass=XS_lhc8_tH_schan_SM(refmass) ! mutab%channel_description(i,1)='tH (s-channel)' else if(p.eq.10) then norm_rate=t%lhc8%XS_qq_hjZ_ratio(j) SMrate=t%lhc8%XS_qq_HZ_SM(j) SMrate_refmass=ZH_cpmix_nnlo_qqbb(refmass,'LHC8 ',1.0D0,1.0D0,1.0D0,0.0D0,0.0D0,.True.) ! mutab%channel_description(i,1)='qq-HZ' else if(p.eq.11) then norm_rate=t%lhc8%XS_gg_hjZ_ratio(j) SMrate=t%lhc8%XS_gg_HZ_SM(j) SMrate_refmass=ZH_cpmix_nnlo_gg(refmass,'LHC8 ',1.0D0,1.0D0,1.0D0,0.0D0,0.0D0,.True.) else if(p.eq.0) then norm_rate=1.0D0 SMrate=1.0D0 SMrate_refmass=1.0D0 ! STXSobs%channel_description(i,1)='none' else write(*,*) "WARNING: Unknown production mode id p=",p," for STXS observable id = ",STXSobs%id endif else if(abs(STXSobs%energy-13.0D0).le.small) then if(p.eq.1) then norm_rate=t%lhc13%XS_hj_ratio(j) SMrate=t%lhc13%XS_H_SM(j) SMrate_refmass=XS_lhc13_gg_H_SM(refmass)+XS_lhc13_bb_H_SM(refmass) ! STXSobs%channel_description(i,1)='singleH' else if(p.eq.2) then norm_rate=t%lhc13%XS_vbf_ratio(j) SMrate=t%lhc13%XS_vbf_SM(j) SMrate_refmass=XS_lhc13_vbf_SM(refmass) ! STXSobs%channel_description(i,1)='VBF' else if(p.eq.3) then norm_rate=t%lhc13%XS_hjW_ratio(j) SMrate=t%lhc13%XS_HW_SM(j) SMrate_refmass=XS_lhc13_HW_SM(refmass) ! STXSobs%channel_description(i,1)='HW' else if(p.eq.4) then norm_rate=t%lhc13%XS_hjZ_ratio(j) SMrate=t%lhc13%XS_HZ_SM(j) SMrate_refmass=ZH_cpmix_nnlo_ggqqbb(refmass,'LHC13',1.0D0,1.0D0,1.0D0,0.0D0,0.0D0,.True.) ! STXSobs%channel_description(i,1)='HZ' else if(p.eq.5) then norm_rate=t%lhc13%XS_tthj_ratio(j) SMrate=t%lhc13%XS_ttH_SM(j) SMrate_refmass=XS_lhc13_ttH_SM(refmass) ! STXSobs%channel_description(i,1)='ttH' else if(p.eq.6) then norm_rate=t%lhc13%XS_gg_hj_ratio(j) SMrate=t%lhc13%XS_gg_H_SM(j) SMrate_refmass=XS_lhc13_gg_H_SM(refmass) ! mutab%channel_description(i,1)='ggH' else if(p.eq.7) then norm_rate=t%lhc13%XS_bb_hj_ratio(j) SMrate=t%lhc13%XS_bb_H_SM(j) SMrate_refmass=XS_lhc13_bb_H_SM(refmass) ! mutab%channel_description(i,1)='bbH' else if(p.eq.8) then norm_rate=t%lhc13%XS_thj_tchan_ratio(j) SMrate=t%lhc13%XS_tH_tchan_SM(j) SMrate_refmass=XS_lhc13_tH_tchan_SM(refmass) ! mutab%channel_description(i,1)='tH (t-channel)' else if(p.eq.9) then norm_rate=t%lhc13%XS_thj_schan_ratio(j) SMrate=t%lhc13%XS_tH_schan_SM(j) SMrate_refmass=XS_lhc13_tH_schan_SM(refmass) ! mutab%channel_description(i,1)='tH (s-channel)' else if(p.eq.10) then norm_rate=t%lhc13%XS_qq_hjZ_ratio(j) SMrate=t%lhc13%XS_qq_HZ_SM(j) SMrate_refmass=ZH_cpmix_nnlo_qqbb(refmass,'LHC13',1.0D0,1.0D0,1.0D0,0.0D0,0.0D0,.True.) ! mutab%channel_description(i,1)='qq-HZ' else if(p.eq.11) then norm_rate=t%lhc13%XS_gg_hjZ_ratio(j) SMrate=t%lhc13%XS_gg_HZ_SM(j) SMrate_refmass=ZH_cpmix_nnlo_gg(refmass,'LHC13',1.0D0,1.0D0,1.0D0,0.0D0,0.0D0,.True.) else if(p.eq.0) then norm_rate=1.0D0 SMrate=1.0D0 SMrate_refmass=1.0D0 ! STXSobs%channel_description(i,1)='none' else write(*,*) "WARNING: Unknown production mode id p=",p," for STXS observable id = ",STXSobs%id endif endif else if(STXSobs%collider.eq.'TEV') then if(p.eq.1) then norm_rate=t%tev%XS_hj_ratio(j) SMrate=t%tev%XS_H_SM(j) SMrate_refmass=XS_tev_gg_H_SM(refmass)+XS_tev_bb_H_SM(refmass) ! STXSobs%channel_description(i,1)='singleH' else if(p.eq.2) then norm_rate=t%tev%XS_vbf_ratio(j) SMrate=t%tev%XS_vbf_SM(j) SMrate_refmass=XS_tev_vbf_SM(refmass) ! STXSobs%channel_description(i,1)='VBF' else if(p.eq.3) then norm_rate=t%tev%XS_hjW_ratio(j) SMrate=t%tev%XS_HW_SM(j) SMrate_refmass=XS_tev_HW_SM(refmass) ! STXSobs%channel_description(i,1)='HW' else if(p.eq.4) then norm_rate=t%tev%XS_hjZ_ratio(j) SMrate=t%tev%XS_HZ_SM(j) SMrate_refmass=ZH_cpmix_nnlo_ggqqbb(refmass,'TEV ',1.0D0,1.0D0,1.0D0,0.0D0,0.0D0,.True.) ! STXSobs%channel_description(i,1)='HZ' else if(p.eq.5) then norm_rate=t%tev%XS_tthj_ratio(j) SMrate=t%tev%XS_ttH_SM(j) SMrate_refmass=XS_tev_ttH_SM(refmass) ! STXSobs%channel_description(i,1)='ttH' else if(p.eq.10) then norm_rate=t%tev%XS_qq_hjZ_ratio(j) SMrate=t%tev%XS_qq_HZ_SM(j) SMrate_refmass=ZH_cpmix_nnlo_qqbb(refmass,'TEV ',1.0D0,1.0D0,1.0D0,0.0D0,0.0D0,.True.) ! mutab%channel_description(i,1)='qq-HZ' else if(p.eq.11) then norm_rate=t%tev%XS_gg_hjZ_ratio(j) SMrate=t%tev%XS_gg_HZ_SM(j) SMrate_refmass=ZH_cpmix_nnlo_gg(refmass,'TEV ',1.0D0,1.0D0,1.0D0,0.0D0,0.0D0,.True.) else if(p.eq.0) then norm_rate=1.0D0 SMrate=1.0D0 SMrate_refmass=1.0D0 ! STXSobs%channel_description(i,1)='none' else write(*,*) "WARNING: Unknown production mode id p=",p," for STXS observable id = ",STXSobs%id endif else if(STXSobs%collider.eq.'ILC') then !--n.B.: As a first attempt, we use the LHC8 normalized cross sections for ZH, VBF, ttH. ! In order to do this properly, a separate input for the ILC cross sections ! has to be provided! It works only for single production mode observables (no ! correct weighting of channels included!)Then, at least in the effective coupling ! approximation, there is no difference to a full implementation. ! The theoretical uncertainty of the ILC production modes will are defined in ! usefulbits_HS.f90. if(p.eq.1.or.p.eq.2) then write(*,*) 'Warning: Unknown ILC production mode (',p,') in table ',STXSobs%id norm_rate=0.0D0 SMrate=1.0D0 SMrate_refmass=1.0D0 ! STXSobs%channel_description(i,1)='unknown' else if(p.eq.3) then norm_rate=t%lhc8%XS_hjW_ratio(j) SMrate=t%lhc8%XS_HW_SM(j) SMrate_refmass=XS_lhc8_HW_SM(refmass) ! STXSobs%channel_description(i,1)='WBF' else if(p.eq.4) then norm_rate=t%lhc8%XS_hjZ_ratio(j) SMrate=t%lhc8%XS_HZ_SM(j) SMrate_refmass=XS_lhc8_HZ_SM(refmass) ! STXSobs%channel_description(i,1)='HZ' else if(p.eq.5) then norm_rate=t%lhc8%XS_tthj_ratio(j) SMrate=t%lhc8%XS_ttH_SM(j) SMrate_refmass=XS_lhc8_ttH_SM(refmass) ! STXSobs%channel_description(i,1)='ttH' else if(p.eq.0) then norm_rate=1.0D0 SMrate=1.0D0 SMrate_refmass=1.0D0 ! STXSobs%channel_description(i,1)='none' else write(*,*) "WARNING: Unknown production mode id p=",p," for STXS observable id = ",STXSobs%id endif endif !--Multiply now by the decay rate if(d.eq.1) then norm_rate=norm_rate*div(t%BR_hjgaga(j),t%BR_Hgaga_SM(j),0.0D0,1.0D0) SMrate=SMrate*t%BR_Hgaga_SM(j) SMrate_refmass = SMrate_refmass*BRSM_Hgaga(refmass) ! STXSobs%channel_description(i,2)='gammagamma' else if(d.eq.2) then norm_rate=norm_rate*div(t%BR_hjWW(j),t%BR_HWW_SM(j),0.0D0,1.0D0) SMrate=SMrate*t%BR_HWW_SM(j) SMrate_refmass = SMrate_refmass*BRSM_HWW(refmass) ! STXSobs%channel_description(i,2)='WW' else if(d.eq.3) then norm_rate=norm_rate*div(t%BR_hjZZ(j),t%BR_HZZ_SM(j),0.0D0,1.0D0) SMrate=SMrate*t%BR_HZZ_SM(j) SMrate_refmass = SMrate_refmass*BRSM_HZZ(refmass) ! STXSobs%channel_description(i,2)='ZZ' else if(d.eq.4) then norm_rate=norm_rate*div(t%BR_hjtautau(j),t%BR_Htautau_SM(j),0.0D0,1.0D0) SMrate=SMrate*t%BR_Htautau_SM(j) SMrate_refmass = SMrate_refmass*BRSM_Htautau(refmass) ! STXSobs%channel_description(i,2)='tautau' else if(d.eq.5) then norm_rate=norm_rate*div(t%BR_hjbb(j),t%BR_Hbb_SM(j),0.0D0,1.0D0) SMrate=SMrate*t%BR_Hbb_SM(j) SMrate_refmass = SMrate_refmass*BRSM_Hbb(refmass) ! STXSobs%channel_description(i,2)='bb' else if(d.eq.6) then norm_rate=norm_rate*div(t%BR_hjZga(j),t%BR_HZga_SM(j),0.0D0,1.0D0) SMrate=SMrate*t%BR_HZga_SM(j) SMrate_refmass = SMrate_refmass*BRSM_HZga(refmass) ! STXSobs%channel_description(i,2)='Zgamma' else if(d.eq.7) then norm_rate=norm_rate*div(t%BR_hjcc(j),t%BR_Hcc_SM(j),0.0D0,1.0D0) SMrate=SMrate*t%BR_Hcc_SM(j) SMrate_refmass = SMrate_refmass*BRSM_Hcc(refmass) ! STXSobs%channel_description(i,2)='cc' else if(d.eq.8) then norm_rate=norm_rate*div(t%BR_hjmumu(j),t%BR_Hmumu_SM(j),0.0D0,1.0D0) SMrate=SMrate*t%BR_Hmumu_SM(j) SMrate_refmass = SMrate_refmass*BRSM_Hmumu(refmass) ! STXSobs%channel_description(i,2)='mumu' else if(d.eq.9) then norm_rate=norm_rate*div(t%BR_hjgg(j),t%BR_Hgg_SM(j),0.0D0,1.0D0) SMrate=SMrate*t%BR_Hgg_SM(j) SMrate_refmass = SMrate_refmass*BRSM_Hgg(refmass) ! STXSobs%channel_description(i,2)='gg' else if(d.eq.10) then norm_rate=norm_rate*div(t%BR_hjss(j),t%BR_Hss_SM(j),0.0D0,1.0D0) SMrate=SMrate*t%BR_Hss_SM(j) SMrate_refmass = SMrate_refmass*BRSM_Hss(refmass) ! mutab%channel_description(i,2)='ss' else if(d.eq.11) then norm_rate=norm_rate*div(t%BR_hjtt(j),t%BR_Htt_SM(j),0.0D0,1.0D0) SMrate=SMrate*t%BR_Htt_SM(j) SMrate_refmass = SMrate_refmass*BRSM_Htoptop(refmass) ! mutab%channel_description(i,2)='tt' else if(d.eq.0) then norm_rate=norm_rate*1.0D0 SMrate=SMrate*1.0D0 SMrate_refmass = SMrate_refmass*1.0D0 ! STXSobs%channel_description(i,2)='none' endif !------------------------- ! NEW FEATURE (since HB-5.2): Enable to set channelrates directly. if(p.ne.0.and.d.ne.0) then select case(d) case(1) BR_SMref = t%BR_Hgaga_SM(j) ! BR_SMref_mpeak = BRSM_Hgaga(refmass) case(2) BR_SMref = t%BR_HWW_SM(j) ! BR_SMref_mpeak = BRSM_HWW(refmass) case(3) BR_SMref = t%BR_HZZ_SM(j) ! BR_SMref_mpeak = BRSM_HZZ(refmass) case(4) BR_SMref = t%BR_Htautau_SM(j) ! BR_SMref_mpeak = BRSM_Htautau(refmass) case(5) BR_SMref = t%BR_Hbb_SM(j) ! BR_SMref_mpeak = BRSM_Hbb(refmass) case(6) BR_SMref = t%BR_HZga_SM(j) ! BR_SMref_mpeak = BRSM_HZga(refmass) case(7) BR_SMref = t%BR_Hcc_SM(j) ! BR_SMref_mpeak = BRSM_Hcc(refmass) case(8) BR_SMref = t%BR_Hmumu_SM(j) ! BR_SMref_mpeak = BRSM_Hmumu(refmass) case(9) BR_SMref = t%BR_Hgg_SM(j) ! BR_SMref_mpeak = BRSM_Hgg(refmass) case(10) BR_SMref = t%BR_Hss_SM(j) ! BR_SMref_mpeak = BRSM_Hgg(refmass) case(11) BR_SMref = t%BR_Htt_SM(j) ! BR_SMref_mpeak = BRSM_Htoptop(refmass) end select if(STXSobs%collider.eq.'LHC') then if(abs(STXSobs%energy-7.0D0).le.small) then if(t%lhc7%channelrates(j,p,d).ge.0.0d0) then norm_rate=div(t%lhc7%channelrates(j,p,d),BR_SMref,0.0D0,1.0D0) endif else if(abs(STXSobs%energy-8.0D0).le.small) then if(t%lhc8%channelrates(j,p,d).ge.0.0d0) then norm_rate=div(t%lhc8%channelrates(j,p,d),BR_SMref,0.0D0,1.0D0) endif else if(abs(STXSobs%energy-13.0D0).le.small) then if(t%lhc13%channelrates(j,p,d).ge.0.0d0) then norm_rate=div(t%lhc13%channelrates(j,p,d),BR_SMref,0.0D0,1.0D0) endif endif else if(STXSobs%collider.eq.'TEV') then if(t%tev%channelrates(j,p,d).ge.0.0d0) then norm_rate=div(t%tev%channelrates(j,p,d),BR_SMref,0.0D0,1.0D0) endif endif endif !------------------------- if(abs(t%particle(Hneut)%M(j) - STXSobs%mass).le.(assignmentrange_STXS * & sqrt(t%particle(Hneut)%dM(j)**2.0D0+STXSobs%dmass**2.0D0))) then ! if(STXSobs%rate_SM_normalized.eq.1) then if(normalize_rates_to_reference_position) then STXSobs%model_rate_per_Higgs(j,i)=norm_rate*SMrate/(SMrate_refmass) else STXSobs%model_rate_per_Higgs(j,i)=norm_rate !! OLD WAY endif if(normalize_rates_to_reference_position_outside_dmtheo) then if(abs(STXSobs%mass-t%particle(Hneut)%M(j)).ge.t%particle(Hneut)%dM(j)) then STXSobs%model_rate_per_Higgs(j,i)=norm_rate*SMrate/(SMrate_refmass) endif endif ! else ! if(use_SMrate_at_reference_position_for_STXS) then !--- ! n.B.: Need to use officially quoted SM prediction here, because HB/HS do not contain ! SM predictions for exclusive STXS bins (but only inclusive SM rates). !--- ! STXSobs%model_rate_per_Higgs(j,i)=norm_rate*STXSobs%SMrate ! else ! STXSobs%model_rate=norm_rate*SMrate ! endif ! endif else STXSobs%model_rate_per_Higgs(j,i) = 0.0D0 ! STXSobs%inclusive_SM_rate(j,i) = 0.0D0 endif ! Inclusive SM rate must always be evaluated at the mass position of the measurement! STXSobs%inclusive_SM_rate(i) = SMrate_refmass * STXSobs%channel_efficiency(i) ! write(*,*) "j, i, STXSobs%model_rate_per_Higgs(j,i) = ",j,i, STXSobs%model_rate_per_Higgs(j,i) ! Turn normalized rate into absolute rate (per Higgs per channel) STXSobs%model_rate_per_Higgs(j,i) = STXSobs%model_rate_per_Higgs(j,i) * & & STXSobs%relative_efficiency(j,i) * & & STXSobs%inclusive_SM_rate(i) !& SMrate * STXSobs%channel_efficiency(i) ! write(*,*) "j, i, absolute STXSobs%model_rate_per_Higgs(j,i), STXSobs%inclusive_SM_rate(i) = ",& ! & j, i, STXSobs%model_rate_per_Higgs(j,i), STXSobs%inclusive_SM_rate(j,i) !--- ! Take into account model-dependent signal efficiency (relative to SM). ! These have to be given by the user for each observable using the subroutine ! assign_modelefficiencies_to_STXS: !--- enddo enddo ! write(*,*) "STXSobs%id = " , STXSobs%id ! write(*,*) " model_rate_per_Higgs = ",STXSobs%model_rate_per_Higgs ! write(*,*) " inclusive SM rate = ",STXSobs%inclusive_SM_rate if(sum(STXSobs%inclusive_SM_rate).ge.vsmall) then STXSobs%model_total_rate = sum(STXSobs%model_rate_per_Higgs)/sum(STXSobs%inclusive_SM_rate) ! Mistake: Don't divide by the sum over SM!! else STXSobs%model_total_rate = 0.0D0 endif ! write(*,*) "STXSobs%model_total_rate (SM norm)= ", STXSobs%model_total_rate if(STXSobs%rate_SM_normalized.eq.0) then STXSobs%model_total_rate = STXSobs%model_total_rate * STXSobs%SMrate ! write(*,*) "STXSobs%model_total_rate (absolute)= ", STXSobs%model_total_rate endif ! write(*,*) "#--------------- ", STXSobs%id, " ---------------#" ! do i=1,STXSobs%Nc ! write(*,*) "channel id = ", STXSobs%channel_id(i), " rate = ", & ! & sum(STXSobs%model_rate_per_Higgs(:,i))/sum(STXSobs%inclusive_SM_rate)*STXSobs%SMrate ! enddo ! STXSobs%model_total_rate + STXSobs%relative_efficiency(j) * STXSobs%model_rate_per_Higgs(j) ! write(*,*) "Total rate: ", STXSobs%model_total_rate end subroutine evaluate_model_for_STXS !------------------------------------------------------------------------------------ subroutine print_STXS() !------------------------------------------------------------------------------------ implicit none integer :: i character(LEN=100) :: formatter do i=lbound(STXSlist,dim=1), ubound(STXSlist,dim=1) write(*,*) "#--------------------------------------------------#" write(*,*) "#- STXS observable ",i," -#" write(*,*) "#--------------------------------------------------#" write(*,'(A,I10)') " ID = ", STXSlist(i)%id write(*,'(A,A)') " Label = ", STXSlist(i)%label write(*,'(A,A)') " Description = ", STXSlist(i)%desc write(*,'(A,A)') " Experiment = ", STXSlist(i)%expt write(*,'(A,2F6.2)') " Energy, Luminosity = ", STXSlist(i)%energy, STXSlist(i)%lumi write(*,'(A,F10.5,A,F10.5,A,F10.5)') " Obs Signal rate [pb] = ",& & STXSlist(i)%rate, " + ", STXSlist(i)%drate_up, " - ", STXSlist(i)%drate_low write(*,'(A,F10.5,A,F10.5,A,F10.5)') " SM Signal rate [pb] = ",& & STXSlist(i)%SMrate, " + ", STXSlist(i)%dSMrate_up, " - ", STXSlist(i)%dSMrate_low write(*,'(A,F10.5)') " Pred. Signal rate [pb] = ", STXSlist(i)%model_total_rate write(formatter,*) "(A,",STXSlist(i)%Nc,"I10)" formatter = trim(adjustl(formatter)) write(*,'(A)') " Channels = ", STXSlist(i)%channel_id_str write(formatter,*) "(A,",STXSlist(i)%Nc,"F10.5)" formatter = trim(adjustl(formatter)) write(*,formatter) " Channel efficiency = ", STXSlist(i)%channel_efficiency enddo write(*,*) "#--------------------------------------------------#" end subroutine print_STXS !------------------------------------------------------------------------------------ subroutine print_STXS_to_file !------------------------------------------------------------------------------------ use usefulbits, only : file_id_common3 use usefulbits_hs, only : StrCompress implicit none character(LEN=100) :: formatspec integer :: i formatspec='(I3,7X,I10,1X,F6.2,1X,6F10.6,1X,A3,1X,F6.2,1X,F6.2,1X,A,5X,A)' open(file_id_common3,file="STXS_information.txt") write(file_id_common3,*) "#HiggsSignals-"//trim(adjustl(HSvers))// & & " with experimental dataset '"//trim(adjustl(Exptdir))//"'" write(file_id_common3,*) "#Number STXS-ID mass-pos rate_obs drate_low drate_high ", & & "rate_SM dSMrate_low dSMrate_high collaboration energy luminosity description reference" write(file_id_common3,*) "#" do i=lbound(STXSlist,dim=1),ubound(STXSlist,dim=1) write(file_id_common3,formatspec) i ,STXSlist(i)%id,STXSlist(i)%mass, & & STXSlist(i)%rate, STXSlist(i)%drate_low,STXSlist(i)%drate_up, & & STXSlist(i)%SMrate, STXSlist(i)%dSMrate_low,STXSlist(i)%dSMrate_up, & & STXSlist(i)%collaboration, STXSlist(i)%energy, & & STXSlist(i)%lumi, trim(strcompress(STXSlist(i)%desc)), STXSlist(i)%label enddo close(file_id_common3) end subroutine print_STXS_to_file !------------------------------------------------------------------------------------ subroutine clear_STXS() !------------------------------------------------------------------------------------ implicit none integer :: i do i=lbound(STXSlist,dim=1), ubound(STXSlist,dim=1) deallocate(STXSlist(i)%model_rate_per_Higgs) deallocate(STXSlist(i)%inclusive_SM_rate) deallocate(STXSlist(i)%channel_id_str) deallocate(STXSlist(i)%channel_p_id) deallocate(STXSlist(i)%channel_d_id) deallocate(STXSlist(i)%channel_efficiency) deallocate(STXSlist(i)%relative_efficiency) enddo deallocate(STXSlist) if(allocated(STXScorrlist)) deallocate(STXScorrlist) end subroutine clear_STXS !------------------------------------------------------------------------------------ end module STXS !------------------------------------------------------------------------------------ \ No newline at end of file Index: trunk/HiggsSignals-2/datatables.f90 =================================================================== --- trunk/HiggsSignals-2/datatables.f90 (revision 599) +++ trunk/HiggsSignals-2/datatables.f90 (revision 600) @@ -1,733 +1,738 @@ !-------------------------------------------------------------------- ! This file is part of HiggsSignals (OS and TS 04/03/2013) !-------------------------------------------------------------------- module datatables !-------------------------------------------------------------------- use usefulbits_HS, only : mutable, mupeak, Exptdir, obs, analyses,corrlist !, withcorrexpsyst !x mutables, peaks use usefulbits, only : np implicit none integer :: ntable1,ntable2 contains !-------------------------------------------------------------------- subroutine initializemutable_blank(table) !-------------------------------------------------------------------- type(mutable) :: table table%id = -1 table%nx = -1 table%particle_x = -1 table%label = '' table%desc = '' table%expt = '' table%lumi = -1.0D0 table%dlumi = -1.0D0 table%energy = -1.0D0 table%Nc = -1 table%xmax = -1.0D0 table%xmin = -1.0D0 table%sep = -1.0D0 table%deltam = -1.0D0 table%deltax = -1.0D0 table%mhchisq = 0 end subroutine initializemutable_blank !-------------------------------------------------------------------- subroutine initialize_observables !-------------------------------------------------------------------- use store_pathname_HS use usefulbits, only: np,Hneut,Hplus,file_id_common3,file_id_common2 ! implicit none character(LEN=150) :: filename character(LEN=100) :: datafile(500) integer n_datafiles, n_correlations, n_correlations_tmp !--------------------------------------input ! type(mutable),allocatable :: tables(:) !-----------------------------------internal logical :: newtables integer :: x, n, xbeg, xend, i, m, int1, int2 double precision :: db1 character(LEN=pathname_length+150) :: fullfilename character(LEN=200) :: comment character(LEN=1) :: firstchar character(LEN=100) :: line integer :: col integer :: ios integer, allocatable :: skip(:) integer :: id, posperiod !------------------------------------------- ! open(file_id_common3, file=trim(adjustl(pathname_HS))// & ! & "Expt_tables/analyses.txt", form='formatted') ! read(file_id_common3,'(A)',iostat=ios) datafile(1) ! if(ios.ne.0)then ! close(file_id_common3) !! call system('ls -1 -p '//trim(adjustl(pathname_HS))// & !!& 'Expt_tables/'//trim(adjustl(Exptdir))//' > '//trim(adjustl(pathname_HS))// & !!& 'Expt_tables/analyses.txt') - call system('basename -a `ls -1 -p '//trim(adjustl(pathname_HS))// & -& 'Expt_tables/'//trim(adjustl(Exptdir))//'/*.txt` > HS_analyses.txt') +! call system('basename -a `ls -1 -p '//trim(adjustl(pathname_HS))// & +! & 'Expt_tables/'//trim(adjustl(Exptdir))//'/*.txt` > HS_analyses.txt') + call system('ls -1 -p '//trim(adjustl(pathname_HS))// & +& 'Expt_tables/'//trim(adjustl(Exptdir))//'/*.txt | xargs -L 1 basename > HS_analyses.txt') !! open(file_id_common3, file=trim(adjustl(pathname_HS))//"Expt_tables/analyses.txt",& !!& form='formatted') open(file_id_common3, file="HS_analyses.txt",form='formatted') ! else ! rewind(file_id_common3) ! endif print *, "Reading in measurements from the following datafiles from analysis-set "// & trim(adjustl(Exptdir))//":" n = 0 n_datafiles = 0 do n = n+1 read(file_id_common3,'(A)', iostat=ios) datafile(n) if(ios.ne.0) exit write(*,'(I4,2X,A)') n, datafile(n) enddo n_datafiles = n - 1 close(file_id_common3) allocate(obs(n_datafiles),skip(n_datafiles)) do i=lbound(obs,dim=1), ubound(obs,dim=1) call initializemutable_blank(obs(i)%table) enddo do n=1,n_datafiles skip(n)=11 !! if(withcorrexpsyst) skip(n)=12 open(file_id_common3, file=trim(adjustl(pathname_HS)) //'Expt_tables/'// & & trim(adjustl(Exptdir))//'/' // datafile(n)) do read(file_id_common3,'(A)') comment comment = trim(adjustl(comment)) write(firstchar,'(A1)') comment if(firstchar.ne.'#') then exit else skip(n)=skip(n)+1 endif enddo backspace(file_id_common3) read(file_id_common3,*) obs(n)%id, obs(n)%table%id, obs(n)%obstype read(file_id_common3,'(A)') obs(n)%table%label read(file_id_common3,*) obs(n)%table%collider,obs(n)%table%collaboration, & & obs(n)%table%expt read(file_id_common3,'(A)') obs(n)%table%desc read(file_id_common3,*) obs(n)%table%energy, obs(n)%table%lumi, obs(n)%table%dlumi !--TESTING correlated experimental systematics: !! if(withcorrexpsyst) read(file_id_common3,*) (obs(n)%table%correxpsyst(i),i=1,4) !! write(*,*) "Systematics: ",obs(n)%table%correxpsyst !--END read(file_id_common3,*) obs(n)%table%particle_x, obs(n)%table%mhchisq !--CHECK FOR ASSIGNMENT GROUP AS SECOND COLUMN: read(file_id_common3,'(A)') line call read_in_mass_resolution_and_assignment_group(line, obs(n)%table%deltam,& & obs(n)%table%assignmentgroup) ! write(*,*) "dm, group = ",obs(n)%table%deltam, obs(n)%table%assignmentgroup ! read(file_id_common3,*) obs(n)%table%deltam read(file_id_common3,*) obs(n)%table%xmin, obs(n)%table%xmax, obs(n)%table%sep read(file_id_common3,*) obs(n)%table%Nc, obs(n)%table%eff_ref_mass allocate(obs(n)%table%channel_id_str(obs(n)%table%Nc)) ! allocate(obs(n)%table%channel_id(obs(n)%table%Nc)) allocate(obs(n)%table%channel_p_id(obs(n)%table%Nc)) allocate(obs(n)%table%channel_d_id(obs(n)%table%Nc)) ! read(file_id_common3,*) (obs(n)%table%channel_id(i),i=1,obs(n)%table%Nc) read(file_id_common3,*) (obs(n)%table%channel_id_str(i),i=1,obs(n)%table%Nc) do i=1,obs(n)%table%Nc posperiod = index(obs(n)%table%channel_id_str(i),'.') if(posperiod.eq.0) then if(len(trim(adjustl(obs(n)%table%channel_id_str(i)))).eq.2) then read(obs(n)%table%channel_id_str(i),*) id obs(n)%table%channel_p_id(i) = int((id-modulo(id,10))/dble(10)) obs(n)%table%channel_d_id(i) = modulo(id,10) else write(*,*) " For observable ID = ",obs(n)%id stop " Error: Cannot handle channel IDs!" endif else read(obs(n)%table%channel_id_str(i)(:posperiod-1),*) obs(n)%table%channel_p_id(i) read(obs(n)%table%channel_id_str(i)(posperiod+1:),*) obs(n)%table%channel_d_id(i) endif enddo ! write(*,*) "Production channels = ",obs(n)%table%channel_p_id ! write(*,*) "Decay channels = ",obs(n)%table%channel_d_id ! write(*,*) obs(n)%table%channel_id_str allocate(obs(n)%table%channel_eff(obs(n)%table%Nc)) allocate(obs(n)%table%channel_eff_ratios(obs(n)%table%Nc)) if(obs(n)%table%eff_ref_mass.ge.0D0) then read(file_id_common3,*) (obs(n)%table%channel_eff(i),i=1,obs(n)%table%Nc) else do i=1,obs(n)%table%Nc obs(n)%table%channel_eff(i)=1.0D0 enddo read(file_id_common3,*) endif obs(n)%table%channel_eff_ratios=1.0D0 allocate(obs(n)%table%channel_description(obs(n)%table%Nc,2)) close(file_id_common3) enddo col = 3 do x=1,n_datafiles if(obs(x)%table%deltam.le.obs(x)%table%sep) then write(*,*) "Warning: Mass resolution for " ,trim(obs(x)%table%label), & & "observable number ",x," with ID = ",obs(x)%id, & & " very small - may lead to unreliable results" write(*,*) "Mass resolution = ",obs(x)%table%deltam write(*,*) "Separation = ",obs(x)%table%sep endif obs(x)%table%nx=nint((obs(x)%table%xmax-obs(x)%table%xmin)/obs(x)%table%sep)+1 allocate(obs(x)%table%mu(obs(x)%table%nx,col)) allocate(obs(x)%table%mass(obs(x)%table%nx)) allocate(obs(x)%table%channel_mu(obs(x)%table%Nc,np(obs(x)%table%particle_x))) allocate(obs(x)%table%channel_systSM(obs(x)%table%Nc,np(obs(x)%table%particle_x))) allocate(obs(x)%table%channel_syst(obs(x)%table%Nc,np(obs(x)%table%particle_x))) allocate(obs(x)%table%channel_w(obs(x)%table%Nc,np(obs(x)%table%particle_x))) allocate(obs(x)%table%channel_w_corrected_eff(obs(x)%table%Nc,np(obs(x)%table%particle_x))) enddo do x=1, n_datafiles fullfilename=trim(adjustl(pathname_HS))//'Expt_tables/'//trim(adjustl(Exptdir))//'/'& & //trim(datafile(x)) call read_mutable(obs(x)%table,skip(x),col,fullfilename) enddo deallocate(skip) !NEW: - call system('basename -a `ls -1 -p '//trim(adjustl(pathname_HS))// & -& 'Expt_tables/'//trim(adjustl(Exptdir))//'/*.corr 2>/dev/null` > HS_correlations.txt 2>/dev/null') - call system('rm -rf HS_ncorrelations.txt') + call system('ls -1 -p '//trim(adjustl(pathname_HS))// & +& 'Expt_tables/'//trim(adjustl(Exptdir))//'/*.corr | xargs -L 1 basename > HS_correlations.txt') + +! call system('basename -a `ls -1 -p '//trim(adjustl(pathname_HS))// & +! & 'Expt_tables/'//trim(adjustl(Exptdir))//'/*.corr 2>/dev/null` > HS_correlations.txt 2>/dev/null') +! call system('rm -rf HS_ncorrelations.txt') open(file_id_common3, file="HS_correlations.txt",form='formatted') print *, "Reading in correlations from the following datafiles in analysis-set "// & trim(adjustl(Exptdir))//":" n = 0 n_datafiles = 0 n_correlations = 0 do n = n+1 read(file_id_common3,'(A)', iostat=ios) datafile(n) if(ios.ne.0) exit fullfilename=trim(adjustl(pathname_HS))//'Expt_tables/'//trim(adjustl(Exptdir))//'/'& & //trim(datafile(n)) call system('cat '//trim(adjustl(fullfilename))//' | wc -l > HS_ncorrelations.txt') open(file_id_common2,file="HS_ncorrelations.txt",form='formatted') read(file_id_common2,'(I10)') n_correlations_tmp close(file_id_common2) write(*,'(2I4,2X,A)') n, n_correlations_tmp, datafile(n) n_correlations = n_correlations + n_correlations_tmp enddo n_datafiles = n - 1 close(file_id_common3) allocate(corrlist(n_correlations)) m=0 do n=1,n_datafiles fullfilename=trim(adjustl(pathname_HS))//'Expt_tables/'//trim(adjustl(Exptdir))//'/'& & //trim(datafile(n)) open(file_id_common3,file=fullfilename) do m= m+1 read(file_id_common3,*,iostat=ios) int1, int2, db1 ! write(*,*) m, int1, int2, db1, ios if(ios.ne.0) exit corrlist(m)%obsID1 = int1 corrlist(m)%obsID2 = int2 corrlist(m)%corr = db1 enddo m=m-1 close(file_id_common3) enddo ! allocate(obs(n_datafiles),skip(n_datafiles)) ! do i=lbound(obs,dim=1), ubound(obs,dim=1) ! call initializemutable_blank(obs(i)%table) ! enddo ! ! do n=1,n_datafiles ! skip(n)=11 ! !! if(withcorrexpsyst) skip(n)=12 ! open(file_id_common3, file=trim(adjustl(pathname_HS)) //'Expt_tables/'// & ! & trim(adjustl(Exptdir))//'/' // datafile(n)) ! do ! read(file_id_common3,'(A)') comment ! comment = trim(adjustl(comment)) ! write(firstchar,'(A1)') comment ! if(firstchar.ne.'#') then ! exit ! else ! skip(n)=skip(n)+1 ! endif ! enddo ! backspace(file_id_common3) ! read(file_id_common3,*) obs(n)%id, obs(n)%table%id, obs(n)%obstype ! read(file_id_common3,'(A)') obs(n)%table%label end subroutine initialize_observables !-------------------------------------------------------------------- subroutine read_in_mass_resolution_and_assignment_group(line, dm, group) !-------------------------------------------------------------------- character(LEN=100), intent(in) :: line character(LEN=100), intent(out) :: group double precision, intent(out) :: dm integer :: i, indx, prev, beginning integer :: j, indxstr, prevstr, beginningstr prev = 0 beginning = 1 prevstr = 0 beginningstr = 1 group = "" do i=1,len(line) indx = index('0123456789.', line(i:i)) if (indx.eq.0 .and. prev.gt.0) then read(line(beginning:i-1), *) dm else if (indx.gt.0 .and. prev.eq.0) then beginning = i end if prev = indx indxstr = index('ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz_', line(i:i)) if (indxstr.eq.0 .and. prevstr.gt.0) then read(line(beginningstr:i-1), *) group else if (indxstr.gt.0 .and. prevstr.eq.0) then beginningstr = i end if prevstr = indxstr end do ! group=trim(adjustl(group_tmp)) end subroutine read_in_mass_resolution_and_assignment_group !-------------------------------------------------------------------- subroutine setup_peak_observables !-------------------------------------------------------------------- implicit none integer :: i,j do i=lbound(obs,dim=1),ubound(obs,dim=1) if(obs(i)%obstype.eq.1) then obs(i)%peak%id = obs(i)%id obs(i)%peak%mpeak = obs(i)%table%mass(1) obs(i)%peak%dm = obs(i)%table%deltam obs(i)%peak%mu = obs(i)%table%mu(1,2) obs(i)%peak%mu_original = obs(i)%table%mu(1,2) obs(i)%peak%dmuup = obs(i)%table%mu(1,3) - obs(i)%table%mu(1,2) obs(i)%peak%dmulow = obs(i)%table%mu(1,2) - obs(i)%table%mu(1,1) obs(i)%peak%dlumi = obs(i)%table%dlumi obs(i)%peak%scale_mu = 1.0D0 obs(i)%peak%undo_assignment = 0 obs(i)%peak%assignmentgroup = obs(i)%table%assignmentgroup ! obs(i)%peak%scale_mh = 1.0D0 if(.not.allocated(obs(i)%peak%Higgs_comb)) then allocate(obs(i)%peak%Higgs_comb(np(obs(i)%table%particle_x))) endif if(.not.allocated(obs(i)%peak%Higgses)) then allocate(obs(i)%peak%Higgses(np(obs(i)%table%particle_x))) endif obs(i)%peak%Nc = obs(i)%table%Nc ! have to allocate channel info and set them initially to zero! ! if(.not.allocated(obs(i)%peak%channel_id)) then ! allocate(obs(i)%peak%channel_id(obs(i)%peak%Nc)) ! endif ! obs(i)%peak%channel_id(:) = obs(i)%table%channel_id(:) if(.not.allocated(obs(i)%peak%channel_p_id)) then allocate(obs(i)%peak%channel_p_id(obs(i)%peak%Nc)) endif obs(i)%peak%channel_p_id(:) = obs(i)%table%channel_p_id(:) if(.not.allocated(obs(i)%peak%channel_d_id)) then allocate(obs(i)%peak%channel_d_id(obs(i)%peak%Nc)) endif obs(i)%peak%channel_d_id(:) = obs(i)%table%channel_d_id(:) if(.not.allocated(obs(i)%peak%channel_eff)) then allocate(obs(i)%peak%channel_eff(obs(i)%peak%Nc)) endif obs(i)%peak%channel_eff(:) = obs(i)%table%channel_eff(:) if(.not.allocated(obs(i)%peak%channel_mu)) then allocate(obs(i)%peak%channel_mu(obs(i)%peak%Nc)) endif if(.not.allocated(obs(i)%peak%channel_w_model)) then allocate(obs(i)%peak%channel_w_model(obs(i)%peak%Nc)) endif if(.not.allocated(obs(i)%peak%channel_w)) then allocate(obs(i)%peak%channel_w(obs(i)%peak%Nc)) endif if(.not.allocated(obs(i)%peak%channel_w_corrected_eff)) then allocate(obs(i)%peak%channel_w_corrected_eff(obs(i)%peak%Nc)) endif if(.not.allocated(obs(i)%peak%channel_systSM)) then allocate(obs(i)%peak%channel_systSM(obs(i)%peak%Nc)) endif if(.not.allocated(obs(i)%peak%channel_syst)) then allocate(obs(i)%peak%channel_syst(obs(i)%peak%Nc)) endif do j=1,obs(i)%peak%Nc obs(i)%peak%channel_mu(j)=0.0D0 obs(i)%peak%channel_w(j)=0.0D0 obs(i)%peak%channel_w_corrected_eff(j)=0.0D0 obs(i)%peak%channel_w_model(j)=0.0D0 obs(i)%peak%channel_systSM(j)=0.0D0 obs(i)%peak%channel_syst(j)=0.0D0 enddo endif enddo end subroutine setup_peak_observables !-------------------------------------------------------------------- subroutine read_mutable(t1,skip,col,fullfilename) !-------------------------------------------------------------------- !--------------------------------------input type(mutable) :: t1 integer :: skip,col character(LEN=*) :: fullfilename !-----------------------------------internal integer :: i,n , file_id_1 double precision :: xdummy,xdummy_store !------------------------------------------- file_id_1 = 666 t1%mu=-1.0D0 open(file_id_1, file=(trim(fullfilename))) do i=1,skip read(file_id_1,*) !skip lines enddo xdummy_store = t1%xmin-t1%sep do i=1,t1%nx read(file_id_1,*)xdummy,(t1%mu(i,n),n=1,3) ! checks that x are evenly spaced as expected if((abs(xdummy-xdummy_store-t1%sep).gt.1.0D-7) & & .or.(abs(xdummy-(t1%xmin+dble(i-1)*t1%sep)).gt.1.0D-7))then !! write(*,*)i,t1%id,xdummy,t1%xmin+dble(i-1)*t1%sep write(*,*) "Problem with observable ",t1%id stop 'error in read_mutable (a1)' endif t1%mass(i) = xdummy !! write(*,*) xdummy xdummy_store=xdummy enddo if(abs(xdummy-t1%xmax).gt.1.0D-7)stop 'error in read_mutable (a2)' close(file_id_1) end subroutine read_mutable !-------------------------------------------------------------------- subroutine setup_tablelist !-------------------------------------------------------------------- implicit none integer,allocatable :: tableid(:) integer :: i,j,k,n, ntab, firstpeak_index logical :: tabused, tab_taken_from_mpredobs allocate(tableid(size(obs,dim=1))) do i=lbound(tableid,dim=1),ubound(tableid,dim=1) tableid(i)=-1 enddo ntab=0 do i=lbound(obs,dim=1),ubound(obs,dim=1) tabused=.False. do k=lbound(tableid,dim=1),ubound(tableid,dim=1) if(obs(i)%table%id.eq.tableid(k)) tabused=.True. enddo if(.not.tabused) then ntab=ntab+1 !--Fill first element which is -1 with the table id do k=lbound(tableid,dim=1),ubound(tableid,dim=1) if(tableid(k).eq.-1) then tableid(k)=obs(i)%table%id exit endif enddo endif enddo if(allocated(analyses)) deallocate(analyses) allocate(analyses(ntab)) do i=lbound(analyses,dim=1),ubound(analyses,dim=1) analyses(i)%id = tableid(i) !--Find observables based on this table n=0 do j=lbound(obs,dim=1),ubound(obs,dim=1) if(analyses(i)%id.eq.obs(j)%table%id) then !----Check if it is a peak observable if(obs(j)%obstype.eq.1) then n=n+1 endif endif enddo allocate(analyses(i)%peaks(n)) analyses(i)%Npeaks=n n=0 firstpeak_index=0 tab_taken_from_mpredobs=.False. do j=lbound(obs,dim=1),ubound(obs,dim=1) if(analyses(i)%id.eq.obs(j)%table%id) then !----Check if it is a peak observable if(obs(j)%obstype.eq.1) then n=n+1 analyses(i)%peaks(n) = obs(j)%peak analyses(i)%peaks(n)%Higgses = obs(j)%Higgses allocate(analyses(i)%peaks(n)%channel_w_allH( & & size(obs(j)%table%channel_w,dim=1),size(obs(j)%table%channel_w,dim=2))) allocate(analyses(i)%peaks(n)%channel_w_corrected_eff_allH( & & size(obs(j)%table%channel_w_corrected_eff,dim=1), & & size(obs(j)%table%channel_w_corrected_eff,dim=2))) allocate(analyses(i)%peaks(n)%channel_systSM_allH( & & size(obs(j)%table%channel_systSM,dim=1),size(obs(j)%table%channel_systSM,dim=2))) allocate(analyses(i)%peaks(n)%channel_syst_allH( & & size(obs(j)%table%channel_syst,dim=1),size(obs(j)%table%channel_syst,dim=2))) allocate(analyses(i)%peaks(n)%channel_mu_allH( & & size(obs(j)%table%channel_mu,dim=1),size(obs(j)%table%channel_mu,dim=2))) analyses(i)%peaks(n)%channel_w_allH = obs(j)%table%channel_w analyses(i)%peaks(n)%channel_w_corrected_eff_allH = obs(j)%table%channel_w_corrected_eff analyses(i)%peaks(n)%channel_systSM_allH = obs(j)%table%channel_systSM analyses(i)%peaks(n)%channel_syst_allH = obs(j)%table%channel_syst analyses(i)%peaks(n)%channel_mu_allH = obs(j)%table%channel_mu if(n.eq.1) firstpeak_index=j elseif(obs(j)%obstype.eq.2) then analyses(i)%table = obs(j)%table if(.not.allocated(analyses(i)%Higgses)) then allocate(analyses(i)%Higgses(size(obs(j)%Higgses,dim=1))) endif analyses(i)%Higgses = obs(j)%Higgses do k=lbound(analyses(i)%Higgses,dim=1),ubound(analyses(i)%Higgses,dim=1) analyses(i)%Higgses(k)%mp_test=1 enddo tab_taken_from_mpredobs = .True. endif endif enddo !-If we did not find a mpred observable, take the table from the first peak observable if(.not.tab_taken_from_mpredobs) then if(firstpeak_index.gt.0) then analyses(i)%table=obs(firstpeak_index)%table if(.not.allocated(analyses(i)%Higgses)) then allocate(analyses(i)%Higgses(size(obs(firstpeak_index)%Higgses,dim=1))) endif analyses(i)%Higgses=obs(firstpeak_index)%Higgses endif endif enddo deallocate(tableid) end subroutine setup_tablelist !-------------------------------------------------------------------- subroutine setup_observables !-------------------------------------------------------------------- use usefulbits, only : debug,np,not_a_particle implicit none !-----------------------------------internal integer :: i,c, nobs call initialize_observables call setup_peak_observables !n.b.: setup of mpred observables will be done during the run. do i=lbound(obs,dim=1),ubound(obs,dim=1) if( count(obs%id.eq.obs(i)%id).ne.1)then write(*,*)'the observable id',obs(i)%id,'is repeated' stop 'error in setup_observables (c1)' endif enddo !check to make sure that obs(i)%mutable%particle_x are all particles do i=lbound(obs,dim=1),ubound(obs,dim=1) if(obs(i)%table%particle_x.eq.not_a_particle)then write(*,*)obs(i)%id,'particle_x=not_a_particle.' stop 'error in setup_observables (d1)' endif enddo end subroutine setup_observables !-------------------------------------------------------------------- subroutine setup_LHC_Run1_combination !-------------------------------------------------------------------- use usefulbits, only : file_id_common3 use usefulbits_HS, only : LHCrun1_rates, LHCrun1_correlationmatrix use store_pathname_HS integer :: ios, n, skip, i,j, id double precision :: r, r_low, r_up, dble_tmp character(LEN=100) :: string_tmp skip =3 open(file_id_common3, file=trim(adjustl(pathname_HS))//'Expt_tables/'// & & 'LHC_combination_Run1/LHC_comb_Run1.exp') n=0 do n = n+1 if(n.le.skip) then read(file_id_common3,*, iostat=ios) else read(file_id_common3,*, iostat=ios) id, r_low, r, r_up if(ios.ne.0) exit LHCrun1_rates(n-skip)%channel_id = id LHCrun1_rates(n-skip)%r_low = r_low LHCrun1_rates(n-skip)%r = r LHCrun1_rates(n-skip)%r_up = r_up LHCrun1_rates(n-skip)%dr_low = LHCrun1_rates(n-skip)%r - LHCrun1_rates(n-skip)%r_low LHCrun1_rates(n-skip)%dr_up = LHCrun1_rates(n-skip)%r_up - LHCrun1_rates(n-skip)%r endif enddo close(file_id_common3) skip =0 n = 0 open(file_id_common3, file=trim(adjustl(pathname_HS))//'Expt_tables/'// & & 'LHC_combination_Run1/LHC_comb_Run1.corr') do n = n+1 if(n.le.skip) then read(file_id_common3,*, iostat=ios) else read(file_id_common3,*, iostat=ios) i, j, dble_tmp if(ios.ne.0) exit LHCrun1_correlationmatrix(i,j) = dble_tmp endif enddo do i=1,20 do j=1,20 if(LHCrun1_correlationmatrix(i,j).ne.LHCrun1_correlationmatrix(j,i)) then write(*,*) i, j, LHCrun1_correlationmatrix(i,j),LHCrun1_correlationmatrix(j,i) endif enddo enddo end subroutine setup_LHC_Run1_combination !-------------------------------------------------------------------- subroutine check_available_Higgses(ii) !-------------------------------------------------------------------- implicit none integer, intent(in) :: ii integer :: i,j if(.not.allocated(analyses)) then stop "Error in subroutine check_available_Higgses: analyses not allocated." endif do i=lbound(analyses(ii)%peaks,dim=1),ubound(analyses(ii)%peaks,dim=1) do j=lbound(analyses(ii)%peaks(i)%Higgs_comb,dim=1), & & ubound(analyses(ii)%peaks(i)%Higgs_comb,dim=1) if(analyses(ii)%peaks(i)%Higgs_comb(j).ne.0) then analyses(ii)%Higgses(analyses(ii)%peaks(i)%Higgs_comb(j))%mp_test=0 endif enddo enddo end subroutine check_available_Higgses !-------------------------------------------------------------------- subroutine interpolate_mutable(value, mh_in, mutab, steps) ! Returns interpolated values for the signal strength (lower, central, upper value) ! for the given mass mh_in of mutable mutab. "steps" defines the number of interpolation ! steps. !-------------------------------------------------------------------- implicit none double precision, dimension(3), intent(out) :: value double precision, intent(in) :: mh_in type(mutable), intent(in) :: mutab integer, intent(in) :: steps integer :: i,j, imin, imax double precision :: subsep, mh, closestdistance, mh_closest double precision, dimension(3) :: mu if(mh_in.ge.mutab%xmax) then value(:) = mutab%mu(ubound(mutab%mu,dim=1),:) else if(mh_in.le.mutab%xmin) then value(:) = mutab%mu(lbound(mutab%mu,dim=1),:) else imin = int((mh_in-mutab%xmin)/mutab%sep)+1 subsep = mutab%sep/steps closestdistance = 1000000.0D0 mh = mutab%mass(imin) do j=1, steps mh = mh + subsep mu(:) = mutab%mu(imin,:) + (mutab%mu(imin+1,:)-mutab%mu(imin,:))*j/steps if(abs(mh-mh_in).le.closestdistance) then closestdistance = abs(mh-mh_in) mh_closest = mh value = mu endif enddo endif end subroutine interpolate_mutable !-------------------------------------------------------------------- subroutine deallocate_observables !-------------------------------------------------------------------- implicit none integer :: x if(allocated(obs)) then do x=lbound(obs,dim=1),ubound(obs,dim=1) deallocate(obs(x)%table%mu) deallocate(obs(x)%table%mass) deallocate(obs(x)%table%channel_mu) deallocate(obs(x)%table%channel_id_str) ! deallocate(obs(x)%table%channel_id) deallocate(obs(x)%table%channel_p_id) deallocate(obs(x)%table%channel_d_id) deallocate(obs(x)%table%channel_systSM) deallocate(obs(x)%table%channel_syst) deallocate(obs(x)%table%channel_w) deallocate(obs(x)%table%channel_w_corrected_eff) deallocate(obs(x)%table%channel_eff) deallocate(obs(x)%table%channel_eff_ratios) deallocate(obs(x)%table%channel_description) if(allocated(obs(x)%table%Toys_mhobs)) deallocate(obs(x)%table%Toys_mhobs) if(allocated(obs(x)%table%Toys_muobs)) deallocate(obs(x)%table%Toys_muobs) if(allocated(obs(x)%peak%Higgs_comb)) then deallocate(obs(x)%peak%Higgs_comb) deallocate(obs(x)%peak%Higgses) ! deallocate(obs(x)%peak%channel_id) deallocate(obs(x)%peak%channel_p_id) deallocate(obs(x)%peak%channel_d_id) deallocate(obs(x)%peak%channel_eff) deallocate(obs(x)%peak%channel_mu) deallocate(obs(x)%peak%channel_w_model) deallocate(obs(x)%peak%channel_w) deallocate(obs(x)%peak%channel_w_corrected_eff) deallocate(obs(x)%peak%channel_systSM) deallocate(obs(x)%peak%channel_syst) endif if(allocated(obs(x)%Higgses)) deallocate(obs(x)%Higgses) enddo deallocate(obs) endif end subroutine deallocate_observables !-------------------------------------------------------------------- end module datatables !-------------------------------------------------------------------- \ No newline at end of file