I want to give some data to someone who’s not a root user but an Origin one.
What I am doing is subtracting a background from a histogram, by fitting it using a user defined function.
What I want is 3 dat files
- One for the total histogram
- One for the fitting function(that will actually be like a histogram in the sense of channel-value)
- One for the subtracted histogram
I’ve found a code that exports histogram data in an ascii file, but I don’t know how to apply this in a function, so as to have channel for x and value for y. Any idea will be more than welcome!
My code is the following
[code] #include “Riostream.h”
Bool_t reject;
Double_t exclude(Double_t *x, Double_t *par)
{
if (reject && x[0] > 610 && x[0] < 650) {//<-------------------------------------------
TF1::RejectPoint();
return 0;
}
//return par[0] + par[1]*x[0];
return par[0] + par[1]*x[0] + par[2]*x[0]*x[0];
//return par[0] + par[1]*x[0] + par[2]*x[0]*x[0] + par[3]*x[0]*x[0]*x[0];
//return par[0] + par[1]*x[0] + par[2]*x[0]*x[0] + par[3]*x[0]*x[0]*x[0] + par[4]*x[0]*x[0]*x[0]*x[0];
}
void testFitExclude(){//(char * file_c) {
gROOT->SetStyle("Plain");
gStyle->SetOptStat(0000);
gStyle->SetOptFit(1111);
gStyle->SetOptTitle(0);
TCanvas *c = new TCanvas("c", "c");
c->SetFillColor(5);
c->SetFrameFillColor(10);
TString dir = gSystem->UnixPathName(gInterpreter->GetCurrentMacroName());
dir.ReplaceAll("testFitExclude.C","");
dir.ReplaceAll("/./","/");
ifstream in;
in.open(Form("%sB_p_3300_120.dat",dir.Data()));//<-------------------------------------------
Float_t x,y;//,z,w,k,l,m;
Int_t nlines = 0;
TH1F *h1 = new TH1F("h1","Experimental",1024,1,1024);
while (1) {
in >> x >> y;
if (!in.good()) break;
if (nlines < 5) printf("x=%8f, y=%8f\n");
h1->Fill(x,y);
nlines++;
}
printf("Found %d points\n",nlines);
in.close();
h1->SetLineColor(kBlack);
h1->SetLineWidth(1);
h1->Draw();
h1->GetXaxis()->SetTitle("Channels");
h1->GetYaxis()->SetTitle("Events");
reject = kTRUE;
TF1 *fit = new TF1("fit",exclude,0,1024,3);
fit->SetParameters(0,0,0);
TH1F *hnew = (TH1F*)h1->Clone("hnew");
h1->Fit("fit","W",NULL,550,850);//<----------------------------------------
reject = kFALSE;
//TF1 *fit = new TF1("fit","pol2",0,1024);
// h1->Fit("fit","W",NULL,465,850);//465 850
//fit->Draw("same");
//fit->SetLineColor(kRed);
h1->GetFunction("fit")->SetLineWidth(4);
h1->GetFunction("fit")->SetLineColor(kRed);
TF1 *extraPol = new TF1("extraPol",exclude,0,1024,5);
//TF1 *extraPol = new TF1("extraPol","pol2",0,1024);
extraPol->SetParameters(fit->GetParameters());
extraPol->Draw("same");
h1->GetFunction("fit")->Draw("same");
//extraPol->SetLineStyle(7);
extraPol->SetLineWidth(2);
//TH1F *hnew = (TH1F*)h1->Clone("hnew");
hnew->Add(extraPol, -1);
hnew->Draw("same");
hnew->SetLineColor(kMagenta);
TString histfilename150 = "B_p_3300_120_Subtracted.dat";//<------------------------------TString::Format("%s_Projected_150.dat",file.Data());
SingleExportAscii(hnew,histfilename150);
//Legend *leg = c->BuildLegend();
//print(fline->GetParameters());
}
/**
- \brief Export Single Histogram into ASCII file
/
Bool_t SingleExportAscii(TH1 hist, TString &filename, TString folder="", TString separator="\t")
{
Int_t i,j;
Double_t xcenter, xwidth;
Bool_t success=kFALSE;
//filename = folder + hist->GetName() + “.dat”;
ofstream file_out(filename);
//file_out << "# Output " << hist->ClassName() << “: " << hist->GetName() << " (” << hist->GetTitle() << “)\n”;
if (hist->GetDimension()==1)
{
//file_out << “# BinCenter” << separator << “Content\n”;
for (i=1; i<=hist->GetNbinsX(); i++)
file_out << int(hist->GetBinCenter(i)) << separator << hist->GetBinContent(i) << endl;
if (i>1)
success=kTRUE;
}
else if (hist->GetDimension()==2)
{
file_out << “# xBinCenter” << separator << “yBinCenter” << separator << “Content” << separator << “xBinHalfWidth” << separator << “yBinHalfWidth” << separator << “Error” << endl;
for (i=1; i <= hist->GetNbinsX(); i++)
{
xcenter = hist->GetXaxis()->GetBinCenter(i);
xwidth = hist->GetXaxis()->GetBinWidth(i)/2;
for (j=1; j <= hist->GetNbinsY(); j++)
file_out << xcenter << separator << hist->GetYaxis()->GetBinCenter(j) << separator << hist->GetBinContent(i,j) << separator << xwidth << separator << hist->GetYaxis()->GetBinWidth(j)/2 << separator << hist->GetBinError(i,j) << endl;
if (j>1)
file_out << endl; // produce a blank line after each set of Ybins for a certain Xbin. Gnuplot likes this.
}
if (i>1)
success=kTRUE;
}
file_out.close();
if (success == kTRUE)
cout << "*** TRemiHistExport: Histogram " << hist->GetName() << " written to " << filename << endl;
return success;
}[/code]