00001 {
00002 const char* files[] = { "add1.root", "add2.root", "add3.root", "add4.root", "add5.root",
00003 "add6.root", "add7.root", "add8.root", "add9.root", "add10.root",
00004 "add11.root", "add12.root", "add13.root", "add14.root", "add15.root",
00005 "add16.root", "add17.root", "add18.root", "add19.root", "add20.root",
00006 "add21.root", "add22.root", "add23.root", "add24.root", "add25.root",
00007 "add26.root", "add27.root", "add28.root", "add29.root", "add30.root" };
00008 const unsigned int nfiles = 30;
00009 const unsigned int npeaks = 6;
00010 const Double_t peak_search[npeaks] = { 1453., 1325., 967., 933., 606., 508.};
00011 const Double_t search_window = 20.;
00012 const Double_t sigma = 1.;
00013 Double_t peaks[nfiles][npeaks];
00014 Double_t sig[nfiles][npeaks];
00015 Double_t peaks_err[nfiles][npeaks];
00016 Double_t sig_err[nfiles][npeaks];
00017 unsigned int t[nfiles] = { 1387604374, 1387606953, 1387609839, 1387612731, 1387615575,
00018 1387618372, 1387621403, 1387624574, 1387627327, 1387630073,
00019 1387633340, 1387636613, 1387639373, 1387642307, 1387645236,
00020 1387648020, 1387650798, 1387662449, 1387667155, 1387669933,
00021 1387672708, 1387675469, 1387678264, 1387681099, 1387683923,
00022 1387686724, 1387689536, 1387692364, 1387697564, 1387700883 };
00023 const unsigned int first_run_start = 1387603224;
00024 const double total_livetime = 104090.;
00025 for (unsigned int i = 0; i < nfiles; ++i)
00026 t[i] -= first_run_start;
00027
00028
00029 TF1* func = new TF1("pkfunc", "gaus(0)+[3]*x+[4]", 400., 1600.);
00030 for (unsigned int i = 0; i < nfiles; ++i) {
00031 printf("File %d\n", i+1);
00032 char path[256];
00033 TFile f(strcat(strcpy(path, "/home/jrquirk/tmp/ge/out/v37/"), files[i]), "READ");
00034 TH1* h = (TH1*)f.Get("GeSpectrum/hEnergy");
00035 for (unsigned int j = 0; j < npeaks; ++j) {
00036 h->GetXaxis()->SetRangeUser(peak_search[j] - search_window, peak_search[j] + search_window);
00037 func->SetParameters(h->GetBinContent(h->FindBin(peak_search[j])),
00038 peak_search[j],
00039 sigma,
00040 0.,
00041 (h->GetBinContent(h->FindBin(peak_search[j] - 20.)) + h->GetBinContent(h->FindBin(peak_search[j] + 20.)))/2.);
00042 TFitResultPtr res = h->Fit(func, "SNQ");
00043 peaks[i][j] = func->GetParameter(1);
00044 peaks_err[i][j] = func->GetParError(1);
00045 sig[i][j] = func->GetParameter(2);
00046 sig_err[i][j] = func->GetParError(2);
00047 }
00048 }
00049
00050
00051 FILE* ofile = fopen("peaks.csv", "w");
00052 fprintf(ofile, "Time Peak +/-Peak Sigma +/-Sigma");
00053 for (unsigned int i = 0; i < nfiles; ++i) {
00054 fprintf(ofile, "%u ", t[i]);
00055 for (unsigned int j = 0; j < npeaks; ++j)
00056 fprintf(ofile, " %5g %5.5g %5g %5.5g", peaks[i][j], peaks_err[i][j], sig[i][j], sig_err[i][j]);
00057 fprintf(ofile, "\n");
00058 }
00059 fclose(ofile);
00060
00061
00062 TLegend* leg = new TLegend(0.2, 0.53, 0.9, 0.755, "Gain Drifts");
00063 TLegend* leg_norm = new TLegend(0.6, 0.1, 0.8, 0.4, "Gain Drifts");
00064 TLegend* leg_shift = new TLegend(0.6, 0.1, 0.8, 0.4, "Gain Drifts");
00065 TLegend* leg_sig = new TLegend(0.5, 0.1, 0.95, 0.4, "Sigma Drifts");
00066 TLegend* leg_res = new TLegend(0.5, 0.1, 0.65, 0.4, "Resolution Drifts");
00067 TMultiGraph* mg = new TMultiGraph("peakpos", "Germanium Peaks");
00068 TMultiGraph* mg_norm = new TMultiGraph("peakpos_norm", "Germanium Peaks (Normalized)");
00069 TMultiGraph* mg_shift = new TMultiGraph("peakpos_shift", "Germanium Peaks (Shifted)");
00070 TMultiGraph* mg_sig = new TMultiGraph("peaksigma", "Germanium Peak Sigma");
00071 TMultiGraph* mg_res = new TMultiGraph("peakresolution", "Germanium Peak Resolution");
00072 for (unsigned int i = 0; i < npeaks; ++i) {
00073 char fmt_pk[256]; memset(fmt_pk, '\0', 256); strcat(fmt_pk, "%lg");
00074 char fmt_sig[256]; memset(fmt_sig, '\0', 256); strcat(fmt_sig, "%lg");
00075 char lab[256]; memset(lab, '\0', 256);
00076 for (unsigned int j = 0; j < npeaks; ++j) {
00077 if (j==i) {
00078 strcat(fmt_pk, " %lg %lg %*lg %*lg");
00079 strcat(fmt_sig, " %*lg %*lg %lg %lg");
00080 } else {
00081 strcat(fmt_pk, " %*lg %*lg %*lg %*lg");
00082 strcat(fmt_sig, " %*lg %*lg %*lg %*lg");
00083 }
00084 }
00085 TGraphErrors* gr = new TGraphErrors("peaks.csv", fmt_pk);
00086 TGraphErrors* gr_norm = (TGraphErrors*)gr->Clone();
00087 TGraphErrors* gr_shift = (TGraphErrors*)gr->Clone();
00088 TGraphErrors* gr_sig = new TGraphErrors("peaks.csv", fmt_sig);
00089 TGraphErrors* gr_res = (TGraphErrors*)gr_sig->Clone();
00090 Double_t norm = gr_norm->GetY()[0];
00091 for (unsigned int j = 0; j < gr_norm->GetN(); ++j) {
00092 gr_norm->GetY()[j] = gr_norm->GetY()[j] / norm;
00093 gr_norm->GetEY()[j] = gr_norm->GetEY()[j] / norm;
00094 gr_shift->GetY()[j] = gr_shift->GetY()[j] - norm;
00095 gr_res->GetY()[j] = gr_res->GetY()[j] / peak_search[i];
00096 gr_res->GetEY()[j] = gr_res->GetEY()[j] / peak_search[i];
00097 }
00098
00099 TFitResultPtr fit = gr->Fit("pol1", "SQ", "AL");
00100 gr_norm->Fit("pol1", "SQ", "AL");
00101 gr_shift->Fit("pol1", "SQ", "AL");
00102 sprintf(lab, "%u keV: %2f +/- %2f eV/s (%3f eV Total)", (unsigned int)peak_search[i], fit->Value(1)*1000., fit->ParError(1)*1000., total_livetime*fit->Value(1)*1000.);
00103 leg->AddEntry(gr, lab);
00104 sprintf(lab, "%u keV", (unsigned int)peak_search[i]);
00105 leg_norm->AddEntry(gr_norm, lab);
00106 leg_shift->AddEntry(gr_shift, lab);
00107 fit = gr_sig->Fit("pol1", "SQ");
00108 gr_res->Fit("pol1", "SQ");
00109 sprintf(lab, "%u keV: %2f +/- %2f eV/s", (unsigned int)peak_search[i], fit->Value(1)*1000., fit->ParError(1)*1000.);
00110 leg_sig->AddEntry(gr_sig, lab);
00111 sprintf(lab, "%u keV", (unsigned int)peak_search[i]);
00112 leg_res->AddEntry(gr_res, lab);
00113
00114 gr->SetMarkerColor(i+1); gr_norm->SetMarkerColor(i+1); gr_shift->SetMarkerColor(i+1);
00115 gr->SetMarkerStyle(21); gr_norm->SetMarkerStyle(21); gr_shift->SetMarkerStyle(21);
00116 gr->SetFillStyle(0); gr_norm->SetFillStyle(0); gr_shift->SetFillStyle(0);
00117 gr->GetFunction("pol1")->SetLineColor(i+1); gr_norm->GetFunction("pol1")->SetLineColor(i+1); gr_shift->GetFunction("pol1")->SetLineColor(i+1);
00118 gr_sig->SetMarkerColor(i+1); gr_res->SetMarkerColor(i+1);
00119 gr_sig->SetMarkerStyle(21); gr_res->SetMarkerStyle(21);
00120 gr_sig->SetFillStyle(0); gr_res->SetFillStyle(0);
00121 gr_sig->GetFunction("pol1")->SetLineColor(i+1); gr_res->GetFunction("pol1")->SetLineColor(i+1);
00122 mg->Add(gr);
00123 mg_norm->Add(gr_norm);
00124 mg_shift->Add(gr_shift);
00125 mg_sig->Add(gr_sig);
00126 mg_res->Add(gr_res);
00127 }
00128
00129
00130 new TCanvas();
00131 leg->SetTextSize(0.03);
00132 mg->Draw("AP");
00133 mg->GetXaxis()->SetTitle("Time since start of first run (s)");
00134 mg->GetYaxis()->SetTitle("Peak Position (keV)");
00135 leg->Draw("SAME");
00136 new TCanvas();
00137 leg_norm->SetTextSize(0.03);
00138 mg_norm->Draw("AP");
00139 mg_norm->GetXaxis()->SetTitle("Time since start of first run (s)");
00140 mg_norm->GetYaxis()->SetTitle("Normalized Peak Position (keV)");
00141 leg_norm->Draw("SAME");
00142 new TCanvas();
00143 leg_shift->SetTextSize(0.03);
00144 mg_shift->Draw("AP");
00145 mg_shift->GetXaxis()->SetTitle("Time since start of first run (s)");
00146 mg_shift->GetYaxis()->SetTitle("Shifted Peak Position (keV)");
00147 leg_shift->Draw("SAME");
00148 new TCanvas();
00149 leg_sig->SetTextSize(0.03);
00150 mg_sig->Draw("AP");
00151 mg_sig->GetXaxis()->SetTitle("Time since start of first run (s)");
00152 mg_sig->GetYaxis()->SetTitle("#sigma_{E}");
00153 leg_sig->Draw("SAME");
00154 new TCanvas();
00155 leg_res->SetTextSize(0.03);
00156 mg_res->Draw("AP");
00157 mg_res->GetXaxis()->SetTitle("Time since start of first run (s)");
00158 mg_res->GetYaxis()->SetTitle("#sigma_{E}/E");
00159 leg_res->Draw("SAME");
00160 }
