/* Author: PengJi Ding (丁鹏基) this program is writed for background-removing with SNIP algorithm */
{
/* input spectrum data and draw the original spectrum histogram h1 */
FILE *fp=fopen("/home/dinghao/dinghao/Spectrum/Co60-1.txt","r"); // 所用的数据是Co-60的实验数据
Int_t nbins=8193;
Float_t m,n;
TH1F *h1=new TH1F("h1","Original spectrum",8193,0,8192);
TCanvas *Background = gROOT->GetListOfCanvases()->FindObject("Background");
if (!Background) Background =
new TCanvas("Background","Estimation of background with decreasing window",
10,10,800,500);
Background->Divide(1,2);
Background->cd(1);
for(Int_t i=0;i<nbins;i++) {
fscanf(fp,"%f %f",&(m[i]),&(n[i]));
h1->Fill(i,n[i]);
h1->Draw();
}
fclose(fp);
TFile f("/home/dinghao/dinghao/Spectrum/demo.root","recreate"); // save root file
h1->Write();
/* remove the background */
int i;
Double_t xmin = 0;
Double_t xmax = 8192;
Float_t *source = new float[nbins];
TH1F *back = new TH1F("back","Original spectrum and Background",nbins,xmin,xmax);
TH1F *d = new TH1F("d","",nbins,xmin,xmax);
TH1F *d1 = new TH1F("d","",nbins,xmin,xmax);
back=(TH1F*)f->Get("h1;1");
TH1F *h3 = h1->Clone();
Background->cd(2);
/* add searching peaks method and draw the result histogram h3 */
Float_t fPositionX[100];
Float_t fPositionY[100];
Int_t fNPeaks = 0;
Int_t i,nfound,bin;
Double_t a;
Float_t *dest = new float[nbins];
TSpectrum *s = new TSpectrum();
for (i = 0; i < nbins; i++) source[i]=back->GetBinContent(i + 1);
nfound = s->SearchHighRes(source, dest, nbins, 2.3, 10, kTRUE, 10, kTRUE, 3);
Float_t *xpeaks = s->GetPositionX();
for (i = 0; i < nfound; i++) {
a=xpeaks[i];
bin = 1+Int_t(a + 0.5);
fPositionX[i] = h3->GetBinCenter(bin);
fPositionY[i] = h3->GetBinContent(bin);
}
TPolyMarker * pm = (TPolyMarker*)h3->GetListOfFunctions()->FindObject("TPolyMarker");
if (pm) {
h3->GetListOfFunctions()->Remove(pm);
delete pm;
}
pm = new TPolyMarker(nfound, fPositionX, fPositionY);
h3->GetListOfFunctions()->Add(pm);
pm->SetMarkerStyle(23);
pm->SetMarkerColor(kRed);
pm->SetMarkerSize(1.3);
h3->Draw();
for (i = 0; i < nbins; i++) d->SetBinContent(i + 1,dest[i]);
d->SetLineColor(kMagenta);
d->Draw("SAME");
printf("Found %d candidate peaksn",nfound);
}