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