本文整理汇总了C++中TSpectrum类的典型用法代码示例。如果您正苦于以下问题:C++ TSpectrum类的具体用法?C++ TSpectrum怎么用?C++ TSpectrum使用的例子?那么恭喜您, 这里精选的类代码示例或许可以为您提供帮助。
在下文中一共展示了TSpectrum类的20个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于我们的系统推荐出更棒的C++代码示例。
示例1: DeconvolutionRL_wide
void DeconvolutionRL_wide() {
Int_t i;
const Int_t nbins = 256;
Double_t xmin = 0;
Double_t xmax = nbins;
Double_t source[nbins];
Double_t response[nbins];
gROOT->ForceStyle();
TString dir = gROOT->GetTutorialsDir();
TString file = dir+"/spectrum/TSpectrum.root";
TFile *f = new TFile(file.Data());
h = (TH1F*) f->Get("decon3");
h->SetTitle("Deconvolution of closely positioned overlapping peaks using Richardson-Lucy deconvolution method");
d = (TH1F*) f->Get("decon_response_wide");
for (i = 0; i < nbins; i++) source[i]=h->GetBinContent(i + 1);
for (i = 0; i < nbins; i++) response[i]=d->GetBinContent(i + 1);
h->SetMaximum(30000);
h->Draw("L");
TSpectrum *s = new TSpectrum();
s->DeconvolutionRL(source,response,256,10000,1,1);
for (i = 0; i < nbins; i++) d->SetBinContent(i + 1,source[i]);
d->SetLineColor(kRed);
d->Draw("SAME L");
}
开发者ID:Y--,项目名称:root,代码行数:28,代码来源:DeconvolutionRL_wide.C
示例2: while
int e428ana1::MakePara(Double_t thr){
int npeak,nlap;
double threshol,thrmin,thrmax;
TSpectrum *spec;
for (int i=0; i<4; i++) {
for (int j=0; j<127; j++) {
npeak=0;nlap=0;
threshol=thr;
thrmax=1;thrmin=0;
while (npeak!=4 && nlap<100){
spec = new TSpectrum(10);
npeak = spec->Search(hdssd[i][j],2,"",threshol);
if (npeak >3) {
thrmin=threshol;
}else thrmax=threshol;
threshol=(thrmax-thrmin)/2;
nlap++;
// printf("Npeak %d %d %d %f\n",i,j,npeak,threshol);
}
float *xpeak=spec->GetPositionX();
float *ypeak=spec->GetPositionY();
printf("Npeak %d %d %d",i,j,npeak);
for (int i1=0; i1<npeak; i1++) {
printf(" %f %f",xpeak[i1],ypeak[i1]);
}
printf("\n");
}
}
return 0;
}
开发者ID:tdtrongiop,项目名称:E428,代码行数:30,代码来源:e428ana1.C
示例3: Deconvolution_wide_boost
void Deconvolution_wide_boost() {
Int_t i;
const Int_t nbins = 256;
Double_t xmin = 0;
Double_t xmax = nbins;
Double_t source[nbins];
Double_t response[nbins];
gROOT->ForceStyle();
TH1F *h = new TH1F("h","Deconvolution",nbins,xmin,xmax);
TH1F *d = new TH1F("d","",nbins,xmin,xmax);
TString dir = gROOT->GetTutorialDir();
TString file = dir+"/spectrum/TSpectrum.root";
TFile *f = new TFile(file.Data());
h = (TH1F*) f->Get("decon3");
h->SetTitle("Deconvolution of closely positioned overlapping peaks using boosted Gold deconvolution method");
d = (TH1F*) f->Get("decon_response_wide");
for (i = 0; i < nbins; i++) source[i]=h->GetBinContent(i + 1);
for (i = 0; i < nbins; i++) response[i]=d->GetBinContent(i + 1);
h->SetMaximum(200000);
h->Draw("L");
TSpectrum *s = new TSpectrum();
s->Deconvolution(source,response,256,200,50,1.2);
for (i = 0; i < nbins; i++) d->SetBinContent(i + 1,source[i]);
d->SetLineColor(kRed);
d->Draw("SAME L");
}
开发者ID:davidlt,项目名称:root,代码行数:31,代码来源:Deconvolution_wide_boost.C
示例4: autogain152
TGraph* autogain152(TH1 *hist) {
hist->GetXaxis()->SetRangeUser(200.,16000.);
TSpectrum *s = new TSpectrum();
Int_t nfound = s->Search(hist,6,"",0.08); //This will be dependent on the source used.
printf("Found %d candidate peaks to fit\n",nfound);
if(nfound > 6)
nfound = 6;
std::vector<float> vec;
for(int x=0;x<nfound;x++)
vec.push_back(s->GetPositionX()[x]);
std::sort(vec.begin(),vec.end());
Float_t energies[] = {121.7830, 244.6920, 344.276, 778.903, 964.131, 1408.011};
TGraph* slopefit = new TGraph(nfound, &(vec[0]), energies);
printf("Now fitting: Be patient\n");
slopefit->Fit("pol1");
if(slopefit->GetFunction("pol1")->GetChisquare() > 10.) {
slopefit->RemovePoint(slopefit->GetN()-1);
slopefit->Fit("pol1");
}
TChannel *chan = 0;
slopefit->Draw("AC*");
return slopefit;
}
开发者ID:atlaffoley,项目名称:GRSISort,代码行数:29,代码来源:autoeffic.C
示例5: peaks
void peaks(Int_t np=10) {
npeaks = TMath::Abs(np);
TH1F *h = new TH1F("h","test",500,0,1000);
//generate n peaks at random
Double_t par[3000];
par[0] = 0.8;
par[1] = -0.6/1000;
Int_t p;
for (p=0;p<npeaks;p++) {
par[3*p+2] = 1;
par[3*p+3] = 10+gRandom->Rndm()*980;
par[3*p+4] = 3+2*gRandom->Rndm();
}
TF1 *f = new TF1("f",fpeaks,0,1000,2+3*npeaks);
f->SetNpx(1000);
f->SetParameters(par);
TCanvas *c1 = new TCanvas("c1","c1",10,10,1000,900);
c1->Divide(1,2);
c1->cd(1);
h->FillRandom("f",200000);
h->Draw();
TH1F *h2 = (TH1F*)h->Clone("h2");
//Use TSpectrum to find the peak candidates
TSpectrum *s = new TSpectrum(2*npeaks);
Int_t nfound = s->Search(h,2,"",0.10);
printf("Found %d candidate peaks to fit\n",nfound);
//Estimate background using TSpectrum::Background
TH1 *hb = s->Background(h,20,"same");
if (hb) c1->Update();
if (np <0) return;
//estimate linear background using a fitting method
c1->cd(2);
TF1 *fline = new TF1("fline","pol1",0,1000);
h->Fit("fline","qn");
//Loop on all found peaks. Eliminate peaks at the background level
par[0] = fline->GetParameter(0);
par[1] = fline->GetParameter(1);
npeaks = 0;
Double_t *xpeaks = s->GetPositionX();
for (p=0;p<nfound;p++) {
Double_t xp = xpeaks[p];
Int_t bin = h->GetXaxis()->FindBin(xp);
Double_t yp = h->GetBinContent(bin);
if (yp-TMath::Sqrt(yp) < fline->Eval(xp)) continue;
par[3*npeaks+2] = yp;
par[3*npeaks+3] = xp;
par[3*npeaks+4] = 3;
npeaks++;
}
printf("Found %d useful peaks to fit\n",npeaks);
printf("Now fitting: Be patient\n");
TF1 *fit = new TF1("fit",fpeaks,0,1000,2+3*npeaks);
//we may have more than the default 25 parameters
TVirtualFitter::Fitter(h2,10+3*npeaks);
fit->SetParameters(par);
fit->SetNpx(1000);
h2->Fit("fit");
}
开发者ID:Y--,项目名称:root,代码行数:59,代码来源:peaks.C
示例6: Background_order
void Background_order() {
Int_t i;
const Int_t nbins = 4096;
Double_t xmin = 0;
Double_t xmax = 4096;
Double_t source[nbins];
gROOT->ForceStyle();
TH1F *d1 = new TH1F("d1","",nbins,xmin,xmax);
TH1F *d2 = new TH1F("d2","",nbins,xmin,xmax);
TH1F *d3 = new TH1F("d3","",nbins,xmin,xmax);
TH1F *d4 = new TH1F("d4","",nbins,xmin,xmax);
TString dir = gROOT->GetTutorialsDir();
TString file = dir+"/spectrum/TSpectrum.root";
TFile *f = new TFile(file.Data());
TH1F *back = (TH1F*) f->Get("back2");
back->SetTitle("Influence of clipping filter difference order on the estimated background");
back->SetAxisRange(1220,1460);
back->SetMaximum(3000);
back->Draw("L");
TSpectrum *s = new TSpectrum();
for (i = 0; i < nbins; i++) source[i]=back->GetBinContent(i + 1);
s->Background(source,nbins,40,TSpectrum::kBackDecreasingWindow,
TSpectrum::kBackOrder2,kFALSE,
TSpectrum::kBackSmoothing3,kFALSE);
for (i = 0; i < nbins; i++) d1->SetBinContent(i + 1,source[i]);
d1->SetLineColor(kRed);
d1->Draw("SAME L");
for (i = 0; i < nbins; i++) source[i]=back->GetBinContent(i + 1);
s->Background(source,nbins,40,TSpectrum::kBackDecreasingWindow,
TSpectrum::kBackOrder4,kFALSE,
TSpectrum::kBackSmoothing3,kFALSE);
for (i = 0; i < nbins; i++) d2->SetBinContent(i + 1,source[i]);
d2->SetLineColor(kBlue);
d2->Draw("SAME L");
for (i = 0; i < nbins; i++) source[i]=back->GetBinContent(i + 1);
s->Background(source,nbins,40,TSpectrum::kBackDecreasingWindow,
TSpectrum::kBackOrder6,kFALSE,
TSpectrum::kBackSmoothing3,kFALSE);
for (i = 0; i < nbins; i++) d3->SetBinContent(i + 1,source[i]);
d3->SetLineColor(kGreen);
d3->Draw("SAME L");
for (i = 0; i < nbins; i++) source[i]=back->GetBinContent(i + 1);
s->Background(source,nbins,40,TSpectrum::kBackDecreasingWindow,
TSpectrum::kBackOrder8,kFALSE,
TSpectrum::kBackSmoothing3,kFALSE);
for (i = 0; i < nbins; i++) d4->SetBinContent(i + 1,source[i]);
d4->SetLineColor(kMagenta);
d4->Draw("SAME L");
}
开发者ID:Y--,项目名称:root,代码行数:56,代码来源:Background_order.C
示例7: autogain
TGraph* autogain(TH1 *hist,TNucleus *nuc) { //Display The fits on a TPad
if(!hist || !nuc)
return 0;
nuc->SetSourceData();
if(nuc->GetA() == 152) {
return autogain152(hist);
}
// Search
hist->GetXaxis()->SetRangeUser(200.,16000.);
TSpectrum *s = new TSpectrum();
Int_t nfound = s->Search(hist,6,"",0.1); //This will be dependent on the source used.
printf("Found %d candidate peaks to fit\n",nfound);
// Match
nuc->TransitionList.Sort();
std::vector<float> engvec;
TIter iter(&(nuc->TransitionList));
TObject* obj;
while(obj = iter.Next()) {
if(!obj->InheritsFrom("TGRSITransition"))
continue;
TGRSITransition *tran = (TGRSITransition*)obj;
engvec.push_back(static_cast<float>(tran->energy));
if(engvec.size() == nfound)
break;
}
if(nfound != engvec.size())
return 0;
Float_t *posPeaks = s->GetPositionX();
Float_t *energies = &(engvec[0]);
for(int x=0;x<nfound;x++) {
printf("posPeaks[%i] = %f\t\tenrgies[%i] = %f\n",x,posPeaks[x],x,energies[x]);
}
TGraph *slopefit = new TGraph(nfound,posPeaks,energies );
printf("Now fitting: Be patient\n");
slopefit->Fit("pol1");
slopefit->Draw("AC*");
return slopefit;
}
开发者ID:atlaffoley,项目名称:GRSISort,代码行数:52,代码来源:autoeffic.C
示例8: Fit
//______________________________________________________________________________
void Fit(Int_t run)
{
// Perform fit.
Char_t tmp[256];
// delete old function
if (gFitFunc) delete gFitFunc;
// create fitting function
if (gFitFunc) delete gFitFunc;
sprintf(tmp, "fGauss_%d", run);
gFitFunc = new TF1(tmp, "expo(0)+gaus(2)");
gFitFunc->SetLineColor(2);
// estimate peak position
TSpectrum s;
s.Search(gH, 10, "goff nobackground", 0.05);
Double_t peak = TMath::MaxElement(s.GetNPeaks(), s.GetPositionX());
// prepare fitting function
gFitFunc->SetRange(gMin, gMax);
gFitFunc->SetParameter(2, gH->GetXaxis()->FindBin(peak));
gFitFunc->SetParLimits(2, 0, 100000);
gFitFunc->SetParameter(3, peak);
gFitFunc->SetParameter(4, 20);
gFitFunc->SetParLimits(4, 18, 100);
for (Int_t i = 0; i < 10; i++)
if (!gH->Fit(gFitFunc, "RBQ0")) break;
// get peak
peak = gFitFunc->GetParameter(3);
// indicator line
gLine->SetX1(peak);
gLine->SetX2(peak);
gLine->SetY1(0);
gLine->SetY2(gH->GetMaximum());
// draw
gCFit->cd();
gH->GetXaxis()->SetRangeUser(0, 2000);
gH->Draw();
gFitFunc->Draw("same");
gLine->Draw("same");
// fill overview histogram
gHOverview->SetBinContent(run+1, peak);
gHOverview->SetBinError(run+1, 0.0001);
}
开发者ID:A2-Collaboration,项目名称:acqu,代码行数:52,代码来源:PIDEnergy.C
示例9: getValByKeyword
//Given an input histogram and TSpectrum returns a numeric value based on the input keyword; supported keywords are "AMPLITUDE,MEAN,PEAK,SIGMA"
float AnalyzeResponseUniformity::getValByKeyword(string strInputKeyword, shared_ptr<TH1F> hInput, TSpectrum &specInput){
//Try to automatically assign a value
if (0 == strInputKeyword.compare("AMPLITUDE") ) { //Case: Histo Amplitude
return hInput->GetBinContent( hInput->GetMaximumBin() );
} //End Case: Histo Amplitude
else if (0 == strInputKeyword.compare("MEAN") ) { //Case: Histo Mean
return hInput->GetMean();
} //End Case: Histo Mean
else if ( 0 == strInputKeyword.compare("PEAK") ){ //Case: Histo Peak
Double_t *dPeakPos = specInput.GetPositionX();
return dPeakPos[0];
} //End Case: Histo Peak
else if (0 == strInputKeyword.compare("SIGMA") ) { //Case: Histo RMS
return hInput->GetRMS();
} //End Case: Histo RMS
else{ //Case: manual assignment
printClassMethodMsg("AnalyzeResponseUniformity","getValByKeyword","Error! Input Keyword Not Recognized");
printClassMethodMsg("AnalyzeResponseUniformity","getValByKeyword", ("\tGiven: " + strInputKeyword ).c_str() );
printClassMethodMsg("AnalyzeResponseUniformity","getValByKeyword","\tRecognized Keywords:\n");
for (int i=0; i < vec_strSupportedKeywords.size(); ++i) {
printClassMethodMsg("AnalyzeResponseUniformity","getValByKeyword", vec_strSupportedKeywords[i].c_str() );
}
printClassMethodMsg("AnalyzeResponseUniformity","getValByKeyword","\tUndefined Behavior May Occur");
return -1e12;
} //End Case: manual assignment
} //End AnalyzeResponseUniformity::getValByKeyword()
开发者ID:jhallan,项目名称:CMS_GEM_Analysis_Framework,代码行数:32,代码来源:AnalyzeResponseUniformity.cpp
示例10: write_spectrum
void MSSM_spectrum_plotter::write_spectrum(const TSpectrum& spectrum, std::ofstream& filestr) const
{
for (std::size_t s = 0; s < spectrum.size(); ++s) {
if (!filestr.good()) {
ERROR("MSSM_spectrum_plotter::write_spectrum: "
"file stream is corrupted");
break;
}
const std::string& name = spectrum[s].name;
const std::string& latex_name = spectrum[s].latex_name;
const std::valarray<double>& masses = spectrum[s].masses;
const std::size_t multiplicity = masses.size();
filestr << std::left << "# " << name << '\n';
for (std::size_t i = 0; i < multiplicity; ++i) {
std::string lname("$" + latex_name + "$");
std::stringstream lname_with_index;
lname_with_index << "$" << latex_name;
if (multiplicity > 1)
lname_with_index << "_{" << (i+1) << "}";
lname_with_index << "$";
filestr << std::left << std::setw(width) << s
<< std::left << std::setw(width) << masses[i]
<< std::left << std::setw(width) << name
<< std::left << std::setw(2*width) << lname
<< std::left << std::setw(2*width) << lname_with_index.str()
<< '\n';
}
}
}
开发者ID:dylan-harries,项目名称:CNE6SSM-Spectrum,代码行数:33,代码来源:MSSM_utilities.cpp
示例11: FindPeak
void FindPeak(TH1 *hm, int * i, char * namefile){
int np =5, p, max = 0;
int npeaks = TMath::Abs(np);
double par[3000];
par[0] = 0.8;
par[1] = -0.6/1000;
for (p=0;p<npeaks;p++) {
par[3*p+2] = 1;
par[3*p+3] = 10+gRandom->Rndm()*980;
par[3*p+4] = 3+2*gRandom->Rndm();
}
TSpectrum *s = new TSpectrum(2*npeaks,1);
int nfound = s->Search(hm,2,"",0.10);
printf("Found %d candidate peaks to fit\n",nfound);
TH1 *hb = s->Background(hm,20,"same");
if (np <0) return;
// loope over peaks
TF1 *fline = new TF1("fline","pol1",0,1000);
hm->Fit("fline","qn");
par[0] = fline->GetParameter(0);
par[1] = fline->GetParameter(1);
npeaks = 0;
float *xpeaks = s->GetPositionX();
for (p=0;p<nfound;p++) {
float xp = xpeaks[p];
int bin = hm->GetXaxis()->FindBin(xp);
float yp = hm->GetBinContent(bin);
if (yp-TMath::Sqrt(yp) < fline->Eval(xp)) continue;
par[3*npeaks+2] = yp;
par[3*npeaks+3] = xp;
par[3*npeaks+4] = 3;
npeaks++;
}
printf("Found %d useful peaks to fit\n",npeaks);
printf("Now fitting: Be patient\n");
if (max < npeaks) max = npeaks;
TF1 *fit = new TF1("fit",fpeaks,0,1000,2+3*npeaks);
TVirtualFitter::Fitter(hm,10+3*npeaks);
fit->SetParameters(par);
fit->SetNpx(1000);
hm->Fit("fit");
}
开发者ID:bt3gl,项目名称:Calculating_the_Diameter_of_Sun,代码行数:46,代码来源:DataRedux.C
示例12: Load
//// This version needs thought
void SiCalibrator::FindPeaks() {
Load();
int iterate = 25;
int nbins;
TSpectrum* spec = new TSpectrum(10);
avesigma = 0;
int nsigma = 0;
for (int src=0; src < CalData.size(); src++) {
if (CalData[src].hSource != 0) {
for (int ch=0; ch<CalData[src].sourcedata.size(); ch++) {
TH1D* hbi = CalData[src].hSource->ProjectionY("hb",ch+1,ch+1);
nbins = hbi->GetNbinsX();
TH1D* hbk1 = (TH1D*) spec->Background(hbi,iterate);
hbi->Add(hbk1,-1);
//Estimate parameters
double lim = 5;
double max = 0;
double peak = 0;
int i = nbins - 1;
while (hbi->GetBinContent(i) < lim && i > 0) {
max = i;
i--;
}
double sigma = (max/CalibSource::sourcelist[src].betas.back().E)*2; //2 keV
if (sigma > 1) {
avesigma += sigma;
nsigma++;
}
gErrorIgnoreLevel = kError;
int npeaks = spec->Search(hbi,sigma,"nodraw ",0.001);
float* adc = spec->GetPositionX();
float* amp = spec->GetPositionY();
for (int i=0; i<npeaks; i++) {
CalData[src].sourcedata[ch].ADC.push_back(adc[i]);
CalData[src].sourcedata[ch].Amp.push_back(amp[i]);
}
}
}
}
avesigma = avesigma/nsigma;
}
开发者ID:losalamos,项目名称:UCNB_Analyzer,代码行数:42,代码来源:SiCalibrator.cpp
示例13: DiffPeak
/*
Finds peak in timing difference spectrum
*/
double DiffPeak(TTree *&tree, TString pidcut = "PID>1000")
{
TSpectrum *spec = new TSpectrum(10); // peak finder
TH1D *diff = new TH1D("diff","diff",200,-5,5);
tree->Project("diff","diff",pidcut.Data());
spec->Search(diff,2,"nodraw",0.7);
double diffpeak = spec->GetPositionX()[0];
// testing
//TCanvas *c_tmp = new TCanvas();
//diff->Draw();
//c_tmp->WaitPrimitive();
//delete c_tmp;
delete diff;
delete spec;
return diffpeak;
} // end DiffPeak()
开发者ID:hyperbolee,项目名称:prtdirc,代码行数:24,代码来源:analysistools.C
示例14: draw_ToF_Calib_Tower
void draw_ToF_Calib_Tower() {
ofstream fout("ToF_Calib_Tower.txt");
TFile* f = new TFile("/phenix/plhf/zji/taxi/Run13pp510MinBias/8328/data/total.root");
TH2I* h_hitmap_tof_energy_PbSc = (TH2I*)f->Get("hitmap_tof_energy_PbSc");
TH2I* h_hitmap_tof_energy_PbGl = (TH2I*)f->Get("hitmap_tof_energy_PbGl");
TH2I* h_hitmap_tof_tower = (TH2I*)f->Get("hitmap_tof_tower");
TCanvas* c = new TCanvas("c", "Canvas", 1800, 600);
gStyle->SetOptStat(1);
gStyle->SetOptFit(1);
c->Divide(3,1);
c->cd(1);
h_hitmap_tof_energy_PbSc->Draw("colz");
c->cd(2);
h_hitmap_tof_energy_PbGl->Draw("colz");
c->cd(3);
fout << "towerid ToF peak" << endl;
for(Int_t id=0; id<25000; id++) {
TH1D* hp_hitmap_tof_tower = (TH1D*)h_hitmap_tof_tower->ProjectionY("_py",id+1,id+1);
TSpectrum* peak = new TSpectrum();
Int_t nfound = peak->Search(hp_hitmap_tof_tower,2.,"nodraw");
if(nfound) {
Float_t peakX = *peak->GetPositionX();
hp_hitmap_tof_tower->Fit("gaus","Q0","",peakX-3.,peakX+3.);
TF1* f_hitmap_tof_tower = hp_hitmap_tof_tower->GetFunction("gaus");
Double_t mean_hitmap_tof_tower = f_hitmap_tof_tower->GetParameter(1);
fout << id << " " << mean_hitmap_tof_tower << endl;
f_hitmap_tof_tower->Delete();
}
else {
fout << id << " " << "-100." << endl;
}
delete peak;
hp_hitmap_tof_tower->Delete();
}
fout.close();
TH1D* hp_hitmap_tof_tower = (TH1D*)h_hitmap_tof_tower->ProjectionY("_py",100,100);
TSpectrum* peak = new TSpectrum();
peak->Search(hp_hitmap_tof_tower,2.,"nodraw");
Float_t peakX = *peak->GetPositionX();
hp_hitmap_tof_tower->Fit("gaus","","",peakX-3.,peakX+3.);
TF1* f_hitmap_tof_tower = hp_hitmap_tof_tower->GetFunction("gaus");
Double_t mean_hitmap_tof_tower = f_hitmap_tof_tower->GetParameter(1);
cout << "peakX=" << peakX << endl;
cout << "ToF mean=" << mean_hitmap_tof_tower << endl;
c->Print("ToF_Calib_Tower.pdf");
}
开发者ID:nfeege,项目名称:phenix-directphotons-pp,代码行数:54,代码来源:draw_ToF_Calib_Tower.C
示例15: peaks
//________________________________________________________________________________
void peaks(TH1* h, Int_t np=10, Int_t ij=0) {
if (! h) return;
npeaks = TMath::Abs(np);
if (! c1) c1 = new TCanvas();
else c1->Clear();
if (c2 && ij > 0) {c2->cd(ij); h->Draw(); c2->Update();}
c1->Divide(1,2);
c1->cd(1);
h->Draw();
TH1F *h2 = (TH1F*)h->Clone("h2");
//Use TSpectrum to find the peak candidates
TSpectrum *s = new TSpectrum(2*npeaks);
Int_t nfound = s->Search(h,5,"",0.05);
printf("Found %d candidate peaks to fit\n",nfound);
if (! nfound) return;
//Estimate background using TSpectrum::Background
TH1 *hb = s->Background(h,20,"same");
hb->Draw("same");
c1->Update();
if (c2 && ij > 0) {c2->cd(ij); h->Draw(); c2->Update();}
if (np <0) return;
//estimate linear background using a fitting method
c1->cd(2);
TF1 *fline = new TF1("fline","pol1",0,1000);
fline->FixParameter(1,0.);
h->Fit("fline","qnlw");
if (c2 && ij > 0) {c2->cd(ij+1); h->Draw(); c2->Update(); c1->cd(2);}
//Loop on all found peaks. Eliminate peaks at the background level
Double_t par[3000];
par[0] = fline->GetParameter(0);
par[1] = fline->GetParameter(1);
npeaks = 0;
Float_t *xpeaks = s->GetPositionX();
Float_t ymax = 0;
for (Int_t p=0;p<nfound;p++) {
Float_t xp = xpeaks[p];
Int_t bin = h->GetXaxis()->FindBin(xp);
Float_t yp = h->GetBinContent(bin);
if (yp-3*TMath::Sqrt(yp) < fline->Eval(xp)) continue;
par[3*npeaks+2] = yp;
if (yp > ymax) ymax = yp;
par[3*npeaks+3] = xp;
par[3*npeaks+4] = 3;
npeaks++;
}
cout << "Found " << npeaks << " useful peaks to fit" << endl;
Int_t NP = 0;
Int_t nbad = 0;
TString LineH(" {\""); LineH += h->GetName(); LineH += "\"";
TString Line("");
struct ParErr_t {Double_t par, err;};
ParErr_t parErr[10];
TF1 *fit = 0;
if (ymax > 2*par[0]) {
cout << "Now fitting: Be patient" << endl;
fit = new TF1("fit",fpeaks,0,1000,2+3*npeaks);
TVirtualFitter::Fitter(h2,10+3*npeaks); //we may have more than the default 25 parameters
fit->SetParameter(0,par[0]);
fit->FixParameter(1,0.);
for (Int_t p = 0; p < npeaks; p++) {
fit->SetParName(3*p+2,Form("A%i",p));
fit->SetParLimits(3*p+2,0,1e6);
fit->SetParameter(3*p+2,par[3*p+2]);
fit->SetParName(3*p+3,Form("#mu%i",p));
fit->SetParLimits(3*p+3,TMath::Max(0.,par[3*p+3]-2), TMath::Min(240.,par[3*p+3]+2));
fit->SetParameter(3*p+3,par[3*p+3]);
fit->SetParName(3*p+4,Form("#sigma%i",p));
fit->SetParLimits(3*p+4,0.01,20);
fit->SetParameter(3*p+4,par[3*p+4]);
}
fit->SetNpx(1000);
h2->SetStats(0);
h2->Fit("fit","l");
if (c2 && ij > 0) {c2->cd(ij); h2->Draw("same"); c2->Update(); c1->cd(2);}
fit->GetParameters(par);
for (Int_t p = 0; p<np;p++) {
if (p < npeaks && par[3*p+2] > 5*fit->GetParError(3*p+2) &&
par[3*p+2] > par[0]) {
if (TMath::Abs(par[3*p+4]) > 5.0) nbad++;
// Line += Form(",%f,%f,%7.2f,%5.2f",par[3*p+2],fit->GetParError(3*p+2),par[3*p+3],TMath::Abs(par[3*p+4]));
parErr[NP].par = par[3*p+3];
parErr[NP].err = TMath::Abs(par[3*p+4]);
for (Int_t i = 0; i < NP; i++) {
if (parErr[i].par > parErr[NP].par) {
ParErr_t temp = parErr[i];
parErr[i] = parErr[NP];
parErr[NP] = temp;
}
}
NP++;
}
}
}
for (Int_t p = 0; p < np; p++) {
if (p < NP) Line += Form(",%7.2f,%5.2f",parErr[p].par,parErr[p].err);
else Line += ",0,0";
}
Line += "},";
//.........这里部分代码省略.........
开发者ID:star-bnl,项目名称:star-svt,代码行数:101,代码来源:FindPeaks.C
示例16: fit_double_gaussian
void fit_double_gaussian(TH1F*& hrsp)
{
if (0==hrsp) {
cout<<"ERROR: Empty pointer to fit_double_gaussian()"<<endl;return;
}
string histname = hrsp->GetName();
double mean = hrsp->GetMean();
double rms = hrsp->GetRMS();
TSpectrum *spec = new TSpectrum(10);
int nfound = spec->Search(hrsp,1,"new");
if(nfound!=2)return;
double *xpeaks = spec->GetPositionX();
double peak1 = xpeaks[0];
double bin1 = hrsp->FindBin(peak1);
double norm1 = hrsp->GetBinContent(bin1);
double sigma1 = 0.2;
double peak2 = xpeaks[1];
double bin2 = hrsp->FindBin(peak2);
double norm2 = hrsp->GetBinContent(bin2);
double sigma2 = 0.3;
cout << " Mean : " << mean << " \t RMS : " << rms << endl;
cout << " norm1 : " << norm1 << " \t norm2 : " << norm2 << endl;
cout << " peak1 : " << peak1 << " \t sig1 : " << sigma1 << endl;
cout << " peak2 : " << peak2 << " \t sig2 : " << sigma2 << endl;
double fitrange_min = 0.2;
double fitrange_max = 1.7;
TF1* fitfnc(0); int fitstatus(-1);
TF1 *fitg1(0), *fitg2(0);
fitfnc = new TF1("fdgaus","gaus(0)+gaus(3)",fitrange_min,fitrange_max);
fitfnc->SetLineColor(1);
fitfnc->SetLineStyle(2);
fitfnc->SetLineWidth(2);
fitfnc->SetParNames("N_{1}", "#mu_{1}", "#sigma_{1}",
"N_{2}", "#mu_{2}", "#sigma_{2}");
fitfnc->SetParameters(norm1, peak1, sigma1,
norm2, peak2, sigma2);
//fitfnc->SetParLimits(0,0.0,2.0*norm1);
//fitfnc->SetParLimits(1,peak1-3*sigma1,peak1+3*sigma1);
//fitfnc->SetParLimits(2,0.1,1.0);
//fitfnc->SetParLimits(3,0.0,1.5*norm2);
//fitfnc->SetParLimits(4,peak2-3*sigma2,peak2+3*sigma2);
//fitfnc->SetParLimits(5,0.01,1.0);
fitstatus = hrsp->Fit(fitfnc,"RQ");
// if (0!=fitstatus){
// fitfnc->SetParLimits(4,0.2,1.7);
// fitfnc->SetParLimits(5,2.0,10.0);
// //cout <<" Not able to Fit this pt bin " << hrsp->GetName() << endl;
// }
fitstatus = hrsp->Fit(fitfnc,"RQ");
hrsp->SetMaximum(norm1+0.2*norm1);
fitg1 = new TF1("fg1","gaus(0)",fitrange_min,fitrange_max);
fitg1->SetParameters(fitfnc->GetParameter(0),
fitfnc->GetParameter(1),
fitfnc->GetParameter(2));
fitg1->SetLineColor(2);
fitg1->SetLineStyle(2);
hrsp->GetListOfFunctions()->Add(fitg1);
fitg2 = new TF1("fg2","gaus(0)",fitrange_min,fitrange_max);
fitg2->SetParameters(fitfnc->GetParameter(3),
fitfnc->GetParameter(4),
fitfnc->GetParameter(5));
fitg2->SetLineColor(4);
fitg2->SetLineStyle(4);
hrsp->GetListOfFunctions()->Add(fitg2);
if(hrsp->GetFunction("fdgaus")==0){
cout << "No function recorded in histogram " << hrsp->GetName() << endl;
}
if (0!=fitstatus){
cout<<"fit_double_gaussian() to "<<hrsp->GetName()
<<" failed. Fitstatus: "<<fitstatus
<<" - FNC deleted."<<endl;
hrsp->GetListOfFunctions()->Delete();
}
}
开发者ID:pawannetrakanti,项目名称:UserCode,代码行数:91,代码来源:CaloL3.C
示例17: I2GFmainLoop
I2GFvalues I2GFmainLoop(TH1F *htemp, int N_iter, float N_sigma_range, bool ShowFit)
//Arguments: (histo to be fit, N iterations to find peak using gaus fit, fit range param., do or do not plot fit on histo)
{
I2GFvalues myI2GFvalues;
//Set initial values...(in case fit fails)
myI2GFvalues.rchi2 = -100;
myI2GFvalues.mean = -100;
myI2GFvalues.mean_err = -100;
myI2GFvalues.sigma = -100;
myI2GFvalues.sigma_err = -100;
TSpectrum *s = new TSpectrum(); //TSpectrum(1,1)->Argument: (Number of peaks to find, Distance to neighboring peak: "1"-->3sigma)
Int_t NPeaks;
Float_t *Peak; //TSpectrum *s = new TSpectrum(); --> No warning message
Float_t *PeakAmp;
float peak_pos = 0;
float peak_pos_amp = 0; //Initial value, assuming positive going peaks
int peak_pos_bin = 0;
int binMaxCnt = 0;
float binMaxCnt_value = 0;
int binMaxCnt_counts = 0;
int Nbins = 0;
Int_t zero_value_bin = 0;
float low_limit = 0;
float high_limit = 0;
float peak_pos_count = 0;
float zero_bin_value = 0;
float max_bin_value = 0;
TF1 *func;
TF1 *func1;
TF1 *func2;
TF1 *func3;
float Chi2;
int NDF = 1;
float f_RChi2 = 1;
float f_const = 1;
float f_mean = 1;
float f_mean_err = 1;
float f_sigma = 1;
float f_sigma_err = 1;
float peak = 1;
float f_const2 = 1;
float f_mean2 = 1;
float f_mean_err2 = 1;
float f_sigma2 = 1;
float f_sigma_err2 = 1;
//---------Basic Histo Peak Finding Parameters----------
binMaxCnt = htemp->GetMaximumBin(); //Finds bin w/ max counts (can be dubious, so use TSpectrum)
binMaxCnt_value = (Double_t) htemp->GetXaxis()->GetBinCenter(binMaxCnt); //if the bin number is known and the bin value (in x-axis units) is wanted
binMaxCnt_counts = (Int_t) htemp->GetBinContent(binMaxCnt); //Finds counts within particular bin
//---------TSpectrum Peak Finding Parameters--------
if (ShowFit) NPeaks = s->Search(htemp, 2, "goff", 0.5); //opens a canvas (one time in a loop), even with: s->Search(htemp, 2, "nodraw", 0.9);
//else NPeaks = s->Search(htemp, 2, "", 0.5); //s->Search(htemp, 2, "nodraw", 0.9);
//Npeaks = s->GetNPeaks(); //If using this, need pointer in declaration above: Int_t *NPeaks
//N_peaks = Npeaks[0];
Peak = s->GetPositionX();
PeakAmp = s->GetPositionY();
for (int i=0; i<NPeaks; i++)
{
if (peak_pos_amp < PeakAmp[i])
{
peak_pos_amp = PeakAmp[i]; //TSpectrum finds peak counts
peak_pos = Peak[i]; //TSpectrum finds pos. of peak in x-axis units
}
}
peak_pos_bin = htemp->GetXaxis()->FindBin(peak_pos); //if the bin value (in x-axis units) is known and the bin number is wanted
peak_pos_count = htemp->GetBinContent(peak_pos_bin); //counts in particular bin
//------------------------------------------------------------------------------------------------------------------
zero_value_bin = htemp->GetXaxis()->FindBin(0.0); //if the bin value (in x-axis units) is known and the bin number is wanted
Nbins = htemp->GetSize() - 2; //total number of bins in histo
zero_bin_value = htemp->GetXaxis()->GetBinCenter(0); //if the bin number is known and the bin value (in x-axis units) is wanted.
max_bin_value = htemp->GetXaxis()->GetBinCenter(Nbins); //if the bin number is known and the bin value (in x-axis units) is wanted.
int TS = 0;
if (peak_pos >= zero_bin_value && peak_pos <= max_bin_value) //Make sure that TSpectrum peak is within histo range
{ //if not, use Par initial values from Basic Histo Peak Find
TS=1; //for cout below
low_limit = peak_pos - (0.1 * abs(max_bin_value-zero_bin_value)); //peakpos-10% of full range of histo
high_limit = peak_pos + (0.1 * abs(max_bin_value-zero_bin_value)); //peakpos+10% of full range of histo
func = new TF1("func", "gaus");
//func->FixParameter(1,0);
//.........这里部分代码省略.........
开发者ID:anmehta,项目名称:FNAL-Beam-Test-Scripts,代码行数:101,代码来源:doubleGausFit_withHistParameter.C
示例18: main
//.........这里部分代码省略.........
//Rebinning the histograms based on the mean value...
for (int p=0; p<nPMT; p++) {
for (int i=0; i<nBinsXY; i++) {
for (int j=0; j<nBinsXY; j++) {
hisxy[p][i][j]->GetXaxis()->SetRange(1,nBinHist);
Double_t mean = hisxy[p][i][j]->GetMean();
hisxy[p][i][j]->GetXaxis()->SetRange(0,nBinHist);
double nGroup = 4.*mean/200.;
nGroup = nGroup>1. ? nGroup : 1.;
hisxy[p][i][j]->Rebin((int)nGroup);
}
}
}
// Extracting mean from 200 keV peak
// Define fit ranges
double xLowBin[nPMT][nBinsXY][nBinsXY];
double xHighBin[nPMT][nBinsXY][nBinsXY];
double xLow[nPMT][nBinsXY][nBinsXY];
double xHigh[nPMT][nBinsXY][nBinsXY];
int maxBin[nPMT][nBinsXY][nBinsXY];
double maxCounts[nPMT][nBinsXY][nBinsXY];
double binCenterMax[nPMT][nBinsXY][nBinsXY];
double meanVal[nPMT][nBinsXY][nBinsXY]; // Holds the mean of the distribution as defined by first fitting the peak
double fitMean[nPMT][nBinsXY][nBinsXY]; // Holds the mean of the low energy peak
double fitSigma[nPMT][nBinsXY][nBinsXY]; // Holds the sigma of the low energy peak
double endpoint[nPMT][nBinsXY][nBinsXY]; // Holds the endpoint
/////// First determine roughly where the low energy peak is
TSpectrum *spec;
for (int p=0; p<nPMT; p++) {
for (int i=0; i<nBinsXY; i++) {
for (int j=0; j<nBinsXY; j++) {
double r = sqrt(power(posmap.getBinCenter(j),2)+power(posmap.getBinCenter(i),2));
// Find bin with maximum content
hisxy[p][i][j]->GetXaxis()->SetRange(2,nBinHist);
maxBin[p][i][j] = hisxy[p][i][j]->GetMaximumBin();
maxCounts[p][i][j] = hisxy[p][i][j]->GetBinContent(maxBin[p][i][j]);
binCenterMax[p][i][j] = hisxy[p][i][j]->GetBinCenter(maxBin[p][i][j]);
if (r<=(50.+2*xyBinWidth))
{
spec = new TSpectrum(20);
Int_t npeaks = spec->Search(hisxy[p][i][j],1.5,"",0.5);
if (npeaks==0)
{
cout << "No peaks identified at PMT" << p << " position " << posmap.getBinCenter(i) << ", " << posmap.getBinCenter(j) << endl;
}
else
{
Float_t *xpeaks = spec->GetPositionX(); // Note that newer versions of ROOT return a pointer to double...
TAxis *xaxis = (TAxis*)(hisxy[p][i][j]->GetXaxis());
Int_t peakBin=0;
Double_t BinSum=0.;
Double_t BinSumHold = 0.;
Int_t maxPeak=0.;
for (int pk=0;pk<npeaks;pk++) {
开发者ID:mabrow05,项目名称:ParallelAnalyzer-1,代码行数:67,代码来源:position_map.cpp
示例19: fit
void fit(){
int range = 1e4;
TFile *f = new TFile("hgrroot.root");
TCanvas *ca = new TCanvas("ca","ca",0,0,1400,900);
ca->Divide(4,4);
TH1F* h[7][4];
TF1* fu[7][4][2];
TF1* fus[7][4][2][3];
double res[7*4];
double det[7*4];
double slope[7*4];
double offset[7*4];
double resolution =5;
TEnv* output = new TEnv("corecal.dat");
for(int d=0;d<7;d++){
for(int c=0;c<4;c++){
h[d][c] = (TH1F*)f->Get(Form("htraceen_b%d_c%d_cr%d_d%d",0,9,c,d+1));
ca->cd(1+c*4);
//h[d][c]->GetXaxis()->SetRangeUser(100,4000);
TSpectrum *sp = new TSpectrum(2,resolution);
sp->SetResolution(resolution);
Int_t nfound = 0;
nfound = sp->Search(h[d][c],resolution,"nobackground",0.5);
Float_t *xpeaks = sp->GetPositionX();
Float_t *ypeaks = sp->GetPositionY();
// for(int p=0;p<nfound;p++){
// cout << xpeaks[p] << "\t" << ypeaks[p] << endl;
// }
if(nfound!=2){
cout << "Found " << nfound << " peaks in spectrum, too many, aborting" << endl;
continue;
}
h[d][c]->DrawCopy();
//check if first peak is lower in energy, otherwise swap them
if(xpeaks[0]>xpeaks[1]){
Float_t temp = xpeaks[1];
xpeaks[1] = xpeaks[0];
xpeaks[0] = temp;
temp = ypeaks[1];
ypeaks[1] = ypeaks[0];
ypeaks[0] = temp;
}
for(int p=0;p<nfound;p++){
ca->cd(1+c*4+1+p);
h[d][c]->GetXaxis()->SetRangeUser(xpeaks[p]-range,xpeaks[p]+range);
h[d][c]->DrawCopy();
fu[d][c][p] = new TF1(Form("fcore_d%d_c%d_p%d",d,c,p),fgammagaussbg,xpeaks[p]-range,xpeaks[p]+range,6);
fu[d][c][p]->SetLineColor(3);
fu[d][c][p]->SetLineWidth(1);
fu[d][c][p]->SetParameter(0,0);//bg const
fu[d][c][p]->SetParameter(1,0);//bg slope
fu[d][c][p]->SetParameter(2,h[d][c]->Integral(xpeaks[p]-range,xpeaks[p]+range));//norm
fu[d][c][p]->SetParameter(3,xpeaks[p]);//mean
fu[d][c][p]->SetParLimits(3,xpeaks[p]-500,xpeaks[p]+500);//mean
fu[d][c][p]->SetParameter(4,100);//sigma
fu[d][c][p]->SetParLimits(4,0.001,1000);//sigma
fu[d][c][p]->SetParameter(5,h[d][c]->GetBinContent(h[d][c]->FindBin(xpeaks[p]-range)));//step
h[d][c]->Fit(fu[d][c][p],"Rqn");
fu[d][c][p]->Draw("same");
fus[d][c][p][0] = new TF1(Form("fcore_d%d_c%d_p%d_bg",d,c,p),fgammabg,xpeaks[p]-range,xpeaks[p]+range,6);
fus[d][c][p][1] = new TF1(Form("fcore_d%d_c%d_p%d_st",d,c,p),fgammastep,xpeaks[p]-range,xpeaks[p]+range,
|
请发表评论