NcBlocks.cxx

Go to the documentation of this file.

 
#include "NcBlocks.h"
#include "Riostream.h"
 
ClassImp(NcBlocks); // Class implementation to enable ROOT I/O
 
NcBlocks::NcBlocks()
{
 
 fMode=0;
}

NcBlocks::~NcBlocks()
{
}

NcBlocks::NcBlocks(const NcBlocks& q)
{
 
 fMode=0;
}

Double_t NcBlocks::GetPrior(Int_t n,Double_t fpr)
{
 
 Double_t prior=0;
 
 if (fMode<1 || fMode>3 || n<1 || fpr<0 || fpr>1)
 {
  cout << " *NcBlocks::GetPrior* Inconsistent input : mode=" << fMode << " n=" << n << " fpr=" << fpr << endl;
  return 0;
 }
 
 Double_t rn=n;
 
 if (fMode==3)
 {
  Double_t c=0;
  Double_t s=0;
  if (fpr>0.045 && fpr<0.055) // Fit of J.D. Scargle for fpr=0.05
  {
   c=1.32;
   s=0.577;
  }
  else
  {
   c=51.29*TMath::Landau(fpr,-0.152,0.03167);
   s=0.5807+0.2317*fpr;
  }
  prior=c+s*log10(rn);
 }
 else // Data Mode 1 and 2
 {
  prior=4.-log(73.53*fpr*pow(rn,-0.478));
 }
 
 prior=-prior;
 
 return prior;
}

Double_t NcBlocks::GetBlockFitness(Double_t n,Double_t len)
{
 
 Double_t fb=0;
 
 if ((fMode!=1 && fMode!=2) || (fMode==1 && n<1) || len<=0)
 {
  cout << " *NcBlocks::GetBlockFitness* Inconsistent input : mode=" << fMode << " n=" << n << " len=" << len << endl;
  return 0;
 }
 
 if (n<=0) return 0;
 
 fb=n*log(n/len);
 
 return fb;
}

Double_t NcBlocks::GetBlocks(TH1* hin,Double_t fpr,TH1* hout,Int_t ntrig)
{
 
 if (!hin || !hout) return 0;
 
 Int_t n=hin->GetNbinsX();
 
 if (!n) return 0;
 
 // Set the Data Mode for binned data, if it was not set already
 // Also initialize the output histogram title for binned data
 if (!fMode)
 {
  fMode=2;
  TString title;
  title.Form("Bayesian Block representation with FPR=%-.3g",fpr);
  hout->SetTitle(title);
 }
 
 TArrayD best(n);    // Array with optimal fitness values
 TArrayI last(n);    // Array with indices of optimal last change points
 TArrayD lengths(n); // Array with the lengths of the optimal partition blocks
 TArrayD counts(n);  // Array with the event counts of the optimal partition blocks
 
 Double_t blen=0;   // Length of a certain block
 Double_t bcount=0; // Event count in a certain block
 Double_t prior=0; // Contribution of each block to the prior fitness for a certain partition
 Double_t bfit=0;  // Fitness of a certain block
 Double_t pfit=0;  // Fitness of a certain partition
 
 Double_t xlow=0;
 Double_t xup=0;
 Int_t index=0;
 
 // Variables for measurement of an observable (i.e. Data Mode 3)
 Double_t yk=0;   // Measured value
 Double_t sigk=0; // Error on the measured value
 Double_t a=0; 
 Double_t b=0;
 
 // The parameters for the optimal partition
 Double_t optfit=0;
 Int_t optj=0;
 Double_t optlen=0;
 Double_t optcount=0;
 
 // Parameters and counter for trigger mode
 Int_t oldoptj=0;
 Double_t oldoptlen=0;
 Int_t ncp=0;
 Double_t xtrig=0;
 Double_t ytrig=0;
 Double_t oldytrig=0;
 
 // Add Data Cells one by one to the sample to be partioned
 Int_t ncells=0;
 Int_t first=1;
 for (Int_t i=1; i<=n; i++)
 {
  ncells=i;
  prior=GetPrior(i,fpr);
  xup=hin->GetBinLowEdge(i)+hin->GetBinWidth(i);
  // Loop over all possible block partitions for this Data Cell sample
  first=1;
  for (Int_t j=1; j<=i; j++)
  {
   xlow=hin->GetBinLowEdge(j);
   blen=xup-xlow;
   bcount=0;
   bfit=0;
   a=0;
   b=0;
   if (fMode==1 || fMode==2)
   {
    bcount=hin->Integral(j,i);
    bfit=GetBlockFitness(fabs(bcount),blen);
   }
   if (fMode==3)
   {
    for (Int_t k=j; k<=i; k++)
    {
     yk=hin->GetBinContent(k);
     sigk=hin->GetBinError(k);
     sigk=sigk*sigk;
     if (sigk)
     {
      a+=1./(sigk);
      b+=yk/sigk;
     }
    }
    if (a)
    {
     bcount=b/a; // Weighted mean of the y-values in the block
     bfit=b*b/(2.*a); 
    }
   }
   pfit=prior+bfit;
   index=j-1; // Array index reduced by 1 because of C++ array index convention
   if (index>0) pfit+=best.At(index-1);
 
   // Record attributes for the optimal partition 
   if (first)
   {
    optfit=pfit;
    optj=j;
    optlen=blen;
    optcount=bcount;
    first=0;
   }
   else
   {
    if (pfit>optfit)
    {
     optfit=pfit;
     optj=j;
     optlen=blen;
     optcount=bcount;
    }
   }
  } // End of loop over possible block partitions
 
  // Store the attributes of the optimal partition
  best.SetAt(optfit,i-1);
  last.SetAt(optj,i-1);
  lengths.SetAt(optlen,i-1);
  counts.SetAt(optcount,i-1);
 
  if (!ntrig || optj==1) continue;
 
  // Check for triggering on a new change point
  oldoptj=last.At(i-2);
  oldoptlen=lengths.At(i-2);
  oldytrig=0;
  if ((fMode==1 || fMode==2) && oldoptlen) oldytrig=counts.At(i-2)/oldoptlen;
  if (fMode==3) oldytrig=counts.At(i-2);
  ytrig=0;
  if ((fMode==1 || fMode==2) && optlen) ytrig=counts.At(i-1)/optlen;
  if (fMode==3) ytrig=counts.At(i-1);
 
  if (optj>oldoptj && ((ntrig>0 && ytrig>oldytrig) || (ntrig<0 && ytrig<oldytrig)))
  {
   ncp++;
   xtrig=hin->GetBinLowEdge(optj);
  }
 
  // Stop when the requested number of triggers has been reached
  if (ncp>=abs(ntrig)) break;
 
 } // End of adding Data Cells one by one
 
 // Obtain the change points and corresponding partition information
 Int_t jcp=0;
 index=ncells;
 Double_t x=0;
 Double_t y=0;
 ncp=0;
 TArrayD xarr(ncells+1);
 TArrayD yarr(ncells+1);
 while (index>0)
 {
  index--; // Reduce array index by 1 because of C++ array convention
  jcp=last.At(index);
  ncp++;
  x=hin->GetBinLowEdge(jcp)+lengths.At(index);
  y=0;
  if (fMode==1 || fMode==2)
  {
   if (lengths.At(index)) y=counts.At(index)/lengths.At(index);
  }
  if (fMode==3) y=counts.At(index);
 
  xarr.SetAt(x,ncp-1);
  yarr.SetAt(y,ncp-1);
 
  // Also  mark the start of the first bin
  if (jcp==1) xarr.SetAt(hin->GetBinLowEdge(jcp),ncp);
 
  index=jcp-1;  
 }
 
 // Create the corresponding variable binned histogram
 Double_t* xbins=new Double_t[ncp+1];
 Double_t* yvals=new Double_t[ncp+1];
 
 for (Int_t i=0; i<=ncp; i++)
 {
  xbins[i]=xarr.At(ncp-i);
  yvals[i]=yarr.At(ncp-i);
 }
 
 hout->SetBins(ncp,xbins);
 for (Int_t i=1; i<=ncp; i++)
 {
  hout->SetBinContent(i,yvals[i]);
 }
 
 // Determine the first trigger point after the full sample analysis
 if (!ntrig) xtrig=xbins[1];
 
 hout->SetLineWidth(2);
 hout->SetLineColor(kBlue);
 hout->SetStats(kFALSE);
 
 // Set the output histogram and axes titles
 TString title,str;
 title="Bayesian Block representation for histogram ";
 title+=hin->GetName();
 title+=" with FPR= %-.3g";
 title+=";";
 str=hin->GetXaxis()->GetTitle();
 if (str=="") str="Recordings (e.g. time)";
 title+=str;
 title+=";";
 str=hin->GetYaxis()->GetTitle();
 if (str=="") str="Counts";
 title+=str;
 str=title.Format(title.Data(),fpr);
 hout->SetTitle(str.Data());
 
 // Indicate the requested trigger in a legend
 if (ntrig)
 {
  str.Form("Requested trigger at : %-.3g",xtrig);
 
  TLegend* leg=new TLegend(0.5,0.85,0.7,0.9,str);
  leg->SetFillColor(0);
  leg->SetTextColor(kBlue);
  leg->SetTextAlign(22);
 
  TList* hlist=hout->GetListOfFunctions();
  hlist->Add(leg);
 }
 
 // Reset the Data Mode for subsequent invokations
 fMode=0;
 
 delete [] xbins;
 delete [] yvals;
 
 return xtrig;
}

