145 std::vector<std::vector<std::vector<std::vector<std::vector<double>>>>> cut_m;
146 std::vector<std::vector<std::vector<std::vector<std::vector<double>>>>> cut_M;
148 for (
int par = 0; par <
c_nbinpar; par++) {
149 std::vector<std::vector<std::vector<std::vector<double>>>> cut_M_par;
150 std::vector<std::vector<std::vector<std::vector<double>>>> cut_m_par;
151 for (
int lay1 = 0; lay1 <
c_nbinlay; lay1++) {
152 std::vector<std::vector<std::vector<double>>> cut_M_lay1;
153 std::vector<std::vector<std::vector<double>>> cut_m_lay1;
154 for (
int lay2 = 0; lay2 <
c_nbinlay; lay2++) {
155 std::vector<std::vector<double>> cut_M_lay2;
156 std::vector<std::vector<double>> cut_m_lay2;
157 for (
int theta = 0; theta <
c_nbint; theta++) {
158 std::vector<double> cut_M_theta;
159 std::vector<double> cut_m_theta;
160 for (
int p = 0; p <
c_nbinp; p++) {
164 TF1* fit_2gaus =
new TF1(
"fit_2gaus",
"[0]*TMath::Gaus(x,[1],[2], kTRUE)+[3]*TMath::Gaus(x,[4],[5],kTRUE)+[6]",
168 int nbin0 =
m_histo.at(par).at(lay1).at(lay2).at(theta).at(p)->GetBinContent(bin0);
169 double moltSigma1 = 0.5;
170 double moltSigma2 = 2;
171 double sigma1 = (double)moltSigma1 * (
m_histo.at(par).at(lay1).at(lay2).at(theta).at(p)->GetStdDev());
172 double sigma2 = (double)moltSigma2 * (
m_histo.at(par).at(lay1).at(lay2).at(theta).at(p)->GetStdDev());
173 double norm1 = sigma1 *
sqrt(2 * M_PI) * 0.9 * nbin0;
174 double norm2 = sigma2 *
sqrt(2 * M_PI) * 0.1 * nbin0;
175 double mean1 = (double)
m_histo.at(par).at(lay1).at(lay2).at(theta).at(p)->GetMean();
176 double mean2 = (double)
m_histo.at(par).at(lay1).at(lay2).at(theta).at(p)->GetMean();
177 double bkg = (double)
m_histo.at(par).at(lay1).at(lay2).at(theta).at(p)->GetBinContent(2);
179 fit_2gaus->SetParameters(norm1, mean1, sigma1, norm2, mean2, sigma2, bkg);
181 cut_M_theta.push_back(3 * (
sqrt(fit_2gaus->GetParameter(2)*fit_2gaus->GetParameter(2))));
182 cut_m_theta.push_back(3 * (-
sqrt(fit_2gaus->GetParameter(2)*fit_2gaus->GetParameter(2))));
188 double integral =
m_histo.at(par).at(lay1).at(lay2).at(theta).at(p)->Integral();
189 double sum_M =
m_histo.at(par).at(lay1).at(lay2).at(theta).at(p)->GetBinContent(
c_nbin + 1);
190 double sum_m =
m_histo.at(par).at(lay1).at(lay2).at(theta).at(p)->GetBinContent(0);
195 while (sum_m < integral * percent / 2) {
197 sum_m = sum_m +
m_histo.at(par).at(lay1).at(lay2).at(theta).at(p)->GetBinContent(bin_m);
199 while (sum_M < integral * percent / 2) {
201 sum_M = sum_M +
m_histo.at(par).at(lay1).at(lay2).at(theta).at(p)->GetBinContent(bin_M);
203 if (
m_histo.at(par).at(lay1).at(lay2).at(theta).at(p)->GetEntries() < 100) {
205 int filledBin_M =
c_nbin + 1;
206 while (
m_histo.at(par).at(lay1).at(lay2).at(theta).at(p)->GetBinContent(filledBin_m) == 0 && filledBin_m < (
double)
c_nbin / 2) {
209 while (
m_histo.at(par).at(lay1).at(lay2).at(theta).at(p)->GetBinContent(filledBin_M) == 0 && filledBin_M > (
double)
c_nbin / 2) {
215 cut_M_theta.push_back(
m_histo.at(par).at(lay1).at(lay2).at(theta).at(p)->GetBinCenter(bin_M));
216 cut_m_theta.push_back(
m_histo.at(par).at(lay1).at(lay2).at(theta).at(p)->GetBinCenter(bin_m));
221 cut_M_lay2.push_back(cut_M_theta);
223 cut_m_lay2.push_back(cut_m_theta);
226 cut_M_lay1.push_back(cut_M_lay2);
228 cut_m_lay1.push_back(cut_m_lay2);
231 cut_M_par.push_back(cut_M_lay1);
233 cut_m_par.push_back(cut_m_lay1);
236 cut_M.push_back(cut_M_par);
238 cut_m.push_back(cut_m_par);
245 std::vector<std::vector<std::vector<TH2F*>>> cut_M_histo;
246 std::vector<std::vector<std::vector<TH2F*>>> cut_m_histo;
248 for (
int par = 0; par <
c_nbinpar; par++) {
249 std::vector<std::vector<TH2F*>> cut_M_histo_par;
250 std::vector<std::vector<TH2F*>> cut_m_histo_par;
251 for (
int lay1 = 0; lay1 <
c_nbinlay; lay1++) {
252 std::vector<TH2F*> cut_M_histo_lay1;
253 std::vector<TH2F*> cut_m_histo_lay1;
254 for (
int lay2 = 0; lay2 <
c_nbinlay; lay2++) {
255 cut_M_histo_lay1.push_back(
new TH2F(
"CUTS_M_" +
m_namePar.at(par) + Form(
"layer%d_%d", lay1, lay2),
257 cut_m_histo_lay1.push_back(
new TH2F(
"CUTS_m_" +
m_namePar.at(par) + Form(
"layer%d_%d", lay1, lay2),
259 for (
int theta = 1; theta <=
c_nbint; theta++) {
260 for (
int p = 1; p <=
c_nbinp; p++) {
261 cut_M_histo_lay1.at(lay2)->SetBinContent(p, theta, cut_M.at(par).at(lay1).at(lay2).at(theta - 1).at(p - 1));
262 cut_m_histo_lay1.at(lay2)->SetBinContent(p, theta, cut_m.at(par).at(lay1).at(lay2).at(theta - 1).at(p - 1));
266 cut_M_histo_par.push_back(cut_M_histo_lay1);
267 cut_M_histo_lay1.clear();
268 cut_m_histo_par.push_back(cut_m_histo_lay1);
269 cut_m_histo_lay1.clear();
271 cut_M_histo.push_back(cut_M_histo_par);
272 cut_M_histo_par.clear();
273 cut_m_histo.push_back(cut_m_histo_par);
274 cut_m_histo_par.clear();
278 std::vector<std::vector<std::vector<std::vector<double>>>> cut_out_norm;
279 std::vector<std::vector<std::vector<std::vector<double>>>> cut_out_pow;
280 std::vector<std::vector<std::vector<std::vector<double>>>> cut_out_bkg;
282 for (
int minmax = 0; minmax < 2; minmax++) {
283 std::vector<std::vector<std::vector<double>>> cut_out_norm_minmax;
284 std::vector<std::vector<std::vector<double>>> cut_out_pow_minmax;
285 std::vector<std::vector<std::vector<double>>> cut_out_bkg_minmax;
286 for (
int par = 0; par <
c_nbinpar; par++) {
287 std::vector<std::vector<double>> cut_out_norm_par;
288 std::vector<std::vector<double>> cut_out_pow_par;
289 std::vector<std::vector<double>> cut_out_bkg_par;
290 for (
int lay1 = 0; lay1 <
c_nbinlay; lay1++) {
291 std::vector<double> cut_out_norm_lay1;
292 std::vector<double> cut_out_pow_lay1;
293 std::vector<double> cut_out_bkg_lay1;
294 for (
int lay2 = 0; lay2 <
c_nbinlay; lay2++) {
296 TF2* fit_MS =
new TF2(
"fit_MS",
"[0]*1/(TMath::Power(x,[1])*TMath::Sqrt(TMath::Sin(y)))+[2]",
c_pmin,
c_pmax,
c_tmin,
c_tmax);
298 double norm = cut_M.at(par).at(lay1).at(lay2).at(1).at(1);
300 double bkg = cut_M.at(par).at(lay1).at(lay2).at(1).at(
c_nbinp - 1);
301 fit_MS->SetParameters(norm, Pow, bkg);
304 cut_m_histo.at(par).at(lay1).at(lay2) -> Fit(
"fit_MS");
308 cut_M_histo.at(par).at(lay1).at(lay2) -> Fit(
"fit_MS");
311 cut_out_norm_lay1.push_back(fit_MS->GetParameter(0));
312 cut_out_pow_lay1.push_back(fit_MS->GetParameter(1));
313 cut_out_bkg_lay1.push_back(fit_MS->GetParameter(2));
315 cut_out_norm_par.push_back(cut_out_norm_lay1);
316 cut_out_norm_lay1.clear();
317 cut_out_pow_par.push_back(cut_out_pow_lay1);
318 cut_out_pow_lay1.clear();
319 cut_out_bkg_par.push_back(cut_out_bkg_lay1);
320 cut_out_bkg_lay1.clear();
322 cut_out_norm_minmax.push_back(cut_out_norm_par);
323 cut_out_norm_par.clear();
324 cut_out_pow_minmax.push_back(cut_out_pow_par);
325 cut_out_pow_par.clear();
326 cut_out_bkg_minmax.push_back(cut_out_bkg_par);
327 cut_out_bkg_par.clear();
329 cut_out_norm.push_back(cut_out_norm_minmax);
330 cut_out_norm_minmax.clear();
331 cut_out_pow.push_back(cut_out_pow_minmax);
332 cut_out_pow_minmax.clear();
333 cut_out_bkg.push_back(cut_out_bkg_minmax);
334 cut_out_bkg_minmax.clear();
341 for (
int g = 0; g <
c_nbinp; g++) {
344 B2DEBUG(30,
"momentum=" << p_out_mom);
345 B2DEBUG(30,
"d0, 3-4, Min: " << cut_m.at(1).at(3).at(4).at(1).at(g));
346 B2DEBUG(30,
"min cut (TH2F):" << cut_m_histo.at(1).at(3).at(4)->GetBinContent(g + 1, 2));
347 double norm_min = cut_out_norm.at(0).at(1).at(3).at(4);
348 double pow_min = cut_out_pow.at(0).at(1).at(3).at(4);
349 double bkg_min = cut_out_bkg.at(0).at(1).at(3).at(4);
350 B2DEBUG(30,
"norm par min:" << norm_min);
351 B2DEBUG(30,
"pow par min:" << pow_min);
352 B2DEBUG(30,
"bkg par min:" << bkg_min);
353 B2DEBUG(30,
"evaluate min cut:" << norm_min / (
sqrt(sin(t_out_theta)) * pow(p_out_mom, pow_min)) + bkg_min);
354 B2DEBUG(30,
"d0, 3-4, Max: " << cut_M.at(1).at(3).at(4).at(1).at(g));
355 B2DEBUG(30,
"max cut (TH2F):" << cut_M_histo.at(1).at(3).at(4)->GetBinContent(g + 1, 2));
356 double norm_max = cut_out_norm.at(1).at(1).at(3).at(4);
357 double pow_max = cut_out_pow.at(1).at(1).at(3).at(4);
358 double bkg_max = cut_out_bkg.at(1).at(1).at(3).at(4);
359 B2DEBUG(30,
"norm par max:" << norm_max);
360 B2DEBUG(30,
"pow par max:" << pow_max);
361 B2DEBUG(30,
"bkg par max:" << bkg_max);
362 B2DEBUG(30,
"evaluate max cut:" << norm_max / (
sqrt(sin(t_out_theta)) * pow(p_out_mom, pow_max)) + bkg_max);
363 B2DEBUG(30,
"----------------------------------------");
371 for (
int par = 0; par <
c_nbinpar; par++) {
372 for (
int lay1 = 0; lay1 <
c_nbinlay; lay1++) {
373 for (
int lay2 = 0; lay2 <
c_nbinlay; lay2++) {
374 for (
int theta = 0; theta <
c_nbint; theta++) {
375 for (
int p = 0; p <
c_nbinp; p++) {
376 double layerdiff = lay2 - lay1;
377 if (layerdiff >= 0 && layerdiff < 3) {
378 if (
m_histo.at(par).at(lay1).at(lay2).at(theta).at(p)->GetEntries() > 0) {
379 m_histo.at(par).at(lay1).at(lay2).at(theta).at(p)->Write();
389 for (
int par = 0; par <
c_nbinpar; par++) {
390 for (
int lay1 = 0; lay1 <
c_nbinlay; lay1++) {
391 for (
int lay2 = 0; lay2 <
c_nbinlay; lay2++) {
392 for (
int minmax = 0; minmax < 2; minmax++) {
394 int layerdiff = lay2 - lay1;
395 if (layerdiff >= 0 && layerdiff < 3) {
396 cut_m_histo.at(par).at(lay1).at(lay2)->Write();
400 int layerdiff = lay2 - lay1;
401 if (layerdiff >= 0 && layerdiff < 3) {
402 cut_M_histo.at(par).at(lay1).at(lay2)->Write();
413 TH3F* output_norm_m =
new TH3F(
"output_norm_m",
"output_norm_m",
c_nbinpar, 0, 4,
c_nbinlay, 0, 4,
c_nbinlay, 0, 4);
414 TH3F* output_pow_m =
new TH3F(
"output_pow_m",
"output_pow_m",
c_nbinpar, 0, 4,
c_nbinlay, 0, 4,
c_nbinlay, 0, 4);
415 TH3F* output_bkg_m =
new TH3F(
"output_bkg_m",
"output_bkg_m",
c_nbinpar, 0, 4,
c_nbinlay, 0, 4,
c_nbinlay, 0, 4);
417 TH3F* output_norm_M =
new TH3F(
"output_norm_M",
"output_norm_M",
c_nbinpar, 0, 4,
c_nbinlay, 0, 4,
c_nbinlay, 0, 4);
418 TH3F* output_pow_M =
new TH3F(
"output_pow_M",
"output_pow_M",
c_nbinpar, 0, 4,
c_nbinlay, 0, 4,
c_nbinlay, 0, 4);
419 TH3F* output_bkg_M =
new TH3F(
"output_bkg_M",
"output_bkg_M",
c_nbinpar, 0, 4,
c_nbinlay, 0, 4,
c_nbinlay, 0, 4);
422 for (
int par = 0; par <
c_nbinpar; par++) {
423 for (
int lay1 = 0; lay1 <
c_nbinlay; lay1++) {
424 for (
int lay2 = 0; lay2 <
c_nbinlay; lay2++) {
425 output_norm_m->SetBinContent(par, lay1, lay2, cut_out_norm.at(0).at(par).at(lay1).at(lay2));
426 output_norm_M->SetBinContent(par, lay1, lay2, cut_out_norm.at(1).at(par).at(lay1).at(lay2));
428 output_pow_m->SetBinContent(par, lay1, lay2, cut_out_pow.at(0).at(par).at(lay1).at(lay2));
429 output_pow_M->SetBinContent(par, lay1, lay2, cut_out_pow.at(1).at(par).at(lay1).at(lay2));
431 output_bkg_m->SetBinContent(par, lay1, lay2, cut_out_bkg.at(0).at(par).at(lay1).at(lay2));
432 output_bkg_M->SetBinContent(par, lay1, lay2, cut_out_bkg.at(1).at(par).at(lay1).at(lay2));
438 output_norm_m->Write();
439 output_norm_M->Write();
440 output_pow_m->Write();
441 output_pow_M->Write();
442 output_bkg_m->Write();
443 output_bkg_M->Write();
447 B2INFO(
"number of spacepoint with theta out of limits=" <<
m_tCounter);
448 B2INFO(
"number of spacepoint with momentum out of limits=" <<
m_pCounter);
449 B2INFO(
"number of tracks cut by global cuts=" <<
m_globCounter);