/// \file /// \ingroup tutorial_tree /// \notebook -nodraw /// Read data from an ascii file and create a root file with an histogram and an ntuple. /// See a variant of this macro in basic2.C. /// /// \macro_code /// /// \author Rene Brun #include "TH1.h" #include "TF1.h" #include "Riostream.h" void basic() { // read file $ROOTSYS/tutorials/tree/basic.dat // this file has 3 columns of float data ifstream in; in.open("eu_gu_1.txt"); auto f = TFile::Open("basic.root","RECREATE"); Int_t i, npeaks = 5; TCanvas *c1 = new TCanvas("c1","c1"); TGraph *g = new TGraph("eu_gu_1.txt"); // g->Sort(); // just a precaution Int_t n = g->GetN();//get the no. of points on the TGraph *g Double_t source[n]; Double_t *x = new Double_t[(n + 1)]; x[0] = (3.0 * g->GetX()[0] - g->GetX()[1]) / 2.0; x[n] = (3.0 * g->GetX()[(n - 1)] - g->GetX()[(n - 2)]) / 2.0; for (Int_t i = 1; i < n; i++) { x[i] = (g->GetX()[(i - 1)] + g->GetX()[i]) / 2.0; } TH1D *h = new TH1D("h", "spectrum;x-value(Channels);#counts", n, x); h->FillN(n, g->GetX(), g->GetY()); h->Sumw2(kFALSE); // restore proper errors delete [] x; // no longer needed delete g; // no longer needed TF1 *f1 = new TF1("f1", "gausn(0) + pol1(3)", 410., 422.); // linear background f1->SetParameters(2800., 415., 2., 0., 0.); h->Fit("f1", "LBR+"); /*std::cout << "Peak Area 121.74 = " << f1->Integral(410., 422.) << std::endl; std::cout << "background = " << f1->Integral(410., 422.) - f1->GetParameter(0) << std::endl;*/ //f1->Print("V"); TF1 *f2 = new TF1("f2", "gausn(0) + pol1(3)", 831., 843.); // linear background f2->SetParameters(1000., 838., 2., 0., 0.); h->Fit("f2", "LBR+"); std::cout << "Peak Area 121.74 = " << f1->Integral(410., 422.) << std::endl; std::cout << "background = " << f1->Integral(410., 422.) - f1->GetParameter(0) << std::endl; TF1 *f3 = new TF1("f3", "gausn(0) + pol1(3)", 1172., 1184.); f3->SetParameters(2700., 1178., 2., 0., 0.); h->Fit("f3", "LBR+"); std::cout << "Peak Area 121.74 = " << f1->Integral(410., 422.) << std::endl; std::cout << "background = " << f1->Integral(410., 422.) - f1->GetParameter(0) << std::endl; TF1 *f4 = new TF1("f4", "gausn(0) + pol1(3)", 1512., 1526.); // clean peak, no background f4->SetParameters(280., 1519., 3., 0., 0.); h->Fit("f4", "LBR+"); std::cout << "Peak Area 121.74 = " << f1->Integral(410., 422.) << std::endl; std::cout << "background = " << f1->Integral(410., 422.) - f1->GetParameter(0) << std::endl; TF1 *f5 = new TF1("f5", "gausn(0) + pol1(3)", 2656., 2673.); // clean peak, no background f5->SetParameters(580., 2664., 3., 0., 0.); h->Fit("f5", "LBR+"); TF1 *f6 = new TF1("f6", "gausn(0) + pol1(3)", 2958., 2976.); // clean peak, no background f6->SetParameters(200., 2966., 3., 0., 0.); h->Fit("f6", "LBR+"); TF1 *f7 = new TF1("f7", "gausn(0) + pol1(3)", 3288., 3307.); // clean peak, no background f7->SetParameters(540., 3296., 3., 0., 0.); h->Fit("f7", "LBR+"); TF1 *f8 = new TF1("f8", "gausn(0) + pol1(3)", 3794., 3814.); // clean peak, no background f8->SetParameters(420., 3805., 3., 0., 0.); h->Fit("f8", "LBR+"); TF1 *f9 = new TF1("f9", "gausn(0) + pol1(3)", 4804., 4826.); // clean peak, no background f9->SetParameters(500., 4815., 3., 0., 0.); h->Fit("f9", "LBR+"); // e.g. "BR+" or "LBR+"*/ h->SetAxisRange(0., 8191.,"X"); h->SetAxisRange(-100., 3000.,"Y"); //Area1 /* std::cout << h->Integral(h->FindFixBin(410.), h->FindFixBin(422.)) << std::endl; std::cout << h->Integral(h->FindFixBin(831.), h->FindFixBin(843.)) << std::endl; std::cout << h->Integral(h->FindFixBin(1172.), h->FindFixBin(1184.)) << std::endl; std::cout << h->Integral(h->FindFixBin(1512.), h->FindFixBin(1526.)) << std::endl; std::cout << h->Integral(h->FindFixBin(2656.), h->FindFixBin(2673.)) << std::endl; std::cout << h->Integral(h->FindFixBin(2958.), h->FindFixBin(2976.)) << std::endl; std::cout << h->Integral(h->FindFixBin(3288.), h->FindFixBin(3307.)) << std::endl; std::cout << h->Integral(h->FindFixBin(3794.), h->FindFixBin(3814.)) << std::endl; std::cout << h->Integral(h->FindFixBin(4804.), h->FindFixBin(4826.)) << std::endl; TSpectrum *s = new TSpectrum(2*npeaks); Int_t nfound = s->Search(h,2,"",0.06); printf("Found %d candidate peaks to fit\n",nfound); s->Background(h->GetArray() + 1, n, 10, TSpectrum::kBackDecreasingWindow, TSpectrum::kBackOrder8, kTRUE, TSpectrum::kBackSmoothing7, kTRUE); TH1 *hb = s->Background(h,20,"same"); hb->SetLineColor(kBlack); gPad->SetLogy(1); h->SetMaximum(300.);//along Y h->Draw(); h->SetAxisRange(0., 150.,"Y"); // create a TF1 with the range from 10 to 1900 and 4 parameters TF1 *RT = new TF1("RT",Fit1,0,2000,3); RT->FixParameter(0,1900.0); //<==== RT->FixParameter(1,0.0028); RT->FixParameter(2,0.); h1->Draw(); h1->Fit("RT","R"); // Writing fuction fitted data in an histogram TH1F *h2 = new TH1F("h2","background",1900,-0.5,1899.5); h2->SetDirectory(0); for (int i=0;i<=1899;i++) //<====== { float x=h2->GetBinCenter(i); float y=RT->Eval(x); h2->SetStats(0); h2->Fill(x,y); //<====== h2->GetXaxis()->CenterTitle(); h2->GetYaxis()->CenterTitle(); h2->GetXaxis()->SetTitle("Energy (keV)"); h2->GetYaxis()->SetTitle("Counts"); h2->SetLineWidth(2); h2->SetLineColor(2); } TString resultfile="background.root"; TFile* fresult=TFile::Open( resultfile.Data(), "RECREATE" ); //fresult->cd(); h2->Write(); h2->SetStats(0); h2->GetXaxis()->CenterTitle(); h2->GetYaxis()->CenterTitle(); h2->GetXaxis()->SetTitle("Energy (keV)"); h2->GetYaxis()->SetTitle("Counts"); //MyC->cd(3); h2->Draw(); cout <<"Wrote result to "<Clone(); diff->SetDirectory(0); diff->SetName("background_subtracted "); diff->Add(h2,-1.0); diff->SetStats(0); diff->GetXaxis()->CenterTitle(); diff->GetYaxis()->CenterTitle(); diff->GetXaxis()->SetTitle("Energy (keV)"); diff->GetYaxis()->SetTitle("Counts"); diff->GetXaxis()->SetRangeUser(180,2000); diff->GetYaxis()->SetRangeUser(0,3000); TString resultfile1="background_subtracted.root"; TFile* fresult1=TFile::Open( resultfile1.Data(), "RECREATE" ); diff->Write(); diff->Draw();*/ // h->Print("all"); in.close(); f->Write(); }