diff --git a/DESCRIPTION b/DESCRIPTION
index ac23a157..2f64adf2 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -55,6 +55,7 @@ Description: Toolkit to perform statistical analyses of correlation
    <https://inspirehep.net/literature/1792113>) are implemented.
    In addition, input/output and plotting routines are available.
 Imports:
+    dplyr,
     abind,
     boot,
     dplyr,
diff --git a/R/cf.R b/R/cf.R
index fd496e4e..8ad139d9 100644
--- a/R/cf.R
+++ b/R/cf.R
@@ -60,7 +60,8 @@ cf_meta <- function (.cf = cf(), nrObs = 1, Time = NA, nrStypes = 1, symmetrised
   return (.cf)
 }
 
-#' Bootstrapped CF mixin constructor
+
+#' Bootstrapped CF mixin constructor 
 #'
 #' @param .cf `cf` object to extend.
 #' @param boot.R Integer, number of bootstrap samples used.
@@ -126,6 +127,7 @@ cf_boot <- function (.cf = cf(), boot.R, boot.l, seed, sim, endcorr, cf.tsboot,
   return (.cf)
 }
 
+
 #' Estimates error from jackknife samples
 #'
 #' Currently this uses the mean over the jackknife samples in order to compute
@@ -497,10 +499,78 @@ gen.block.array <- function(n, R, l, endcorr=TRUE) {
   return(list(starts = st, lengths = lens))
 }
 
+#' Computes the samples for reweighted correlation function
+#'
+#' @param cf `cf` object.
+#' @param rw `rw` object.
+#' @param boot.R Integer
+#' @param boot.l Integer
+#' @param seed Integer
+#' @param sim string
+#' @param endcorr boolean
+#' @export
+bootstrap_rw.cf <- function(cf, rw, boot.R=400, boot.l=2, seed=1234, sim="geom", endcorr=TRUE) {
+  stopifnot(inherits(cf, 'cf_orig'))
+  stopifnot(inherits(rw, 'rw_orig'))
+  stopifnot(inherits(rw, 'rw_meta'))
+  stopifnot(inherits(cf, 'cf_indexed'))
 
 
-#' bootstrap a set of correlation functions
-#' 
+  ##We should also check that the cf object and the rw object contains the same gauge configurations
+
+  stopifnot(rw$conf.index == cf$conf.index)
+ 
+  stopifnot( nrow(cf$cf) == length(rw$conf.index) )
+
+  boot.l <- ceiling(boot.l)
+  boot.R <- floor(boot.R)
+
+  stopifnot(boot.l >= 1)
+  stopifnot(boot.l <= nrow(cf$cf))
+  stopifnot(boot.R >= 1)
+
+  ##Construct correlation function for the reweighting samples
+  rw_cf <- cf
+  rw_cf$cf <- replicate(ncol(cf$cf), rw$rw)
+
+  ## we set the seed for reproducability and correlation
+  old_seed <- swap_seed(seed)
+
+  ## now we bootstrap the correlators*reweighting factor
+  rwcf.tsboot <- boot::tsboot(cf$cf*rw_cf$cf, statistic = function(x){ return(apply(x, MARGIN=2L, FUN=mean))},
+                            R = boot.R, l=boot.l, sim=sim, endcorr=endcorr)
+
+
+  restore_seed(old_seed)
+
+  ## we set the seed for reproducability and correlation
+  old_seed <- swap_seed(seed)
+
+  ## now we bootstrap the reweighting factor
+  rw.tsboot <- boot::tsboot(rw_cf$cf, statistic = function(x){ return(apply(x, MARGIN=2L, FUN=mean))},
+                            R = boot.R, l=boot.l, sim=sim, endcorr=endcorr)
+
+
+  rwcf.tsboot$t0<- rwcf.tsboot$t0/rw.tsboot$t0
+  rwcf.tsboot$t <- rwcf.tsboot$t/rw.tsboot$t
+
+  cf <- cf_boot(cf,
+                boot.R = boot.R,
+                boot.l = boot.l,
+                seed = seed,
+                sim = sim,
+                cf.tsboot = rwcf.tsboot)
+
+  class(cf) <- append(class(cf), 'cfrw_boot')
+
+  class(cf) <- setdiff(class(cf), 'cf_orig')
+
+
+  restore_seed(old_seed)
+
+  return(invisible(cf))
+}
+
 #' bootstrap a set of correlation functions
 #' 
 #' 
@@ -574,9 +644,6 @@ bootstrap.cf <- function(cf, boot.R=400, boot.l=2, seed=1234, sim="geom", endcor
 }
 
 
-
-#' jackknife a set of correlation functions
-#' 
 #' jackknife a set of correlation functions
 #' 
 #' 
@@ -666,8 +733,83 @@ jackknife.cf <- function(cf, boot.l = 1) {
                 resampling_method = 'jackknife')
 
   return (invisible(cf))
+
 }
 
+#' Computes the jackknife samples for reweighted correlation function
+#'
+#' @param cf `cf` object.
+#' @param rw `rw` object.
+#' @param boot.l Integer
+#' @export
+jackknife_rw.cf <- function(cf, rw, boot.l = 1) {
+  stopifnot(inherits(cf, 'cf_orig'))
+  stopifnot(inherits(rw, 'rw_orig'))
+  stopifnot(inherits(rw, 'rw_meta'))
+  stopifnot(inherits(cf, 'cf_indexed'))
+
+  stopifnot(rw1$conf.index == rw2$conf.index)
+
+  ##We should also check that the cf object and the rw object contains the same gauge configurations
+
+  stopifnot( nrow(cf$cf) == length(rw$conf.index) )
+
+
+  stopifnot(boot.l >= 1)
+  boot.l <- ceiling(boot.l)
+
+  ##Construct correlation function for the reweighting samples
+  rw_cf <- cf
+  rw_cf$cf <- replicate(ncol(cf$cf), rw$rw)
+
+
+  ## blocking with fixed block length, but overlapping blocks
+  ## number of observations
+  n <- nrow(cf$cf)
+  ## number of overlapping blocks
+  N <- n-boot.l+1
+
+  
+  numerator <- apply(cf$cf*rw_cf$cf, 2, mean)
+  denominator <- apply(rw_cf$cf, 2, mean)
+  t0 <- numerator/denominator
+
+  t <- array(NA, dim = c(N, ncol(cf$cf)))
+  for (i in 1:N) {
+    ## The measurements that we are going to leave out.
+    ii <- c(i:(i+boot.l-1))
+    ## jackknife replications of the mean
+    t[i, ] <- 
+
+    numerator <- apply(cf$cf[-ii, ]* rw_cf$cf[ ii, ], 2L, mean)
+    denominator <- apply( rw_cf$cf[ ii, ] , 2L, mean )
+    t[i, ] < numerator/denominator 
+  }
+
+
+  cf <- invalidate.samples.cf(cf)
+
+  cf.tsboot <- list(t = t,
+                    t0 = t0,
+                    R = N,
+                    l = boot.l)
+
+
+  cf <- cf_boot(cf,
+                boot.R = cf.tsboot$R,
+                boot.l = cf.tsboot$l,
+                seed = 0,
+                sim = 'geom',
+                cf.tsboot = cf.tsboot,
+                resampling_method = 'jackknife')
+
+  class(cf) <- append(class(cf), 'cfrw_boot')
+
+  class(cf) <- setdiff(class(cf), 'cf_orig')
+
+  return (invisible(cf))
+}
+# Gamma method analysis on all time-slices in a 'cf' object
 #' uwerr.cf
 #' @description
 #' Gamma method analysis on all time-slices in a 'cf' object
@@ -738,26 +880,20 @@ addConfIndex2cf <- function(cf, conf.index) {
   if(is.null(cf$conf.index)) {
     cf$conf.index <- conf.index
   }
+  class(cf) <- append(class(cf), 'cf_indexed')
   return(cf)
 }
 
-
-
-#' Combine statistics of two cf objects
-#' 
-#' \code{addStat.cf} takes the raw data of two \code{cf} objects and combines
-#' them into one
+#' Combine correlation function from different replicas
 #' 
-#' Note that the two \code{cf} objects to be combined need to be compatible.
-#' Otherwise, \code{addStat.cf} will abort with an error.
-#' 
-#' @param cf1 the first of the two \code{cf} objects to be combined
-#' @param cf2 the second of the two \code{cf} objects to be combined
-#' @return an object of class \code{cf} with the statistics of the two input
-#' \code{cf} objects combined
-#' @author Carsten Urbach, \email{curbach@@gmx.de}
-#' @seealso \code{\link{cf}}
-#' @keywords correlation function
+#' @param cf1 `cf` object: correlation function for replicum A
+#' @param cf2 `cf` object: correlation function for replicum B
+#' @param reverse1 `boolean` After the bifurcation point one of
+#'                           the replicas (chain of correlation 
+#'                           functions in simulation time) has  
+#'                           to be reversed.
+#' @param reverse2 `boolean`
+#'
 #' @examples
 #'
 #' data(samplecf)
