NcSample.cxx

Go to the documentation of this file.

 
#include "NcSample.h"
#include "Riostream.h"
 
ClassImp(NcSample); // Class implementation to enable ROOT I/O
 
NcSample::NcSample(const char* name,const char* title) : TNamed(name,title)
{
 
 fDim=fMaxdim;
 fN=0;
 fRemove=0;
 fStore=0;
 fNmax=0;
 fMoveVariable=0;
 fX=0;
 fY=0;
 fZ=0;
 fT=0;
 fArr=0;
 fIndices=0;
 fOrdered=0;
 fGraphT=0;
 fCanvas=0;
 fAnimObject=0;
 Reset();
 SetNames();
}

NcSample::~NcSample()
{
 
 if (fX)
 {
  delete fX;
  fX=0;
 }
 if (fY)
 {
  delete fY;
  fY=0;
 }
 if (fZ)
 {
  delete fZ;
  fZ=0;
 }
 if (fT)
 {
  delete fT;
  fT=0;
 }
 if (fArr)
 {
  delete fArr;
  fArr=0;
 }
 if (fIndices)
 {
  delete fIndices;
  fIndices=0;
 }
 if (fGraphT)
 {
  delete fGraphT;
  fGraphT=0;
 }
 if (fCanvas)
 {
  delete fCanvas;
  fCanvas=0;
 }
 if (fAnimObject)
 {
  delete fAnimObject;
  fAnimObject=0;
 }
}

NcSample::NcSample(const NcSample& s) : TNamed(s)
{
 
 for (Int_t i=0; i<s.fMaxdim; i++)
 {
  fNames[i]=s.fNames[i];
 }
 
 fDim=s.fDim;
 fStore=s.fStore;
 fNmax=s.fNmax;
 fMoveVariable=s.fMoveVariable;
 fN=s.fN;
 fRemove=s.fRemove;
 
 fX=0;
 fY=0;
 fZ=0;
 fT=0;
 fArr=0;
 fIndices=0;
 fOrdered=0;
 fGraphT=0;
 fCanvas=0;
 fAnimObject=0;
 
 if (s.fX) fX=new TArrayD(fN);
 if (s.fY) fY=new TArrayD(fN);
 if (s.fZ) fZ=new TArrayD(fN);
 if (s.fT) fT=new TArrayD(fN);
 
 Double_t val=0;
 for (Int_t i=0; i<fDim; i++)
 {
  fSum[i]=s.fSum[i];
  fMean[i]=s.fMean[i];
  fRMSdev[i]=s.fRMSdev[i];
  fSigma[i]=s.fSigma[i];
  fMin[i]=s.fMin[i];
  fMax[i]=s.fMax[i];
  for (Int_t j=0; j<fDim; j++)
  {
   fSum2[i][j]=s.fSum2[i][j];
   fCov[i][j]=s.fCov[i][j];
   fCor[i][j]=s.fCor[i][j];
  }
 }
 
 if (!fStore) return;
 
 for (Int_t i=0; i<fN; i++)
 {
  if (fX)
  {
   val=s.fX->At(i);
   fX->AddAt(val,i);
  }
  if (fY)
  {
   val=s.fY->At(i);
   fY->AddAt(val,i);
  }
  if (fZ)
  {
   val=s.fZ->At(i);
   fZ->AddAt(val,i);
  }
  if (fT)
  {
   val=s.fT->At(i);
   fT->AddAt(val,i);
  }
 }
}

void NcSample::Reset()
{
 
 fN=0;
 fRemove=0;
 fOrdered=0;
 for (Int_t i=0; i<fDim; i++)
 {
  fSum[i]=0;
  fMean[i]=0;
  fRMSdev[i]=0;
  fSigma[i]=0;
  fMin[i]=0;
  fMax[i]=0;
  for (Int_t j=0; j<fDim; j++)
  {
   fSum2[i][j]=0;
   fCov[i][j]=0;
   fCor[i][j]=0;
  }
 }
 
 // Set storage arrays to initial size
 if (fX) fX->Set(10);
 if (fY) fY->Set(10);
 if (fZ) fZ->Set(10);
 if (fT) fT->Set(10);
 
 // Resetting the variable names to their defaults
 SetNames();
 
 // Delete the temp. storage arrays for ordering
 if (fArr)
 {
  delete fArr;
  fArr=0;
 }
 if (fIndices)
 {
  delete fIndices;
  fIndices=0;
 }
 
 // Delete the temp. TGraphTime pointer
 if (fGraphT)
 {
  delete fGraphT;
  fGraphT=0;
 }
 
 // Delete the multi-purpose canvas
 if (fCanvas)
 {
  delete fCanvas;
  fCanvas=0;
 }
 
 // Delete the multi-purpose animation object
 if (fAnimObject)
 {
  delete fAnimObject;
  fAnimObject=0;
 }
}

void NcSample::SetNames(TString name1,TString name2,TString name3,TString name4)
{
 
 if (!fN)
 {
  fNames[0]=name1;
  fNames[1]=name2;
  fNames[2]=name3;
  fNames[3]=name4;
 }
 else
 {
  cout << " *NcSample::SetNames* Names not modified, since data were already present." << endl;
 }
}

void NcSample::Enter(Double_t x)
{
 
 if (!fN)
 {
  fDim=1;
  fNames[1]="";
  fNames[2]="";
  fNames[3]="";
 }
 
 if (fDim != 1)
 {
  cout << " *NcSample::Enter* Error : Not a 1-dim sample." << endl;
  return;
 }
 
 // Remove (after re-ordering) the first entry if needed
 if (fNmax && fN==fNmax) RemoveEntry(1,fMoveVariable,abs(fMoveVariable));
 
 fN+=1;
 fSum[0]+=x;
 fSum2[0][0]+=x*x;
 
 if (fN==1)
 {
  fMin[0]=x;
  fMax[0]=x;
 }
 else
 {
  if (x<fMin[0]) fMin[0]=x;
  if (x>fMax[0]) fMax[0]=x;
 }
 
 // Store all entered data when storage mode has been selected 
 if (fStore)
 {
  if (!fX) fX=new TArrayD(10);
  if (fX->GetSize() < fN) fX->Set(fN+10);
  fX->AddAt(x,fN-1);
 }
 
 // Compute the various statistics
 Compute();
}

void NcSample::Remove(Double_t x)
{
 
 if (!fN) return;
 
 if (fDim != 1)
 {
  cout << " *NcSample::Remove* Error : Not a 1-dim sample." << endl;
  return;
 }
 
 fRemove=1;
 fN-=1;
 fSum[0]-=x;
 fSum2[0][0]-=x*x;
 
 // Remove data entry from the storage
 if (fStore && fX)
 {
  Double_t val=0;
  for (Int_t i=0; i<=fN; i++)
  {
   val=fX->At(i);
   if (fabs(x-val)>1.e-10) continue;
 
   for (Int_t j=i+1; j<=fN; j++)
   {
    val=fX->At(j);
    fX->AddAt(val,j-1);
   }
  }
 }
 
 // Compute the various statistics
 Compute();
}

void NcSample::Enter(Double_t x,Double_t y)
{
 
 if (!fN)
 {
  fDim=2;
  fNames[2]="";
  fNames[3]="";
 }
 
 if (fDim != 2)
 {
  cout << " *NcSample::Enter* Error : Not a 2-dim sample." << endl;
  return;
 }
 
 // Remove (after re-ordering) the first entry if needed
 if (fNmax && fN==fNmax) RemoveEntry(1,fMoveVariable,abs(fMoveVariable));
 
 fN+=1;
 fSum[0]+=x;
 fSum[1]+=y;
 fSum2[0][0]+=x*x;
 fSum2[0][1]+=x*y;
 fSum2[1][0]+=y*x;
 fSum2[1][1]+=y*y;
 
 if (fN==1)
 {
  fMin[0]=x;
  fMax[0]=x;
  fMin[1]=y;
  fMax[1]=y;
 }
 else
 {
  if (x<fMin[0]) fMin[0]=x;
  if (x>fMax[0]) fMax[0]=x;
  if (y<fMin[1]) fMin[1]=y;
  if (y>fMax[1]) fMax[1]=y;
 }
 
 // Store all entered data when storage mode has been selected 
 if (fStore)
 {
  if (!fX) fX=new TArrayD(10);
  if (fX->GetSize() < fN) fX->Set(fN+10);
  fX->AddAt(x,fN-1);
  if (!fY) fY=new TArrayD(10);
  if (fY->GetSize() < fN) fY->Set(fN+10);
  fY->AddAt(y,fN-1);
 }
 
 // Compute the various statistics
 Compute();
}

void NcSample::Remove(Double_t x,Double_t y)
{
 
 if (!fN) return;
 
 if (fDim != 2)
 {
  cout << " *NcSample::Remove* Error : Not a 2-dim sample." << endl;
  return;
 }
 
 fRemove=1;
 fN-=1;
 fSum[0]-=x;
 fSum[1]-=y;
 fSum2[0][0]-=x*x;
 fSum2[0][1]-=x*y;
 fSum2[1][0]-=y*x;
 fSum2[1][1]-=y*y;
 
 // Remove data entry from the storage
 if (fStore && fX && fY)
 {
  Double_t val=0;
  for (Int_t i=0; i<=fN; i++)
  {
   val=fX->At(i);
   if (fabs(x-val)>1.e-10) continue;
   val=fY->At(i);
   if (fabs(y-val)>1.e-10) continue;
 
   for (Int_t j=i+1; j<=fN; j++)
   {
    val=fX->At(j);
    fX->AddAt(val,j-1);
    val=fY->At(j);
    fY->AddAt(val,j-1);
   }
  }
 }
 
 // Compute the various statistics
 Compute();
}

void NcSample::Enter(Double_t x,Double_t y,Double_t z)
{
 
 if (!fN)
 {
  fDim=3;
  fNames[3]="";
 }
 
 if (fDim != 3)
 {
  cout << " *NcSample::Enter* Error : Not a 3-dim sample." << endl;
  return;
 }
 
 // Remove (after re-ordering) the first entry if needed
 if (fNmax && fN==fNmax) RemoveEntry(1,fMoveVariable,abs(fMoveVariable));
 
 fN+=1;
 fSum[0]+=x;
 fSum[1]+=y;
 fSum[2]+=z;
 fSum2[0][0]+=x*x;
 fSum2[0][1]+=x*y;
 fSum2[0][2]+=x*z;
 fSum2[1][0]+=y*x;
 fSum2[1][1]+=y*y;
 fSum2[1][2]+=y*z;
 fSum2[2][0]+=z*x;
 fSum2[2][1]+=z*y;
 fSum2[2][2]+=z*z;
 
 if (fN==1)
 {
  fMin[0]=x;
  fMax[0]=x;
  fMin[1]=y;
  fMax[1]=y;
  fMin[2]=z;
  fMax[2]=z;
 }
 else
 {
  if (x<fMin[0]) fMin[0]=x;
  if (x>fMax[0]) fMax[0]=x;
  if (y<fMin[1]) fMin[1]=y;
  if (y>fMax[1]) fMax[1]=y;
  if (z<fMin[2]) fMin[2]=z;
  if (z>fMax[2]) fMax[2]=z;
 }
 
 // Store all entered data when storage mode has been selected 
 if (fStore)
 {
  if (!fX) fX=new TArrayD(10);
  if (fX->GetSize() < fN) fX->Set(fN+10);
  fX->AddAt(x,fN-1);
  if (!fY) fY=new TArrayD(10);
  if (fY->GetSize() < fN) fY->Set(fN+10);
  fY->AddAt(y,fN-1);
  if (!fZ) fZ=new TArrayD(10);
  if (fZ->GetSize() < fN) fZ->Set(fN+10);
  fZ->AddAt(z,fN-1);
 }
 
 // Compute the various statistics
 Compute();
}

void NcSample::Remove(Double_t x,Double_t y,Double_t z)
{
 
 if (!fN) return;
 
 if (fDim != 3)
 {
  cout << " *NcSample::Remove* Error : Not a 3-dim sample." << endl;
  return;
 }
 
 fRemove=1;
 fN-=1;
 fSum[0]-=x;
 fSum[1]-=y;
 fSum[2]-=z;
 fSum2[0][0]-=x*x;
 fSum2[0][1]-=x*y;
 fSum2[0][2]-=x*z;
 fSum2[1][0]-=y*x;
 fSum2[1][1]-=y*y;
 fSum2[1][2]-=y*z;
 fSum2[2][0]-=z*x;
 fSum2[2][1]-=z*y;
 fSum2[2][2]-=z*z;
 
 // Remove data entry from the storage
 if (fStore && fX && fY && fZ)
 {
  Double_t val=0;
  for (Int_t i=0; i<=fN; i++)
  {
   val=fX->At(i);
   if (fabs(x-val)>1.e-10) continue;
   val=fY->At(i);
   if (fabs(y-val)>1.e-10) continue;
   val=fZ->At(i);
   if (fabs(z-val)>1.e-10) continue;
 
   for (Int_t j=i+1; j<=fN; j++)
   {
    val=fX->At(j);
    fX->AddAt(val,j-1);
    val=fY->At(j);
    fY->AddAt(val,j-1);
    val=fZ->At(j);
    fZ->AddAt(val,j-1);
   }
  }
 }
 
 // Compute the various statistics
 Compute();
}

