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