70 Gaudi::svcLocator() -> service(
"MessageSvc",
msgSvc);
71 MsgStream log(
msgSvc,
"XtMdcCalib");
72 log << MSG::INFO <<
"XtMdcCalib::initialize()" << endreq;
75 m_mdcGeomSvc = mdcGeomSvc;
76 m_mdcFunSvc = mdcFunSvc;
77 m_mdcUtilitySvc = mdcUtilitySvc;
90 m_fdXt =
new TFolder(
"fdXt",
"fdXt");
91 m_hlist ->
Add(m_fdXt);
93 m_nlayer = m_mdcGeomSvc -> getLayerSize();
96 for(lay=0; lay<m_nlayer; lay++){
97 for(iEntr=0; iEntr<m_nEntr[lay]; iEntr++){
98 for(iLR=0; iLR<m_nLR; iLR++){
100 sprintf(hname,
"Hxt%02d_E%02d_LR%01d_B%02d", lay, iEntr, iLR,
bin);
101 hist =
new TH1D(hname,
"", 300, -1.5, 1.5);
102 m_hxt.push_back( hist );
103 m_fdXt ->
Add( hist );
116 Gaudi::svcLocator() -> service(
"MessageSvc",
msgSvc);
117 MsgStream log(
msgSvc,
"XtMdcCalib");
118 log << MSG::DEBUG <<
"XtMdcCalib::fillHist()" << endreq;
123 bool esCutFg =
event->getEsCutFlag();
124 if( ! esCutFg )
return -1;
154 for(iEntr=0; iEntr<m_nEntr[lay]; iEntr++){
157 if(0 == iEntr) m_mdcFunSvc -> getXtpar(lay, 8, iLR, xtpar);
158 else if(1 == iEntr) m_mdcFunSvc -> getXtpar(lay, 9, iLR, xtpar);
159 m_tm[lay][iEntr][iLR] = xtpar[6];
170 int ntrk =
event -> getNTrk();
171 if((ntrk < m_param.
nTrkCut[0]) || (ntrk > m_param.
nTrkCut[1]))
return -1;
172 for(i=0; i<ntrk; i++){
173 rectrk =
event->getRecTrk(i);
174 nhit = rectrk -> getNHits();
177 dr = rectrk->
getDr();
178 if(fabs(dr) > m_param.
drCut)
continue;
181 dz = rectrk->
getDz();
182 if(fabs(dz) > m_param.
dzCut)
continue;
185 double p = rectrk -> getP();
186 if((fabs(p) < m_param.
pCut[0]) || (fabs(p) > m_param.
pCut[1]))
continue;
189 double phi0 = rectrk -> getPhi0();
190 double phiTrk = phi0 + CLHEP::halfpi;
191 if(phiTrk > CLHEP::twopi) phiTrk -= CLHEP::twopi;
192 if(m_param.
cosmicDwTrk && (phiTrk<CLHEP::pi))
continue;
195 fgHitLay[lay] =
false;
198 for(k=0; k<nhit; k++){
199 rechit = rectrk -> getRecHit(k);
200 lay = rechit -> getLayid();
201 fgHitLay[lay] =
true;
206 if(fgHitLay[lay]) nhitlay++;
213 for(k=0; k<nhit; k++){
214 rechit = rectrk -> getRecHit(k);
215 lay = rechit -> getLayid();
216 doca = rechit -> getDocaInc();
217 iLR = rechit -> getLR();
218 entrance = rechit -> getEntra();
219 resi = rechit -> getResiInc();
220 tdr = rechit -> getTdrift();
221 stat = rechit -> getStat();
223 if(1 != stat)
continue;
225 if( (fabs(doca) > m_docaMax[lay]) ||
226 (fabs(resi) > m_param.
resiCut[lay]) ){
231 if( ! fgHitLay[1] )
continue;
232 }
else if(42 == lay){
233 if( ! fgHitLay[41] )
continue;
235 if( (!fgHitLay[lay-1]) && (!fgHitLay[lay+1]) )
continue;
238 iEntr = m_mdcFunSvc -> getXtEntrIndex(entrance);
239 if(1 == m_nEntr[lay]){
241 }
else if(2 == m_nEntr[lay]){
242 if(entrance > 0.0) iEntr = 1;
246 bin = (int)(tdr / m_tbinw);
250 if((
bin < m_nbin) && (1 == m_hxtmap.count(
key))){
252 m_hxt[hid] ->
Fill( resi );
257 if((
bin < m_nbin) && (1 == m_hxtmap.count(
key))){
259 m_hxt[hid] ->
Fill( resi );
262 if(fabs(tdr - m_tm[lay][iEntr][iLR]) < 5){
264 if( 1 == m_hxtmap.count(
key) ){
266 m_hxt[hid] ->
Fill( resi );
270 if( 1 == m_hxtmap.count(
key) ){
272 m_hxt[hid] ->
Fill( resi );
287 Gaudi::svcLocator() -> service(
"MessageSvc",
msgSvc);
288 MsgStream log(
msgSvc,
"XtMdcCalib");
289 log << MSG::INFO <<
"XtMdcCalib::updateConst()" << endreq;
319 Double_t arglist[10];
321 TMinuit* gmxt =
new TMinuit(6);
322 gmxt -> SetPrintLevel(-1);
323 gmxt -> SetFCN(
fcnXT);
324 gmxt -> SetErrorDef(1.0);
325 gmxt -> mnparm(0,
"xtpar0", 0, 0.1, 0, 0, ierflg);
326 gmxt -> mnparm(1,
"xtpar1", 0, 0.1, 0, 0, ierflg);
327 gmxt -> mnparm(2,
"xtpar2", 0, 0.1, 0, 0, ierflg);
328 gmxt -> mnparm(3,
"xtpar3", 0, 0.1, 0, 0, ierflg);
329 gmxt -> mnparm(4,
"xtpar4", 0, 0.1, 0, 0, ierflg);
330 gmxt -> mnparm(5,
"xtpar5", 0, 0.1, 0, 0, ierflg);
332 gmxt -> mnexcm(
"SET NOW", arglist, 0, ierflg);
334 TMinuit* gmxtEd =
new TMinuit(1);
335 gmxtEd -> SetPrintLevel(-1);
337 gmxtEd -> SetErrorDef(1.0);
338 gmxtEd -> mnparm(0,
"xtpar0", 0, 0.1, 0, 0, ierflg);
340 gmxtEd -> mnexcm(
"SET NOW", arglist, 0, ierflg);
342 ofstream fxtlog(
"xtlog");
343 for(lay=0; lay<m_nlayer; lay++){
344 if(0 == m_param.
fgCalib[lay])
continue;
346 for(iEntr=0; iEntr<m_nEntr[lay]; iEntr++){
347 for(iLR=0; iLR<m_nLR; iLR++){
349 fxtlog <<
"Layer " << setw(3) << lay << setw(3) << iEntr
350 << setw(3) << iLR << endl;
352 for(ord=0; ord<m_nxtpar; ord++){
354 if(0 == iEntr) xtpar = calconst -> getXtpar(lay, 8, iLR, ord);
355 else if(1 == iEntr) xtpar = calconst -> getXtpar(lay, 9, iLR, ord);
367 entry = m_hxt[hid] -> GetEntries();
369 deltx = m_hxt[hid] -> GetMean();
370 xerr = m_hxt[hid]->GetRMS();
376 tbcen = ( (double)
bin + 0.5 ) * m_tbinw;
377 else tbcen = m_tm[lay][iEntr][iLR];
378 xcor =
xtFun(tbcen, xtini) - deltx;
380 if((tbcen <=
Tmax) || (
bin == m_nbin)){
382 XMEAS.push_back( xcor );
389 fxtlog << setw(3) <<
bin
397 if(
XMEAS.size() < 12 ){
409 for(ord=0; ord<=5; ord++){
410 arglist[0] = ord + 1;
411 arglist[1] = xtini[ord];
412 gmxt -> mnexcm(
"SET PARameter", arglist, 2, ierflg);
419 gmxt -> mnexcm(
"SET PARameter", arglist, 2, ierflg);
420 gmxt -> mnexcm(
"FIX", arglist, 1, ierflg);
425 gmxt -> mnexcm(
"MIGRAD", arglist, 2, ierflg);
426 gmxt -> mnstat(fmin, edm, errdef, nvpar, nparx, istat);
428 fxtlog <<
"Xtpar: " << endl;
429 if( (0 == ierflg) && (istat >= 2) ){
430 for(ord=0; ord<=5; ord++){
431 gmxt -> GetParameter(ord, xtpar, xterr);
435 if(1 == m_nEntr[lay]){
437 calconst -> resetXtpar(lay, i, iLR, ord, xtpar);
438 }
else if(2 == m_nEntr[lay]){
441 calconst->
resetXtpar(lay, i, iLR, ord, xtpar);
444 calconst->
resetXtpar(lay, i, iLR, ord, xtpar);
447 fxtlog << setw(15) << xtpar << setw(15) << xterr << endl;
450 for(ord=0; ord<=5; ord++){
451 fxtlog << setw(15) << xtini[ord] << setw(15) <<
"0" << endl;
454 fxtlog << setw(15) <<
Tmax << setw(15) <<
"0" << endl;
459 gmxt -> mnexcm(
"REL", arglist, 1, ierflg);
467 arglist[1] = xtini[7];
468 gmxtEd -> mnexcm(
"SET PARameter", arglist, 2, ierflg);
472 gmxtEd -> mnexcm(
"MIGRAD", arglist, 2, ierflg);
473 gmxtEd -> mnstat(fmin, edm, errdef, nvpar, nparx, istat);
475 if( (0 == ierflg) && (istat >=2) ){
476 gmxtEd -> GetParameter(0, xtpar, xterr);
477 if(xtpar < 0.0) xtpar = 0.0;
481 if(1 == m_nEntr[lay]){
483 calconst -> resetXtpar(lay, i, iLR, 7, xtpar);
484 }
else if(2 == m_nEntr[lay]){
493 fxtlog << setw(15) << xtpar << setw(15) << xterr << endl;
495 fxtlog << setw(15) << xtini[7] << setw(15) <<
"0" << endl;
498 fxtlog << setw(15) << xtini[7] << setw(15) <<
"0" << endl;
500 fxtlog <<
"Tm " << setw(15) <<
Tmax
501 <<
" Dmax " << setw(15) <<
Dmax << endl;
516 cout <<
"Xt update finished" << endl;