
/**************************************
  *         Dataset Version 1.2        *
  *            26.4.1999               *
  * Left to do:                        *
  *         -- get rid of direct       *
  *            datapoints manipulation *
  *         -- routine to add errors   *
  *         -- routine to resample     *
  *         -- improve operators       *
  **************************************/
 
 
 
 #include <iostream>
 #include <dataset.h>
 #include <fstream>
 #include <math.h>
 #include <functional>
 #include <algorithm>
 #include <datapoint.h>
 #include <functions.h>
 
 Dataset::Dataset()
 {
   numberpoints=0;
         comment="No comment";
 }
 
 Dataset::Dataset(const Dataset &data)
 {
   *this=data;
 }
 
 Dataset::Dataset(const int& np)
 {
   int i;
 
   numberpoints=np;
   for (i=0;i<np;i++) {
     dplist.push_back(Datapoint(0.0,0.0));
   }
         comment="No comment";
 }
 Dataset::Dataset(const int& np,const double& min,const double& max)
 {
   int i;
   double d=(max-min)/(np-1);
   dplist.clear();
 
   numberpoints=np;
   for (i=0;i<np;i++) {
     dplist.push_back(Datapoint(i*d+min,0.0));
   }
   comment="No comment";
 }
 Dataset::Dataset(const int& np,const double* xin,const double* yin)
 {
   int i;
   dplist.clear();
 
   numberpoints=np;
   for (i=0;i<np;i++) {
     dplist.push_back(Datapoint(xin[i],yin[i]));
   }
         comment="No comment";
 }
 Dataset::Dataset(const vector<Datapoint>& dplistin)
 {
   int i;
   int np=dplistin.size();
 
   numberpoints=np;
   for (i=0;i<np;i++) {
     dplist.push_back(dplistin[i]);
   }
         comment="No comment";
 }
 Dataset::Dataset(const vector<double>& dplistin)
 {
   int i;
   int np=dplistin.size();
 
   numberpoints=np;
   for (i=0;i<np;i++) {
     dplist.push_back(Datapoint(double(i),dplistin[i]));
   }
         comment="No comment";
 }
 Dataset::Dataset(const valarray<double>& dplistin)
 {
   int i;
   int np=dplistin.size();
 
   numberpoints=np;
   for (i=0;i<np;i++) {
     dplist.push_back(Datapoint(double(i),dplistin[i]));
   }
         comment="No comment";
 }
 Dataset::Dataset(const valarray<double>& dplistinx,const valarray<double>& dplistiny)
 {
   int i;
   int np=dplistinx.size();
 
   numberpoints=np;
   for (i=0;i<np;i++) {
     dplist.push_back(Datapoint(dplistinx[i],dplistiny[i]));
   }
         comment="No comment";
 }
 Dataset::Dataset(const char* filename)
 {
   int i;
         comment="No comment";
   ifstream in_file(filename);
   dplist.clear();
 
   in_file >> (*this);
 }
 
 Dataset::Dataset(const string& filename)
 {
   int i;
         comment="No comment";
   ifstream in_file(filename.c_str());
   dplist.clear();
 
   in_file >> (*this);  
 }
 
 Dataset::~Dataset()
 {
   dplist.clear();
 }
 
 int Dataset::Init(const int& np)
 {
   int i;
 
   numberpoints=np;
         comment="No comment";
   dplist.clear();
   for (i=0;i<np;i++) {
     dplist.push_back(Datapoint(0.0,0.0));
   }
   return ok;
 }
 int Dataset::SetDatapoint(const int& i,const Datapoint& dpin) {
   if (i<0 || i>numberpoints-1) return out_of_bounds; // Error occured -- out of bounds.
   dplist[i]=dpin;
   return ok;
 }
 
 void Dataset::Cumulative()
 {
   int i;
   Datapoint sum(0.0,0.0);
   for (i=numberpoints-1;i>=0;i--) {
     sum.add_y(dplist[i]);
     dplist[i].SetY(sum.GetY());
     dplist[i].SetYErr(sum.GetYErr());
   }
 }
 void Dataset::InvCumulative()
 {
   int i;
   Datapoint sum(0.0,0.0);
   for (i=0;i<numberpoints;i++) {
     sum.add_y(dplist[i]);
     dplist[i].SetY(sum.GetY());
     dplist[i].SetYErr(sum.GetYErr());
   }
 }
 int Dataset::AddTo(const double& xval,const double& addval) {
   return (AddTo(GetIndex(xval),addval));
 }
 int Dataset::AddTo(const int& i,const double& addval) {
   if (i<0 || i>numberpoints-1) return out_of_bounds;
   dplist[i].add_y(addval);
   return ok;
 }
 int Dataset::SetTo(const double& xval,const double& addval) {
   return (SetTo(GetIndex(xval),addval));
 }
 int Dataset::SetTo(const int& i,const double& addval) {
   if (i<0 || i>numberpoints-1) return out_of_bounds;
   dplist[i].SetY(addval);
   dplist[i].SetYErr(0.0);
   return ok;
 }
 
 Datapoint Dataset::Integrate(const double& a,const double& b) const
 {
   Datapoint result(b,0.0);
   int i;
 
   for (i=1;i<numberpoints;i++) {
     if (dplist[i].GetX()<=b && dplist[i].GetX()>=a) result.add_y(dplist[i]*(dplist[i].GetX()-dplist[i-1].GetX())); 
   }
   return result;
 }
 Datapoint Dataset::Sum(const double& a,const double& b) const
 {
   Datapoint result(0.0,0.0);
   int i;
   for (i=0;i<numberpoints;i++) {
     if (dplist[i].GetX()<=b && dplist[i].GetX()>=a) result.add_y(dplist[i]); 
   }
   return result;
 }
 
 int Dataset::ReBin(const int& newbin) {
         if (newbin>=numberpoints) return 1;
         Dataset res(newbin,GetX(0),GetX(numberpoints-1));
         Dataset nums(newbin,GetX(0),GetX(numberpoints-1));
   for (int i=0;i<numberpoints;i++) {
     res.AddTo(dplist[i].GetX(),dplist[i].GetY()); 
     nums.AddTo(dplist[i].GetX(),1.0); 
   }
         res=res/nums;
         *this=res;
 }
 int Dataset::Smooth(int numin)
 {
   int i;
   int ii;
   for (ii=0;ii<numin;ii++) {
     for (i=0;i<numberpoints;i++) {
       dplist[i].SetY(0.33333*(ringacc(i-1).GetY()+ringacc(i).GetY()+ringacc(i+1).GetY()));
       dplist[i].SetYErr(0.33333*(ringacc(i-1).GetYErr()+ringacc(i).GetYErr()+ringacc(i+1).GetYErr()));
     }
   }
   return ok;
 }
 int Dataset::IncreaseDensity(int numin)
 {
   int i;
   int ii;
   vector<Datapoint> insertlist;
   for (i=0;i<numberpoints-1;i++) {
     for (ii=0;ii<numin;ii++) {
       Datapoint apoint(dplist[i].GetX()+(double(ii)/double(numin))*(dplist[i+1].GetX()-dplist[i].GetX()),dplist[i].GetY()+(double(ii)/double(numin))*(dplist[i+1].GetY()-dplist[i].GetY())); // New point calculated using linear interpolation.*/
       insertlist.push_back(apoint);
     }
     insertlist.push_back(dplist[i+1]);
   }
   (*this)=Dataset(insertlist);
   return ok;
 }
 int Dataset::Close()
 {
   if (numberpoints<=0) return empty;
   dplist.push_back(dplist[0]);
   numberpoints++;
   return ok;
 }
 
 int Dataset::RemoveDatapoint(const int& no) {
   dplist.erase(vector<Datapoint>::iterator(&dplist[no]));
   numberpoints--;
   return ok;
 }
 
 int Dataset::TakeLogX()
 {
   int i,flag=0;
   for (i=0;i<numberpoints;i++) {
     if (dplist[i].GetX()<0) {
       flag=flag | log_neg;
       RemoveDatapoint(i);
     }
     if (dplist[i].GetX()==0) {
       flag=flag | log_zero;
       RemoveDatapoint(i);
     }
    if (dplist[i].GetX()>0) dplist[i].TakeLogX();
   }
   return flag;
 }
 int Dataset::TakeLogY()
 {
   int i,flag=0;
   vector<Datapoint>::iterator datpoint;
   for (datpoint=dplist.begin();datpoint<dplist.end();datpoint++) {
     if (datpoint->GetY()<0) {
       flag=flag | log_neg;
       //dplist.erase(datpoint,datpoint+1);
       datpoint->SetY(-1e99);
     }
     if (datpoint->GetY()==0) {
       flag=flag | log_zero;
       datpoint->SetY(-1e99);
       //dplist.erase(datpoint,datpoint+1);
     }
     if (datpoint->GetY()>0) datpoint->TakeLogY();
   }
   return flag;
 }
 Dataset Dataset::Subset(const double& xmin,const double& xmax) const
 {
   int i,count=0;
   for (i=0;i<numberpoints;i++) {
     if (dplist[i].GetX()<xmax && dplist[i].GetX()>xmin) count++;
   }
   Dataset result(count,xmin,xmax);
   count=0;
   for (i=0;i<numberpoints;i++) {
     if (dplist[i].GetX()<xmax && dplist[i].GetX()>xmin) {
       result.dplist[count]=dplist[i];
       count++;
     }
   }
   return result;
 }
 int Dataset::Save(const string& filename) const
 { 
   return (Save(filename.c_str()));
 }
 int Dataset::Save(const char* filename) const
 {  
   int i;
   ofstream out_file(filename);
   out_file << (*this);
   return ok; // Still need to add proper error management.
 }
 int Dataset::Load(const string& filename)
 {  
   return (Load(filename.c_str()));
 }
 int Dataset::Load(const char* filename)
 {  
   int i,np;
   ifstream in_file(filename);
   dplist.clear();
   numberpoints=0;
   in_file >> (*this);
   return ok;
 }
 
 
 void Dataset::Multiply(const double& factor)
 {
   int np,i;
 
   np=numberpoints;
   for (i=0;i<np;i++) {
     dplist[i].mul_y(factor);
   }
   return;
 }
 int Dataset::Normalize()
 {
   int np,i;
   Datapoint norm;
 
   np=numberpoints;
   norm=Integrate(dplist[0].GetX(),dplist[np-1].GetX());
   if (norm.GetY()==0) return div_zero; 
   for (i=0;i<np;i++) {
     dplist[i].div_y(norm);
   }
   return ok;
 }
 double Dataset::Sum() const
 {
   int np,i;
   double result=0.0;
 
   np=numberpoints;
   for (i=0;i<np;i++) {
     result+=dplist[i].GetY();
   }
   return result;
 }
 void Dataset::TextDump()
 {
   int i;
   cout << "Datapoints: "<< numberpoints << "\n";
   for (i=0;i<numberpoints;i++) {
     cout << "Data Element " << i << " has values (x,xerr,y,yerr)= (" << dplist[i] << ").\n";
   }
 }
 
 int Dataset::InitFunction(double (*function)(double x))
 {
   int i;
   
   for (i=0;i<numberpoints;i++) {
     dplist[i].SetY((*function)(dplist[i].GetX()));
   }
   return ok;
 }
 
 double Dataset::Xi(double (*function)(double x),const double& min,const double& max) const
 {
   int i;
   double z,xi;
 
   xi=0;z=0;
   for (i=0;i<numberpoints;i++) {
     if (dplist[i].GetX()<max && dplist[i].GetX()>min) {
       z+=pow(dplist[i].GetY()-(*function)(dplist[i].GetX()),2); // should also divide by sigma here.
     }
   }
   xi=sqrt(z);
   return xi;
 }
 double Dataset::GetValue(const double& xval) const
 {
   double result;
   int i,i1=0,i2=numberpoints-1;
   if (xval>dplist[i2].GetX()) return dplist[i2].GetY();
   while (abs(i1-i2)>1) {
     i=int((i1+i2)/2);
     if (dplist[i].GetX()>=xval) i2=i;
     else i1=i;
     //cout << "xva=" << xval << "datx=" << dplist[i].GetX() << "i1=" << i1 << "i2=" << i2 << "\n";
   }
   double dx=dplist[i2].GetX()-dplist[i1].GetX();
   if (dx==0.0) cerr << "Errrrrrror in dataset getvalue!!\n";
   result=(xval-dplist[i1].GetX())/dx*(dplist[i2].GetY()-dplist[i1].GetY())+dplist[i1].GetY();
   return result;
 }
 int Dataset::GetIndex(const double& xval) const
 {
   int result;
   int i,i1=0,i2=numberpoints-1;
   if (xval>dplist[i2].GetX()) return i2;
   while (abs(i1-i2)>1) {
     i=int((i1+i2)/2);
     if (dplist[i].GetX()>=xval) i2=i;
     else i1=i;
     //cout << "xva=" << xval << "datx=" << dplist[i].GetX() << "i1=" << i1 << "i2=" << i2 << "\n";
   }
   double dx=dplist[i2].GetX()-dplist[i1].GetX();
   if (dx==0.0) cerr << "Errrrrrror in dataset getvalue!!\n";
   if (fabs(xval-dplist[i1].GetX())<fabs(xval-dplist[i2].GetX())) result=i1;
   else result=i2;
   return result;
 }
 Datapoint Dataset::GetMax() const
 {
   Datapoint result=*(max_element(dplist.begin(),dplist.end(),mem_fun_ref(&Datapoint::less_y)));
   return result;
 }
 Datapoint Dataset::GetMin() const
 {
   Datapoint result=*(min_element(dplist.begin(),dplist.end(),mem_fun_ref(&Datapoint::less_y)));
   return result;
 }
 
 Datapoint Dataset::GetMaxX() const
 {
         Datapoint result(*max_element(dplist.begin(),dplist.end(),mem_fun_ref(&Datapoint::less_x)));
   return result;
 }
 Datapoint Dataset::GetMinX() const
 {
   Datapoint result=*(min_element(dplist.begin(),dplist.end(),mem_fun_ref(&Datapoint::less_x)));
   return result;
 }
 
 
 int Dataset::BinValue(const double& value,double v)
 {
   int bin;
   
   if (value<dplist[0].GetX() || value>=dplist[numberpoints-1].GetX()) return out_of_bounds;
   //the above line gave trouble with histogram of cross section plottings -- high end cut off.
   bin=0;
   while (dplist[bin].GetX()<=value && bin<numberpoints) bin++;
   if (bin) dplist[bin-1].add_y(v);
   return ok;
 }
 int Dataset::NormaliseToValue(const double& xval,const double& value)
 {
   int i;
   double ratio;
   double nval;
 
   if (xval<dplist[0].GetX() || xval>dplist[numberpoints-1].GetX()) return out_of_bounds;
   i=0;
   while (dplist[i].GetX()<=xval) i++;
   if (i) nval=dplist[i-1].GetY();
   else nval=dplist[i].GetY();
   if (nval==0) return -1;
         ratio=value/nval;
         Multiply(ratio);
   return ok;
 }
 int Dataset::ShiftToValue(const double& xval,const double& value)
 {
   int i;
   double ratio;
   double nval;
 
   if (xval<dplist[0].GetX() || xval>dplist[numberpoints-1].GetX()) return out_of_bounds;
   i=0;
   while (dplist[i].GetX()<=xval) i++;
   if (i) nval=dplist[i-1].GetY();
   else nval=dplist[i].GetY();
   double dd=value-nval;
   for (i=0;i<dplist.size();i++) {
     dplist[i].add_y(value-nval);
   }
   return ok;
 }
 
 
 void Dataset::ZeroAll()
 {
   int i;
   for (i=0;i<numberpoints;i++) dplist[i]=0.0;
 }
 int Dataset::fill(const double& what,const double& with)
 {
   int i;
   int count=0;
   for (i=0;i<numberpoints;i++) {
     if (dplist[i].GetY()==what) {
       dplist[i].SetY(with);
       dplist[i].SetYErr(with);
       count++;
     }
   }
   return count;
 }
 
 vector<Dataset> Dataset::Split(double t) const
 {  
   double* xdata;
   double* ydata;
   int i,count=0;
   bool closed=false;
   bool ok=false;
   int ls=0;
   double dist=0.0;
   double ldist=0.0;
   vector<Dataset> result;
   vector<Datapoint> dpl;
   xdata=new double[numberpoints];
   ydata=new double[numberpoints];
   double thresh=0.0;
   for (i=0;i<numberpoints-1;i++) {
     xdata[0]=dplist[i].GetX();
     ydata[0]=dplist[i].GetY();
     xdata[1]=dplist[i+1].GetX();
     ydata[1]=dplist[i+1].GetY();
     thresh+=sqrt(pow(xdata[1]-xdata[0],2)+pow(ydata[1]-ydata[0],2));
   }
   thresh=thresh/numberpoints*t; 
   while (count<numberpoints) {
     i=0;closed=false;ok=true;
     dist=0.0;ldist=0.0;dpl.clear();
     while (!closed && i<numberpoints && count<numberpoints) {
       xdata[i]=dplist[count].GetX();
       ydata[i]=dplist[count].GetY();
       if (i>0) dist=sqrt(pow(xdata[i]-xdata[i-1],2)+pow(ydata[i]-ydata[i-1],2));
       if (dist>thresh && ok) {closed=true;i--;count--;}
       else dpl.push_back(dplist[count]);
       i++;count++;
       if (dist>thresh*2.0) ok=true;
     }
     Dataset tempdat(dpl);
     if (i>2 && numberpoints>1 && dpl.size()>0) {
       xdata[0]=dpl[0].GetX();
       ydata[0]=dpl[0].GetY();
       xdata[1]=dpl[dpl.size()-1].GetX();
       ydata[1]=dpl[dpl.size()-1].GetY();
       dist=sqrt(pow(xdata[1]-xdata[0],2)+pow(ydata[1]-ydata[0],2));
       if (dist<thresh) tempdat.Close();
     }
     if (i>2) result.push_back(tempdat);
   }
   delete[] xdata;xdata=NULL;
   delete[] ydata;ydata=NULL;
   return result;
 }
 
 int Dataset::Split(vector<Dataset>& result,double t) const
 {  
   if (numberpoints<2) return empty;
   double* xdata;
   double* ydata;
   int i,count=0;
   bool closed=false;
   int ls=0;
   double dist=0.0;
   double ldist=0.0;
   vector<Datapoint> dpl;
   xdata=new double[numberpoints];
   ydata=new double[numberpoints];
   double thresh=0.0;
   for (i=0;i<numberpoints-1;i++) {
     xdata[0]=dplist[i].GetX();
     ydata[0]=dplist[i].GetY();
     xdata[1]=dplist[i+1].GetX();
     ydata[1]=dplist[i+1].GetY();
     thresh+=sqrt(pow(xdata[1]-xdata[0],2)+pow(ydata[1]-ydata[0],2));
   }
   thresh=thresh/numberpoints*t; 
   while (count<numberpoints) {
     i=0;closed=false;
     dist=0.0;ldist=0.0;dpl.clear();
     while (!closed && i<numberpoints && count<numberpoints) {
       xdata[i]=dplist[count].GetX();
       ydata[i]=dplist[count].GetY();
       if (i>0) dist=sqrt(pow(xdata[i]-xdata[i-1],2)+pow(ydata[i]-ydata[i-1],2));
       if (dist>thresh) {closed=true;i--;count--;}
       else dpl.push_back(dplist[count]);
       i++;count++;
     }
     Dataset tempdat(dpl);
     if (i>2) result.push_back(tempdat);
   }
   delete[] xdata;xdata=NULL;
   delete[] ydata;ydata=NULL;
   return ok;
 }
 
 
 int Dataset::FindIndex(const double& xval) const {
   int result=0;
 
   if (xval<dplist[0].GetX() || xval>dplist[numberpoints-1].GetX()) return (-1);
   if (xval==dplist[numberpoints-1].GetX()) return (numberpoints-1);
   while (dplist[result].GetX()<=xval) result++;
   result--;
   return result;
 }
 int Dataset::NumberPeaks() const {
   Dataset diff=Differential();
   int sign=int(diff.dplist[0].GetY()/fabs(diff.dplist[0].GetY()));
   int lsign=sign;
   int npeaks=0,i;
   for (i=0;i<numberpoints;i++) {
     if (diff.dplist[i].GetY()!=0.0) sign=int(diff.dplist[i].GetY()/fabs(diff.dplist[i].GetY()));
     if (sign==-1 && lsign==1) {
       npeaks++;
     }
     lsign=sign;
   }
   return npeaks;
 }
 Dataset Dataset::Differential() const {
   Dataset result=*this;
   double dx;
   
   int i;
   for (i=1;i<numberpoints;i++) {
     Datapoint dum=dplist[i];
     dx=dum.GetX()-dplist[i-1].GetX();
     if (dx) result.dplist[i]=dum.sub_y(dplist[i-1])/dx;
   }
   result.dplist[0]=result.dplist[1];
   return result;
 }
 Dataset Dataset::Integrate() const
 {
   Dataset result=*this;
   double dx;
   
   int i;
   for (i=1;i<numberpoints;i++) {
     dx=dplist[i].GetX();
     result.dplist[i]=Integrate(dplist[0].GetX(),dx);
   }
   result.dplist[0]=0.0;
   return result;
 }
 Dataset Dataset::IntegrateCum(const double& a,const double& b) const
 {
   Dataset result=*this;
   double dx;
   
   int i;
   for (i=0;i<numberpoints;i++) {
     dx=dplist[i].GetX();
     if (dx>a && dx<b) result.dplist[i]=Integrate(a,dx);
     else result.dplist[i]=0.0;
   }
   return result;
 }
 Dataset Dataset::IntegrateInvCum(const double& a,const double& b) const
 {
   Dataset result=*this;
   double dx;
   
   int i;
   for (i=0;i<numberpoints;i++) {
     dx=dplist[i].GetX();
     if (dx>a && dx<b) result.dplist[i]=Integrate(dx,b);
     else result.dplist[i]=0.0;
     result.dplist[i].SetX(dx);
   }
   return result;
 }
 
 Dataset& Dataset::operator=(const Dataset &data)
 {
   int i;
 
   if (this!=&data) {
     Init(data.numberpoints);
     for (i=0;i<numberpoints;i++) {
       dplist[i]=data.dplist[i];
     }
   }
   comment=data.comment;
         return (*this);
 }
 
 Dataset operator/(const Dataset &data1,const double &number)
 {
   Dataset result;
   Datapoint dummy;
   int i,np;
 
   np=data1.numberpoints;
   result.Init(np);
   for (i=0;i<np;i++) {
     dummy=data1.dplist[i];
     result.dplist[i]=dummy.div_y(number);
   }
   return result;
 }
 Dataset operator/(const double& number,const Dataset &data1)
 {
   Dataset result;
   Datapoint dummy;
   int i,np;
 
   np=data1.numberpoints;
   result.Init(np);
   for (i=0;i<np;i++) {
     dummy=data1.dplist[i];
     double res=dummy.GetY();
     result.dplist[i].SetY(number/res);
   }
   return result;
 }
 Dataset operator*(const Dataset &data1,const double &number)
 {
   Dataset result;
   Datapoint dummy;
   int i,np;
 
   np=data1.numberpoints;
   result.Init(np);
   for (i=0;i<np;i++) {
     dummy=data1.dplist[i];
     result.dplist[i]=dummy.mul_y(number);
   }
   return result;
 }
 Dataset operator+(const Dataset &data1,const Dataset &data2)
 {
   double min,max;
   int i,np,j,k;
 
   np=data1.numberpoints;
   Dataset result(data1);
   for (i=0;i<np;i++) {
     k=data2.FindIndex(result.dplist[i].GetX());
     if (k!=-1) {
       result.dplist[i].add_y(data2.dplist[k]);
     }
   }
   return result;
 }
 
 Dataset Dataset::operator+(const Dataset &data2)
 {
   double min,max;
   int i,np,j,k;
         np=numberpoints;
         Dataset result(*this);
   for (i=0;i<np;i++) {
     k=data2.FindIndex(result.dplist[i].GetX());
     if (k!=-1) {
       result.dplist[i].add_y(data2.dplist[k]);
     }
   }
   return result;
 }
 
 Dataset Dataset::operator-(const Dataset &data2)
 {
   double min,max;
   int i,np,j,k;
         np=numberpoints;
         Dataset result(*this);
   for (i=0;i<np;i++) {
     k=data2.FindIndex(result.dplist[i].GetX());
     if (k!=-1) {
       result.dplist[i].sub_y(data2.dplist[k]);
     }
   }
   return result;
 }
 
 
 Dataset Dataset::operator/(const Dataset &data2)
 {
   double min,max;
   int i,np,j,k;
   np=numberpoints;
   Dataset result(*this);
   for (i=0;i<np;i++) {
     k=data2.FindIndex(result.dplist[i].GetX());
     if (k!=-1) {
       result.dplist[i].div_y(data2.dplist[k]);
     }
   }
   return result;
 }
 
 
 
 Dataset Dataset::operator*(const Dataset &data2)
 {
   double min,max;
   int i,np,j,k;
   np=numberpoints;
   Dataset result(*this);
   for (i=0;i<np;i++) {
     k=data2.FindIndex(result.dplist[i].GetX());
     if (k!=-1) {
       result.dplist[i].mul_y(data2.dplist[k]);
     }
   }
   return result;
 }
 
 Dataset Dataset::operator+(const double &number)
 {
   double min,max;
   int i,np,j,k;
         np=numberpoints;
         Dataset result(*this);
   for (i=0;i<np;i++) {
       result.dplist[i].add_y(number);
   }
   return result;
 }
 Dataset Dataset::operator-(const double &number)
 {
   double min,max;
   int i,np,j,k;
         np=numberpoints;
         Dataset result(*this);
   for (i=0;i<np;i++) {
       result.dplist[i].sub_y(number);
   }
   return result;
 }
 Dataset Dataset::operator/(const double &number)
 {
   double min,max;
   int i,np,j,k;
         np=numberpoints;
         Dataset result(*this);
   for (i=0;i<np;i++) {
       result.dplist[i].div_y(number);
   }
   return result;
 }
 Dataset Dataset::operator*(const double &number)
 {
   double min,max;
   int i,np,j,k;
         np=numberpoints;
         Dataset result(*this);
   for (i=0;i<np;i++) {
       result.dplist[i].mul_y(number);
   }
   return result;
 }
 
 valarray<double> Dataset::getY(void) const {
         int i;
         valarray<double> result(numberpoints);
         for (i=0;i<numberpoints;i++) {
                 result[i]=dplist[i].GetY();
         }
         return result;
 }
                 
 valarray<double> Dataset::getX(void) const {
         int i;
         valarray<double> result(numberpoints);
         for (i=0;i<numberpoints;i++) {
                 result[i]=dplist[i].GetX();
         }
         return result;
 }
                 
 Dataset operator-(const Dataset &data1,const Dataset &data2)
 {
   double min,max;
   int i,np,j,k;
 
   np=data1.numberpoints;
   Dataset result(data1);
   for (i=0;i<np;i++) {
     k=data2.FindIndex(result.dplist[i].GetX());
     if (k!=-1) {
       result.dplist[i].sub_y(data2.dplist[k]);
     }
   }
   return result;
 }
 Dataset operator/(const Dataset &data1,const Dataset &data2)
 {
   double min,max;
   int i,np,j,k;
 
   np=data1.numberpoints;
   Dataset result(data1);
   for (i=0;i<np;i++) {
     k=data2.FindIndex(result.dplist[i].GetX());
     if (k!=-1) {
       if (data2.dplist[k].GetY()!=0.0) result.dplist[i].div_y(data2.dplist[k]);
     }
   }
   return result;
 }
 Dataset operator*(const Dataset &data1,const Dataset &data2)
 {
   double min,max;
   int i,np,j,k;
 
   np=data1.numberpoints;
   Dataset result(data1);
   for (i=0;i<np;i++) {
     k=data2.FindIndex(result.dplist[i].GetX());
     if (k!=-1) {
       result.dplist[i].mul_y(data2.dplist[k]);
     }
   }
   return result;
 }
 
 double SigDiff(const Dataset &data1,const double& val)
 {
   double result=0.0;
   int i,np;
 
   np=data1.numberpoints;
   for (i=0;i<np;i++) {
     result+=data1.dplist[i].GetY()-val;
   }
   result=fabs(result)/np;
   return result;
 }
 
 int Dataset::SwapXY() {
   int i,np;
   double temp;
   np=numberpoints;
   for (i=0;i<np;i++) {
     dplist[i].SwapXY();
   }
   return ok;
 }
 int Dataset::SwapX() {
   int i,np;
   Datapoint temp;
   np=numberpoints;
   for (i=0;i<int(np/2);i++) {
     temp=dplist[i];
     dplist[i].SetX(dplist[np-1-i].GetX());
     dplist[i].SetXErr(dplist[np-1-i].GetXErr());
     dplist[np-1-i].SetX(temp.GetX());
     dplist[np-1-i].SetXErr(temp.GetXErr());
   }
   return ok;
 }
 int Dataset::SwapY() {
   int i,np;
   Datapoint temp;
   np=numberpoints;
   for (i=0;i<int(np/2);i++) {
     temp=dplist[i];
     dplist[i].SetY(dplist[np-1-i].GetY());
     dplist[i].SetYErr(dplist[np-1-i].GetYErr());
     dplist[np-1-i].SetY(temp.GetY());
     dplist[np-1-i].SetYErr(temp.GetYErr());
   }
   return ok;
 }
 Dataset Dataset::getXHistogram(const int& np,const double& x1,const double& x2,const valarray<double>& weights) const {
         Dataset result(np,x1,x2);
         int i;
         for (i=0;i<numberpoints;i++) {
                 result.AddTo(dplist[i].GetX(),weights[i]);
         }
         return result;
 }
 Dataset Dataset::getYHistogram(const int& np,const double& x1,const double& x2,const valarray<double>& weights) const {
         Dataset result(np,x1,x2);
         int i;
         for (i=0;i<numberpoints;i++) {
                 result.AddTo(dplist[i].GetY(),weights[i]);
         }
         return result;
 }
 Dataset Dataset::getXHistogram(const int& np,const double& x1,const double& x2) const {
         Dataset result(np,x1,x2);
         int i;
         for (i=0;i<numberpoints;i++) {
                 result.AddTo(dplist[i].GetX(),1.0);
         }
         return result;
 }
 Dataset Dataset::getYHistogram(const int& np,const double& x1,const double& x2) const {
         Dataset result(np,x1,x2);
         int i;
         for (i=0;i<numberpoints;i++) {
                 result.AddTo(dplist[i].GetY(),1.0);
         }
         return result;
 }
 
 int Dataset::loadAscii(const char* infile,const int& xi,const int& yi) {
         ifstream inf(infile);
         loadAscii(inf,xi,yi);
         return 0;
 }
 int Dataset::loadAscii(const string& infile,const int& xi,const int& yi) {
         ifstream inf(infile.c_str());
         loadAscii(inf,xi,yi);
         return 0;
 }
 istream& Dataset::loadAscii(istream& infile,const int& xi,const int& yi) {
         string line;
         int nrows=0;
         vector<double> dumvals;
         while (infile) {
                 getline(infile,line);
                 int lpos=0;
                 int npos=line.find(' ');
                 //cout << "npos=" << npos << flush;
                 dumvals.clear();
                 while (npos!=string::npos) {
                         string dum(line,lpos,npos-lpos);
                         dumvals.push_back(atof(dum.c_str()));
                         nrows++;
                         lpos=npos;
                         npos=line.find(' ',npos+1);
                         //cout << "npos=" << npos << flush;
                 }               
                 Datapoint dum(dumvals[xi],dumvals[yi]);
     //cout << dum << endl;
                 if (infile) dplist.push_back(dum);
         }
   numberpoints=dplist.size();
   return infile;
 }
 
 istream& Dataset::loadAscii(istream& infile) {
         while (infile) {
                 Datapoint dum;
                 infile >> dum;
     if (infile) dplist.push_back(dum);
         }
   numberpoints=dplist.size();
   return infile;
 }
 
 ostream& operator<<(ostream& out,const Dataset& dat) {
   int i;
   out << dat.comment << "\n";
   out << dat.numberpoints << "\n";
   for (i=0;i<dat.numberpoints;i++) {
     out << dat.dplist[i] << "\n";
   }
   return out;
 }
 istream& operator>>(istream& in_file,Dataset& dat) {
         string dummy;
         int i,np;
   getline(in_file,dummy);
         cout << dummy << "\n";
         cout << dat.GetComment() << "\n";
         dat.SetComment(dummy);
   in_file >> np;
   for (i=0;i<np;i++) {
     Datapoint dum;
     in_file >> dum;
     (dat.dplist).push_back(dum);
   }  
   dat.numberpoints=dat.dplist.size();
   return in_file;
 }
 Dataset operator*(const Dataset &data1,Function &funct)
 {
   Dataset result;
   Datapoint dummy;
   int i,np;
 
   np=data1.numberpoints;
   result.Init(np);
   for (i=0;i<np;i++) {
     dummy=data1.dplist[i];
     result.dplist[i]=dummy.mul_y(funct.function(dummy.GetX()));
   }
   return result;
 }
 Dataset operator/(const Dataset &data1,Function &funct)
 {
   Dataset result;
   Datapoint dummy;
   int i,np;
 
   np=data1.numberpoints;
   result.Init(np);
   for (i=0;i<np;i++) {
     dummy=data1.dplist[i];
     result.dplist[i]=dummy.div_y(funct.function(dummy.GetX()));
   }
   return result;
 }
 
 Dataset& Sqrt(Dataset& datin) {
     int i;
     for (i=0;i<datin.GetSize();i++) {
         datin.SetTo(i,sqrt(datin.GetValue(i)));
     }
     return datin;
 }
 
 
 /*
 
 *********************** Definitions of all routines that are using PGplot here .......
 
 void Dataset::Draw(PGGraph &output,int style,int width)
 {  
   float* xdata;
   float* ydata;
   int i;
 
   //output.SetLine(style,width);
   xdata=new float[numberpoints];
   ydata=new float[numberpoints];
   for (i=0;i<numberpoints;i++) {
     xdata[i]=dplist[i].GetX();
     ydata[i]=dplist[i].GetY();
   }
   cpgline(numberpoints,xdata,ydata);
   delete[] xdata;xdata=NULL;
   delete[] ydata;ydata=NULL;
 }
 void Dataset::DrawBinned(PGGraph &output,int style,int width)
 { 
   float* xdata;
   float* ydata;
   int i;
 
   //output.SetLine(style,width);
   xdata=new float[numberpoints];
   ydata=new float[numberpoints];
   for (i=0;i<numberpoints;i++) {
     xdata[i]=dplist[i].GetX();
     ydata[i]=dplist[i].GetY();
   }
   cpgbin(numberpoints,xdata,ydata,1);
   delete[] xdata;xdata=NULL;
   delete[] ydata;ydata=NULL;
 }
 
 void Dataset::DrawPoints(PGGraph &output,int marker)
 { 
   float* xdata;
   float* ydata;
   int i;
 
   xdata=new float[numberpoints];
   ydata=new float[numberpoints];
   for (i=0;i<numberpoints;i++) {
     xdata[i]=dplist[i].GetX();
     ydata[i]=dplist[i].GetY();
   }
   cpgpt(numberpoints,xdata,ydata,marker);
   delete[] xdata;xdata=NULL;
   delete[] ydata;ydata=NULL;
 }
 void Dataset::DrawXHistogram(PGGraph &output,int nbins)
 { 
   float* data;
   double dmin,dmax;
   int i,flag=1;
 
   if (nbins<0) {flag=0;nbins=20;}
   data=new float[numberpoints];
   dmin=dplist[0].GetX();
   dmax=dplist[0].GetX();
   for (i=0;i<numberpoints;i++) {
     data[i]=dplist[i].GetX();
     if (data[i]>dmax) dmax=data[i];
     if (data[i]<dmin) dmin=data[i];
   }
   cpghist(numberpoints,data,dmin,dmax,nbins,flag);
   delete[] data;data=NULL;
 }
 void Dataset::DrawYHistogram(PGGraph &output,int nbins)
 { 
   float* data;
   double dmin,dmax;
   int i,flag=1;
 
   if (nbins<0) {flag=0;nbins=20;}
   data=new float[numberpoints];
   dmin=dplist[0].GetY();
   dmax=dplist[0].GetY();
   for (i=0;i<numberpoints;i++) {
     data[i]=dplist[i].GetY();
     if (data[i]>dmax) dmax=data[i];
     if (data[i]<dmin) dmin=data[i];
   }
   cpghist(numberpoints,data,dmin,dmax,nbins,flag);
   delete[] data;data=NULL;
 }
 void Dataset::Display(PGGraph &output)
 {
   output.NewPage(dplist[0].GetX(),dplist[numberpoints-1].GetX(),GetMin().GetY(),GetMax().GetY(),0,0);
   Draw(output);
 }
 void Dataset::Display(PGGraph &output,int way)
 {
   int j,w;
   
   j=0;w=0;
   if (way & LOGX) w=10;
   if (way & LOGY) w=20;
   if ((way & LOGX) && (way & LOGY)) w=30;
   if (way & JUST) j=1;
   if (way & XHISTOGRAM) DrawXHistogram(output,-1);
   if (way & YHISTOGRAM) DrawYHistogram(output,-1);
   if (way & POINT) {
       output.NewPage(dplist[0].GetX(),dplist[numberpoints-1].GetX(),GetMin().GetY(),GetMax().GetY(),j,w);
     DrawPoints(output);
     return;
   }
   if (way & BINNED) {
       output.NewPage(dplist[0].GetX(),dplist[numberpoints-1].GetX(),GetMin().GetY(),GetMax().GetY(),j,w);
     DrawBinned(output);
   }
   else {
       output.NewPage(dplist[0].GetX(),dplist[numberpoints-1].GetX(),GetMin().GetY(),GetMax().GetY(),j,w);
     Draw(output);
   }
 }
 void Dataset::Display(PGGraph &output,int way,char* title,char *labelx,char* labely,int font,float height,int style,int width) {
   int j,w;
   
   j=0;w=0;
   if (way & LOGX) w=10;
   if (way & LOGY) w=20;
   if ((way & LOGX) && (way & LOGY)) w=30;
   if (way & JUST) j=1;
   output.SetFont(font,height);
   output.SetLine(FULL,1);
   output.NewPage(dplist[0].GetX(),dplist[numberpoints-1].GetX(),GetMin().GetY(),GetMax().GetY(),j,w,title,labelx,labely);
   Draw(output,style,width);
 }  
 
 int Dataset::DrawStdDev(double stdev) {
   int i,np;
   float* xdata;
   float* topdata;
   float* botdata;
 
   //output.SetLine(style,width);
   xdata=new float[numberpoints];
   topdata=new float[numberpoints];
   botdata=new float[numberpoints];
   for (i=0;i<numberpoints;i++) {
     xdata[i]=dplist[i].GetX();
     topdata[i]=dplist[i].GetY()+stdev;
     botdata[i]=dplist[i].GetY()-stdev;
   }
   cpgerry(numberpoints,xdata,topdata,botdata,1.0);
   delete[] xdata;xdata=NULL;
   delete[] topdata;topdata=NULL;
   delete[] botdata;botdata=NULL;
 }
 
 int Dataset::DrawStdDev(Function& fin) {
   int i,np;
   float* xdata;
   float* topdata;
   float* botdata;
 
   //output.SetLine(style,width);
   xdata=new float[numberpoints];
   topdata=new float[numberpoints];
   botdata=new float[numberpoints];
   for (i=0;i<numberpoints;i++) {
     xdata[i]=dplist[i].GetX();
     topdata[i]=dplist[i].GetY()+fin.function(dplist[i].GetY());
     botdata[i]=dplist[i].GetY()-fin.function(dplist[i].GetY());
   }
   cpgerry(numberpoints,xdata,topdata,botdata,1.0);
   delete[] xdata;xdata=NULL;
   delete[] topdata;topdata=NULL;
   delete[] botdata;botdata=NULL;
 }
 
 
 */