Double_t NcBlocks::GetBlocks(NcSample s,Int_t i,Double_t fpr,TH1* hout,Int_t ntrig)
{
 
 if (!hout)
 {
  cout << " *NcBlocks::GetBlocks* Error : Output histogram not specified." << endl;
  return 0;
 }
 
 Int_t n=s.GetN();
 Int_t store=s.GetStoreMode();
 Int_t dim=s.GetDimension();
 
 if (n<2 || !store || dim<1 || i<1 || i>dim || fpr<0 || fpr>1)
 {
  cout << " *NcBlocks::GetBlocks* Inconsistent input for NcSample treatment." << endl;
  cout << " Store Mode:" << store << " Entries:" << n << " Dimension:" << dim << " i:" << i << " fpr:" << fpr << endl; 
  return 0;
 }
 
 // Set Data Mode for unbinned event data
 fMode=1;
 
 // Represent each observation as 1 count in a variable binned histogram
 Double_t* xbins=new Double_t[n];
 Double_t val=0;
 for (Int_t idx=1; idx<=n; idx++)
 {
  val=s.GetEntry(idx,i,1,i);
  xbins[idx-1]=val;
 }
 
 TH1F hin("","",n-1,xbins);
 for (Int_t j=1; j<n; j++)
 {
  hin.SetBinContent(j,1);
 }
 
 Double_t xtrig=GetBlocks(&hin,fpr,hout,ntrig);
 
 // Set the output histogram and axes titles
 TString title,str;
 title="Bayesian Block representation for NcSample ";
 title+=s.GetName();
 title+=" with FPR=%-.3g";
 title+=";Recordings of variable ";
 title+=i;
 title+=" (";
 title+=s.GetVariableName(i);
 title+=")";
 title+=";Count rate";
 str=title.Format(title.Data(),fpr);
 hout->SetTitle(str.Data());
 
 delete [] xbins;
 
 return xtrig;
}

Double_t NcBlocks::GetBlocks(NcSample s,TString name,Double_t fpr,TH1* hout,Int_t ntrig)
{
 
 Int_t i=s.GetIndex(name);
 
 Double_t xtrig=GetBlocks(s,i,fpr,hout,ntrig);
 
 return xtrig;
}