void NcSample::Enter(Double_t x,Double_t y,Double_t z,Double_t t)
{
 
 if (!fN) fDim=4;
 
 if (fDim != 4)
 {
  cout << " *NcSample::Enter* Error : Not a 4-dim sample." << endl;
  return;
 }
 
 // Remove (after re-ordering) the first entry if needed
 if (fNmax && fN==fNmax) RemoveEntry(1,fMoveVariable,abs(fMoveVariable));
 
 fN+=1;
 fSum[0]+=x;
 fSum[1]+=y;
 fSum[2]+=z;
 fSum[3]+=t;
 fSum2[0][0]+=x*x;
 fSum2[0][1]+=x*y;
 fSum2[0][2]+=x*z;
 fSum2[0][3]+=x*t;
 fSum2[1][0]+=y*x;
 fSum2[1][1]+=y*y;
 fSum2[1][2]+=y*z;
 fSum2[1][3]+=y*t;
 fSum2[2][0]+=z*x;
 fSum2[2][1]+=z*y;
 fSum2[2][2]+=z*z;
 fSum2[2][3]+=z*t;
 fSum2[3][0]+=t*x;
 fSum2[3][1]+=t*y;
 fSum2[3][2]+=t*z;
 fSum2[3][3]+=t*t;
 
 if (fN==1)
 {
  fMin[0]=x;
  fMax[0]=x;
  fMin[1]=y;
  fMax[1]=y;
  fMin[2]=z;
  fMax[2]=z;
  fMin[3]=t;
  fMax[3]=t;
 }
 else
 {
  if (x<fMin[0]) fMin[0]=x;
  if (x>fMax[0]) fMax[0]=x;
  if (y<fMin[1]) fMin[1]=y;
  if (y>fMax[1]) fMax[1]=y;
  if (z<fMin[2]) fMin[2]=z;
  if (z>fMax[2]) fMax[2]=z;
  if (t<fMin[3]) fMin[3]=t;
  if (t>fMax[3]) fMax[3]=t;
 }
 
 // Store all entered data when storage mode has been selected 
 if (fStore)
 {
  if (!fX) fX=new TArrayD(10);
  if (fX->GetSize() < fN) fX->Set(fN+10);
  fX->AddAt(x,fN-1);
  if (!fY) fY=new TArrayD(10);
  if (fY->GetSize() < fN) fY->Set(fN+10);
  fY->AddAt(y,fN-1);
  if (!fZ) fZ=new TArrayD(10);
  if (fZ->GetSize() < fN) fZ->Set(fN+10);
  fZ->AddAt(z,fN-1);
  if (!fT) fT=new TArrayD(10);
  if (fT->GetSize() < fN) fT->Set(fN+10);
  fT->AddAt(t,fN-1);
 }
 
 // Compute the various statistics
 Compute();
}

void NcSample::Remove(Double_t x,Double_t y,Double_t z,Double_t t)
{
 
 if (!fN) return;
 
 if (fDim != 4)
 {
  cout << " *NcSample::Remove* Error : Not a 4-dim sample." << endl;
  return;
 }
 
 fRemove=1;
 fN-=1;
 fSum[0]-=x;
 fSum[1]-=y;
 fSum[2]-=z;
 fSum[3]-=t;
 fSum2[0][0]-=x*x;
 fSum2[0][1]-=x*y;
 fSum2[0][2]-=x*z;
 fSum2[0][3]-=x*t;
 fSum2[1][0]-=y*x;
 fSum2[1][1]-=y*y;
 fSum2[1][2]-=y*z;
 fSum2[1][3]-=y*t;
 fSum2[2][0]-=z*x;
 fSum2[2][1]-=z*y;
 fSum2[2][2]-=z*z;
 fSum2[2][3]-=z*t;
 fSum2[3][0]-=t*x;
 fSum2[3][1]-=t*y;
 fSum2[3][2]-=t*z;
 fSum2[3][3]-=t*t;
 
 // Remove data entry from the storage
 if (fStore && fX && fY && fZ && fT)
 {
  Double_t val=0;
  for (Int_t i=0; i<=fN; i++)
  {
   val=fX->At(i);
   if (fabs(x-val)>1.e-10) continue;
   val=fY->At(i);
   if (fabs(y-val)>1.e-10) continue;
   val=fZ->At(i);
   if (fabs(z-val)>1.e-10) continue;
   val=fT->At(i);
   if (fabs(t-val)>1.e-10) continue;
 
   for (Int_t j=i+1; j<=fN; j++)
   {
    val=fX->At(j);
    fX->AddAt(val,j-1);
    val=fY->At(j);
    fY->AddAt(val,j-1);
    val=fZ->At(j);
    fZ->AddAt(val,j-1);
    val=fT->At(j);
    fT->AddAt(val,j-1);
   }
  }
 }
 
 // Compute the various statistics
 Compute();
}

void NcSample::RemoveEntry(Int_t i,Int_t mode,Int_t j)
{
 
 if (fDim<1 || fN<=0 || i<1 || i>fN) return;
 
 if (mode && (j<1 || j>fDim))
 {
  cout << " *NcSample::RemoveEntry* Error : Invalide argument j=" << j << endl;
  return;
 }
 
 if (!fStore)
 {
  cout << " *NcSample::RemoveEntry* Error : Storage of data entries was not activated." << endl;
  return;
 }
 
 Order(mode,j);
 
 // Get the corresponding original entry index
 Int_t idx=fIndices->At(i-1);
 
 Double_t x=0;
 Double_t y=0;
 Double_t z=0;
 Double_t t=0;
 
 if (fX) x=GetEntry(idx+1,1);
 if (fY) y=GetEntry(idx+1,2);
 if (fZ) z=GetEntry(idx+1,3);
 if (fT) t=GetEntry(idx+1,4);
 
 if (fDim==1) Remove(x);
 if (fDim==2) Remove(x,y);
 if (fDim==3) Remove(x,y,z);
 if (fDim==4) Remove(x,y,z,t);
}

void NcSample::RemoveEntry(Int_t i,Int_t mode,TString name)
{
 
 Int_t j=GetIndex(name);
 RemoveEntry(i,mode,j);
}

void NcSample::RemoveEntry(Int_t i,Int_t j,Int_t mode,Int_t k)
{
 
 for (Int_t m=i; m<=j; m++)
 {
  RemoveEntry(i,mode,k);
 }
}

void NcSample::RemoveEntry(Int_t i,Int_t j,Int_t mode,TString name)
{
 
 for (Int_t m=i; m<=j; m++)
 {
  RemoveEntry(i,mode,name);
 }
}

void NcSample::Order(Int_t mode,Int_t i)
{
 
 if (mode && (i<1 || i>fDim))
 {
  cout << " *NcSample::Order* Error : Invalide argument i=" << i << endl;
  return;
 }
 
 if (fDim<1)
 {
  cout << " *NcSample::Order* Error : Dimension less than 1." << endl;
  return;
 }
 
 if (!fStore)
 {
  cout << " *NcSample::Order* Error : Storage of data entries was not activated." << endl;
  return;
 }
 
 if (fN<=0)
 {
  fOrdered=0;
  return;
 }
 
 // Use i=1 for mode=0 to obtain a correct setting of the ordering status word.
 if (!mode) i=1;
 
 // Set the corresponding ordering status word
 Int_t iword=10*abs(mode)+i;
 if (mode<0) iword*=-1;
 
 // No new ordering needed if the ordering status word hasn't changed
 if (iword==fOrdered) return;
 
 // Store the new ordering status word
 fOrdered=iword;
 
 // Prepare temp. array to hold the ordered values
 if (!fArr)
 {
  fArr=new TArrayD(fN);
 }
 else
 {
  if (fArr->GetSize() < fN) fArr->Set(fN);
 }
 
 // Prepare temp. array to hold the ordered indices
 if (!fIndices)
 {
  fIndices=new TArrayI(fN);
 }
 else
 {
  if (fIndices->GetSize() < fN) fIndices->Set(fN);
 }
 
 Double_t val=0;
 
 // Just the order in which the entries were entered
 if (!mode)
 {
  for (Int_t j=0; j<fN; j++)
  {
   if (i==1) val=fX->At(j);
   if (i==2) val=fY->At(j);
   if (i==3) val=fZ->At(j);
   if (i==4) val=fT->At(j);
   fArr->AddAt(val,j);
   fIndices->AddAt(j,j);
  }
  return;
 }
 
 // Order the values of the specified variable
 Int_t iadd=0;
 for (Int_t j=0; j<fN; j++)
 {
  if (i==1) val=fX->At(j);
  if (i==2) val=fY->At(j);
  if (i==3) val=fZ->At(j);
  if (i==4) val=fT->At(j);
 
  iadd=0;
  if (j==0)
  {
   fArr->AddAt(val,j);
   fIndices->AddAt(j,j);
   iadd=1;
  }
  else
  {
   for (Int_t k=0; k<j; k++)
   {
    if (mode>0 && val>=fArr->At(k)) continue; // Increasing ordering
    if (mode<0 && val<=fArr->At(k)) continue; // Decreasing ordering
    // Put value in between the existing ones
    for (Int_t m=j-1; m>=k; m--)
    {
     fArr->AddAt(fArr->At(m),m+1);
     fIndices->AddAt(fIndices->At(m),m+1);
    }
    fArr->AddAt(val,k);
    fIndices->AddAt(j,k);
    iadd=1;
    break;
   }
 
   if (!iadd)
   {
    fArr->AddAt(val,j);
    fIndices->AddAt(j,j);
   }
  }
 }
}

Int_t NcSample::GetIndex(TString name) const
{
 
 Int_t idx=0;
 
 for (Int_t i=0; i<fMaxdim; i++)
 {
  if (fNames[i]==name) idx=i+1;
 }
 
 return idx;
}

void NcSample::Compute()
{
 
 // Reset the ordering status word
 fOrdered=0;
 
 if (fN<=0) return;
 
 Double_t rn=fN;
 Double_t rn1=fN-1;
 Double_t var=0;
 Double_t rmsij;
 for (Int_t k=0; k<fDim; k++)
 {
  fMean[k]=fSum[k]/rn;
  var=(fSum2[k][k]/rn)-(fMean[k]*fMean[k]);
  if (var<0) var=0;
  fRMSdev[k]=sqrt(var);
  var=0;
  if (rn1)
  {
   var=(fSum2[k][k]-(rn*fMean[k]*fMean[k]))/rn1;
   if (var<0) var=0;
  }
  fSigma[k]=sqrt(var);
 }
 for (Int_t i=0; i<fDim; i++)
 {
  for (Int_t j=0; j<fDim; j++)
  {
   fCov[i][j]=(fSum2[i][j]/rn)-(fMean[i]*fMean[j]);
   fCor[i][j]=0;
   rmsij=fRMSdev[i]*fRMSdev[j];
   if (rmsij != 0) fCor[i][j]=fCov[i][j]/rmsij;
  }
 }
}

Int_t NcSample::GetDimension() const
{
 
 return fDim;
}

Int_t NcSample::GetN() const
{
 
 return fN;
}

TString NcSample::GetVariableName(Int_t i) const
{
 
 TString name="";
 
 if (i<1 || i>fDim)
 {
  cout << " *NcSample::GetVariableName* Error : Invalid index " << i << endl;
 }
 else
 {
  name=fNames[i-1];
 }
 
 return name;
}

Double_t NcSample::GetSum(Int_t i) const
{
 
 if (i<1 || i>fDim)
 {
  cout << " *NcSample::GetSum* Error : Invalid index " << i << endl;
  return 0;
 }
 else
 {
  return fSum[i-1];
 }
}

Double_t NcSample::GetSum(TString name) const
{
 
 Int_t i=GetIndex(name);
 return GetSum(i);
}

Double_t NcSample::GetMean(Int_t i) const
{
 
 if (i<1 || i>fDim)
 {
  cout << " *NcSample::GetMean* Error : Invalid index " << i << endl;
  return 0;
 }
 
 return fMean[i-1];
}

Double_t NcSample::GetMean(TString name) const
{
 
 Int_t i=GetIndex(name);
 return GetMean(i);
}

Double_t NcSample::GetRMS(Int_t i) const
{
 
 if (i<1 || i>fDim)
 {
  cout << " *NcSample::GetRMS* Error : Invalid index " << i << endl;
  return 0;
 }
 
 Double_t rms=0;
 Double_t val=fSum2[i-1][i-1]/float(fN);
 if (val>=0) rms=sqrt(val);
 
 return rms;
}

Double_t NcSample::GetRMS(TString name) const
{
 
 Int_t i=GetIndex(name);
 return GetRMS(i);
}

Double_t NcSample::GetVar(Int_t i,Int_t model) const
{
 
 if (i<1 || i>fDim)
 {
  cout << " *NcSample::GetVar* Error : Invalid index " << i << endl;
  return 0;
 }
 
 Double_t var=fRMSdev[i-1]*fRMSdev[i-1];
 if (model) var=fSigma[i-1]*fSigma[i-1];
 
 return var;
}

Double_t NcSample::GetVar(TString name,Int_t model) const
{
 
 Int_t i=GetIndex(name);
 return GetVar(i);
}

Double_t NcSample::GetSigma(Int_t i,Int_t model) const
{
 
 if (i<1 || i>fDim)
 {
  cout << " *NcSample::GetSigma* Error : Invalid index " << i << endl;
  return 0;
 }
 
 Double_t val=fRMSdev[i-1];
 if (model) val=fSigma[i-1];
 
 return val;
}

Double_t NcSample::GetSigma(TString name,Int_t model) const
{
 
 Int_t i=GetIndex(name);
 return GetSigma(i,model);
}

