301 {
303 Gaudi::svcLocator()->service(
"MessageSvc",
msgSvc );
304 MsgStream log(
msgSvc,
"XtMdcCalib" );
305 log << MSG::INFO << "XtMdcCalib::updateConst()" << endmsg;
306
308
309 int i;
310 int hid;
311 int lay;
312 int iLR;
313 int iEntr;
315 int ord;
317
318 Stat_t entry;
319 double xtpar;
320 double xterr;
321 double tbcen;
322 double deltx;
323 double xcor;
324 double xerr;
325 double xtini[8];
326 double xtfit[8];
327
328 Int_t ierflg;
329 Int_t istat;
330 Int_t nvpar;
331 Int_t nparx;
332 Double_t fmin;
333 Double_t edm;
334 Double_t errdef;
335 Double_t arglist[10];
336
337 TMinuit* gmxt = new TMinuit( 6 );
338 gmxt->SetPrintLevel( -1 );
339 gmxt->SetFCN(
fcnXT );
340 gmxt->SetErrorDef( 1.0 );
341 gmxt->mnparm( 0, "xtpar0", 0, 0.1, 0, 0, ierflg );
342 gmxt->mnparm( 1, "xtpar1", 0, 0.1, 0, 0, ierflg );
343 gmxt->mnparm( 2, "xtpar2", 0, 0.1, 0, 0, ierflg );
344 gmxt->mnparm( 3, "xtpar3", 0, 0.1, 0, 0, ierflg );
345 gmxt->mnparm( 4, "xtpar4", 0, 0.1, 0, 0, ierflg );
346 gmxt->mnparm( 5, "xtpar5", 0, 0.1, 0, 0, ierflg );
347 arglist[0] = 0;
348 gmxt->mnexcm( "SET NOW", arglist, 0, ierflg );
349
350 TMinuit* gmxtEd = new TMinuit( 1 );
351 gmxtEd->SetPrintLevel( -1 );
353 gmxtEd->SetErrorDef( 1.0 );
354 gmxtEd->mnparm( 0, "xtpar0", 0, 0.1, 0, 0, ierflg );
355 arglist[0] = 0;
356 gmxtEd->mnexcm( "SET NOW", arglist, 0, ierflg );
357
359 for ( lay = 0; lay < m_nlayer; lay++ )
360 {
361 if ( 0 == m_param.fgCalib[lay] ) continue;
362
363 for ( iEntr = 0; iEntr < m_nEntr[lay]; iEntr++ )
364 {
365 for ( iLR = 0; iLR < m_nLR; iLR++ )
366 {
367
368 fxtlog << "Layer " << setw( 3 ) << lay << setw( 3 ) << iEntr << setw( 3 ) << iLR
369 << endl;
370
371 for ( ord = 0; ord < m_nxtpar; ord++ )
372 {
373
374 if ( 0 == iEntr ) xtpar = calconst->
getXtpar( lay, 8, iLR, ord );
375 else if ( 1 == iEntr ) xtpar = calconst->
getXtpar( lay, 9, iLR, ord );
376
377 xtini[ord] = xtpar;
378 xtfit[ord] = xtpar;
379 }
380
382
384 {
387
388 entry = m_hxt[hid]->GetEntries();
389 if ( entry > 100 )
390 {
391 deltx = m_hxt[hid]->GetMean();
392 xerr = m_hxt[hid]->GetRMS();
393 }
394 else { continue; }
395
396 if (
bin < m_nbin ) tbcen = ( (double)
bin + 0.5 ) * m_tbinw;
397 else tbcen = m_tm[lay][iEntr][iLR];
398 xcor =
xtFun( tbcen, xtini ) - deltx;
399
400 if ( ( tbcen <=
Tmax ) || (
bin == m_nbin ) )
401 {
403 XMEAS.push_back( xcor );
404 ERR.push_back( xerr );
405 }
406 else
407 {
410 ERRED.push_back( xerr );
411 }
412 fxtlog << setw( 3 ) <<
bin << setw( 15 ) << deltx << setw( 15 ) << xcor << setw( 15 )
413 << tbcen << setw( 15 ) << xerr << endl;
414 }
415
416 if (
XMEAS.size() < 12 )
417 {
421
425
426 continue;
427 }
428
429 for ( ord = 0; ord <= 5; ord++ )
430 {
431 arglist[0] = ord + 1;
432 arglist[1] = xtini[ord];
433 gmxt->mnexcm( "SET PARameter", arglist, 2, ierflg );
434 }
435
436
437 if ( 1 == m_param.fixXtC0 )
438 {
439 arglist[0] = 1;
440 arglist[1] = 0.0;
441 gmxt->mnexcm( "SET PARameter", arglist, 2, ierflg );
442 gmxt->mnexcm( "FIX", arglist, 1, ierflg );
443 }
444
445 arglist[0] = 1000;
446 arglist[1] = 0.1;
447 gmxt->mnexcm( "MIGRAD", arglist, 2, ierflg );
448 gmxt->mnstat( fmin, edm, errdef, nvpar, nparx, istat );
449
450 fxtlog << "Xtpar: " << endl;
451 if ( ( 0 == ierflg ) && ( istat >= 2 ) )
452 {
453 for ( ord = 0; ord <= 5; ord++ )
454 {
455 gmxt->GetParameter( ord, xtpar, xterr );
456
457 xtfit[ord] = xtpar;
458
459 if ( 1 == m_nEntr[lay] )
460 {
461 for ( i = 0; i < 18; i++ ) calconst->
resetXtpar( lay, i, iLR, ord, xtpar );
462 }
463 else if ( 2 == m_nEntr[lay] )
464 {
465 if ( 0 == iEntr )
466 {
467 for ( i = 0; i < 9; i++ )
468 calconst->
resetXtpar( lay, i, iLR, ord, xtpar );
469 }
470 else
471 {
472 for ( i = 9; i < 18; i++ )
473 calconst->
resetXtpar( lay, i, iLR, ord, xtpar );
474 }
475 }
476 fxtlog << setw( 15 ) << xtpar << setw( 15 ) << xterr << endl;
477 }
478 }
479 else
480 {
481 for ( ord = 0; ord <= 5; ord++ )
482 { fxtlog << setw( 15 ) << xtini[ord] << setw( 15 ) << "0" << endl; }
483 }
484 fxtlog << setw( 15 ) <<
Tmax << setw( 15 ) <<
"0" << endl;
485
486
487 if ( 1 == m_param.fixXtC0 )
488 {
489 arglist[0] = 1;
490 gmxt->mnexcm( "REL", arglist, 1, ierflg );
491 }
492
494
496 {
497
498 arglist[0] = 1;
499 arglist[1] = xtini[7];
500 gmxtEd->mnexcm( "SET PARameter", arglist, 2, ierflg );
501
502 arglist[0] = 1000;
503 arglist[1] = 0.1;
504 gmxtEd->mnexcm( "MIGRAD", arglist, 2, ierflg );
505 gmxtEd->mnstat( fmin, edm, errdef, nvpar, nparx, istat );
506
507 if ( ( 0 == ierflg ) && ( istat >= 2 ) )
508 {
509 gmxtEd->GetParameter( 0, xtpar, xterr );
510 if ( xtpar < 0.0 ) xtpar = 0.0;
511 if ( m_param.fixXtEdge ) xtpar = xtini[7];
512
513
514 if ( 1 == m_nEntr[lay] )
515 {
516 for ( i = 0; i < 18; i++ ) calconst->
resetXtpar( lay, i, iLR, 7, xtpar );
517 }
518 else if ( 2 == m_nEntr[lay] )
519 {
520 if ( 0 == iEntr )
521 {
522 for ( i = 0; i < 9; i++ ) calconst->
resetXtpar( lay, i, iLR, 7, xtpar );
523 }
524 else
525 {
526 for ( i = 9; i < 18; i++ ) calconst->
resetXtpar( lay, i, iLR, 7, xtpar );
527 }
528 }
529 fxtlog << setw( 15 ) << xtpar << setw( 15 ) << xterr << endl;
530 }
531 else { fxtlog << setw( 15 ) << xtini[7] << setw( 15 ) << "0" << endl; }
532 }
533 else { fxtlog << setw( 15 ) << xtini[7] << setw( 15 ) << "0" << endl; }
534 fxtlog <<
"Tm " << setw( 15 ) <<
Tmax <<
" Dmax " << setw( 15 ) <<
Dmax << endl;
535
542 }
543 }
544 }
545
546 fxtlog.close();
547 delete gmxt;
548 delete gmxtEd;
549 cout << "Xt update finished" << endl;
550 return 1;
551}
void resetXtpar(int lay, int entr, int lr, int order, double val)
double getXtpar(int lay, int entr, int lr, int order)
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)