// root.cern.ch/doc/v608/classTSpectrum.html
// root.cern.ch/root/htmldoc/guides/spectrum/Spectrum.html

// Illustrates how to find peaks in histograms.
//
// This script generates a random number of gaussian peaks
// on top of a linear background.
// The positions of the peaks are found via TSpectrum
// and injected as initial values of parameters to make a global fit.
// The background is computed and drawn on top of the original histogram.
//
// To execute this example, do (in ROOT 5 or ROOT 6):
//
// ~~~{.cpp}
//  root > .x peaks.C  (generate 10 peaks by default)
//  root > .x peaks.C++ (use the compiler)
//  root > .x peaks.C++(30) (generates 30 peaks)
// ~~~
//
// Author: Rene Brun


#include "TCanvas.h"
#include "TMath.h"
#include "TH1.h"
#include "TF1.h"
#include "TRandom.h"
#include "TSpectrum.h"
#include "TVirtualFitter.h"


Int_t NPeaks = 10;


Double_t funPeaks (Double_t *x, Double_t *par)
{
   Double_t Result = par[0] + par[1]*x[0];

   for (Int_t p = 0; p < NPeaks; p++) {
      Double_t Norm  = par [3*p + 2];
      Double_t Mean  = par [3*p + 3];
      Double_t Sigma = par [3*p + 4];
      Result += Norm * TMath::Gaus (x[0], Mean, Sigma);
   }
   return Result;
}


void peaks (Int_t NPeaks = 10)
{
  TH1F* hspect = new TH1F ("hspect", "test", 500, 0, 1000);

  // Generate [NPeaks] peaks at random

  TF1* funModel = new TF1 ("funModel", funPeaks, 0, 1000, 2+3*NPeaks);
  funModel->SetNpx (10000);

  Double_t par [3000];
  par[0] =  0.8;
  par[1] = -0.6/1000;

  for (Int_t p = 0; p < NPeaks; p++) {
    par [3*p + 2] =  1;                         // "norm"
    par [3*p + 3] = 40 +   gRandom->Rndm()*920; // "mean"
    par [3*p + 4] =  3 + 2*gRandom->Rndm();     // "sigma"
  }
  funModel->SetParameters (par);

  TCanvas* c1 = new TCanvas ("c1", "c1", 10, 10, 1000, 900);
  c1->Divide (1,2);
  c1->cd (1);
  hspect->FillRandom ("funModel", 200000);
  hspect->Draw ();
  TH1F* h2 = (TH1F*) hspect->Clone ("h2");

  /***                                                       ***/
  /*** From here we will treat hspect as unknown spectrum.   ***/
  /*** We intend to find Bkgnd, Peaks, and later - Fit peaks ***/
  /***                                                       ***/


  // Use TSpectrum to find the peak candidates

  TSpectrum* Analyser1dim = new TSpectrum (2*NPeaks);

  Int_t NFound = Analyser1dim->Search (hspect, 2, "", 0.10);

/*
  Int_t TSpectrum::Search (
   const TH1* hin,             : pointer to the histogram of source spectrum
   Double_t   sigma = 2,       : sigma of searched peaks, see manual for details
   Option_t*  option = "",     :
   Double_t   threshold = 0.05 : [0,1]. Peaks with Amp <= threshold*HighestPeak are discarded. 
  )

  The number of found peaks and their positions are written into members:
      fNpeaks, fPositionX and fPositionY.
  Positions are retrievable via TSpectrum::GetPositionX() and GetPositionY() methods,
  which return the arrays of resp. X and Y positions of found peaks.

  The search is performed in the current histogram range.

  Options:

  - Background is removed by default.          Otherwise: "nobackground"
  - Smoothing is done by default.              Otherwise: "nomarkov"
  - Polymarker with peaks is added by default. Otherwise: "goff"
  - The resulting TH1 is drawn by default.     Otherwise: "nodraw"

  See: root.cern.ch/doc/master/classTSpectrum.html#a5f8b7b208f813670259b51a8bcdcbfc5
*/

  printf ("Found %d candidate peaks to fit\n", NFound);

  // Estimate background using TSpectrum::Background

  TH1* hbkgnd = Analyser1dim->Background (hspect, 20, "same");

/*
  TH1* TSpectrum::Background (
   const     TH1* h,     : input 1-d histogram
   Int_t     niter = 20, : numberIterations. Rising niter makes the result smoother and lower
   Option_t* option = "" : Options, see below.
  )

  Options:

  - To set the direction parameter "BackIncreasingWindow". (Default: BackDecreasingWindow)
  - Order of clipping filter (default "BackOrder2"). Possible values= BackOrder4, ..6, ..8
  - "nosmoothing" : if selected, the background is not smoothed (Default: smoothed)
  - Width of smoothing window, (Default: BackSmoothing3). Possible: 3, ..., 15
  - "Compton"     : if selected, the Compton edges will be estimated.
  - "same"        : if selected, the background histogram is superimposed

  This function calculates the background spectrum in the input histogram h.
  The background is returned as a TH1* histogram.

  See root.cern.ch/doc/master/classTSpectrum.html#a56e15d1c36c7557b17d5f6d68fa315a7
*/

  if (hbkgnd) c1->Update();


  // Estimate the linear background using the very simple fit

  c1->cd (2);
  TF1* funLine = new TF1 ("funline", "pol1", 0, 1000);
  hspect->Fit ("funline", "QN");


  /***  For every found peak:                             ***/
  /***  - eliminate too weak peaks                        ***/
  /***  - use found (X,Y) as inits for further fitting    ***/

  par[0] = funLine->GetParameter (0);  			  // Estimate initial fit parameters
  par[1] = funLine->GetParameter (1);
  NPeaks = 0;
  Double_t* xpeaks = Analyser1dim->GetPositionX ();       // Returns X-positions of found peaks
  Double_t* ypeaks = Analyser1dim->GetPositionY ();

  for (Int_t p = 0; p < NFound; p++)
  {
    Double_t xp = xpeaks[p];
    Int_t   bin = hspect->GetXaxis()->FindBin (xp);
    Double_t yp = hspect->GetBinContent (bin);
    if (yp - sqrt(yp) < funLine->Eval(xp)) continue;      // Simple 'significance' test
    par [3*NPeaks + 2] = yp;
    par [3*NPeaks + 3] = xp;
    par [3*NPeaks + 4] = 3;
    NPeaks++;
  }
  printf ( "Found %d useful peaks to fit\n", NPeaks);

  printf ( "Now fitting: Be patient\n" );

  TF1* funFit = new TF1 ("funfit", funPeaks, 0, 1000, 2+3*NPeaks);
  funFit->SetParName (0, "Bg const");
  funFit->SetParName (1, "Bg lin"  );

  for (Int_t p = 0; p < NPeaks; p++) {
    funFit->SetParName (3*p+2 , Form ("P%d Norm" , p) );
    funFit->SetParName (3*p+3 , Form ("P%d Mean" , p) );
    funFit->SetParName (3*p+4 , Form ("P%d Sigma", p) );
  }

  // We may have more than the default 25 parameters

  TVirtualFitter::Fitter (h2, 10+3*NPeaks);
  funFit->SetParameters (par);
  funFit->SetNpx (10000);

  h2->Fit ("funfit");					// Now fitting.
}