Double_t NcSample::GetCov(Int_t i,Int_t j) const
{
 
 if (i<1 || i>fDim || j<1 || j>fDim)
 {
  cout << " *NcSample::GetCov* Error : Invalid index encountered i=" << i << " j=" << j << endl;
  return 0;
 }
 
 return fCov[i-1][j-1];
}

Double_t NcSample::GetCov(TString nameA, TString nameB) const
{
 
 Int_t i=GetIndex(nameA);
 Int_t j=GetIndex(nameB);
 return GetCov(i,j);
}

Double_t NcSample::GetCor(Int_t i,Int_t j) const
{
 
 if (i<1 || i>fDim || j<1 || j>fDim)
 {
  cout << " *NcSample::GetCor* Error : Invalid index encountered i=" << i << " j=" << j << endl;
  return 0;
 }
 
  return fCor[i-1][j-1];
}

Double_t NcSample::GetCor(TString nameA, TString nameB) const
{
 
 Int_t i=GetIndex(nameA);
 Int_t j=GetIndex(nameB);
 return GetCor(i,j);
}

void NcSample::Data(Int_t i,Int_t j)
{
 
 const char* name=GetName();
 const char* title=GetTitle();
 
 cout << " *" << ClassName() << "::Data*";
 if (strlen(name))  cout << " Info for sample Name : " << name;
 if (strlen(title)) cout << " Title : " << title;
 cout << endl;
 cout << " Access to the individual data entries is";
 if (!fStore) cout << " not";
 cout << " possible (StoreMode=" << fStore << ")" << endl;
 
 if (fStore && fNmax)
 {
  if (fNmax==fN && fRemove)
  {
   cout << " *** At entering data the dynamic FIFO storage was limited to " << fNmax << " entries";
  }
  else
  {
   cout << " *** At entering new data the dynamic FIFO storage will be limited to " << fNmax << " entries";
  }
  if (!fMoveVariable) cout << " using the order of original entering.";
  if (fMoveVariable>0) cout << " after increasing";
  if (fMoveVariable<0) cout << " after decreasing";
  if (fMoveVariable) cout << " ordering w.r.t. variable " << abs(fMoveVariable);
  cout << endl;
 }
 
 if (i<0 || i>fDim || j<0 || j>fDim)
 {
  cout << " Inconsistent input i=" << i << " and j=" << j << " for dimension " << fDim << endl;
  return;
 }
 if (!i && !j) cout << " Statistics and correlations of all variables";
 if (i && !j)  cout << " Statistics of variable " << i;
 if (!i && j)  cout << " Statistics of variable " << j;
 if (i && j)   cout << " Correlation statistics of the variables " << i << " and " << j;
 cout << endl;
 
 if (!fN)
 {
  cout << " No data has been entered." << endl;
  return;
 }
 
 // Statistics and correlations of all variables
 if (!i && !j)
 {
  // Statistics of individual variables
  for (Int_t iv=1; iv<=fDim; iv++)
  {
   List(iv);
  }
 
  if (fDim<2) return;
 
  // Covariance(s) and correlation coefficient(s)
  for (Int_t iv=1; iv<=fDim; iv++)
  {
   for (Int_t jv=iv+1;jv<=fDim; jv++)
   {
    List(iv,jv);
   }
  }
 }
 
 // Statistics of the i-th variable
 if (i && !j) List(i);
 
 // Statistics of the j-th variable
 if (!i && j) List(j);
 
 // Correlation statistics of the variables i and j
 if (i && j) List(i,j);
}

void NcSample::Data(TString nameA,TString nameB)
{
 
 Int_t i=GetIndex(nameA);
 Int_t j=GetIndex(nameB);
 Data(i,j);
}

void NcSample::ListOrdered(Int_t mode,Int_t i)
{
 
 if (!fStore)
 {
  cout << " *NcSample::ListOrdered* Error : Storage of data entries was not activated." << endl;
  return;
 }
 
 if (fN<=0)
 {
  cout << " *NcSample::ListOrdered* No entries were stored." << endl;
  return;
 }
 
 if (fDim<1)
 {
  cout << " *NcSample::ListOrdered* Error : Dimension less than 1." << endl;
  return;
 }
 
 if (mode && (i<1 || i>fDim))
 {
  cout << " *NcSample::ListOrdered* Error : Invalid argument i=" << i << endl;
  return;
 }
 
 Order(mode,i);
 
 const char* name=GetName();
 const char* title=GetTitle();
 
 cout << " *NcSample::ListOrdered* Ordered listing for the sample";
 if (strlen(name))  cout << " Name : " << name;
 if (strlen(title)) cout << " Title : " << title;
 cout << endl;
 if (fNmax==fN && fRemove)
 {
  cout << " *** At entering data the dynamic FIFO storage was limited to " << fNmax << " entries";
  if (!fMoveVariable) cout << " using the order of original entering.";
  if (fMoveVariable>0) cout << " after increasing";
  if (fMoveVariable<0) cout << " after decreasing";
  if (fMoveVariable) cout << " ordering w.r.t. variable " << abs(fMoveVariable);
  cout << endl;
 }
 cout << " Listing of the stored entries in";
 if (!mode) cout << " order of original entering." << endl;
 if (mode<0) cout << " decreasing order of variable : " << i << " (" << fNames[i-1] << ")" << endl;
 if (mode>0) cout << " increasing order of variable : " << i << " (" << fNames[i-1] << ")" << endl;
 if (mode) cout << " The number between brackets indicates the original data entry number." << endl; 
 
 Int_t index=0;
 for (Int_t j=0; j<fN; j++)
 {
  index=fIndices->At(j);
  
  if (index<0 || index>=fN) continue;
 
  cout << " Index : " << (j+1);
  if (mode) cout << " (" << (index+1) << ") ";
  cout << " " << fNames[0] << "=" << fX->At(index);
  if (fDim>1) cout << " " << fNames[1] << "=" << fY->At(index);
  if (fDim>2) cout << " " << fNames[2] << "=" << fZ->At(index);
  if (fDim>3) cout << " " << fNames[3] << "=" << fT->At(index);
  cout << endl;
 }
}

void NcSample::ListOrdered(Int_t mode,TString name)
{
 
 Int_t i=GetIndex(name);
 ListOrdered(mode,i);
}

void NcSample::List(Int_t i)
{
 
 if (i<1 || i>fDim)
 {
  cout << " *NcSample::List(i)* Error : Invalid index " << i << endl;
 }
 else
 {
  cout << " Variable " << fNames[i-1] << " : N=" << fN;
  cout << " Sum=" << fSum[i-1] << " Mean=" << fMean[i-1];
  cout << " Deviation(rms)=" << fRMSdev[i-1] << " Sigma=" << fSigma[i-1];
  if (!fRemove || fStore)
  {
   cout << endl;
   cout << "          Minimum=" << GetMinimum(i) << " Maximum=" << GetMaximum(i);
  }
  if (fStore)
  {
   cout << " Median=" << GetMedian(i);
   cout << " Spread(w.r.t. median)=" << GetSpread(i,0) << " Spread(w.r.t. mean)=" << GetSpread(i,1);
  }
  cout << endl;
 }
}

void NcSample::List(Int_t i,Int_t j) const
{
 
 if (i<1 || i>fDim || j<1 || j>fDim)
 {
  cout << " *NcSample::List(i,j)* Error : Invalid index encountered i=" << i << " j=" << j << endl;
 }
 else
 {
  cout << " " << fNames[i-1] << "-" << fNames[j-1];
  cout << " Correlation(coefficient)=" << fCor[i-1][j-1] << " Covariance=" << fCov[i-1][j-1] << endl;
 }
}

void NcSample::SetStoreMode(Int_t mode,Int_t nmax,Int_t i)
{
 
 if (fN)
 {
  cout << "  Storage mode can only be set before first data." << endl;
 }
 else if ((!mode && nmax) || mode<0 || mode>1 || nmax<0 || abs(i)>fDim)
 {
  cout << "  Inconsistent input : mode=" << mode << " nmax=" << nmax << " i=" << i << endl;
  fNmax=0;
  fMoveVariable=0;
 }
 else
 {
  fStore=mode;
  fNmax=nmax;
  fMoveVariable=i;
 }
}

Int_t NcSample::GetStoreMode() const
{
 
 return fStore;
}

Int_t NcSample::GetStoreNmax() const
{
 
 return fNmax;
}

Double_t NcSample::GetQuantile(Double_t f,Int_t i)
{
 
 if (i<1 || i>fDim)
 {
  cout << " *NcSample::GetQuantile* Error : Invalid index " << i << endl;
  return 0;
 }
 
 if (!fStore)
 {
  cout << " *NcSample::GetQuantile* Error : Storage of data entries was not activated." << endl;
  return 0;
 }
 
 if (fN<=0) return 0;
 
 if (f<0 || f>1) return 0;
 
 Double_t quantile=0;
 
 if (fN==1)
 {
  if (i==1) quantile=fX->At(0);
  if (i==2) quantile=fY->At(0);
  if (i==3) quantile=fZ->At(0);
  if (i==4) quantile=fT->At(0);
  return quantile;
 }
 
 if (f==0) return GetMinimum(i);
 if (f==1) return GetMaximum(i);
 
 // Order the values of the i-th variable in increasing order
 Order(1,i);
 
 quantile=0;
 Double_t rindex=double(fN)*f;
 Int_t index=int(rindex);
 if (fN%2) // Odd number of entries
 {
  quantile=fArr->At(index);
 }
 else // Even number of entries
 {
  quantile=(fArr->At(index-1)+fArr->At(index))/2.;
 }
 return quantile;
}

Double_t NcSample::GetQuantile(Double_t f,TString name)
{
 
 Int_t i=GetIndex(name);
 return GetQuantile(f,i);
}

Double_t NcSample::GetMedian(Int_t i)
{
 
 Double_t median=GetQuantile(0.5,i);
 return median;
}

Double_t NcSample::GetMedian(TString name)
{
 
 Double_t median=GetQuantile(0.5,name);
 return median;
}

Double_t NcSample::GetSpread(Int_t i,Int_t model,Double_t vref)
{
 
 if (model<0 || model>2)
 {
  cout << " *NcSample::GetSpread* Error : Unsupported parameter model=" << model << endl;
  return -1;
 }
 
 if (i<1 || i>fDim)
 {
  cout << " *NcSample::GetSpread* Error : Invalid index " << i << endl;
  return -1;
 }
 
 if (!fStore)
 {
  cout << " *NcSample::GetSpread* Error : Storage of data entries was not activated." << endl;
  return -1;
 }
 
 if (fN<=1) return -1;
 
 Double_t central=GetMedian(i);
 if (model==1) central=GetMean(i);
 if (model==2) central=vref;
 
 Double_t spread=0;
 for (Int_t j=0; j<fN; j++)
 {
  spread+=fabs(central-(fArr->At(j)));
 }
 
 spread=spread/float(fN);
 
 return spread;
}

Double_t NcSample::GetSpread(TString name,Int_t model,Double_t vref)
{
 
 Int_t i=GetIndex(name);
 return GetSpread(i,model,vref);
}

Double_t NcSample::GetMinimum(Int_t i) const
{
 
 if (i<1 || i>fDim)
 {
  cout << " *NcSample::GetMinimum* Error : Invalid index " << i << endl;
  return 0;
 }
 
 if (!fRemove) return fMin[i-1];
 
 if (!fStore)
 {
  cout << " *NcSample::GetMinimum* Error : Storage of data entries was not activated." << endl;
  return 0;
 }
 
 if (fN<=0) return 0;
 
 Double_t min=0;
 
 if (i==1) min=fX->At(0);
 if (i==2) min=fY->At(0);
 if (i==3) min=fZ->At(0);
 
 for (Int_t k=1; k<fN; k++)
 {
  if (i==1 && fX->At(k)<min) min=fX->At(k);
  if (i==2 && fY->At(k)<min) min=fY->At(k);
  if (i==3 && fZ->At(k)<min) min=fZ->At(k);
 }
 
 return min;
}

Double_t NcSample::GetMinimum(TString name) const
{
 
 Int_t i=GetIndex(name);
 return GetMinimum(i);
}

Double_t NcSample::GetMaximum(Int_t i) const
{
 
 if (i<1 || i>fDim)
 {
  cout << " *NcSample::GetMaximum* Error : Invalid index " << i << endl;
  return 0;
 }
 
 if (!fRemove) return fMax[i-1];
 
 if (!fStore)
 {
  cout << " *NcSample::GetMaximum* Error : Storage of data entries was not activated." << endl;
  return 0;
 }
 
 if (fN<=0) return 0;
 
 Double_t max=0;
 
 if (i==1) max=fX->At(0);
 if (i==2) max=fY->At(0);
 if (i==3) max=fZ->At(0);
 
 for (Int_t k=1; k<fN; k++)
 {
  if (i==1 && fX->At(k)>max) max=fX->At(k);
  if (i==2 && fY->At(k)>max) max=fY->At(k);
  if (i==3 && fZ->At(k)>max) max=fZ->At(k);
 }
 
 return max;
}

Double_t NcSample::GetMaximum(TString name) const
{
 
 Int_t i=GetIndex(name);
 return GetMaximum(i);
}

