// Richardson-Lucy deconvolution algorithm
//
// root.cern.ch/root/htmldoc/guides/spectrum/Spectrum.html#deconvolution---unfolding
// root.cern.ch/doc/master/DeconvolutionRL__wide_8C.html

// Example to illustrate deconvolution function (class TSpectrum).
//
// Authors: Miroslav Morhac, Olivier Couet


void DeconvolutionRL_wide () 
{
  Int_t i;
  const Int_t nbins = 256;
  Double_t  xmin  = 0;
  Double_t  xmax  = nbins;
  Double_t* source   = new Double_t [nbins];
  Double_t* response = new Double_t [nbins];
  gROOT->ForceStyle();

  TFile* f = new TFile ("TSpectrum.root");
  TH1F*  hdata = (TH1F*) f->Get ("decon3");
  hdata->SetTitle ("Deconvolution of closely positioned overlapping peaks using Richardson-Lucy deconvolution method");
  TH1F* hresponse = (TH1F*) f->Get ("decon_response_wide");

  for (i = 0; i < nbins; i++)  source  [i] = hdata    ->GetBinContent (i+1);
  for (i = 0; i < nbins; i++)  response[i] = hresponse->GetBinContent (i+1);
  hdata->SetMaximum (30000);
  hdata->Draw ("L");

  TSpectrum *Analyser1dim = new TSpectrum ();
  Analyser1dim->DeconvolutionRL (source, response, 256, 10000, 1, 1);

/*
  const char* TSpectrum::DeconvolutionRL (
     Double_t*       source,   // pointer to the vector of source   spectrum
     const Double_t* response, // pointer to the vector of response spectrum
     Int_t ssize,              // length of source and response spectra
     Int_t numberIterations,   // for details we refer to documentation
     Int_t numberRepetitions,  // for repeated boosted deconvolution
     Double_t boost            // boosting coefficient, [xi] --> [xi]^boost
  )

  This function calculates deconvolution from source spectrum according to response
  spectrum using Richardson-Lucy deconvolution algorithm. The result is placed in 
  the vector pointed by source pointer. On successful completion it returns 0. 
  On error it returns pointer to the string describing error. If desired after 
  every numberIterations one can apply boosting operation (exponential function 
  with exponent given by boost coefficient) and repeat it numberRepetitions times 
  (see Gold deconvolution).

  See: root.cern.ch/doc/v608/classTSpectrum.html#ac0a2c2b867d5e962f93cf4ddd458fe80
*/

  for (i = 0; i < nbins; i++) cout << source[i] << endl;

  TH1F* hresult = (TH1F*) hresponse->Clone ("hresult");
  hresult->Reset ();
  for (i = 0; i < nbins; i++)  hresult->SetBinContent (i+1, source[i]);
  hresult->SetLineColor (kRed);
  hresult->Draw ("SAME L");
}