24 #include "TSpectrum.h"
25 #include "TVirtualFitter.h"
28 Double_t
fpeaks(Double_t *
x, Double_t *par) {
29 Double_t result = par[0] + par[1]*x[0];
31 Double_t
norm = par[3*
p+2];
32 Double_t
mean = par[3*
p+3];
33 Double_t
sigma = par[3*
p+4];
34 result += norm*TMath::Gaus(x[0],mean,sigma);
40 TH1F *
h =
new TH1F(
"h",
"test",500,0,1000);
48 par[3*p+3] = 10+gRandom->Rndm()*980;
49 par[3*p+4] = 3+2*gRandom->Rndm();
51 TF1 *
f =
new TF1(
"f",
fpeaks,0,1000,2+3*npeaks);
53 f->SetParameters(par);
54 TCanvas *c1 =
new TCanvas(
"c1",
"c1",10,10,1000,900);
57 h->FillRandom(
"f",200000);
59 TH1F *
h2 = (TH1F*)h->Clone(
"h2");
61 TSpectrum *
s =
new TSpectrum(2*npeaks);
62 Int_t nfound = s->Search(h,2,
"",0.10);
63 printf(
"Found %d candidate peaks to fit\n",nfound);
65 TH1 *hb = s->Background(h,20,
"same");
71 TF1 *fline =
new TF1(
"fline",
"pol1",0,1000);
74 par[0] = fline->GetParameter(0);
75 par[1] = fline->GetParameter(1);
77 Double_t *xpeaks = s->GetPositionX();
78 for (p=0;p<nfound;p++) {
79 Double_t xp = xpeaks[
p];
80 Int_t bin = h->GetXaxis()->FindBin(xp);
81 Double_t yp = h->GetBinContent(bin);
82 if (yp-TMath::Sqrt(yp) < fline->Eval(xp))
continue;
88 printf(
"Found %d useful peaks to fit\n",npeaks);
89 printf(
"Now fitting: Be patient\n");
90 TF1 *
fit =
new TF1(
"fit",
fpeaks,0,1000,2+3*npeaks);
92 TVirtualFitter::Fitter(h2,10+3*npeaks);
93 fit->SetParameters(par);
117 Double_t result = par[0] + par[1]*x[0];
119 result += par[3]*TMath::Gaus(x[0],par[4],par[5]);
121 Double_t
norm = par[6 +
p*2];
122 Double_t
mean = par[4] + par[2]*(
p+1);
123 Double_t
sigma = par[8 +
p*2];
124 result += norm*TMath::Gaus(x[0],mean,sigma);
133 TH1F *
data = (TH1F*)(gSystem->FindObject(hName));
135 cout<<
"spmCalib Data Histogramm "<<hName<<
" not found"<<endl;
139 Double_t range(150.);
140 Double_t Double_t norm0(1.), slope0(-1./range);
141 Double_t peak_to_peak(10.);
142 Double_t pLoc(0.), pRMS(1.5);
143 Double_t par[6+maxPeaks*2];
147 par[2] = peak_to_peak;
151 for (
int p=0;
p< maxPeaks;
p++) {
152 par[6 +
p*3] = norm0;
153 par[7 +
p*3] = par[4] + par[2]*(
p+1);
156 TString scId =
"sc_"; scId += hName;
157 TF1 * scF = (TF1*) (gROOT->FindObject(scId));
if(scF)
delete scF;
158 scF =
new TF1(scId,&
sipmpeaks,0,range,6+maxPeaks*2);
159 TH1F * dataCl = (TH1F*)(data->Clone(
"dataCl"));
160 dataCl -> SetAxisRange(0., range);
161 Int_t binsToUse = dataCl->GetXaxis()->FindBin(range);
162 TCanvas *
peaks =
new TCanvas(
"peaks",
"peaks",10,10,1000,900);
163 peaks -> Divide(1,2);
166 TH1F * dataCl2 = (TH1F*)data->Clone(
"dataCl2");
169 TSpectrum *
s =
new TSpectrum(maxPeaks+1);
170 Int_t nfound = s->Search(dataCl, 2,
"", 0.10);
171 printf(
"TSpectrum: Found %d candidate peaks to fit\n",nfound);
174 TH1 *hb = s->Background(
h,20,
"same");
175 if (hb) peaks->Update();
179 TF1 *fline =
new TF1(
"fline",
"pol1",0,range);
180 dataCl -> Fit(
"fline",
"qn");
182 par[0] = fline->GetParameter(0);
183 par[1] = fline->GetParameter(1);
185 Double_t *xpeaks = s->GetPositionX();
186 Double_t spacing = 0.;
187 for (
p=0;
p<nfound;
p++) {
188 Double_t xp = xpeaks[
p];
189 Int_t bin = dataCl->GetXaxis()->FindBin(xp);
190 Double_t yp = dataCl->GetBinContent(bin);
191 cout<<
"Testing Peak "<<
p<<
" xp "<<xp<<
" bin "<<bin<<
" yP "<<yp<<
" BCKG "<<fline->Eval(xp)<<endl;
192 if (yp-TMath::Sqrt(yp) < fline->Eval(xp))
continue;
193 if(npeaks) spacing += xp - par[4 + (npeaks-1)+3];
194 par[3 + npeaks*3] = yp;
195 par[4 + npeaks*3] = xp;
196 par[5 + npeaks*3] = 3.;
199 printf(
"Found %d useful peaks to fit\n",npeaks);
200 printf(
"Now fitting: Be patient\n");
203 TVirtualFitter::Fitter(
h2,10+3*npeaks);
204 fit->SetParameters(par);
218 Double_t rcpar[2+3+2*nPeaks];
219 Double_t fcpar[2+3+2*nPeaks];
223 TString cSC =
"Single Cell Calibration (black - raw, red - fit) Run "; cSC +=
runId.Data();
224 scCalibDisplay =
new TCanvas(
"scCalibDisplay",cSC,200*nx_c,100*ny_c);
226 gStyle->SetOptStat(1110);
227 gStyle->SetOptFit(111);
228 Int_t binsToUse(200); Double_t slope0(-1./(Double_t))binsToUse); Double_t norm0(1.);
229 Double_t peak_to_peak(10.);
230 for (
int ich = 0; ich<
NCH; ich++){
231 TH1F *hrv = (TH1F*)rpeak[ich];
233 hrv -> SetAxisRange(0., binsToUse);
234 TH1F *fmv = (TH1F*)
fm[ich];
236 fmv -> SetAxisRange(0., binsToUse);
241 rcpar[2] = peak_to_peak;
246 for (
int p=1;
p<= nPeaks;
p++) {
247 rcpar[2*
p+4] = norm0;
252 fcpar[2] = peak_to_peak;
257 for (
int p=1;
p<= nPeaks;
p++) {
258 fcpar[2*
p+4] = norm0;
262 TString rvscId =
"rvsc_"; rvscId += ich;
263 TString fmscId =
"fmsc_"; fmscId += ich;
264 rvsc[ich] = (TF1*) (gROOT->FindObject(rvscId));
if(rvsc[ich])
delete rvsc[ich];
265 fmsc[ich] = (TF1*) (gROOT->FindObject(fmscId));
if(fmsc[ich])
delete fmsc[ich];
266 rvsc[ich] =
new TF1(rvscId,&scpeaks,0,500,5+2*nPeaks); fmsc[ich] =
new TF1(fmscId,&scpeaks,0,500,5+2*nPeaks);
267 rvsc[ich] ->
SetNpx(binsToUse); fmsc[ich] ->
SetNpx(binsToUse);
268 rvsc[ich] -> SetParameters(rcpar); fmsc[ich] -> SetParameters(fcpar);
270 TH1F * rh = (TH1F*)(hrv->Clone(
"rh")); TH1F * fh = (TH1F*)(fmv->Clone(
"fh"));
273 TF1 *fline =
new TF1(
"fline",
"pol1",0,binsToUse);
275 TSpectrum *
s =
new TSpectrum(2*nPeaks);
277 Int_t rfound = s->Search(hrv,1,
"new");
278 printf(
"Found %d candidate peaks in Raw hist ",rfound);
280 hrv -> Fit(
"fline",
"qn");
281 rcpar[0] = fline->GetParameter(0);
282 rcpar[1] = fline->GetParameter(1);
285 Float_t *xrpeaks = s->GetPositionX();
286 for (
p=0;
p<nfound;
p++) {
287 Float_t xp = xpeaks[
p];
288 Int_t bin = hrv->GetXaxis()->FindBin(xp);
289 Float_t yp =
h->GetBinContent(bin);
290 if (yp-TMath::Sqrt(yp) < fline->Eval(xp))
continue;
301 int_t ffound = s->Search(fmv,1,
"new");
302 printf(
"Found %d candidate peaks in Fit hist ",ffound);
304 fmv -> Fit(
"fline",
"qn");
305 fcpar[0] = fline->GetParameter(0);
306 fcpar[1] = fline->GetParameter(1);
310 Float_t *xpeaks = s->GetPositionX();
311 for (
p=0;
p<nfound;
p++) {
312 Float_t xp = xpeaks[
p];
313 Int_t bin =
h->GetXaxis()->FindBin(xp);
314 Float_t yp =
h->GetBinContent(bin);
315 if (yp-TMath::Sqrt(yp) < fline->Eval(xp))
continue;
330 scCalibDisplay ->
cd(ich+1); hrv->SetAxisRange(0., 200); hrv ->Draw(); fmv->Draw(
"same");
331 hrv -> Fit(rvsc[ich],
"rnQ"); fmv -> Fit(fmsc[ich],
"rnq");
332 rvsc[ich] ->
Draw(
"same"); fmsc[ich] ->
Draw(
"same");