.packageName <- "mscalib"
#Copyright 2004, W. Wolski, all rights reserved.
##not neccessary - its the same for each peaklist
##



calibspline <- function(ispl,mv,error,theo)
  {
    ##t Constructor
    ##- Returns object of class calibspline. It is used by the function \code{getextcalib.massvector}.
    ##n Used by by method getextcalib.
    ##+ ispl : object of class "smooth.spline"
    ##+ mv : massvector
    ##+ error : vector with errors.
    ##+ theo : vector with theoretical masses
    ##sa getextcalib.massvector
    ##ex data(ppg)
    ##ex getextcalib(ppg)
    ##ex if(!inherits(ppg,"calibspline")){stop("its not a calibspline")}
    ##ex plot(ppg)
    
    if(!inherits(ispl,"smooth.spline") || !(inherits(mv,"massvectorlist")||inherits(mv,"massvector") ))
      {
        print(class(ispl))
        print(class(mv))
        stop("param ispl not of class smooth.spline or mv not of class massvector!\n")
      }
    ispl$allow <- c("info","error","theo")
    ispl$info <- info(mv)
    if(length(error) != length(theo))
      stop("arg error and arg theo must have the same length\n")
    
    ispl$error <- error
    ispl$theo <- theo
    class(ispl)<-c("calibspline","smooth.spline","mylistobj")
    return(ispl)
  }



print.calibspline <- function(x,...)
  {
    ##t Print Calibspline Object
    ##- `print' prints its argument and returns it invisibly (via `invisible(x)')
    ##+ x : calibspline
    ##v info : the id of the calibspline
    ##v smooth.spline : smooth.spline
    ##e data(ppg)
    ##e tmp <- getextcalib(ppg,error=200)
    ##e print(tmp)
    
    cat("info :",x$info,"\n")
    res <- NextMethod("print")
    invisible(x)
  }

summary.calibspline <- function(object,...)
  {
    ##t Calibspline Summaries
    ##- Generates a summary, min, max etc.
    ##+ object : massvector
    ##sa summary
    ##e data(ppg)
    ##e tmp <- getextcalib(ppg,error=200)
    ##e summary(tmp)
    NextMethod(object)
  }


##MASSVECTOR

plot.calibspline <- function(x,...)
  {
    ##t Calibspline Plotting
    ##- Function for plotting object of class calibspline.
    ##+ x : calibspline
    ##+ ... : graphical parameters can be given as arguments to `plot'.
    ##e data(ppg)
    ##e tmp <- getextcalib(ppg,error=200)
    ##e plot(tmp)
    cS<-x
    x<-NULL
    plot(cS$theo,cS$error,pch="*",xlab="thoretical mass [m/z]",ylab="error [ppm]",main=mget(cS,"info"),...)
    lines(predict(cS,cS$x),col=2)
  }



applycalib.calibspline<-function(object,mv,...)
  {
    ##t External Calbiration
    ##- Applys object of class calibspline to massvector or massvectorlist to correct for measurment errors.
    ##d In case of \bold{external calibration} some sample spots are only dedicated
    ##d to calibration. Calibration samples which produces equidistant
    ##d peaks, which exact masses are known, can be used to precisely
    ##d estimate the mass dependent error function.
    ##+ object : massvector
    ##+ cS : calibspline
    ##v massvector : calibrated massvector
    ##sa applyextcalib.massvectorlist, getextcalib.massvector, getextcalib.massvectorlist
    ##r Gobom J, Mueller M, Egelhofer V, Theiss D, Lehrach H, Nordhoff E, 2002. A calibration method that simplifies and improves accurate determination of peptide molecular masses by MALDI-TOF MS. Anal Chem. 74(15):3915-23.
    ##r Wolski http://www.molgen.mpg.de/~wolski/mscalib
    ##e data(mv1)
    ##e data(ppg)
    ##e res<- getextcalib(ppg,getPPGmasses(),error=150)
    ##e plot(res)
    ##e mv2<-applycalib(res,mv1)
    ##e compare(mv1,mv2,error=300)
    ##e rm(mv1,mv2)
    ##e data(mvl)
    ##e mvl<-mvl[1:100]
    ##e res<-applycalib(res,mvl)
    
    if(!(inherits(mv,"massvectorlist")||inherits(mv,"massvector")))
      {
        stop(as.character(substitute(mv))," should be an object of class massvector or massvectorlist but is of class : ",class(mv),"\n")
      }
                                        #peaklist- array of peaks masses.
                                        #spline - spline to predict the error.
    if(inherits(mv,"massvector"))
      {
        error <- predict(object,mass(mv))
        error <- error$y
        masspred <- mv[,1]/(1+error/1e6)
        mv[,1] <- masspred
        return(mv)
      }
    if(inherits(mv,"massvectorlist"))
      {
        res<-lapply(mv,applyextcalib,object)
        mv <- massvectorlist(experiment(mv),res,project(mv))
        rm(res)
        return(mv)
      }
  }



#
#caliblist
#
#
#calibextlist  function(experiment,calib,...)
#  {
#    
#    res<-list()
#    tmp<-list(...)
#    allow<-c("experiment","project","data","calib")
#    attr(res,"allow") <- allow
#    if(!missing(experiment))
#      attr(res,"experiment") <- experiment
#    if(!missing(calib))
#      attr(res,"calib") <- calib
#    for(x in names(tmp))
#      {
#        if(x %in% allow)
#          {
#            attr(res,x)<-tmp[[x]]
#          }
#      }
#    class(res) <- c("calibextlist" , "caliblist" , "myobj" , "list")
#    res
#  }

#Copyright 2004, W. Wolski, all rights reserved.
############################################
##calibstat
#
summary.calibstat<-function(object,...)
  {
    ##t Calibstat summary
    ##- Generates Summary.
    ##+ object : object of class calibstat.
    ##e data(mv1)
    ##e data(cal)
    ##e test<-getintcalib(mv1,cal)
    ##e summary(test)
    cat("info :",info(object),"\n")
    print(as.vector(object))
    invisible(c(info(object),as.vector(object)))
  }

hist.calibstat <- function(x,...)
  {
    ##t Histogram
    ##- Histogram
    ##+ x : object of class calibstat
    warning("Not implemented!!\n")
  }
print.calibstat <- function(x,...)
  {
    ##t print
    ##- print
    ##+ x : object of class calibstat
    warning("Not implemented!!\n")
  }
plot.calibstat <- function(x,...)
  {
    ##t plot
    ##- plot
    ##+ x : object of class calibstat
    warning("Not implemented!!\n")
  }


################################################
#caliblist
#


plot.caliblist <- function(x,...)
  {
    ##t plot
    ##- plot
    ##+ x : object of class caliblist
    warning("not implemented yet!!\n")
  }





print.caliblist <- function(x,...)
  {
    ##t Print caliblist
    ##- `print' prints its argument and returns it invisibly (via `invisible(x)')
    ##+ x : caliblist
    ##v list : list
    ##e data(mvl)
    ##e mvl<-mvl[1:100]
    ##e data(cal)
    ##e test <- getintcalib(mvl,cal,error=500)
    ##e print(test)
    cl<-class(x)[1]
    cat("Class  :", cl ,"\n")
    exp <- mget(x,"experiment")
    cat("Experiment :", exp ,"\n")
    pro <- mget(x,"project")
    cat("Project    :", pro ,"\n")
    ll<- length(x)
    cat("Calibstat List lenght : ",ll ,"\n")
     if(length(x)>0)
      {
        cmv<- class(x[[1]])[1]
        cat("Class calibstat object: ", cmv ,"\n")
        allow<- x[[1]][["allow"]]
        cat("Fields in calibstat objects: ", join(allow,sep=" ") ,"\n")
      }
    else{
      cmv<-NULL
      allow<-NULL
    }
    invisible(list(class=cl,experiment=exp,project=pro,length=ll,classcalibstat=cmv,allow=allow))
  }

"[[<-.caliblist" <- function(x,i,value)
  {
    ##t Replace calibstat object in the caliblist.
    ##- Replace a calibstat object in the caliblist with a different one.
    ##+ x : caliblist
    ##+ i : index or name (info) of calibstat object to replace
    ##+ value : object of class calibstat (e.g. calibrestat, calibextstat).
    ##e data(mvl)
    ##e data(cal)
    ##e res<-getintcalib(mvl,cal,error=300)
    ##e res[[1]]<-res[[10]]
    
    if(!inherits(value,"calibstat"))
      stop("only calibstats allowed to assing")
    x <- NextMethod("[[<-")
    if(is.numeric(i))
      {
        names(x)[i]<- mget(value,"info")
      }
    x
  }

"[.caliblist"<-function(x,i)
  {
    ##t Extract Parts of a Caliblist.
    ##- The caliblist extends list. The calibstat objects in the list can therefore be accessed like list elements.
    ##+ x : caliblist
    ##+ i : indices of calibstat object to extract
    ##v caliblist : caliblist
    ##sa [<-.caliblist, \link[base]{[.list}
    ##e data(mvl)
    ##e data(cal)
    ##e print(cal)
    ##e res<-getintcalib(mvl,cal,error=300)
    ##e res<-res[1:10]
    ##e class(res)
    ##e plot(res)
    
    tmp<-NextMethod("[")
    res<- caliblist(class(x)[1],mget(x,"experiment"),tmp,project=mget(x,"project"))
    res
  }


caliblist <- function(class,experiment,data,...)
  {
    ##t Constructor
    ##- Returns object of class caliblist.
    ##a calibrelist, calibintlist
    ##d Returns objects of class calibrelist and calibintlist. The name of the class
    ##d have to be specified as the first argument in the constructor.
    ##d Constructor is used by function getrecalib.massvectorlist & getintcalib.massvectorlist
    ##+ class : string with class name "calibrestat","calibintstat".
    ##+ experiment : string with experiment name.
    ##+ data : a list with calibstat objects
    ##v calibrelist : If called with first (\code{class}) argument set to "calibrelist".
    ##v calibintlist : If called with first (\code{class}) argument set to "calibintlist".
    ##sa getrecalib.massvectorlist, getintcalib.massvectorlist
    ##e #Example calibrelist class:
    ##e data(mvl)
    ##e mvl<-mvl[1:10]
    ##e res <- getrecalib(mvl)
    ##e print(res)
    ##e summary(res)
    ##e image(res,what="Coef.Intercept")
    ##e image(res,what="Coef.Slope")
    ##e plot(res)
    ##e hist(res)
    ##e dres<-as.data.frame(res)
    ##e plot(dres$Coef.Intercept,dres$PQM,xlab="Coef.Intercept",ylab="PQM")
    ##e #greate subset.
    ##e res2<-subset(res,PQM>10)
    ##e length(res2)
    ##e plot(res2)
    ##e test<-applyrecalib(mvl, res2)

    if(missing(data))
      {
        res<-list()
      }
    else
      {
        if(!inherits(data[[1]],"calibstat"))
          stop(as.character(substitute(data))," must be a list of calibstat objects")
        res<-data
      }
    allow <- c("experiment","project","data")
    tmp<-list(...)
    for(x in names(tmp))
      {
        if(x %in% allow)
          {
            attr(res,x) <- tmp[[x]]
          }
      }
    attr(res,"allow") <- allow
    if(!missing(experiment))
      attr(res,"experiment") <- experiment

    if(missing(class))
      {
        class(res) <- c("caliblist","mlist","list","myobj")
      }
    else
      {
        class(res) <- c(class,"caliblist","mlist","list","myobj")
      }
    res
  }
#Copyright 2001, W. Wolski, all rights reserved.
##calibintstat - object with calibration statistics
##calibintlist - collection with calibration data.
#this is an calibration object.

calibintstat<-function(info , lm , ...)
  {
    ##t Constructor
    ##- Returns an object of class calibintstat. Is used by function \code{getintcalib}.
    ##+ info : identifier of the calibstat object. Links it with the massvector.
    ##+ lm : object of class 'lm' (linear model)
    ##+ ... : further parameters Coef.Intercept, Coef.Slope, lengthmv, nrmatch, error.mean, error.stdv, ppm.
    ##sa \link[base]{lm}
    ##v calibintstat : object of class calibintstat.
    ##e data(mv1)
    ##e data(cal)
    ##e res<-getintcalib(mv1,cal,error=2,ppm=FALSE)
    ##e class(res)
    ##e plot(res)
    ##e res<-getintcalib(mv1,cal,error=400,ppm=TRUE)
    ##e class(res)
    ##e plot(res)
    allow<-c("info","Coeff.Intercept", "Coeff.Slope","info","lengthmv","nrmatch","ppm","tcoor")

    if(missing(lm))
      {
        res<-list()
      }
    else
      {
        res<-lm
        res$Coeff.Intercept<-coef(lm)[1]
        res$Coeff.Slope<-coef(lm)[2]
      }
    tmp<-list(...)
    res$allow<- allow 
    if(!missing(info))
      res$info<-info
    for(x in names(tmp))
      {
        if(x %in% allow)
          {
            res[[x]]<-tmp[[x]]
          }
      }
    class(res)<-c("calibintstat","calibstat","lm","mylistobj")
    res
  }


as.lm.calibintstat<-function(object,...)
  {
    ##t Fitting Linear Models
    ##- attempts to coerce its argument into a linear model.
    ##+ object of class calibintstat
    ##v object of class "lm".
    al <- object$allow
    for(i in al)
      {
        object[[i]]<-NULL
      }
    object$allow<-NULL
    class(object)<-"lm"
    return(object)
  }

print.calibintstat<-function(x,...)
  {
    ##t Print calibintstat
    ##- `print' prints its argument and returns it invisibly (via `invisible(x)')
    ##+ x : calibintstat object
    ##v info : the id of the massvector
    ##v type : type of error
    ##v Intercept : Intercept of the mass dependent error
    ##v Slope : The slope of the mass dependent error
    ##v Lenght pl : Lenght of the peaklist
    ##v mean : mean of the error
    ##v stdv : stdv of the error
    ##v nrmatch : nr matches.
    ##v Xcoor : x coordinate on sample support
    ##v Ycoor : y coordinate on sample support
    ##e data(mv1)
    ##e data(cal)
    ##e res<-getintcalib(mv1,cal,error=500,ppm=TRUE)
    ##e print(res)

    inf <- mget(x,"info")
    cat("info         :",inf,"\n")
    typ <- ifelse(mget(x,"ppm"),"ppm","abs")
    cat("type         :",typ,"\n")
    int <- mget(x,"Coeff.Intercept")
    cat("Intercept    :",int,"\n")
    slop <- mget(x,"Coeff.Slope")
    cat("Slope        :",slop,"\n")
    ll <- mget(x,"lengthmv")
    cat("Length pl    :",ll,"\n")
    nrmat<- mget(x,"nrmatch")
    cat("nr match     :",nrmat,"\n")
    tcoor<-mget(x,"tcoor")
    cat("tcoor        :",tcoor,"\n")
    invisible(x)
  }

summary.calibintstat<-function(object,...)
  {
    ##t Calibintstat Summaries
    ##- Generates a summary: info, Intercept, Slope etc. for object of class calibintastat
    ##+ object : calibinstat.
    ##+ ... : further parameters.
    ##sa summary
    ##e data(mv1)
    ##e data(cal)
    ##e res<-getintcalib(mv1,cal,error=500,ppm=TRUE)
    ##e summary(res)

    res<-c(info(object),as.vector(object))
    names(res)[1]<-"info"
    class(res) <- "table"
    res
  }

as.vector.calibintstat<-function(x,mode="any")
  {
    ##t Coerces to vectors
    ##- 'as.vector', a generic, attempts to coerce its argument into a
    ##-  vector of mode 'mode' (the default is to coerce to whichever mode
    ##- is most convenient).  The attributes of 'x' are removed.
    ##+ x : object of class calibintstat
    ##+ mode : any
    ##e data(mv1)
    ##e data(cal)
    ##e res<-getintcalib(mv1,cal,error=500,ppm=TRUE)
    ##e as.vector(res)
    temp<-summary(as.lm(x))
    res <- c(ifelse(is.null(x$lengthmv),NA, x$lengthmv ),
             ifelse(is.null(x$Coeff.Intercept),NA , x$Coeff.Intercept ),
             ifelse(is.null(x$Coeff.Slope),NA , x$Coeff.Slope ),
             ifelse(is.null(x$nrmatch),NA,x$nrmatch),
             if(is.null(x$tcoor)){c(NA,NA)}else{x$tcoor},
             temp$r.squared,
             temp$adj.r.squared,
             if(!is.null(temp$fstatistic[1]))
             {
               pf(temp$fstatistic[1],temp$fstatistic[2],temp$fstatistic[3],lower.tail=FALSE)
             }
             else
             {
               NA
             }
             )
    names(res) <- c("lengthmv","Coef.Intercept","Coef.Slope","nrmatch","Xcoor","Ycoor","R.Squared","Adjusted.R.squared","p.value")
    res
  }

applycalib.calibintstat<-function(object,mv,...)
  {
    ##t Internal Calibration
    ##- Corrects the massvector for the error model stored in calibintstat object.
    ##d \bold{Internal calibration} aligns masses of
    ##d peaks to known masses and determines by linear regression a affine
    ##d function that describing the relative error. The internal
    ##d correction fails when no calibration peaks can be found.
    ##+ object :calibintstat
    ##+ mv : massvector
    ##+ ... : further parameters
    ##v massvector : calibrated massvector. 
    ##sa applyintcalib.massvector, getintcalib.massvector, correctinternal.massvector
    ##r Wolski
    ##e data(mv1)
    ##e data(cal)
    ##e res<-getintcalib(mv1,cal,error=300)
    ##e mv2<- applycalib(res,mv1)
    ##e plot(mv1[,1],mv2[,1]-mv1[,1])
    if(!inherits(mv,"massvector"))
      {
        print(class(mv))
        stop(as.character(substitute(mv)),"should be of class massvector!\n")
      }
    if(inherits(object,"lm"))
      {
        errp <- predict(object,data.frame(masstheo=mass(mv) ) )
        if(mget(object,"ppm"))
          mv[,1] <- mv[,1]/(1 - errp/1e6)
        else
          mv[,1] <- mv[,1] + errp
      }
    mv
  }


plot.calibintstat<-function(x,...)
  {
    y<-x
    rm(x)
    ##t Calibinstat Plotting
    ##- Plots a line with intercept and slope given by the error model.
    ##+ x : object of class calibintstat
    ##+ ... : graphical parameters can be given as arguments to plot.
    ##e data(mv1)
    ##e data(cal)
    ##e res<-getintcalib(mv1,cal,error=500,ppm=TRUE)
    ##e plot(res)
    ylab <- ifelse(mget(y,"ppm"),"error ppm","error abs")
    curve(y$Coeff.Intercept + y$Coeff.Slope*x,xlim=c(800,4000),xlab="m/z",ylab=ylab,main="mass dependent error functon",...)
  }


#################################################
## calibintlist
#

