const Int_t     kNumMaxPeaks   = 30;
const Int_t     kNumMaxPars    = 3*kNumMaxPeaks;
const Int_t     kNumMaxHalfWid = 30;         // for FWHM estimate

class TTest {

private:
  TH1D      *fHist;   // a Gaussian spectrum 
  TH1D      *fHsub;
  TSpectrum *fTSpec;  // a TSpectrum object
  TF1       *fFitFunc;

public:
  TCanvas   *c;
  virtual ~TTest(); // Destructor
  TTest(TCanvas*c);  // Constructor
  void FindAPeak(void);

};


TH1D      *fHback;  // estimated background spectrum

//______________________________________________________________________________
//
// Peak-fit function
//
Double_t
fpeaks(Double_t *x, Double_t *par) {
  Double_t xx = x[0];
  Int_t bin = fHback->GetXaxis()->FindBin(xx);
  Double_t result = fHback->GetBinContent(bin);
  for (Int_t p=0; p <kNumMaxPeaks; p++) {
    Double_t norm  = par[3*p+0];
    Double_t mean  = par[3*p+1];
    Double_t sigma = par[3*p+2]/2.35482;  // fPar[3*p+2]: FWHM
    result += norm*TMath::Gaus(x[0],mean,sigma);
  }
  return result;
}

//______________________________________________________________________________

TTest::TTest (TCanvas* c) : fHist(0) {
  TFile *f = new TFile("test.root","read");
  fHist = (TH1D*)f->Get("hH00");
  fHist->GetXaxis()->SetRangeUser(600.,800.);
  fHist->Draw();
  c->Update();
}

TTest::~TTest() {
  if (fHist)    delete fHist;
  if (fTSpec)   delete fTSpec;
  if (fFitFunc) delete fFitFunc;
}

void
TTest::FindAPeak(void) {
  fTSpec = new TSpectrum(1,1.0); // maxPos=1, resolution=1.0
  Int_t nPeaks = fTSpec->Search(fHist,2,"",0.10); //sigma=2,thresh.=0.10
  cout << "nPeaks = " << nPeaks << endl;
  fHback = (TH1D*)fTSpec->Background(fHist,20,"Compton,same");
  //fHsub = (TH1D*)fHist->Clone(fHist->GetName());
  //fHsub->Add(fHback,-1.);
  //fHsub->Draw("same");
  //c->Update(); // <--- leads to segmentation violation
# if 0
  fHist->Fit("gaus","+","same",682.,702.);
# else
  fFitFunc = new TF1("fit", fpeaks, 682.,702.,3);
  fFitFunc->SetParameter(0,385712);
  fFitFunc->SetParameter(1,690.);
  fFitFunc->SetParameter(2,9.12);
  //fFitFunc->SetNpx(1000); // Set number of points along X
  fFitFunc->SetLineColor(4); // Blue
  fFitFunc->SetLineWidth(3);
  fHist->Fit("fit","","SAME");
  fFitFunc->Print("all");
# endif // use Gaussian
}