Double_t NcBlocks::GetBlocks(Int_t n,Double_t* arr,Double_t fpr,TH1* hout,Int_t ntrig)
{
 
 if (!hout)
 {
  cout << " *NcBlocks::GetBlocks* Error : Output histogram not specified." << endl;
  return 0;
 }
 
 if (!arr || n<2 || fpr<0 || fpr>1)
 {
  cout << " *NcBlocks::GetBlocks* Inconsistent input for array treatment." << endl;
  if (!arr)
  {
   cout << " Array pointer is 0." << endl;
  }
  else
  {
   cout << " Entries:" << n <<  " fpr:" << fpr << endl;
  }
  return 0;
 }
 
 NcSample s;
 s.SetStoreMode();
 
 for (Int_t i=0; i<n; i++)
 {
  s.Enter(arr[i]);
 }
 
 Double_t xtrig=GetBlocks(s,1,fpr,hout,ntrig);
 
 // Set the output histogram and axes titles
 TString title;
 title.Form("Bayesian Block representation for unbinned array data with FPR=%-.3g",fpr);
 title+=";Recordings (e.g. time);Count rate";
 hout->SetTitle(title);
 
 return xtrig;
}

Double_t NcBlocks::GetBlocks(Int_t n,Float_t* arr,Double_t fpr,TH1* hout,Int_t ntrig)
{
 
 if (!hout)
 {
  cout << " *NcBlocks::GetBlocks* Error : Output histogram not specified." << endl;
  return 0;
 }
 
 if (!arr || n<2 || fpr<0 || fpr>1)
 {
  cout << " *NcBlocks::GetBlocks* Inconsistent input for array treatment." << endl;
  if (!arr)
  {
   cout << " Array pointer is 0." << endl;
  }
  else
  {
   cout << " Entries:" << n <<  " fpr:" << fpr << endl;
  }
  return 0;
 }
 
 NcSample s;
 s.SetStoreMode();
 
 for (Int_t i=0; i<n; i++)
 {
  s.Enter(arr[i]);
 }
 
 Double_t xtrig=GetBlocks(s,1,fpr,hout,ntrig);
 
 // Set the output histogram and axes titles
 TString title;
 title.Form("Bayesian Block representation for unbinned array data with FPR=%-.3g",fpr);
 title+=";Recordings (e.g. time);Count rate";
 hout->SetTitle(title);
 
 return xtrig;
}

Double_t NcBlocks::GetBlocks(TGraphErrors gr,Double_t fpr,TH1* hout,Int_t ntrig)
{
 
 if (!hout)
 {
  cout << " *NcBlocks::GetBlocks* Error : Output histogram not specified." << endl;
  return 0;
 }
 
 Int_t n=gr.GetN();
 
 if (n<2 || fpr<0 || fpr>1)
 {
  cout << " *NcBlocks::GetBlocks* Inconsistent input for TGraphErrors treatment." << endl;
  cout << " Entries:" << n << " fpr:" << fpr << endl; 
  return 0;
 }
 
 // Set Data Mode for measurements of a continuous observable
 fMode=3;
 
 // Sort the data points with increasing x-value
 gr.Sort();
 
 // Represent each observation as a value in a variable binned histogram
 Double_t* xbins=new Double_t[n+1];
 Double_t x=0;
 Double_t y=0;
 Double_t err=0;
 Double_t dist=0;
 Double_t dmin=-1;
 for (Int_t i=0; i<n; i++)
 {
  gr.GetPoint(i,x,y);
  err=fabs(gr.GetErrorX(i));
  xbins[i]=x-err;
  if (i>0)
  {
   dist=xbins[i]-xbins[i-1];
   if (dmin<0 || dist<dmin) dmin=dist;
  }
 }
 // Add an extra bin to contain the last measurement
 xbins[n]=xbins[n-1]+dmin;
 
 TH1F hin("","",n,xbins);
 for (Int_t j=1; j<=n; j++)
 {
  gr.GetPoint(j-1,x,y);
  err=fabs(gr.GetErrorY(j-1));
  hin.SetBinContent(j,y);
  hin.SetBinError(j,err);
 }
 
 Double_t xtrig=GetBlocks(&hin,fpr,hout,ntrig);
 
 // Set the output histogram and axes titles
 TString title,str;
 TString xtitle="Samplings (e.g. time)";
 TString ytitle="Measured value";
 TAxis* ax=gr.GetXaxis();
 if (ax)
 {
  str=ax->GetTitle();
  if (str != "") xtitle=str;
 }
 ax=gr.GetYaxis();
 if (ax)
 {
  str=ax->GetTitle();
  if (str != "") ytitle=str;
 }
 title="Bayesian Block representation for TGraphErrors ";
 title+=gr.GetName();
 title+=" with FPR=%-.3g;";
 title+=xtitle;
 title+=";";
 title+=ytitle;
 str=title.Format(title,fpr);
 hout->SetTitle(str);
 
 delete [] xbins;
 
 return xtrig;
}

Double_t NcBlocks::GetBlocks(TGraph gr,TF1 f,Double_t fpr,TH1* hout,Int_t ntrig)
{
 
 NcSample s;
 TGraphErrors gre=s.GetGraphErrors(&gr,0,0,0,&f);
 
 Double_t xtrig=GetBlocks(gre,fpr,hout,ntrig);
 
 TString str=" and input errors : ";
 str+=f.GetExpFormula("p");
 str.ReplaceAll("x","y");
 
 TString title=hout->GetTitle();
 title+=str;
 hout->SetTitle(title);
 
 return xtrig; 
}