applycalib.calibintlist<-function(object,mvl,...)
  {
    ##t Internal Calibration
    ##- Corrects the massvectors in the massvectorlist \code{mvl} using the error model
    ##- of the \code{calibintstat} objects stored in the \code{calibintlist}.
    ##d \bold{Internal calibration} aligns masses of
    ##d peaks to known masses and determines by linear regression a affine
    ##d function that describing the relative error. The internal
    ##d correction fails when no calibration peaks can be found.
    ##+ object: calibintlist
    ##+ mvl : massvectorlist
    ##+ ... : further params
    ##v massvectorlist : calibrated massvectorlist. 
    ##sa applyintcalib.massvectorlist, getintcalib.massvectorlist, correctinternal.massvectorlist, calibintstat, caliblist
    ##r Wolski
    ##e data(mvl)
    ##e data(cal)
    ##e res<-getintcalib(mvl,cal,error=300)
    ##e mvl2<-applycalib(res,mvl)
    

    if(!inherits(mvl,"massvectorlist"))
      stop(as.character(substitute(mvl)),"have to be a object of class massvectorlist!!!\n")
                                       #calconstants -  output of function recalibrate
    for(x in 1:length(object))
      {
        tmp <- mvl[[names(object)[x]]]
        mvl[[names(object)[x]]] <- applyintcalib(tmp,object[[x]])
        if(x%%10==0)
          cat(formatC(x,width=3)," ",sep="")
        if(x%%100==0)
          cat("\n")
      }
    cat("\n")
    mvl
  }

plot.calibintlist<-function(x,...)
  {
    ##t Calibintlist Plotting
    ##- A matrix of scatterplots is produced.
    ##+ x : object of class calibintlist
    ##+ ... : graphical parameters can be given as arguments to plot.
    ##e data(mvl)
    ##e data(cal)
    ##e mvl <- mvl[1:100]
    ##e ires <- getintcalib(mvl,cal,error=250)
    ##e plot(ires)

    dat<-as.matrix(x)
    colnames(dat)<-names(as.vector(x[[1]]))
    plot(data.frame(dat),pch="*",...)
    invisible(dat)
  }


hist.calibintlist<-function(x,...)
  {
    ##t Histogram Plot
    ##- Computes histograms of the lengthmv, Coef.Intercept, Coef.Slope ...
    ##+ x : calibrelist
    ##+ ... : further graphical parameters.
    ##e data(mvl)
    ##e data(cal)
    ##e mvl <- mvl[1:100]
    ##e data(cal)
    ##e ires <- getintcalib(mvl,cal,error=250)
    ##e hist(ires)

    dat<-as.matrix(x)
    colnames(dat)<-names(as.vector(x[[1]]))
    par(mfrow=c(3,2))
    hist(dat[,1],main=colnames(dat)[1])
    hist(dat[,2],main=colnames(dat)[2])
    hist(dat[,3],main=colnames(dat)[3])
    hist(dat[,4],main=colnames(dat)[4])
    hist(dat[,5],main=colnames(dat)[5])
    hist(dat[,6],main=colnames(dat)[6])
    par(mfrow=c(1,1))
    invisible(dat)
  }




#Copyright 2004, W. Wolski, all rights reserved.
                                        #generic methods provided by the package
                                        #function used at least by more than one object.
                                        #plot
                                        #summary
                                        #hist
                                        #print

#comparing two massvectors.
compare<-function(object,...)
  {
    UseMethod("compare")
  }

as.lm <- function(object,...)
  {
    UseMethod("as.lm")
  }

#for setting and getting masses
mass <- function(object,...)
  {
    UseMethod("mass")
  }

#each object has a unique human readable field
info<-function(object,...)
  {
    UseMethod("info")
  }
#is the unique identifier either readable or not readable
#id<-function(object,...)
#  {
#    UseMethod("id")
#  }
#is either the id, or the info
#key<-function(object,...)
#  {
#    UseMethod("key")
#  }


##as.data.frame

##(replaces correct affine)
#internalCalib<-function(object,...)
#  {
#    UseMethod("internalCalib")
#  }


##searches the mascot server with the list.#
#maskotSearch<-function(object,...)
#  {
#    UseMethod("maskotSearch")
#  }

###################################################
##
## the Wool and Smilanski filtering method
##

wsFilter<-function(object,...)
  {
    UseMethod("wsFilter")
  }


wsiFilter<-function(object,...)
  UseMethod("wsiFilter")

wsdist<-function(object,...)
    UseMethod("wsdist")


#p1 -have to be a double
#p2 -can be a vector
distance <- function(p1,p2)
  {
    ##t Intra massvector mass distance
    ##- Distance between two peaks o a massvector defined as deviation from the peptide rule.
    ##+ p1 : mass
    ##+ p2 : mass
    ##v distance : distance of the two masses, as deviation from the peptide rule.
    ##r Wool A, Smilansky Z 2002. Precalibration of matrix-assisted laser desorption/ionization-time of flight spectra for peptide mass fingerprinting. {\em Proteomics.} 2(10):1365-73.
    ##r Wolski http://www.molgen.mpg.de/~wolski/mscalib
    d12<-abs(p1-p2)
    m<-1.000495
    return(mmod(d12,m))
  }

mmod <- function(x,m)
  {
    mn <- mod(x,m)
    mnt<-mn
    mns<-mnt[mnt > m/2]
    mn[mnt > m/2] <- (m - mns)
    return(mn*2)
  }

mod <-function(x,m)
  {
    t1<-floor(x/m)
    return(x-t1*m)
  }

#
#calibEternal functions
#

getextcalib <- function(object,...)
  {
    UseMethod("getextcalib")
  }


##calibrates a peaklist by a spline
applyextcalib<-function(object,...)
  {
    UseMethod("applyextcalib")
  }

calibexternal<-function(object,...)
  {
    UseMethod("calibexternal")
  }
###################################################
#
#Recalibration functions
#

recalibrate<-function(object,...)
  {
    UseMethod("recalibrate")
  }
##returns recalibobject.
getrecalib<-function(object,...)
  {
    UseMethod("getrecalib")
  }

applyrecalib<-function(object,...)
  {
    UseMethod("applyrecalib")
  }

applycalib<-function(object,...)
  {
    UseMethod("applycalib")
  }


FullWidthatHalfMaximum<-function(sumcol)
{
   minmax<-range(sumcol[,2])
                                        #   compute the half maximum
   abl1<-diff(sumcol[,2])
   ta1a<-c(abl1,1)
   ta1b<-c(1,abl1)
   minima <- which(ta1a>0 & ta1b<0)
   mmax <- which(sumcol[,2]==max(sumcol[,2]))
   if(mmax==1 | mmax==length(sumcol[,2]))
     {
          return(list(PQM = 0,hm=c(0,0,0),width=c(0,0)))
     }
   else
     {
       ll <- which(minima<mmax)
       rr <- which(minima>mmax)
       if(length(ll)>0)
         {
           lmin<- max(minima[ll])
         }
       else
         {
           lmin<-1
         }
       
       if(length(rr)>0)
         {
           rmin<- min(minima[rr])
         }
       else
         {
           rmin<-length(ta1a)
         }
       
       int<-sumcol[,2]
                                        #peak hight.
       pmax <- int[mmax] # peak maximum
       pmin <- max(c(int[lmin],int[rmin])) #peak minimum
       ph <- pmax - pmin
       if(FALSE)
         {
           plot(sumcol[,2],type="l")
           
           abline(h=pmax)
           abline(h=pmin)
           abline(v=mmax,col=2)
           abline(v=rmin,col=3)
                                        #   maxima <- which(tt<0 & tt2>0)
           abline(v=lmin,col=4)
         }
                                        #peak width
                                        #halfmax
       hm<-ph/2
       sumcc<-sumcol[lmin:rmin,]
       
       
                                        #hm <- (minmax[1] + minmax[2])/2
       dumm<-0.05
       while(TRUE){
         hi<-(pmin + hm + hm*dumm)
         lo<-(pmin + hm - hm*dumm)
         
         width <-sumcc[,1][sumcc[,2]<hi &sumcc[,2]>lo]
         if(length(width)>2) break;
         dumm<-dumm+0.05
       }
       if(FALSE)
         {
           plot(sumcc,type="l")
           abline(h=hi)
           abline(h=lo)
         }

       ind <- which(diff(width)==max(diff(width)))
       width <- c(mean(width[1:ind]),mean(width[(ind+1):length(width)]))
       
       hm<-c(pmin + hm,pmax,pmin)
       names(hm)<-c("HalfMaximum","max","min")
       width<-c(width[2]-width[1],width)
     }
   return(list(PQM=-log(as.numeric(width[1]/ph)),hm=hm,width=width))
}


####################################################
##Affine calibration.
#


##expects an lm object as second param
correctinternal <- function(object,...)
  UseMethod("correctinternal")

getintcalib <- function(object,...)
  UseMethod("getintcalib")

applyintcalib<-function(object,...)
  UseMethod("applyintcalib")

####################################################
##global calibration
#
getglobalcalib<-function(object,...)
  UseMethod("getglobalcalib")

globalcalib<-function(object,...)
    UseMethod("globalcalib")

globalcalib.default<-function(object,...)
    warning("need massvectorlist\n")


#
#get abundant masses
#
gamasses<-function(object,...)
  UseMethod("gamasses")

#
#getting and setting calibobjects
#
getcalib<-function(object,...)
  UseMethod("getcalib")

setcalib<-function(object,...)
  UseMethod("setcalib")

getdiff<-function(object,...)
  UseMethod("getdiff")

diffFilter<-function(object,...)
  UseMethod("diffFilter")


##reads a peaklist from a file.
readBruker<-function(object,...)
    UseMethod("readBruker")

peaks<-function(object,...)
  UseMethod("peaks")


peaks.default<-function(object, max=TRUE,na.rm=FALSE,...){
  ##t Find peaks
  ##- Finds peaks - neighborpeaks smaller than central.
  ##+ object : numeric array.
  ##+ max : TRUE find maxima, FALSE find minima
  ##+ na.rm : handling of na.rm values.
  ##v index : index of the peaks   
  x<-object
  if (na.rm)
   omit<-is.na(x)
  else
   omit<-FALSE
  if (max){
   rval<-1+which(diff(sign(diff(x[!omit])))<0)
  }else{
   rval<-1+which(diff(sign(diff(x[!omit])))>0)
  }
  if (na.rm)
  {
   rval<-rval+cumsum(omit)[rval]
  }
  rval
}


#method for filtering massvectors and massvectorlists
mvFilter<-function(object,...)
  UseMethod("mvFilter")

join<-function(object,...)
    UseMethod("join")

join.default<-function(object,sep,...)
{
  if(missing(sep))
    sep<-""
  res <- object[1]
  if(length(object)>1)
    {
      for(y in 2:length(object))
        {
          res<-paste(res,object[y],sep=sep)
        }
    }
  return(res)
}

"setParms<-"<-function(object,value)
  {
    UseMethod("setParms<-")
  }


experiment<-function(object, ... )
  UseMethod("experiment")


project<-function(object, ... )
  UseMethod("project")

                                        #gets an object field
mget<-function(object,...)
    UseMethod("mget")

writeF<-function(object,...)
  UseMethod("writeF")


readF<-function(object,...)
  UseMethod("readF")
#Copyright 2004, W. Wolski, all rights reserved.
getaccC<- function(pl,cal,error=500,ppm=TRUE,uniq=FALSE)
{
  ##t Find Matching Masses
  ##- Returns the indices of masses matching to each other.
  ##+ pl : vector with masses
  ##+ cal : vector with masses
  ##+ error : error in Da or ppm.
  ##+ ppm : \code{TRUE} - then error has to be given as relative error in ppm. \code{FALSE} - error are absolute error in dalton.
  ##+ uniq : \code{TRUE} - return only closest match to mass in cal. \code{FALSE} - return all matches in the error range.
  ##e getaccC(1001:1010,1001:1010,error=300,ppm=TRUE,uniq=TRUE)
  lpl <- length(pl)
  lcal <- length(cal)
  #cat("lpl ",lpl, " lcal ",lcal," lmods ",lmods,"\n")
  tmp <- max(lpl,lcal) # there can not be more matches than tmp.
  plind <- rep(0,lpl*lcal)
  calind <- rep(0,lpl*lcal)
  ind <- 0
  if(length(pl)>0 & length(cal)>0)
  {
    if(uniq)
      {
      res <- .C("getaccU", 
                as.double(pl), 
                as.integer(lpl),  
                as.double(cal), 
                as.integer(lcal),
                as.double(error),
                as.integer(plind),
                as.integer(calind),
                as.integer(ind),
                as.integer(ppm),
                PACKAGE="mscalib"
                )
    }
    else
      {
        res <- .C("getaccD", 
                  as.double(pl), 
                  as.integer(lpl),  
                  as.double(cal), 
                  as.integer(lcal),
                  as.double(error),
                  as.integer(plind),
                  as.integer(calind),
                  as.integer(ind),
                  as.integer(ppm),
                  PACKAGE="mscalib"
                  )
      }
   }
   else{
    return(list(plind=NULL,calind=NULL))
   }     
  if(res[[8]]>0)
   {
      test <- list(plind=res[[6]][1:res[[8]]]+1,calind=res[[7]][1:res[[8]]]+1)
   }
  else
    {
      test <- list(plind=NULL,calind=NULL)
    }
  return(test)
}
  

#Copyright 2004, W. Wolski, all rights reserved.
as.data.frame.mlist<-function(x,row.names=NULL,optional = FALSE)
  {
    ##t Data Frames
    ##- Turns the mlist object into a data.frame
    ##d These functions create a data frame, tightly coupled
    ##d collections of variables which share many of the properties of
    ##d matrices and of lists, used as the fundamental data structure by
    ##d most of R's modeling software.
    ##+ x : object of class massvectorlist
    ##+ ... : further arguments
    ##sa as.matrix.mlist, as.matrix
    ##e data(mvl)
    ##e tmp <- as.data.frame(mvl)
    ##e names(mvl)
    ##e plot(tmp$lengthmv,tmp$mass.Min.)
    ##e data(mvl)
    ##e mvl<-mvl[1:100]
    ##e data(cal)
    ##e test<-getintcalib(mvl,cal,error=500)
    ##e tmp<-as.data.frame(test)
    ##e names(tmp)
    res <- NULL
    tmp <- as.matrix(x)
    ntmp <- names(x)
    res <- data.frame(info=as.character(ntmp),tmp)
    res
  }


summary.mlist<-function(object,...)
  {
    ##t mist Summaries
    ##- Generates a summary for the data.frame generated by as.data.frame.mlist
    ##+ object : mlist
    ##+ ... : further arguments
    ##sa summary.massvector,as.data.frame.mlist
    ##e data(mvl)
    ##e summary(mvl)
    ##e data(mvl)
    ##e mvl<-mvl[1:100] 
    ##e data(cal)
    ##e test <- getintcalib(mvl,cal,error=500)
    ##e summary(test)

    res<-list(info = info(object))
    res<-c(res,list(summary=summary(as.matrix(object))))
    return(res)
  }


as.matrix.mlist<-function(x)
{
  ##t Matrices
  ##- Turns the caliblist into a matrix
  ##+ x : caliblist
  ##v matrix : matrix
  ##sa as.data.frame.mlist
  ##e data(mvl)
  ##e mvl<-mvl[1:100]
  ##e data(cal)
  ##e test<-getintcalib(mvl,cal,error=500)
  ##e tmp<-as.matrix(test)
  ##e colnames(tmp)
  ##e dim(tmp)
  ##e data(mvl)
  ##e tmp<-as.matrix(mvl)
  ##e print(colnames(mvl))
  ##e plot(tmp[,"lengthmv"],tmp[,"mass.Min."])
  dat<-t(sapply(x,as.vector))
  dat
}


image.mlist<-function(x,what="",col=terrain.colors(100),...)
  {
    ##t Display a Color Image
    ##- Creates a grid of colored or gray-scale rectangles with colors
    ##- corresponding to the values in 'z'.  This can be used to display
    ##- three-dimensional or spatial data aka "images". This is a generic
    ##- function.
    ##+ x : object of class mlist. (e.g: caliblist or massvectorlist)
    ##+ what : what value to display on the image.
    ##+ col : a list of colors such as that generated by 'rainbow', 'heat.colors', 'topo.colors', 'terrain.colors' or similar functions.
    ##e data(mvl)
    ##e image(mvl,what="lengthmv")
    cal<-x
    rm(x) # do not like to use x.
    if(length(cal)==0){
      warning("List has length 0")
      return()
    }
    res <- as.matrix(cal)
    if(! what %in% colnames(res))
      stop("Only following fields can be shown for ", class(cal)[1]," : \n", join(colnames(res),sep=" "),"\n pass one to the what paramter.")
    nam <- mget(cal,"tcoor")
    if(!is.null(names(nam$coorX)) & !is.null(names(nam$coorY)))
      {
        X <- nam$coorX[unique(names(nam$coorX))]
        Y <- nam$coorY[unique(names(nam$coorY))]
        XX <- 1:max(X)
        names(XX)<-rep("",max(X))
        names(XX)[X]<-names(X)
        X<-XX
        YY <- 1:max(Y)
        names(YY)<-rep("",max(Y))
        names(YY)[Y]<-names(Y)
        Y<-YY
      }
    else
      {
        X <- unique(nam$coorX)
        Y <- unique(nam$coorY)
      }
    hello <- matrix(NA,max(X),max(Y))
   


    for(z in 1:length(cal))
      {
        hello[ nam$coorX[z] , nam$coorY[z] ] <- res[z,what]
      }

    if(!is.null(names(X)) & !is.null(names(Y)))
       {
         rownames(hello) <- names(X)
         colnames(hello) <- names(Y)
       }

    par(bg="gray")
    tmar<-par()$mar
                                        #define layout
    nf <- layout(matrix(c(1,2),1,2),widths=c(5,1), TRUE)
    par(mar=c(3,3,2,0.5))
                                        #2.03.2004
    image(t(hello) , main=what,xaxs="i",yaxs="i",axes=FALSE,col=col,...)
    #image(hello , main=what,xaxs="i",yaxs="i",axes=FALSE,col=col,...)
    if((length(Y)-1)>0)
      {
        axis( 1 , at=seq(0,1,1/(length(Y)-1)) , labels=names(sort(Y)))
      }
    else
      {
         axis(1,at=0.5,labels=names(Y))
      }
    if((length(X)-1)>0)
      {
        axis( 2 , at=seq(0,1,1/(length(X)-1)) , labels=names(sort(X)))
      }
    else
      {
        axis(2,at=6,labels=names(X))
      }
    tres<-na.omit(c(hello))
    if(min(tres)!=max(tres))
      {
        scale<-seq(min(tres),max(tres),(max(tres)-min(tres))/9)
      }
    else
      {
        scale<-rep(min(tres),10)
      }
    scale<-matrix(scale,nrow=1)
    par(mar=c(3,0,2,0.5))
                                        #check for scalig factor
    #pp<-mget(cal[[1]],"ppm")
    lable<-""
    #if(!is.null(pp))
    #  {
    #    lable<-ifelse(pp,"* -1e4","* -1e6")
    #  }
    image(1,1:10,scale,axes=FALSE,xlab="",ylab="",col=col,main=lable )
    scalet<-format(scale,digits=1)
    for(x in 1:length(scalet))
      {
        text(0.5,x,scalet[x])
      }
    par(mar=tmar)
    layout(matrix(1))
    invisible(t(hello))
  }