Double_t NcSample::GetQuantile(Double_t f,TH1* histo,Int_t mode) const
{
 
 if (!histo) return 0;
 
 if (f<0 || f>1) return 0;
 
 Double_t xlow=0;
 Double_t xup=0;
 Double_t quantile=0;
 
 if (mode==2) // Quantile of Y values
 {
  NcSample temp;
  temp.SetStoreMode(1);
  Double_t val=0;
  for (Int_t i=1; i<=histo->GetNbinsX(); i++)
  {
   val=histo->GetBinContent(i);
   temp.Enter(val);
  }
  quantile=temp.GetQuantile(f,1);
 }
 else // Quantile of X values
 {
  // Take the average of two quantiles closely around f
  // This will enhance the precision for low statistics
  Double_t q[2];
  Double_t p[2];
  p[0]=f-0.01;
  p[1]=f+0.01;
  if (p[0]<0) p[0]=0;
  if (p[0]>1) p[0]=1;
  histo->ComputeIntegral();
  Int_t nq=histo->GetQuantiles(2,q,p);
 
  if (!nq) return 0;
 
  xlow=q[0];
  xup=q[1];
  if (mode==1)
  {
   Int_t mbin=histo->FindBin(q[0]);
   xlow=histo->GetBinCenter(mbin);
   mbin=histo->FindBin(q[1]);
   xup=histo->GetBinCenter(mbin);
  }
 }
 
 quantile=(xlow+xup)/2.;
 
 return quantile;
}

Double_t NcSample::GetMedian(TH1* histo,Int_t mode) const
{
 
 Double_t median=GetQuantile(0.5,histo,mode);
 return median;
}

Double_t NcSample::GetSpread(TH1* histo,Int_t mode,Int_t model,Double_t vref) const
{
 
 if (model<0 || model>2)
 {
  cout << " *NcSample::GetSpread* Error : Unsupported parameter model=" << model << endl;
  return -1;
 }
 
 if (!histo) return -1;
 
 Int_t nbins=histo->GetNbinsX();
 
 if (nbins<1) return -1;
 
 Double_t spread=0;
 
 if (mode==2) // Spread in Y values
 {
  NcSample temp;
  temp.SetStoreMode(1);
  Double_t val=0;
  for (Int_t i=1; i<=nbins; i++)
  {
   val=histo->GetBinContent(i);
   temp.Enter(val);
  }
  spread=temp.GetSpread(1,model,vref);
 }
 else // Spread in X values
 {
  Double_t central=GetMedian(histo,mode);
  if (model==1) central=histo->GetMean();
  if (model==2) central=vref;
  Double_t x=0;
  Double_t y=0;
  Double_t ysum=0;
  for (Int_t jbin=1; jbin<=nbins; jbin++)
  {
   x=histo->GetBinCenter(jbin);
   y=histo->GetBinContent(jbin);
   if (y>0)
   {
    spread+=fabs(x-central)*y;
    ysum+=y;
   }
  }
  if (ysum>0) spread=spread/ysum;
 }
 
 return spread;
}

Double_t NcSample::GetEntry(Int_t i,Int_t j,Int_t mode,Int_t k)
{
 
 TString name=GetName();
 
 if (!fStore)
 {
  cout << " *NcSample::GetEntry* Error : Storage mode not activated." << endl;
  return 0;
 }
 
 if (i<1 || i>fN)
 {
  cout << " *NcSample::GetEntry* Error : Invalid index number i=" << i;
  if (name!="") printf(" for NcSample %-s",name.Data());
  printf("\n");
  return 0;
 }
 
 if (j<1 || j>fDim)
 {
  cout << " *NcSample::GetEntry* Error : Invalid variable number j=" << j;
  if (name!="") printf(" for NcSample %-s",name.Data());
  printf("\n");
  return 0;
 }
 
 if (mode && (k<1 || k>fDim))
 {
  cout << " *NcSample::GetEntry* Error : Invalid argument k=" << k;
  if (name!="") printf(" for NcSample %-s",name.Data());
  printf("\n");
  return 0;
 }
 
 Double_t value=0;
 
 // Determine the entry index in the main arrays.
 Int_t index=i-1;
 
 if (mode)
 {
  Order(mode,k);
  index=fIndices->At(i-1);
 }
 
 if (j==1) value=fX->At(index);
 if (j==2) value=fY->At(index);
 if (j==3) value=fZ->At(index);
 if (j==4) value=fT->At(index);
 
 return value;
}

Double_t NcSample::GetEntry(Int_t i,TString nameA,Int_t mode,TString nameB)
{
 
 Int_t j=GetIndex(nameA);
 Int_t k=GetIndex(nameB);
 return GetEntry(i,j,mode,k);
}

void NcSample::GetSubset(NcSample* s,Int_t ifirst,Int_t ilast,Int_t mode,Int_t k)
{
 
 if (!s)
 {
  cout << " *NcSample::GetSubset* Error : Pointer for output NcSample object is zero." << endl;
  return;
 }
 
 s->Reset();
 
 if (!fStore)
 {
  cout << " *NcSample::GetSubset* Error : Storage mode not activated." << endl;
  return;
 }
 
 if (ifirst<1 || ifirst>fN || ilast>fN || ilast<ifirst)
 {
  cout << " *NcSample::GetSubset* Error : Invalid input ifirst=" << ifirst << " ilast=" << ilast << endl;
  return;
 }
 
 if (mode && (k<1 || k>fDim))
 {
  cout << " *NcSample::GetSubset* Error : Invalid argument k=" << k << endl;
  return;
 }
 
 s->SetStoreMode();
 s->SetNames(fNames[0],fNames[1],fNames[2],fNames[3]);
 
 Double_t* values=new Double_t[fDim];
 for (Int_t i=ifirst; i<=ilast; i++)
 {
  for (Int_t j=1; j<=fDim; j++)
  {
   values[j-1]=GetEntry(i,j,mode,k);
  }
  if (fDim==1) s->Enter(values[0]);
  if (fDim==2) s->Enter(values[0],values[1]);
  if (fDim==3) s->Enter(values[0],values[1],values[2]);
  if (fDim==4) s->Enter(values[0],values[1],values[2],values[3]);
 }
 
 delete [] values;
}

NcSample NcSample::GetDtSample(Int_t i,Int_t nc,Int_t store,Int_t nmax,Int_t order)
{
 
 // Make "order" value to comply with a single variable sample
 if (order>0) order=1;
 if (order<0) order=-1; 
 
 NcSample sdt;
 sdt.SetStoreMode(store,nmax,order);
 sdt.SetName("DtSample");
 sdt.SetNames("Dt");
 
 if (i<1 || i>GetDimension()) return sdt;
 
 if (nc<1) return sdt;
 
 Int_t nen=GetN();
 if (nen<=nc) return sdt;
 
 TString varname=GetVariableName(i);
 TString str;
 str.Form("Interval sampling for variable %-s between %-i consecutive entries (nc=%-i)",varname.Data(),nc+1,nc);
 sdt.SetNameTitle("DtSample",str);
 
 // Fill the dt sample
 Double_t t1=0;
 Double_t t2=0;
 Double_t deltat=0;
 for (Int_t ien=1; ien<=nen; ien++)
 {
  if (ien+nc>nen) break;
 
  t1=GetEntry(ien,i,1,i);
  t2=GetEntry(ien+nc,i,1,i);
  deltat=fabs(t2-t1);
 
  sdt.Enter(deltat);
 }
 
 return sdt;
}

NcSample NcSample::GetDtSample(TString name,Int_t nc,Int_t store,Int_t nmax,Int_t order)
{
 
 Int_t i=GetIndex(name);
 NcSample sdt=GetDtSample(i,nc,store,nmax,order);
 
 return sdt;
}

TH1D NcSample::GetSamplingHistogram(Int_t i,TF1* f)
{
 
 TString s="Sampling Histogram for NcSample ";
 s+=GetName();
 s+=";Sampling number;";
 if (f)
 {
  TString sf=f->GetExpFormula("p");
  sf.ReplaceAll("x",fNames[i-1]);   
  s+=sf;
 }
 else
 {
  s+=fNames[i-1];
 }
 
 TH1D hist("",s,fN,0,fN);
 
 if (!fStore)
 {
  cout << " *NcSample::GetSamplingHistogram* Error : Storage mode has not been activated." << endl;
  return hist;
 }
 
 if (i<1 || i>fDim)
 {
  cout << " *NcSample::GetSamplingHistogram* Error : Invalid index encountered i=" << i << endl;
  return hist;
 }
 
 Double_t y=0;
 for (Int_t ip=1; ip<=fN; ip++)
 {
  y=GetEntry(ip,i);
  if (f) y=f->Eval(y);
  hist.SetBinContent(ip,y);
 }
 
 return hist;
}

TH1D NcSample::GetSamplingHistogram(Int_t i,TString f)
{
 
 TF1 func("func",f);
 
 return GetSamplingHistogram(i,&func);
}

TH1D NcSample::GetSamplingHistogram(TString nameA,TF1* f)
{
 
 Int_t i=GetIndex(nameA);
 return GetSamplingHistogram(i,f);
}

TH1D NcSample::GetSamplingHistogram(TString nameA,TString f)
{
 
 TF1 func("func",f);
 
 return GetSamplingHistogram(nameA,&func);
}

TH1D NcSample::Get1DHistogram(Int_t i,Int_t j,Bool_t sumw2,Int_t nbx,TF1* f)
{
 
 TString s="1D Histogram for NcSample ";
 s+=GetName();
 s+=";";
 s+=fNames[i-1];
 s+=";Counts";
 if (j>0)
 {
  s+=" weighted by : ";
  if (f)
  {
   TString sf=f->GetExpFormula("p");
   sf.ReplaceAll("x",fNames[j-1]);   
   s+=sf;
  }
  else
  {
   s+=fNames[j-1];
  }
 }
 
 Double_t xlow=GetMinimum(i);
 Double_t xup=GetMaximum(i);
 
 Double_t dx=0;
 if (nbx) dx=(xup-xlow)/double(nbx);
 
 // Extension to include the maximum value
 xup+=1e-6*fabs(dx);
 
 TH1D hist("",s.Data(),nbx,xlow,xup);
 hist.Sumw2(sumw2);
 
 if (!fStore)
 {
  cout << " *NcSample::Get1DHistogram* Error : Storage mode has not been activated." << endl;
  return hist;
 }
 
 if (i<1 || i>fDim || j>fDim)
 {
  cout << " *NcSample::Get1DHistogram* Error : Invalid index encountered i=" << i << " j=" << j << endl;
  return hist;
 }
 
 Double_t x=0;
 Double_t y=0;
 for (Int_t ip=1; ip<=fN; ip++)
 {
  x=GetEntry(ip,i);
  if (j>0)
  {
   y=GetEntry(ip,j);
   if (f) y=f->Eval(y);
   hist.Fill(x,y);
  }
  else
  {
   hist.Fill(x);
  }
 }
 
 return hist;
}

TH1D NcSample::Get1DHistogram(TString nameA,TString nameB,Bool_t sumw2,Int_t nbx,TF1* f)
{
 
 Int_t i=GetIndex(nameA);
 Int_t j=GetIndex(nameB);
 return Get1DHistogram(i,j,sumw2,nbx,f);
}

TH2D NcSample::Get2DHistogram(Int_t i,Int_t j,Int_t k,Bool_t sumw2,Int_t nbx,Int_t nby,TF1* f)
{
 
 TString s="2D Histogram for NcSample ";
 s+=GetName();
 s+=";";
 s+=fNames[i-1];
 s+=";";
 s+=fNames[j-1];
 s+=";Counts";
 if (k>0)
 {
  s+=" weighted by : ";
  if (f)
  {
   TString sf=f->GetExpFormula("p");
   sf.ReplaceAll("x",fNames[k-1]);
   s+=sf;
  }
  else
  {
   s+=fNames[k-1];
  }
 }
 
 Double_t xlow=GetMinimum(i);
 Double_t xup=GetMaximum(i);
 Double_t ylow=GetMinimum(j);
 Double_t yup=GetMaximum(j);
 
 Double_t dx=0;
 if (nbx) dx=(xup-xlow)/double(nbx);
 Double_t dy=0;
 if (nby) dy=(yup-ylow)/double(nby);
 
 // Extension to include the maximum values
 xup+=1e-6*fabs(dx);
 yup+=1e-6*fabs(dy);
 
 TH2D hist("",s.Data(),nbx,xlow,xup,nby,ylow,yup);
 hist.Sumw2(sumw2);
 hist.SetMarkerStyle(20);
 hist.SetMarkerSize(1);
 
 if (!fStore)
 {
  cout << " *NcSample::Get2DHistogram* Error : Storage mode has not been activated." << endl;
  return hist;
 }
 
 if (i<1 || i>fDim || j<1 || j>fDim || k>fDim)
 {
  cout << " *NcSample::Get2DHistogram* Error : Invalid index encountered i=" << i << " j=" << j << " k=" << k << endl;
  return hist;
 }
 
 Double_t x=0;
 Double_t y=0;
 Double_t z=0;
 for (Int_t ip=1; ip<=fN; ip++)
 {
  x=GetEntry(ip,i);
  y=GetEntry(ip,j);
  if (k>0)
  {
   z=GetEntry(ip,k);
   if (f) z=f->Eval(z);
   hist.Fill(x,y,z);
  }
  else
  {
   hist.Fill(x,y);
  }
 }
 
 return hist;
}

TH2D NcSample::Get2DHistogram(TString nameA,TString nameB,TString nameC,Bool_t sumw2,Int_t nbx,Int_t nby,TF1* f)
{
 
 Int_t i=GetIndex(nameA);
 Int_t j=GetIndex(nameB);
 Int_t k=GetIndex(nameC);
 return Get2DHistogram(i,j,k,sumw2,nbx,nby,f);
}