Double_t NcBlocks::GetBlocks(TGraph gr,TString f,Double_t fpr,TH1* hout,Int_t ntrig)
{
 
 TF1 func("func",f);
 NcSample s;
 TGraphErrors gre=s.GetGraphErrors(&gr,0,0,0,&func);
 
 Double_t xtrig=GetBlocks(gre,fpr,hout,ntrig);
 
 TString str=" and input errors : ";
 str+=func.GetExpFormula("p");
 str.ReplaceAll("x","y");
 
 TString title=hout->GetTitle();
 title+=str;
 hout->SetTitle(title);
 
 return xtrig; 
}

Double_t NcBlocks::GetBlocks(TGraph gr,Double_t nrms,Double_t fpr,TH1* hout,Int_t ntrig)
{
 
 // Obtain the RMS deviation of all the y-values
 Double_t rms=gr.GetRMS(2);
 
 // Determine the error for each y-value and convert into a function format
 Double_t err=fabs(nrms*rms);
 
 TString f;
 f.Form("%-.5g",err);
 
 Double_t xtrig=GetBlocks(gr,f,fpr,hout,ntrig);
 
 TString str;
 str.Form(" from nrms=%-.3g",fabs(nrms));
 
 TString title=hout->GetTitle();
 title+=str;
 hout->SetTitle(title);
 
 return xtrig; 
}

Int_t NcBlocks::GetBlocks(TH1* hin,TH1* hout,Int_t n,Int_t mode)
{
 
 if (!hin)
 {
  cout << " *NcBlocks::GetBlocks* Error : Input histogram not specified." << endl;
  return 0;
 }
 
 if (!hout)
 {
  cout << " *NcBlocks::GetBlocks* Error : Output histogram not specified." << endl;
  return 0;
 }
 
 Int_t nbins=hin->GetNbinsX();
 
 if (!nbins || n<1 || n>nbins || mode<0 || mode>2)
 {
  cout << " *NcBlocks::GetBlocks* Inconsistent input nbins=" << nbins << " n=" << n << " mode=" << mode << endl;
  return 0;
 }
 
 // Retrieve the various sets of "n" bins from the input histogram
 Int_t jbin=0;
 Double_t x=0;
 Double_t y=0;
 Double_t xlow=0;
 Double_t xup=0;
 Double_t binwidth=0;
 NcSample s;
 if (mode==1) s.SetStoreMode(1);
 Int_t nblocks=0;
 Double_t average=0;
 TArrayD xarr(nbins);
 TArrayD yarr(nbins);
 while (jbin<nbins)
 {
  for (Int_t i=0; i<n; i++)
  {
   jbin++;
 
   if (jbin>nbins) break;
 
   x=hin->GetBinCenter(jbin);
   y=hin->GetBinContent(jbin);
   binwidth=hin->GetBinWidth(jbin);
   if (i==0) xlow=hin->GetBinLowEdge(jbin);
   xup=x+0.5*binwidth;
   s.Enter(x,y);
  }
  if (mode==0) average=s.GetMean(2);
  if (mode==1) average=s.GetMedian(2);
  if (mode==2) average=s.GetRMS(2);
  nblocks++;
  xarr.SetAt(xlow,nblocks-1);
  yarr.SetAt(average,nblocks-1);
  s.Reset();
 } // End of while loop
 
 // Create the corresponding variable binned output histogram
 Double_t* xbins=new Double_t[nblocks+1];
 Double_t* yvals=new Double_t[nblocks+1];
 
 for (Int_t i=0; i<nblocks; i++)
 {
  xbins[i]=xarr.At(i);
  yvals[i]=yarr.At(i);
 }
 // Add an extra bin to contain the last data
 xbins[nblocks]=(1.+1e-6)*xup;
 yvals[nblocks]=yvals[nblocks-1];
 
 hout->SetBins(nblocks,xbins);
 for (Int_t i=1; i<=nblocks; i++)
 {
  hout->SetBinContent(i,yvals[i-1]);
 }
 
 hout->SetLineWidth(2);
 hout->SetLineColor(kBlue);
 hout->SetStats(kFALSE);
 
 // Set the output histogram and axes titles
 TString title,str;
 title="Block representation for histogram ";
 title+=hin->GetName();
 title+=" grouped in %-d consecutive bins";
 title+=";";
 str=hin->GetXaxis()->GetTitle();
 if (str=="") str="Recordings (e.g. time)";
 title+=str;
 if (mode==0) title+=";Mean ";
 if (mode==1) title+=";Median ";
 if (mode==2) title+=";RMS ";
 str=hin->GetYaxis()->GetTitle();
 if (str=="") str="Counts";
 title+=str;
 str=title.Format(title.Data(),n);
 hout->SetTitle(str.Data());
 
 delete [] xbins;
 delete [] yvals;
 
 return nblocks;
}

Int_t NcBlocks::GetBlocks(NcSample s,Int_t i,TH1* hout,Int_t n,Int_t mode)
{
 
 if (!hout)
 {
  cout << " *NcBlocks::GetBlocks* Error : Output histogram not specified." << endl;
  return 0;
 }
 
 Int_t nen=s.GetN();
 Int_t store=s.GetStoreMode();
 Int_t dim=s.GetDimension();
 
 if (!store || dim<1 || i<1 || i>dim || n<1 || n>nen || mode<0 || mode>2)
 {
  cout << " *NcBlocks::GetBlocks* Inconsistent input for NcSample treatment." << endl;
  cout << " Store Mode:" << store << " Entries:" << nen << " Dimension:" << dim << " i:" << i << " n:" << n << " mode:" << mode << endl; 
  return 0;
 }
 
 TGraph gr=s.GetGraph(i);
 
 Int_t nblocks=GetBlocks(&gr,hout,n,mode);
 
 // Set the output histogram and axes titles
 TString title,str;
 title="Block representation for NcSample ";
 title+=s.GetName();
 title+=" grouped in %-d consecutive samples";
 title+=";Sampling number";
 if (mode==0) title+=";Mean ";
 if (mode==1) title+=";Median ";
 if (mode==2) title+=";RMS ";
 title+="of variable ";
 title+=i;
 title+=" (";
 title+=s.GetVariableName(i);
 title+=")";
 str=title.Format(title.Data(),n);
 hout->SetTitle(str.Data());
 
 return nblocks;
}

