3#include "GaudiKernel/MsgStream.h"
4#include "GaudiKernel/IMessageSvc.h"
5#include "GaudiKernel/StatusCode.h"
6#include "GaudiKernel/ISvcLocator.h"
7#include "GaudiKernel/Bootstrap.h"
9#include "GaudiKernel/SmartDataPtr.h"
10#include "GaudiKernel/IDataProviderSvc.h"
11#include "GaudiKernel/PropertyMgr.h"
44 delete m_hTrec[lay][lr];
46 delete m_hTrecZ[lay][lr][
bin];
62 Gaudi::svcLocator() -> service(
"MessageSvc",
msgSvc);
63 MsgStream log(
msgSvc,
"PreT0MdcCalib");
64 log << MSG::INFO <<
"PreT0MdcCalib::initialize()" << endreq;
67 m_mdcGeomSvc = mdcGeomSvc;
68 m_mdcFunSvc = mdcFunSvc;
69 m_mdcUtilitySvc = mdcUtilitySvc;
78 m_fdTrec =
new TFolder(
"Trec",
"Trec");
79 m_hlist->Add(m_fdTrec);
81 m_fdTrecZ =
new TFolder(
"TrecZ",
"TrecZ");
82 m_hlist->Add(m_fdTrecZ);
86 if(0 == lr)
sprintf(hname,
"Trec%02d_L", lay);
87 else if(1 == lr)
sprintf(hname,
"Trec%02d_R", lay);
88 else sprintf(hname,
"Trec%02d_C", lay);
90 if(lay < 8) m_hTrec[lay][lr] =
new TH1F(hname,
"", 100, 0, 400);
91 else m_hTrec[lay][lr] =
new TH1F(hname,
"", 125, 0, 500);
92 m_fdTrec ->
Add(m_hTrec[lay][lr]);
97 for(
int iud=0; iud<2; iud++){
98 if(0 == iud)
sprintf(hname,
"TrecCosm%02d_Up", lay);
99 else sprintf(hname,
"TrecCosm%02d_Dw", lay);
100 if(lay < 8) m_hTrecCosm[lay][iud] =
new TH1F(hname,
"", 100, 0, 400);
101 else m_hTrecCosm[lay][iud] =
new TH1F(hname,
"", 125, 0, 500);
102 m_fdTrec ->
Add(m_hTrecCosm[lay][iud]);
109 if(0 == lr)
sprintf(hname,
"Trec%02d_z%02d_L", lay,
bin);
110 else if(1 == lr)
sprintf(hname,
"Trec%02d_z%02d_R", lay,
bin);
111 else sprintf(hname,
"Trec%02d_z%02d_C", lay,
bin);
113 if(lay < 8) m_hTrecZ[lay][lr][
bin] =
new TH1F(hname,
"", 100, 0, 400);
114 else m_hTrecZ[lay][lr][
bin] =
new TH1F(hname,
"", 125, 0, 500);
115 m_fdTrecZ ->
Add(m_hTrecZ[lay][lr][
bin]);
125 m_zwid[lay] = (zeast - zwest) / (
double)m_nzbin;
127 if(lay < 8) m_vp[lay] = 220.0;
128 else m_vp[lay] = 240.0;
130 if( 0 == (lay % 2) ){
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];
461Double_t PreT0MdcCalib::funTmin(Double_t* x, Double_t* par){
463 fitval = par[0] + par[1]*
exp( -par[2]*(x[0]-par[3]) ) /
464 ( 1 +
exp( -(x[0]-par[4])/par[5] ));
468Double_t PreT0MdcCalib::funTmax(Double_t* x, Double_t* par){
470 fitval = par[0] + par[1] / (1 +
exp((x[0]-par[2])/par[3]));
sprintf(cut,"kal_costheta0_em>-0.93&&kal_costheta0_em<0.93&&kal_pxy0_em>=0.05+%d*0.1&&kal_pxy0_em<0.15+%d*0.1&&NGch>=2", j, j)
EvtComplex exp(const EvtComplex &c)
*******INTEGER m_nBinMax INTEGER m_NdiMax !No of bins in histogram for cell exploration division $ !Last vertex $ !Last active cell $ !Last cell in buffer $ !No of sampling when dividing cell $ !No of function total $ !Flag for random ceel for $ !Flag for type of for WtMax $ !Flag which decides whether vertices are included in the sampling $ entire domain is hyp !Maximum effective eevents per bin
const double MdcCalTdcCnv
virtual double getTimeWalk(int layid, double Q) const =0
virtual const MdcGeoWire *const Wire(unsigned id)=0
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)
HepPoint3D Forward(void) const
HepPoint3D Backward(void) const
int updateConst(MdcCalibConst *calconst)
int fillHist(MdcCalEvent *event)
void initialize(TObjArray *hlist, IMdcGeomSvc *mdcGeomSvc, IMdcCalibFunSvc *mdcFunSvc, IMdcUtilitySvc *mdcUtilitySvc)