TH3D NcSample::Get3DHistogram(Int_t i,Int_t j,Int_t k,Int_t m,Bool_t sumw2,Int_t nbx,Int_t nby,Int_t nbz,TF1* f)
{
 
 TString s="3D Histogram for NcSample ";
 s+=GetName();
 if (m>0)
 {
  s+=" with ";
  if (f)
  {
   TString sf=f->GetExpFormula("p");
   sf.ReplaceAll("x",fNames[m-1]);
   s+=sf;
  }
  else
  {
   s+=fNames[m-1];
  }
  s+=" as weight";
 }
 s+=";";
 s+=fNames[i-1];
 s+=";";
 s+=fNames[j-1];
 s+=";";
 s+=fNames[k-1];
 
 Double_t xlow=GetMinimum(i);
 Double_t xup=GetMaximum(i);
 Double_t ylow=GetMinimum(j);
 Double_t yup=GetMaximum(j);
 Double_t zlow=GetMinimum(k);
 Double_t zup=GetMaximum(k);
 
 Double_t dx=0;
 if (nbx) dx=(xup-xlow)/double(nbx);
 Double_t dy=0;
 if (nby) dy=(yup-ylow)/double(nby);
 Double_t dz=0;
 if (nbz) dz=(zup-zlow)/double(nbz);
 
 // Extension to include the maximum values
 xup+=1e-6*fabs(dx);
 yup+=1e-6*fabs(dy);
 zup+=1e-6*fabs(dz);
 
 TH3D hist("",s.Data(),nbx,xlow,xup,nby,ylow,yup,nbz,zlow,zup);
 hist.Sumw2(sumw2);
 hist.SetMarkerStyle(20);
 hist.SetMarkerSize(1);
 
 if (!fStore)
 {
  cout << " *NcSample::Get3DHistogram* Error : Storage mode has not been activated." << endl;
  return hist;
 }
 
 if (i<1 || i>fDim || j<1 || j>fDim || k<1 || k>fDim || m>fDim)
 {
  cout << " *NcSample::Get2DHistogram* Error : Invalid index encountered i=" << i << " j=" << j << " k=" << k << " m=" << m << endl;
  return hist;
 }
 
 Double_t x=0;
 Double_t y=0;
 Double_t z=0;
 Double_t t=0;
 for (Int_t ip=1; ip<=fN; ip++)
 {
  x=GetEntry(ip,i);
  y=GetEntry(ip,j);
  z=GetEntry(ip,k);
  if (m>0)
  {
   t=GetEntry(ip,m);
   if (f) t=f->Eval(t);
   hist.Fill(x,y,z,t);
  }
  else
  {
   hist.Fill(x,y,z);
  }
 }
 
 return hist;
}

TH3D NcSample::Get3DHistogram(TString nameA,TString nameB,TString nameC,TString nameD,Bool_t sumw2,Int_t nbx,Int_t nby,Int_t nbz,TF1* f)
{
 
 Int_t i=GetIndex(nameA);
 Int_t j=GetIndex(nameB);
 Int_t k=GetIndex(nameC);
 Int_t m=GetIndex(nameD);
 return Get3DHistogram(i,j,k,m,sumw2,nbx,nby,nbz,f);
}

TGraph NcSample::GetGraph(Int_t i,TF1* f)
{
 
 TGraph gr;
 gr.SetName("NcSample");
 
 if (!fStore)
 {
  cout << " *NcSample::GetGraph* Error : Storage mode has not been activated." << endl;
  return gr;
 }
 
 if (i<1 || i>fDim)
 {
  cout << " *NcSample::GetGraph* Error : Invalid index encountered i=" << i << endl;
  return gr;
 }
 
 Double_t x=0;
 Double_t y=0;
 for (Int_t ip=0; ip<fN; ip++)
 {
  x=ip+1;
  y=GetEntry(ip+1,i);
  if (f) y=f->Eval(y);
  gr.SetPoint(ip,x,y);
 }
 
 TString s="TGraph for NcSample ";
 s+=GetName();
 s+=" : X-axis=Sampling number ";
 s+="  Y-axis=";
 TString sf=fNames[i-1];
 if (f)
 {
  sf=f->GetExpFormula("p");
  sf.ReplaceAll("x",fNames[i-1]);
 }
 s+=sf;
 s+=";Sampling number;";
 s+=sf;
 
 gr.SetTitle(s.Data());
 
 gr.SetMarkerStyle(20);
 gr.SetMarkerSize(1);
 gr.SetDrawOption("AP");
 
 return gr;
}

TGraph NcSample::GetGraph(Int_t i,TString f)
{
 
 TF1 func("func",f);
 
 return GetGraph(i,&func);
}

TGraph NcSample::GetGraph(TString nameA,TF1* f)
{
 
 Int_t i=GetIndex(nameA);
 return GetGraph(i,f);
}

TGraph NcSample::GetGraph(Int_t i,Int_t j,TF1* f)
{
 
 TGraph gr;
 gr.SetName("NcSample");
 
 if (!fStore)
 {
  cout << " *NcSample::GetGraph* Error : Storage mode has not been activated." << endl;
  return gr;
 }
 
 if (i<1 || i>fDim || j<1 || j>fDim)
 {
  cout << " *NcSample::GetGraph* Error : Invalid index encountered i=" << i << " j=" << j << endl;
  return gr;
 }
 
 Double_t x=0;
 Double_t y=0;
 for (Int_t ip=0; ip<fN; ip++)
 {
  x=GetEntry(ip+1,i);
  y=GetEntry(ip+1,j);
  if (f) y=f->Eval(y);
  gr.SetPoint(ip,x,y);
 }
 
 TString s="TGraph for NcSample ";
 s+=GetName();
 s+=" : X-axis=";
 s+=fNames[i-1];
 s+="  Y-axis=";
 TString sf=fNames[j-1];
 if (f)
 {
  sf=f->GetExpFormula("p");
  sf.ReplaceAll("x",fNames[j-1]);
 }
 s+=sf;
 s+=";";
 s+=fNames[i-1];
 s+=";";
 s+=sf;
 
 gr.SetTitle(s.Data());
 
 gr.SetMarkerStyle(20);
 gr.SetMarkerSize(1);
 gr.SetDrawOption("AP");
 
 return gr;
}

TGraph NcSample::GetGraph(Int_t i,Int_t j,TString f)
{
 
 TF1 func("func",f);
 
 return GetGraph(i,j,&func);
}

TGraph NcSample::GetGraph(TString nameA,TString nameB,TF1* f)
{
 
 Int_t i=GetIndex(nameA);
 Int_t j=GetIndex(nameB);
 return GetGraph(i,j,f);
}

TGraphErrors NcSample::GetGraphErrors(TGraph* gr,Int_t ix,Int_t iy,TF1* fx,TF1* fy)
{
 
 TGraphErrors gre;
 
 if (!gr)
 {
  cout << " *NcSample::GetGraphErrors* Error : No input Graph is specified." << endl;
  return gre;
 }
 
 TString name=gr->GetName();
 name+=".E";
 
 gre.SetName(name);
 
 if ((ix || iy) && !fStore)
 {
  cout << " *NcSample::GetGraphErrors* Error : Storage mode has not been activated." << endl;
  return gre;
 }
 
 if (ix<0 || ix>fDim || iy<0 || iy>fDim)
 {
  cout << " *NcSample::GetGraphErros* Error : Invalid index encountered ix=" << ix << " iy=" << iy << endl;
  return gre;
 }
 
 // The number of data points of the input Graph
 Int_t np=gr->GetN();
 
 if (!np)
 {
  cout << " *NcSample::GetGraphErrors* Error : Input Graph does not contain any data points." << endl;
  return gre;
 }
 
 if ((ix || iy) && np!=fN)
 {
  cout << " *NcSample::GetGraphErrors* Error : Entries don't match Nsample=" << fN << " Ngraph=" << np << endl;
  return gre;
 }
 
 Double_t x=0;
 Double_t y=0;
 Double_t ex=0;
 Double_t ey=0;
 for (Int_t ip=0; ip<np; ip++)
 {
  gr->GetPoint(ip,x,y);
  if (ix) ex=GetEntry(ip+1,ix);
  if (fx)
  {
   if (!ix)
   {
    ex=fx->Eval(x);
   }
   else
   {
    ex=fx->Eval(ex);
   }
  }
  if (iy) ey=GetEntry(ip+1,iy);
  if (fy)
  {
   if (!iy)
   {
    ey=fy->Eval(y);
   }
   else
   {
    ey=fy->Eval(ey);
   }
  }
  ex=fabs(ex);
  ey=fabs(ey);
  gre.SetPoint(ip,x,y);
  gre.SetPointError(ip,ex,ey);
 }
 
 TString sfx="0";
 if (ix || fx)
 {
  sfx="|";
  if (!fx)
  {
   sfx+=fNames[ix-1];
  }
  else
  {
   sfx+=fx->GetExpFormula("p");
   if (ix) sfx.ReplaceAll("x",fNames[ix-1]);
  }
  sfx+="|";
 }
 
 TString sfy="0";
 if (iy || fy)
 {
  sfy="|";
  if (!fy)
  {
   sfy+=fNames[iy-1];
  }
  else
  {
   sfy+=fy->GetExpFormula("p");
   if (iy)
   {
    sfy.ReplaceAll("x",fNames[iy-1]);
   }
   else
   {
    sfy.ReplaceAll("x","y");
   }
  }
  sfy+="|";
 }
 
 TString s="TGraphErrors for Graph : ";
 s+=gr->GetName();
 s+="   X-error=";
 s+=sfx;
 s+=" Y-error=";
 s+=sfy;
 
 // Copy the axes titles of the input Graph (if present)
 TAxis* axis=0;
 axis=gr->GetXaxis();
 if (axis)
 {
  s+=";";
  s+=axis->GetTitle();
 } 
 axis=gr->GetYaxis();
 if (axis)
 {
  s+=";";
  s+=axis->GetTitle();
 } 
 
 gre.SetTitle(s.Data());
 
 gre.SetMarkerStyle(20);
 gre.SetMarkerSize(0.5);
 gre.SetDrawOption("AP");
 
 return gre;
}

TGraphErrors NcSample::GetGraphErrors(TGraph* gr,TString nameA,TString nameB,TF1* fx,TF1* fy)
{
 
 Int_t ix=GetIndex(nameA);
 Int_t iy=GetIndex(nameB);
 
 return GetGraphErrors(gr,ix,iy,fx,fy);
}

TGraphTime* NcSample::GetGraph(Int_t i,Int_t j,Int_t mode,Int_t k,Bool_t smp)
{
 
 if (!fStore)
 {
  cout << " *NcSample::GetGraph* Error : Storage mode has not been activated." << endl;
  return 0;
 }
 
 if (fN<1)
 {
  cout << " *NcSample::GetGraph* Error : No data entries are stored." << endl;
  return 0;
 }
 
 if (i<1 || i>fDim || j<1 || j>fDim || (mode && (k<1 || k>fDim)))
 {
  cout << " *NcSample::GetGraph* Error : Invalid index encountered i=" << i << " j=" << j << " k=" << k << endl;
  return 0;
 }
 
 Double_t xlow=GetMinimum(i);
 Double_t xup=GetMaximum(i);
 Double_t ylow=GetMinimum(j);
 Double_t yup=GetMaximum(j);
 
 // Extensions to display the points well within the frame
 xlow-=0.1*fabs(xlow);
 xup+=0.1*fabs(xup);
 ylow-=0.1*fabs(ylow);
 yup+=0.1*fabs(yup);
 
 if (fGraphT)
 {
  delete fGraphT;
  fGraphT=0;
 }
 
 fGraphT=new TGraphTime(fN,xlow,ylow,xup,yup);
 fGraphT->SetName("NcSample");
 
 Double_t x=0;
 Double_t y=0;
 TMarker* m=0;
 Int_t jpstart=0;
 for (Int_t istep=0; istep<fN; istep++)
 {
  if (!smp) jpstart=istep;
  for (Int_t jp=jpstart; jp<=istep; jp++)
  {
   x=GetEntry(jp+1,i,mode,k);
   y=GetEntry(jp+1,j,mode,k);
   m=new TMarker(x,y,20);
   m->SetMarkerSize(1);
   fGraphT->Add(m,istep);
  }  
 }
 
 TString s="TGraphTime for NcSample ";
 s+=GetName();
 s+=";";
 s+=fNames[i-1];
 s+=";";
 s+=fNames[j-1];
 fGraphT->SetNameTitle("",s.Data());
 
 return fGraphT;
}

TGraphTime* NcSample::GetGraph(TString nameA,TString nameB,Int_t mode,TString nameC,Bool_t smp)
{
 
 Int_t i=GetIndex(nameA);
 Int_t j=GetIndex(nameB);
 Int_t k=GetIndex(nameC);
 return GetGraph(i,j,mode,k,smp);
}

TGraph2D NcSample::GetGraph(Int_t i,Int_t j,Int_t k,TF1* f)
{
 
 TGraph2D gr;
 gr.SetName("NcSample");
 
 if (!fStore)
 {
  cout << " *NcSample::GetGraph* Error : Storage mode has not been activated." << endl;
  return gr;
 }
 
 if (fN<1)
 {
  cout << " *NcSample::GetGraph* Error : No data entries are stored." << endl;
  return gr;
 }
 
 if (i<1 || i>fDim || j<1 || j>fDim || k<1 || k>fDim)
 {
  cout << " *NcSample::GetGraph* Error : Invalid index encountered i=" << i << " j=" << j << " k=" << k << endl;
  return gr;
 }
 
 Double_t x=0;
 Double_t y=0;
 Double_t z=0;
 for (Int_t ip=0; ip<fN; ip++)
 {
  x=GetEntry(ip+1,i);
  y=GetEntry(ip+1,j);
  z=GetEntry(ip+1,k);
  if (f) z=f->Eval(z);
  gr.SetPoint(ip,x,y,z);
 }
 
 TString s="TGraph2D for NcSample ";
 s+=GetName();
 s+=" : X-axis=";
 s+=fNames[i-1];
 s+="  Y-axis=";
 s+=fNames[j-1];
 s+="  Z-axis=";
 TString sf=fNames[k-1];
 if (f)
 {
  sf=f->GetExpFormula("p");
  sf.ReplaceAll("x",fNames[k-1]);
 }
 s+=sf;
 s+=";";
 s+=fNames[i-1];
 s+=";";
 s+=fNames[j-1];
 s+=";";
 s+=sf;
 
 gr.SetTitle(s.Data());
 
 gr.SetMarkerStyle(20);
 gr.SetMarkerSize(1);
 gr.SetDrawOption("P");
 
 return gr;
}

