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 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 } else if (fitResKFactorMethod == 3) {
234 TF1* kgsFitFcn =
new TF1(
"kgsFitFcn",
scaleFactorFitFcn, minFitRange, maxFitRange, 3);
235 kgsFitFcn->SetParameter(0, minFitRange);
236 kgsFitFcn->SetParameter(1, 1.0);
237 kgsFitFcn->SetParameter(2, 1.0);
238 for (int ifit = 0; ifit < 1; ifit++)
239 kgs_z_vs_ntrk->Fit(kgsFitFcn, "Q");
240 kgs_z_vs_ntrk->Fit(kgsFitFcn, "Q");
241 kgs_z_ntrk_fit = kgsFitFcn;
242 kgs_z_ntrk_fit_er = kgsFitFcn->GetParErrors();
243
244
245
246
247 } else if (fitResKFactorMethod == 4) {
248
249 kgs_z_vs_ntrk->Fit("pol0", "Q", "", minFitRange, maxFitRange);
250 kgs_z_vs_ntrk->GetFunction("pol0")->SetLineColor(kRed);
251 kgs_z_ntrk_fit = kgs_z_vs_ntrk->GetFunction("pol0");
252 kgs_z_ntrk_fit_er = kgs_z_ntrk_fit->GetParErrors();
253
254
255 }
256
257
258 int nbins_z_err_ntrk = h_Vrt_err_vs_Something->GetNbinsX();
259
260 std::vector<float> av_err_z;
261 std::vector<float> av_err_z_er;
262
263
264 std::vector<float> err_bins_z_nt;
265 std::vector<float> err_bins_z_nt_er;
266 std::vector<float> res_z;
267 std::vector<float> res_z_er;
268
269
270
271 for (int bin_count = 1; bin_count <= nbins_z_err_ntrk; ++bin_count) {
272 err_bins_z_nt.push_back(h_Vrt_err_vs_Something->GetXaxis()->GetBinCenter(bin_count));
273 err_bins_z_nt_er.push_back(h_Vrt_err_vs_Something->GetXaxis()->GetBinWidth(bin_count) / 2.);
274
275 TH1D* profileY = h_Vrt_err_vs_Something->ProjectionY("projectionErrors", bin_count, bin_count, "e");
276
277
278
279
280 float mean = profileY->GetMean();
281 float mean_error = profileY->GetMeanError();
282
283
284
285
286
287
288 delete profileY;
289
290
291 av_err_z.push_back(
mean);
292 av_err_z_er.push_back(mean_error);
293
294
295
296
297 double pr_er = 0.0;
299 if (fitResKFactorMethod == 1) {
300 pr_er =
error_func(bin_count, kgs_z_ntrk_fit_er);
301 } else if (fitResKFactorMethod == 2) {
302 val = h_Vrt_err_vs_Something->GetXaxis()->GetBinCenter(bin_count);
303 pr_er = TMath::Power(kgs_z_ntrk_fit_er[1] * val, 2) + TMath::Power(kgs_z_ntrk_fit_er[0], 2);
304 pr_er = TMath::Sqrt(pr_er);
305
306
307 } else if (fitResKFactorMethod == 3) {
308 val = h_Vrt_err_vs_Something->GetXaxis()->GetBinCenter(bin_count);
309
310 pr_er = kgs_z_ntrk_fit_er[2];
311 } else if (fitResKFactorMethod == 4) {
312 pr_er = kgs_z_ntrk_fit_er[0];
313 }
314
315 res_z.push_back(
mean * kgs_z_ntrk_fit->Eval(h_Vrt_err_vs_Something->GetXaxis()->GetBinCenter(bin_count)));
316 res_z_er.push_back(TMath::Sqrt(TMath::Power(mean_error *
317 kgs_z_ntrk_fit->Eval(h_Vrt_err_vs_Something->GetXaxis()->GetBinCenter(
318 bin_count)),
319 2) + TMath::Power(pr_er *
mean, 2)));
320
321
322
323
324 }
325 TGraphErrors* res_z_vs_ntrk =
326 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]));
327 res_z_vs_ntrk->GetYaxis()->SetTitle(coordinate + " Vertex Resolution [mm]");
328 res_z_vs_ntrk->GetXaxis()->SetTitle(xAxisLabel);
329 res_z_vs_ntrk->SetTitle(coordinate + " Vertex Resolution");
330 res_z_vs_ntrk->SetName("resolution_" + coordinate + "_" + versus);
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
347
349 if (versus == "Ntrk") res_z_vs_ntrk->GetXaxis()->SetRangeUser(0., 100.);
350 else res_z_vs_ntrk->GetXaxis()->SetRangeUser(0., 20.);
351 res_z_vs_ntrk->GetYaxis()->SetRangeUser(0.0025, 1.);
352 res_z_vs_ntrk->Write("", TObject::kOverwrite);
353 delete res_z_vs_ntrk;
354
355 if (versus == "Ntrk") krms_z_vs_ntrk->GetXaxis()->SetRangeUser(0., 100.);
356 else krms_z_vs_ntrk->GetXaxis()->SetRangeUser(0., 20.);
357 krms_z_vs_ntrk->GetYaxis()->SetRangeUser(0.5, 1.3);
358 krms_z_vs_ntrk->Write("", TObject::kOverwrite);
359 delete krms_z_vs_ntrk;
360
361 if (versus == "Ntrk") kgs_z_vs_ntrk->GetXaxis()->SetRangeUser(0., 100.);
362 else kgs_z_vs_ntrk->GetXaxis()->SetRangeUser(0., 20.);
363 kgs_z_vs_ntrk->GetYaxis()->SetRangeUser(0.5, 1.3);
364 kgs_z_vs_ntrk->Write("", TObject::kOverwrite);
365 delete kgs_z_vs_ntrk;
366
367
368
369
370
371 return;
372 }
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)