mget.mlist<-function(object,attrn,...)
  {
    ##t Field Access
    ##- Acces fields in object of class myobj
    ##+ object : object of class  mlist
    ##+ attrn : name of field (Attribute)
    ##e data(mvl)
    ##e mget(mvl)
    ##e mget(mvl,"info")
    ##e mget(mvl,"tcoor")
    if(missing(attrn))
      return(attr(object,"allow"))
    if(attrn %in% "tcoor")
      {
        res <- NULL
        nres<-NULL
        for(x in object)
          {
            res <- rbind( res , mget(x,"tcoor") )
            nres <- rbind( nres,names(mget(x,"tcoor")) )
          }
        coorX <- res[,1]
        names(coorX) <- nres[,1]
        coorY<-res[,2]
        names(coorY)<-nres[,2]
        return(list(coorX=coorX,coorY=coorY))
      }
    else
      NextMethod("mget")
  }

experiment.mlist<-function(object,experiment,...)
  {
    ##t Info Acces
    ##- Access to the experiment field of the mlist. Can be used for setting or getting the experiment field.
    ##a info.mlist
    ##+ object : object of class mlist
    ##+ experiment : New experiment name. If not missing function returns massvector with new info field content.
    ##e data(mvl)
    ##e experiment(mvl)
    ##e mvl<-experiment(mvl,"newname")
    if(missing(experiment))
      return(mget(object,"experiment"))
    else
      {
        setParms(object)<-list(experiment=experiment)
      }
    object
  }

info.mlist<-experiment.mlist

subset.mlist <- function(x,subset,...)
  {
    ##t Subset mlist
    ##- Return subsets of list elements which meet conditions.
    ##+ x : object of class mlist
    ##+ subset : logical expression.
    ##e data(mvl)
    ##e mvl<-subset(mvl,lengthmv>30)
    u <- as.data.frame(x)
    if (missing(subset))
      {
        r <- TRUE
        cat("For subsetting use comparison on : \n", join(names(u),sep=" ")  ,"\n")
        return()
      }
    else {
        e <- substitute(subset)
        r <- eval(e, u, parent.frame())
        r <- r & !is.na(r)
    }
    vars <- TRUE
    u <- u[r, vars, drop = FALSE]
    x[as.character(u$info)]
}
#Copyright 2004, W. Wolski, all rights reserved.

                                        #constructor
                                        #massvector
                                        #infofield
                                       #coordinates on target
massvector <- function(info,masses,tcoor)
{
  ##t Constructor
  ##- massvector extends matrix
  ##d The class massvector keeps the masses and their intensities in a matrix.
  ##d It stores additional attributes like the coordinates of the mass spectrometric sample on the support.
  ##d It also stores a identifier of the massvector.
  ##+ masses : a matrix or a array with masses. First column must contain masses.
  ##+ info : a unique identifier for the list.
  ##+ tcoor : a array of length two if integer sample support coordinates.
  ##e massvector("hello march",NULL)
  ##e massvector("hello march",1:100)
  ##e massvector("hello march",cbind(1:10,10:1))
  ##e tmp<-cbind(1:10,1:10)
  ##e colnames(tmp)<-c("mass","test")
  ##e rr<-massvector("hello march",tmp)
  ##e rr2<-massvector("hello bart",cbind(1:12,1:12))
  ##e # plot functions for massvector
  ##e plot(rr)
  ##e hist(rr)
  ##e summary(rr)
  ##e image(rr)
  ##e info(rr)
  ##e # setting new masses
  ##e mass(rr)
  ##e mass(rr,1:10)
  ##e peaks(rr,cbind(1:10,11:20))
  ##e # plotting with masses
  ##e data(mv1)
  ##e data(mv2)
  ##e plot(mv1,mv2)
  ##e image(mv1  ,mv2)
  ##e summary(mv1)
  ##e print(mv1)
  ##e hist(mv1)
  ##e plot(mv1)
  ##e image(mv1,mv2,error=199,ppm=FALSE)


  
  res<-matrix()
  if(!missing(masses))
    {
      if(is.null(masses))
        {
          masses <- matrix(nrow=0,ncol=2)
          colnames(masses)<-c("mass","C1")
        }
      else
        {
          #if not a matrix.
          if(!inherits(masses,"matrix"))
            {
              masses<-cbind(masses,rep(1,length(masses)))
              colnames(masses)<-c("mass",paste("C",2:length(masses[1,]),sep=""))
            }
          if(dim(masses)[2]<2)
            masses<-cbind(masses,rep(1,length(masses)))
                                        #order the matrix
          or<-order(masses[,1])
          for(i in 1:dim(masses)[2])
            {
              masses[,i]<-masses[,i][or]
            }
          if(length(masses[,1])>0)
            rownames(masses) <- 1:length(masses[,1])

          if(is.null(colnames(masses)))
            {
              colnames(masses) <- c("mass",paste("C",2:length(masses[1,]),sep=""))
            }
          else
            colnames(masses)[1] <- "mass"
        }
    }
  else
    {
      masses<-matrix(nrow=0,ncol=2)
      colnames(masses)<-c("mass","C1")
    }
  if(!missing(info))
    {
      attr(masses,"info")<-info
    }
  if(!missing(tcoor))
    {
      if(length(tcoor)!=2){stop("There must be two coordinates!\n")}
      attr(masses,"tcoor") <- tcoor
    }
  attr(masses,"allow")<-c("info","tcoor","gelcoor")
  class(masses)<-c("massvector","myobj","matrix")
  return(masses)
}



summary.massvector<-function(object,...)
  {
    ##t Massvector Summaries
    ##- Generates a summary: min, max etc.
    ##+ object : massvector
    ##sa summary
    tmp <- dim(object)[2]
    res <- list( lengthmv=length(object) )
    for(i in 1:tmp)
      {
        xx<-object[,i]
        res <- c(res,list(summary(ifelse(length(xx)>0,xx,0))))
        names(res)[i+1]<-colnames(object)[i]
      }
    res
  }

as.vector.massvector<-function(x, mode="any")
  {
    ##t Vector
    ##- Casts massvector into vector. It does it by calling the summary.massvector method first and unlisting the result.
    ##+ x : massvector
    ##e data(mv1)
    ##e as.vector(mv1)
    unlist(summary(x))
  }

mass.massvector<-function(object,mas,...)
  {
    ##t Mass Access.
    ##- Access to the mass field of a massvector.
    ##+ object : massvector
    ##+ mas : A array with masses. If missing function returns masses.
    ##e data(mv1)
    ##e mass(mv1)
    ##e mass(mv1,1:10)
    
    if(missing(mas))
      return(object[,1])
    else
      {
        if(inherits(mas,"matrix"))
           stop(substitute(object) ," : has to be an array, to set a matrix use peaks instead!")
        masses<-cbind(sort(mas),rep(1,length(mas)))
        return(masses)
        rownames(masses)<-1:length(mas)
        colnames(masses)<-c("mass","C1")
        tmp <- attributes(object)
        at <- names(attributes(masses))
        for(l in 1:length(tmp))
          {
            if(! names(tmp)[l] %in% at)
              {
                attr(masses,names(tmp)[l])<-tmp[[l]]
              }
          }
        return(masses)
      }
  }


hist.massvector<-function(x,accur = 0.1,abund = 0, main=info(x) ,xlab="m/z",xlim=c(min(mass(x)),max(mass(x))),add=FALSE,col=1,...)
  {
    ##t Histograms
    ##- Histograms
    ##+ accur : sets the bin width of the histogramm.
    ##+ abund : draws a horizontal line at the frequency given by abund.
    ##+ xlab : sets the xlabels.
    ##+ xlim : sets the min and max value to be displayed.
    ##+ add : T-adds the histogram to an existing image.
    ##+ col : the color of the histogram.
    ##+ ... : further plotting arguments.
    ##sa hist
    ##e data(mv1)
    ##e hist(mv1)
    mhist<-list(NULL)
                                        #assign indices of bins with high peak abundance
    mhist[[1]] <- hist(mass(x),breaks=seq(min(mass(x))-accur/2,max(mass(x))+1.5*accur,accur),plot=TRUE,main=main,xlab=xlab,xlim=xlim,add=add,col=col,border=col,...)
    mhist[[2]] <- hist(mass(x),breaks=seq(min(mass(x))-accur,max(mass(x))+accur,accur),add=TRUE,col=col,border=col)
    abline(h=abund,lty=2)
 
  }




#setting and getting masses

peaks.massvector<-function(object,masses,...)
  {
    ##t Data Access
    ##- Access to the mass field of the massvector.
    ##+ object : massvector
    ##+ masses : matrix with masses in first column. If missing the matrix of the massvector is returned.
    ##sa mass.massvector
    ##e mv1<-massvector()
    ##e mv1<- peaks(mv1,cbind(1:10,1:10))
    ##e peaks(mv1)
    if(missing(masses))
      {
        for(y in attr(object,"allow"))
          {
            attr(object,y)<-NULL
          }
        attr(object,"allow")<-NULL
        class(object)<-"matrix"
        return(object)
      }
    else
      {

        if(!inherits(masses,"matrix"))
          {
            stop("The second arg (masses) must be a matrix\n")
          }
        or<-order(masses[,1])
        for(i in 1:length(dim(masses)[2]))
          {
            masses[,i]<-masses[,i][or]
          }
        if(is.null(rownames(masses)))
          rownames(masses)<-1:length(masses[,1])
        if(is.null(colnames(masses)))
          colnames(masses)<-c("mass",paste("C",2:length(masses[1,]),sep=""))
        else
          colnames(masses)[1]<-c("mass")

        tmp <- attributes(object)
        at <- names(attributes(masses))

        for(l in 1:length(tmp))
          {
            if(! names(tmp)[l] %in% at)
              attr(masses,names(tmp)[l])<-tmp[[l]]
          }
        return(masses)
      }
  }

image.massvector<-function(x,mv2,error=NULL,ppm=FALSE,col=topo.colors(100),...)
  {
    ##t Display a Color Image
    ##- Creates a grid of colored or gray-scale rectangles with colors
    ##- corresponding to the mass differences within the peaklist
    ##- or within two peaklists.
    ##+ x : massvector
    ##+ mv2 : massvector
    ##+ error : up to which mass difference display the differences.
    ##+ col : a list of colors such as that generated by `rainbow',`heat.colors', `topo.colors', `terrain.colors' or similar functions.
    ##... : graphical parameters for `plot' may also be passed as arguments to this function.
    ##sa plot.massvector, hist.massvector
    ##e data(mv2)
    ##e data(mv1)
    ##e image(mv1,mv2)
    ##e image(mv1,mv2,error=500)
    
    if(!missing(mv2))
      {
        if(!inherits(mv2,"massvector"))
          stop("Second arg should be a massvector too!\n")
      }
    else
      {
        mv2<-x
      }
    res<-NULL
    for(u in 1:length(x[,1]))
      {
        if(ppm)
          tmp<-abs(mv2[,1]-x[u,1])/mv2[,1]*10e6
        else
          tmp<-abs(mv2[,1]-x[u,1])
        tmp[tmp>error]<-NA
        res<-rbind(res,tmp)
      }
    par(bg="gray")
    nf <- layout(matrix(c(1,2),1,2),widths=c(4,1), TRUE)
    tmar<-par()$mar
    par(mar=c(5,5,1,1))
    image(1:length(x[,1]),1:length(mv2[,1]),res,col=col,xlab=info(x),ylab=info(mv2))
    tres<-na.omit(c(res))
    if(min(tres)!=max(tres))
      {
        scale<-seq(min(tres),max(tres),(max(tres)-min(tres))/9)
      }
    else
      {
        scale<-rep(min(tres),10)
      }
    scale<-matrix(scale,nrow=1)
    par(mar=c(5,1,1,1))
    image(1,1:10,scale,axes=FALSE,xlab="",ylab="",col=col)
    scalet<-format(scale,digits=1)
    for(u in 1:length(scalet))
      {
        text(25,u,scalet[u])
      }
    layout(matrix(1))
    par(mar=tmar)
    invisible(res)
  }

min.massvector<-function(mv,...)
  {
    ##t Maxima and Minima
    ##- Returns the maxima and minima of the masses and intensities
    ##sa max.massvector
    ##+ mv : massvector
    ##v named array with minima of the columns of mv.
    return(apply(mv,2,min))
  }

max.massvector<-function(mv,...)
  {
    ##t Maxima and Minima
    ##- Returns the maxima and minima of the masses and intensities
    ##sa min.massvector
    ##+ mv : massvector
    ##v named array with maxima of the columns of mv.
    return(apply(mv,2,max))
  }

length.massvector<-function(x,...)
  {
    ##t Length of Massvector
    ##- get the Length of a Massvector
    ##+ x : massvector
    ##e data(mv1)
    ##e length(mv1)
    
    return(length(x[,1]))
  }


plot.massvector<-function(x,...)
  {
    ##t Massvector Plotting
    ##- Function for plotting massvectors. If one massvector are given it shows a stick masspectrum.
    ##- If two massvectors are given their masses are plotted against each other.
    ##+ x : massvector
    ##+ ... : a second massvector and graphical parameters can be given as arguments to `plot'.
    ##e data(mv1)
    ##e data(mv2)
    ##e plot(mv1)
    ##e plot(mv1,mv2)
    mv<-x
    rm(x)
    if(length(mv)==0)
      return(paste(as.character(substitute(x))," has length 0!",sep=""))
    pars<-list(...)
    tmp<-FALSE
    if(length(pars)>0)
      {
        tmp<- inherits(pars[[1]],"massvector")
      }
    if(tmp)
      {
        plot(1,1
             ,xlim = if(is.null(pars[["xlim"]])){c(min(mv[,1]),max(mv[,1]))}else{pars[["xlim"]]}
             ,ylim = if(is.null(pars[["ylim"]])){c(min(pars[[1]][,1]),max(pars[[1]][,1]))}else{pars[["ylim"]]}
             ,ylab = info(pars[[1]]),xlab=info(mv))
        abline(v=mv[,1])
        abline(h=pars[[1]][,1])
        if(!is.null(pars[["error"]]))
          {
            abline(coef=c(0,1),col=2)
            abline(coef=c(pars[["error"]],1),col=3)
            abline(coef=c(-pars[["error"]],1),col=3)
          }
        else
          {
            abline(coef=c(0,1),col=2)
            abline(coef=c(1,1),col=3)
            abline(coef=c(-1,1),col=3)
          }
      }
    else
      {
        parsn<-names(pars)
        if("main" %in% parsn)
          {
            m <- pars[["main"]]
            pars <- subset(pars,select=-main)
          }
        else
          m <- info(mv)
        x<-mv[,1]
        y<-mv[,2]
        ylim<-c(0,max(y))
        xlim<-c(min(x),max(x))
                                        #        xlim<-c(min(x),max(x))
        sw<-FALSE
        if(!is.null(pars[["add"]]))
          sw<-pars[["add"]]
        
        if("add" %in% parsn)
          {
            pars <- as.list(subset(as.data.frame(pars),TRUE,select=-add))
          }
        parj<-join(paste(names(pars),pars,sep="="),sep=",")
        if(is.na(parj))
          {
            parj<-""
          }
        if(sw)
          {
            test <- parse(text=paste("points(x,y,type=\"h\",",parj,")",sep="\n"))
          }
        else
          {
            test <- parse(text=paste("plot.default(x,y,type=\"h\",main=\"",m,"\",axes=TRUE,xlab=\"m/z\",ylab=\"",colnames(mv)[2],"\",ylim=\c(",join(ylim,sep=","),"),xlim=\c(",join(xlim,sep=","),"),",parj ,")",sep="") )
          }
        dataf <- eval(test)
        abline(h=0)
      }
  } 


                                        #filtering
mvFilter.massvector<-function(object,fmass,match=FALSE,error=250,ppm=TRUE,uniq=FALSE,...)
  {
    ##t Filtering Massvector
    ##- Filters massvector for masses given in a second massvector.
    ##+ object : massvector
    ##+ abundm : massvector
    ##+ error : mesurment error.
    ##+ ppm : given either in ppm (\code{TRUE}) or as absolut error (F).
    ##+ match: logical; \code{TRUE} - than returns masses matching to the masses in massvector abundant, \code{FALSE} - returns masses not matching.
    ##+ uniq : logical; \code{FALSE} - returns all masses in the range given by error. \code{TRUE} - returns only the closest mass.
    ##v massvector : with matching or not matchin masses.
    ##e data(mv1)
    ##e data(mv2)
    ##e mvFilter(mv1,mv2,error=250,match=FALSE)
    ##e mvFilter(mv1,mv2,error=250,match=TRUE)

    mmatch<-getaccC(mass(object),mass(fmass),error=error,ppm=ppm,uniq=uniq)
    if(length(mmatch$plind)>0)
      {
        if(!match)
          {
            return(object[-mmatch$plind,])
          }
        else
          {
           return(object[mmatch$plind,])
          }
      }
    else
      {
        if(!match)
          return(object)
        else
          return(object[NULL,])
      }
  }

print.massvector <- function(x,quote=FALSE,...)
  {
    ##t Print massvector
    ##- `print' prints its argument and returns it invisibly (via `invisible(x)')
    ##+ x : massvector
    ##sa \link[base]{print}
    ##e data(mv1)
    ##e print(mv1)
        
    print(paste("info   : (",attr(x,"info"),")",sep=""),quote=quote,...)
    inf<-info(x)
    mcor<-attr(x,"tcoor")
    ncor<-names(mcor)
    if(!is.null(ncor))
      print(paste("coor   :  ",ncor[1],"=",mcor[1]," ; ",ncor[2],"=",mcor[2] ,sep=""),quote=quote,...)
    else
      print(paste("coor   :  ",mcor[1]," ; ",mcor[2] ,sep=""),quote=quote,...)
    attr(x,"allow") <- NULL
    attr(x,"info") <- NULL
    attr(x,"tcoor") <- NULL
    class(x) <- NULL
    print(x,...)
    invisible(x)
  }


compare.massvector <- function(object,mv2,plot=TRUE,error=1000,ppm=TRUE,uniq=FALSE,...)
  {
    ##t Compares massvectors
    ##- Compares the masses in the massvectors. Returns basic statistics about the matching peaks.
    ##- Plots the relative or absolute error of matchin peaks.
    ##+ object : massvector
    ##+ mv2 : massvector
    ##+ plot : True - plot the relatvetor absolute error. \code{FALSE} - no plotting.
    ##+ error : size of the measurment error (default 150 ppm)
    ##+ ppm : \code{TRUE} - relative error in parts per million, \code{FALSE} - absolute error.
    ##+ ... : further parameters.
    ##+ uniq : logical : default \code{FALSE}. Select all peaks matching in a massrange.
    ##v FMSTAT : Fowlkes & Mallows statistik (nr matching)/sqrt(length(object)*length(mv2))
    ##v min : smallest error
    ##v ... : 1st qu. , mean, median, 3rd qu., max and stdv of error.
    ##sa resid
    ##e data(mv1)
    ##e data(mv2)
    ##e compare(mv1,mv2,error=5000,ppm=TRUE,uniq=TRUE)
    ##e compare(mv2,mv1,error=1,ppm=FALSE,uniq=TRUE)

    main<-list(...)$main
    if(!inherits(mv2,"massvector"))
      stop("Second arg are not a massvector!\n")
    match<-getaccC(mass(object) , mass(mv2) , error=error , ppm=ppm , uniq=uniq)
    if(ppm)
      {
        if(plot)
          plot(object[match$plind,1],(object[match$plind,1]-mv2[match$calind,1])*1e6/mv2[match$calind,1],main=ifelse(is.null(main),paste(info(object),info(mv2)),main),xlab="mass",ylab="error[ppm]")
        res <- (object[match$plind,1]-mv2[match$calind,1])*1e6/mv2[match$calind,1]
      }
    else
      {
        if(plot)
          plot(object[match$plind,1],(object[match$plind,1]-mv2[match$calind,1]),main=ifelse(is.null(main),paste(info(object),info(mv2)),main),xlab="mass",ylab="Da")
        res <- (object[match$plind,1]-mv2[match$calind,1])
        
      }
    if(length(res)>0)
      {
        res<-c(length(match$plind)/sqrt(length(object)*length(mv2)),min(res),quantile(res,0.25),mean(res),median(res),quantile(res,0.75),max(res),sqrt(var(res)))
      }
    else
      {
         res<-c(length(match$plind)/sqrt(length(object)*length(mv2)),NA,NA,NA,NA,NA,NA,NA)
      }
    names(res)<-c("FMSTAT","min","1st qu.","mean","median","3rd qu.","max","stdv")
    invisible(res)
  }

