228 {
230 Gaudi::svcLocator() -> service(
"MessageSvc",
msgSvc);
231 MsgStream log(
msgSvc,
"PreT0MdcCalib");
232 log << MSG::DEBUG << "PreT0MdcCalib::updateConst()" << endreq;
233
234
235
236
237 int lay;
238 int wir;
239 int lr;
243 double initT0 = m_param.
initT0;
244
247 double chisq;
248 double ndf;
251 char funname[200];
252
253
256
258
259
262 fitTminFg[lay][lr] = 0;
263 chindfTmin[lay][lr] = -1;
264 sprintf(funname, "ftmin%02d_%d", lay, lr);
265 ftmin[lay][lr] = new TF1(funname, funTmin, 0, 150, 6);
266
268 Stat_t nEntryTot = 0;
269 for(int ibin=1; ibin<=25; ibin++){
270 Stat_t entry = m_hTrec[lay][lr]->GetBinContent(ibin);
271 nEntryTot += entry;
272 }
273 double c0Ini = (double)nEntryTot / 25.0;
274 double c1Ini = (m_hTrec[lay][lr]->GetMaximum()) - c0Ini;
275
276 ftmin[lay][lr] -> SetParameter(0, c0Ini);
277 ftmin[lay][lr] -> SetParameter(1, c1Ini);
278 ftmin[lay][lr] -> SetParameter(2, 0);
279 ftmin[lay][lr] -> SetParameter(4, initT0);
280 ftmin[lay][lr] -> SetParameter(5, 2);
281
282 m_hTrec[lay][lr] ->
Fit(funname,
"Q",
"",
285 gStyle -> SetOptFit(11);
286
287
288 chisq = ftmin[lay][lr]->GetChisquare();
289 ndf = ftmin[lay][lr]->GetNDF();
290 chindfTmin[lay][lr] = chisq / ndf;
291
293 fitTminFg[lay][lr] = 1;
294 t0Fit[lay][lr] = ftmin[lay][lr]->GetParameter(4);
295
296 t0Fit[lay][lr] += m_param.
t0Shift;
297 t0Cal[lay][lr] = t0Fit[lay][lr] - m_param.
timeShift;
298 }
299 }
300
301 if(0 == fitTminFg[lay][lr]){
302 wir = m_mdcGeomSvc->
Wire(lay, 0)->
Id();
303 t0Cal[lay][lr] = calconst->
getT0(wir);
304 t0Fit[lay][lr] = t0Cal[lay][lr] + m_param.
timeShift;
305 }
306 }
307
308 for(int iud=0; iud<2; iud++){
309 sprintf(funname, "ftminCosm_%02d_%d", lay, iud);
310 ftminCosm[lay][iud] = new TF1(funname, funTmin, 0, 150, 6);
311 ftminCosm[lay][iud] -> SetParameter(0, 0);
312 ftminCosm[lay][iud] -> SetParameter(4, initT0);
313 ftminCosm[lay][iud] -> SetParameter(5, 1);
314 m_hTrecCosm[lay][iud] ->
Fit(funname,
"Q",
"",
317 gStyle -> SetOptFit(11);
318 t0FitCosm[lay][iud] += m_param.
t0Shift;
319 t0FitCosm[lay][iud] = ftminCosm[lay][iud]->GetParameter(4);
320 }
321 }
322
323
327 fitTmaxFg[lay][lr] = 0;
328 chindfTmax[lay][lr] = -1;
329 sprintf(funname, "ftmax%02d_%d", lay, lr);
330 ftmax[lay][lr] = new TF1(funname, funTmax, 250, 500, 4);
331
333 ftmax[lay][lr] -> SetParameter(2, m_param.
initTm[lay]);
334 ftmax[lay][lr] -> SetParameter(3, 10);
335 m_hTrec[lay][lr] ->
Fit(funname,
"Q+",
"",
338 gStyle -> SetOptFit(11);
339 chisq = ftmax[lay][lr]->GetChisquare();
340 ndf = ftmax[lay][lr]->GetNDF();
341 chindfTmax[lay][lr] = chisq / ndf;
343 fitTmaxFg[lay][lr] = 1;
344 tmax[lay][lr] = ftmax[lay][lr]->GetParameter(2);
345 }
346 }
347
348 if(0 == fitTmaxFg[lay][lr]){
349 tmax[lay][lr] = (calconst->
getXtpar(lay, 0, lr, 6)) + t0Fit[lay][2];
350 }
351 }
352 }
353
354
355 ofstream ft0("preT0.dat");
357 ft0 << setw(5) << lay << setw(3) << fitTminFg[lay][2]
358 << setw(15) << t0Cal[lay][2] << setw(15) << t0Fit[lay][2]
359 << setw(15) << chindfTmin[lay][2] << endl;
360 }
361 ft0 << endl;
363 ft0 << setw(5) << lay
364 << setw(3) << fitTmaxFg[lay][0] << setw(10) << tmax[lay][0]
365 << setw(10) << chindfTmax[lay][0]
366 << setw(3) << fitTmaxFg[lay][1] << setw(10) << tmax[lay][1]
367 << setw(10) << chindfTmax[lay][1]
368 << setw(3) << fitTmaxFg[lay][2] << setw(10) << tmax[lay][2]
369 << setw(10) << chindfTmax[lay][2]
370 << setw(10) << tmax[lay][0] - t0Fit[lay][2]
371 << setw(10) << tmax[lay][1] - t0Fit[lay][2]
372 << setw(10) << tmax[lay][2] - t0Fit[lay][2]
373 << endl;
374 }
375 ft0.close();
376 cout << "preT0.dat was written." << endl;
377
378
379 ofstream ft0cosm("cosmicT0.dat");
381 ft0cosm << setw(5) << lay << setw(15) << t0Fit[lay][2]
382 << setw(15) << t0FitCosm[lay][0] << setw(15) << t0FitCosm[lay][1] << endl;
383 }
384 ft0cosm.close();
385
386
387 int i;
388 int nwire = m_mdcGeomSvc -> getWireSize();
389 for(i=0; i<nwire; i++){
390 lay = m_mdcGeomSvc -> Wire(i) -> Layer();
392 calconst -> resetT0(i, t0Cal[lay][2]);
393 calconst -> resetDelT0(i, 0.0);
394 }
395 }
396
397
399 int iEntr;
400 double tm;
402 if(1 != m_param.
fgCalib[lay])
continue;
403
406 tm = tmax[lay][lr] - t0Fit[lay][2];
409 calconst -> resetXtpar(lay, iEntr, lr, 6, tm);
410 }
411 }
412 }
413 }
414 }
415
416
421 for(lr=0; lr<2; lr++){
423 calconst -> resetSdpar(lay, iEntr, lr,
bin, sdpar);
424 }
425 }
426 }
427 }
428
429
432 delete ftmin[lay][lr];
433 delete ftmax[lay][lr];
434 }
435 }
436 return 1;
437}
double initTm[MdcCalNLayer]
int fgCalib[MdcCalNLayer]
double tminFitRange[MdcCalNLayer][2]
double tmaxFitRange[MdcCalNLayer][2]
double getT0(int wireid) const
double getXtpar(int lay, int entr, int lr, int order)