00001 #include "TH1.h"
00002 #include "TF1.h"
00003 #include "TAxis.h"
00004 #include "TSQLiteServer.h"
00005
00006 void FitParametersFirstGuess(TH1* h, TF1* f, const Double_t sigma) {
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 }
00026
00027
00028
00029 ge_energy_runbyrun_calib(const char* mergedb_name) {
00030 const unsigned int nsets = 7;
00031 const unsigned int npeaks = 8;
00032 const char* sets[nsets] = { "Al100", "Al50awithNDet2", "Al50awithoutNDet2", "Al50b", "Si16P", "SiR21pct", "SiR23pct"};
00033 const Double_t en[npeaks] = { 351.932, 510.998928, 583.191, 609.312, 911.204, 1173.237, 1332.501, 1460.83 };
00034 const Double_t en_err[npeaks] = { 0.002, 0.000011, 0.002, 0.007, 0.004, 0.004, 0.005, 0.01 };
00035 const Double_t sigma = 8.;
00036 const Double_t nbins = 100.;
00037 const Double_t bounds[npeaks][2] = { { 2800., 2800. + nbins },
00038 { 4100., 4100. + nbins },
00039 { 4700., 4700. + nbins },
00040 { 4900., 4900. + nbins },
00041 { 7360., 7360. + nbins },
00042 { 9500., 9500. + nbins },
00043 { 10800., 10800.+ nbins },
00044 { 11850., 11850. + nbins } };
00045 char dbname[256];
00046 sprintf(dbname, "sqlite://%s", mergedb_name);
00047 TSQLiteServer* mergedb = new TSQLiteServer(dbname);
00048 unsigned int nfiles[nsets] = { 0, 0, 0, 0, 0, 0, 0 };
00049 unsigned int nfiles_tmp = 0;
00050 for (unsigned int i = 0; i < nsets; ++i) {
00051 char cmd[128];
00052 sprintf(cmd, "SELECT COUNT(file) FROM Merge WHERE file LIKE '%s%%'", sets[i]);
00053 TSQLiteResult* res = (TSQLiteResult*)mergedb->Query(cmd);
00054 TSQLiteRow* row = (TSQLiteRow*)res->Next();
00055 nfiles[i] = (unsigned int)atoi(row->GetField(0));
00056 if (nfiles[i] > nfiles_tmp) nfiles_tmp = nfiles[i];
00057 delete row;
00058 delete res;
00059 }
00060 const unsigned int nfiles_max = nfiles_tmp;
00061 Double_t time[nsets][nfiles_max];
00062 Double_t peaks[nsets][nfiles_max][npeaks], sigmas[nsets][nfiles_max][npeaks], peaks_err[nsets][nfiles_max][npeaks], sigmas_err[nsets][nfiles_max][npeaks];
00063 Double_t chi2[nsets][nfiles_max][npeaks], ndf[nsets][nfiles_max][npeaks];
00064 Int_t status[nsets][nfiles_max][npeaks];
00065 for (unsigned int i = 0; i < nsets; ++i) {
00066 char cmd[128];
00067 sprintf(cmd, "SELECT time FROM Merge WHERE file LIKE '%s%%' ORDER BY time ASC", sets[i]);
00068 TSQLiteResult* res = (TSQLiteResult*)mergedb->Query(cmd);
00069 for (unsigned int j = 0; j < nfiles[i]; ++j) {
00070 TSQLiteRow* row = (TSQLiteRow*)res->Next();
00071 time[i][j] = atof(row->GetField(0));
00072 delete row;
00073 }
00074 delete res;
00075 }
00076
00077 char fname[64];
00078 for (unsigned int i = 0; i < nsets; ++i) {
00079 for (unsigned int j = 0; j < nfiles[i]; ++j) {
00080 sprintf(fname, "%s_%u.root", sets[i], j+1);
00081 TFile* f = new TFile(fname, "READ");
00082 TH1* h = (TH1*)f->Get("GeSpectrum/hEnergyFarOOT");
00083
00084 h->GetXaxis()->SetRangeUser(3500., 4500.);
00085 Double_t shift_511 = (bounds[1][0] + bounds[1][1])/2. - (Double_t)h->GetMaximumBin();
00086 for (unsigned int k = 0; k < npeaks; ++k) {
00087 Double_t shift = shift_511*(1.+((bounds[k][0]+bounds[k][1])/2.-4100.)/6000.);
00088 TF1 fit("fitfunc", "gaus(0)+pol1(3)");
00089 h->GetXaxis()->SetRangeUser(bounds[k][0] - shift, bounds[k][1] - shift);
00090 FitParametersFirstGuess(h, &fit, sigma);
00091 TFitResultPtr fitres = h->Fit(&fit, "SEQM");
00092 if (!fitres->IsValid()) {
00093 printf("Peak: %s(%u)\tStatus: %d\n", fname, k, (Int_t)fitres);
00094 return;
00095 }
00096 peaks[i][j][k] = fitres->Parameter(1); peaks_err[i][j][k] = fitres->ParError(1);
00097 sigmas[i][j][k] = fitres->Parameter(2); sigmas_err[i][j][k] = fitres->ParError(2);
00098 status[i][j][k] = (Int_t)fitres; chi2[i][j][k] = fitres->Chi2(); ndf[i][j][k] = (Double_t) fitres->Ndf();
00099 }
00100 }
00101 }
00102 TMultiGraph* mg[nsets];
00103 char grname[64], grtitle[64];
00104 for (unsigned int i = 0; i < nsets; ++i) {
00105 sprintf(grname, "mg%s", sets[i]);
00106 sprintf(grtitle, "Gain (%s)", sets[i]);
00107 mg[i] = new TMultiGraph(grname, grtitle);
00108 }
00109 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];
00110 for (unsigned int i = 0; i < nsets; ++i) {
00111 printf("--------------------------------------------------------------------------------\n");
00112 printf("%s\nFile\tGain\tPedestal\tChi2\tNDF\Chi2/NDF\tGain Pct Err\n", sets[i]);
00113 for (unsigned int j = 0; j < nfiles[i]; ++j) {
00114 TGraphErrors* gr = new TGraphErrors(npeaks, peaks[i][j], en, peaks_err[i][j], en_err);
00115 TFitResultPtr fitres = gr->Fit("pol1", "SQ");
00116 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));
00117 mg[i]->Add(gr); gr->GetFunction("pol1")->SetLineColor(j+1); gr->SetMarkerColor(j+1); gr->SetMarkerStyle(1);
00118 gain[i][j] = fitres->Parameter(1); gain_err[i][j] = fitres->ParError(1);
00119 offset[i][j] = fitres->Parameter(0); offset_err[i][j] = fitres->ParError(0);
00120 chi2_cal[i][j] = fitres->Chi2(); ndf_cal[i][j] = (Double_t)fitres->Ndf();
00121 }
00122 }
00123 TMultiGraph* mg_gain = new TMultiGraph("mg_gain", "Gain Drift");
00124 TLegend* leg_gain = new TLegend(0.7,0.7,0.9,0.9);
00125 for (unsigned int i = 0; i < nsets; ++i) {
00126 new TCanvas();
00127 TGraphErrors* gr = new TGraphErrors(nfiles[i], time[i], gain[i], 0x0, gain_err[i]);
00128 mg_gain->Add(gr); gr->SetMarkerColor(i+1); gr->SetMarkerStyle(1);
00129 mg[i]->Draw("A*");
00130 leg_gain->AddEntry(gr, sets[i]);
00131 }
00132 new TCanvas();
00133 mg_gain->Draw("A*");
00134 mg_gain->GetXaxis()->SetTitle("Unix Time (s)");
00135 mg_gain->GetYaxis()->SetTitle("Gain (E=Gain*ADC+Offset)");
00136 leg_gain->Draw("SAME");
00137
00138 for (unsigned int i = 0; i < nsets; ++i) {
00139 for (unsigned int j = 0; j < nfiles[i]; ++j) {
00140 char cmd[128];
00141 sprintf(cmd, "SELECT runs FROM Merge WHERE file=='%s_%u'", sets[i], j+1);
00142 TSQLiteResult* res = (TSQLiteResult*)mergedb->Query(cmd);
00143 TSQLiteRow* row = (TSQLiteRow*)res->Next();
00144 std::stringstream ss(row->GetField(0));
00145 delete res; delete row;
00146 while (ss.good()) {
00147 char ofname[128];
00148 unsigned int run;
00149 ss >> run;
00150 sprintf(ofname, "calib.run%05u.Energy.csv", run);
00151 std::ofstream ofile(ofname);
00152 ofile << "run,channel,gain,gain_err,offset,offset_err,chi2,ndf" << std::endl;
00153 ofile << run << ",Ge-S," << gain[i][j] << ',' << gain_err[i][j] << ','
00154 << offset[i][j] << ',' << offset_err[i][j] << ','
00155 << chi2_cal[i][j] << ',' << ndf_cal[i][j] << std::endl;
00156 ofile.close();
00157 }
00158 }
00159 }
00160 }