mget.massvector<-function(object,attrn,...)
  {
    ##t Field Access
    ##- Access to the fields in the massvector
    ##+ object : massvector
    ##+ attrn : The value of which field to return. If missing the fields of the objects are returned.
    ##+ ... : further parameters
    ##v xxx : depends which field in the massvector are accessed.
    ##e data(mv1)
    ##e mget(mv1,"info")
    ##e mget(mv1,"peaks")
    
    if(missing(attrn))
      {
        return(attr(object,"allow"))
      }
    else
      {
        if(attrn %in% "peaks")
          {
            al<-attr( object , "allow" )
            for(x in al)
              attr(object,x) <- NULL
            attr(object,"allow")<- NULL
            class(object)<-"matrix"
            return(object)
          }
        else
          NextMethod("mget")
      }
  }


as.matrix.massvector <- function(x)
  {
    ##t Matrices
    ##- Turns the massvector into a matrix
    ##+ x : massvector
    ##v matrix : matrix with masses and intensities.
    ##sa peaks.massvector
    ##e data(mv1)
    ##e res<-as.matrix(mv1)
    ##e class(res)
    return(mget(x,"peaks"))
  }

"[.massvector"<-function(peak,i,j)
  {
    ##t Extract Parts of an Massvector
    ##- The massvector extends matrix. The masses and intensities can be accessed like in case of a matrix.
    ##- Do not use mv[1:10] (not handled properly)
    ##+ peak : object from which to extract elements.
    ##+ i : row access
    ##+ j : column access
    ##v xxx : If rows are selected then a massvector is returned. If columns are accessed than arrays are returned.
    ##sa peaks.massvector, mass.massvector, mget.massvector
    ##e data(mv1)
    ##e ls()
    ##e mv1[1:10,] # returns a massvector of length 10
    ##e mv1[1:10,1] # the first ten masses are returned.
    ##e mv1[,1] # the masses are returned.
    ##e mv1[,2] # the peak area are returned.
 
    
    res<-NextMethod("[")
    if(missing(j) & !missing(i))
                                        #     if(inherits(res,"matrix"))
      {
        if(!inherits(res,"matrix"))
          {
            tt<-matrix(res,ncol=dim(peak)[2])
            colnames(tt)<-colnames(peak)
            rownames(tt)<-1:length(tt[,1])
            res<-tt
          }
        at<-names(attributes(res))
        atr<-attributes(peak)
        for(x in 1:length(atr))
          {
            cur <- names(atr)[x]
            if(!cur %in% at)
              attr(res,cur) <- atr[[x]]
          }
      }
    res
  }

c.massvector<-function(x,...)
  {
    ##t Combine Massvectors into one Massvector.
    ##- Combines Massvectors into one Massvector.
    ##+ x : massvector
    ##+ ... : massvectors to be concatenated.
    ##sa rbind
    ##v massvector : massvector
    ##e data(mv1)
    ##e data(mv2)
    ##e par(mfrow=c(2,1))
    ##e plot(mv1)
    ##e plot(mv2,add=TRUE)
    ##e plot(c(mv1,mv2))
    
    tmp<-list(...)
    for(y in tmp)
      {
        if(inherits(y,"massvector"))
          x <- massvector(paste(info(x),info(y),sep="_"),rbind(x,y))
      }
  }
                                        #------------------------------------------------





getPPGmasses <- function(start=10,end=100)
{
  ##t PPG masses
  ##- Computes poly-(propylene glycol) masses.
  ##+ start : length of shortest polymer.
  ##+ end : length of longest polymer.
  ##v massvector : massvector of ppg masses
  ##e plot(getPPGmasses(start=12,end=100))
  mC  <-  12
  mH <- 1.007825
  mO <- 15.994915
  mNa <- 22.98977
  me <- 0.00054858
  massPPG <- mC*3 + mH*6 + mO
  massOH <- mO+mH

  n <- start:end
  n <- massPPG * n + massOH + mH + mNa - me
  names(n) <- start:end
  tt <-cbind(n,start:end)
  colnames(tt)<-c("mass","n")
  massvector("theoretical ppg masses",tt)
}

##################################################
#
# Wool Smilanski filtering
#

wsiFilter.massvector<-function(object , mdist=0.25 , fraction=0.2 , ... )
  {
    ##t Smilanski Filtering
    ##- The function returns the inidces of masses identified as chemical noise.
    ##d Chemical noise can be removed from the peptide mass lists
    ##d due to the strong clustering of mono-isotopic peptide
    ##d peaks. Following the distance measure and filtering
    ##d method proposed by Wool Smilanski we developed an algorithm to
    ##d classify masses as peptide and non-peptide. The algorithm is based
    ##d on a modified distance measure and hierarchical clustering of all
    ##d intra massvector distances.
    ##+ object : massvector.
    ##+ mdist : minimal distance to branch to be prune. The unit of this distance are Daltons.
    ##+ fraction : maximal fraction (nr masses in branch)/(length of massvector) of branche  to be prune.
    ##+ ... : further arguments.
    ##v indices :  Indices of masses which are identified as being nonpeptide.
    ##r Wool A, Smilansky Z 2002. Precalibration of matrix-assisted laser desorption/ionization-time of flight spectra for peptide mass fingerprinting. \emph{Proteomics.} 2(10):1365-73.
    ##r Wolski
    ##sa wsFilter.massvector,wsdist.massvector, wsFilter.massvectorlist
    ##e data(mv1)
    ##e data(mv2)
    ##e length(mv1)
    ##e length(wsFilter(mv1))
    ##e length(mv2)
    ##e length(wsFilter(mv2))
    
    if(length(object)<2)
      return(NULL)
                                        #die function gibt die indices der contaminanten zurck.
                                        #der aufrufenden function bleibt berlassen was sie mit ihnen macht.
                                        #pl - peaklist
                                        #mdist - distance of the contaminant cluster.
                                        #fraction - the contaminant cluster shouldnt be greater than fraction.
    if(is.null(rownames(object)))
      rownames(object)<-1:length(object)
    hst <- hclust(wsdist(object),method="single")
    tmp <- which(hst$height > mdist)
    if(length(tmp)>0)
      {
         tmp<-cutree(hst,2)
                                        #count how big are the branches.
         c2<-sum(tmp==2)
         c1<-sum(tmp==1)
         if(min(c(c2,c1))/max(c(c2,c1))< fraction)
           {
             if(c2>c1)
               {
                 res<- as.numeric(rownames(object)[tmp==1])
                 res<-c( res , wsiFilter(object[tmp==2,] , mdist=mdist , fraction=fraction))
                 return(as.numeric(res))
               }else{
                 res<-as.numeric(rownames(object)[tmp==2])
                 res<-c(res , wsiFilter(object[tmp==1,] , mdist=mdist , fraction=fraction))
                 return(res)
               }
           }
      }
                                        # if no contaminations where found.
                                        #wenn keine contaminationen gefunden wruden.
    return(NULL)
  }

wsFilter.massvector <- function(object,mdist=0.25,fraction=0.2, peptides=TRUE,...)
  {
    ##t Smilanski Filtering
    ##- Removes chemical noise from the massvectorlist.
    ##d Chemical noise can be removed from the peptide mass lists
    ##d due to the strong clustering of mono-isotopic peptide
    ##d peaks. Following the distance measure and filtering
    ##d method proposed by Wool Smilanski we developed an algorithm to
    ##d classify masses as peptide and non-peptide. The algorithm is based
    ##d on a modified distance measure and hierarchical clustering of all
    ##d intra massvector distances.
    ##+ object : massvector.
    ##+ mdist : minimal distance of branch to be cut.
    ##+ fraction : maximal size of branche (nr masses in branch)/(length of massvector) to be cut.
    ##+ peptides : \code{TRUE} - returns peptides, \code{FALSE} - returns chemical noise.
    ##+ ... : further parameters.
    ##v massvector :  returns either a massvector of nonpeptide masses or massvector of peptide masses.
    ##r Wool A, Smilansky Z 2002. Precalibration of matrix-assisted laser desorption/ionization-time of flight spectra for peptide mass fingerprinting. Proteomics. 2(10):1365-73.
    ##sa wsiFilter.massvector,wsFilter.massvectorlist
    ##e data(mv1)
    ##e tmp <- wsFilter(mv1,peptide=FALSE)
    ##e plot(tmp)
    ##e tmp <- wsFilter(mv1,peptide=TRUE)
    ##e plot(tmp)
    
    
    tmp<-wsiFilter(object,mdist=mdist,fraction=fraction)
    if(peptides)
      {
        if(length(tmp)==0)
          return(object)
        else
          return(object[-tmp,])
      }
    else
      {
        return(object[tmp,])
      }
  }




wsdist.massvector <- function(object,...)
  {
    ##t Wool Smilanski Distance Matrix
    ##- This function computes and returns the distance matrix
    ##- using the intra massvecotor distance  measure.
    ##+ object : massvector
    ##v dist : an object of class distance.
    ##sa wsFilter.massvector, wsiFilter.massvector,
    ##r Wool A, Smilansky Z 2002. Precalibration of matrix-assisted laser desorption/ionization-time of flight spectra for peptide mass fingerprinting. \emph{Proteomics.} 2(10):1365-73.
    ##r Wolski \url{http://www.molgen.mpg.de/~wolski/mscalib}
    ##e data(mv1)
    ##e plot(hclust(wsdist(mv1),method="single"))
    ret<-NULL
    pl<-object[,1]
    for(x in pl)
    {
      ret<-rbind(ret,distance(x,pl))
    }
    rownames(ret) <- formatC(pl,digits=2,format="f")
    colnames(ret) <- format(pl,digits=2,format="f")
    ret <- as.dist(ret)
    return(ret)
  }


###########################################
# External calibration


  
getextcalib.massvector <- function(object,calib,error=300,...)
{
  ##t External Error Model
  ##- Returns the error model obtained from the calibration sample.
  ##d In case of \bold{external calibration} some sample spots are only dedicated
  ##d to calibration. Calibration samples which produces equidistant
  ##d peaks, which exact masses are known, can be used to precisely
  ##d estimate the mass dependent error function.
  ##+ object : massvector
  ##+ calib : massvector with calibration masses
  ##+ error : relative measurment error in ppm.
  ##v calibspline : can be used to calibrate peaklists
  ##sa calibspline, applyextcalib.massvector
  ##r Gobom J, Mueller M, Egelhofer V, Theiss D, Lehrach H, Nordhoff E, 2002. A calibration method that simplifies and improves accurate determination of peptide molecular masses by MALDI-TOF MS. \emph{Anal Chem.} 74(15):3915-23.
  ##r Wolski \url{http://www.molgen.mpg.de/~wolski/mscalib}
  ##e data(mv1)
  ##e data(ppg)
  ##e res<-getextcalib(ppg[[1]])
  ##e mv2 <- applycalib(res,mv1)
  ##e plot(mv1[,1],mv2[,1]-mv1[,1])

  if(missing(calib))
    {
      calib<-getPPGmasses()
    }
  object <- correctinternal(object,calib,error=error)
  theo<-NULL
  expd<-NULL
  exp<-object
  match<-getaccC(mass(exp),mass(calib),error=250,uniq=TRUE)
  expd <- c(expd,exp[match$plind])
  theo <- c(theo,calib[match$calind])

  mord <- order(theo)
  expd <- expd[mord]
  theo <- theo[mord]
  require(modreg)
                                        #calculate mass dependent error
  error <- (expd-theo)*1e6/theo
  ispl <- smooth.spline(theo,error)
  ispl <- calibspline(ispl,object,error,theo)
  return(ispl)
}


applyextcalib.massvector <- function(object,cS,...)
{
  ##t External Calbiration
  ##- Applys object of class calibspline to massvector to correct ofr measurment errors.
  ##- The error model stored in the calibspline are obtained by the function \code{getextcalib}
  ##d In case of \bold{external calibration} some sample spots are only dedicated
  ##d to calibration. Calibration samples which produces equidistant
  ##d peaks, which exact masses are known, can be used to precisely
  ##d estimate the mass dependent error function.
  ##+ object : massvector
  ##+ cS : calibspline
  ##v massvector : calibrated massvector
  ##sa applyextcalib.massvectorlist, applycalib.calibspline,getextcalib.massvector, getextcalib.massvectorlist
  ##r Gobom J, Mueller M, Egelhofer V, Theiss D, Lehrach H, Nordhoff E, 2002. A calibration method that simplifies and improves accurate determination of peptide molecular masses by MALDI-TOF MS. \emph{Anal Chem.} 74(15):3915-23.
  ##r Wolski \url{http://www.molgen.mpg.de/~wolski/mscalib}
  ##e data(ppg)
  ##e data(mv1)
  ##e res<-getextcalib(ppg[[1]],getPPGmasses(),error=200)
  ##e mv2<-applyextcalib(mv1,res)
  ##e compare(mv1,mv2)

  if(!inherits(cS,"calibspline"))
    {
      stop(as.character(substitute(cS)), " should be a calibspline but is of class: ", class(cS),"\n" )
    }
  #peaklist- array of peaks masses.
  #spline - spline to predict the error.
  error <- predict(cS,mass(object))
  error <- error$y

  masspred <- object[,1]/(1+error/1e6)
  object[,1] <- masspred
  return(object)
}


calibexternal.massvector <- function(object,ppg,calib,...)
  {
    ##t External Calbiration
    ##- Perfroms external calibration. Obtains the error model by calling \code{getextcalib}
    ##- and corrects the masses in the massvector for the errror.
    ##d In case of \bold{external calibration} some sample spots are only dedicated
    ##d to calibration. Calibration samples which produces equidistant
    ##d peaks, which exact masses are known, can be used to precisely
    ##d estimate the mass dependent error function.
    ##+ object : massvector
    ##+ ppg : either a massvector or massvectorlist with masses of the calibration sample (e.g. poly-(propylene glycol)).
    ##+ calib : a massvector with exact masses of the calibration sample (ppg).
    ##+ ... : further parameters.
    ##v massvector: calibrated massvector
    ##sa applyextcalib.massvectorlist, getextcalib.massvector, getextcalib.massvectorlist, calibexternal.massvector, applycalib.calibspline
    ##r Gobom J, Mueller M, Egelhofer V, Theiss D, Lehrach H, Nordhoff E, 2002. A calibration method that simplifies and improves accurate determination of peptide molecular masses by MALDI-TOF MS.\emph{Anal Chem.} 74(15):3915-23.
    ##r Wolski
    ##e data(ppg)
    ##e data(mv1)
    ##e mv2<-calibexternal(mv1,ppg)
    ##e compare(mv1,mv2)
    
    if(inherits(ppg,"massvector") | inherits(ppg,"massvectorlist"))
      {}
    else
      {
        warning("The second param should be a massvector or massvectorlist!\n")
        return(FALSE)
      }
    if(missing(calib))
      {
        calib<-getPPGmasses()
      }
    cS <- getextcalib(ppg,calib)
    res <- applyextcalib(object,cS)
    res
  }

#recalibration



recalibrate.massvector <- function(object,PQM=7,...)
  {
    ##t Precalibration
    ##d Precalibration method utilizes the knowledge that masses
    ##d of peptides are in equidistant spaced clusters. The wavelength of
    ##d the massesvector can be determined as described by
    ##d Wool. The comparision of the experimental wavelength with
    ##d the theoretical one, makes possible to find an affine function
    ##d that corrects the masses. Chemical noise in the spectra may hamper
    ##d the determination of mass list frequency. The package provides a
    ##d function to filter chemical noise.
    ##- Obtains the error and performs the calibration in one step.
    ##+ object : massvector
    ##+ PQM : Peak Quality Measure. Indicates how well the wavelenght of the massvector was determined.
    ##+ ... : further parameters
    ##w massvector : recalibrated massvector.
    ##r Wool A, Smilansky Z 2002. Precalibration of matrix-assisted laser desorption/ionization-time of flight spectra for peptide mass fingerprinting. {\em Proteomics.} 2(10):1365-73.
    ##r Wolski \url{http://www.molgen.mpg.de/~wolski/mscalib}
    ##e data(mv1)
    ##e mv2<-recalibrate(mv1)
    ##e plot(mv1[,1],mv1[,1]-mv2[,1],type="l")

    cal<-getrecalib(object)
    if(mget(cal,"PQM")>PQM)
      object<-applyrecalib(object,cal)
    return(object)
  }


applyrecalib.massvector<-function(object,calc,...)
  {
    ##t Precalibration
    ##d \bold{Precalibration} method utilizes the knowledge that masses
    ##d of peptides are in equidistant spaced clusters. The wavelength of
    ##d the massesvector can be determined as described by
    ##d Wool. The comparision of the experimental wavelength with
    ##d the theoretical one, makes possible to find an affine function
    ##d that corrects the masses. Chemical noise in the spectra may hamper
    ##d the determination of mass list frequency. The package provides a
    ##d function to filter chemical noise.
    ##- Uses the model of the measurment error stored in a calibrestat object
    ##- to correct the masses in the massvector.
    ##+ object : massvector.
    ##+ calc: object calibrestat class.
    ##+ ...: further arguments.
    ##w massvector : recalibrated massvector.
    ##sa getrecalib.massvector, getrecalib.massvectorlist, recalibrate.massvector, recalibrate.massvectorlist, calibrestat, calibrelist
    ##r Wool A, Smilansky Z 2002. Precalibration of matrix-assisted laser desorption/ionization-time of flight spectra for peptide mass fingerprinting.\emph{Proteomics.} 2(10):1365-73.
    ##r Wolski \url{http://www.molgen.mpg.de/~wolski/mscalib}
    ##e data(mv1)
    ##e res <- getrecalib(mv1)
    ##e plot(res)
    ##e mv2<-applyrecalib(mv1,res)
    ##e compare(mv1,mv2)

    
    peak<- object[,1]*(mget(calc,"Coeff.Slope")/1e6+1) + mget(calc,"Coeff.Intercept")
    object[,1] <- peak
    return(object)
  }

