251 {
253 Gaudi::svcLocator() -> service(
"MessageSvc",
msgSvc);
254 MsgStream log(
msgSvc,
"PreT0MdcCalib");
255 log << MSG::DEBUG << "PreT0MdcCalib::updateConst()" << endreq;
256
257
258
259
260 int lay;
261 int wir;
262 int lr;
266 double initT0 = m_param.
initT0;
267
270 double chisq;
271 double ndf;
274 char funname[200];
275
276
279
281
282
285 fitTminFg[lay][lr] = 0;
286 chindfTmin[lay][lr] = -1;
287 sprintf(funname,
"ftmin%02d_%d", lay, lr);
288 ftmin[lay][lr] = new TF1(funname, funTmin, 0, 150, 6);
289
291 Stat_t nEntryTot = 0;
292 for(int ibin=1; ibin<=25; ibin++){
293 Stat_t entry = m_hTrec[lay][lr]->GetBinContent(ibin);
294 nEntryTot += entry;
295 }
296 double c0Ini = (double)nEntryTot / 25.0;
297 double c1Ini = (m_hTrec[lay][lr]->GetMaximum()) - c0Ini;
298
299 ftmin[lay][lr] -> SetParameter(0, c0Ini);
300 ftmin[lay][lr] -> SetParameter(1, c1Ini);
301 ftmin[lay][lr] -> SetParameter(2, 0);
302 ftmin[lay][lr] -> SetParameter(4, initT0);
303 ftmin[lay][lr] -> SetParameter(5, 2);
304
305 m_hTrec[lay][lr] ->
Fit(funname,
"Q",
"",
308 gStyle -> SetOptFit(11);
309
310
311 chisq = ftmin[lay][lr]->GetChisquare();
312 ndf = ftmin[lay][lr]->GetNDF();
313 chindfTmin[lay][lr] = chisq / ndf;
314
316 fitTminFg[lay][lr] = 1;
317 t0Fit[lay][lr] = ftmin[lay][lr]->GetParameter(4);
318
319 t0Fit[lay][lr] += m_param.
t0Shift;
320 t0Cal[lay][lr] = t0Fit[lay][lr] - m_param.
timeShift;
321 }
322 }
323
324 if(0 == fitTminFg[lay][lr]){
325 wir = m_mdcGeomSvc->
Wire(lay, 0)->
Id();
326 t0Cal[lay][lr] = calconst->
getT0(wir);
327 t0Fit[lay][lr] = t0Cal[lay][lr] + m_param.
timeShift;
328 }
329 }
330
331 for(int iud=0; iud<2; iud++){
332 sprintf(funname,
"ftminCosm_%02d_%d", lay, iud);
333 ftminCosm[lay][iud] = new TF1(funname, funTmin, 0, 150, 6);
334 ftminCosm[lay][iud] -> SetParameter(0, 0);
335 ftminCosm[lay][iud] -> SetParameter(4, initT0);
336 ftminCosm[lay][iud] -> SetParameter(5, 1);
337 m_hTrecCosm[lay][iud] ->
Fit(funname,
"Q",
"",
340 gStyle -> SetOptFit(11);
341 t0FitCosm[lay][iud] += m_param.
t0Shift;
342 t0FitCosm[lay][iud] = ftminCosm[lay][iud]->GetParameter(4);
343 }
344 }
345
346
350 fitTmaxFg[lay][lr] = 0;
351 chindfTmax[lay][lr] = -1;
352 sprintf(funname,
"ftmax%02d_%d", lay, lr);
353 ftmax[lay][lr] = new TF1(funname, funTmax, 250, 500, 4);
354
356 ftmax[lay][lr] -> SetParameter(2, m_param.
initTm[lay]);
357 ftmax[lay][lr] -> SetParameter(3, 10);
358 m_hTrec[lay][lr] ->
Fit(funname,
"Q+",
"",
361 gStyle -> SetOptFit(11);
362 chisq = ftmax[lay][lr]->GetChisquare();
363 ndf = ftmax[lay][lr]->GetNDF();
364 chindfTmax[lay][lr] = chisq / ndf;
366 fitTmaxFg[lay][lr] = 1;
367 tmax[lay][lr] = ftmax[lay][lr]->GetParameter(2);
368 }
369 }
370
371 if(0 == fitTmaxFg[lay][lr]){
372 tmax[lay][lr] = (calconst->
getXtpar(lay, 0, lr, 6)) + t0Fit[lay][2];
373 }
374 }
375 }
376
377
378 ofstream ft0("preT0.dat");
380 ft0 << setw(5) << lay << setw(3) << fitTminFg[lay][2]
381 << setw(15) << t0Cal[lay][2] << setw(15) << t0Fit[lay][2]
382 << setw(15) << chindfTmin[lay][2] << endl;
383 }
384 ft0 << endl;
386 ft0 << setw(5) << lay
387 << setw(3) << fitTmaxFg[lay][0] << setw(10) << tmax[lay][0]
388 << setw(10) << chindfTmax[lay][0]
389 << setw(3) << fitTmaxFg[lay][1] << setw(10) << tmax[lay][1]
390 << setw(10) << chindfTmax[lay][1]
391 << setw(3) << fitTmaxFg[lay][2] << setw(10) << tmax[lay][2]
392 << setw(10) << chindfTmax[lay][2]
393 << setw(10) << tmax[lay][0] - t0Fit[lay][2]
394 << setw(10) << tmax[lay][1] - t0Fit[lay][2]
395 << setw(10) << tmax[lay][2] - t0Fit[lay][2]
396 << endl;
397 }
398 ft0.close();
399 cout << "preT0.dat was written." << endl;
400
401
402 ofstream ft0cosm("cosmicT0.dat");
404 ft0cosm << setw(5) << lay << setw(15) << t0Fit[lay][2]
405 << setw(15) << t0FitCosm[lay][0] << setw(15) << t0FitCosm[lay][1] << endl;
406 }
407 ft0cosm.close();
408
409
410 int i;
411 int nwire = m_mdcGeomSvc -> getWireSize();
412 for(i=0; i<nwire; i++){
413 lay = m_mdcGeomSvc -> Wire(i) -> Layer();
415 calconst -> resetT0(i, t0Cal[lay][2]);
416 calconst -> resetDelT0(i, 0.0);
417 }
418 }
419
420
422 int iEntr;
423 double tm;
425 if(1 != m_param.
fgCalib[lay])
continue;
426
429 tm = tmax[lay][lr] - t0Fit[lay][2];
432 calconst -> resetXtpar(lay, iEntr, lr, 6, tm);
433 }
434 }
435 }
436 }
437 }
438
439
444 for(lr=0; lr<2; lr++){
446 calconst -> resetSdpar(lay, iEntr, lr,
bin, sdpar);
447 }
448 }
449 }
450 }
451
452
455 delete ftmin[lay][lr];
456 delete ftmax[lay][lr];
457 }
458 }
459 return 1;
460}
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)