TGraph2D NcSample::GetGraph(Int_t i,Int_t j,Int_t k,TString f)
{
 
 TF1 func("func",f);
 
 return GetGraph(i,j,k,&func);
}

TGraph2D NcSample::GetGraph(TString nameA,TString nameB,TString nameC,TF1* f)
{
 
 Int_t i=GetIndex(nameA);
 Int_t j=GetIndex(nameB);
 Int_t k=GetIndex(nameC);
 return GetGraph(i,j,k,f);
}

TGraph2D NcSample::GetGraph(TString nameA,TString nameB,TString nameC,TString f)
{
 
 TF1 func("func",f);
 
 return GetGraph(nameA,nameB,nameC,&func);
}

TGraph2DErrors NcSample::GetGraph2DErrors(TGraph2D* gr,Int_t ix,Int_t iy,Int_t iz,TF1* fx,TF1* fy,TF1* fz)
{
 
 TGraph2DErrors gre;
 
 if (!gr)
 {
  cout << " *NcSample::GetGraph2DErrors* Error : No input Graph is specified." << endl;
  return gre;
 }
 
 TString name=gr->GetName();
 name+=".E";
 
 gre.SetName(name);
 
 if ((ix || iy || iz) && !fStore)
 {
  cout << " *NcSample::GetGraph2DErrors* Error : Storage mode has not been activated." << endl;
  return gre;
 }
 
 if (ix<0 || ix>fDim || iy<0 || iy>fDim || iz<0 || iz>fDim)
 {
  cout << " *NcSample::GetGraph2DErrors* Error : Invalid index encountered ix=" << ix << " iy=" << iy << " iz=" << iz << endl;
  return gre;
 }
 
 // The number of data points of the input Graph
 Int_t np=gr->GetN();
 
 if (!np)
 {
  cout << " *NcSample::GetGraph2DErrors* Error : Input Graph does not contain any data points." << endl;
  return gre;
 }
 
 if ((ix || iy || iz) && np!=fN)
 {
  cout << " *NcSample::GetGraph2DErrors* Error : Entries don't match Nsample=" << fN << " Ngraph=" << np << endl;
  return gre;
 }
 
 Double_t* xarr=gr->GetX();
 Double_t* yarr=gr->GetY();
 Double_t* zarr=gr->GetZ();
 Double_t x=0;
 Double_t y=0;
 Double_t z=0;
 Double_t ex=0;
 Double_t ey=0;
 Double_t ez=0;
 for (Int_t ip=0; ip<np; ip++)
 {
  x=xarr[ip];
  y=yarr[ip];
  z=zarr[ip];
  if (ix) ex=GetEntry(ip+1,ix);
  if (fx)
  {
   if (!ix)
   {
    ex=fx->Eval(x);
   }
   else
   {
    ex=fx->Eval(ex);
   }
  }
  if (iy) ey=GetEntry(ip+1,iy);
  if (fy)
  {
   if (!iy)
   {
    ey=fy->Eval(y);
   }
   else
   {
    ey=fy->Eval(ey);
   }
  }
  if (iz) ez=GetEntry(ip+1,iz);
  if (fz)
  {
   if (!iz)
   {
    ez=fz->Eval(z);
   }
   else
   {
    ez=fz->Eval(ez);
   }
  }
  ex=fabs(ex);
  ey=fabs(ey);
  ez=fabs(ez);
  gre.SetPoint(ip,x,y,z);
  gre.SetPointError(ip,ex,ey,ez);
 }
 
 TString sfx="0";
 if (ix || fx)
 {
  sfx="|";
  if (!fx)
  {
   sfx+=fNames[ix-1];
  }
  else
  {
   sfx+=fx->GetExpFormula("p");
   if (ix) sfx.ReplaceAll("x",fNames[ix-1]);
  }
  sfx+="|";
 }
 
 TString sfy="0";
 if (iy || fy)
 {
  sfy="|";
  if (!fy)
  {
   sfy+=fNames[iy-1];
  }
  else
  {
   sfy+=fy->GetExpFormula("p");
   if (iy)
   {
    sfy.ReplaceAll("x",fNames[iy-1]);
   }
   else
   {
    sfy.ReplaceAll("x","y");
   }
  }
  sfy+="|";
 }
 
 TString sfz="0";
 if (iz || fz)
 {
  sfz="|";
  if (!fz)
  {
   sfz+=fNames[iz-1];
  }
  else
  {
   sfz+=fz->GetExpFormula("p");
   if (iz)
   {
    sfz.ReplaceAll("x",fNames[iz-1]);
   }
   else
   {
    sfz.ReplaceAll("x","z");
   }
  }
  sfz+="|";
 }
 
 TString s="TGraphErrors for Graph : ";
 s+=gr->GetName();
 s+="   X-error=";
 s+=sfx;
 s+=" Y-error=";
 s+=sfy;
 s+=" Z-error=";
 s+=sfz;
 
 // Copy the axes titles of the input Graph (if present)
 TAxis* axis=0;
 axis=gr->GetXaxis();
 if (axis)
 {
  s+=";";
  s+=axis->GetTitle();
 } 
 axis=gr->GetYaxis();
 if (axis)
 {
  s+=";";
  s+=axis->GetTitle();
 } 
 axis=gr->GetZaxis();
 if (axis)
 {
  s+=";";
  s+=axis->GetTitle();
 } 
 
 gre.SetTitle(s.Data());
 
 gre.SetMarkerStyle(20);
 gre.SetMarkerSize(0.5);
 gre.SetDrawOption("err P");
 
 return gre;
}

TGraph2DErrors NcSample::GetGraph2DErrors(TGraph2D* gr,TString nameA,TString nameB,TString nameC,TF1* fx,TF1* fy,TF1* fz)
{
 
 Int_t ix=GetIndex(nameA);
 Int_t iy=GetIndex(nameB);
 Int_t iz=GetIndex(nameC);
 
 return GetGraph2DErrors(gr,ix,iy,iz,fx,fy,fz);
}

TGraphQQ NcSample::GetQQplot(Int_t i,Int_t j,TF1* f)
{
 
 TGraphQQ gr;
 gr.SetName("NcSample");
 
 if (!f && (j<1 || j>fDim)) return gr;
 
 if (!fStore)
 {
  cout << " *NcSample::GetQQplot* Error : Storage mode has not been activated." << endl;
  return gr;
 }
 
 if (fN<1)
 {
  cout << " *NcSample::GetQQplot* Error : No data entries are stored." << endl;
  return gr;
 }
 
 if (i<1 || i>fDim || (!f && (j<1 || j>fDim)))
 {
  cout << " *NcSample::GetQQplot* Error : Invalid index encountered i=" << i << " j=" << j << endl;
  return gr;
 }
 
 Double_t* arri=0;
 if (i==1) arri=fX->GetArray();
 if (i==2) arri=fY->GetArray();
 if (i==3) arri=fZ->GetArray();
 if (i==4) arri=fT->GetArray();
 
 if (!f)
 {
  Double_t* arrj=0;
  if (j==1) arrj=fX->GetArray();
  if (j==2) arrj=fY->GetArray();
  if (j==3) arrj=fZ->GetArray();
  if (j==4) arrj=fT->GetArray();
 
  if (arri && arrj) gr=TGraphQQ(fN,arri,fN,arrj);
 }
 else
 {
  if (arri) gr=TGraphQQ(fN,arri,f);
 }
 
 TString s="QQ-plot (TGraphQQ) for NcSample ";
 s+=GetName();
 if (f)
 {
  s+=" of ";
  s+=fNames[i-1];
  s+=" versus f(a)=";
  TString sf=f->GetExpFormula("p");
  sf.ReplaceAll("x","a");
  s+=sf;
 }
 else
 {
  s+=";";
  s+=fNames[i-1];
  s+=";";
  s+=fNames[j-1];
 }
 
 gr.SetTitle(s.Data());
 
 gr.SetMarkerStyle(20);
 gr.SetMarkerSize(1);
 
 return gr;
}

TGraphQQ NcSample::GetQQplot(TString nameA,TString nameB,TF1* f)
{
 
 Int_t i=GetIndex(nameA);
 Int_t j=GetIndex(nameB);
 return GetQQplot(i,j,f);
}

void NcSample::Load(TGraph* g,Int_t clr)
{
 
 if (!g)
 {
  printf(" *%-s::Load* Error : No TGraph specified ==> No action taken. \n",ClassName());
  return;
 }
 
 if (!clr && fDim!=2)
 {
  printf(" *%-s::Load* Error : Dimension %-i does not match a TGraph. \n",ClassName(),fDim);
  return;
 }
 
 if (clr) Reset();
 
 Int_t npoints=g->GetN();
 Double_t* x=g->GetX();
 Double_t* y=g->GetY();
 
 if (!x || !y || npoints<1) return;
 
 for (Int_t i=0; i<npoints; i++)
 {
  Enter(x[i],y[i]);
 }
}

void NcSample::Load(TGraph2D* g,Int_t clr)
{
 
 if (!g)
 {
  printf(" *%-s::Load* Error : No TGraph2D specified ==> No action taken. \n",ClassName());
  return;
 }
 
 if (!clr && fDim!=3)
 {
  printf(" *%-s::Load* Error : Dimension %-i does not match a TGraph2D. \n",ClassName(),fDim);
  return;
 }
 
 if (clr) Reset();
 
 Int_t npoints=g->GetN();
 Double_t* x=g->GetX();
 Double_t* y=g->GetY();
 Double_t* z=g->GetZ();
 
 if (!x || !y || !z || npoints<1) return;
 
 for (Int_t i=0; i<npoints; i++)
 {
  Enter(x[i],y[i],z[i]);
 }
}

void NcSample::Load(TH1* h,Int_t clr)
{
 
 if (!h)
 {
  printf(" *%-s::Load* Error : No histogram specified ==> No action taken. \n",ClassName());
  return;
 }
 
 Int_t hdim=h->GetDimension();
 
 if (hdim!=1)
 {
  printf(" *%-s::Load* Error : Histograms of dimension %-i are not supported. \n",ClassName(),hdim);
  return;
 }
 
 if (!clr && fDim!=1)
 {
  printf(" *%-s::Load* Error : Dimension %-i does not match a 1-dimensional histogram. \n",ClassName(),fDim);
  return;
 }
 
 if (clr) Reset();
 
 Int_t nbins=h->GetNbinsX();
 
 Double_t val=0;
 for (Int_t i=1; i<=nbins; i++)
 {
  val=h->GetBinContent(i);
  Enter(val);
 }
}

void NcSample::Load(TArray* arr,Int_t clr)
{
 
 if (!arr)
 {
  printf(" *%-s::Load* Error : No data array specified ==> No action taken. \n",ClassName());
  return;
 }
 
 if (!clr && fDim!=1)
 {
  printf(" *%-s::Load* Error : Dimension %-i does not match a 1-dimensional data array. \n",ClassName(),fDim);
  return;
 }
 
 if (clr) Reset();
 
 Int_t n=arr->GetSize();
 
 Double_t val=0;
 for (Int_t i=0; i<n; i++)
 {
  val=arr->GetAt(i);
  Enter(val);
 }
}

Double_t NcSample::GetSNR(Int_t i,Int_t mode,Bool_t dB) const
{
 
 if (i<1 || i>fDim || !mode || abs(mode)>4) return -9999;
 
 Double_t snr=-1;
 Double_t ymax=GetMaximum(i);
 Double_t ymin=GetMinimum(i);
 Double_t mean=GetMean(i);
 Double_t var=GetVar(i);
 Double_t rmsdev=sqrt(var);
 Double_t sigma=GetSigma(i,1);
 Double_t stdev=rmsdev;
 if (mode<0) stdev=sigma;
 
 if (stdev<=0) return -9999;

 // SNR of a transient signal //
 
 if (abs(mode)>2)
 {
  if (abs(mode)==3) // Data are power or RMS amplitudes
  {
   snr=(ymax-ymin)/stdev;
  }
  if (abs(mode)==4) // Data are bipolar amplitudes
  {
   snr=(ymax-ymin)/(2.*stdev);
  }
 
  if (dB && snr>0) snr=10.*log10(snr);
 
  return snr;
 }

 // SNR of a continuous signal //
 
 // Treat the values as observed amplitudes and work in dB scale
 Double_t psignal=mean*mean;
 Double_t pnoise=stdev*stdev;
 
 if (psignal>0 && pnoise>0)
 {
  snr=10.*(log10(psignal)-log10(pnoise));
  if (abs(mode)==1) snr=snr/2.; // Treat values as observed power
  if (!dB) // Convert to the straight ratio
  {
   snr=snr/10.;
   snr=pow(10,snr);
  }
 }
 
 return snr;
}