Int_t NcBlocks::GetBlocks(NcSample s,TString name,TH1* hout,Int_t n,Int_t mode)
{
 
 Int_t i=s.GetIndex(name);
 
 Int_t nblocks=GetBlocks(s,i,hout,n,mode);
 
 return nblocks;
}

Int_t NcBlocks::GetBlocks(Int_t nr,Double_t* arr,TH1* hout,Int_t n,Int_t mode)
{
 
 if (!hout)
 {
  cout << " *NcBlocks::GetBlocks* Error : Output histogram not specified." << endl;
  return 0;
 }
 
 if (!arr || n<1 || n>nr || mode<0 || mode>2)
 {
  cout << " *NcBlocks::GetBlocks* Inconsistent input for array treatment." << endl;
  if (!arr)
  {
   cout << " Array pointer is 0." << endl;
  }
  else
  {
   cout << " Entries:" << nr << " n:" << n << " mode:" << mode << endl;
  }
  return 0;
 }
 
 NcSample s;
 s.SetStoreMode();
 
 for (Int_t i=0; i<nr; i++)
 {
  s.Enter(arr[i]);
 }
 
 Int_t nblocks=GetBlocks(s,1,hout,n,mode);
 
 // Set the output histogram and axes titles
 TString title;
 title.Form("Block representation for array data grouped in %-d consecutive recordings",n);
 title+=";Sampling number";
 if (mode==0) title+=";Mean ";
 if (mode==1) title+=";Median ";
 if (mode==2) title+=";RMS ";
 title+="value";
 hout->SetTitle(title);
 
 return nblocks;
}

Int_t NcBlocks::GetBlocks(Int_t nr,Float_t* arr,TH1* hout,Int_t n,Int_t mode)
{
 
 if (!hout)
 {
  cout << " *NcBlocks::GetBlocks* Error : Output histogram not specified." << endl;
  return 0;
 }
 
 if (!arr || n<1 || n>nr || mode<0 || mode>2)
 {
  cout << " *NcBlocks::GetBlocks* Inconsistent input for array treatment." << endl;
  if (!arr)
  {
   cout << " Array pointer is 0." << endl;
  }
  else
  {
   cout << " Entries:" << nr << " n:" << n << " mode:" << mode << endl;
  }
  return 0;
 }
 
 NcSample s;
 s.SetStoreMode();
 
 for (Int_t i=0; i<nr; i++)
 {
  s.Enter(arr[i]);
 }
 
 Int_t nblocks=GetBlocks(s,1,hout,n,mode);
 
 // Set the output histogram and axes titles
 TString title;
 title.Form("Block representation for array data grouped in %-d consecutive recordings",n);
 title+=";Sampling number";
 if (mode==0) title+=";Mean ";
 if (mode==1) title+=";Median ";
 if (mode==2) title+=";RMS ";
 title+="value";
 hout->SetTitle(title);
 
 return nblocks;
}

Int_t NcBlocks::GetBlocks(TGraph* gr,TH1* hout,Int_t n,Int_t mode)
{
 
 if (!gr)
 {
  cout << " *NcBlocks::GetBlocks* Error : Input TGraph not specified." << endl;
  return 0;
 }
 
 if (!hout)
 {
  cout << " *NcBlocks::GetBlocks* Error : Output histogram not specified." << endl;
  return 0;
 }
 
 if (n<1 || mode<0 || mode>2)
 {
  cout << " *NcBlocks::GetBlocks* Inconsistent input for TGraph treatment : n=" << n << " mode=" << mode << endl;
  return 0;
 }
 
 Int_t npoints=gr->GetN();
 
 if (!npoints) return 0;
 
 // Sort the data points with increasing x-value
 gr->Sort();
 
 // Represent each observation as a value in a variable binned histogram
 Double_t* xbins=new Double_t[npoints+1];
 Double_t x=0;
 Double_t y=0;
 for (Int_t i=0; i<npoints; i++)
 {
  gr->GetPoint(i,x,y);
  xbins[i]=x;
 }
 // Add an extra bin to contain the last measurement
 xbins[npoints]=(1.+1e-6)*xbins[npoints-1];
 
 TH1F hin("","",npoints,xbins);
 for (Int_t j=1; j<=npoints; j++)
 {
  gr->GetPoint(j-1,x,y);
  hin.SetBinContent(j,y);
 }
 
 Int_t nblocks=GetBlocks(&hin,hout,n,mode);
 
 // Set the output histogram and axes titles
 TString title,str;
 TString xtitle="Samplings (e.g. time)";
 TString ytitle="Measured value";
 TAxis* ax=gr->GetXaxis();
 if (ax)
 {
  str=ax->GetTitle();
  if (str != "") xtitle=str;
 }
 ax=gr->GetYaxis();
 if (ax)
 {
  str=ax->GetTitle();
  if (str != "") ytitle=str;
 }
 title="Block representation for TGraph ";
 title+=gr->GetName();
 title+=" grouped in %-d consecutive samples;";
 title+=xtitle;
 if (mode==0) title+=";Mean ";
 if (mode==1) title+=";Median ";
 if (mode==2) title+=";RMS ";
 title+=ytitle;
 str=title.Format(title,n);
 hout->SetTitle(str);
 
 delete [] xbins;
 
 return nblocks;
}

