rootana/scripts/ge/ge_energy_runbyrun_calib.c File Reference

#include "TH1.h"
#include "TF1.h"
#include "TAxis.h"
#include "TSQLiteServer.h"

Functions
void	FitParametersFirstGuess (TH1 h, TF1 f, const Double_t sigma)
	ge_energy_runbyrun_calib (const char *mergedb_name)

Function Documentation

void FitParametersFirstGuess	(	TH1 *	h,
		TF1 *	f,
		const Double_t	sigma
	)

Definition at line 6 of file ge_energy_runbyrun_calib.c.

Referenced by ge_energy_runbyrun_calib().

00006                                                                    {
00007   const unsigned int npar = 5;
00008   Double_t par[npar];
00009 
00010   TAxis* x = h->GetXaxis();
00011   const Double_t n = x->GetLast() - x->GetFirst() + 1.;
00012 
00013   par[4] = 0.;
00014   par[2] = sigma;
00015   par[1] = (Double_t)h->GetMaximumBin();
00016   h->GetXaxis()->SetRangeUser(par[1] - n/2., par[1] + n/2.);
00017   par[3] = h->Integral()/n;
00018   par[0] = h->GetMaximum() - par[3];
00019 
00020   f->SetParameters(par);
00021   f->SetParLimits(0, 0., h->GetMaximum());
00022   f->SetParLimits(1, x->GetFirst(), x->GetLast());
00023   f->SetParLimits(2, 1., 3.*sigma);
00024   return;
00025 }

ge_energy_runbyrun_calib ( const char * mergedb_name )

Definition at line 29 of file ge_energy_runbyrun_calib.c.

References FitParametersFirstGuess().

00030 00031 00032 00033 00034 00035 00036 00037 00038 00039 00040 00041 00042 00043 00044 00045 00046 00047 00048 00049 00050 00051 00052 00053 00054 00055 00056 00057 00058 00059   } 00060 00061 00062 00063 00064 00065 00066 00067 00068 00069 00070 00071 00072 00073     } 00074 00075   } 00076 00077 00078 00079 00080 00081 00082 00083 00084 00085 00086 00087 00088 00089 00090 00091 00092 00093 00094 00095         } 00096 00097 00098 00099       } 00100     } 00101   } 00102 00103 00104 00105 00106 00107 00108   } 00109 00110 00111 00112 00113 00114 00115 00116 00117 00118 00119 00120 00121     } 00122   } 00123 00124 00125 00126 00127 00128 00129 00130 00131   } 00132 00133 00134 00135 00136 00137 00138 00139 00140 00141 00142 00143 00144 00145 00146 00147 00148 00149 00150 00151 00152 00153 00154 00155 00156 00157       } 00158     } 00159   } 00160 }