Double_t NcSample::GetSNR(TString name,Int_t mode,Bool_t dB) const
{
 
 Int_t i=GetIndex(name);
 return GetSNR(i,mode,dB);
}

Double_t NcSample::GetCV(Int_t i,Int_t mode) const
{
 
 if (i<1 || i>fDim || abs(mode)<1 || abs(mode)>2) return -1;
 
 Double_t vmax=GetMaximum(i);
 Double_t vmin=GetMinimum(i);
 Double_t value=GetMean(i);
 if (abs(mode)==2) value=(vmax-vmin)/2.;
 Double_t stdev=GetSigma(i,0);
 if (mode<0) stdev=GetSigma(i,1);
 
 Double_t cv=-1;
 
 if (value) cv=fabs(stdev/value);
 
 return cv;
}

Double_t NcSample::GetCV(TString name,Int_t mode) const
{
 
 Int_t i=GetIndex(name);
 return GetCV(i,mode);
}

void NcSample::Animation(Int_t i,Int_t j,Int_t mode,Int_t k,Int_t delay,TString opt)
{
 
 if (!fStore || fN<1 || i<1 || i>fDim || j<1 || j>fDim || (mode && (k<1 || k>fDim)) || delay<0)
 {
  cout << " *NcSample::Animation* Inconsistent input data." << endl;
  return;
 }
 
 if (fCanvas)
 {
  delete fCanvas;
  fCanvas=0;
 }
 
 if (fAnimObject)
 {
  delete fAnimObject;
  fAnimObject=0;
 }
 
 fCanvas=new TCanvas("fCanvas","Sampling animation");
 
 TGraph* gr=new TGraph();
 fAnimObject=gr;
 
 TString s="Sampling animation for NcSample ";
 s+=GetName();
 gr->SetTitle(s.Data());
 gr->SetMarkerStyle(20);
 gr->SetMarkerSize(1);
 
 Double_t x=0;
 Double_t y=0;
 for (Int_t ip=0; ip<fN; ip++)
 {
  x=GetEntry(ip+1,i,mode,k);
  y=GetEntry(ip+1,j,mode,k);
  gr->SetPoint(ip,x,y);
 
  if (!ip) gr->Draw(opt.Data());
 
  s="Variable ";
  s+=fNames[i-1];
  gr->GetXaxis()->SetTitle(s.Data());
  s="Variable ";
  s+=fNames[j-1];
  gr->GetYaxis()->SetTitle(s.Data());
  
  fCanvas->Modified();
  fCanvas->Update();
 
  gSystem->Sleep(delay);
 }
}

void NcSample::Animation(TString nameA,TString nameB,Int_t mode,TString nameC,Int_t delay,TString opt)
{
 
 Int_t i=GetIndex(nameA);
 Int_t j=GetIndex(nameB);
 Int_t k=GetIndex(nameC);
 Animation(i,j,mode,k,delay,opt);
}

void NcSample::Animation(Int_t i,Int_t j,Int_t k,Int_t mode,Int_t m,Int_t delay,TString opt)
{
 
 if (!fStore || fN<1 || i<1 || i>fDim || j<1 || j>fDim || k<1 || k>fDim || (mode && (m<1 || m>fDim)) || delay<0)
 {
  cout << " *NcSample::Animation* Inconsistent input data." << endl;
  return;
 }
 
 if (fCanvas)
 {
  delete fCanvas;
  fCanvas=0;
 }
 
 if (fAnimObject)
 {
  delete fAnimObject;
  fAnimObject=0;
 }
 
 fCanvas=new TCanvas("fCanvas","Sampling animation");
 
 TGraph2D* gr=new TGraph2D(fN);
 fAnimObject=gr;
 
 TString s="Sampling animation for NcSample ";
 s+=GetName();
 gr->SetTitle(s.Data());
 gr->SetMarkerStyle(20);
 gr->SetMarkerSize(1);
 
 Double_t x=0;
 Double_t y=0;
 Double_t z=0;
 for (Int_t ip=0; ip<fN; ip++)
 {
  x=GetEntry(ip+1,i,mode,m);
  y=GetEntry(ip+1,j,mode,m);
  z=GetEntry(ip+1,k,mode,m);
  gr->SetPoint(ip,x,y,z);
 
  if (!ip) gr->Draw(opt.Data());
 
  s="Variable ";
  s+=fNames[i-1];
  gr->GetXaxis()->SetTitle(s.Data());
  s="Variable ";
  s+=fNames[j-1];
  gr->GetYaxis()->SetTitle(s.Data());
  s="Variable ";
  s+=fNames[k-1];
  gr->GetZaxis()->SetTitle(s.Data());
  
  fCanvas->Modified();
  fCanvas->Update();
 
  gSystem->Sleep(delay);
 }
}

void NcSample::Animation(TString nameA,TString nameB,TString nameC,Int_t mode,TString nameD,Int_t delay,TString opt)
{
 
 Int_t i=GetIndex(nameA);
 Int_t j=GetIndex(nameB);
 Int_t k=GetIndex(nameC);
 Int_t m=GetIndex(nameD);
 Animation(i,j,k,mode,m,delay,opt);
}

Double_t NcSample::Digitize(Int_t i,Int_t nbits,Double_t vcal,Int_t mode)
{
 
 if (mode<0 || mode>3)
 {
  cout << " *NcSample::Digitize* Inconsistent input mode=" << mode << endl;
  return 0;
 }
 
 if (!nbits || nbits>60 || nbits<-60 || fabs(vcal)<=0)
 {
  cout << " *NcSample::Digitize* Inconsistent input nbits=" << nbits << " vcal=" << vcal << endl;
  return 0;
 }
 
 if ((mode==1 || mode==3) && nbits==1)
 {
  cout << " *NcSample::Digitize* Inconsistent input nbits=" << nbits << " mode=" << mode << endl;
  return 0;
 }
 
 Bool_t logmode=kFALSE;
 if (nbits<0) logmode=kTRUE;
 
 nbits=abs(nbits);
 
 Long_t nlevels=pow(2,nbits);
 Double_t range=fabs(vcal);
 if (mode==1 || mode==3) range*=2.;
 Double_t step=fabs(vcal);
 if (mode==0) step=range/double(nlevels-1);
 if (mode==1) step=range/double(nlevels-2);
 
 if (step<=0) return 0;
 
 Long_t nstepsmin=0;
 Long_t nstepsmax=nlevels-1;
 if ((mode==0 || mode==2) && vcal<0)
 {
  nstepsmin=-nlevels+1;
  nstepsmax=0;
 }
 if (mode==1 || mode==3)
 {
  nstepsmin=-nlevels/2;
  nstepsmax=nlevels/2-1;
 }
 
 Double_t digvalmin=double(nstepsmin)*step;
 Double_t digvalmax=double(nstepsmax)*step;
 
 cout << " *NcSample::Digitize* Parameters of the " << nbits << "-bits";
 if (logmode) cout << " Log10";
 cout << " ADC digitization";
 if (i>0 && i<=fDim) cout << " of variable : " << i << " (" << fNames[i-1] << ")"; 
 cout << endl;
 TString s="linear";
 if (logmode) s="Log10";
 if (mode==0 || mode==2)
 {
  cout << " Digitized " << s << " full scale range : [" << digvalmin << "," << digvalmax << "]  Step size : " << step << endl;
 }
 if (mode==1 || mode==3)
 {
  cout << " Digitized " << s << " full scale range : [" << (digvalmin+step) << "," << digvalmax << "]  Step size : " << step;
  cout << "  Actual range : [" << digvalmin << "," << digvalmax << "]" << endl;
 }
 
 if (logmode)
 {
  Double_t linvalmin=pow(10,digvalmin);
  Double_t linvalmax=pow(10,digvalmax);
  if (mode==0 || mode==2)
  {
   cout << " Digitized linear full scale range : [" << linvalmin << "," << linvalmax << "]" << endl;
  }
  if (mode==1 || mode==3)
  {
   cout << " Digitized linear full scale range : [" << (linvalmin*pow(10,step)) << "," << linvalmax << "]";
   cout << "  Actual range : [" << linvalmin << "," << linvalmax << "]" << endl;
  }
 }
 
 Double_t retval=step;
 if (mode==2)
 {
  if (vcal<0) retval=digvalmin;
  if (vcal>0) retval=digvalmax;
 }
 if (mode==3) retval=digvalmax;
 
 if (!fStore)
 {
  cout << endl;
  cout << " *NcSample::Digitize* Error : Storage mode has not been activated." << endl;
  return retval;
 }
 
 Double_t minval=1;
 if (i>0 && i<=fDim) minval=GetMinimum(i);
 
 if (logmode && minval<=0)
 {
  cout << endl;
  cout << " *NcSample::Digitize* Non-positive input value(s) for Log10 ADC i=" << i << " minimum=" << minval << endl;
  return retval;
 }
 
 if (i<1 || i>fDim)
 {
  cout << endl;
  cout << " *NcSample::Digitize* Value i=" << i << " is outside range. Only listing of ADC specs." << endl;
  return retval;
 }
 
 Double_t val=0;
 Long_t nsteps=0;
 Double_t digval=0;
 for (Int_t j=1; j<=fN; j++)
 {
  val=GetEntry(j,i);
 
  if (logmode) val=log10(val);
 
  nsteps=val/step;
 
  if (nsteps<nstepsmin) nsteps=nstepsmin;
  if (nsteps>nstepsmax) nsteps=nstepsmax;
 
  digval=double(nsteps)*step;
 
  if (logmode) digval=pow(10,digval);
 
  if (i==1) fX->AddAt(digval,j-1);
  if (i==2) fY->AddAt(digval,j-1);
  if (i==3) fZ->AddAt(digval,j-1);
  if (i==4) fT->AddAt(digval,j-1);
 }
 
 return retval;
}

Double_t NcSample::Digitize(TString name,Int_t nbits,Double_t vcal,Int_t mode)
{
 
 Int_t i=GetIndex(name);
 
 return Digitize(i,nbits,vcal,mode);
}

TArrayL64 NcSample::ADC(Int_t nbits,Double_t range,Double_t Vbias,TArray* Vsig,TH1* hist,Int_t B,Int_t C) const
{
 
 TArrayL64 arradc(0);
 if (hist) hist->Reset();
 
 Int_t nVsig=0;
 if (Vsig) nVsig=Vsig->GetSize();
 
 if (nbits<=0 || nbits>60 || !range || fabs(Vbias)>fabs(range) || B<0 || (B && C<1))
 {
  printf("\n *%-s::ADC* Inconsistent input nbits=%-i range=%-g Vbias=%-g B=%-i C=%-i \n",ClassName(),nbits,range,Vbias,B,C);
  return arradc;
 }
 
 NcMath math;
 
 Long64_t N=pow(2,nbits); // The number of quantization levels
 Long64_t adcmin=0;
 Long64_t adcmax=N-1;
 Double_t Vref=0;
 Double_t Vfs=0;
 Double_t LSB=0;
 
 // Floating point version of some parameters
 Double_t rN=N;
 Double_t rB=B;
 if (B==1) rB=exp(1);
 Double_t rC=C;
 Double_t radcmax=adcmax;
 
 if (range<0) // |range| represents Vref
 {
  Vref=fabs(range);
  LSB=Vref/rN;
  Vfs=Vref-LSB;
 }
 else // range represents Vfs
 {
  Vfs=range;
  LSB=Vfs/radcmax;
  Vref=Vfs+LSB;
 }
 
 Long64_t ped=0;
 Double_t frac=0;
 if (B) // Log ADC
 {
  if (range<0) Vfs=pow(rB,-rC)*pow(rB,radcmax*rC/rN)*Vref;
  if (range>0) Vref=pow(rB,rC)*pow(rB,-radcmax*rC/rN)*Vfs;
  LSB=Vref*pow(rB,-rC)*(pow(rB,rC/rN)-1.);
  frac=Vbias/Vref;
  ped=0;
  if (frac>0) ped=Long64_t(rN*(math.Log(rB,frac)+rC)/rC);
 }
 
 if (LSB<=0 || Vfs<=0)
 {
  printf("\n *%-s::ADC* Inconsistent parameters : LSB=%-g Vfs=%-g \n",ClassName(),LSB,Vfs);
  return arradc;
 }
 
 if (!B) ped=Vbias/LSB; // Pedestal value for a linear ADC
 
 Double_t DR=20.*log10(Vfs/LSB);
 
 TString mode="linear";
 if (B==1) mode="Log_e";
 if (B>1)
 {
  mode="Log_";
  mode+=B;
 }
 
 if (!nVsig)
 {
  printf(" *%-s::ADC* No input data have been provided --> Only the value of adc(Vbias) is returned. \n",ClassName());
  if (!B)
  {
   printf(" Parameters for a %-i-bits %-s ADC with adc=[%-lli,%-lli] : \n",nbits,mode.Data(),adcmin,adcmax);
   printf(" Vref=%-g Vfs=%-g LSB=%-g DR=%-g (dB) Vbias=%-g adc(Vbias)=%-lli \n",Vref,Vfs,LSB,DR,Vbias,ped);
  }
  else
  {
   printf(" Parameters for a %-i-bits %-s ADC with adc=[%-lli,%-lli] and a code efficiency factor of %-i: \n",nbits,mode.Data(),adcmin,adcmax,C);
   printf(" Vref=%-g Vfs=%-g LSB=%-g DR=%-g (dB) Vbias=%-g adc(Vbias)=%-lli \n",Vref,Vfs,LSB,DR,Vbias,ped);
  }
  arradc.Set(1);
  arradc[0]=ped;
  return arradc;
 }
 
 arradc.Set(nVsig);
 
 if (hist)
 {
  TString title;
  title.Form("%-s %-i-bits %-s ADC with Vfs=%-g LSB=%-g;Sample number;ADC counts",ClassName(),nbits,mode.Data(),Vfs,LSB);
  hist->SetBins(nVsig,0,nVsig);
  hist->SetMarkerStyle(20);
  hist->SetTitle(title);
 }
 
 Double_t val=0;
 Double_t radcval=0;
 Long64_t adcval=0;
 for (Int_t j=0; j<nVsig; j++)
 {
  val=Vbias+Vsig->GetAt(j);
  
  if (B) // Log ADC
  {
   radcval=0;
   frac=val/Vref;
   if (frac>0) radcval=(rN/rC)*(math.Log(rB,frac)+rC);
  }
  else // Linear ADC
  {
   radcval=val/LSB;
  }
 
  adcval=Long64_t(radcval);
 
  // Check for saturation
  if (adcval<adcmin) adcval=adcmin;
  if (adcval>adcmax) adcval=adcmax;
 
  arradc[j]=adcval;
 
  if (hist) hist->SetBinContent(j+1,adcval);
 }
 
 return arradc;
}