Int_t NcBlocks::Add(TH1* h1,TH1* h2,TH1* hout,Bool_t scale,Double_t c,Double_t d)
{
 
 if (hout) hout->Reset();
 
 if (!h1 || !h2 || !hout) return 1;
 
 Int_t ndim1=h1->GetDimension();
 Int_t ndim2=h2->GetDimension();
 Int_t ndimo=hout->GetDimension();
 
 if (ndim1!=1 || ndim2!=1 || ndimo!=1)
 {
  cout << " *NcBlocks::Add* Error : Histograms should all be 1-dimensional." << endl;
  return 1;
 }
 
 // Make the X-axis of "hout" identical to the X-axis of "h1"
 TString name=hout->GetName();
 h1->Copy(*hout);
 hout->Reset();
 
 TString name1=h1->GetName();
 if (name1=="") name1="h1";
 
 TString name2=h2->GetName();
 if (name2=="") name2="h2";
 
 TString sc="+";
 if (c<0) sc="-";
 if (fabs(c)!=1)
 {
  sc.Form("%-+.3g",c);
  sc+="*";
 }
 
 TString sd="";
 sd.Form("%-+.3g",d);
 
 TString title="Resulting histogram of: ";
 title+=name1;
 if (c)
 {
  title+=sc;
  title+=name2;
 }
 if (d) title+=sd;
 if (scale)
 {
  title+=" (scaled w.r.t. bin size)"; 
 }
 else
 {
  title+=" (not scaled w.r.t. bin size)"; 
 }
 title+=";";
 TAxis* axis1=h1->GetXaxis();
 title+=axis1->GetTitle();
 title+=";";
 axis1=h1->GetYaxis();
 title+=axis1->GetTitle();
 
 hout->SetName(name);
 hout->SetTitle(title);
 
 Int_t nb1=h1->GetNbinsX();
 Int_t nb2=h2->GetNbinsX();
 
 if (!nb1 || !nb2) return 1;
 
 // Get the largest bin size of "h1"
 Double_t bwidth1=0;
 Double_t bwmax1=0;
 Int_t imax1=0;
 for (Int_t i=1; i<=nb1; i++)
 {
  bwidth1=h1->GetBinWidth(i);
  if (i==1 || bwidth1>bwmax1)
  {
   bwmax1=bwidth1;
   imax1=i;
  }
 }
 
 // Get the smallest bin size of "h2" 
 Double_t bwidth2=0;
 Double_t bwmin2=0;
 Int_t imin2=0;
 for (Int_t i=1; i<=nb2; i++)
 {
  bwidth2=h2->GetBinWidth(i);
  if (i==1 || bwidth2<bwmin2)
  {
   bwmin2=bwidth2;
   imin2=i;
  }
 }
 
 Double_t ratio=bwmax1/bwmin2;
 if (ratio>1.001)
 {
  printf(" *NcBlocks::Add* Error : Larger bin size encountered in histogram %-s than in %-s \n",name1.Data(),name2.Data());
  printf(" %-s: binsize=%-g for bin=%-i  %-s: binsize=%-g for bin=%-i \n",name1.Data(),bwmax1,imax1,name2.Data(),bwmin2,imin2);
  return 1;
 }
 
 // Loop over all the bins of the input histogram h1
 Double_t x1=0;
 Double_t y1=0;
 Int_t i2=0;
 Double_t y2=0;
 Double_t ynew=0;
 TAxis* axis2=h2->GetXaxis();
 for (Int_t i1=1; i1<=nb1; i1++)
 {
  x1=h1->GetBinCenter(i1);
  y1=h1->GetBinContent(i1);
  bwidth1=h1->GetBinWidth(i1);
  i2=axis2->FindFixBin(x1);
  y2=h2->GetBinContent(i2);
  bwidth2=h2->GetBinWidth(i2);
 
  // Do not consider underflow or overflow bins
  if (i2<1 || i2>nb2) continue;
 
  if (scale) y2*=bwidth1/bwidth2;
  ynew=y1+c*y2+d;
  hout->SetBinContent(i1,ynew);
 }
 
 return 0;
}

Int_t NcBlocks::Add(TGraph* gr,TH1* h,TGraph* gout,Double_t c,Double_t d)
{
 
 if (gout) gout->Set(0);
 
 if (!gr || !h || !gout) return 1;
 
 TString nameg=gr->GetName();
 if (nameg=="") nameg="gr";
 
 TString nameh=h->GetName();
 if (nameh=="") nameh="h";
 
 TString sc="+";
 if (c<0) sc="-";
 if (fabs(c)!=1)
 {
  sc.Form("%-+.3g",c);
  sc+="*";
 }
 
 TString sd="";
 sd.Form("%-+.3g",d);
 
 TString title="Resulting graph of: ";
 title+=nameg;
 if (c)
 {
  title+=sc;
  title+=nameh;
 }
 if (d) title+=sd;
 title+=";";
 TAxis* axis=0;
 axis=gr->GetXaxis();
 title+=axis->GetTitle();
 title+=";";
 axis=gr->GetYaxis();
 title+=axis->GetTitle();
 
 gout->SetTitle(title);
 
 Int_t ndim=h->GetDimension();
 
 if (ndim!=1)
 {
  cout << " *NcBlocks::Add* Error : Histogram " << nameh << " should be 1-dimensional." << endl;
  return 1;
 }
 
 Int_t np=gr->GetN();
 Int_t nb=h->GetNbinsX();
 
 if (!np || !nb) return 1;
 
 // Loop over all the points in the graph
 Double_t x=0;
 Double_t y=0;
 Double_t ex=0;
 Double_t ey=0;
 Int_t hbin=0;
 Double_t hval=0;
 Double_t ynew=0;
 axis=h->GetXaxis();
 for (Int_t i=0; i<np; i++)
 {
  gr->GetPoint(i,x,y);
  hbin=axis->FindFixBin(x);
  hval=h->GetBinContent(hbin);
 
  // Do not consider underflow or overflow bins
  if (hbin<1 || hbin>nb) continue;
 
  ynew=y+c*hval+d;
  gout->SetPoint(i,x,ynew);
 
  if (gr->InheritsFrom("TGraphErrors") && gout->InheritsFrom("TGraphErrors"))
  {
   TGraphErrors* gre=(TGraphErrors*)gr;
   TGraphErrors* goute=(TGraphErrors*)gout;
   ex=gre->GetErrorX(i);
   ey=gre->GetErrorY(i);
   goute->SetPointError(i,ex,ey);
  }
 }
 
 return 0;
}

