271 {
273 Gaudi::svcLocator() -> service(
"MessageSvc",
msgSvc);
274 MsgStream log(
msgSvc,
"XtMdcCalib");
275 log << MSG::INFO << "XtMdcCalib::updateConst()" << endreq;
276
278
279 int i;
280 int hid;
281 int lay;
282 int iLR;
283 int iEntr;
285 int ord;
287
288 Stat_t entry;
289 double xtpar;
290 double xterr;
291 double tbcen;
292 double deltx;
293 double xcor;
294 double xerr;
295 double xtini[8];
296 double xtfit[8];
297
298 Int_t ierflg;
299 Int_t istat;
300 Int_t nvpar;
301 Int_t nparx;
302 Double_t fmin;
303 Double_t edm;
304 Double_t errdef;
305 Double_t arglist[10];
306
307 TMinuit* gmxt = new TMinuit(6);
308 gmxt -> SetPrintLevel(-1);
309 gmxt -> SetFCN(
fcnXT);
310 gmxt -> SetErrorDef(1.0);
311 gmxt -> mnparm(0, "xtpar0", 0, 0.1, 0, 0, ierflg);
312 gmxt -> mnparm(1, "xtpar1", 0, 0.1, 0, 0, ierflg);
313 gmxt -> mnparm(2, "xtpar2", 0, 0.1, 0, 0, ierflg);
314 gmxt -> mnparm(3, "xtpar3", 0, 0.1, 0, 0, ierflg);
315 gmxt -> mnparm(4, "xtpar4", 0, 0.1, 0, 0, ierflg);
316 gmxt -> mnparm(5, "xtpar5", 0, 0.1, 0, 0, ierflg);
317 arglist[0] = 0;
318 gmxt -> mnexcm("SET NOW", arglist, 0, ierflg);
319
320 TMinuit* gmxtEd = new TMinuit(1);
321 gmxtEd -> SetPrintLevel(-1);
323 gmxtEd -> SetErrorDef(1.0);
324 gmxtEd -> mnparm(0, "xtpar0", 0, 0.1, 0, 0, ierflg);
325 arglist[0] = 0;
326 gmxtEd -> mnexcm("SET NOW", arglist, 0, ierflg);
327
329 for(lay=0; lay<m_nlayer; lay++){
330 if(0 == m_param.
fgCalib[lay])
continue;
331
332 for(iEntr=0; iEntr<m_nEntr[lay]; iEntr++){
333 for(iLR=0; iLR<m_nLR; iLR++){
334
335 fxtlog << "Layer " << setw(3) << lay << setw(3) << iEntr
336 << setw(3) << iLR << endl;
337
338 for(ord=0; ord<m_nxtpar; ord++){
339
340 if(0 == iEntr) xtpar = calconst -> getXtpar(lay, 8, iLR, ord);
341 else if(1 == iEntr) xtpar = calconst -> getXtpar(lay, 9, iLR, ord);
342
343 xtini[ord] = xtpar;
344 xtfit[ord] = xtpar;
345 }
346
348
352
353 entry = m_hxt[hid] -> GetEntries();
354 if(entry > 100){
355 deltx = m_hxt[hid] -> GetMean();
356 xerr = m_hxt[hid]->GetRMS();
357 } else{
358 continue;
359 }
360
362 tbcen = ( (double)
bin + 0.5 ) * m_tbinw;
363 else tbcen = m_tm[lay][iEntr][iLR];
364 xcor =
xtFun(tbcen, xtini) - deltx;
365
366 if((tbcen <=
Tmax) || (
bin == m_nbin)){
368 XMEAS.push_back( xcor );
369 ERR.push_back( xerr );
370 } else{
373 ERRED.push_back( xerr );
374 }
375 fxtlog << setw(3) <<
bin
376 << setw(15) << deltx
377 << setw(15) << xcor
378 << setw(15) << tbcen
379 << setw(15) << xerr
380 << endl;
381 }
382
383 if(
XMEAS.size() < 12 ){
387
391
392 continue;
393 }
394
395 for(ord=0; ord<=5; ord++){
396 arglist[0] = ord + 1;
397 arglist[1] = xtini[ord];
398 gmxt -> mnexcm("SET PARameter", arglist, 2, ierflg);
399 }
400
401
403 arglist[0] = 1;
404 arglist[1] = 0.0;
405 gmxt -> mnexcm("SET PARameter", arglist, 2, ierflg);
406 gmxt -> mnexcm("FIX", arglist, 1, ierflg);
407 }
408
409 arglist[0] = 1000;
410 arglist[1] = 0.1;
411 gmxt -> mnexcm("MIGRAD", arglist, 2, ierflg);
412 gmxt -> mnstat(fmin, edm, errdef, nvpar, nparx, istat);
413
414 fxtlog << "Xtpar: " << endl;
415 if( (0 == ierflg) && (istat >= 2) ){
416 for(ord=0; ord<=5; ord++){
417 gmxt -> GetParameter(ord, xtpar, xterr);
418
419 xtfit[ord] = xtpar;
420
421 if(1 == m_nEntr[lay]){
422 for(i=0; i<18; i++)
423 calconst -> resetXtpar(lay, i, iLR, ord, xtpar);
424 } else if(2 == m_nEntr[lay]){
425 if(0 == iEntr){
426 for(i=0; i<9; i++)
427 calconst->
resetXtpar(lay, i, iLR, ord, xtpar);
428 } else{
429 for(i=9; i<18; i++)
430 calconst->
resetXtpar(lay, i, iLR, ord, xtpar);
431 }
432 }
433 fxtlog << setw(15) << xtpar << setw(15) << xterr << endl;
434 }
435 } else{
436 for(ord=0; ord<=5; ord++){
437 fxtlog << setw(15) << xtini[ord] << setw(15) << "0" << endl;
438 }
439 }
440 fxtlog << setw(15) <<
Tmax << setw(15) <<
"0" << endl;
441
442
444 arglist[0] = 1;
445 gmxt -> mnexcm("REL", arglist, 1, ierflg);
446 }
447
449
451
452 arglist[0] = 1;
453 arglist[1] = xtini[7];
454 gmxtEd -> mnexcm("SET PARameter", arglist, 2, ierflg);
455
456 arglist[0] = 1000;
457 arglist[1] = 0.1;
458 gmxtEd -> mnexcm("MIGRAD", arglist, 2, ierflg);
459 gmxtEd -> mnstat(fmin, edm, errdef, nvpar, nparx, istat);
460
461 if( (0 == ierflg) && (istat >=2) ){
462 gmxtEd -> GetParameter(0, xtpar, xterr);
463 if(xtpar < 0.0) xtpar = 0.0;
465
466
467 if(1 == m_nEntr[lay]){
468 for(i=0; i<18; i++)
469 calconst -> resetXtpar(lay, i, iLR, 7, xtpar);
470 } else if(2 == m_nEntr[lay]){
471 if(0 == iEntr){
472 for(i=0; i<9; i++)
474 } else{
475 for(i=9; i<18; i++)
477 }
478 }
479 fxtlog << setw(15) << xtpar << setw(15) << xterr << endl;
480 } else {
481 fxtlog << setw(15) << xtini[7] << setw(15) << "0" << endl;
482 }
483 } else {
484 fxtlog << setw(15) << xtini[7] << setw(15) << "0" << endl;
485 }
486 fxtlog <<
"Tm " << setw(15) <<
Tmax
487 <<
" Dmax " << setw(15) <<
Dmax << endl;
488
495 }
496 }
497 }
498
499 fxtlog.close();
500 delete gmxt;
501 delete gmxtEd;
502 cout << "Xt update finished" << endl;
503 return 1;
504}
int fgCalib[MdcCalNLayer]
void resetXtpar(int lay, int entr, int lr, int order, double val)
virtual int updateConst(MdcCalibConst *calconst)=0
static double xtFun(double t, double xtpar[])
static void fcnXT(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
static void fcnXtEdge(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)