@@ -766,7 +902,7 @@ addConfIndex2cf <- function(cf, conf.index) {
 #' cfnew <- addStat.cf(cf1=samplecf, cf2=samplecf)
 #' 
 #' @export
-addStat.cf <- function(cf1, cf2) {
+addStat.cf <- function(cf1, cf2, reverse1=TRUE,reverse2=FALSE) {
   stopifnot(inherits(cf1, 'cf'))
   stopifnot(inherits(cf2, 'cf'))
 
@@ -780,6 +916,9 @@ addStat.cf <- function(cf1, cf2) {
   stopifnot(inherits(cf1, 'cf_meta'))
   stopifnot(inherits(cf2, 'cf_meta'))
 
+  ##Either both should have an index or none of them
+  stopifnot(inherits(cf1, 'cf_indexed') == inherits(cf1, 'cf_indexed') )
+
   stopifnot(cf1$Time == cf2$Time)
   stopifnot(dim(cf1$cf)[2] == dim(cf2$cf)[2])
   stopifnot(cf1$nrObs == cf2$nrObs )
@@ -787,8 +926,37 @@ addStat.cf <- function(cf1, cf2) {
 
   cf <- cf1
 
-  cf$cf <- rbind(cf1$cf, cf2$cf)
-  cf$icf <- rbind(cf1$icf, cf2$icf)
+  cf1_temp<- cf1$cf
+  icf1_temp <- cf1$icf
+  if (reverse1 == TRUE){
+    apply(cf1_temp,2,rev)
+    if ( has_icf(cf1)){
+      apply(icf1_temp,2,rev)
+    }
+  }
+  cf2_temp <- cf2$cf
+  icf2_temp <- cf2$icf
+  if (reverse2 == TRUE){
+    apply(cf2_temp,2,rev)
+    if ( has_icf(cf2)){
+      apply(icf2_temp,2,rev)
+    }
+  }
+  if (inherits(cf1, 'cf_indexed')){
+    conflist_temp1 <- cf1$conf.index 
+    if (reverse1 == TRUE){
+      conflist_temp1 <- rev(conflist_temp1) 
+    }
+    conflist_temp2 <- cf2$conf.index
+    if (reverse2 == TRUE){
+      conflist_temp2 <- rev(conflist_temp2)
+    }
+    cf$conf.index <- c(conflist_temp1,conflist_temp2)
+  }
+     
+
+  cf$cf <- rbind(cf1_temp, cf2_temp)
+  cf$icf <- rbind(icf1_temp, icf2_temp)
 
   cf <- invalidate.samples.cf(cf)
 
@@ -958,7 +1126,6 @@ apply_elementwise.cf <- function(cf1, cf2, `%op%` = `/`) {
     stopifnot(cf1$Time == cf2$Time)
 
     cf$cf <- cf1$cf %op% cf2$cf
-
     if( has_icf(cf1) | has_icf(cf2) ){
       stopifnot(has_icf(cf1) & has_icf(cf2))
       cf$icf <- cf1$icf %op% cf2$icf
@@ -1317,6 +1484,40 @@ shift.cf <- function(cf, places) {
   return(invisible(cf))
 }
 
+
+avg.sparsify.cf <- function(cf, avg, sparsity){
+  if(!any(class(cf) == "cf")) {
+    stop("avg.sparsify.cf: Input must be of class 'cf'\n")
+  }
+  if( sparsity > 1 ){
+    stop("avg.sparsify.cf: Sparsification has not been implemented yet, only sparsity = 1 supported for now.\n")
+  }
+
+  Lt <- cf$Time
+  if(cf$symmetrised){
+    Lt <- cf$Time/2 + 1
+  }
+  nmeas <- length(cf$cf) / Lt
+  targetstats <- nmeas / avg
+
+  cf2 <- invalidate.samples.cf(cf)
+  cf2$cf <- cf$cf[(1:targetstats),]
+  cf2$icf <- cf$icf[(1:targetstats),]
+  
+  for( m_idx in 1:targetstats ){
+    from <- (m_idx-1)*avg + 1
+    to <- m_idx*avg
+    avg_idx <- c(from:to)
+    cf2$cf[m_idx,] <- apply(X = cf$cf[avg_idx,],
+                            MARGIN = 2,
+                            FUN = sum)/avg
+    #cf2$icf[m_idx,] <- apply(X = cf$icf[avg_idx,],
+    #                         MARGIN = 2,
+    #                         FUN = sum)/avg
+  }
+  return(cf2)
+}
+
 #' Invalidate samples
 #'
 #' When a correlation function is modified, any resampling should be
diff --git a/R/cyprus_readutils.R b/R/cyprus_readutils.R
index d2c76cb5..8f40e60e 100644
--- a/R/cyprus_readutils.R
+++ b/R/cyprus_readutils.R
@@ -28,6 +28,7 @@ cyprus_make_key_scalar <- function(istoch, loop_type, cid = 4, accumulated = FAL
         )
 }
 
+
 #' @title HDF5 key for Cyprus CalcLoops derivative-type loops
 #' @description Generates an HDF5 key (full path) for the derivative
 #'              type loops from the Cyprus CalcLoops application.
@@ -65,6 +66,91 @@ cyprus_make_key_vector <- function(istoch, loop_type, dir, cid = 4, accumulated
 }
 
 
+
+#' @title HDF5 key for Cyprus Calc2pt-and Calc2pt-3pt baryon contractions
+#' @description Generates an HDF5 key (full path) for the baryon 2pt -3pt 
+#'              correlation function produced by plegma postprocessed by 
+#'              extract2.py python filtering routine
+#'              (one usually filters out some momentum combinations)
+#'              and merge.py python merge routine (one merges all the 
+#'              source position and gauge configuration). 
+#'           
+#' @param type_of_correlation_function twopt or threept
+#'              generated by PLEGMA
+#' @param confnumber String, index of the gauge field configuration
+#' 
+#' @param source_location String, sx??sy??sz??st??? format string that 
+#'              contains the position o
+#'      
+#' @param baryon_type String that charachterizes the type of the baryon
+#'              DeltaMn, Delta0, DeltaPl, DeltaPlPl 
+#'
+#' @param interpolator_type: String that determines the type of the 
+#'               interpolator: Comes from PLEGMA, possible values
+#                Pp_Cgi_Cgi
+#cyprus_make_key_baryon <- function( type_of_correlation_function, confnumber, source_location, baryon_type, interpolator_type ){
+
+#  return (sprintf("/%s/%s/%s/%s/%s", type_of_correlation_function,
+#                                                                confnumber,
+#                                                                source_location,
+#                                                                baryon_type,
+#                                                                interpolator_type))
+#}
+#/0280_SS_gN50a4p_aN50a0p5/sx09sy15sz18st37/baryons_u[+2.5e-03]d[-2.5e-03]s[+1.5e-02]/OmegaMn/Pp_Cgi_Cgi Dataset {64}
+
+cyprus_make_key_baryon <- function( confnumber, type_of_correlation_function, fermion_smearing_string,gauge_smearing_string, source_location, baryon_udsc_string, baryon_type, interpolator_type ){
+
+  return (sprintf("%s_%s_%s_%s/%s/%s/%s/%s",confnumber,
+                                     type_of_correlation_function,
+                                     fermion_smearing_string,
+                                     gauge_smearing_string,
+                                     source_location,
+                                     baryon_udsc_string,
+                                     baryon_type,
+                                     interpolator_type))
+}
+
+#' @title HDF5 key for Cyprus piNdiagramms- 2pt baryon contractions
+#' @description Generates an HDF5 key (full path) for the baryon 2pt
+#'              Generated by the piNdiagramms executable
+#'              The key encodes the type of the baryon and the source 
+#'              position, note that no more post-processing is needed,
+#'              the correlation function is already shifted relativ to the
+#'              source. 
+#'           
+#' @param type_of_baryon twopt it can be D or N
+#' 
+#' @param source_location String, sx??sy??sz??st??? format string that 
+#'              contains the position o
+cyprus_make_key_scattering2pt <- function( source_location, type_of_baryon ){
+
+  return (sprintf("%s/%s",source_location,
+                           type_of_baryon))
+}
+
+
+#' @title HDF5 key for Cyprus piNdiagramms- 4pt baryon contractions
+#' @description Generates an HDF5 key (full path) for the baryon 2pt
+#'              Generated by the piNdiagramms executable
+#'              The key encodes the type of the contraction (B,W,Z), the
+#'              source position for the point source propagator and the 
+#'              momentum in the sequential propagators
+#'           
+#' @param type_of_correlation can be B1,B2,W1,W2,W3,W4,Z1,Z2,Z3,Z4
+#' 
+#' @param source_location String, sx??sy??sz??st??? format string 
+#'
+#' @param momentum string specifying the momentum in the sequential 
+#'        propagator for example: pi2=1_-1_1
+cyprus_make_key_scattering4pt <- function( source_location, momentum_i2, type_of_correlation ){
+
+  return (sprintf("%s/pi2=%d_%d_%d/%s",source_location,
+                             momentum_i2[1],momentum_i2[2],momentum_i2[3],
+                             type_of_correlation))
+}
+
+
+
 #' @title read HDF5 loop files in the Cyprus CalcLoops format
 #' @description The CalcLoops code produces HDF5 files which contain
 #'              a matrix of momenta and the data for the loops (without
@@ -202,7 +288,6 @@ cyprus_read_loops <- function(selections, files, Time, nstoch,
     stopifnot( !is.null(project_gamma) )
     stopifnot( all( names(selections) == names(project_gamma) ) )
     selected_vec_loop_types <- which( names(selections) %in% vector_loop_types )
-
     incorrectly_sized_projectors <- NULL
     for( loop_type in selected_vec_loop_types ){
       if( length(project_gamma[[loop_type]] != 4) ){
@@ -554,4 +639,668 @@ cyprus_read_loops <- function(selections, files, Time, nstoch,
   } # for( loop_type )
   return(rval)
 }
+#' @export
+#' @title read HDF5 contraction files in the Cyprus Calc2pt-3pt format
+#' @description The Calc2pt function produces all types of baryon 2 and 3 pt
+#'              functions for different momentum interpolating operators,  
+#'              and source positions.
+#'
+#'              In this routine we read the types of baryon contractions 
+#'              specified by selections.
+#'              For the moment here you can only define the type
+#'              of interpolators and the type of correlation function
+#'              currently the routine also assumes that all the datasets
+#'              for different configurations, different source positions, etc.
+#'              are merged together in one hdf5 file.
+#'              The list of gauge congfigurations and source positions does
+#'              not have to be specified.
+#'              i.e we read all available data from the hdf5 file.
+#'              Currently the routine only reads one momentum
+#'
+#' @param selections Named list with names from the list 'DeltaMn','Delta0',etc.
+#'                   The elements of this list are in turn expected
+#'                   be data frames of the form
+#'                     \tabular{rr}{
+#'                       \strong{interp} \tab \strong{type} \cr
+#'                       "Cgi_Cgi"    \tab "twop"        \cr
+#'                       "Cgi_Cgig0"  \tab "twop"        \cr
+#'                       ...             \tab ...         
+#'                     }
+#'                   specifying the interpolating operators and contractions types
+#'                   for each baryon type.
+#' @param symmetrize bool space time reflection symmetries of the action and the 
+#'                        anti periodic boundary conditions in the temporal direction 
+#'                        for the quark field at zero momentum imply, C_x^+(t)=-C_x^{-}(T-t)
+#'                        in order to decrease errors we average correlators in the forward 
+#'                        backward direction and define
+#'                        C_x(t)=C_x^+(t)-C_x^-(T-t)
+#' @param files Vector of strings, list of HDF5 files to be processed.
+#'                    Currently only one element supported
+#' @param Time Integer, time extent of the lattice.
+#' @param verbose Boolean, output I/O time per file. Requires 'tictoc' package.
+#' @return 
+#'         raw_cf object containing all the correlation functions.
+#'         Find it useful temporary for using the gevp
+
+#' @export
+cyprus_read_baryon_correlation <- function(selections, files, symmetrize, Time, verbose = FALSE){
+  rhdf5_avail <- requireNamespace("rhdf5")
+  dplyr_avail <- requireNamespace("dplyr")
+  stringr_avail <- requireNamespace("stringr")
+  if( !rhdf5_avail | !dplyr_avail | !stringr_avail){
+    stop("The 'dplyr', 'stringr'  and 'rhdf5' packages are required to use this function!\n")
+  }
+  if( verbose ){
+    tictoc_avail <- requireNamespace("tictoc")
+    if( !tictoc_avail ){
+      stop("Time reporting requires the 'tictoc' package!")
+    }
+  }
+ 
+  selected_baryon_types <- names(selections)
+
+  for( ifile in 1:length(files) ){
+    f <- files[ifile]
+    if(verbose) tictoc::tic(sprintf("Reading %s",f))    
+    h5f <- rhdf5::H5Fopen(f, flags = "H5F_ACC_RDONLY")
+    #First we extract the configuration id-s
+    temporary1 <- rhdf5::h5ls(h5f)
+
+    # Extracting the type of the correlation function we want possibility twop , threep, we extract all other information from the hd5f file
+    # for now we have only real support for the twop case
+
+
+    #checking whether the wanted baryon type is available in the file
+    avail_baryon_types <- unlist( lapply( selected_baryon_types, function(x){ x %in% temporary1$name } ) )
+    if( any( !avail_baryon_types ) ){
+        msg <- sprintf("Some selected baryon types could not be found in %s:\n %s",
+                       f,
+                       do.call(paste, as.list( selected_baryon_types[!avail_baryon_types] ) )
+                     )
+        stop(msg)
+    }
+
+    type_of_correlation_function <- as.vector(selections[[avail_baryon_types[1]]]$type)[1]
+
+    filtering_pattern <- "/...."
+    temporary3 <- unique(temporary1$group %>% str_extract(pattern=filtering_pattern))
+    gauge_conf_list <- temporary3[!is.na(temporary3)] 
+
+    interpolator_list_length <- length(selections[[avail_baryon_types[1]]]$interp)
+
+    # In case we symmetrize the correlator we do not store separately
+    # the plus and minus projection, and we store the timeslices only
+    # from 1 to Time/2+1
+    if (symmetrize == TRUE){
+
+      final_val <- array(as.complex(NA), dim=c(length(gauge_conf_list), length(selected_baryon_types)*interpolator_list_length*(Time/2+1)))
+
+    }
+    # On the other hand if we do not symmetrize we store the plus
+    # and minus correlator separately for the full time extent
+    else{
+      final_val <- array(as.complex(NA), dim=c(length(gauge_conf_list), length(selected_baryon_types)*2*interpolator_list_length*Time))
+    }
+    for (gauges in gauge_conf_list){
+
+      #we wanted to extract all source position for a given config
+
+      temporary1 <- rhdf5::h5ls(h5f)
+      temporary1 <- temporary1$group[grepl(gauges,temporary1$group)]
+      temporary2 <- unique(temporary1 %>% str_extract(pattern="sx..sy..sz..st..."))
+      source_position_list <- temporary2[!is.na(temporary2)]
+      
+      for( baryon_type in selected_baryon_types ){
+
+         
+        for(interp in selections[[ baryon_type ]]$interp){ 
+
+          if (symmetrize == TRUE){
+
+            rval <- array(as.complex(NA), dim=c(length(source_position_list),
+                                                                        Time/2+1                                                 ))
+          }
+          else {
+            rval1 <- array(as.complex(NA), dim=c(length(source_position_list),
+                                                                        Time                                                 ))
+            rval2 <- array(as.complex(NA), dim=c(length(source_position_list),
+                                                                        Time                                                 ))
+
+          }
+          sp_num <- 0
+          for (sp in source_position_list){
+            sp_string <- strsplit(sp, "/")[[1]][[1]]
+
+            #Reading the positive
+            interp_pp <- sprintf("Pp_%s", interp)
+
+
+            key <- cyprus_make_key_baryon( gauges,
+                                           as.vector(selections[[baryon_type]]$type)[1], 
+                                           as.vector(selections[[baryon_type]]$fermion_smearing_string)[1], 
+                                           as.vector(selections[[baryon_type]]$gauge_smearing_string)[1],
+                                           sp_string, 
+                                           as.vector(selections[[baryon_type]]$baryon_udsc_string)[1], 
+                                           baryon_type,
+                                           interp_pp )
+            print(key)
+
+            # read the data, which comes in the ordering
+            # time, complex
+            data_pp <- h5_get_dataset(h5f, key)
+            if (dim(data_pp) == 0) 
+              next;
+            sp_num <- sp_num+1
+            #temp<-str_split(sp,"st")[[1]][[2]]
+            #ts <- as.double(str_split(temp,"/")[[1]][[1]])
+            #if (length(dim(data_pp))==1){
+            #             data_pp[(Time-ts):Time] <- -1*data_pp[(Time-ts):Time]
+            #    } else if (length(dim(data_pp))==3){
+            #             data_pp[,1,(Time-ts):Time] <- -1*data_pp[,1,(Time-ts):Time]
+            #    } else {
+            #             data_pp[,(Time-ts):Time] <- -1*data_pp[,(Time-ts):Time]
+            #    }
+
+
+            #Reading the negative
+            interp_pm <- sprintf("Pm_%s", interp)
+            key <- cyprus_make_key_baryon( gauges,
+                                           as.vector(selections[[baryon_type]]$type)[1],
+                                           as.vector(selections[[baryon_type]]$fermion_smearing_string)[1],
+                                           as.vector(selections[[baryon_type]]$gauge_smearing_string)[1],
+                                           sp_string,
+                                           as.vector(selections[[baryon_type]]$baryon_udsc_string)[1],
+                                           baryon_type,
+                                           interp_pm )
+
+
+            # read the data, which comes in the ordering
+            # time, complex
+            data_pm <- h5_get_dataset(h5f, key)
+            #if (length(dim(data_pm))==1){
+            #             data_pm[(Time-ts):Time] <- -1*data_pm[(Time-ts):Time]
+            #} else if (length(dim(data_pm))==3){
+            #             data_pm[,1,(Time-ts):Time] <- -1*data_pm[,1,(Time-ts):Time]
+            #} else {
+            #             data_pm[,(Time-ts):Time] <- -1*data_pm[,(Time-ts):Time]
+            #}
+
+             
+            if (symmetrize == TRUE){
+              reversed <- data_pm[,(dim(data_pm)[length(dim(data_pm))]):1] 
+              
+              for (line in 2:(Time/2+1)){
+                data_pp[,line] <- data_pp[,line]-reversed[,line-1] 
+              }
+              rval[which(sp==source_position_list),1:(Time/2+1)] <- complex(real=data_pp[1,1:(Time/2+1)], imaginary=data_pp[2,1:(Time/2+1)])
+            } else{
+              rval1[which(sp==source_position_list),1:Time] <- complex(real=data_pp[1,1:Time], imaginary=data_pp[2,1:Time])
+              rval2[which(sp==source_position_list),1:Time] <- complex(real=data_pm[1,1:Time], imaginary=data_pm[2,1:Time])    
+            }
+          } # source_positions 
+          if (symmetrize == TRUE){
+
+            for (nt in 1:Time/2+1){
+              tmpr <- 0
+              tmpi <- 0
+              for (line in 1:sp_num){
+
+                tmpr <- tmpr + Re(rval[line,nt])
+                tmpi <- tmpi + Im(rval[line,nt])
+
+              }
+              tmpr <- tmpr/ length(source_position_list)
+              tmpi <- tmpi/ length(source_position_list)
+
+              final_val[which(gauges==gauge_conf_list),(Time/2+1)*interpolator_list_length*(which(baryon_type==selected_baryon_types)-1)+(Time/2+1)*(which(interp==selections[[ baryon_type ]]$interp)-1)+nt] <- complex(real=tmpr, imaginary=tmpi)
+
+            }
+
+          }
+          else{
+            for (nt in 1:Time){
+              tmpr <- 0 
+              tmpi <- 0
+              for (line in 1:sp_num){
+
+                tmpr <- tmpr + Re(rval1[line,nt])
+                tmpi <- tmpi + Im(rval1[line,nt])
+
+              }
+              tmpr <- tmpr/ sp_num
+              tmpi <- tmpi/ sp_num
+
+              final_val[which(gauges==gauge_conf_list),2*Time*interpolator_list_length*(which(baryon_type==selected_baryon_types)-1)+2*Time*(which(interp==selections[[ baryon_type ]]$interp)-1)+nt] <- complex(real=tmpr, imaginary=tmpi)
+
+              tmpr <- 0
+              tmpi <- 0
+              for (line in 1:sp_num){
+
+                tmpr <- tmpr + Re(rval2[line,nt])
+                tmpi <- tmpi + Im(rval2[line,nt])
+
+              }
+              tmpr <- tmpr/ sp_num
+              tmpi <- tmpi/ sp_num
+
+              final_val[which(gauges==gauge_conf_list),Time*2*interpolator_list_length*(which(baryon_type==selected_baryon_types)-1)+2*Time*(which(interp==selections[[ baryon_type ]]$interp)-1)+Time+nt] <- complex(real=tmpr, imaginary=tmpi)
+
+            }# time
+          } #symmetrize
+        } # interpolating operators
+      } # baryon type
+    } # gauge conf
+    H5Fclose(h5f)
+    if(verbose) tictoc::toc()
+  } # ifile
+
+  return(final_val)
+
+}
+#' @export
+#' @title read HDF5 unprojected contraction files in the Cyprus piNdiagramms format
+#' @description The piNdiagramms function produces all types of baryon
+#'              correlation function for isospin 3/2  for different momentum 
+#               combinations, interpolating operators and source positions.
+#'
+#'              In this routine we read the types of baryon contractions 
+#'              specified by selections.
+#'              For the moment here you can only define the 
+#'                     a;type of baryon (N,D)
+#'                     b;gamma structures (Cg1,Cg2,Cg2)
+#'                     c;momentum (0,0,0)
+#' @param selections Named list with names from the list 'D','N'
+#'                   The elements of this list are in turn expected
+#'                   be data frames of the form
+#'                     \tabular{rrrr}{
+#'                       \strong{interp} \tab \strong{px} \tab \strong{py} \tab \stron{pz} \cr
+#'                       "Cg1,Cg2,Cg3"  \tab "0" \tab "1" \tab "1" \cr
+#'                       "Cg1,Cg2,Cg3"  \tab "0" \tab "1" \tab "0" \cr
+#'                       ...        \tab ... \tab ... \tab ...         
+#'                     }
+#'                   specifying the interpolating operators and momenta 
+#'                   for each baryon type.
+#' @param files Vector of strings, list of HDF5 files to be processed.
+#'                    Currently only one element supported
+#' @param Time Integer, time extent of the lattice.
+#' @param verbose Boolean, output I/O time per file. Requires 'tictoc' package.
+#' @return 
+#'         raw_cf object containing all the correlation functions.
+#'         Find it useful temporary for using the gevp
+
+#' @export
+cyprus_read_scattering_2pt_correlation_all <- function(selections, files, Time, verbose = FALSE){
+  rhdf5_avail <- requireNamespace("rhdf5")
+  dplyr_avail <- requireNamespace("dplyr")
+  if( verbose ){
+    tictoc_avail <- requireNamespace("tictoc")
+    if( !tictoc_avail ){
+      stop("Time reporting requires the 'tictoc' package!")
+    }
+  }
+
+  selected_baryon_types <- names(selections)
+  #print(selected_baryon_types)
+
+  #we store the real parts and imaginary parts for
+  #the different source positions separately
+  data_final_real <- NULL
+  data_final_imag <- NULL
+
+  #loop on the list of files provided in the input
+  for( ifile in 1:length(files) ){
+    f <- files[ifile]
+    if(verbose) tictoc::tic(sprintf("Reading %s",f))
+
+    #open the file with read only argument
+    h5f <- rhdf5::H5Fopen(f, flags = "H5F_ACC_RDONLY")
+
+    #checking whether the wanted baryon type is available in the file
+    temporary1 <- rhdf5::h5ls(h5f)
+    avail_baryon_types <- unlist( lapply( selected_baryon_types, function(x){ x %in% temporary1$name } ) )
+    if( any( !avail_baryon_types ) ){
+        msg <- sprintf("Some selected baryon types could not be found in %s:\n %s",
+                       f,
+                       do.call(paste, as.list( selected_baryon_types[!avail_baryon_types] ) )
+                     )
+        stop(msg)
+    }
+
+    #First we extract the source positions contained in the file
+    h5filedatasetcontent <- rhdf5::h5ls(h5f)
+    #print(h5filedatasetcontent)
+
+    filtering_pattern <- "/...."
+    temporary <- unique(h5filedatasetcontent$group %>% str_extract(pattern=filtering_pattern))
+    gauge_conf_list <- temporary[!is.na(temporary)]
+    #print(gauge_conf_list)
+
+    for (gauge  in  gauge_conf_list){
+
+      #loop over all the source positions in the file, 
+      #for each sourceposition we write out seperately "mvec" dataset
+      #that is the reason for the division by 2
+      #temporary <- h5filedatasetcontent$group[grepl(gauge,h5filedatasetcontent$group)]
+      #temporary2 <- temporary[grepl("sx",temporary)]
+      #sourceposition_list <- temporary2[!grepl("pi2",temporary2)]
+      #for (sourceposition in 1:length(sourceposition_list)){
+
+      #key <- cyprus_make_key_scattering2pt( sourceposition_list[sourceposition], "mvec" )
+
+      # read the available momenta
+      # array of 3 integer number
+      #data_mom <- h5_get_dataset(h5f, key)
+      keytag <- sprintf("%s/%s",gauge,
+                                selected_baryon_types)
+      #key <- cyprus_make_key_scattering2pt( sourceposition_list[sourceposition], keytag )
+
+        # read the data
+
+      data_corr <- h5_get_dataset(h5f, keytag) 
+      #str(data_corr)
+        
+      #if by mistake you have not applied the boundary condition
+      #ts<-as.double(str_split(sourceposition_list[[1]],"st")[[1]][[2]])
+      #data_corr[,,,(Time-ts):Time] <- -1*data_corr[,,,(Time-ts):Time]
+      #print(dim(data_corr))
+        
+      # convert it to a two-dimensional format, in such a way
+      # that the first index corresponds to the gauge configuration (and source position)
+      # and the second index corresponds to everything else
+      #res <- matrix(data_corr, prod(dim(data_corr)[1:3]),dim(data_corr)[5])
+      dim(data_corr)<- c(2,16,Time)
+ 
+      #Doing possible spin projection
+      data_projection <- data_corr
+      #Doing the symmetrization
+      #C+-B=<0.25*Tr(1+-g4)*C>
+      #Csymmetried=C+B(t)-C-B(T-t)
+
+
+      partial_array<- array(0,dim=c(2,Time/2+1))
+
+      #Obtaining C(T-t)
+
+      reversed <- data_projection[,,Time:1]
+
+      partial_array[,1] <- 0.25*(data_projection[,1,1]+data_projection[,6,1])
+      for (line in 2:(Time/2+1)){
+        partial_array[,line] <- 0.25*(data_projection[,1,line]+data_projection[,6,line])
+        partial_array[,line] <- partial_array[,line]-0.25*(reversed[,11,line-1]+reversed[,16,line-1])
+      }
+
+      #convert the final array into two dimensional
+      #merge the dimension time, gamma and momentum
+      #leave the complex (real and imag) intact
+      #real part will be equal to cf
+      #complex part will be equal to icf
+      tmpprod <- prod(dim(partial_array))
+      dim(partial_array) <- c(2, tmpprod/2)
+      data_final_real <- rbind(data_final_real, partial_array[1,])
+      data_final_imag <- rbind(data_final_imag, partial_array[2,])
+   }
+
+   rhdf5::H5Fclose(h5f)
+  }
+  attr(data_final_real,"dimnames") <- NULL
+  attr(data_final_imag,"dimnames") <- NULL
+
+  cf <- cf_meta(nrObs=1,Time=Time,nrStypes=1,symmetrised=TRUE)
+  cf$symmetrised=TRUE
+  cf <- cf_orig(cf, cf=data_final_real,icf=data_final_imag)
+  return(invisible(cf))
+
+}
+
+
+
+
+      #looking for the momentum combinations
+      #filtering out the momentum
+      #final_indices <- NULL
+      #sink_source_t_gamma5 <- NULL
+      #works now for only one combination
+      #for (differentcombinations in 1: length(selections[[avail_baryon_types]]$interp)){
+
+          #determining the indices of momentum
+          #for(mom in 1:ncol(data_mom)){
+          #  if (all.equal(data_mom[,mom],
+          #                    c(selections[[selected_baryon_types]]$px[differentcombinations],
+          #                      selections[[selected_baryon_types]]$py[differentcombinations],
+          #                      selections[[selected_baryon_types]]$pz[differentcombinations]))==TRUE){
+          #  momindex <- mom
+          #  break;
+          #  }
+          #}
+          #getting the available gamma structures
+          #descrip <- rhdf5::h5readAttributes(h5f,key,"description")
+          #tmp <- strsplit(descrip$description,"/")
+          #tmp1<- strsplit(tmp[[1]][[4]],"}")
+          #tmp2 <- gsub(",\\{","",tmp1[[1]][[3]])
+          #gamma_source_list <- strsplit(tmp2,",")[[1]]
+          #gammalist <- strsplit(as.character(selections[[selected_baryon_types]]$interp[differentcombinations]),",")[[1]]
+
+          #Determining the indices of the gamma structures
+          #indices <- NULL
+          #error_gamma <- TRUE
+          #list of gamma structures at the source
+          #for (line in 1:length(gammalist)){
+          #  for (line2 in 1:length(gamma_source_list)){
+          #    if (all.equal(gammalist[line], gamma_source_list[line2])==TRUE){
+          #      indices <- c(indices,line2);
+          #      error_gamma <- FALSE
+          #      break;
+          #    }
+          #  }
+          #}
+          #if we have a gamma structure in the input that cannot be found in the
+          # hdf5 file we return with error message
+          #if( error_gamma == TRUE ){
+          #  msg <- sprintf("Some selected gamma structures could not be found in %s:\n %s",
+          #             f,
+          #             do.call(paste, as.list( gamma_source_list ) )
+          #           )
+          #  stop(msg)
+          #}
+
+          #Determining the indices of all the correlation functions
+          #for example if the input is cg1,cg2,cg3 then we need the following
+          #matrix
+          #c( (cg1,cg1)  (cg1,cg2) (cg1,cg3) )
+          #c( (cg2,cg1)  (cg2,cg2) (cg2,cg3) )
+          #c( (cg3,cg1)  (cg3,cg2) (cg3,cg3) )
+          #Note that here we do not project out, we return all the spin components
+          #with real and imaginary part
+
+        #  indlen <- length(indices)
+        #  for (line2 in 1:indlen){
+        #    if (gamma_source_list[indices[line2]]=="C" || gamma_source_list[indices[line2]]=="Cg4"){
+        #      sink_source_t_gamma5 <- c( sink_source_t_gamma5, TRUE)
+        #    }
+        #    else {
+        #      sink_source_t_gamma5 <- c( sink_source_t_gamma5, FALSE)
+        #    }
+        #  }
+        #
+        #  for (line1 in 1:indlen){
+        #    for (line2 in 1:indlen){
+        #      final_indices <- c(final_indices, (indices[line1]-1)*length(gamma_source_list)+(indices[line2]))
+        #    }
+        #  }
+        #}
+        #data_without_external_gamma<- data_corr[,,final_indices,]
+        #str(data_corr[,,1,])
+        #str(data_without_external_gamma)
+        #stop("testing")
+        
+        #data_tmp <- data_corr[,,final_indices,]
+        #indices <- final_indices
+        #Multiplication by external gammas
+        #indlen<- sqrt(length(indices))
+        #for (line1 in 1:(indlen)){
+        #  for (line2 in 1:(indlen)){
+        #    if (length(dim(data_tmp))==4){
+        #      data_tmp[,,(line1-1)*indlen+line2,] <- data_without_external_gamma[,,(line1-1)*indlen+line2,]
+        #    } else if (length(dim(data_tmp))==3){
+        #      data_tmp[,,] <- data_without_external_gamma[,,]
+        #    }
+        #    ##if the gamma structure at the source is either C or Cg4 we have to multiply the correlator in space space
+           # if (line2==1 & line1==1 & length(dim(data_tmp))==4){
+            #  print("Testing")
+            #  print(data_tmp[,,1,])
+            #}
+            #if (line2==1 & line1==1 & length(dim(data_tmp))==3){
+            #  print("Testing")
+            #  print(data_tmp[,,])
+            #}
+
+            ##by gamma5 from the left 
+         #   if (sink_source_t_gamma5[line2]==TRUE){
+         #     if (length(dim(data_tmp))==4){
+         #       data_tmp[,1,(line1-1)*indlen+line2,] <- data_without_external_gamma[,3,(line1-1)*indlen+line2,]
+         #       data_tmp[,2,(line1-1)*indlen+line2,] <- data_without_external_gamma[,4,(line1-1)*indlen+line2,]
+         #       data_tmp[,3,(line1-1)*indlen+line2,] <- data_without_external_gamma[,1,(line1-1)*indlen+line2,]
+         #       data_tmp[,4,(line1-1)*indlen+line2,] <- data_without_external_gamma[,2,(line1-1)*indlen+line2,]
+         #       data_tmp[,5,(line1-1)*indlen+line2,] <- data_without_external_gamma[,7,(line1-1)*indlen+line2,]
+         #       data_tmp[,6,(line1-1)*indlen+line2,] <- data_without_external_gamma[,8,(line1-1)*indlen+line2,]
+          #      data_tmp[,7,(line1-1)*indlen+line2,] <- data_without_external_gamma[,5,(line1-1)*indlen+line2,]
+         #       data_tmp[,8,(line1-1)*indlen+line2,] <- data_without_external_gamma[,6,(line1-1)*indlen+line2,]
+         #       data_tmp[,9,(line1-1)*indlen+line2,] <-  data_without_external_gamma[,11,(line1-1)*indlen+line2,]
+         #       data_tmp[,10,(line1-1)*indlen+line2,] <- data_without_external_gamma[,12,(line1-1)*indlen+line2,]
+         #       data_tmp[,11,(line1-1)*indlen+line2,] <- data_without_external_gamma[,9,(line1-1)*indlen+line2, ]
+         #       data_tmp[,12,(line1-1)*indlen+line2,] <- data_without_external_gamma[,10,(line1-1)*indlen+line2,]
+         #       data_tmp[,13,(line1-1)*indlen+line2,] <- data_without_external_gamma[,15,(line1-1)*indlen+line2,]
+         #       data_tmp[,14,(line1-1)*indlen+line2,] <- data_without_external_gamma[,16,(line1-1)*indlen+line2,]
+         #       data_tmp[,15,(line1-1)*indlen+line2,] <- data_without_external_gamma[,13,(line1-1)*indlen+line2,]
+         #       data_tmp[,16,(line1-1)*indlen+line2,] <- data_without_external_gamma[,14,(line1-1)*indlen+line2,]
+         #     }
+         #     else if (length(dim(data_tmp))==3){
+          #      data_tmp[,1,] <- data_without_external_gamma[,3,]
+          #      data_tmp[,2,] <- data_without_external_gamma[,4,]
+          #      data_tmp[,3,] <- data_without_external_gamma[,1,]
+          #      data_tmp[,4,] <- data_without_external_gamma[,2,]
+          #      data_tmp[,5,] <- data_without_external_gamma[,7,]
+          #      data_tmp[,6,] <- data_without_external_gamma[,8,]
+          #      data_tmp[,7,] <- data_without_external_gamma[,5,]
+          #      data_tmp[,8,] <- data_without_external_gamma[,6,]
+           #     data_tmp[,9,] <-  data_without_external_gamma[,11,]
+           #     data_tmp[,10,] <- data_without_external_gamma[,12,]
+           #     data_tmp[,11,] <- data_without_external_gamma[,9, ]
+           #     data_tmp[,12,] <- data_without_external_gamma[,10,]
+           #     data_tmp[,13,] <- data_without_external_gamma[,15,]
+           #     data_tmp[,14,] <- data_without_external_gamma[,16,]
+           #     data_tmp[,15,] <- data_without_external_gamma[,13,]
+           #     data_tmp[,16,] <- data_without_external_gamma[,14,]
+#
+ #             }
+  #          }
+   #         if (length(dim(data_tmp))==4){
+    #          data_corr[,,(line1-1)*indlen+line2,] <- data_tmp[,,(line1-1)*indlen+line2,]
+     ##       }
+       #     else if (length(dim(data_tmp))==3){
+        #      data_corr <- data_tmp
+         #   }
+            ##if the gamma structure at the sink is either C or Cg4 we have to multiply the correlator in space space
+            ##by gamma5 from the right
+          #  if (sink_source_t_gamma5[line1]==TRUE){
+       #      if (length(dim(data_tmp))==4){
+       #       data_corr[,1,(line1-1)*indlen+line2,] <- data_tmp[,9,(line1-1)*indlen+line2,]
+       #       data_corr[,2,(line1-1)*indlen+line2,] <- data_tmp[,10,(line1-1)*indlen+line2,]
+       #       data_corr[,3,(line1-1)*indlen+line2,] <- data_tmp[,11,(line1-1)*indlen+line2,]
+       #       data_corr[,4,(line1-1)*indlen+line2,] <- data_tmp[,12,(line1-1)*indlen+line2,]
+       #       data_corr[,5,(line1-1)*indlen+line2,] <- data_tmp[,13,(line1-1)*indlen+line2,]
+       #       data_corr[,6,(line1-1)*indlen+line2,] <- data_tmp[,14,(line1-1)*indlen+line2,]
+       #       data_corr[,7,(line1-1)*indlen+line2,] <- data_tmp[,15,(line1-1)*indlen+line2,]
+       #       data_corr[,8,(line1-1)*indlen+line2,] <- data_tmp[,16,(line1-1)*indlen+line2,]
+       #       data_corr[,9,(line1-1)*indlen+line2,] <-  data_tmp[,1,(line1-1)*indlen+line2,]
+       #       data_corr[,10,(line1-1)*indlen+line2,] <- data_tmp[,2,(line1-1)*indlen+line2,]
+       #       data_corr[,11,(line1-1)*indlen+line2,] <- data_tmp[,3,(line1-1)*indlen+line2,]
+       #       data_corr[,12,(line1-1)*indlen+line2,] <- data_tmp[,4,(line1-1)*indlen+line2,]
+       #       data_corr[,13,(line1-1)*indlen+line2,] <- data_tmp[,5,(line1-1)*indlen+line2,]
+       #       data_corr[,14,(line1-1)*indlen+line2,] <- data_tmp[,6,(line1-1)*indlen+line2,]
+        #      data_corr[,15,(line1-1)*indlen+line2,] <- data_tmp[,7,(line1-1)*indlen+line2,]
+         #     data_corr[,16,(line1-1)*indlen+line2,] <- data_tmp[,8,(line1-1)*indlen+line2,]
+         #    }
+         #    else if (length(dim(data_tmp))==3){
+         #     data_corr[,1,1,] <- data_tmp[,9,]
+         #     data_corr[,1,2,] <- data_tmp[,10,]
+         #     data_corr[,1,3,] <- data_tmp[,11,]
+         #     data_corr[,1,4,] <- data_tmp[,12,]
+         #     data_corr[,1,5,] <- data_tmp[,13,]
+         #     data_corr[,1,6,] <- data_tmp[,14,]
+         #     data_corr[,1,7,] <- data_tmp[,15,]
+         #     data_corr[,1,8,] <- data_tmp[,16,]
+         #     data_corr[,1,9,] <-  data_tmp[,1,]
+         #     data_corr[,1,10,] <- data_tmp[,2,]
+         #     data_corr[,1,11,] <- data_tmp[,3,]
+         #     data_corr[,1,12,] <- data_tmp[,4,]
+         #     data_corr[,1,13,] <- data_tmp[,5,]
+         #     data_corr[,1,14,] <- data_tmp[,6,]
+         #     data_corr[,1,15,] <- data_tmp[,7,]
+         #     data_corr[,1,16,] <- data_tmp[,8,]
+      #       }
+      #       
+      #      }
+      ##    }
+      #  }
+
+        #Doing possible spin projection
+       # data_projection <- data_corr
+        #Doing the symmetrization
+        #C+-B=<0.25*Tr(1+-g4)*C>
+        #Csymmetried=C+B(t)-C-B(T-t)
+
+
+        #partial_array<- array(0,dim=c(2,Time/2+1,length(final_indices)))
+
+        #Obtaining C(T-t)
+
+        #if we filter out just one gammastructure and momentum, then
+        #data_projection will be not 4 but three dimensional structure
+        #therefore we have do an if statement here
+        
+       # if (length(dim(data_corr))==3){
+       #   reversed <- data_projection[,,Time:1]
+#
+ #         partial_array[,1,] <- 0.25*(data_projection[,1,1]+data_projection[,6,1])
+ #         for (line in 2:(Time/2+1)){
+ #           partial_array[,line,] <- 0.25*(data_projection[,1,line]+data_projection[,6,line])
+ #           partial_array[,line,] <- partial_array[,line,]-0.25*(reversed[,11,line-1]+reversed[,16,line-1])
+ #         }
+#
+#        }
+#        else{
+#          reversed <- data_projection[,,,Time:1]
+#
+#          partial_array[,1,] <- 0.25*(data_projection[,1,,1]+data_projection[,6,,1])
+#          for (line in 2:(Time/2+1)){
+#            partial_array[,line,] <- 0.25*(data_projection[,1,,line]+data_projection[,6,,line])
+#            partial_array[,line,] <- partial_array[,line,]-0.25*(reversed[,11,,line-1]+reversed[,16,,line-1])
+#
+#          }
+ #       }
+
+        #convert the final array into two dimensional
+        #merge the dimension time, gamma and momentum
+        #leave the complex (real and imag) intact
+        #real part will be equal to cf
+        #complex part will be equal to icf
+  #      tmpprod <- prod(dim(partial_array))
+   #     dim(partial_array) <- c(2, tmpprod/2)
+   #     data_final_real <- rbind(data_final_real, partial_array[1,])
+   #     data_final_imag <- rbind(data_final_imag, partial_array[2,])
+#
+ #     }
+ #   }
+
+  #  rhdf5::H5Fclose(h5f)
+  #}
+  #attr(data_final_real,"dimnames") <- NULL
+  #attr(data_final_imag,"dimnames") <- NULL
+
+  #cf <- cf_meta(nrObs=length(final_indices)*length(selections[[avail_baryon_types]]$interp),Time=Time,nrStypes=1,symmetrised=TRUE)
+  #cf$symmetrised=TRUE
+  #cf <- cf_orig(cf, cf=data_final_real,icf=data_final_imag)
+  #return(invisible(cf))
 
+#}
diff --git a/R/gevp.R b/R/gevp.R
index e7ed455d..8fb506e3 100644
--- a/R/gevp.R
+++ b/R/gevp.R
@@ -122,13 +122,17 @@ gevp <- function(cf, Time, t0 = 1, element.order = 1:cf$nrObs,
   cM <- cM + t(cM)
   ## check for positive definiteness
   ev.cM <- eigen(cM, symmetric=TRUE, only.values = TRUE)
+  symmetric_problem <- TRUE
   if(any(ev.cM$values < 0)) {
-    stop("gevp: matrix at t0 is not positive definite. Aborting...\n")
+    #stop("gevp: matrix at t0 is not positive definite. Aborting...\n")
+    symmetric_problem <- FALSE
+    cM_orig <- cM
+  }
+  else {
+    ## compute Cholesky factorisation
+    L <- chol(cM)
+    invL <- solve(L)
   }
-  ## compute Cholesky factorisation
-  L <- chol(cM)
-  invL <- solve(L)
-  
   ## now the time dependence for t != t0
   ## we need to multiply from the left with t(invL) and from the right with invL
   for(t in c((t0 + 1):(Thalf), (t0-1):0)) {
@@ -150,14 +154,21 @@ gevp <- function(cf, Time, t0 = 1, element.order = 1:cf$nrObs,
     }
     else {
       ## determine eigenvalues and vectors
-      
-      variational.solve <- eigen(t(invL) %*% cM %*% invL,
+      if (symmetric_problem == TRUE){
+         variational.solve <- eigen(t(invL) %*% cM %*% invL,
                                  symmetric=TRUE, only.values = FALSE, EISPACK=FALSE)
+      }
+      else {
+	 variational.solve <- geigen(cM, cM_orig,
+                                 symmetric=FALSE, only.values = FALSE)
+	 variational.solve$values <- sqrt(Re(variational.solve$values)**2+Im(variational.solve$values)**2)
+	 print(variational.solve$values)
+      }
       ## sort depending on input by values or vectors
       sortindex <- integer(matrix.size)
       decreasing <- (t >= t0)
       if(sort.type == "values" || t == t0+1) {
-        sortindex <- order(variational.solve$values, decreasing=decreasing)
+        sortindex <- order(Re(variational.solve$values), decreasing=decreasing)
       }
       else if(sort.type == "vectors") {
         ## compute the scalar product of eigenvectors with those at t.sort
@@ -193,18 +204,24 @@ gevp <- function(cf, Time, t0 = 1, element.order = 1:cf$nrObs,
         sortindex <- Perms[which.max(DMp),]
         ##if(!decreasing) sortindex <- rev(sortindex)
       }
+      print(sortindex)
       evalues[t+1,] <- variational.solve$values[sortindex]
       evectors[t+1,,] <- variational.solve$vectors[, sortindex]
     }
   }
   for(t in c((0:(t0-1)), (t0 + 1):(Thalf))) {
-    evectors[t+1,,] <- invL %*% evectors[t+1,,]
+    if (symmetric_problem == TRUE){
+      evectors[t+1,,] <- invL %*% evectors[t+1,,]
+    }
     tmp <- matrix(Cor[ii+t], nrow=matrix.size, ncol=matrix.size) %*% evectors[t+1,,]
     ## t(evectors[t+1,,]) %*% tmp must be proportional to delta_ij
     ## these are the amplitudes up to a factor sqrt(exp(-mt) \pm exp(-m(T-t)))
     ## diag(t(evectors[t+1,,]) %*% tmp) might get negative due to fluctuations
     ## we set them to NA first
     d <- diag(t(evectors[t+1,,]) %*% tmp)
+    if (symmetric_problem == FALSE){
+      d <- Re(d)
+    }
     d[d < 0] <- NA
     amplitudes[t+1,,] <- t(t(tmp)/sqrt(d))
     rm(tmp)
@@ -218,6 +235,7 @@ gevp <- function(cf, Time, t0 = 1, element.order = 1:cf$nrObs,
   }
   
   return(invisible(list(evalues=evalues, evectors=evectors, amplitudes=amplitudes)))
+  
 }
 
 
diff --git a/R/gl_fse_corr.R b/R/gl_fse_corr.R
new file mode 100644
index 00000000..18f2f3a3
--- /dev/null
+++ b/R/gl_fse_corr.R
@@ -0,0 +1,18 @@
+#' apply parameter-free Gasser-Leutwyler FSE corrections
+#' @details Apply FSE correction of a meson mass and decay constant
+#' based on eqs. 8 and 9 of Colangeo, Duerr, Haefeli [Nucl.Phys.B 721 (2005)]
+#' going back to Gasser and Leutwyler. Note that we employ the normalisation
+#' in which f_pi = 130.4.
+#' @param xi Numeric vector, the ratio \eqn{M_{ps}^2 / (4 \pi f_\pi)^2}. 
+#' @param lambda Numeric, the dimensionless product \eqn{M_{ps} \cdot L}. 
+#' @return Data frame with two columns containing the correction factors \code{K_mps}
+#' and \code{K_fps} for which we have:
+#' \deqn{ M_{ps}(L) = M_{ps}(\infty) \cdot K_{M_{ps} + \mathcal{O}(M_{ps} \xi^2) }
+#' \deqn{ f_{ps}(L) = f_{ps}(\infty) \cdot K_{f_{ps} + \mathcal{O}(f_{ps} \xi^2) }
+#' @export
+gl_fse_corr <- function(xi, lambda){
+  stopifnot( length(xi) == length(lambda) )
+  g1_tilde <- g1(lambda)
+  return(data.frame(K_Mps = (1+0.5*xi*g1_tilde), K_fps = (1-2*xi*g1_tilde)))
+}
+
diff --git a/R/matrixfit.R b/R/matrixfit.R
index cac0612d..e189fae7 100644
--- a/R/matrixfit.R
+++ b/R/matrixfit.R
@@ -359,9 +359,9 @@ matrixfit <- function(cf, t1, t2,
   if(pcmodel) {
     stopifnot(inherits(cf, 'cf_principal_correlator'))
   }
-  if(cf$symmetrised == FALSE){
-    stop('You must symmetrize and bootstrap the function before fitting.')
-  }
+  #if(cf$symmetrised == FALSE){
+  #  stop('You must symmetrize and bootstrap the function before fitting.')
+  #}
 
   t1p1 <- t1 + 1
   t2p1 <- t2 + 1
@@ -443,7 +443,8 @@ matrixfit <- function(cf, t1, t2,
       stop("not all parameters are used in the fit! Aborting\n")
     }
   }
-  
+  print(mSize)
+  print(dim(parlist)[2]) 
   if(dim(parlist)[2] != mSize) {
     warning(mSize, " ", dim(parlist)[2], "\n")
     stop("parlist has not the correct length! Aborting! Use e.g. extractSingleCor.cf or c to bring cf to correct number of observables\n")
diff --git a/R/new_matrixfit.R b/R/new_matrixfit.R
index 9555cb97..fd5e6a27 100644
--- a/R/new_matrixfit.R
+++ b/R/new_matrixfit.R
@@ -76,6 +76,53 @@ MatrixModel <- R6::R6Class(
   )
 )
 
+SingleBaryonModel <- R6::R6Class(
+  'SingleBaryonModel',
+  inherit = MatrixModel,
+  public = list(
+    prediction = function (par, x, ...) {
+      parind <- make_parind(self$parlist, length(x), summands = 1)
+      sign_vec <- make_sign_vec(self$sym_vec, length(x), self$m_size)
+      ov_sign_vec <- make_ov_sign_vec(self$neg_vec, length(x), self$m_size)
+
+      ov_sign_vec * 0.5 * par[parind[, 1]] * par[parind[, 2]] *
+        (exp(-par[1] * x) )#+ sign_vec * exp(-par[1] * (self$time_extent - x)))
+    },
+    prediction_jacobian = function (par, x, ...) {
+      parind <- make_parind(self$parlist, length(x), summands = 1)
+      sign_vec <- make_sign_vec(self$sym_vec, length(x), self$m_size)
+      ov_sign_vec <- make_ov_sign_vec(self$neg_vec, length(x), self$m_size)
+
+      ## Derivative with respect to the mass, `par[1]`.
+      zp <- ov_sign_vec * 0.5 * par[parind[, 1]] * par[parind[, 2]] *
+        (-x * exp(-par[1] * x) )#-
+          # (self$time_extent-x) * sign_vec * exp(-par[1] * (self$time_extent-x)))
+      res <- zp
+      
+      ## Derivatives with respect to the amplitudes.
+      for (i in 2:length(par)) {
+        zp1 <- rep(0, length(zp))
+        j <- which(parind[, 1] == i)
+        zp1[j] <- ov_sign_vec * 0.5 * par[parind[j, 2]] *
+          (exp(-par[1] * x[j]) )#+ sign_vec[j] * exp(-par[1] * (self$time_extent-x[j])))
+        
+        zp2 <- rep(0, length(zp))
+        j <- which(parind[, 2] == i)
+        zp2[j] <- ov_sign_vec * 0.5 * par[parind[j, 1]] *
+          (exp(-par[1] * x[j]) )#+ sign_vec[j] * exp(-par[1] * (self$time_extent-x[j])))
+        
+        res <- cbind(res, zp1 + zp2)
+      }
+      stopifnot(nrow(res) == length(x))
+      stopifnot(ncol(res) == length(par))
+
+      dimnames(res) <- NULL
+      
+      return (res)
+    }
+  )
+)
+
 SingleModel <- R6::R6Class(
   'SingleModel',
   inherit = MatrixModel,
@@ -98,19 +145,19 @@ SingleModel <- R6::R6Class(
         (-x * exp(-par[1] * x) -
            (self$time_extent-x) * sign_vec * exp(-par[1] * (self$time_extent-x)))
       res <- zp
-      
+
       ## Derivatives with respect to the amplitudes.
       for (i in 2:length(par)) {
         zp1 <- rep(0, length(zp))
         j <- which(parind[, 1] == i)
         zp1[j] <- ov_sign_vec * 0.5 * par[parind[j, 2]] *
           (exp(-par[1] * x[j]) + sign_vec[j] * exp(-par[1] * (self$time_extent-x[j])))
-        
+
         zp2 <- rep(0, length(zp))
         j <- which(parind[, 2] == i)
         zp2[j] <- ov_sign_vec * 0.5 * par[parind[j, 1]] *
           (exp(-par[1] * x[j]) + sign_vec[j] * exp(-par[1] * (self$time_extent-x[j])))
-        
+
         res <- cbind(res, zp1 + zp2)
       }
 
@@ -118,12 +165,13 @@ SingleModel <- R6::R6Class(
       stopifnot(ncol(res) == length(par))
 
       dimnames(res) <- NULL
-      
+
       return (res)
     }
   )
 )
 
+
 TwoAmplitudesModel <- R6::R6Class(
   'TwoAmplitudesModel',
   inherit = MatrixModel,
@@ -367,12 +415,57 @@ NParticleModel <- R6::R6Class(
   )
 )
 
-## Model for a single correlator and a simple additive constant.
-##
-## @description
-## This model is just the “single” model plus a simple additive
-## constant. In some way it is a specialization of the “n particles” model as
-## the energy difference is just set to zero.
+#' Model for the n particle correlator fit
+#'
+#' @description
+#' n particle correlator with thermal pollution term(s)
+NBaryonParticleModel <- R6::R6Class(
+  'NBaryonParticleModel',
+  inherit = MatrixModel,
+  public = list(
+    initialize = function (time_extent, parlist, sym_vec, neg_vec, m_size) {
+      super$initialize(time_extent, parlist, sym_vec, neg_vec, m_size)
+      self$sm <- SingleBaryonModel$new(time_extent, parlist, sym_vec, neg_vec, m_size)
+     },
+    prediction = function (par, x, ...) {
+      n <- length(par) / 2
+      y <- 0
+      par_aux <- numeric(2)
+      for (i in 1:n) {
+        par_aux[1] <- par[2*i - 1]
+        par_aux[2] <- par[2*i]
+        corr_part <- self$sm$prediction(par_aux, x, ...)
+        y <- y + corr_part
+      }
+      return(y)
+    },
+    prediction_jacobian = function (par, x, ...) {
+      n <- length(par) / 2
+      par_aux <- numeric(2)
+      par_aux[1] <- par[1]
+      par_aux[2] <- par[2]
+      y <- self$sm$prediction_jacobian(par_aux, x, ...)
+      if (n > 1) {
+        for (i in 2:n) {
+          par_aux[1] <- par[2*i - 1]
+          par_aux[2] <- par[2*i]
+          jacobian_part_aux <- self$sm$prediction_jacobian(par_aux, x, ...)
+          y <- cbind(y, jacobian_part_aux)
+        }
+      }
+      return(y)
+    },
+    sm = NA
+  )
+)
+
+
+#' Model for a single correlator and a simple additive constant.
+#'
+#' @description
+#' This model is just the “single” model plus a simple additive
+#' constant. In some way it is a specialization of the “n particles” model as
+#' the energy difference is just set to zero.
 SingleConstantModel <- R6::R6Class(
   'SingleConstantModel',
   inherit = MatrixModel,
@@ -565,6 +658,8 @@ new_matrixfit <- function(cf,
   
   if (model == 'single') {
     model_object <- SingleModel$new(cf$Time, parlist, sym.vec, neg.vec, mSize)
+  } else if (model == 'single_baryon') {
+    model_object <- SingleBaryonModel$new(cf$Time, parlist, sym.vec, neg.vec, mSize)
   } else if (model == 'two_amplitudes') {
     stopifnot(cf$nrObs == 1)
     model_object <- TwoAmplitudesModel$new(cf$Time, parlist, sym.vec, neg.vec, mSize)
@@ -580,6 +675,9 @@ new_matrixfit <- function(cf,
   } else if (model == 'n_particles') {
     stopifnot(cf$nrObs == 1)
     model_object <- NParticleModel$new(cf$Time, parlist, sym.vec, neg.vec, mSize)
+  } else if (model == 'n_baryon_particles') {
+    stopifnot(cf$nrObs == 1)
+    model_object <- NBaryonParticleModel$new(cf$Time, parlist, sym.vec, neg.vec, mSize)
   } else if (model == 'single_constant') {
     stopifnot(cf$nrObs == 1)
     model_object <- SingleConstantModel$new(cf$Time, parlist, sym.vec, neg.vec, mSize)
diff --git a/R/plotutils.R b/R/plotutils.R
index 06fb9560..8fc962ae 100644
--- a/R/plotutils.R
+++ b/R/plotutils.R
@@ -19,6 +19,7 @@ compute.plotlims <- function(val, logscale, cumul.dval, cumul.mdval){
     tmp <- val+2*apply(X=cumul.mdval,MARGIN=1,FUN=min,na.rm=TRUE)
     tmpp <- val+2*apply(X=cumul.dval,MARGIN=1,FUN=max,na.rm=TRUE)
   }
+  # remove negative values if we're working on a log scale
   if(logscale) {
     tmp <- tmp[ tmp > 0 ]
     tmpp <- tmpp[ tmpp > 0 ]
diff --git a/R/readutils.R b/R/readutils.R
index 305a3080..b256da35 100644
--- a/R/readutils.R
+++ b/R/readutils.R
@@ -556,114 +556,115 @@ readoutputdata <- function(filename) {
   return(invisible(data))
 }
 
-
-
-#' Read correlator data from single file
-#' 
-#' Reads arbitrary number of samples for a complex correlation function from a
-#' text file.
-#' 
-#' 
-#' @param file filename of file to read from.
-#' @param Time time extent of the correlation function
-#' @param sym if \code{TRUE} average C(+t) and C(-t), otherwise C(+t) and
-#' -C(-t). Averaging can be switched off using the \code{symmetrise} option.
-#' @param skip number of lines to skip at beginning of file
-#' @param check.t if set to an integer value larger than zero the function will
-#' assume that in the corresponding column of the file the Euclidean time is
-#' counted and it will check whether the maximum in this column is identical to
-#' Time-1.
-#' @param ind.vector index vector of length 2 with the indices of real and
-#' imaginary values of correlator, respectivley.
-#' @param symmetrise if set to \code{TRUE}, the correlation function will be
-#' averaged for \code{t} and \code{Time-t}, with the sign depending on the value
-#' of \code{sym}. Note that currently the correlator with t-values larger than
-#' \code{Time/2} will be discarded.
-#' @param path the path to the files.
-#' @param autotruncate Boolean. Whether to autotruncate or not
-#' @param avg Integer. Average over successive number samples
-#' @param stride Integer. Read only subset of files with corresponding stride.
-#' @param Nmin Integer. Minimal number of measurements that must remain after
-#' sparsification and averaging. Default equals to 4.
-#' @return returns a list with two arrays \code{cf} and \code{icf} with real
-#' and imaginary parts of the correlator, and integers \code{Time},
-#' \code{nrStypes=1} and \code{nrObs=1}. Both of the arrays have dimension
-#' \code{c(N, (Time/2+1))}, where \code{N} is the number of measurements
-#' (gauges).  \code{Time} is the time extent, \code{nrStypes} the number of
-#' smearing levels and \code{nrObs} the number of operators, both of which are
-#' currently fixed to 1.
-#' @author Carsten Urbach, \email{curbach@@gmx.de}
-#' @seealso \code{\link{readcmidatafiles}}, \code{\link{readbinarydisc}},
-#' \code{\link{readcmidisc}}, \code{\link{readcmicor}},
-#' \code{\link{readbinarycf}}
-#' @keywords file
-#' @export readtextcf
-readtextcf <- function(file, Time=48, sym=TRUE, path="", skip=1, check.t=0, ind.vector=c(2,3), symmetrise=TRUE,
-                       stride=1, avg=1, Nmin=4, autotruncate=TRUE) {
-  stopifnot(!missing(file))
-  stopifnot(Time >= 1)
-
-  tmp <- read.table(paste(path, file, sep=""), skip=skip)
-  stopifnot(!((length(ind.vector) < 2) || (max(ind.vector) > length(tmp)) || (min(ind.vector) < 1)))
-  
-  if(check.t > 0 && max(tmp[[check.t]]) != Time-1) {
-    stop("Time in function call does not match the one in the file, aborting...\n")
-  }
-
-  if(length(tmp[[ind.vector[1]]]) %% Time != 0) {
-    stop("Time does not devide the number of rows in file, aborting... check value of paramter skip to readtextcf!\n")
+processcf <- function(dat, Lt, sym, ind.vector, symmetrise, sparsity, avg, Nmin, autotruncate, Nmax){
+  i1 <- c(2:(Lt/2))
+  i2 <- c(T:(Lt/2+2))
+  ii <- c(1:(Lt))
+  sign <- +1
+  if(!sym) sign <- -1
+
+  dat <- array(dat[[ ind.vector[1] ]] + 1i*dat[[ ind.vector[2] ]], dim=c(Lt, length(dat[[ ind.vector[1] ]])/Lt))
+
+  ## remove unwanted measurements
+  if( Nmax != -1 ){
+    if( ncol(dat) > Nmax ){
+      cat(sprintf("[processcf] Reducing the number of total measurements to %d from %d\n",Nmax,ncol(dat)))
+      dat <- dat[,1:Nmax]
+    }
   }
-  
-  ii <- c(1:(Time/2+1))
 
-  tmp <- array(tmp[[ind.vector[1]]] + 1i*tmp[[ind.vector[2]]], dim=c(Time, length(tmp[[ind.vector[1]]])/Time))
-  if( (stride > 1 | avg > 1) & (ncol(tmp) %% (stride*avg) != 0) ){
+  ## make sure that the number of measurements is compatible with any sparsification
+  ## and/or averaging requested
+  if( (sparsity > 1 | avg > 1) & (ncol(dat) %% (sparsity*avg) != 0) ){
     if(autotruncate){
-      warning(sprintf("stride=%d, avg=%d, ncol=%d\n",stride,avg,ncol(tmp)))
-      warning("readtextcf: Sparsification and/or averaging requested, but their product does not divide the number of measurements!\n")
-      warning("readtextcf: Reducing the number of total measurements to fit!\n")
-      nmeas <- as.integer( (stride*avg)*floor( ncol(tmp)/(stride*avg) ))
-      if( nmeas/(stride*avg) >= Nmin ){
-        tmp <- tmp[,1:nmeas]
+      cat(sprintf("[processcf] sparsity=%d, avg=%d, ncol=%d\n",sparsity,avg,ncol(dat)))
+      cat("[processcf] Sparsification and/or averaging requested, but their product does not divide the number of measurements!\n")
+      cat("[processcf] Reducing the number of total measurements to fit!\n")
+      nmeas <- as.integer( (sparsity*avg)*floor( ncol(dat)/(sparsity*avg) ))
+      if( nmeas/(sparsity*avg) >= Nmin ){
+        dat <- dat[,1:nmeas]
       } else {
-        warning(sprintf("readtextcf: After sparsification and averaging, less than %d measurements remain, disabling sparsification and averaging!\n",Nmin))
-        stride <- 1
+        cat(sprintf("[processcf] After sparsification and averaging, less than %d measurements remain, disabling sparsification and averaging!\n",Nmin))
+        sparsity <- 1
         avg <- 1
       }
     } else {
-      stop("readtextcf: Sparsification and/or averaging requested, but their product does not divide the number of measurements!\n")
+      stop("[processcf] Sparsification and/or averaging requested, but their product does not divide the number of measurements!\n")
     }
   }
 
-  ## sparsify data
-  if(stride > 1){
-    sp.idx <- seq(from=1,to=ncol(tmp),by=stride)
-    tmp <- tmp[,sp.idx]
+  ## sparsify dat
+  if(sparsity > 1){
+    sp.idx <- seq(from=1,to=ncol(dat),by=sparsity)
+    dat <- dat[,sp.idx]
   } 
   # average over 'avg' measurements sequentially
   if(avg > 1){
-    tmp2 <- tmp
-    tmp <- array(0, dim=c(Time,ncol(tmp2)/avg))
-    for( i in c(1:ncol(tmp)) ){
-      tmp[,i] <- (1.0/avg)*apply(X=tmp2[,((i-1)*avg+1):(i*avg)],
-                                 MARGIN=1,
-                                 FUN=sum)
+    dat2 <- dat
+    dat <- array(0, dim=c(Lt,ncol(dat2)/avg))
+    for( i in c(1:ncol(dat)) ){
+      dat[,i] <- (1.0/avg)*apply(X=dat2[,((i-1)*avg+1):(i*avg)],
+                                  MARGIN=1,
+                                  FUN=sum)
     }
   }
-  
-  ret <- cf_meta(nrObs = 1, Time=Time, nrStypes = 1)
-  ret <- cf_orig(ret, cf = t(Re(tmp)), icf = t(Im(tmp)))
+
+  ret <- cf_meta(nrObs = 1, Time=Lt, nrStypes = 1)
+  ret <- cf_orig(ret, cf = t(Re(dat[ii,])), icf = t(Im(dat[ii,])))
 
   if (symmetrise) {
     sign <- +1
     if (!sym) sign <- -1
-
+    ## average +-t
     ret <- symmetrise.cf(ret, sign)
   }
-
-  return (invisible(ret))
+  
+  return(invisible(ret))
 }
+#' @title reading reweighting factors for a list of gauge configuration 
+#'        and random samples from ASCII files
+#' @param file_names_to_read list of filenames for the reweighting factors
+#' @param gauge_conf_list <- a list of integers with the indices of the gauge configs
+#' @param nsamples number of stochastic samples used for computing the reweighting factors
+read.rw <- function( file_names_to_read, gauge_conf_list, nsamples, monomial_id ) 
+{
+  stopifnot(length(gauge_conf_list)==length(file_names_to_read)) 
+  ret <- rw_meta(conf.index=gauge_conf_list)
+  tmp <- readcmidatafiles(files=file_names_to_read,skip=0,verbose=TRUE,colClasses=c("integer","integer","numeric","numeric","numeric","numeric","numeric"))
+  names(tmp)[1] <- "monomialid"
+  names(tmp)[2] <- "stochastic_index"
+  names(tmp)[3] <- "kappa_target"
+  names(tmp)[4] <- "kappa_original"
+  names(tmp)[5] <- "light_quark_mass_target"
+  names(tmp)[6] <- "light_quark_mass_original"
+  names(tmp)[7] <- "reweightingfactor"
+
+# Select the reweighting factor for a particular monomial
+
+  dplyr_avail <- requireNamespace("dplyr")
+  stopifnot(dplyr_avail)
+  
+  tmp <- dplyr::filter(tmp,monomialid==monomial_id)
+
+# Number of reweighted determinants for each gauge configuration
+
+  n_rew_factors <- length(tmp$reweightingfactor)/(nsamples*length(gauge_conf_list))
+  stopifnot(n_rew_factors == 1)
+  
 
+# Exponentianing and Averaging over the stochastic samples
+  
+  tmp2 <- matrix(tmp$reweightingfactor,nrow=nsamples,ncol=length(gauge_conf_list)*n_rew_factors)
+  tmp3 <- apply(exp(-tmp2),2,mean)
+
+# Normalize the largest reweighting factor to be one and storing this factor
+# this is neccessary due to the large value of the reweighting factor 
+# after exponentiating
+  tmp4 <- tmp3/max(tmp3)
+  
+  ret <- rw_orig(ret, rw = tmp4, conf.index=gauge_conf_list, max_value = max(tmp3))
+
+}
 #' @title reader for Nissa text format correlation functions
 #' @param file_basenames_to_read Character vector of file names without the
 #'                              smearing combination suffixes (such as 'll', 'ls', 'sl', 'ss')
@@ -786,6 +787,7 @@ readnissatextcf <- function(file_basenames_to_read,
 #' X$cf
 #' 
 #' @export readbinarycf
+
 readbinarycf <- function(files, 
                          Time, 
                          obs=5, 
@@ -1227,3 +1229,151 @@ readgradflow <- function(path, skip=0, basename="gradflow", col.names) {
   return(invisible(ldata))
 }
 
+
+read_bsm_temporal_phifield <- function(file, Lt, sym, path, scalars,
+                                       symmetrise=FALSE,
+                                       sparsity=1, avg=1, Nmin=4, autotruncate=TRUE,
+                                       Nmax=-1)
+{
+  if(missing(file)){
+    stop("files must be given!")
+  }
+  if(missing(Lt)){
+    stop("Lt time extent must be passed!")
+  }
+  if(missing(sym)){
+    stop("Symmetry must be passed!")
+  }
+  
+  cat(sprintf("[read_bsm_temporal_phifield] Attempting to read spatially averaged phi field %s\n",
+              file))
+
+  tmp <- read.table(file = sprintf("%s/%s", path, file),
+                    header = FALSE,
+                    skip = 0,
+                    col.names = c("t","re","sidx")
+                    )
+
+  # use only the requested scalars in the right order 
+  tmp <- tmp[ which( tmp$sidx %in% scalars ), ]
+  if( any( unique(tmp$sidx) != scalars ) ){
+    stop(sprintf("[read_bsm_temporal_phifield] It was not possible to extract the requested set of scalars from %s\n",file))
+  }
+  
+  tmp <- cbind(tmp,tmp[,3])
+  tmp[,3] <- 0
+
+  ret <- processcf(dat = tmp,
+                   Lt = Lt,
+                   sparsity = sparsity,
+                   avg = avg,
+                   Nmin = Nmin,
+                   Nmax = Nmax,
+                   sym = sym,
+                   symmetrise = symmetrise,
+                   ind.vector = c(2,3),
+                   autotruncate = autotruncate)
+
+  ret$scalars <- scalars
+
+  return(invisible(ret))
+}
+
+readbsmcf <- function(file, Lt, sym, path, nscalars, 
+                      symmetrise=FALSE, skip=1,
+                      sparsity=1, avg=1, Nmin=4, autotruncate=TRUE,
+                      Nmax=-1){
+  if(missing(file)){
+    stop("files must be given!")
+  }
+  if(missing(Lt)){
+    stop("Lt time extent must be passed!")
+  }
+  if(missing(sym)){
+    stop("Symmetry must be passed!")
+  }
+  if(missing(nscalars)){
+    stop("nscalars must be passed!")
+  }
+
+  tmp <- read.table(file = sprintf("%s/%s", path, file),
+                    header = FALSE,
+                    skip = skip,
+                    col.names = c("t","re","im","gidx","sidx")
+                    )
+
+  gidx <- unique(tmp$gidx)
+  sidx <- unique(tmp$sidx)
+
+  meas_gauges            <- length(unique(tmp$gidx))
+  meas_total             <- nrow(tmp)/Lt
+  meas_scalars_per_gauge <- meas_total / meas_gauges
+
+  cat(sprintf("[readbsmcf] Data has %d scalars per gauge, %d gauges and %d scalars per gauge were requested\n",
+              meas_scalars_per_gauge, meas_gauges, nscalars))
+  if ( meas_scalars_per_gauge > nscalars ){
+    if( meas_scalars_per_gauge %% nscalars != 0 ){
+      stop(sprintf("The requested number of scalars per gauge %d does not divide the actual %d in %s",
+                   nscalars,
+                   meas_scalars_per_gauge))
+    } else {
+      # we use every 'scalar_step' measurement to get the correct number of scalars
+      cat(sprintf("[readbsmcf] Adjusting the number of scalars per gauge from %d to %d as requested\n",meas_scalars_per_gauge,nscalars))
+      scalar_step <- meas_scalars_per_gauge / nscalars
+      sidx <- sidx[ seq(1, length(sidx), scalar_step) ]
+      tmp <- tmp[ which( tmp$sidx %in% sidx ), ]
+    }
+  } else if ( meas_scalars_per_gauge < nscalars ){
+    stop(sprintf("It seems that the number of scalars per gauge in %s does not match the requested %d",
+                 file,
+                 nscalars))
+  }
+  ret <- processcf(dat = tmp,
+                   Lt = Lt, 
+                   ind.vector = c(2,3), 
+                   sym = sym,
+                   symmetrise = symmetrise,
+                   sparsity = sparsity,
+                   avg = avg,
+                   Nmin = Nmin,
+                   autotruncate = autotruncate,
+                   Nmax = Nmax
+                   )
+  ret$scalars <- sidx[1:Nmax]
+  ret$gauges <- gidx[1:as.integer(Nmax/avg)]
+  return(invisible(ret))
+}
+
+readtextcf <- function(file, T=48, sym=TRUE, path="", skip=1, check.t=0, ind.vector=c(2,3), symmetrise=TRUE,
+                       sparsity=1, avg=1, Nmin=4, autotruncate=TRUE, Nmax=-1) {
+  if(missing(file)) {
+    stop("files must be given! Aborting...\n")
+  }
+  if(T < 1) {
+    stop("T must be larger than 0 and integer, aborting...\n")
+  }
+
+  tmp <- read.table(paste(path, file, sep=""), skip=skip)
+  if((length(ind.vector) < 2) || (max(ind.vector) > length(tmp)) || (min(ind.vector) < 1)){
+    stop("index vector too short or out of range\n")
+  }
+  
+  if(check.t > 0 && max(tmp[[check.t]]) != T-1) {
+    stop("T in function call does not match the one in the file, aborting...\n")
+  }
+
+  if(length(tmp[[ind.vector[1]]]) %% T != 0) {
+    stop("T does not devide the number of rows in file, aborting... check value of paramter skip to readtextcf!\n")
+  }
+  
+  return(invisible(processcf(dat = tmp,
+                             Lt = T,
+                             sym = sym,
+                             ind.vector = ind.vector,
+                             symmetrise = symmetrise,
+                             sparsity = sparsity,
+                             avg = avg,
+                             Nmin = Nmin,
+                             autotruncate = autotruncate,
+                             Nmax = Nmax)))
+}
diff --git a/R/removeTemporal.cf.R b/R/removeTemporal.cf.R
index 016fa335..9a081689 100644
--- a/R/removeTemporal.cf.R
+++ b/R/removeTemporal.cf.R
@@ -151,16 +151,34 @@ old_removeTemporal.cf <- function(cf,
   # hold and that no new fields have been added that we are not aware of.
   ret <- cf_meta(nrObs = cf$nrObs, Time = cf$Time, nrStypes = cf$nrStypes,
                  symmetrised = cf$symmetrised)
-  ret <- cf_orig(ret,
-                 cf = cf$cf)
-  ret <- cf_boot(ret,
+  if (has_icf(cf)){
+    ret <- cf_orig(ret,
+                 cf = cf$cf,
+                 icf= cf$icf)
+    ret <- cf_boot(ret,
                  boot.R = cf$boot.R,
                  boot.l = cf$boot.l,
                  seed = cf$seed,
                  sim = cf$sim,
                  endcorr = cf$endcorr,
                  cf.tsboot = cf$cf.tsboot,
+                 icf.tsboot = cf$icf.tsboot,
                  resampling_method = cf$resampling_method)
+
+  }
+  else{
+    ret <- cf_orig(ret,
+                 cf = cf$cf)  
+  
+    ret <- cf_boot(ret,
+                 boot.R = cf$boot.R,
+                 boot.l = cf$boot.l,
+                 seed = cf$seed,
+                 sim = cf$sim,
+                 endcorr = cf$endcorr,
+                 cf.tsboot = cf$cf.tsboot,
+                 resampling_method = cf$resampling_method)
+  }
   ret <- cf_shifted(ret,
                     deltat = cf$deltat,
                     forwardshift = cf$forwardshift)
@@ -227,6 +245,11 @@ takeTimeDiff.cf <- function (cf, deltat = 1, forwardshift = FALSE) {
   if( inherits(cf, 'cf_orig') ) {
     cf$cf[,tlhs] <- cf$cf[,trhs1]-cf$cf[,trhs2]
     cf$cf[,-tlhs] <- NA
+    print(has_icf(cf))
+    if (has_icf(cf)){
+      cf$icf[,tlhs] <- cf$icf[,trhs1]-cf$icf[,trhs2]
+      cf$icf[,-tlhs] <- NA
+    }
   }
 
   # now the bootstrap samples
@@ -235,16 +258,44 @@ takeTimeDiff.cf <- function (cf, deltat = 1, forwardshift = FALSE) {
     cf$cf.tsboot$t0[-tlhs] <- NA
     cf$cf.tsboot$t[,tlhs] <- cf$cf.tsboot$t[,trhs1]-cf$cf.tsboot$t[,trhs2]
     cf$cf.tsboot$t[,-tlhs] <- NA
+    if (has_icf(cf)){
+      cf$icf.tsboot$t0[tlhs] <- cf$icf.tsboot$t0[trhs1]-cf$icf.tsboot$t0[trhs2]
+      cf$icf.tsboot$t0[-tlhs] <- NA
+      cf$icf.tsboot$t[,tlhs] <- cf$icf.tsboot$t[,trhs1]-cf$icf.tsboot$t[,trhs2]
+      cf$icf.tsboot$t[,-tlhs] <- NA
+
+    }
   }
 
   # We perform a new construction in order to have only the fields defined that
   # we want and also to make sure that invariants are holding.
   ret <- cf_meta(nrObs = cf$nrObs, Time = cf$Time, nrStypes = cf$nrStypes,
                  symmetrised = cf$symmetrised)
-  ret <- cf_orig(ret,
+  if (has_icf(cf)){
+    ret <- cf_orig(ret,
+                 cf = cf$cf,
+                 icf= cf$icf)
+  }
+  else{
+    ret <- cf_orig(ret,
                  cf = cf$cf)
+  }
+  
   if (inherits(cf, 'cf_boot')) {
-    ret <- cf_boot(ret,
+    if (has_icf(cf)){
+      ret <- cf_boot(ret,
+                   boot.R = cf$boot.R,
+                   boot.l = cf$boot.l,
+                   seed = cf$seed,
+                   sim = cf$sim,
+                   endcorr = cf$endcorr,
+                   cf.tsboot = cf$cf.tsboot,
+                   icf.tsboot = cf$icf.tsboot,
+                   resampling_method = cf$resampling_method)
+
+    }
+    else {
+      ret <- cf_boot(ret,
                    boot.R = cf$boot.R,
                    boot.l = cf$boot.l,
                    seed = cf$seed,
@@ -252,6 +303,7 @@ takeTimeDiff.cf <- function (cf, deltat = 1, forwardshift = FALSE) {
                    endcorr = cf$endcorr,
                    cf.tsboot = cf$cf.tsboot,
                    resampling_method = cf$resampling_method)
+    }
   }
   ret <- cf_shifted(ret,
                     deltat = deltat,
diff --git a/R/rw.R b/R/rw.R
new file mode 100644
index 00000000..1be37888
--- /dev/null
+++ b/R/rw.R
@@ -0,0 +1,251 @@
+#' Reweighting factor container
+#'
+#' This function `rw()` creates containers for reweighting factors
+#' of class `rw`. This class is particularly designed to store reweighting 
+#' factors for each gauge configuration, that can be applied on 
+#' correlation functions emerging in statistical and quantum field theory
+#' simulations. Multiplication operation is defined for this class,
+#' as well as for increasing statistics and \link{is.rw}.
+#'
+#' @details
+#'#'
+#' @family rw constructors
+#'
+#' @export
+rw <- function () {
+  rw <- list()
+  class(rw) <- append(class(rw), 'rw')
+  return (rw)
+}
+
+#' rw metadata mixin constructor
+#'
+#' @param .rw `rw` object to extend.
+#' @param conf.index list of Integers, containing the index of gauge configurations
+#'
+#' @family rw constructors
+#'
+#' @export
+rw_meta <- function (.rw = rw(), conf.index ) {
+   stopifnot(inherits(.rw, 'rw'))
+
+   .rw$conf.index <- conf.index
+
+   class(.rw) <- append(class(.rw), 'rw_meta')
+   return (.rw)
+}
+
+
+#' @export
+rw_orig <- function (.rw = rw(), rw, conf.index, max_value) {
+  stopifnot(inherits(.rw, 'rw'))
+
+  .rw$rw <- rw
+  .rw$max_value <- max_value
+  .rw$conf.index <- conf.index
+
+  class(.rw) <- append(class(.rw), 'rw_orig')
+  class(.rw) <- append(class(.rw), 'rw_meta')
+
+  ## Norm is needed in order to be able to extend statistics
+  ## If the object does not have `rw_norm` it cannot be extended
+  class(.rw) <- append(class(.rw), 'rw_norm')
+  return (.rw)
+}
+
+#' @export
+rw_unit <- function (.rw = rw(), conf.index) {
+  stopifnot(inherits(.rw, 'rw'))
+
+  .rw$conf.index <- conf.index
+  .rw$rw <- rep(1.0,length(conf.index))
+  .rw$max_value <- 1.0
+
+  class(.rw) <- append(class(.rw), 'rw_orig')
+  class(.rw) <- append(class(.rw), 'rw_meta')
+  class(.rw) <- append(class(.rw), 'rw_norm')
+  return (.rw)
+}
+
+#' Checks whether the cf object contains no data
+#'
+#' @param .rw `rw` object.
+#'
+#' @examples
+#' # The empty rw object must be empty:
+#' is_empty.rw(rw())
+#'
+is_empty.rw <- function (.rw) {
+  setequal(class(.rw), class(rw())) &&
+    is.null(names(.rw))
+}
+
+
+#' Arithmetically multiplies two reweighting factors 
+#'
+#' @param rw1,rw2 `rw_orig` objects to be muplitplied
+#' @param nf1,nf2 Integer. Factors that determines the number of flavours we
+#'        have for reweighting factor 1 and 2, the default value is 1, because
+#'        usually we compute the reweighting factors for Q: the product of 
+#'        up and down or strange and charm determinant, there are no additional 
+#'        terms in the sea determinant that have to be reweighted with the same
+#'        rw factor. In case we compute the rw1 factor only for the up quark, we 
+#'        for example, we have to set nf1=2 to obtain the rw factor
+#'        for the light determinant.  
+#'
+#' @return
+#' The value is
+#' \deqn{rw_1*nf1*rw_2*nf2) \,.}
+#'
+#' @export
+multiply.rw <- function(rw1, rw2, nf1 = 1, nf2 = 1) {
+
+  stopifnot(inherits(rw1, 'rw'))
+  stopifnot(inherits(rw2, 'rw'))
+ 
+  stopifnot(inherits(rw1, 'rw_orig'))
+  stopifnot(inherits(rw2, 'rw_orig'))
+
+  stopifnot(inherits(rw1, 'rw_meta'))
+  stopifnot(inherits(rw2, 'rw_meta'))
+ 
+  stopifnot(identical( rw1$conf.index, rw2$conf.index ))
+
+  rw <- rw()
+  rw$rw <- nf1*rw1$rw * nf2*rw2$rw
+  rw$max_value <- NA
+  rw$conf.index <- rw1$conf.index 
+
+  class(rw) <- append(class(rw), 'rw_orig')
+  class(rw) <- append(class(rw), 'rw_meta')
+  
+  return(rw)
+}
+
+#' Arithmetically multiply reweighting factors
+#'
+#' @param rw1,rw2 `rw_orig` objects.
+#'
+#' @export
+'*.rw' <- function (rw1, rw2) {
+  multiply.rw(cf1, cf2, nf1 = 1, nf2 = 1)
+}
+
+
+#' Checks whether an object is a rw
+#'
+#' @param x Object, possibly of class `rw`.
+#'
+#' @export
+is.rw <- function (x) {
+  inherits(x, "rw")
+}
+
+meanval <- function(d, w) {return(mean(d$rw[w]))}
+
+#' Plot a reweighting factor function
+#'
+#' @param x `rw` object
+#' @param rep See \code{\link{plotwitherror}}.
+#'
+#' @inheritParams plotwitherror
+#'
+#' @export
+plot.rw <- function(x, rep = FALSE, ...) {
+  rw <- x
+  stopifnot(inherits(rw, 'rw'))
+  negs <- which(rw$rw < 0)
+  neg.vec <- rep(1,times=length(rw$rw))
+  neg.vec[negs] <- -1
+ 
+  colneg <- rep("black",times=length(rw$rw))
+  colneg[negs] <- "red"
+
+# Estimating the statistical error on the reweighting factor
+
+  rw_boot <- boot(rw, meanval, R=1500)
+  rw_error <- sd(rw_boot)
+
+  val <- rw$rw/mean(rw$rw)
+  err <- rep(x=rw_error,length(rw$rw))
+
+  df <- data.frame(t = c(1:length(val)),
+                   CF = val,
+                   Err = err)
+
+  plotwitherror(x = c(1:length(val)), y = neg.vec * df$CF, dy = df$Err, col=colneg, rep = rep, ...)
+
+  return(invisible(df))
+}
+
+#' Combine reweighting factors from different replicas
+#' 
+#' @param rw1 `rw` object: reweighting factor for replicum A
+#' @param rw2 `cf` object: reweighting factor for replicum B
+#' @param reverse1 `boolean` After the bifurcation point one of
+#'                           the replicas (chain of reweighting 
+#'                           factors in simulation time) has  
+#'                           to be reversed.
+#' @param reverse2 `boolean`
+#'
+#' @examples
+#' # Suppose we have reweighting factors in replicum A from 0 to 500
+#' # in steps of 4 and in replicum B from 4 to 500 in steps of 4.
+#' # To combined the two replicas we have to use
+#'
+#' #addStat.rw(rw_replicumB, rw_replicumA, TRUE, FALSE)
+#'
+#' # which means
+#' # combined=(rw500 from B, rw496 from B,...,rw004 from B, rw000 from A, ..
+#' # rw500 from A) 
+#' @export
+addStat.rw <- function(rw1, rw2,reverse1=FALSE, reverse2=FALSE) {
+  stopifnot(inherits(rw1, 'rw'))
+  stopifnot(inherits(rw2, 'rw'))
+
+  if (is_empty.rw(rw1)) {
+    return (invisible(rw2))
+  }
+  if (is_empty.cf(rw2)) {
+    return (invisible(rw1))
+  }
+
+  stopifnot(inherits(rw1, 'rw_orig'))
+  stopifnot(inherits(rw2, 'rw_orig'))
+
+
+  stopifnot(inherits(rw1, 'rw_meta'))
+  stopifnot(inherits(rw2, 'rw_meta'))
+ 
+  stopifnot(inherits(rw1, 'rw_norm'))
+  stopifnot(inherits(rw2, 'rw_norm'))
+
+  rw    <- rw1
+
+  if (reverse1 == TRUE){
+    rw_1   <- rev(rw1$rw)
+    conf_1 <- rev(rw1$conf.index)
+  }
+  else{
+    rw_1   <- rw1$rw
+    conf_1 <- rw1$conf.index 
+  }
+  if (reverse2 == TRUE){
+    rw_2   <- rev(rw2$rw)
+    conf_2 <- rev(rw2$conf.index)
+  }
+  else{
+    rw_2   <- rw2$rw
+    conf_2 <- rw2$conf.index 
+  }
+  rw_1 <- rw_1*rw1$max_value
+  rw_2 <- rw_2*rw2$max_value
+  
+  rw$rw <- c(rw_1, rw_2)
+  max_value <- max(rw$rw)
+  rw$rw <- rw$rw/max_value
+  rw$conf.index <- c(conf_1, conf_2) 
+  rw$max_value  <- max_value 
+
+  return (invisible(rw))
+}