class="fragment">00029 { const unsigned int nsets = 7; const unsigned int npeaks = 8; const char* sets[nsets] = { "Al100", "Al50awithNDet2", "Al50awithoutNDet2", "Al50b", "Si16P", "SiR21pct", "SiR23pct"}; const Double_t en[npeaks] = { 351.932, 510.998928, 583.191, 609.312, 911.204, 1173.237, 1332.501, 1460.83 }; const Double_t en_err[npeaks] = { 0.002, 0.000011, 0.002, 0.007, 0.004, 0.004, 0.005, 0.01 }; const Double_t sigma = 8.; const Double_t nbins = 100.; const Double_t bounds[npeaks][2] = { { 2800., 2800. + nbins }, { 4100., 4100. + nbins }, { 4700., 4700. + nbins }, { 4900., 4900. + nbins }, { 7360., 7360. + nbins }, { 9500., 9500. + nbins }, { 10800., 10800.+ nbins }, { 11850., 11850. + nbins } }; char dbname[256]; sprintf(dbname, "sqlite://%s", mergedb_name); TSQLiteServer* mergedb = new TSQLiteServer(dbname); unsigned int nfiles[nsets] = { 0, 0, 0, 0, 0, 0, 0 }; unsigned int nfiles_tmp = 0; for (unsigned int i = 0; i < nsets; ++i) { char cmd[128]; sprintf(cmd, "SELECT COUNT(file) FROM Merge WHERE file LIKE '%s%%'", sets[i]); TSQLiteResult* res = (TSQLiteResult*)mergedb->Query(cmd); TSQLiteRow* row = (TSQLiteRow*)res->Next(); nfiles[i] = (unsigned int)atoi(row->GetField(0)); if (nfiles[i] > nfiles_tmp) nfiles_tmp = nfiles[i]; delete row; delete res; const unsigned int nfiles_max = nfiles_tmp; Double_t time[nsets][nfiles_max]; Double_t peaks[nsets][nfiles_max][npeaks], sigmas[nsets][nfiles_max][npeaks], peaks_err[nsets][nfiles_max][npeaks], sigmas_err[nsets][nfiles_max][npeaks]; Double_t chi2[nsets][nfiles_max][npeaks], ndf[nsets][nfiles_max][npeaks]; Int_t status[nsets][nfiles_max][npeaks]; for (unsigned int i = 0; i < nsets; ++i) { char cmd[128]; sprintf(cmd, "SELECT time FROM Merge WHERE file LIKE '%s%%' ORDER BY time ASC", sets[i]); TSQLiteResult* res = (TSQLiteResult*)mergedb->Query(cmd); for (unsigned int j = 0; j < nfiles[i]; ++j) { TSQLiteRow* row = (TSQLiteRow*)res->Next(); time[i][j] = atof(row->GetField(0)); delete row; delete res; char fname[64]; for (unsigned int i = 0; i < nsets; ++i) { for (unsigned int j = 0; j < nfiles[i]; ++j) { sprintf(fname, "%s_%u.root", sets[i], j+1); TFile* f = new TFile(fname, "READ"); TH1* h = (TH1*)f->Get("GeSpectrum/hEnergyFarOOT"); // Get first guess at shift h->GetXaxis()->SetRangeUser(3500., 4500.); Double_t shift_511 = (bounds[1][0] + bounds[1][1])/2. - (Double_t)h->GetMaximumBin(); for (unsigned int k = 0; k < npeaks; ++k) { Double_t shift = shift_511*(1.+((bounds[k][0]+bounds[k][1])/2.-4100.)/6000.); TF1 fit("fitfunc", "gaus(0)+pol1(3)"); h->GetXaxis()->SetRangeUser(bounds[k][0] - shift, bounds[k][1] - shift); FitParametersFirstGuess(h, &fit, sigma); TFitResultPtr fitres = h->Fit(&fit, "SEQM"); if (!fitres->IsValid()) { printf("Peak: %s(%u)\tStatus: %d\n", fname, k, (Int_t)fitres); return; peaks[i][j][k] = fitres->Parameter(1); peaks_err[i][j][k] = fitres->ParError(1); sigmas[i][j][k] = fitres->Parameter(2); sigmas_err[i][j][k] = fitres->ParError(2); status[i][j][k] = (Int_t)fitres; chi2[i][j][k] = fitres->Chi2(); ndf[i][j][k] = (Double_t) fitres->Ndf(); TMultiGraph* mg[nsets]; char grname[64], grtitle[64]; for (unsigned int i = 0; i < nsets; ++i) { sprintf(grname, "mg%s", sets[i]); sprintf(grtitle, "Gain (%s)", sets[i]); mg[i] = new TMultiGraph(grname, grtitle); Double_t gain[nsets][nfiles_max], offset[nsets][nfiles_max], gain_err[nsets][nfiles_max], offset_err[nsets][nfiles_max], chi2_cal[nsets][nfiles_max], ndf_cal[nsets][nfiles_max]; for (unsigned int i = 0; i < nsets; ++i) { printf("--------------------------------------------------------------------------------\n"); printf("%s\nFile\tGain\tPedestal\tChi2\tNDF\Chi2/NDF\tGain Pct Err\n", sets[i]); for (unsigned int j = 0; j < nfiles[i]; ++j) { TGraphErrors* gr = new TGraphErrors(npeaks, peaks[i][j], en, peaks_err[i][j], en_err); TFitResultPtr fitres = gr->Fit("pol1", "SQ"); printf("%u\t%g\t%g\t%g\t%g\t%g\t%g\n", j, fitres->Value(1), fitres->Value(0), fitres->Chi2(), (Double_t)fitres->Ndf(), fitres->Chi2()/(Double_t)fitres->Ndf(), fitres->ParError(1)/fitres->Value(1)); mg[i]->Add(gr); gr->GetFunction("pol1")->SetLineColor(j+1); gr->SetMarkerColor(j+1); gr->SetMarkerStyle(1); gain[i][j] = fitres->Parameter(1); gain_err[i][j] = fitres->ParError(1); offset[i][j] = fitres->Parameter(0); offset_err[i][j] = fitres->ParError(0); chi2_cal[i][j] = fitres->Chi2(); ndf_cal[i][j] = (Double_t)fitres->Ndf(); TMultiGraph* mg_gain = new TMultiGraph("mg_gain", "Gain Drift"); TLegend* leg_gain = new TLegend(0.7,0.7,0.9,0.9); for (unsigned int i = 0; i < nsets; ++i) { new TCanvas(); TGraphErrors* gr = new TGraphErrors(nfiles[i], time[i], gain[i], 0x0, gain_err[i]); mg_gain->Add(gr); gr->SetMarkerColor(i+1); gr->SetMarkerStyle(1); mg[i]->Draw("A*"); leg_gain->AddEntry(gr, sets[i]); new TCanvas(); mg_gain->Draw("A*"); mg_gain->GetXaxis()->SetTitle("Unix Time (s)"); mg_gain->GetYaxis()->SetTitle("Gain (E=Gain*ADC+Offset)"); leg_gain->Draw("SAME"); for (unsigned int i = 0; i < nsets; ++i) { for (unsigned int j = 0; j < nfiles[i]; ++j) { char cmd[128]; sprintf(cmd, "SELECT runs FROM Merge WHERE file=='%s_%u'", sets[i], j+1); TSQLiteResult* res = (TSQLiteResult*)mergedb->Query(cmd); TSQLiteRow* row = (TSQLiteRow*)res->Next(); std::stringstream ss(row->GetField(0)); delete res; delete row; while (ss.good()) { char ofname[128]; unsigned int run; ss >> run; sprintf(ofname, "calib.run%05u.Energy.csv", run); std::ofstream ofile(ofname); ofile << "run,channel,gain,gain_err,offset,offset_err,chi2,ndf" << std::endl; ofile << run << ",Ge-S," << gain[i][j] << ',' << gain_err[i][j] << ',' << offset[i][j] << ',' << offset_err[i][j] << ',' << chi2_cal[i][j] << ',' << ndf_cal[i][j] << std::endl; ofile.close();

rootana/scripts/ge/ge_energy_runbyrun_calib.c File Reference

Functions

Function Documentation