getrecalib.massvector<-function(object,plot=FALSE,...)
  {
    ##t Precalibration
    ##d \bold{Precalibration} method utilizes the knowledge that masses
    ##d of peptides are in equidistant spaced clusters. The wavelength of
    ##d the massesvector can be determined as described by
    ##d Wool. The comparision of the experimental wavelength with
    ##d the theoretical one, makes possible to find an affine function
    ##d that corrects the masses. Chemical noise in the spectra may hamper
    ##d the determination of mass list frequency. The package provides a
    ##d function to filter chemical noise.
    ##- Obtains the an error model using the wavelength analysis of the peaklist.
    ##+ mv : massvector
    ##+ plot: \code{TRUE} - A graphic showing the \eqn{F(\omega)}{F(omega)} function. default \code{FALSE}
    ##v calibrestat : object of class calibrestat.
    ##sa applyrecalib.massvector, massvector, calibrestat, wsFilter.massvector
    ##r Wool A, Smilansky Z 2002. Precalibration of matrix-assisted laser desorption/ionization-time of flight spectra for peptide mass fingerprinting.\emph{ Proteomics.} 2(10):1365-73.
    ##r Wolski \url{http://www.molgen.mpg.de/~wolski/mscalib}
    ##e data(mv1)
    ##e res <- getrecalib(mv1)
    ##e print(res)
    ##e as.vector(res)
    ##e summary(res)
    ##e image(res)
    ##e plot(res)


    mv<-mass(object)
                                        #returns calibration constants obtained by the wools smilanski method.
    if(length(object)<3)
      {
        res<-calibrestat(info(object))
        setParms(res)<-list(Coeff.Intercept=0)
        setParms(res)<-list(Coeff.Slope=0)
        setParms(res)<-list(lengthmv = length(object))
        setParms(res)<-list(tcoor=mget(object,"tcoor"))
        setParms(res) <- list(PQM = 0)
        return(res)
      }
                                        #   lambda<-seq(1,1.001,0.000001)
    lambda<-seq(0.999495,1.001495,0.000005)
    omega<-2*pi/lambda
    test<-mv%*%t(omega)
    testsin<-sin(test)
    testcos<-cos(test)
    colsumsin <-apply(testsin,2,sum)
    colsumcos <-apply(testcos,2,sum)
    sumcol<-sqrt(colsumsin^2+colsumcos^2)
    sumcol <- cbind(lambda,sumcol)
                                        #determine the maximum and get its index
    mmax <- max(sumcol[,2])
    index<-sumcol[,2] == mmax
    if(FALSE)
      {
        return(sumcol)
      }
                                        #temporary only for testing purposes.
    if(plot){
      len<-length(mv)
      main <- paste("Massvector length : ", len, sep="")
      plot(sumcol,type="l",ylab="Amplitude",xlab="Wavelength",xlim=c(0.9995,1.0015),las=2, main = main)
      points(sumcol[,1][index],sumcol[,2][index],col=2,pch="*")
      stat<-FullWidthatHalfMaximum(sumcol)
      abline(h=stat$hm)
      abline(v=stat$width)
                                        #end temporary
    }
                                        #by which wavelength the maximum occure
    lambdamax<-sumcol[,1][index]
                                        #calculate the phase shift.
    phimax<-atan(colsumsin[index]/colsumcos[index])
                                        #we now see that the peak centers lie on the line
                                        # M= lambdamax*N + bmaxx
    bmax<-lambdamax*phimax/(2*pi)
                                        #to coorect the peaklist apply
                                        #peak<-(peaklist-bzero)/alpha
    alpha <- 1.000495/lambdamax
    res<-calibrestat(info(object))
    tmp<-c(bmax,alpha)
    names(tmp)<- c("Intercept","Slope")
    setParms(res)<-list(Coeff.Intercept=- (tmp[1]*tmp[2]))
    setParms(res)<-list(Coeff.Slope=(tmp[2]*1e6-1e6))
    setParms(res)<-list(lengthmv = length(object))
    setParms(res)<-list(tcoor=mget(object,"tcoor"))
    tmp <- FullWidthatHalfMaximum(sumcol)
    setParms(res) <- list(PQM = tmp$PQM)
                                        #    setParms(res)<-list(stat=tmp)
    return(res)
  }


##############################################
## affine calibration
#


applyintcalib.massvector <- function(object,cal,...)
  {
    ##t Internal Calibration
    ##- Corrects the massvector for the error model stored in calibintstat object.
    ##d \bold{Internal calibration} aligns masses of
    ##d peaks to known masses and determines by linear regression a affine
    ##d function that describing the relative error. The internal
    ##d correction fails when no calibration peaks can be found.
    ##+ object : massvector
    ##+ cal : object of class calibintstat
    ##+ ... : further parameters.
    ##v massvector : calibrated massvector. 
    ##sa applycalib.calibintstat,getintcalib.massvector, correctinternal.massvector
    ##r Wolski \url{http://www.molgen.mpg.de/~wolski/mscalib}
    ##e data(mv1)
    ##e data(cal)
    ##e res<-getintcalib(mv1,cal,error=300)
    ##e mv2<- applyintcalib(mv1,res)
    ##e plot(mv1[,1],mv2[,1]-mv1[,1])
    if(length(object)==0)
      {
        return(object)
      }
    if(inherits(cal,"lm"))
      {
        errp <- predict(cal,data.frame(masstheo=mass(object) ) )
        if(mget(cal,"ppm"))
          {
            object[,1] <- object[,1]/(1 - errp/1e6)
          }
        else
          {
            object[,1] <- object[,1] + errp
          }
      }
    object
  }


correctinternal.massvector<-function(object,calib,error=500,uniq=FALSE,ppm=TRUE,...)
  {
    ##t Internal Calibration
    ##- Corrects the masses of the massvector. It first obtains the model of the
    ##- measrument error by calling \code{getintcalib.massvector}. It than corrects the masses
    ##- by a call to \code{applyintcalib.massvector}.
    ##d \bold{Internal calibration} aligns masses of
    ##d peaks to known masses and determines by linear regression a affine
    ##d function that describing the relative error. The internal
    ##d correction fails when no calibration peaks can be found.
    ##+ object : massvector
    ##+ calib : massvector with calibration masses
    ##v massvector : calibrated massvector. 
    ##sa getintcalib.massvector, calibintstat, applycalib.calibintstat
    ##r Wolski \url{http://www.molgen.mpg.de/~wolski/mscalib}
    ##e data(mv1)
    ##e data(cal)
    ##e res <- correctinternal(mv1,cal,error=200)
    
    tmp <- getintcalib(object,calib,error=error,uniq=uniq,ppm=ppm)
    tmp <- applyintcalib(object,tmp,ppm=ppm)
    tmp
  }


getintcalib.massvector <- function(object,calib,error=500,uniq=FALSE,ppm=TRUE,...)
  {
    ##t Internal Calbiration.
    ##- Obtains error model by alingning masses in massvector to known masses (calibration list).
    ##d \bold{Internal calibration} aligns masses of
    ##d peaks to known masses and determines by linear regression a affine
    ##d function that describing the relative error. The internal
    ##d correction fails when no calibration peaks can be found.
    ##d 
    ##+ object : massvector
    ##+ calib : massvector with calibration masses
    ##+ error : assumed measurment error.
    ##+ uniq : \code{TRUE}- use only mass closest to calibration mass. \code{FALSE}- use all masses closer to the calibration mass then given error.
    ##+ ppm : \code{TRUE}- describe the error as relative error. \code{FALSE}- describe the error as absolute error.
    ##v calibintstat : object of class calibintstat.
    ##sa applyintcalib.massvector, getintcalib.massvector, correctinternal.massvectorlist, calibintstat
    ##r Wolski \url{http://www.molgen.mpg.de/~wolski/mscalib}
                                        #peaklist - peaklist.
                                        #calib - list with calibrants.
    massmv<-mass(object)                    #use mass(.. because you still dont shure about implementation.
    masscalib<-mass(calib)
    match<-NULL
    match<-getaccC(massmv,masscalib,error=error,ppm=ppm,uniq=uniq)
    if(length(match$plind)>1)
      {
                                        #calculate mass dependent error function.
                                        #get peaklist matching
        smallerrpl<-massmv[match$plind]
                                        #get calibrants matching
        masstheo<-masscalib[match$calind]
        if(ppm)
          {
            err  <-  (masstheo - smallerrpl) * 1e6 /masstheo
          }
        else
          {
            err  <-  masstheo - smallerrpl
          }
                                        #finally calculate error
                                        #distinguishing 2 cases. if peaks quite close to each other only correct for offset.
        if(abs(diff(range(masstheo)))<200)
          {
            err<-mean(err)
            masstheo <- mean(masstheo)
            mymod <- lm(err~masstheo)
          }else{
            mymod <- lm(err~masstheo)
          }
        coraf <- calibintstat(info(object),mymod)
        stat <- c(mean(abs(err)),sqrt(var(err)))
        names(stat) <- c("mean","stdv")
        setParms(coraf) <- list(Coeff.Intercept= mymod$coefficients[1],
                                Coeff.Slope =  ifelse(ppm,mymod$coefficients[2]*1e4,mymod$coefficients[2]*1e6) ,
                                lengthmv = length(object),
                                nrmatch = length(match$plind),
                                ppm=ppm,
                                tcoor = mget(object,"tcoor")
                                )
      }
    else
      {
        err <- rep(0,10)
        masstheo <- 1:10
        mymod <- lm(err~masstheo)
        coraf <- calibintstat(info(object),mymod) 
        setParms(coraf) <- list( Coeff.Intercept = mymod$coefficients[1] ,
                                Coeff.Slope =  ifelse(ppm,mymod$coefficients[2]*1e4,mymod$coefficients[2]*1e6),
                                lengthmv = length(object),
                                nrmatch = 0 ,
                                ppm=ppm,
                                tcoor = mget(object,"tcoor")
                               )
      }
    return(coraf)
  }


################################################
# gamasses
#


gamasses.massvector <- function(object,accur = 0.1,abund = 50,...)
  {
    ##t Abundant masses
    ##- Determines abundant masses in a massvector.
    ##d Abundant masses are masses that occur in a large
    ##d fraction of the massvectors. Typical abundant mass
    ##d are derived from tryptic autoproteolysis products.
    ##d Abundant masses can also often be assigned to keratin
    ##d isoforms (human hair- and skin proteins).
    ##d Many of the abundant masses cannot be assigned to any protein.
    ##d Abundant masses can be used for calibration.
    ##d Removing them may increase the identification specificity.
    ##sa gamasses.massvectorlist
    ##+ object : massvector
    ##+ accur : measurment accuracy in m/z
    ##+ abund : how many times a mass have to occur in a mass bin to be considered as an abundant mass.
    ##+ ... : further parameters.
    ##v massvector : massvector with abundant masses.
    ##r Wolski 
    ##e data(mvl)
    ##e mv<-unlist(mvl) 
    ##e res<-gamasses(mv,abund=30,accur=0.3)
    ##e plot(res)
    
    dat<-object[,1]
    attributes(dat)<-NULL
    mhist<-list(NULL)
    mhist[[1]] <-hist(dat,breaks=seq(min(dat)-accur/2,max(dat)+1.5*accur,accur),plot=FALSE)#,main=main,xlab=xlab,xlim=xlim)
    mhist[[2]]<-hist(dat,breaks=seq(min(dat)-accur,max(dat)+accur,accur),plot=FALSE)
    wh1<-which(mhist[[1]]$counts > abund)
    wh2<-which(mhist[[2]]$counts > abund)
    if(length(wh1)==0 & length(wh2)==0)
      {
        res<-matrix(ncol=3)
        colnames(res)<-c("mass","number","wm")
        return(massvector(paste(info(object),"abundant",sep="_"),res))
      }
                                        #    cat("wh1 : ", wh1 ," \nwh2 : ",wh2,"\n")
                                        #bestimme die maxima.
    z1<- rep(0,length(mhist[[1]]$mids))
    z2<- rep(0,length(mhist[[2]]$mids))
                                        #    cat("z1 " , z1 , " z2 " , z2 , "\n" )
                                        #bernehme nur die die hufig sind.
    z1[wh1]<- mhist[[1]]$counts[wh1]
    z2[wh2]<- mhist[[2]]$counts[wh2]
                                        #find the peaks.returns the indexes    
    p1 <- peaks(z1,max=TRUE)
    p2 <- peaks(z2,max=TRUE)
    names(p1)<-rep(1,length(p1))        #name the indexes
    names(p2)<-rep(2,length(p2))        #name the indexes
                                        #Determin which indexes are close to each other.
    p12 <- sort(c(p1,p2))
                                        #if difference between peaks are 0 or 1 than they are close to each other.
    dp <- diff(p12)
                                        #find the indexes of this cases
    ipl <- which(dp==1|dp==0)
    ipu <- ipl+1
                                        #now i have the indices of the peaks in array p12.
                                        #Array p12 by themselves gives me the indices of the peaks.
    p12close <- cbind(p12[ipl],p12[ipu])
    p12closeN <- cbind(as.numeric(names(p12)[ipl]),as.numeric(names(p12)[ipu]))
                                        #now in each row are the indices of the peaks that are close to each other.
                                        #using this indices i can retrieve the mids and the counts
                                        #for this bins out of the histogram and calculate a weighted
                                        #average.
    wm <- NULL

    if(length(p12close)>0)
      {
        for(x in 1:length(p12close[,1]))
          {
            w1 <- mhist[[ p12closeN[x,1] ]]$counts[p12close[x,1]]
            m1 <- mhist[[ p12closeN[x,1] ]]$mids[p12close[x,1]]
            w2 <- mhist[[ p12closeN[x,2] ]]$counts[p12close[x,2]]
            m2 <- mhist[[ p12closeN[x,2] ]]$mids[p12close[x,2]]
            wm <- c( wm , (w1*m1+w2*m2)/(w1+w2) )
          }
      }

    if(length(c(ipl,ipu)>0))
      {
        uniquepeaks <- p12[-c(ipl,ipu)]
      }
    else
      {
        uniquepeaks <- p12
      }
    uN <- as.numeric(names(uniquepeaks))

    if(length(uniquepeaks)>0)
      {
        for(x in 1:length(uniquepeaks))
          {
            wm<-c(wm,mhist[[ uN[x] ]]$mids[uniquepeaks[x]]) 
          }
      }
    rm(mhist)
    res <- NULL
    num<-NULL
                                        #calculating exact mass
    for(x in wm)
      {
        #print(x)
        res <- c(res, mean( dat[(x+accur*0.8) > dat & (x-accur*0.8) < dat ]))
        num<-c(num,length(dat[(x+accur*0.8) > dat & (x-accur*0.8) < dat ]))
      }
                                        #    print(paste(info(object),"abundant",sep="_"))

    res<-cbind(res,num,wm)
    colnames(res)<-c("mass","number","wm")
    return(massvector(paste(info(object),"abundant",sep="_"),res))
  }


diffFilter.massvector<-function(object,listofdiffs,higher=TRUE,error=0.05,uniq=TRUE,prune=TRUE,...)
  {
    ##t Abundant Differences
    ##- Removes mass differences from the massvector.
    ##d Removes one of the masses contributing to a mass difference given in the list of differences.
    ##d Can be used if a variable modification are present in the massvector but can not be considered by the identification software. It also can be used to return the modified masses.
    ##+ object : massvector
    ##+ listoffdiffs : massvector with mass differences
    ##+ higher : logical;\code{TRUE} - remove higher mass, \code{FALSE} = remove lower mass.
    ##+ prune : logical;\code{TRUE} - remove mass, \code{FALSE} = return modified mass.
    ##+ error : How much the differences can diviate from the differences given in listofdiffs.
    ##+ ... : further parameters.
    ##v massvector : filtered massvector.
    ##sa getdiff.massvector, getdiff.massvectorlist, diffFilter.massvectorlist
    ##r Wolski
    ##e data(mv1)
    ##e res<-getdiff(mv1,range=c(0,100))
    ##e diffFilter(mv1,res,higher=TRUE)
    ##e diffFilter(mv1,res,higher=FALSE)
    if(!inherits(listofdiffs,"massvector"))
      stop("second arg has to be a massvector\n")
    if(length(object)<=1)
      {
        return(object)
      }
    pl <- object[,1]
    names(pl)<-1:length(pl)
    pl <- sort(pl) # should be sorted anyway.
    res <- NULL
    ldiff<-listofdiffs[,1] # get the difference masses.
    for(x in 1:(length(pl)-1))
      {
        diffpl <- sort(diff(pl,x)) # computes differences
        tmp <- getaccC(diffpl,ldiff,error=error,ppm=FALSE,uniq=uniq)
        ind <- as.numeric(names(diffpl))[tmp$plind]
        if(higher)
          {
            res <- c(res,ind)
          }
        else
          {
            res<-c(res , ind-x )
          }
      }
    if(prune)
      {
        if(length(res)>0)
          {
            #cat("length : ",length(object),"\n")
            #print(res)
            return(object[-res,])
          }
        else
          return(object)
      }
    else
      {
        if(length(res)>0)
          return(object[res,])
        else
          return(object[NULL,])
      }
  }

getdiff.massvector<-function(object,range=c(0,100),...)
  {
    ##t Massdifferences
    ##- Computes mass differences in the object.
    ##d Removes one of the masses contributing to a mass difference given in the list of diffs.
    ##d Can be used if a variable modification are present in the massvector but can not be considered by the identification software.
    ##+ object : massvector
    ##+ listoffdiffs : massvector with mass differences
    ##+ higher : \code{TRUE} - remove higher mass, \code{FALSE} = remove lower mass.
    ##+ error : How much the differences can diviate from the differences given in listofdiffs
    ##v massvector : filtered massvector.
    ##r Wolski

    rrr<-range
    res <- NULL
    rarea <- NULL
    if(!length(object)>1)
      {
        res<-massvector(info(object),NULL)
        setParms(res)<-list(tcoor=mget(object,"tcoor"))
        return(res)
      }
    else
      {
        m <- mass(object)
        a<-object[,2]
        for(x in 1:(length(m)-1))
          {
            diffpl<-m - m[x]
            ratarea<-a / a[x]
            names(diffpl)<-NULL
            names(ratarea)<-NULL
            res <- c(res,diffpl[diffpl>rrr[1] & diffpl<rrr[2]])
            rarea <- c(rarea,ratarea[diffpl>rrr[1] & diffpl<rrr[2]])
          }

      }
    if(length(res)==0)
      {
        res<-massvector(info(object),NULL)
        setParms(res)<-list(tcoor=mget(object,"tcoor"))
        return(res)
      }
    res<-cbind(res,rarea)
    colnames(res)<-c("massd","arear")
    res<-massvector(info(object),res)
    setParms(res)<-list(tcoor=mget(object,"tcoor"))
    return(res)
  }

