Hi all,
I am trying to deconvolve a gamma spectra with a linear system of changing response. I successfully implemented the Deconvolution function, but have a little trouble with the Unfolding function.
When filling resp[f][i] = hinr->GetBinContent(f+1,i+1);
, after several loops, the program Breaks:
" *** Break *** segmentation violation"
In addition, when I feel it with 0’s, ROOT starts the Unfolding Function but also breaks:
" *** Break *** segmentation violation"
Please let me know anyone has come through this and how to remedy is.
Thanks.
//ROOT STUFF~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
TSpectrum pfinder;
TFile *f = new TFile("Spectrum.root");
TH1F *hins = new TH1F("h1","my histogram",2220,0,3000);
TH2F *hinr = new TH2F("h2","Detector Response Function",50,0,11.1, 50, 0, 5500 );
//FILL IN OBSERVED SPECTRUM
for (int i = 0; i < bins2; i = i+1){
hins->Fill(energySpec[i], frequency[i]);
} //end "for" to print results
//FILL IN RESPONSE MATRIX
for (int f = 0; f < xbins ; f = f+1){ // f is for differen response functions
for (int i = 0; i < ybins; i = i+1){
hinr->Fill(energy[i], f, response_matrix[f][i]);
} //end "for" to print results
}
Int_t size = hins->GetXaxis()->GetNbins();
cout << size << endl;
Double_t *source;
source = new Double_t [size];
Double_t ** resp = new Double_t * [xbins];
Int_t i;
//FILL IN SOURCE MATRIX
for (i=0;i<size;i++) {
source[i] = hins->GetBinContent(i+1);
}
//FILL IN RESPONSE MATRIX
for (int f = 0; f < xbins ; f = f+1){ // f is for differen response functions
cout << f << endl;
for (int i = 0; i < ybins; i = i+1){
resp[f][i] = hinr->GetBinContent(f+1,i+1);
}
}
cout << "Commencing Deconvolution now" << endl;
pfinder.Unfolding(source, const_cast<const double**>(resp), xbins, ybins, 100,1,1);
cout << "Deconvolution Complete" << endl;
TH1F *hout = new TH1F("hdec1_out","hdec1_out",size,0,size);
for (i=0;i<nbins;i++) {
hout->Fill(i+1,source[i]);
}
gStyle->SetFrameFillColor(18);
TCanvas *c1 = new TCanvas("c1","Bacground estimator",600,800);
c1->SetFillColor(38);
c1->Divide(1,2);
c1->cd(1);
// draw original spectrum and response signal
hins->Draw("HIST");
hinr->SetLineColor(kRed);
hinr->Draw("same");
c1->SetLogy();
// draw deconvolution
c1->cd(2);
hout->Draw("HIST");
c1->cd();
//delete [] resp;
delete [] source;
//END OF ROOT STUFF~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~