210 for (
int TCId = 1; TCId < 577; TCId++) {
212 mean_FWD +=
TCID[TCId - 1];
213 }
else if (TCId < 513) {
214 mean_BAR +=
TCID[TCId - 1];
216 mean_BWD +=
TCID[TCId - 1];
220 mean_BAR /= (512 - 80);
221 mean_BWD /= (576 - 512);
223 for (
int TCId = 1; TCId < 577; TCId++) {
225 if (
TCID[TCId - 1] < mean_FWD * 0.1)
m_FWD++;
226 }
else if (TCId < 513) {
227 if (
TCID[TCId - 1] < mean_BAR * 0.1)
m_BAR++;
229 if (
TCID[TCId - 1] < mean_BWD * 0.1)
m_BWD++;
234 const int etot_size =
etot.size();
235 for (
int i = 0; i < etot_size; i++) {
240 RooRealVar* E_tot =
new RooRealVar(
"E_tot",
"E_tot", 0, 4000);
241 RooDataHist* data_E_tot =
new RooDataHist(
"data_E_tot",
"dataset with E_tot", *E_tot,
h_etot);
242 RooRealVar* mean_check =
new RooRealVar(
"mean_check",
"Mean for checking", 2000, 1000, 3000);
243 RooRealVar* sigma_check =
new RooRealVar(
"sigma_check",
"Sigma for checking", 100, 0, 300);
244 RooGaussian* gauss_check =
new RooGaussian(
"gauss_check",
"Gaussian for checking", *E_tot, *mean_check, *sigma_check);
245 gauss_check->fitTo(*data_E_tot);
246 double mean_ch = mean_check->getVal();
248 RooRealVar* mean_E_tot =
new RooRealVar(
"mean_E_tot",
"Mean of Gaussian", mean_ch, mean_ch - 200, mean_ch + 200);
251 RooRealVar* sigma1_E_tot =
new RooRealVar(
"sigma1_E_tot",
"Width of Gaussian", 30, 0, 60);
252 RooGaussian* gauss1_E_tot =
new RooGaussian(
"gauss1_E_tot",
"gauss(x,mean,sigma)", *E_tot, *mean_E_tot, *sigma1_E_tot);
253 RooRealVar* sigma2_E_tot =
new RooRealVar(
"sigma2_E_tot",
"Width of Gaussian", 100, 30, 250);
254 RooGaussian* gauss2_E_tot =
new RooGaussian(
"gauss2_E_tot",
"gauss(x,mean,sigma)", *E_tot, *mean_E_tot, *sigma2_E_tot);
256 RooRealVar* frac_core_E_tot =
new RooRealVar(
"frac_core",
"core fraction", 0.85, 0.0, 1.0);
257 RooAddModel* gaussm_E_tot =
new RooAddModel(
"gaussm",
"core+tail gauss", RooArgList(*gauss1_E_tot, *gauss2_E_tot),
258 RooArgList(*frac_core_E_tot));
262 gaussm_E_tot->fitTo(*data_E_tot, RooFit::Range(mean_ch - 300, mean_ch + 300));
264 RooPlot* xframe_E_tot = E_tot->frame(RooFit::Title(
"energy_tot"));
266 data_E_tot->plotOn(xframe_E_tot, RooFit::MarkerStyle(7));
268 gaussm_E_tot->plotOn(xframe_E_tot, RooFit::Name(
"Cal_energy"), RooFit::LineColor(2));
269 double v_m_E_tot = mean_E_tot->getVal();
270 double e_m_E_tot = mean_E_tot->getError();
272 double v1_s1_E_tot = sigma1_E_tot->getVal();
274 double v2_s1_E_tot = sigma2_E_tot->getVal();
277 double frac_E_tot = frac_core_E_tot->getVal();
278 double v_s1_E_tot = v1_s1_E_tot * (frac_E_tot) + v2_s1_E_tot * (1 - frac_E_tot);
282 for (
long unsigned int ii = 0; ii <
caltime.size(); ii++) {
288 RooRealVar T_c(
"T_c",
"T_c", -100, 100);
289 RooDataHist data_T_c(
"data_T_c",
"dataset with T_c", T_c,
h_caltime);
291 RooRealVar mean_T_c(
"mean_T_c",
"Mean of Gaussian", 0, -100, 100);
292 RooRealVar sigma1_T_c(
"sigma1_T_c",
"Width of Gaussian", 10, 0, 50);
293 RooGaussian gauss1_T_c(
"gauss1_T_c",
"gauss(x,mean,sigma)", T_c, mean_T_c, sigma1_T_c);
294 RooRealVar sigma2_T_c(
"sigma2_T_c",
"Width of Gaussian", 30, 10, 100);
295 RooGaussian gauss2_T_c(
"gauss2_T_c",
"gauss(x,mean,sigma)", T_c, mean_T_c, sigma2_T_c);
297 RooRealVar frac_core_T_c(
"frac_core",
"core fraction", 0.85, 0.0, 1.0);
298 RooAddModel gaussm_T_c(
"gaussm",
"core+tail gauss", RooArgList(gauss1_T_c, gauss2_T_c), RooArgList(frac_core_T_c));
300 gaussm_T_c.fitTo(data_T_c, RooFit::Range(-45, 45));
301 RooPlot* xframe_T_c = T_c.frame(RooFit::Title(
"Caltime"));
302 data_T_c.plotOn(xframe_T_c, RooFit::MarkerStyle(7));
303 gaussm_T_c.plotOn(xframe_T_c, RooFit::Name(
"Caltime"), RooFit::LineColor(2));
304 gaussm_T_c.plotOn(xframe_T_c, RooFit::Components(
"gauss1_T_c"), RooFit::LineColor(kGreen));
305 gaussm_T_c.plotOn(xframe_T_c, RooFit::Components(
"gauss2_T_c"), RooFit::LineColor(kMagenta));
307 xframe_T_c->GetXaxis()->SetTitle(
"Caltime[ns]");
308 xframe_T_c->GetYaxis()->SetTitle(
"Entries");
310 double v_m_T_c = mean_T_c.getVal();
311 double e_m_T_c = mean_T_c.getError();
313 double v1_s1_T_c = sigma1_T_c.getVal();
315 double v2_s1_T_c = sigma2_T_c.getVal();
318 double frac_T_c = frac_core_T_c.getVal();
319 double v_s1_T_c = v1_s1_T_c * (frac_T_c) + v2_s1_T_c * (1 - frac_T_c);
324 const int cluster_size =
cluster.size();
325 for (
int ii = 0; ii < cluster_size; ii++) {
329 RooRealVar x(
"x",
"esum", 0, 4000);
330 RooDataHist dh(
"dh",
"ecltrg", x,
h_clusterE);
332 RooRealVar mean(
"mean",
"mean", 2000, 1000, 2500);
333 RooRealVar sigma(
"sigma",
"sigma", 150, 0, 300);
334 RooRealVar tail(
"tail",
"tail", 0.45, 0, 1);
335 RooNovosibirsk novo(
"novo",
"", x, mean, sigma, tail);
336 novo.fitTo(dh, RooFit::Extended(0), RooFit::Range(1500, 2200));
339 double clusterE_vm, clusterE_em, clusterE_vs;
341 clusterE_vm = mean.getVal();
342 clusterE_em = mean.getError();
343 clusterE_vs = sigma.getVal();
347 TFile file(outputfile,
"RECREATE");
348 TTree* tree =
new TTree(
"tree",
"tree");
369 tree->Branch(
"m_nRun", &
m_nRun,
"m_nRun/I");
370 tree->Branch(
"m_nevent", &nevent,
"m_nevent/I");
371 tree->Branch(
"m_FWD", &FWD,
"m_FWD/D");
372 tree->Branch(
"m_BAR", &BAR,
"m_BAR/D");
373 tree->Branch(
"m_BWD", &BWD,
"m_BWD/D");
374 tree->Branch(
"m_ALL", &ALL,
"m_ALL/D");
375 tree->Branch(
"m_etot_mean", &
m_etot_mean,
"m_etot_mean/D");
376 tree->Branch(
"m_etot_error", &
m_etot_error,
"m_etot_error/D");
377 tree->Branch(
"m_etot_sigma", &
m_etot_sigma,
"m_etot_sigma/D");
378 tree->Branch(
"m_caltime_mean", &
m_caltime_mean,
"m_caltime_mean/D");
379 tree->Branch(
"m_caltime_error", &
m_caltime_error,
"m_caltime_error/D");
380 tree->Branch(
"m_caltime_sigma", &
m_caltime_sigma,
"m_caltime_sigma/D");
381 tree->Branch(
"m_clusterE_mean", &clusterE_vm,
"m_clusterE_mean/D");
382 tree->Branch(
"m_clusterE_error", &clusterE_em,
"m_clusterE_error/D");
383 tree->Branch(
"m_clusterE_sigma", &clusterE_vs,
"m_clusterE_sigma/D");
391 B2INFO(
"ECL Trigger QAM result "
393 <<
LogVar(
"Total Event", nevent)
394 <<
LogVar(
"The # of Low Hit TC in Forward Endcap", FWD)
395 <<
LogVar(
"The # of Low Hit TC in Barrel", BAR)
396 <<
LogVar(
"The # of Low Hit TC in Backward Endcap", BWD)
401 <<
LogVar(
"Cluster Energy Peak", clusterE_vm)
402 <<
LogVar(
"Cluster Energy Peak width", clusterE_vs));
405 B2DEBUG(100,
"TRGECLQAMModule ... endRun called ");