TArrayD NcSample::DAC(Int_t nbits,Double_t range,Double_t Vbias,TArray* adcs,TArray* peds,TH1* hist,Int_t B,Int_t C) const
{
 
 TArrayD arrdac(0);
 if (hist) hist->Reset();
 
 Int_t nadcs=0;
 if (adcs) nadcs=adcs->GetSize();
 
 Int_t npeds=0;
 if (peds) npeds=peds->GetSize();
 
 if (nbits<=0 || nbits>60 || !range || fabs(Vbias)>fabs(range) || (peds && npeds<nadcs) || B<0 || (B && C<1))
 {
  printf("\n *%-s::DAC* Inconsistent input nbits=%-i range=%-g Vbias=%-g nadcs=%-i npeds=%-i B=%-i C=%-i \n",ClassName(),nbits,range,Vbias,nadcs,npeds,B,C);
  return arrdac;
 }
 
 NcMath math;
 
 Long64_t N=pow(2,nbits); // The number of quantization levels
 Long64_t adcmin=0;
 Long64_t adcmax=N-1;
 Double_t Vref=0;
 Double_t Vfs=0;
 Double_t LSB=0;
 
 // Floating point version of some parameters
 Double_t rN=N;
 Double_t rB=B;
 if (B==1) rB=exp(1);
 Double_t rC=C;
 Double_t radcmax=adcmax;
 
 if (range<0) // |range| represents Vref
 {
  Vref=fabs(range);
  LSB=Vref/rN;
  Vfs=Vref-LSB;
 }
 else // range represents Vfs
 {
  Vfs=range;
  LSB=Vfs/radcmax;
  Vref=Vfs+LSB;
 }
 
 Long64_t ped=0;
 Double_t frac=0;
 if (B) // Log DAC
 {
  if (range<0) Vfs=pow(rB,-rC)*pow(rB,radcmax*rC/rN)*Vref;
  if (range>0) Vref=pow(rB,rC)*pow(rB,-radcmax*rC/rN)*Vfs;
  LSB=Vref*pow(rB,-rC)*(pow(rB,rC/rN)-1.);
  frac=Vbias/Vref;
  ped=0;
  if (frac>0) ped=Long64_t(rN*(math.Log(rB,frac)+rC)/rC);
 }
 
 if (LSB<=0 || Vfs<=0)
 {
  printf("\n *%-s::DAC* Inconsistent parameters : LSB=%-g Vfs=%-g \n",ClassName(),LSB,Vfs);
  return arrdac;
 }
 
 if (!B) ped=Vbias/LSB; // Pedestal value for a linear DAC
 
 Double_t DR=20.*log10(Vfs/LSB);
 
 TString mode="linear";
 if (B==1) mode="Log_e";
 if (B>1)
 {
  mode="Log_";
  mode+=B;
 }
 
 if (!nadcs)
 {
  printf(" *%-s::DAC* No input data have been provided --> Only the value of adc(Vbias) is returned. \n",ClassName());
  if (!B)
  {
   printf(" Parameters for a %-i-bits %-s DAC with adc=[%-lli,%-lli] : \n",nbits,mode.Data(),adcmin,adcmax);
   printf(" Vref=%-g Vfs=%-g LSB=%-g DR=%-g (dB) Vbias=%-g adc(Vbias)=%-lli \n",Vref,Vfs,LSB,DR,Vbias,ped);
  }
  else
  {
   printf(" Parameters for a %-i-bits %-s DAC with adc=[%-lli,%-lli] and a code efficiency factor of %-i: \n",nbits,mode.Data(),adcmin,adcmax,C);
   printf(" Vref=%-g Vfs=%-g LSB=%-g DR=%-g (dB) Vbias=%-g adc(Vbias)=%-lli \n",Vref,Vfs,LSB,DR,Vbias,ped);
  }
  arrdac.Set(1);
  arrdac[0]=ped;
  return arrdac;
 }
 
 arrdac.Set(nadcs);
 
 if (hist)
 {
  TString title;
  title.Form("%-s %-i-bits %-s DAC with Vfs=%-g LSB=%-g;Sample number;Amplitude",ClassName(),nbits,mode.Data(),Vfs,LSB);
  hist->SetBins(nadcs,0,nadcs);
  hist->SetMarkerStyle(20);
  hist->SetTitle(title);
 }
 
 Long64_t adc=0;
 Double_t radc=0;
 Double_t dacval=0;
 for (Int_t j=0; j<nadcs; j++)
 {
  adc=adcs->GetAt(j);
  radc=adc;
 
  if (B) // Log DAC
  {
   dacval=Vref*pow(rB,-rC)*pow(rB,radc*rC/rN);
  }
  else // Linear DAC
  {
   if (peds)
   {
    ped=peds->GetAt(j);
    adc=adc-ped;
    radc=adc;
   }
   dacval=radc*LSB;
  }
 
  // Correct for bias voltage
  if (B || (!B && !peds))
  {
   dacval=dacval-Vbias;
  }
 
  arrdac[j]=dacval;
 
  if (hist) hist->SetBinContent(j+1,dacval);
 }
 
 return arrdac;
}

Long64_t NcSample::Transmit(Int_t i,Int_t nbits,Double_t range,Double_t Vbias,TArray* peds,TH1* hist,Int_t B,Int_t C)
{
 
 TArrayL64 adcarr;
 TArrayD dacarr;
 if (hist) hist->Reset();
 
 if (!fStore)
 {
  printf("\n *%-s::Transmit* Error : Storage mode has not been activated. \n",ClassName());
  printf(" No digitization is performed. \n");
  return -1;
 }
 
 if (i<1 || i>fDim)
 {
  printf("\n *%-s::Transmit* Invalid variable name or index (i=%-i) --> Only listing of ADC specs. \n",ClassName(),i);
  adcarr=ADC(nbits,range,Vbias,0,0,B,C);
  return adcarr[0];
 }
 
 if (!fN)
 {
  printf("\n *%-s::Transmit* No data available for variable %-s (index=%-i) --> Only listing of ADC specs. \n",ClassName(),fNames[i-1].Data(),i);
  adcarr=ADC(nbits,range,Vbias,0,0,B,C);
  return adcarr[0];
 }
 
 Double_t minval=GetMinimum(i);
 Double_t Vin=minval+Vbias;
 
 if (B && Vin<=0)
 {
  TString mode="Log_e";
  if (B>1)
  {
   mode="Log_";
   mode+=B;
  }
  printf("\n *%-s::Transmit* Non-positive input value(s) for %-s ADC --> No digitization is performed. \n",ClassName(),mode.Data());
  printf(" Minimum+Vbias=%-g+%-g=%-g for variable with name %-s (index i=%-i). \n",minval,Vbias,Vin,fNames[i-1].Data(),i);
  return -1;
 }
 
 // Get the corresponding data array
 TArrayD* arrin=0;
 if (i==1) arrin=fX;
 if (i==2) arrin=fY;
 if (i==3) arrin=fZ;
 if (i==4) arrin=fT;
 
 // Perform the digitization via the ADC processeor
 adcarr=ADC(nbits,range,Vbias,arrin,0,B,C);
 
 // Convert the digital data into analog signals via the DAC processor
 dacarr=DAC(nbits,range,Vbias,&adcarr,peds,hist,B,C);
 
 if (hist)
 {
  TString title=hist->GetTitle();
  title.ReplaceAll("DAC","ADC-DAC (Transmit)");
  hist->SetTitle(title);
 }
 
 Int_t ndac=dacarr.GetSize();
 if (ndac<fN)
 {
  printf("\n *%-s::Transmit* Inconsistent number of entries : Nsamplings=%-i Ndac=%-i --> No digitization is performed. \n",ClassName(),fN,ndac);
  return -1;
 }
 
 // Replace the values of the i-th variable
 Double_t val=0;
 for (Int_t j=0; j<fN; j++)
 {
  val=dacarr[j];
  arrin->AddAt(val,j);
 }
 
 return -1;
}

Long64_t NcSample::Transmit(TString name,Int_t nbits,Double_t range,Double_t Vbias,TArray* peds,TH1* hist,Int_t B,Int_t C)
{
 
 Int_t i=GetIndex(name);
 
 return Transmit(i,nbits,range,Vbias,peds,hist,B,C);
}

NcSample NcSample::SampleAndHold(TF1 f,Double_t step,Double_t vmin,Double_t vmax,Int_t loc) const
{
 
 NcSample s("SampleAndHold","Error occurred in SampleAndHold");
 s.SetStoreMode(1);
 
 if (step<=0 || vmax<=vmin)
 {
  cout << " *NcSample::SampleAndHold* Error : Invalid input step=" << step << " vmin=" << vmin << " vmax=" << vmax << endl;
  return s;
 }
 
 TString str="For Function ";
 str+=f.GetExpFormula("p");
 str+=" with step=%-10.3g";
 TString str2=str.Format(str.Data(),step);
 s.SetTitle(str2.Data());
 
 // Enter the sampled data into the new sample
 Double_t xlow=vmin;
 Double_t xval=0;
 Double_t yval=0;
 while (xlow<=vmax)
 {
  yval=f.Eval(xlow);
  if (loc<0) xval=xlow;
  if (!loc) xval=xlow+step/2.;
  if (loc>0) xval=xlow+step;
  if (xval>vmax) xval=vmax;
  s.Enter(xval,yval);
 
  xlow+=step;
 }
 
 return s;
}

NcSample NcSample::SampleAndSum(Int_t i,Double_t step,Int_t loc,Int_t j,Double_t vmin,Double_t vmax)
{
 
 NcSample s("SampleAndSum","Error occurred in SampleAndSum");
 s.SetStoreMode(1);
 
 if (!fStore)
 {
  cout << " *NcSample::SampleAndSum* Error : Storage mode has not been activated." << endl;
  return s;
 }
 
 if (i<1 || i>fDim || j>fDim)
 {
  cout << " *NcSample::SampleAndSum* Error : Invalid index encountered i=" << i << " j=" << j << endl;
  return s;
 }
 
 if (step<=0)
 {
  cout << " *NcSample::SampleAndSum* Error : Invalid step size step=" << step << endl;
  return s;
 }
 
 TString xname=fNames[i-1];
 TString yname="Counts";
 if (j>0) yname=fNames[j-1];
 
 s.SetNames(xname,yname);
 
 TString str="For Variable ";
 str+=xname;
 str+=" with step=%-10.3g";
 TString str2=str.Format(str.Data(),step);
 s.SetTitle(str2.Data());
 
 // Define a histogram with "step" as binsize
 Double_t xmin=GetMinimum(i);
 Double_t xmax=GetMaximum(i);
 
 if (vmax>vmin)
 {
  xmin=vmin;
  xmax=vmax;
 }
 
 Int_t nbx=(xmax-xmin)/step;
 nbx++;
 
 Double_t* bins=new Double_t[nbx+1];
 Double_t xlow=xmin;
 for (Int_t idx=0; idx<=nbx; idx++)
 {
  bins[idx]=xlow;
  xlow+=step;
 }
 
 TH1D hist("","",nbx,bins);
 
 // Fill the histogram with data of the original sample
 Double_t xval=0;
 Double_t yval=1;
 for (Int_t ient=1; ient<=fN; ient++)
 {
  xval=GetEntry(ient,i);
  if (j>0) yval=GetEntry(ient,j);
  hist.Fill(xval,yval);
 }
 
 // Enter the histgram data into the new sample
 Int_t nbins=hist.GetNbinsX();
 Double_t binwidth=hist.GetBinWidth(1);
 for (Int_t k=1; k<=nbins; k++)
 {
  if (loc) xval=hist.GetBinLowEdge(k);
  if (!loc) xval=hist.GetBinCenter(k);
  if (loc>0) xval+=binwidth;
  yval=hist.GetBinContent(k);
  s.Enter(xval,yval);
 }
 
 delete [] bins;
 return s;
}

NcSample NcSample::SampleAndSum(TString nameA,Double_t step,Int_t loc,TString nameB,Double_t vmin,Double_t vmax)
{
 
 Int_t i=GetIndex(nameA);
 Int_t j=GetIndex(nameB);
 
 return SampleAndSum(i,step,loc,j,vmin,vmax);
}

TObject* NcSample::Clone(const char* name) const
{
 
 NcSample* q=new NcSample(*this);
 if (name)
 {
  if (strlen(name)) q->SetName(name);
 }
 return q;
}