Int_t NcBlocks::Divide(TH1* h1,TH1* h2,TH1* hout,Bool_t scale,Double_t c,Double_t d)
{
 
 if (hout) hout->Reset();
 
 if (!c)
 {
  printf(" *NcBlocks::Divide* Error : Invalid value c=0. \n");
  return 1;
 }
 
 if (!h1 || !h2 || !hout) return 1;
 
 Int_t ndim1=h1->GetDimension();
 Int_t ndim2=h2->GetDimension();
 Int_t ndimo=hout->GetDimension();
 
 if (ndim1!=1 || ndim2!=1 || ndimo!=1)
 {
  cout << " *NcBlocks::Divide* Error : Histograms should all be 1-dimensional." << endl;
  return 1;
 }
 
 // Make the X-axis of "hout" identical to the X-axis of "h1"
 TString name=hout->GetName();
 h1->Copy(*hout);
 hout->Reset();
 
 TString name1=h1->GetName();
 if (name1=="") name1="h1";
 
 TString name2=h2->GetName();
 if (name2=="") name2="h2";
 
 TString sc="/";
 if (fabs(c)!=1) sc.Form("%-.3g*",fabs(c));
 
 TString sd;
 if (c>0)
 {
  sd="";
  if (d) sd.Form("%-.3g+",d);
 }
 else
 {
  sd="-";
  if (d) sd.Form("%-.3g-",d);
 }
 
 TString title="Resulting histogram of: ";
 if (sd!="") title+=sd;
 if (c)
 {
  title+=name1;
  if (sc.Contains("*")) title+="/(";
  title+=sc;
  title+=name2;
  if (sc.Contains("*")) title+=")";
 }
 if (scale)
 {
  title+=" (scaled w.r.t. bin size)"; 
 }
 else
 {
  title+=" (not scaled w.r.t. bin size)"; 
 }
 title+=";";
 TAxis* axis1=h1->GetXaxis();
 title+=axis1->GetTitle();
 title+=";";
 axis1=h1->GetYaxis();
 title+=axis1->GetTitle();
 
 hout->SetName(name);
 hout->SetTitle(title);
 
 Int_t nb1=h1->GetNbinsX();
 Int_t nb2=h2->GetNbinsX();
 
 if (!nb1 || !nb2) return 1;
 
 // Get the largest bin size of "h1"
 Double_t bwidth1=0;
 Double_t bwmax1=0;
 Int_t imax1=0;
 for (Int_t i=1; i<=nb1; i++)
 {
  bwidth1=h1->GetBinWidth(i);
  if (i==1 || bwidth1>bwmax1)
  {
   bwmax1=bwidth1;
   imax1=i;
  }
 }
 
 // Get the smallest bin size of "h2" 
 Double_t bwidth2=0;
 Double_t bwmin2=0;
 Int_t imin2=0;
 for (Int_t i=1; i<=nb2; i++)
 {
  bwidth2=h2->GetBinWidth(i);
  if (i==1 || bwidth2<bwmin2)
  {
   bwmin2=bwidth2;
   imin2=i;
  }
 }
 
 Double_t ratio=bwmax1/bwmin2;
 if (ratio>1.001)
 {
  printf(" *NcBlocks::Divide* Error : Larger bin size encountered in histogram %-s than in %-s \n",name1.Data(),name2.Data());
  printf(" %-s: binsize=%-g for bin=%-i  %-s: binsize=%-g for bin=%-i \n",name1.Data(),bwmax1,imax1,name2.Data(),bwmin2,imin2);
  return 1;
 }
 
 // Loop over all the bins of the input histogram h1
 Double_t x1=0;
 Double_t y1=0;
 Int_t i2=0;
 Double_t y2=0;
 Double_t val=0;
 Double_t ynew=0;
 TAxis* axis2=h2->GetXaxis();
 for (Int_t i1=1; i1<=nb1; i1++)
 {
  x1=h1->GetBinCenter(i1);
  y1=h1->GetBinContent(i1);
  bwidth1=h1->GetBinWidth(i1);
  i2=axis2->FindFixBin(x1);
  y2=h2->GetBinContent(i2);
  bwidth2=h2->GetBinWidth(i2);
 
  // Do not consider underflow or overflow bins
  if (i2<1 || i2>nb2) continue;
 
  val=c*y2;
  if (!val) continue;
 
  if (scale) y2*=bwidth1/bwidth2;
  ynew=d+y1/val;
  hout->SetBinContent(i1,ynew);
 }
 
 return 0;
}