writeF.massvector<-function(object,path,file=info(object),ext="txt",...)
  {
    ##t Write massvector
    ##- Write massvector to File.
    ##d The read and write functions for different peak-list formats are not provided by the package. This is because
    ##d there are a oodles of different formats. I will try to collect read-write functions for as many as possible peak-list format's in an add on package
    ##d which you can find at \url{http://www.molgen.mpg.de/~wolski/mscalib/IO/}.
    ##+ object : massvector
    ##+ path : path to folder.
    ##+ file : file name. defualt info(object)
    ##+ ext : file extension.
    ##sa writeF.massvectorlist, readF.massvector
    ##e data(mv1)
    ##e writeF(mv1,".") # writes the file in the home directory.
    ##e readF(massvector(info(mv1)),".")
    ##e file.remove(paste(info(mv1),".txt",sep=""))
    if(missing(path))
      {
        path<-"."
      }
    
    filep <- file.path(path,paste(file,".",ext,sep=""),fsep = .Platform$file.sep)
    con <- file(filep, "w")  # open an output file connection
    firstline<-paste(">",file,":",join(mget(object,"tcoor"),sep=","),sep="")
    writeLines(firstline,con=con)
    close(con)
    write.table(object,file = filep, append = TRUE, quote = FALSE, sep = "\t",
            eol = "\n", na = "NA", dec = ".", row.names = FALSE,
            col.names = FALSE, qmethod = c("escape", "double"))
  }


readF.massvector<-function(object,path,file=info(object),ext="txt",...)
  {
    ##t Read Massvector
    ##- Reads massvector written with the function \code{writeF.massvector}
    ##d The read and write functions for all the different peak-list formats are not provided by the package. This is because
    ##d there are oodles of different formats. I will try to collect read-write functions for as many as possible peak-list format's in an add on package
    ##d which you can find at \url{http://www.molgen.mpg.de/~wolski/mscalib/IO/}.
    ##d The file format of the file to be read must be:
    ##d \code{>0_A1_1SRef:1,1}
    ##d \code{842.4257236	650.66}
    ##d \code{987.3931319	180.3}
    ##d \code{...  ...}
    ##d The char between > and : is read into the info field and used as a key in the massvectorlist. So it must be unique.
    ##d After the double colon the coordinates have to written.
    ##d The first column are the masses. The second column are the intensities.
    ##+ object : object of class massvector. Use constructor \code{massvector()}
    ##+ path : path to file.
    ##+ file : file name; default : info(object).
    ##+ ext : file extension; default : txt.
    ##sa writeF.massvector
    ##e data(mv1)
    ##e writeF(mv1,".")
    ##e readF(massvector(info(mv1)),".")
    ##e file.remove(paste(info(mv1),".txt",sep=""))
    filep <- file.path(path,paste(file,".",ext,sep="") , fsep = .Platform$file.sep)
    con<-file(filep,"r")
    res <- readLines(con=filep,n=-1)
    close(con)
    res1 <- unlist(strsplit(res[1],":"))
    object<-info(object,sub(">","",res1[1]))
    setParms(object) <- list(tcoor= unlist(strsplit(res1[2],",")))
    intern<-function(x){as.numeric(unlist(strsplit(x,"\t")))}
    res <- t(sapply(res[2:length(res)],intern))
    res<-res[order(res[,1]),]
    rownames(res)<-1:length(res[,1])
    object <- peaks(object,res)
    return(object)
  }



residuals.massvector<-function(object,fmass,error=250,ppm=TRUE,uniq=TRUE,...)
{
  ##t Residues of 2 massvectors.
  ##- Computes Residues of 2 massvectors.
  ##d Computes the mass differences between aligned masses of two massvectorlists.
  ##+ object : object of class massvector. Use constructor \code{massvector()}
  ##+ fmass : object of class massvector.
  ##+ error : the assumed difference between the aligned peaks either in ppm or in Da.
  ##+ ppm : logical; if TRUE then the error must be specified in ppm.
  ##sa getaccC,compare
  ##e data(mv1)
  ##e data(mv2)
  ##e plot(resid(mv1,mv2,error=250,uniq=TRUE))

  if(!inherits(fmass,"massvector"))
    stop("second arg must be massvector to\n")
  mmatch<-getaccC(mass(object),mass(fmass),error=error,ppm=ppm,uniq=uniq)
  
  if(length(mmatch$plind)>0)
    {
      ttmp <- object[mmatch$plind,1]-fmass[mmatch$calind,1]
      res<-cbind(fmass[mmatch$calind,1], ttmp)
      return(massvector(info(object),res))
    }
  else
    {
      return(object[NULL,])
    }
}

                                        #Copyright 2004, W. Wolski, all rights reserved.
                                        #massvector list.
massvectorlist <- function(experiment,data,project,...)
  {
    ##t Constructor
    ##- constructor of object massvectorlist (extends list).
    ##+ experiment : name of experiment
    ##+ data : list of massvectors
    ##+ project : name of project.
    ##v massvectorlist : object of class massvector.
    ##e # testing constructor.
    ##e massvectorlist("my1experiment")
    ##e data(mvl)
    ##e massvectorlist("my2experiment",mvl,"hello project")
    ##e plot(mvl)
    ##e summary(mvl)
    ##e hist(mvl)
    ##e hist(mvl)
    ##e image(mvl,what="lengthmv")
    ##e mvl2<-mvl[1:100]
    ##e plot(mvl2)
    ##e summary(mvl2)
    ##e hist(mvl2)
    ##e image(mvl2,what="lengthmv")
    ##e #testing assingments
    ##e mvl2[[11]]<-mvl2[[1]]
    ##e plot(mvl2[[11]],mvl2[[1]])
    ##e image(mvl2[[11]],mvl2[[1]])
    ##e #make one massvector out of the peaklist
    ##e tt<-unlist(mvl)
    ##e plot(tt)
    
    if(!missing(data))
      {
        res <- data
      }
    else
      res <- list()
    attr(res,"allow") <- c("data","experiment","project")
    if(!missing(project))
      attr(res,"project") <- project #project.
    if(!missing(experiment))
      attr(res,"experiment") <- experiment # the big project may consists of several experiments
    class(res) <- c("massvectorlist","mlist","list","myobj")
    return(res)
  }

c.massvectorlist <- function(mvl,...)
  {
    ##t Combine 
    ##- Combines Massvectorlist into one Massvectorlist.
    ##d It does not check if the massvectors in the list are unique.
    ##+ mvl : massvector
    ##+ ... : massvectorlists to be concatenated.
    ##sa c
    ##v massvectorlist : massvectorlist
    ##e data(mvl)
    ##e mvl2<-c(mvl[1:100],mvl[1:50])
    ##e mvl2
    
    tmp <- list(...)
    for(x in tmp)
      {
        if(!inherits(x,"massvectorlist"))
          stop("All Arguments have to be a massvectorlist!!!")
      }
    res <- NextMethod("c")
    res<-massvectorlist(info(mvl),res,mget(mvl,"project"))
    res
  }

  
print.massvectorlist<-function(x,...)
{
  ##t Print massvectorlist
  ##- `print' prints its argument and returns it invisibly (via `invisible(x)')
  ##+ x : massvectorlist
  ##sa \link[base]{print}
  ##e data(mvl)
  ##e print(mvl)
  exp <- mget(x,"experiment")
  print(paste("Experiment : ", ifelse(is.null(exp),"()",paste("(",exp,")",sep="")),sep=""),quote=FALSE,...)
  proj<-mget(x,"project")
  print(paste("Project    : ", ifelse(is.null(proj),"()",paste("(",proj,")",sep="")),sep=""),quote=FALSE,...)
  print(paste("Nr (mv)    : ", paste("(",length(x),")",sep=""),sep=""),quote=FALSE,...)
  invisible(x)
}





as.list.massvectorlist<-function(x,...)
  {
    ##t List
    ##- Turns the massvectorlist into a list
    ##+ x : massvectorlist.
    ##+ ... : further parameter.
    ##e data(mvl)
    ##e res<-as.list(mvl)
    ##e class(res)
    al <- attr(x,"allow")
    for(u in al)
      {
        attr(x,u)<-NULL
      }
    class(x)<-"list"
    x
  }



project.massvectorlist <- function(object,project,...)
  {
    ##t Acces
    ##- Access to the project field of the massvectorlist.
    ##d Can be used for setting or getting the project field.
    ##+ object : massvectorlist
    ##+ project : info character. If missing function returns the current info. If not missing function returns massvector with new project field content.
    ##+ ... : further arguments
    ##sa info.mlist, experiment.mlist
    ##e data(mvl)
    ##e project(mvl)
    ##e mvl <- project(mvl,"newprojectname")

    if(missing(project))
      return(mget(object,"project"))
    else
      {
        setParms(object)<-list(project=project)
      }
    object
  }

histLengths <- function(x,main=info(x),xlab="length of massvectors",...)
  {
    ##t Histograms
    ##- Histogram of the massvector lengths.
    ##+ x : massvectorlist
    ##+ main : info(x)
    ##+ xlab : length of massvectors.
    ##+ ... : further parameters.
    ##sa \link[base]{hist}
    ll <- lapply(x,length)
    res <- hist.default(unlist(ll),xlab="massvector length",main=main,border=1,...)
  }


hist.massvectorlist<-function(x,accur = 0.1, main=info(mvl) ,xlab="m/z",xlim=c(700,4500),add=FALSE,col=1,...)
  {
    ##t Histograms
    ##- Histogram of mass frequencies in the massvectorlist.
    ##+ x : massvectorlist
    ##+ accur : bin width for plotting mass frequencies.
    ##+ col : color of the histogram,
    ##+ main : title of graph
    ##+ xlim : the range to be encompassed by the x axis.
    ##+ xlab : a title for the x axis.
    ##+ add : logical; If \code{TRUE} add to already existing plot.
    ##+ ... : further parameters.
    ##sa hist.massvector
    ##e data(mvl)
    ##e hist(mvl)
    mvl<-x
    rm(x)
    dat<-unlist(mvl)
                                        #attributes(dat)<-NULL
    mhist<-list(NULL)
                                        #assign indices of bins with high peak abundance
    mhist[[1]] <- hist(mass(dat),breaks=seq(min(mass(dat))-accur/2,max(mass(dat))+1.5*accur,accur),plot=TRUE,main=main,xlab=xlab,xlim=xlim,add=add,col=col,border=col,...)
    mhist[[2]] <- hist(mass(dat),breaks=seq(min(mass(dat))-accur,max(mass(dat))+accur,accur),add=TRUE,col=col,border=col)
  }

plot.massvectorlist <-function(x, main=info(mvl) ,xlab="m/z",xlim=c(700,4500),add=FALSE,col=1,cex=0.5,...)
{
  ##t Massvectorlist Plotting
  ##- Plots masses (m/z) in massvector against sample in massvectorlist.
  ##+ x : massvectorlist
  ##+ main : an overall title for the plot.
  ##+ xlab : a title for the x axis.
  ##+ xlim : the range to be encompassed by the x axis.
  ##+ add : \code{TRUE} - masses of a new massvectorlist are added to existing plot.
  ##+ col : color of the points denoting the masses.
  ##+ ... : graphical parameters can be given as arguments to `plot'.
  ##sa hist.massvectorlist, plot.massvector, image.mlist, plot, par
  ##e data(mvl)
  ##e plot(mvl,col= 3 )
  ##e tt<- gamasses(mvl,abund=50)
  ##e plot(mvl,col=1, xlim=c(tt[1,1]-0.4,tt[1,1]+0.4))
  mvl<-x
  rm(x)
  par(bg="lightgray")
  if(add)
    {
      for( x in 1:length(mvl))
        {
          points(mass(mvl[[x]]), rep(x,length(mass(mvl[[x]]))),pch=15,cex=cex,col=col)
        }
    }
  else
    {
      
      allm  <- unlist(mvl)
      mmin  <- min(allm)[1]
      mmax <- max(allm)[1]
      lod<-mvl
      if(is.null(xlim))
        {

          plot.default(mass(lod[[1]]), rep(1,length(mass(lod[[1]]))) , xlim = c(mmin,mmax) , ylim = c(1,length(lod)),pch=15,cex=cex,xlab=xlab,ylab="sample",main=main,col=col,...)
        }
      else
        {
          plot.default(mass(lod[[1]]), rep(1,length(mass(lod[[1]]))) , xlim = xlim , ylim = c(1,length(lod)) , pch=15, cex=cex , xlab=xlab , ylab="sample" , main=main,col=col,...)
        }
      for(x in 2:length(lod))
        {
          points(mass(lod[[x]]), rep(x,length(mass(lod[[x]]))),pch=15,cex=cex)
        }
      abline( h = seq(0,length(lod),10) , lty=3)
    }
}


"[.massvectorlist"<-function(mvl,i)
  {
    ##t Extract Parts of an Massvectorlist.
    ##- The massvectorlist extends list. The massvectors in the massvectorlist can therefore be accessed like list elements.
    ##+ mvl : massvectorlist
    ##+ i : indices of massvectors to extract
    ##v massvectorlist : massvectorlist
    ##sa [<-.massvectorlist, [<-.list
    ##e data(mvl)
    ##e mvl<-mvl[1:10] # returns a massvectorlist of length 10
    ##e mvl
    ##e class(mvl)

 
    tmp<-NextMethod("[");
    res<- massvectorlist(experiment(mvl),tmp,project(mvl))
    res
  }


"[<-.massvectorlist"<-function(mvl,i,value)
  {
    ##t Replace Parts of Massvectorlist
    ##- The massvectorlist extends list. The massvectors in the list can therefore be accessed like list elements.
    ##+ mvl : massvectorlist.
    ##+ i : elements to replace.
    ##+ value : replace by value.
    ##sa [.massvector,[.list
    ##v massvectorlist : massvectorlist
    ##e data(mvl)
    ##e mvl2<-mvl
    ##e mvl2[11:20]<-mvl[1:10]
    ##e compare(mvl[[11]],mvl[[1]])
    
    if(!inherits(value,"massvectorlist"))
      {
        stop(as.character(substitute(value)),"not a massvectorlist")
      }
    NextMethod("[<-")
  }

mvFilter.massvectorlist <-  function(object,fmass,match=FALSE,error=250,ppm=TRUE,uniq=FALSE,...)
{
  ##t Filtering Massvector
  ##- Filters massvector for masses given in a second massvector.
  ##+ object : massvectorlist
  ##+ abund : massvector
  ##+ error : mesurment error.
  ##+ ppm : given either in ppm (\code{TRUE}) or as absolut error (\code{FALSE}).
  ##+ match : logical;\code{TRUE} - than returns masses matching to the masses in massvector abundant, \code{FALSE} - returns masses not matching.
  ##+ uniq : logical; \code{FALSE} - returns or removes all masses in the range given by error. \code{TRUE} - returns ore removes only the closest mass.
  ##+ ... : further parameters.
  ##v massvectorlist : with massvectors with matching or not matching masses.
  ##sa mvFilter.massvector
  ##e data(mvl)
  ##e mvl<-mvFilter(mvl,mvl[[1]],match=FALSE,error=250)
  ##e length(mvl[[1]])
   res<-lapply(object,mvFilter,fmass,match=match,error=error,ppm=ppm,uniq=uniq)
   res<-  massvectorlist(experiment(object),res,project(object))
   res
}
  




unlist.massvectorlist<-function(x,...)
  {
    ##t Flatten massvectorlist
    ##-  Given a list structure `x', `unlist' simplifies it to produce a
    ##-  massvector which contains all the atomic components which occur in x.
    ##+ x : massvectorlist.
    ##+ ... : further arguments.
    ##sa \link[base]{unlist}
    ##v massvector : massvector
    ##e data(mvl)
    ##e amv <- unlist(mvl)
    ##e length(amv)
                                        #     tmp<-NextMethod("unlist")
    tmp<-NULL
    for(u in x)
      {
        if(length(u)!=0)
          tmp<-rbind(tmp,u)
      }
    tmp<-massvector(info(x),tmp)
    tmp
  }



"[[<-.massvectorlist"<-function(x,i,value)
  {
    ##t Replace Parts of an Object
    ##- Replace a massvector in the massvectorlist with a different one.
    ##+ x : massvectorlist
    ##+ i : index or name (info) of massvector to replace
    ##+ value : massvector
    ##e data(mvl)
    ##e data(mv1)
    ##e mvl[[10]]<-mv1
    
    
    if(!inherits(value,"massvector"))
      stop("only massvectors allowed to assing")
    x <- NextMethod("[[<-")
    if(is.numeric(i))
      names(x)[i]<-mget(value,"info")
    x
  }

  
wsFilter.massvectorlist<-function(object,mdist=0.25,fraction=0.2, peptides=TRUE,... )
  {
    ##t Smilanski Filtering
    ##- Removes chemical noise from massvectors in the massvectorlist (if \code{peptides} argument \code{TRUE}) or returns it.
    ##d Chemical noise can be removed from the peptide mass lists
    ##d due to the strong clustering of mono-isotopic peptide
    ##d peaks. Following the distance measure and filtering
    ##d method proposed by Wool Smilanski we developed an algorithm to
    ##d classify masses as peptide and non-peptide. The algorithm is based
    ##d on a modified distance measure and hierarchical clustering of all
    ##d intra massvector distances.
    ##+ object : massvectorlist.
    ##+ mdist : Minimal distance to branch to prune. The unit of the distance is Dalton.
    ##+ fraction : Maximal fraction (nr masses in branch)/(length of massvector) of branche to be prune.
    ##+ peptides : logical; \code{TRUE} - returns peptides, \code{FALSE} - returns chemical noise.
    ##+ ... : further parameters.
    ##v massvectorlist :  Returns a massvectorlist where the massvectors either contain the peptides or the non-peptides.
    ##r Wool A, Smilansky Z 2002. Precalibration of matrix-assisted laser desorption/ionization-time of flight spectra for peptide mass fingerprinting. Proteomics. 2(10):1365-73.
    ##sa wsiFilter.massvector, wsFilter.massvector
    ##e data(mvl)
    ##e res <- wsFilter(mvl,peptides = TRUE)
    ##e plot(res)
    ##e res2 <- wsFilter(mvl,peptides = FALSE)
    ##e plot(res2,col=2,add=TRUE)
    ##e image(res2,what="lengthmv")
    ##e hist(res)
    ##e hist(res2,col=2,add=TRUE)
    res<-lapply(object,wsFilter,mdist=mdist,fraction=fraction,peptides=peptides)
    res<- massvectorlist(experiment(object),res,project(object))
    res
  }


#----------------------------------------------
##
#calculates the spline function from the return value of getListforCalib
#it returns an spline object this spline object can be used for calibration of ppg spectra by the function
#returns a list with the spline prediction used for calibration the theoretical time and the averaged time
##

getextcalib.massvectorlist <- function(object,calib,error=250,...)
{
  ##t External Error Model
  ##- Returns the error model obtained from the calibration sample.
  ##+ object : massvectorlist with calibration sample (ppg) masses
  ##+ calib : massvector with calibration masses (ppg)
  ##+ error : relative measurment error in ppm. (peaks in this range are assumed as matching)
  ##+ ... : further parameters.
  ##v calibspline : can be used to calibrate peaklists
  ##sa calibspline, applyextcalib.massvectorlist ,getextcalib.massvector
  ##r Gobom J, Mueller M, Egelhofer V, Theiss D, Lehrach H, Nordhoff E, 2002. A calibration method that simplifies and improves accurate determination of peptide molecular masses by MALDI-TOF MS. Anal Chem. 74(15):3915-23.
  ##r Wolski \url{http://www.molgen.mpg.de/~wolski/mscalib}

  if(missing(calib))
    {
      calib <- getPPGmasses()
    }
                                        #ppg<-getPPGmasses()
  ppg<-calib
  theo<-NULL
  expd<-NULL
  for(x in 1:length(object))
     {

      exp <- object[[x]]      
      exp <- correctinternal(exp,calib,error=error,uniq=TRUE)
      match <- getaccC(mass(exp),mass(calib),error=error,uniq=TRUE)
      expd <- c(expd,exp[match$plind,1])
      theo <- c(theo,ppg[match$calind,1])
    }
  mord <- order(theo)
  expd <- expd[mord]
  theo <- theo[mord]
  require(modreg)
                                        #calculate mass dependent error
  error <- (expd-theo)*1e6/theo
  ispl <- smooth.spline(theo,error)
  ispl <- calibspline(ispl,object,error,theo)
  return(ispl)
}



