// root.cern.ch/root/htmldoc/guides/spectrum/Spectrum.html#smoothing

// Attention: ROOT implements the Markov smoothing, not the simple convolution.
//
// See: Z.K. Silagadze, Nucl. Instr. Meth. A 376, 451 (1996).

// Example to illustrate smoothing using Markov algorithm (class TSpectrum).
//
// Author: Miroslav Morhac


TH1F* generate_spectrum (int NBins)
{
  TH1F *hspect = new TH1F ("hspect", "Input data", NBins, 0., NBins);

  for (int i = 0; i < NBins*10; i++)
    hspect->Fill ( gRandom->Uniform (NBins) );

  TF1 fgauss ("fgauss", "TMath::Gaus (x,  [0] ,  [1]  , 1 )", 0., NBins);
  fgauss.SetParameter (1, NBins/100);

  for (int PeakNum = 1; PeakNum<= 5; PeakNum++) 
  {
    int NeventsPeak = gRandom->Uniform (0.3*NBins, 1.0*NBins);

    fgauss.SetParameter (0, (PeakNum-0.5)/5. * NBins);

    for (int EventNum = 0; EventNum < NeventsPeak; EventNum++)
      hspect->Fill ( fgauss.GetRandom () );
  }
  return hspect;
}



void Smoothing (Int_t averWindow = 5)
{
  Int_t i;
  Int_t nbins = 1024;
  Double_t xmin  = 0;
  Double_t xmax  = nbins;
  Double_t* source = new Double_t [nbins];
  gROOT->ForceStyle();

  TH1F* hInput  = generate_spectrum (nbins);

  for (i = 0; i < nbins; i++) source[i] = hInput->GetBinContent (i+1);
   
  hInput->SetAxisRange (1, 1024);

  TSpectrum* Analyser1dim = new TSpectrum ();

  Analyser1dim->SmoothMarkov (source, 1024, averWindow);  //3, 7, 10, ...

/*
   const char* TSpectrum::SmoothMarkov (
        Double_t* source, : pointer to the array of source spectrum
	Int_t ssize,      : length of source array
	Int_t averWindow  : width of averaging smoothing window
   )	
  
  root.cern.ch/doc/v608/classTSpectrum.html#a2767b767d39ee7826e9a50b9c0359beb

  "This function calculates smoothed spectrum from source spectrum based on Markov chain
   method. The result is placed in the array pointed by source pointer. On successful
   completion it returns 0. On error it returns pointer to the string describing error."
*/

  TH1F* hSmoMark = new TH1F ("hsmomarkov", Form ("Smoothed via Markov for averWindow = %d", averWindow) , nbins, 0., nbins);
  hSmoMark->SetLineColor (kRed);

  for (i = 0; i < nbins; i++) hSmoMark->SetBinContent (i+1, source[i]);

  TCanvas* c1 = new TCanvas ("c1", "TSpectrum -> Markov smoothing", 15, 10, 600, 480);

  hSmoMark->SetMaximum ( TMath::Max (hInput->GetMaximum() , hSmoMark->GetMaximum()) );
  hSmoMark->Draw ("L");
  hInput  ->Draw ("L SAME");


  /*****              Another smoothing procedure:                 *****/
  /*                                                                   */
  /* TH1::Smooth implementing the 353QH algorithm. From HBOOK's hsmoof */
  /*                                                                   */
  /*********************************************************************/
  
  TCanvas* c2 = new TCanvas ("c2", "TH1 -> 353QH smoothing", 600, 10, 600, 480);

  TH1F* hSmo353QH = (TH1F*) hInput->Clone ();
  hSmo353QH->SetNameTitle ("hsmo353qh", "Smoothed via 353QH algorithm");

  hSmo353QH->Smooth ( 1 );
  
/*
  void TH1::Smooth (
    Int_t      ntimes = 1,  : how many times to repeat
    Option_t*  option = ""  : for "R", smoothing is applied only to the bins defined in the X axis range
  )	

  See: root.cern.ch/doc/master/classTH1.html#a0d08651c37b622f4bcc0e1a0affefb33
*/
  
  hInput->SetMaximum ( TMath::Max (hInput->GetMaximum() , hSmo353QH->GetMaximum()) );
  hSmo353QH->SetLineColor (2);
  hSmo353QH->SetLineWidth (2);

  hSmo353QH->Draw ("L");
  hInput   ->Draw ("L SAME");
  hSmo353QH->Draw ("L SAME");
}