Int_t NcBlocks::Divide(TGraph* gr,TH1* h,TGraph* gout,Double_t c,Double_t d)
{
 
 if (gout) gout->Set(0);
 
 if (!c)
 {
  printf(" *NcBlocks::Divide* Error : Invalid value c=0. \n");
  return 1;
 }
 
 if (!gr || !h || !gout) return 1;
 
 TString nameg=gr->GetName();
 if (nameg=="") nameg="gr";
 
 TString nameh=h->GetName();
 if (nameh=="") nameh="h";
 
 TString sc="/";
 if (fabs(c)!=1) sc.Form("%-.3g*",fabs(c));
 
 TString sd;
 if (c>0)
 {
  sd="";
  if (d) sd.Form("%-.3g+",d);
 }
 else
 {
  sd="-";
  if (d) sd.Form("%-.3g-",d);
 }
 
 TString title="Resulting graph of: ";
 if (sd!="") title+=sd;
 if (c)
 {
  title+=nameg;
  if (sc.Contains("*")) title+="/(";
  title+=sc;
  title+=nameh;
  if (sc.Contains("*")) title+=")";
 }
 title+=";";
 TAxis* axis=0;
 axis=gr->GetXaxis();
 title+=axis->GetTitle();
 title+=";";
 axis=gr->GetYaxis();
 title+=axis->GetTitle();
 
 gout->SetTitle(title);
 
 Int_t ndim=h->GetDimension();
 
 if (ndim!=1)
 {
  cout << " *NcBlocks::Divide* Error : Histogram " << nameh << " should be 1-dimensional." << endl;
  return 1;
 }
 
 Int_t np=gr->GetN();
 Int_t nb=h->GetNbinsX();
 
 if (!np || !nb) return 1;
 
 // Loop over all the points in the graph
 Double_t x=0;
 Double_t y=0;
 Double_t ex=0;
 Double_t ey=0;
 Int_t hbin=0;
 Double_t hval=0;
 Double_t val=0;
 Double_t ynew=0;
 axis=h->GetXaxis();
 Int_t j=0;
 for (Int_t i=0; i<np; i++)
 {
  gr->GetPoint(i,x,y);
  hbin=axis->FindFixBin(x);
  hval=h->GetBinContent(hbin);
 
  // Do not consider underflow or overflow bins
  if (hbin<1 || hbin>nb) continue;
 
  val=c*hval;
  if (!val) continue;
 
  ynew=d+y/val;
  gout->SetPoint(j,x,ynew);
 
  if (gr->InheritsFrom("TGraphErrors") && gout->InheritsFrom("TGraphErrors"))
  {
   TGraphErrors* gre=(TGraphErrors*)gr;
   TGraphErrors* goute=(TGraphErrors*)gout;
   ex=gre->GetErrorX(i);
   ey=gre->GetErrorY(i);
   goute->SetPointError(j,ex,ey);
  }
 
  j++;
 }
 
 return 0;
}

Int_t NcBlocks::Rebin(TH1* hin,TH1* hout,Bool_t scale,Int_t nbins,Double_t xmin,Double_t xmax)
{
 
 if (hout) hout->Reset();
 
 if (!hin || !hout) return 0;
 
 Int_t ndimi=hin->GetDimension();
 Int_t ndimo=hout->GetDimension();
 
 if (ndimi!=1 || ndimo!=1)
 {
  cout << " *NcBlocks::Rebin* Error : Histograms should both be 1-dimensional." << endl;
  return 0;
 }
 
 Int_t nb1=hin->GetNbinsX();
 
 if (!nb1) return 0;
 
 TAxis* xaxis=hin->GetXaxis();
 TAxis* yaxis=hin->GetYaxis();
 
 if (xmax<xmin)
 {
  xmin=xaxis->GetXmin();
  xmax=xaxis->GetXmax();
 }
 
 Double_t bwidth1=0;
 Double_t xlow1=0;
 Double_t xup1=0;
 
 // Automatic binwidth setting
 if (nbins<=0)
 {
  Double_t bwmin=-1;
  for (Int_t i=1; i<=nb1; i++)
  {
   bwidth1=hin->GetBinWidth(i);
   xlow1=hin->GetBinLowEdge(i);
   xup1=xlow1+bwidth1;
 
   if (xlow1>=xmax || xup1<=xmin) continue;
 
   if (bwmin<0 || bwidth1<bwmin) bwmin=bwidth1;
  }
 
  if (bwmin<0)
  {
   cout << " *NcBlocks::Rebin* Error : Input histogram had no data in requested interval [xmin,xmax]." << endl;
   return 0;
  }
 
  Double_t rnbins=(xmax-xmin)/bwmin;
  nbins=int(rnbins);
  Double_t diff=rnbins-double(nbins);
  if (diff) nbins=nbins+1;
 }
 
 hout->SetBins(nbins,xmin,xmax);
 Double_t bwidth=hout->GetBinWidth(1);
 
 TString name=hin->GetName();
 if (name=="") name="hin";
 
 TString snb="";
 snb.Form("   nbins=%-i",nbins);
 
 TString smin="";
 smin.Form(" xmin=%-.3g",xmin);
 
 TString smax="";
 smax.Form(" xmax=%-.3g",xmax);
 
 TString title="Uniformly binned version of histogram: ";
 title+=name;
 title+=snb;
 title+=smin;
 title+=smax;
 if (scale)
 {
  title+=" (scaled w.r.t. bin size)"; 
 }
 else
 {
  title+=" (not scaled w.r.t. bin size)"; 
 }
 title+=";";
 title+=xaxis->GetTitle();
 title+=";";
 title+=yaxis->GetTitle();
 
 hout->SetTitle(title);
 
 // Check if the binning of "hout" is not coarser than that of "hin" within [xmin,xmax]
 for (Int_t i=1; i<=nb1; i++)
 {
  bwidth1=hin->GetBinWidth(i);
  xlow1=hin->GetBinLowEdge(i);
  xup1=xlow1+bwidth1;
 
  if (xlow1>=xmax || xup1<=xmin) continue;
 
  if (bwidth1<bwidth)
  {
  printf(" *NcBlocks::Rebin* Error : Input histogram had finer binning than output histogram. \n");
  printf(" Input: binwidth=%-g for bin=%-i  Output: uniform binwidth=%-g \n",bwidth1,i,bwidth);
  return 0;
  }
 }
 
 // Loop over all the bins of the output histogram hout
 Double_t x=0;
 Int_t i1=0;
 Double_t y1=0;
 for (Int_t i=1; i<=nbins; i++)
 {
  x=hout->GetBinCenter(i);
  i1=xaxis->FindFixBin(x);
 
  // Do not consider underflow or overflow bins
  if (i1<1 || i1>nb1) continue;
 
  y1=hin->GetBinContent(i1);
  bwidth1=hin->GetBinWidth(i1);
  if (scale) y1=y1*bwidth/bwidth1;
  hout->SetBinContent(i,y1);
 }
 
 return nbins;
}