applyextcalib.massvectorlist <- function(object,cS,...)
  {
    ##t External Calbiration
    ##- Corrects the massvectorlist for the measurment error stored calibspline object.
    ##d In case of external calibration some sample spots are only dedicated
    ##d to calibration. Calibration samples which produces equidistant
    ##d peaks, which exact masses are known, can be used to precisely
    ##d estimate the mass dependent error function.
    ##+ object : massvectorlist
    ##+ cS : object of class calibspline
    ##+ ... : further parameters
    ##v massvectorlist : calibrated massvectorlist
    ##sa applyextcalib.massvector, applyextcalib.massvectorlist, getextcalib.massvector, getextcalib.massvectorlist
    ##r Gobom J, Mueller M, Egelhofer V, Theiss D, Lehrach H, Nordhoff E, 2002. A calibration method that simplifies and improves accurate determination of peptide molecular masses by MALDI-TOF MS. Anal Chem. 74(15):3915-23.
    ##r Wolski \url{http://www.molgen.mpg.de/~wolski/mscalib}
    ##e data(ppg)
    ##e data(mvl)
    ##e res <- getextcalib(ppg,getPPGmasses(),error=150)
    ##e mvl2 <- applyextcalib(mvl,res)
    
    if(!inherits(cS,"calibspline"))
      {
        print("second param should be a calibspline")
        return()
      }
    res<-lapply(object,applyextcalib,cS)
    res<- massvectorlist(experiment(object),res,project(object))
    res
  }


calibexternal.massvectorlist <- function(object,ppg,calib,error=300,...)
{
  ##t External Calbiration
  ##- Perfroms external calibration of massvectorlist. Obtains the model of the measurment error by \code{getextcalib} and corrects the masses for this error.
  ##d In case of external calibration some sample spots are only dedicated
  ##d to calibration. Calibration samples which produces equidistant
  ##d peaks, which exact masses are known, can be used to precisely
  ##d estimate the mass dependent error function.
  ##+ object : massvectorlist
  ##+ ppg : either a massvector or massvectorlist with masses of the calibration sample (e.g poly-(propylene glycol) ppg)
  ##+ calib : a massvector with exact masses for the massvectors of the calibration samples (ppg).
  ##+ error : relative error of the measurment in (ppm).
  ##+ ... : further parameters.
  ##v massvectorlist: calibrated massvector
  ##sa calibexternal.massvector, applyextcalib.massvectorlist, getextcalib.massvector, getextcalib.massvectorlist, calibexternal.massvector, applycalib.calibspline
  ##r Gobom J, Mueller M, Egelhofer V, Theiss D, Lehrach H, Nordhoff E, 2002. A calibration method that simplifies and improves accurate determination of peptide molecular masses by MALDI-TOF MS. Anal Chem. 74(15):3915-23.
  ##r Wolski \url{http://www.molgen.mpg.de/~wolski/mscalib}
  ##e data(mvl)
  ##e data(ppg)
  ##e res <- calibexternal(mvl,ppg,getPPGmasses(),error=150)
  cS <- getextcalib(ppg,calib,error=error)
  object <- applyextcalib(object,cS)
  object
}

##########################################################
##Recalibration

recalibrate.massvectorlist <- function(object,PQM=7,...)
  {
    ##t Precalibration
    ##d Precalibration method utilizes the knowledge that masses
    ##d of peptides are in equidistant spaced clusters. The wavelength of
    ##d the massesvector can be determined as described by
    ##d Wool. The comparision of the experimental wavelength with
    ##d the theoretical one, makes possible to find an affine function
    ##d that corrects the masses. Chemical noise in the spectra may hamper
    ##d the determination of mass list frequency. The package provides a
    ##d function to filter chemical noise.
    ##sa recalibrate.massvector
    ##- Obtains the error and performs the calibration in one step.
    ##+ object : massvectorlist.
    ##+ PQM : Peak Quality Measure. Indicates how well the wavelenght of the massvector was determined.
    ##+ ... : further parameters.
    ##v massvectorlist : recalibrated massvectorlist.
    ##r Wool A, Smilansky Z 2002. Precalibration of matrix-assisted laser desorption/ionization-time of flight spectra for peptide mass fingerprinting. Proteomics. 2(10):1365-73.
    ##r Wolski  \url{http://www.molgen.mpg.de/~wolski/mscalib}
    ##e data(mvl)
    ##e mvl<-mvl[1:100]
    ##e pp <- recalibrate(mvl,PQM=8)
    ##e compare(pp[[1]],mvl[[1]],error=1,ppm=FALSE)

    res<-lapply(object,recalibrate,PQM=PQM)
    res<- massvectorlist(experiment(object),res,project(object))
    res
  }


getrecalib.massvectorlist<-function(object,plot=FALSE,...)
  {
    ##t Precalibration
    ##d Precalibration method utilizes the knowledge that masses
    ##d of peptides are in equidistant spaced clusters. The wavelength of
    ##d the massesvector can be determined as described by
    ##d Wool. The comparision of the experimental wavelength with
    ##d the theoretical one, makes possible to find an affine function
    ##d that corrects the masses. Chemical noise in the spectra may hamper
    ##d the determination of mass list frequency. The package provides a
    ##d function to filter chemical noise.
    ##- Obtain the error model using the wavelength analysis of the peaklist.
    ##+ object : massvectorlist
    ##+ plot: logical; \code{TRUE} - A graphic showing the \eqn{F(\omega)}{F(omega)}. default \code{FALSE}
    ##+ ... : further parameters.
    ##v caliblist : caliblist with objects of class calibrestat.
    ##sa getrecalib.massvector,  applyrecalib.massvector, massvector, calibrestat, wsFilter.massvectorlist
    ##r Wool A, Smilansky Z 2002. Precalibration of matrix-assisted laser desorption/ionization-time of flight spectra for peptide mass fingerprinting. Proteomics. 2(10):1365-73.
    ##r Wolski \url{http://www.molgen.mpg.de/~wolski/mscalib}
    ##e data(mvl)
    ##e mvl<-mvl[1:100]
    ##e res <- getrecalib(mvl)
    ##e print(res)
    ##e summary(res)
    ##e image(res,what="Coef.Intercept")
    ##e image(res,what="Coef.Slope")
    ##e plot(res)
    ##e hist(res)
    ##e dres<-as.data.frame(res)
    ##e plot(dres$Coef.Intercept,dres$PQM,xlab="Coef.Intercept",ylab="PQM")
    ##e #create subset.
    ##e res2<-subset(res,PQM>10)
    ##e length(res2)
    ##e plot(res2)
    ##e test<-applyrecalib(mvl, res2)
    ##e plot(test)
    res<-lapply(object,getrecalib,plot=FALSE)
    res<-caliblist("calibrelist",experiment(object),res,project=project(object))
    #res
  }



applyrecalib.massvectorlist <- function(object,calc,...)
  {
    ##t Recalibration
    ##- Corrects the massvector for model of the measurment error stored in the calibrestat object.
    ##d Precalibration method utilizes the knowledge that masses
    ##d of peptides are in equidistant spaced clusters. The wavelength of
    ##d the massesvector can be determined as described by
    ##d Wool. The comparision of the experimental wavelength with
    ##d the theoretical one, makes possible to find an affine function
    ##d that corrects the masses. Chemical noise in the spectra may hamper
    ##d the determination of mass list frequency. The package provides a
    ##d function to filter chemical noise.
    ##+ object : massvectorlist.
    ##+ calc : caliblist with objects of class calibretstat.
    ##+ ... : further parameters.
    ##v massvectorlist : calibrated massvectorlist. 
    ##sa recalibrate.massvectorlist, getrecalib.massvectorlist, correctinternal.massvectorlist
    ##r Wolski \url{http://www.molgen.mpg.de/~wolski/mscalib}
    ##e data(mvl)
    ##e res<-getrecalib(mvl)
    ##e summary(res)
    ##e mvl2<-applyrecalib(mvl,res)


    
    if(!inherits(calc,"calibrelist"))
      stop("Second argument have to be a calibrelist object!!! (use getrecalib)\n")
    for(x in 1:length(calc))
      {
        nami<-names(calc)[x]
        tmp <- object[[nami]]
        tmp<-applyrecalib(tmp,calc[[x]])
        object[[nami]]<-tmp
        #if(x%%10==0)
        #  cat(formatC(x,width=3)," ",sep="")
        #if(x%%100==0)
        #  cat("\n")
      }
    cat("\n")
    object
  }


#########################################
#affine calibration
#

applyintcalib.massvectorlist <- function(object,calc,...)
  {
    ##t Internal Calibration
    ##- Corrects the massvectors in the list for the error model stored in calibintstat object.
    ##d Internal calibration aligns masses of
    ##d peaks to known masses and determines by linear regression a affine
    ##d function that describing the relative error. The internal
    ##d correction fails when no calibration peaks can be found.
    ##+ object : massvectorlist
    ##+ calc : caliblist with objects of class calibintstat
    ##+ ... : further params
    ##v massvectorlist : calibrated massvectorlist. 
    ##sa applyintcalib.massvector,getintcalib.massvectorlist, correctinternal.massvectorlist, calibintstat, caliblist
    ##r Wolski
    ##e data(mvl)
    ##e mvl<-mvl[1:100]
    ##e data(cal)
    ##e res <- getintcalib(mvl,cal,error=300,ppm=FALSE)
    ##e mvl2<-applyintcalib(mvl,res)
    if(!inherits(calc,"calibintlist"))
      stop("Second argument have to be a calibintlist object!!! (use getintcalib)\n")
                                       #calconstants -  output of function recalibrate
    for(x in 1:length(calc))
      {
        tmp <- object[[names(calc)[x]]]
        tmp <- applyintcalib(tmp,calc[[x]])
        object[[names(calc)[x]]] <- tmp
                                        #        if(x%%10==0)
                                        #          cat(formatC(x,width=3)," ",sep="")
                                        #        if(x%%100==0)
                                        #          cat("\n")
      }
    cat("\n")
    object
  }


correctinternal.massvectorlist <- function(object,calib,error=500,ppm=TRUE,...)
  {
    ##t Internal Calibration
    ##- Determines the measurment error of the masses using \code{getintcalib.massvectorlist}.
    ##- It refines the error model and applies it to the massvector by \code{applyintcalib.massvectorlist}.
    ##d \bold{Internal calibration} aligns masses of
    ##d peaks to known masses and determines by linear regression a affine
    ##d function that describing the relative error. The internal
    ##d correction fails when no calibration peaks can be found.
    ##+ object : massvectorlist
    ##+ calib : massvector with calibration masses
    ##+ error : assumed measurment error.
    ##+ uniq : logical;\code{TRUE}- use only mass closest to calibration mass. \code{FALSE}- use all masses closer to the calibration mass then given error.
    ##+ ppm : logical;\code{TRUE}- describe the error as relative error. \code{FALSE}- describe the error as absolute error.
    ##v massvectorlist : calibrated massvectorlist. 
    ##sa getintcalib.massvectorlist, correctinternal.massvector, getintcalib.massvector,calibintstat,caliblist
    ##r wolski
    ##e data(mvl)
    ##e data(cal)
    ##e mvl2 <- correctinternal( mvl,cal, error=500 , ppm=TRUE )
    tmp<-getintcalib(object,calib,error=error,ppm=ppm)
    object<-applyintcalib(object,tmp)
    object
  }

getintcalib.massvectorlist <- function(object,calib,error=500,ppm=TRUE, ... )
  {
    ##t Internal Calbiration
    ##- Obtains error model using massvector with known masses.
    ##d Internal calibration aligns masses of
    ##d peaks to known masses and determines by linear regression a affine
    ##d function that describing the relative error. The internal
    ##d correction fails when no calibration peaks can be found.
    ##+ object : massvectorlist
    ##+ calib : massvector with calibration masses
    ##+ error : assumed measurment error.
    ##+ uniq : \code{TRUE}- use only mass closest to calibration mass. \code{FALSE}- use all masses closer to the calibration mass then given error.
    ##+ ppm : \code{TRUE}- describe the error as relative error. \code{FALSE}- describe the error as absolute error.
    ##v calibintstat : object of class calibintstat. 
    ##sa applyintcalib.massvector, getintcalib.massvector, correctinternal.massvectorlist, calibintstat, caliblist
    ##r Wolski
    ##e data(mvl)
    ##e data(cal)
    ##e res<-getintcalib(mvl,cal,error=400,ppm=TRUE)
    ##e plot(res)
    
    res <- lapply(object,getintcalib,calib,error=error,ppm=ppm)
    res <- caliblist("calibintlist",experiment(object) , res, project=project(object))
    res
  }



##################################################
#Global calibration
#
getglobalcalib.massvectorlist<-function(object , calib , error=500,ppm=TRUE, labund=12, abund = length(object)/5 , accur = ifelse(ppm,error/2000,error) ,...)
  {
    ##t Set Based Internal Calbiration
    ##- Obtains error model using massvector with known masses.
    ##d Set based Calibration copes with the problem of missing calibration
    ##d masses. It first extracts about 15 most abundant masses of the
    ##d massvectorlist, then they are internally calibrated and used
    ##d as new calibration masses. In this fashion more massvectors
    ##d can be internally calibrated.
    ##d Internal calibration aligns masses of
    ##d peaks to known masses and determines by linear regression a affine
    ##d function that describing the relative error. The internal
    ##d correction fails when no calibration peaks can be found.
    ##+ object : massvectorlist
    ##+ calib : massvector with calibration masses
    ##+ error : assumed measurment error.
    ##+ uniq : \code{TRUE}- use only mass closest to calibration mass. \code{FALSE}- use all masses closer to the calibration mass then given error.
    ##+ ppm : logical; \code{TRUE}- describe the error as relative error. \code{FALSE}- describe the error as absolute error.
    ##+ labund : how many abundant masses to use for calibration. default=12.
    ##+ abund : default =  length(object)/5
    ##+ accur : default =  ifelse(ppm,error,error/2000)
    ##+ ... : further parameters.
    ##v calibintstat : object of class calibintstat. 
    ##sa applyintcalib.massvector, getintcalib.massvector, correctinternal.massvectorlist, calibintstat, caliblist
    ##r Wolski
    ##e data(mvl)
    ##e data(cal)
    ##e res<-getglobalcalib(mvl,cal,error=500,ppm=TRUE)
    ##e hist(res)
    if(accur<0.01) stop("error to small")
    mabund <- gamasses(object , accur=accur , abund=abund , ...)
    while(length(mabund)<15)
      {
        print("finding abundant masses")
        abund <- abund - 5
        print(abund)
        mabund <- gamasses(object , accur=accur , abund=abund , ... )
      }
    if(length(mabund)>labund)
      {
        ord<-order(mabund[,2])
        mabund <- mabund[ord,]
        mabund<- mabund[c((length(mabund)-(labund-1)):length(mabund)),]
        ord<-order(mabund[,1])
        mabund <- mabund[ord,]
      }
    abund <- correctinternal(mabund,calib,error=error,ppm=ppm)
    res <- getintcalib(object,abund,error=error,ppm=ppm)
    res
  }



globalcalib.massvectorlist<-function(object , calib , error=500, labund = 12, ppm=TRUE , abund=length(object)/5,accur = ifelse(ppm,error/2000,error) ,...)
  {
    ##t Set Based Internal Calibration
    ##- Determines the error and corrects for it.
    ##d Set based Calibration copes with the problem of missing calibration
    ##d masses. It first extracts about 15 most abundant masses of the
    ##d massvectorlist, then they are internally calibrated and used
    ##d as new calibration masses. In this fashion more massvectors
    ##d can be internally calibrated.
    ##d Internal calibration aligns masses of
    ##d peaks to known masses and determines by linear regression a affine
    ##d function that describing the relative error. The internal
    ##d correction fails when no calibration peaks can be found.
    ##+ object : massvectorlist
    ##+ calib : massvector with calibration masses
    ##+ error : assumed measurment error.
    ##+ labund : how many abundant masses use for calibration. (8-12 masses are sufficient).
    ##+ accur : default := ifelse(ppm,error/2000,error), used to determine abundant masses.
    ##+ abund : default  := length(object)/5
    ##+ ppm : \code{TRUE}- describe the error as relative error. \code{FALSE}- describe the error as absolute error.
    ##+ ... : further parameters.
    ##v massvectorlist : calibrated massvectorlist. 
    ##sa getintcalib.massvectorlist, correctinternal.massvectorlist, getintcalib.massvector,calibintstat,caliblist
    ##r wolski
    ##e data(mvl)
    ##e data(cal)
    ##e mvl2<-globalcalib(mvl,cal,error=500,labund=12)

    mabund <- gamasses(object , accur=accur , abund=abund , ...)
    print(mabund)
    while(length(mabund)<10)
      {
        abund <- abund - 10
        mabund <- gamasses(object , accur=accur , abund=abund , ... )
      }
    if(length(mabund)>labund)
      {
        ord<-order(mabund[,2])
        mabund <- mabund[ord,]
        mabund<- mabund[c((length(mabund)- (labund-1)):length(mabund)),]
      }
    print(mabund)
    abund <- correctinternal(mabund,calib,error=error)
    object <- correctinternal(object,mabund,error=error)
    object
  }

#############################################
#getting abundant masses
#

gamasses.massvectorlist <- function(object,accur=0.1,abund=50,...)
  {
    ##t Abundant masses
    ##- Determines abundant masses in a massvector.
    ##d Abundant masses are masses that occur in a large
    ##d fraction of the massvectors. Typical abundant mass
    ##d are derived from tryptic autoproteolysis products.
    ##d Abundant masses can also often be assigned to keratin
    ##d isoforms (human hair- and skin proteins).
    ##d Many of the abundant masses cannot be assigned to any protein.
    ##d Abundant masses can be used for calibration.
    ##d Removing them may increase the identification specificity.
    ##sa gamasses.massvector
    ##+ object : massvectorlist
    ##+ accur : measurment accuracy
    ##+ abund : how many times a mass have to occur to be an abundant mass.
    ##v massvector : massvector with abundant masses.
    ##r Wolski
    ##e data(mvl)
    ##e #Filtering for abundant masses.
    ##e res<-gamasses(mvl,abund=50)
    ##e plot(res)
    ##e mvFilter(mvl[[1]],res)
    ##e res2<-mvFilter(mvl,res,abundant=TRUE)
    ##e image(res2,what="lengthmv")
    ##e image(mvl,what="lengthmv")
    ##e image(image(res2,what="lengthmv")/image(mvl,what="lengthmv"))
    ##e 
    ##e hist(mvl,accur=0.3)
    ##e hist(res2,add=TRUE,col=2,accur=0.3)
    tmp<-unlist(object)
    res <- gamasses(tmp,accur=accur,abund=abund,main=main,xlab=xlab,xlim=xlim,...)
    res
  }




