106 {
107
108
109 TH2F* h_Vrt_pullVsSomething_split(0);
110 TH2F* h_Vrt_err_vs_Something(0);
111
112 std::string xAxisLabel("");
113
114 if (versus == "Ntrk") {
115 h_Vrt_pullVsSomething_split = (TH2F*) gDirectory->Get(("Vrt_" + coordinate + "pullVsNtrkAverage_split").c_str());
116 h_Vrt_err_vs_Something = (TH2F*) gDirectory->Get(("Vrt_" + coordinate + "err_vs_ntrk").c_str());
117
118 xAxisLabel = "Number of fitted tracks";
119 } else if (versus == "SumPt2") {
120 h_Vrt_pullVsSomething_split = (TH2F*) gDirectory->Get(("Vrt_" + coordinate + "pullVsPt2Average_split").c_str());
121 h_Vrt_err_vs_Something = (TH2F*) gDirectory->Get(("Vrt_" + coordinate + "err_vs_pt2").c_str());
122
123 xAxisLabel = "#sqrt{#sum p_{T}^{2}} [GeV]";
124 } else return;
125
126
127
128 if (h_Vrt_pullVsSomething_split == 0 or h_Vrt_err_vs_Something == 0) return;
129
130 int n_bins = h_Vrt_pullVsSomething_split->GetNbinsX();
131 std::vector<float> rms_z;
132 std::vector<float> rms_z_er;
133 std::vector<float> sigma_z;
134 std::vector<float> sigma_z_er;
135 std::vector<float> bins_z_nt;
136 std::vector<float> bins_z_nt_er;
137
138
139
140
141
142
143
144 Int_t startBin = 0;
145 TH1D* profileZ = 0;
146 const Int_t minEntriesForKFactorBin = 1000;
147 for (int bin_count = 1; bin_count < n_bins + 1; bin_count++) {
148
149
150
151
152
153
154 TH1D* profileZTmp = h_Vrt_pullVsSomething_split->ProjectionY("projectionPulls", bin_count, bin_count, "e");
155
156 if (profileZ == 0) {
157 startBin = bin_count;
158 profileZ = (TH1D*) profileZTmp->Clone("projectionPulls_Integrated");
159
160 } else {
161 profileZ->Add(profileZTmp);
162 }
163 delete profileZTmp;
164 profileZTmp = 0;
165 if ((profileZ->GetEntries() < minEntriesForKFactorBin) && (bin_count < n_bins))
166 continue;
167
168 Double_t lowEdge = h_Vrt_pullVsSomething_split->GetXaxis()->GetBinLowEdge(startBin);
169 Double_t highEdge = h_Vrt_pullVsSomething_split->GetXaxis()->GetBinLowEdge(bin_count) +
170 h_Vrt_pullVsSomething_split->GetXaxis()->GetBinWidth(bin_count);
171 Double_t binCenter = (lowEdge + highEdge) / 2;
172 Double_t
binWidth = (highEdge - lowEdge) / 2;
173
174
175
176 bins_z_nt.push_back(binCenter);
178
179 rms_z.push_back(profileZ->GetRMS());
180 rms_z_er.push_back(profileZ->GetRMSError());
181
182
183 if (profileZ->GetEntries() > 100.) {
185 sigma_z.push_back(fit_res[0]);
186 sigma_z_er.push_back(fit_res[1]);
187 } else {
188 sigma_z.push_back(0.);
189 sigma_z_er.push_back(0.);
190 }
191
192 delete profileZ;
193 profileZ = 0;
194 }
195
196 TGraphErrors* krms_z_vs_ntrk = new TGraphErrors(
197 bins_z_nt.size(), &(bins_z_nt[0]), &(rms_z[0]), &(bins_z_nt_er[0]), &(rms_z_er[0]));
198 krms_z_vs_ntrk->GetYaxis()->SetTitle((coordinate + " scale factor from RMS").c_str());
199 krms_z_vs_ntrk->GetXaxis()->SetTitle(xAxisLabel.c_str());
200 krms_z_vs_ntrk->SetTitle(("scaleFactor" + coordinate + "_RMS").c_str());
201 krms_z_vs_ntrk->SetName(("scaleFactor" + coordinate + "_" + versus + "_RMS").c_str());
202
203 TGraphErrors* kgs_z_vs_ntrk = new TGraphErrors(
204 bins_z_nt.size(), &(bins_z_nt[0]), &(sigma_z[0]), &(bins_z_nt_er[0]), &(sigma_z_er[0]));
205 kgs_z_vs_ntrk->GetYaxis()->SetTitle((coordinate + " scale factor from gauss fit").c_str());
206 kgs_z_vs_ntrk->GetXaxis()->SetTitle(xAxisLabel.c_str());
207 kgs_z_vs_ntrk->SetTitle(("scaleFactor" + coordinate + "_Fit").c_str());
208 kgs_z_vs_ntrk->SetName(("scaleFactor_" + coordinate + "_" + versus + "_Fit").c_str());
209
210
211 float maxFitRange(100.);
212 float minFitRange(2.);
213 if (versus == "SumPt2") {
214 minFitRange = 0.5;
215 maxFitRange = 20.;
216 }
217 TF1* kgs_z_ntrk_fit;
218 const Double_t* kgs_z_ntrk_fit_er;
219 int fitResKFactorMethod = 2;
220 if (fitResKFactorMethod == 1) {
221
222
223 kgs_z_vs_ntrk->Fit("pol2", "Q", "", minFitRange, maxFitRange);
224 kgs_z_vs_ntrk->GetFunction("pol2")->SetLineColor(kRed);
225 kgs_z_ntrk_fit = kgs_z_vs_ntrk->GetFunction("pol2");
226 kgs_z_ntrk_fit_er = kgs_z_ntrk_fit->GetParErrors();
227 } else if (fitResKFactorMethod == 2) {
228
229 kgs_z_vs_ntrk->Fit("pol1", "Q", "", minFitRange, maxFitRange);
230 kgs_z_vs_ntrk->GetFunction("pol1")->SetLineColor(kRed);
231 kgs_z_ntrk_fit = kgs_z_vs_ntrk->GetFunction("pol1");
232 kgs_z_ntrk_fit_er = kgs_z_ntrk_fit->GetParErrors();
233
234 } else if (fitResKFactorMethod == 3) {
235 TF1* kgsFitFcn =
new TF1(
"kgsFitFcn",
scaleFactorFitFcn, minFitRange, maxFitRange, 3);
236 kgsFitFcn->SetParameter(0, minFitRange);
237 kgsFitFcn->SetParameter(1, 1.0);
238 kgsFitFcn->SetParameter(2, 1.0);
239 for (int ifit = 0; ifit < 1; ifit++)
240 kgs_z_vs_ntrk->Fit(kgsFitFcn, "Q");
241 kgs_z_vs_ntrk->Fit(kgsFitFcn, "Q");
242 kgs_z_ntrk_fit = kgsFitFcn;
243 kgs_z_ntrk_fit_er = kgsFitFcn->GetParErrors();
244
245
246
247
248 } else if (fitResKFactorMethod == 4) {
249
250 kgs_z_vs_ntrk->Fit("pol0", "Q", "", minFitRange, maxFitRange);
251 kgs_z_vs_ntrk->GetFunction("pol0")->SetLineColor(kRed);
252 kgs_z_ntrk_fit = kgs_z_vs_ntrk->GetFunction("pol0");
253 kgs_z_ntrk_fit_er = kgs_z_ntrk_fit->GetParErrors();
254
255
256 }
257
258
259 int nbins_z_err_ntrk = h_Vrt_err_vs_Something->GetNbinsX();
260
261 std::vector<float> av_err_z;
262 std::vector<float> av_err_z_er;
263
264
265 std::vector<float> err_bins_z_nt;
266 std::vector<float> err_bins_z_nt_er;
267 std::vector<float> res_z;
268 std::vector<float> res_z_er;
269
270
271
272 for (int bin_count = 1; bin_count <= nbins_z_err_ntrk; ++bin_count) {
273 err_bins_z_nt.push_back(h_Vrt_err_vs_Something->GetXaxis()->GetBinCenter(bin_count));
274 err_bins_z_nt_er.push_back(h_Vrt_err_vs_Something->GetXaxis()->GetBinWidth(bin_count) / 2.);
275
276 TH1D* profileY = h_Vrt_err_vs_Something->ProjectionY("projectionErrors", bin_count, bin_count, "e");
277
278
279
280
281 float mean = profileY->GetMean();
282 float mean_error = profileY->GetMeanError();
283
284
285
286
287
288
289 delete profileY;
290
291
292 av_err_z.push_back(
mean);
293 av_err_z_er.push_back(mean_error);
294
295
296
297
298 double pr_er = 0.0;
300 if (fitResKFactorMethod == 1) {
301
302 pr_er =
error_func(bin_count, kgs_z_ntrk_fit_er);
303 } else if (fitResKFactorMethod == 2) {
304 val = h_Vrt_err_vs_Something->GetXaxis()->GetBinCenter(bin_count);
305 pr_er = TMath::Power(kgs_z_ntrk_fit_er[1] * val, 2) + TMath::Power(kgs_z_ntrk_fit_er[0], 2);
306 pr_er = TMath::Sqrt(pr_er);
307
308
309
310 } else if (fitResKFactorMethod == 3) {
311 val = h_Vrt_err_vs_Something->GetXaxis()->GetBinCenter(bin_count);
312
313 pr_er = kgs_z_ntrk_fit_er[2];
314 } else if (fitResKFactorMethod == 4) {
315 pr_er = kgs_z_ntrk_fit_er[0];
316 }
317
318 res_z.push_back(
mean * kgs_z_ntrk_fit->Eval(h_Vrt_err_vs_Something->GetXaxis()->GetBinCenter(bin_count)));
319 res_z_er.push_back(TMath::Sqrt(TMath::Power(mean_error *
320 kgs_z_ntrk_fit->Eval(h_Vrt_err_vs_Something->GetXaxis()->GetBinCenter(
321 bin_count)),
322 2) + TMath::Power(pr_er *
mean, 2)));
323
324
325
326
327 }
328 TGraphErrors* res_z_vs_ntrk =
329 new TGraphErrors(err_bins_z_nt.size(), &(err_bins_z_nt[0]), &(res_z[0]), &(err_bins_z_nt_er[0]), &(res_z_er[0]));
330 res_z_vs_ntrk->GetYaxis()->SetTitle((coordinate + " Vertex Resolution [mm]").c_str());
331 res_z_vs_ntrk->GetXaxis()->SetTitle(xAxisLabel.c_str());
332 res_z_vs_ntrk->SetTitle((coordinate + " Vertex Resolution").c_str());
333 res_z_vs_ntrk->SetName(("resolution_" + coordinate + "_" + versus).c_str());
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
350
352 if (versus == "Ntrk") res_z_vs_ntrk->GetXaxis()->SetRangeUser(0., 100.);
353 else res_z_vs_ntrk->GetXaxis()->SetRangeUser(0., 20.);
354 res_z_vs_ntrk->GetYaxis()->SetRangeUser(0.0025, 1.);
355 res_z_vs_ntrk->Write("", TObject::kOverwrite);
356 delete res_z_vs_ntrk;
357
358 if (versus == "Ntrk") krms_z_vs_ntrk->GetXaxis()->SetRangeUser(0., 100.);
359 else krms_z_vs_ntrk->GetXaxis()->SetRangeUser(0., 20.);
360 krms_z_vs_ntrk->GetYaxis()->SetRangeUser(0.5, 1.3);
361 krms_z_vs_ntrk->Write("", TObject::kOverwrite);
362 delete krms_z_vs_ntrk;
363
364 if (versus == "Ntrk") kgs_z_vs_ntrk->GetXaxis()->SetRangeUser(0., 100.);
365 else kgs_z_vs_ntrk->GetXaxis()->SetRangeUser(0., 20.);
366 kgs_z_vs_ntrk->GetYaxis()->SetRangeUser(0.5, 1.3);
367 kgs_z_vs_ntrk->Write("", TObject::kOverwrite);
368 delete kgs_z_vs_ntrk;
369
370
371
372
373
374 return;
375 }
void mean(std::vector< double > &bins, std::vector< double > &values, const std::vector< std::string > &files, const std::string &histname, const std::string &tplotname, const std::string &label="")
std::vector< float > stableGaussianFit(TH1 *histo)
double error_func(float x, const Double_t *par)
double scaleFactorFitFcn(double *x, double *par)