140 Gaudi::svcLocator() -> service(
"MessageSvc",
msgSvc);
141 MsgStream log(
msgSvc,
"PreT0MdcCalib");
142 log << MSG::DEBUG <<
"PreT0MdcCalib::fillHist()" << endreq;
166 IDataProviderSvc* eventSvc =
NULL;
167 Gaudi::svcLocator()->service(
"EventDataSvc", eventSvc);
170 double tes =
event->getTes();
172 bool esCutFg =
event->getEsCutFlag();
173 if( ! esCutFg )
return -1;
176 int ntrk =
event -> getNTrk();
177 for(i=0; i<ntrk; i++){
178 rectrk =
event->getRecTrk(i);
181 double dr = rectrk->
getDr();
182 if(fabs(dr) > m_param.
drCut)
continue;
185 double dz = rectrk->
getDz();
186 if(fabs(dz) > m_param.
dzCut)
continue;
189 double p = rectrk -> getP();
190 if((fabs(p) < m_param.
pCut[0]) || (fabs(p) > m_param.
pCut[1]))
continue;
193 double phi0 = rectrk -> getPhi0();
194 double phiTrk = phi0 + CLHEP::halfpi;
195 if(phiTrk > CLHEP::twopi) phiTrk -= CLHEP::twopi;
196 if(m_param.
cosmicDwTrk && (phiTrk<CLHEP::pi))
continue;
199 nhit = rectrk -> getNHits();
200 for(k=0; k<nhit; k++){
201 rechit = rectrk -> getRecHit(k);
202 lay = rechit -> getLayid();
203 cel = rechit -> getCellid();
204 lr = rechit -> getLR();
205 zhit = rechit -> getZhit();
206 tdc = rechit -> getTdc();
207 tdr = rechit -> getTdrift();
212 zprop = fabs(zhit - m_zst[lay]);
213 bin = (int)(zprop / m_zwid[lay]);
214 tp = zprop / m_vp[lay];
215 t0 = m_mdcFunSvc -> getT0(lay, cel);
216 tof = rechit -> getTof();
218 double adc = rechit -> getQhit();
222 const MdcGeoWire* pWire = m_mdcGeomSvc -> Wire(lay, cel);
223 double y1 = pWire -> Backward().y();
224 double y2 = pWire -> Forward().y();
225 double ymid = (y1+y2)*0.5;
232 trec = traw - tp - tof - tes + m_param.
timeShift - tw;
234 m_hTrec[lay][lr] ->
Fill(trec);
235 m_hTrec[lay][2] ->
Fill(trec);
236 if(ymid > 0) m_hTrecCosm[lay][0] ->
Fill(trec);
237 else m_hTrecCosm[lay][1] ->
Fill(trec);
240 m_hTrecZ[lay][lr][
bin] ->
Fill(trec);
241 m_hTrecZ[lay][2][
bin] ->
Fill(trec);
253 Gaudi::svcLocator() -> service(
"MessageSvc",
msgSvc);
254 MsgStream log(
msgSvc,
"PreT0MdcCalib");
255 log << MSG::DEBUG <<
"PreT0MdcCalib::updateConst()" << endreq;
266 double initT0 = m_param.
initT0;
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);
291 Stat_t nEntryTot = 0;
292 for(
int ibin=1; ibin<=25; ibin++){
293 Stat_t entry = m_hTrec[lay][lr]->GetBinContent(ibin);
296 double c0Ini = (double)nEntryTot / 25.0;
297 double c1Ini = (m_hTrec[lay][lr]->GetMaximum()) - c0Ini;
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);
305 m_hTrec[lay][lr] ->
Fit(funname,
"Q",
"",
308 gStyle -> SetOptFit(11);
311 chisq = ftmin[lay][lr]->GetChisquare();
312 ndf = ftmin[lay][lr]->GetNDF();
313 chindfTmin[lay][lr] = chisq / ndf;
316 fitTminFg[lay][lr] = 1;
317 t0Fit[lay][lr] = ftmin[lay][lr]->GetParameter(4);
319 t0Fit[lay][lr] += m_param.
t0Shift;
320 t0Cal[lay][lr] = t0Fit[lay][lr] - m_param.
timeShift;
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;
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);
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);
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);
371 if(0 == fitTmaxFg[lay][lr]){
372 tmax[lay][lr] = (calconst->
getXtpar(lay, 0, lr, 6)) + t0Fit[lay][2];
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;
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]
399 cout <<
"preT0.dat was written." << endl;
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;
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);
425 if(1 != m_param.
fgCalib[lay])
continue;
429 tm = tmax[lay][lr] - t0Fit[lay][2];
432 calconst -> resetXtpar(lay, iEntr, lr, 6, tm);
444 for(lr=0; lr<2; lr++){
446 calconst -> resetSdpar(lay, iEntr, lr,
bin, sdpar);
455 delete ftmin[lay][lr];
456 delete ftmax[lay][lr];