diffFilter.massvectorlist <- function(object,listofdiffs,higher=TRUE,error=0.05,uniq=TRUE,prune=FALSE,...)
{
  ##t Abundant Differences
  ##- Removes masses from the massvector.
  ##d Removes one of the masses contributing to a mass difference given in the list of diffs.
  ##d Can be used if a variable modification are present in the massvector but can not be considered by the identification software.
  ##d Abundant intra massvector mass differences indicate the
  ##d presence of variable modifications in the data set. This
  ##d information can be used to optimize the search strategy.
  ##sa diffFilter.massvector,getdiff.massvector, getdiff.massvectorlist
  ##+ mv : massvector
  ##+ listoffdiffs : massvector with mass differences
  ##+ higher : \code{TRUE} - remove higher mass, \code{FALSE} = remove lower mass.
  ##+ error : How much the differences can diviate from the differences given in listofdiffs
  ##+ prune : logical;default = \code{TRUE} - removes modified mass; \code{FALSE} - return modified mass.
  ##v massvector : filtered massvector.
  ##r Wolski
  ##e data(mvl)
  ##e res<-getdiff(mvl,range=c(0,100))
  ##e hist(res)
  ##e res<-gamasses(res,abund=100)
  ##e test<-diffFilter(mvl,res,higher=TRUE,error=0.1,uniq=TRUE)
  ##e test
  res <- lapply(object,diffFilter, listofdiffs, higher=higher, error=error,uniq=uniq,prune=prune)
  res <- massvectorlist(experiment(object),res,project=project(object))
  res
}

getdiff.massvectorlist<-function(object,masst="massrea",range=c(0,100),...)
{
  ##t Massdifferences
  ##- Computes mass differences in the mv.
  ##d Removes one of the masses contributing to a mass difference given in the list of diffs.
  ##d Can be used if a variable modification are present in the massvector but can not be considered by the identification software.
  ##d Abundant intra massvector mass differences indicate the
  ##d presence of variable modifications in the data set. This
  ##d information can be used to optimize the search strategy.
  ##sa getdiff.massvector
  ##+ mv : massvector
  ##+ listoffdiffs : massvector with mass differences
  ##+ higher : \code{TRUE} - remove higher mass, \code{FALSE} - remove lower mass.
  ##+ error : How much the differences can diviate from the differences given in listofdiffs
  ##v massvector : filtered massvector.
  ##r Wolski
  ##e data(mvl)
  ##e res<-getdiff(mvl,range=c(0,100))
  ##e plot(res)
  ##e tt<-gamasses(res,abund=40)
  ##e plot(tt)
  res<-massvector(NULL)
  res<-lapply(object,getdiff,range=c(0,100))
  res<-massvectorlist(experiment(object),res,project=project(object))
  res<-unlist(res)
  res<-info(res,paste(info(object),"diffs",sep="_"))
  return(res)
}


writeF.massvectorlist<-function(object,path,file=experiment(object),ext="txt",...)
  {
    ##t Write massvectorlist
    ##- Write massvectorlist to File
    ##d The read and write functions for all the different peak-list formats are not provided by the package. This is because
    ##d there are oodles of different formats. I will try to collect read-write functions for as many as possible peak-list format's in an add on package
    ##d which you can find at \url{http://www.molgen.mpg.de/~wolski/mscalib/IO/}.
    ##+ object : massvectorlist.
    ##+ path : path to directory.
    ##+ file : file name. default experiment(object)
    ##+ ext : file extension. default txt.
    ##sa readF.massvectorlist, readF.massvector

    
    filep <- file.path(path,paste(file,".",ext,sep=""),fsep = .Platform$file.sep)
    con <- file(filep, "w")  # open an output file connection

    intern <-function(lobject,con)
      {
        firstline<-paste(">",info(lobject),":",join(mget(lobject,"tcoor"),sep=","),sep="")
        writeLines(firstline,con=con)
        write.table(lobject,file = con, append = TRUE, quote = FALSE, sep = "\t",
                    eol = "\n", na = "NA", dec = ".", row.names = FALSE,
                    col.names = FALSE, qmethod = c("escape", "double"))
      }
    res<-lapply(object,intern,con)
    close(con)
  }

readF.massvectorlist<-function(object,path,file=experiment(object),ext="txt",...)
  {
    ##t Reads massvectorlist
    ##- written to disk with \code{writeF.massvectorlist}
    ##d The read and write functions for all the different peak-list formats are not provided by the package. This is because
    ##d there are oodles of different formats. I will try to collect read-write functions for as many as possible peak-list format's in an add on package
    ##d which you can find at \url{http://www.molgen.mpg.de/~wolski/mscalib/IO/}.
    ##+ object :massvectorlist
    ##+ path : path to directory.
    ##+ file : file to read. default =experiment(object)
    ##sa readF.massvector,writeF.massvectorlist
    ##e data( mvl )
    ##e mvl
    ##e writeF( mvl , "." )
    ##e test <- readF(massvectorlist(info(mvl)),".")
    ##e test
    ##e unlink(paste(info(test),".txt",sep=""))
    
    filep <- file.path(path,paste(file,".",ext,sep="") , fsep = .Platform$file.sep)
    con<-file(filep,"r")
    mm <- readLines(con=filep,n=-1)
    starts <- grep(">",mm)
    starts <- c(starts,length(mm))
    res<-massvectorlist(file)
    setParms(res)<-list(project=project(object))
    intern<-function(xx)
      {
        object<-massvector()
        res1 <- unlist(strsplit(xx[1],":"))
        object<-info(object,sub(">","",res1[1]))
        setParms(object) <- list(tcoor= unlist(strsplit(res1[2],",")))
        intern2<-function(x){as.numeric(unlist(strsplit(x,"\t")))}
        tt<-xx[2:length(xx)]
        rxx <- t(sapply(tt,intern2))
        rxx<-rxx[order(rxx[,1]),]
        rownames(rxx)<-1:length(rxx[,1])
        object <- peaks(object,rxx)
        return(object)
      }

    for(x in 1:(length(starts)-1))
      {
        res[[x]]<-intern(mm[starts[x]:(starts[x+1]-1)])
      }
    close(con)
    return(res)
  }
#Copyright 2004, W. Wolski, all rights reserved.
hist.mylistobj<-function(x,...)
  {
    ##t hist
    ##- hist
    ##d mylistobj provides basic functionality for objects implemented using list to store the attributes.
    ##+ x : object of class mylistobj
    warning("hist not implemented for object of class mylistobj!\n")
  }

image.mylistobj <- function(x,...)
  {
    ##t image
    ##- image
    ##d mylistobj provides basic functionality for objects implemented using list to store the attributes.
    ##+ x : object of class mylistobj
    warning("image not implemented for object of calss mylistobj!\n")
  }


mget.mylistobj <- function(object,attrn,...)
  {
    ##t Field Access
    ##- Access to the fields in the mylistobj
    ##d mylistobj provides basic functionality for objects implemented using list to store the attributes.
    ##+ object : mylistobj
    ##+ attrn : The of the field to access. If missing the fields of the objects are returned.
    ##+ ... : further parameters
    ##v xxx : depends which field in the massvector are accessed.
    ##e data(mv1)
    ##e res <- getrecalib(mv1)
    ##e class(res)
    ##e mget(mv1,"lenghtmv")
    ##e mget(mv1,"Coef.Intercept")
    
    allow <- object$allow
    if(attrn %in% allow)
      {
          return(object[[attrn]])
      }
    else
      {
        stop( attrn , " not in the attributes list :", join(allow,sep=" ") ,"\n")
      }
  }


info.mylistobj <- function(object,info , ...)
  {
    ##t Info Acces
    ##- access to the info field of the mylistobj
    ##d mylistobj provides basic functionality for objects implemented using list to store the attributes.
    ##+ object : mylistobj
    ##+ info : info character. If missing function returns the current info. If not missing function returns mylistobj with new info field content.
    ##e data(mv1)
    ##e res <- getrecalib(mv1)
    ##e info(res)
    ##e res<-info(res,"testname")

    if(missing(info))
      mget(object,"info")
    else
      setParms(object)<-list(info=info)
  }

plot.mylistobj<-function(x,...)
  {
    ##t plot
    ##- plot
    ##d mylistobj provides basic functionality for objects implemented using list to store the attributes.
    ##+ x : object of class mylistobj
    stop(" plot not implemented for object of class mylistobj!\n")
  }


"setParms<-.mylistobj"<-function(object,value)
  {
    ##t Field Access
    ##- Set attributes in object of class mylistobj
    ##d mylistobj provides basic functionality for objects implemented using list to store the attributes.
    ##+ object : object of class myobj
    ##+ value : a list where list names are attributes names.
    ##sa "setParms<-.myobj",mget.mylistobj
    ##e data(mv1)
    ##e res<-getrecalib(mv1)
    ##e setParms(res) <- list(info="test")
    ##e res

    tmp<-value
    allow<-object$allow
    for(x in names(tmp))
      {
        if(x %in% allow)
          {
            object[[x]]<-tmp[[x]]
          }
        else{
          stop(x," are not a allowed attribute list:",join(object$allow,sep=" "),"\n")
        }
      }
    object
  }
#Copyright 2004, W. Wolski, all rights reserved.
mget.myobj <- function(object,attrn,...)
  {
    ##t Field Access
    ##- Acces fields in object of class myobj
    ##+ object : object of class myobj
    ##+ attrn : name of field (Attribute)
    ##e data(mv1)
    ##e mget(mv1)
    ##e mget(mv1,"info")

    allow <- attr(object,"allow")
    if(attrn %in% allow)
      {
          return(attr(object,attrn))
      }
    else
      {
        warning( attrn , " not in the attributes list :", join(allow,sep=" ") ,"\n")
      }
  }

"setParms<-.myobj" <- function(object,value)
  {
    ##t Field Access
    ##- Set attributes in object of class myobj
    ##+ object : object of class myobj
    ##+ value : a list where list names are attributes names.
    ##e data(mv1)
    ##e mv1
    ##e setParms(mv1) <- list(info="test")
    ##e mv1
    
    tmp<-value
    allow<-attr(object,"allow")
    for(u in names(tmp))
      {
        if(u %in% allow)
          {
            if(!(u %in% "mass"))
              {
                attr(object,u)<-tmp[[u]]
              }
            else
              {
                tp<-attributes(object)
                object<-sort(tmp[[u]])
                attributes(object)<-tp
                return(object)
              }
          }
        else
          {
            warning(u," are not a allowed attribute list:", print(object$allow) ,"\n")
          }
      }
    object
  }


info.myobj <- function(object,info,...)
  {
    ##t Info Acces
    ##- access to the info field of the massvector.
    ##+ object : massvector
    ##+ info : info character. If missing function returns the current info. If not missing function returns massvector with new info field content.
    ##e data(mv1)
    ##e info(mv1)
    ##e mv1<-info(mv1,"testname")
    if(!missing(info))
      {
        attr(object,"info")<-info
        return(object)
      }
    return(attr(object,"info"))
  }
#Copyright 2004, W. Wolski, all rights reserved.
##calibrestat

##calibrelist

calibrestat<-function(info,...)
  {
    ##t Constructor
    ##- Returns object of class calibrestat. Used by function \code{getrecalib.massvector}.
    ##- A calibrestat object is returned by the \code{getrecalib.massvector} method.
    ##+ info : unique identifier.
    ##+ ... : can be Coeff.Intercept,Coeff.Slope,lengthvm,PQM,tcoor.
    ##sa getrecalib.massvector
    ##e data(mv1)
    ##e res <- getrecalib(mv1)
    ##e print(res)
    ##e as.vector(res)
    ##e summary(res)
    ##e image(res)
    ##e plot(res)

    
    res<-list()
    tmp<-list(...)
    allow<-c("info","Coeff.Intercept","Coeff.Slope","lengthmv","PQM","tcoor")
    res$allow<-allow
    if(!missing(info))
      res$info<-info
    for(x in names(tmp))
      {
        if(x %in% allow)
          {
            res[[x]]<-tmp[[x]]
          }
      }
    class(res)<-c("calibrestat","calibstat","mylistobj")
    res
  }


print.calibrestat<-function(x,...)
  {
    ##t Print Calibrestat Object
    ##- Prints the fields of the calibrestat object.
    ##+ x : object of class calibrestat
    ##sa print
    ##e data(mv1)
    ##e res<-getrecalib(mv1)
    ##e print(res)
    
    tmp<-x$allow
    cat("class           :", class(x),"\n")
    cat("info            :",mget(x,"info"),"\n")
    cat("lengthmv        :",mget(x,"lengthmv"),"\n")
    cat("Coeff.Intercept :",mget(x,"Coeff.Intercept"),"\n")
    cat("Coeff.Slope     :",mget(x,"Coeff.Slope"),"\n")
    cat("PQM             :",mget(x,"PQM"),"\n")
    cat("tcoor           :",mget(x,"tcoor"),"\n")
    invisible(x)
  }


as.vector.calibrestat<-function(x,...)
  {
    ##t Coerces to vector
    ##- coerces the calibrestat object into a vector.
    ##+ x : calibrestat
    ##e data(mv1)
    ##e res<-getrecalib(mv1)
    ##e as.vector(res)
    res<-c(
           ifelse(is.null(x$lengthmv),NA,x$lengthmv),
           ifelse(is.null(x$Coeff.Intercept),NA,x$Coeff.Intercept),
           ifelse(is.null(x$Coeff.Slope),NA,x$Coeff.Slope),
           ifelse(is.null(x$PQM),NA,x$PQM),
           if(is.null(x$tcoor)){c(NA,NA)}else{x$tcoor}
           )
    names(res)<-c("lengthmv","Coef.Intercept","Coef.Slope","PQM","Xcoor","Ycoor")
    res
  }



applycalib.calibrestat<-function(object,mv,...)
  {
    ##t Precalibration
    ##d \bold{Precalibration} method utilizes the knowledge that masses
    ##d of peptides are in equidistant spaced clusters. The wavelength of
    ##d the \emph{massesvector} can be determined as described by
    ##d Wool. The comparision of the experimental wavelength with
    ##d the theoretical one, makes possible to find an affine function
    ##d that corrects the masses. Chemical noise in the spectra may hamper
    ##d the determination of mass list frequency. The package provides a
    ##d function to filter chemical noise.
    ##- Uses the error model obtained by the method \code{getrecalib} to correct masses in the massvector.
    ##+ object : calibrestat
    ##+ mv : massvector
    ##+ ...: further arguments
    ##w massvector : recalibrated massvector.
    ##sa applycalib.calibintstat, getrecalib.massvector, getrecalib.massvectorlist, recalibrate.massvector, recalibrate.massvectorlist, calibrestat, calibrelist
    ##r Wool A, Smilansky Z 2002. Precalibration of matrix-assisted laser desorption/ionization-time of flight spectra for peptide mass fingerprinting. \bold{Proteomics.} 2(10):1365-73.
    ##r Wolski \url{http://www.molgen.mpg.de/~wolski/mscalib}
    ##e data(mv1)
    ##e res <- getrecalib(mv1)
    ##e plot(res)
    ##e mv2<-applycalib(res,mv1)
    ##e plot(mv1[,1],mv2[,1]-mv1[,1])
    peak<- (mv[,1]-mget(object,"Coeff.Intercept"))*(mget(object,"Coeff.Slope")/1e6+1)
    mv[,1] <- peak
    return(mv)
  }

###############################################
## calibrelist
#


applycalib.calibrelist<-function(object,mvl,...)
{
  ##t Precalibration
  ##- Uses the error model obtained by the method \code{getrecalib} and stored in a \code{calibrelist}
  ##- to correct masses in the massvectorlist.
  ##d \bold{Precalibration} method utilizes the knowledge that masses
  ##d of peptides are in equidistant spaced clusters. The wavelength of
  ##d the \emph{massesvector} can be determined as described by
  ##d Wool. The comparision of the experimental wavelength with
  ##d the theoretical one, makes possible to find an affine function
  ##d that corrects the masses. Chemical noise in the spectra may hamper
  ##d the determination of mass list frequency. The package provides a
  ##d function to filter chemical noise.
  ##+ object : calibrelist
  ##+ mvl : massvectorlist
  ##+ ... : further parameters
  ##v massvectorlist : calibrated massvectorlist. 
  ##sa recalibrate.massvectorlist, getrecalib.massvectorlist, correctinternal.massvectorlist
  ##r Wolski \url{http://www.molgen.mpg.de/~wolski/mscalib}
  ##e data(mvl)
  ##e mvl<-mvl[1:100]
  ##e res<-getrecalib(mvl)
  ##e plot(res)
  ##e image(res,what="PQM")
  ##e mvlr<-applycalib(res,mvl)
  
   if(!inherits(mvl,"massvectorlist"))
      stop(as.character(substitute(mvl)),"have to be a object of class massvectorlist!\n")
    for(x in 1:length(object))
      {
        nami<-names(object)[x]
        tmp <- mvl[[nami]]
        tmp<-applyrecalib(tmp,object[[x]])
        mvl[[nami]]<-tmp
        if(x%%10==0)
          cat(formatC(x,width=3)," ",sep="")
        if(x%%100==0)
          cat("\n")
      }
    cat("\n")
    mvl
}

plot.calibrelist<-function(x,...)
  {
    ##t Plot
    ##- Shows a matrix of scatterplots for different varialbes..
    ##+ x : calibrelist
    ##+ ... : graphical parameters can be given as arguments to plot.
    ##e data(mvl)
    ##e mvl<-mvl[1:100]
    ##e res2 <- getrecalib(mvl) # get recalibration model for not filtered data
    ##e plot(res2)

    dat<-as.matrix(x)[,1:4]
    colnames(dat)<-names(as.vector(x[[1]]))[1:4]
    plot(data.frame(dat),pch="*",...)
    par(mfrow=c(1,1))
  }
  
hist.calibrelist<-function(x,...)
  {
    ##t Histogram Plot
    ##- Computes Histograms of the lengthmv, Coef.Intercept, Coef.Slope, PQM.
    ##+ x : calibrelist
    ##+ ... : further graphical parameters.
    ##e data(mvl)
    ##e mvl<-mvl[1:100]
    ##e res2 <- getrecalib(mvl) # get recalibration model for not filtered data
    ##e hist(res2)
    
    dat<-as.matrix(x)
    colnames(dat)<-names(as.vector(x[[1]]))
    par(mfrow=c(2,2))
    par(cex.main=0.6)
    hist(dat[,1],main="Recalibration",xlab="lengthmv",...)
    hist(dat[,2],main="Recalibration",xlab="Coef.Intercept",...)
    hist(dat[,3],main="Recalibration",xlab="Coef.Slope",...)
    hist(dat[,4],main="Recalibration",xlab="PQM",...)
    par(cex.main=1)
    par(mfrow=c(1,1))
  }
        


        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
#Copyright 2004, W. Wolski, all rights reserved.
.First.lib <- function(lib, pkg) library.dynam("mscalib",pkg,lib)
.Last.lib <- function(libpath) library.dynam.unload("mscalib", libpath)

