10 #include "TFitResult.h"
21 namespace ExtractTimeOffsetsAndCEffNS {
27 TH1I *
Get1DHistogram(
const char * plugin,
const char * directoryName,
const char * name,
bool print =
true){
29 TString fullName = TString(plugin) +
"/" + TString(directoryName) +
"/" + TString(name);
30 thisFile->GetObject(fullName, histogram);
32 if (print) cout <<
"Unable to find histogram " << fullName.Data() << endl;
40 TString fullName = TString(name);
41 thisFile->GetObject(fullName, histogram);
43 if (print) cout <<
"Unable to find histogram " << fullName.Data() << endl;
49 TH2I *
Get2DHistogram(
const char * plugin,
const char * directoryName,
const char * name){
51 TString fullName = TString(plugin) +
"/" + TString(directoryName) +
"/" + TString(name);
52 thisFile->GetObject(fullName, histogram);
54 cout <<
"Unable to find histogram " << fullName.Data() << endl;
67 printf(
"ExtractTimeOffsets\n");
71 char outfilename[255];
72 sprintf(outfilename,
"output/timing/rootfiles/BCALTimeOffset_Results_%s_%i.root",tag.Data(),
run);
73 TFile *outputFile = TFile::Open(outfilename,
"RECREATE");
74 outputFile->mkdir(
"Fits");
75 outputFile->mkdir(
"ResultOverview");
81 cout <<
"Unable to open file " << filename.Data() <<
"...Exiting" << endl;
86 ofstream channel_global_offset_File;
87 sprintf(outfilename,
"output/timing/correction/%s/channel_offset_%i.txt",tag.Data(),
run);
88 channel_global_offset_File.open(outfilename);
90 TH1I *OldTimeOffsets =
new TH1I(
"OldTimeOffsets",
"Total time offset (pre)", 2000, -4, 4);
91 TH1F *OldTimeOffsetsVsChannel =
new TH1F(
"OldTimeOffsetsVsChannel",
"Total time offset (pre) vs channel;CCDB channel", 768, 0.5, 768.5);
92 TH1I *NewTimeOffsets =
new TH1I(
"NewTimeOffsets",
"Total time offset (post)", 2000, -4, 4);
93 TH1F *NewTimeOffsetsVsChannel =
new TH1F(
"NewTimeOffsetsVsChannel",
"Total time offset (post) vs channel;CCDB channel", 768, 0.5, 768.5);
94 TH1I *NewTimeOffsetsFit =
new TH1I(
"NewTimeOffsetsFit",
"Total time offset fit (post)", 2000, -4, 4);
95 TH1F *NewTimeOffsetsFitVsChannel =
new TH1F(
"NewTimeOffsetsFitVsChannel",
"Total time offset fit (post) vs channel;CCDB channel", 768, 0.5, 768.5);
96 TH1I *RelativeTimeOffsetsFit =
new TH1I(
"RelativeTimeOffsetsFit",
"Change in time offset", 2000, -4, 4);
97 TH1F *RelativeTimeOffsetsFitVsChannel =
new TH1F(
"RelativeTimeOffsetsFitVsChannel",
"Change in time offset vs channel;CCDB channel", 768, 0.5, 768.5);
99 TH1D *priorOffsetHist = NULL;
103 if (priorOffsetHist != NULL) {
108 printf(
"Failed to load %s\n",
"CCDB_raw_channel_global_offset");
111 if(residualHist != NULL && priorOffsetHist != NULL){
113 double NumberOfFiles = priorOffsetHist->GetBinContent(769);
114 if (NumberOfFiles!=1) priorOffsetHist->Scale(1/NumberOfFiles);
116 int nBinsX = residualHist->GetNbinsX();
117 int nBinsY = residualHist->GetNbinsY();
119 TString newtag = tag +
"_channel";
121 TH1D *meanhist =
GetGausFitMean(residualHist,tag.Data(),makeplots);
123 for (
int i = 1 ; i <= nBinsX; i++){
125 TH1D *projY = residualHist->ProjectionY(
"temp", i, i);
127 float nsPerBin = (projY->GetBinCenter(projY->GetNbinsX()) - projY->GetBinCenter(1)) / projY->GetNbinsX();
128 float timeWindow = 0.5;
129 int binWindow = int(timeWindow / nsPerBin);
130 double maxEntries = 0;
132 for (
int j = 1 ; j <= projY->GetNbinsX();j++){
134 int maxBin = (j + binWindow) <= projY->GetNbinsX() ? (j + binWindow) : projY->GetNbinsX();
135 double sum = 0, nEntries = 0;
136 for (
int bin = minBin; bin <= maxBin; bin++){
137 sum += projY->GetBinContent(bin) * projY->GetBinCenter(bin);
138 nEntries += projY->GetBinContent(bin);
140 if (nEntries > maxEntries) {
141 maxMean = sum / nEntries;
142 maxEntries = nEntries;
147 float priorOffset = priorOffsetHist->GetBinContent(i);
148 float residualOffset = maxMean;
149 float residualOffsetFit = meanhist->GetBinContent(i);
150 float residualOffsetFitErr = meanhist->GetBinError(i);
151 float newOffset = priorOffset + residualOffset;
152 float newOffsetFit = priorOffset + residualOffsetFit;
154 OldTimeOffsets->Fill(priorOffset);
155 OldTimeOffsetsVsChannel->SetBinContent(i,priorOffset);
156 NewTimeOffsets->Fill(newOffset);
157 NewTimeOffsetsVsChannel->SetBinContent(i,newOffset);
158 NewTimeOffsetsFit->Fill(newOffsetFit);
159 NewTimeOffsetsFitVsChannel->SetBinContent(i,newOffsetFit);
160 RelativeTimeOffsetsFit->Fill(residualOffsetFit);
161 RelativeTimeOffsetsFitVsChannel->SetBinContent(i,residualOffsetFit);
163 if (residualOffsetFitErr<0.1) {
164 sprintf(outfilename,
"output/timing/channels_fit/%s_%03i.txt",tag.Data(),chanint);
165 outfile = fopen(outfilename,
"a");
167 fprintf(outfile,
"%i %.3f %.3f\n",run,newOffsetFit,residualOffsetFitErr);
171 sprintf(outfilename,
"output/timing/channels_roll/%s_%03i.txt",tag.Data(),chanint);
172 outfile = fopen(outfilename,
"a");
174 fprintf(outfile,
"%i %.3f\n",run,newOffset);
177 printf(
"%4i rolling max: %6.3f %+.3f = %6.3f gaus fit: %6.3f %+.3f = %6.3f diff: %6.3f\n",
178 i,priorOffset,residualOffset,newOffset,priorOffset,residualOffsetFit,newOffsetFit,residualOffset-residualOffsetFit);
179 channel_global_offset_File << newOffsetFit << endl;
182 channel_global_offset_File.close();
190 double outermintime=-2, outermaxtime=2;
192 TAxis * xaxis = hist2d->GetXaxis();
193 TAxis * yaxis = hist2d->GetYaxis();
194 int nBinsX = hist2d->GetNbinsX();
195 int nBinsY = hist2d->GetNbinsY();
197 printf(
"bins (%i,%i)\n",nBinsX,nBinsY);
199 TH1D *xmean = hist2d->ProjectionX(
"xmean");
202 gStyle->SetOptStat(2210);
203 gStyle->SetOptFit(1);
206 canvas_fit =
new TCanvas(
"canvas_fit",
"canvas_fit");
209 for (
int i = 1 ; i <= nBinsX; i++){
210 TH1D *singlebin_projY = hist2d->ProjectionY(
"temp", i, i);
211 float center = yaxis->GetBinCenter(singlebin_projY->GetMaximumBin());
213 float mintime = center-width;
214 float maxtime = center+width;
216 TF1 *f_gaus =
new TF1(
"f_gaus",
"gaus",mintime,maxtime);
219 singlebin_projY->SetAxisRange(mintime,maxtime);
220 double mean=0, meanerr=0;
222 singlebin_projY->Fit(f_gaus,
"RQ");
223 mean = f_gaus->GetParameter(1);
224 meanerr = f_gaus->GetParError(1);
225 xmean->SetBinContent(i,mean);
226 xmean->SetBinError(i,meanerr);
227 if (makeplots || meanerr>0.5) {
229 sprintf(plotname,
"plots/gausfitformean/fit_%s_%i.png",tag.Data(),i);
230 canvas_fit->Print(plotname);
232 singlebin_projY->Reset();
sprintf(text,"Post KinFit Cut")
TH1I * Get1DHistogramBasedir(const char *name, bool print=true)
TH2I * Get2DHistogram(const char *plugin, const char *directoryName, const char *name)
TH1I * Get1DHistogram(const char *plugin, const char *directoryName, const char *name, bool print=true)
printf("string=%s", string)