00001 #include "MdcCalibAlg/XtMdcCalib.h"
00002
00003 #include "GaudiKernel/MsgStream.h"
00004 #include "GaudiKernel/IMessageSvc.h"
00005 #include "GaudiKernel/StatusCode.h"
00006 #include "GaudiKernel/ISvcLocator.h"
00007 #include "GaudiKernel/Bootstrap.h"
00008
00009 #include <iostream>
00010 #include <fstream>
00011 #include <iomanip>
00012 #include <string>
00013 #include <cstring>
00014 #include <cmath>
00015
00016 #include "TF1.h"
00017 #include "TMinuit.h"
00018 #include "TRandom.h"
00019
00020 using namespace std;
00021
00022 typedef map<int, int>::value_type valType2;
00023
00024 vector<double> XtMdcCalib::XMEAS;
00025 vector<double> XtMdcCalib::TBINCEN;
00026 vector<double> XtMdcCalib::ERR;
00027 double XtMdcCalib::Tmax;
00028 double XtMdcCalib::Dmax;
00029 vector<double> XtMdcCalib::XMEASED;
00030 vector<double> XtMdcCalib::TBINCENED;
00031 vector<double> XtMdcCalib::ERRED;
00032
00033
00034 XtMdcCalib::XtMdcCalib(){
00035 m_nbin = 50;
00036 m_nLR = 3;
00037 m_nxtpar = 8;
00038
00039 for(int lay=0; lay<43; lay++){
00040 if(lay < 15) m_nEntr[lay] = 1;
00041 else m_nEntr[lay] = 2;
00042 }
00043
00044 m_tbinw = 10.0;
00045 m_fgIni = false;
00046 }
00047
00048 XtMdcCalib::~XtMdcCalib(){
00049 }
00050
00051 void XtMdcCalib::clear(){
00052 cout << "~XtMdcCalib" << endl;
00053
00054 unsigned int i;
00055 for(i=0; i<m_hxt.size(); i++){
00056 delete m_hxt[i];
00057 if(0 == (i % 1000)) cout << "delete hxt[" << i << "]" << endl;
00058 }
00059 delete m_fdXt;
00060 m_hxt.clear();
00061 m_hxtmap.clear();
00062
00063 MdcCalib::clear();
00064 }
00065
00066
00067 void XtMdcCalib::initialize(TObjArray* hlist, IMdcGeomSvc* mdcGeomSvc,
00068 IMdcCalibFunSvc* mdcFunSvc, IMdcUtilitySvc* mdcUtilitySvc) {
00069 IMessageSvc* msgSvc;
00070 Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
00071 MsgStream log(msgSvc, "XtMdcCalib");
00072 log << MSG::INFO << "XtMdcCalib::initialize()" << endreq;
00073
00074 m_hlist = hlist;
00075 m_mdcGeomSvc = mdcGeomSvc;
00076 m_mdcFunSvc = mdcFunSvc;
00077 m_mdcUtilitySvc = mdcUtilitySvc;
00078
00079 MdcCalib::initialize(m_hlist, m_mdcGeomSvc, m_mdcFunSvc, m_mdcUtilitySvc);
00080
00081 int lay;
00082 int iLR;
00083 int iEntr;
00084 int bin;
00085 int hid;
00086 int key;
00087 char hname[200];
00088 TH1D* hist;
00089
00090 m_fdXt = new TFolder("fdXt", "fdXt");
00091 m_hlist -> Add(m_fdXt);
00092
00093 m_nlayer = m_mdcGeomSvc -> getLayerSize();
00094
00095 hid = 0;
00096 for(lay=0; lay<m_nlayer; lay++){
00097 for(iEntr=0; iEntr<m_nEntr[lay]; iEntr++){
00098 for(iLR=0; iLR<m_nLR; iLR++){
00099 for(bin=0; bin<=m_nbin; bin++){
00100 sprintf(hname, "Hxt%02d_E%02d_LR%01d_B%02d", lay, iEntr, iLR, bin);
00101 hist = new TH1D(hname, "", 300, -1.5, 1.5);
00102 m_hxt.push_back( hist );
00103 m_fdXt -> Add( hist );
00104
00105 key = getHxtKey(lay, iEntr, iLR, bin);
00106 m_hxtmap.insert( valType2( key, hid ) );
00107 hid++;
00108 }
00109 }
00110 }
00111 }
00112 }
00113
00114 int XtMdcCalib::fillHist(MdcCalEvent* event){
00115 IMessageSvc* msgSvc;
00116 Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
00117 MsgStream log(msgSvc, "XtMdcCalib");
00118 log << MSG::DEBUG << "XtMdcCalib::fillHist()" << endreq;
00119
00120 MdcCalib::fillHist(event);
00121
00122
00123 bool esCutFg = event->getEsCutFlag();
00124 if( ! esCutFg ) return -1;
00125
00126 int i;
00127 int k;
00128 int lay;
00129 int iLR;
00130 int iEntr;
00131 int bin;
00132 int key;
00133 int hid;
00134
00135 double dr;
00136 double dz;
00137 double tdr;
00138 double doca;
00139 double resi;
00140 double entrance;
00141 int stat;
00142
00143 double xtpar[MdcCalXtNPars];
00144 int nhitlay;
00145 bool fgHitLay[MdcCalNLayer];
00146
00147 if(! m_fgIni){
00148 for(lay=0; lay<MdcCalNLayer; lay++){
00149 if(lay < 8) m_docaMax[lay] = m_param.maxDocaInner;
00150 else m_docaMax[lay] = m_param.maxDocaOuter;
00151 }
00152
00153 for(lay=0; lay<MdcCalNLayer; lay++){
00154 for(iEntr=0; iEntr<m_nEntr[lay]; iEntr++){
00155 for(iLR=0; iLR<MdcCalLR; iLR++){
00156
00157 if(0 == iEntr) m_mdcFunSvc -> getXtpar(lay, 8, iLR, xtpar);
00158 else if(1 == iEntr) m_mdcFunSvc -> getXtpar(lay, 9, iLR, xtpar);
00159 m_tm[lay][iEntr][iLR] = xtpar[6];
00160 }
00161 }
00162 }
00163 m_fgIni = true;
00164 }
00165
00166 MdcCalRecTrk* rectrk;
00167 MdcCalRecHit* rechit;
00168
00169 int nhit;
00170 int ntrk = event -> getNTrk();
00171 if((ntrk < m_param.nTrkCut[0]) || (ntrk > m_param.nTrkCut[1])) return -1;
00172 for(i=0; i<ntrk; i++){
00173 rectrk = event->getRecTrk(i);
00174 nhit = rectrk -> getNHits();
00175
00176
00177 dr = rectrk->getDr();
00178 if(fabs(dr) > m_param.drCut) continue;
00179
00180
00181 dz = rectrk->getDz();
00182 if(fabs(dz) > m_param.dzCut) continue;
00183
00184 for(lay=0; lay<MdcCalNLayer; lay++){
00185 fgHitLay[lay] = false;
00186 }
00187
00188 for(k=0; k<nhit; k++){
00189 rechit = rectrk -> getRecHit(k);
00190 lay = rechit -> getLayid();
00191 fgHitLay[lay] = true;
00192 }
00193
00194 nhitlay = 0;
00195 for(lay=0; lay<MdcCalNLayer; lay++){
00196 if(fgHitLay[lay]) nhitlay++;
00197 }
00198 if(nhitlay < m_param.nHitLayCut) continue;
00199
00200
00201
00202
00203 for(k=0; k<nhit; k++){
00204 rechit = rectrk -> getRecHit(k);
00205 lay = rechit -> getLayid();
00206 doca = rechit -> getDocaInc();
00207 iLR = rechit -> getLR();
00208 entrance = rechit -> getEntra();
00209 resi = rechit -> getResiInc();
00210 tdr = rechit -> getTdrift();
00211 stat = rechit -> getStat();
00212
00213 if(1 != stat) continue;
00214
00215 if( (fabs(doca) > m_docaMax[lay]) ||
00216 (fabs(resi) > m_param.resiCut[lay]) ){
00217 continue;
00218 }
00219
00220 if(0 == lay){
00221 if( ! fgHitLay[1] ) continue;
00222 } else if(42 == lay){
00223 if( ! fgHitLay[41] ) continue;
00224 } else{
00225 if( (!fgHitLay[lay-1]) && (!fgHitLay[lay+1]) ) continue;
00226 }
00227
00228 iEntr = m_mdcFunSvc -> getXtEntrIndex(entrance);
00229 if(1 == m_nEntr[lay]){
00230 iEntr = 0;
00231 } else if(2 == m_nEntr[lay]){
00232 if(entrance > 0.0) iEntr = 1;
00233 else iEntr = 0;
00234 }
00235
00236 bin = (int)(tdr / m_tbinw);
00237
00238
00239 key = getHxtKey(lay, iEntr, iLR, bin);
00240 if((bin < m_nbin) && (1 == m_hxtmap.count(key))){
00241 hid = m_hxtmap[key];
00242 m_hxt[hid] -> Fill( resi );
00243 }
00244
00245
00246 key = getHxtKey(lay, iEntr, 2, bin);
00247 if((bin < m_nbin) && (1 == m_hxtmap.count(key))){
00248 hid = m_hxtmap[key];
00249 m_hxt[hid] -> Fill( resi );
00250 }
00251
00252 if(fabs(tdr - m_tm[lay][iEntr][iLR]) < 5){
00253 key = getHxtKey(lay, iEntr, iLR, m_nbin);
00254 if( 1 == m_hxtmap.count(key) ){
00255 hid = m_hxtmap[key];
00256 m_hxt[hid] -> Fill( resi );
00257 }
00258
00259 key = getHxtKey(lay, iEntr, 2, m_nbin);
00260 if( 1 == m_hxtmap.count(key) ){
00261 hid = m_hxtmap[key];
00262 m_hxt[hid] -> Fill( resi );
00263 }
00264 }
00265
00266 }
00267 }
00268 return 1;
00269 }
00270
00271 int XtMdcCalib::updateConst(MdcCalibConst* calconst){
00272 IMessageSvc* msgSvc;
00273 Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
00274 MsgStream log(msgSvc, "XtMdcCalib");
00275 log << MSG::INFO << "XtMdcCalib::updateConst()" << endreq;
00276
00277 MdcCalib::updateConst(calconst);
00278
00279 int i;
00280 int hid;
00281 int lay;
00282 int iLR;
00283 int iEntr;
00284 int bin;
00285 int ord;
00286 int key;
00287
00288 Stat_t entry;
00289 double xtpar;
00290 double xterr;
00291 double tbcen;
00292 double deltx;
00293 double xcor;
00294 double xerr;
00295 double xtini[8];
00296 double xtfit[8];
00297
00298 Int_t ierflg;
00299 Int_t istat;
00300 Int_t nvpar;
00301 Int_t nparx;
00302 Double_t fmin;
00303 Double_t edm;
00304 Double_t errdef;
00305 Double_t arglist[10];
00306
00307 TMinuit* gmxt = new TMinuit(6);
00308 gmxt -> SetPrintLevel(-1);
00309 gmxt -> SetFCN(fcnXT);
00310 gmxt -> SetErrorDef(1.0);
00311 gmxt -> mnparm(0, "xtpar0", 0, 0.1, 0, 0, ierflg);
00312 gmxt -> mnparm(1, "xtpar1", 0, 0.1, 0, 0, ierflg);
00313 gmxt -> mnparm(2, "xtpar2", 0, 0.1, 0, 0, ierflg);
00314 gmxt -> mnparm(3, "xtpar3", 0, 0.1, 0, 0, ierflg);
00315 gmxt -> mnparm(4, "xtpar4", 0, 0.1, 0, 0, ierflg);
00316 gmxt -> mnparm(5, "xtpar5", 0, 0.1, 0, 0, ierflg);
00317 arglist[0] = 0;
00318 gmxt -> mnexcm("SET NOW", arglist, 0, ierflg);
00319
00320 TMinuit* gmxtEd = new TMinuit(1);
00321 gmxtEd -> SetPrintLevel(-1);
00322 gmxtEd -> SetFCN(fcnXtEdge);
00323 gmxtEd -> SetErrorDef(1.0);
00324 gmxtEd -> mnparm(0, "xtpar0", 0, 0.1, 0, 0, ierflg);
00325 arglist[0] = 0;
00326 gmxtEd -> mnexcm("SET NOW", arglist, 0, ierflg);
00327
00328 ofstream fxtlog("xtlog");
00329 for(lay=0; lay<m_nlayer; lay++){
00330 if(0 == m_param.fgCalib[lay]) continue;
00331
00332 for(iEntr=0; iEntr<m_nEntr[lay]; iEntr++){
00333 for(iLR=0; iLR<m_nLR; iLR++){
00334
00335 fxtlog << "Layer " << setw(3) << lay << setw(3) << iEntr
00336 << setw(3) << iLR << endl;
00337
00338 for(ord=0; ord<m_nxtpar; ord++){
00339
00340 if(0 == iEntr) xtpar = calconst -> getXtpar(lay, 8, iLR, ord);
00341 else if(1 == iEntr) xtpar = calconst -> getXtpar(lay, 9, iLR, ord);
00342
00343 xtini[ord] = xtpar;
00344 xtfit[ord] = xtpar;
00345 }
00346
00347 Tmax = xtini[6];
00348
00349 for(bin=0; bin<=m_nbin; bin++){
00350 key = getHxtKey(lay, iEntr, iLR, bin);
00351 hid = m_hxtmap[key];
00352
00353 entry = m_hxt[hid] -> GetEntries();
00354 if(entry > 100){
00355 deltx = m_hxt[hid] -> GetMean();
00356 xerr = m_hxt[hid]->GetRMS();
00357 } else{
00358 continue;
00359 }
00360
00361 if(bin < m_nbin)
00362 tbcen = ( (double)bin + 0.5 ) * m_tbinw;
00363 else tbcen = m_tm[lay][iEntr][iLR];
00364 xcor = xtFun(tbcen, xtini) - deltx;
00365
00366 if((tbcen <= Tmax) || (bin == m_nbin)){
00367 TBINCEN.push_back( tbcen );
00368 XMEAS.push_back( xcor );
00369 ERR.push_back( xerr );
00370 } else{
00371 TBINCENED.push_back( tbcen );
00372 XMEASED.push_back( xcor );
00373 ERRED.push_back( xerr );
00374 }
00375 fxtlog << setw(3) << bin
00376 << setw(15) << deltx
00377 << setw(15) << xcor
00378 << setw(15) << tbcen
00379 << setw(15) << xerr
00380 << endl;
00381 }
00382
00383 if( XMEAS.size() < 12 ){
00384 TBINCEN.clear();
00385 XMEAS.clear();
00386 ERR.clear();
00387
00388 TBINCENED.clear();
00389 XMEASED.clear();
00390 ERRED.clear();
00391
00392 continue;
00393 }
00394
00395 for(ord=0; ord<=5; ord++){
00396 arglist[0] = ord + 1;
00397 arglist[1] = xtini[ord];
00398 gmxt -> mnexcm("SET PARameter", arglist, 2, ierflg);
00399 }
00400
00401
00402 if(1 == m_param.fixXtC0){
00403 arglist[0] = 1;
00404 arglist[1] = 0.0;
00405 gmxt -> mnexcm("SET PARameter", arglist, 2, ierflg);
00406 gmxt -> mnexcm("FIX", arglist, 1, ierflg);
00407 }
00408
00409 arglist[0] = 1000;
00410 arglist[1] = 0.1;
00411 gmxt -> mnexcm("MIGRAD", arglist, 2, ierflg);
00412 gmxt -> mnstat(fmin, edm, errdef, nvpar, nparx, istat);
00413
00414 fxtlog << "Xtpar: " << endl;
00415 if( (0 == ierflg) && (istat >= 2) ){
00416 for(ord=0; ord<=5; ord++){
00417 gmxt -> GetParameter(ord, xtpar, xterr);
00418
00419 xtfit[ord] = xtpar;
00420
00421 if(1 == m_nEntr[lay]){
00422 for(i=0; i<18; i++)
00423 calconst -> resetXtpar(lay, i, iLR, ord, xtpar);
00424 } else if(2 == m_nEntr[lay]){
00425 if(0 == iEntr){
00426 for(i=0; i<9; i++)
00427 calconst->resetXtpar(lay, i, iLR, ord, xtpar);
00428 } else{
00429 for(i=9; i<18; i++)
00430 calconst->resetXtpar(lay, i, iLR, ord, xtpar);
00431 }
00432 }
00433 fxtlog << setw(15) << xtpar << setw(15) << xterr << endl;
00434 }
00435 } else{
00436 for(ord=0; ord<=5; ord++){
00437 fxtlog << setw(15) << xtini[ord] << setw(15) << "0" << endl;
00438 }
00439 }
00440 fxtlog << setw(15) << Tmax << setw(15) << "0" << endl;
00441
00442
00443 if(1 == m_param.fixXtC0){
00444 arglist[0] = 1;
00445 gmxt -> mnexcm("REL", arglist, 1, ierflg);
00446 }
00447
00448 Dmax = xtFun(Tmax, xtfit);
00449
00450 if( XMEASED.size() >= 3 ){
00451
00452 arglist[0] = 1;
00453 arglist[1] = xtini[7];
00454 gmxtEd -> mnexcm("SET PARameter", arglist, 2, ierflg);
00455
00456 arglist[0] = 1000;
00457 arglist[1] = 0.1;
00458 gmxtEd -> mnexcm("MIGRAD", arglist, 2, ierflg);
00459 gmxtEd -> mnstat(fmin, edm, errdef, nvpar, nparx, istat);
00460
00461 if( (0 == ierflg) && (istat >=2) ){
00462 gmxtEd -> GetParameter(0, xtpar, xterr);
00463 if(xtpar < 0.0) xtpar = 0.0;
00464 if(m_param.fixXtEdge) xtpar = xtini[7];
00465
00466
00467 if(1 == m_nEntr[lay]){
00468 for(i=0; i<18; i++)
00469 calconst -> resetXtpar(lay, i, iLR, 7, xtpar);
00470 } else if(2 == m_nEntr[lay]){
00471 if(0 == iEntr){
00472 for(i=0; i<9; i++)
00473 calconst->resetXtpar(lay, i, iLR, 7, xtpar);
00474 } else{
00475 for(i=9; i<18; i++)
00476 calconst->resetXtpar(lay, i, iLR, 7, xtpar);
00477 }
00478 }
00479 fxtlog << setw(15) << xtpar << setw(15) << xterr << endl;
00480 } else {
00481 fxtlog << setw(15) << xtini[7] << setw(15) << "0" << endl;
00482 }
00483 } else {
00484 fxtlog << setw(15) << xtini[7] << setw(15) << "0" << endl;
00485 }
00486 fxtlog << "Tm " << setw(15) << Tmax
00487 << " Dmax " << setw(15) << Dmax << endl;
00488
00489 TBINCEN.clear();
00490 XMEAS.clear();
00491 ERR.clear();
00492 TBINCENED.clear();
00493 XMEASED.clear();
00494 ERRED.clear();
00495 }
00496 }
00497 }
00498
00499 fxtlog.close();
00500 delete gmxt;
00501 delete gmxtEd;
00502 cout << "Xt update finished" << endl;
00503 return 1;
00504 }
00505
00506 int XtMdcCalib::getHxtKey(int lay, int entr, int lr, int bin) const {
00507 int key;
00508
00509 key = ( (lay << HXTLAYER_INDEX) & HXTLAYER_MASK ) |
00510 ( (entr << HXTENTRA_INDEX) & HXTENTRA_MASK ) |
00511 ( (lr << HXTLR_INDEX) & HXTLR_MASK ) |
00512 ( (bin << HXTBIN_INDEX) & HXTBIN_MASK );
00513
00514 return key;
00515 }
00516
00517 void XtMdcCalib::fcnXT(Int_t &npar, Double_t *gin, Double_t &f,
00518 Double_t *par, Int_t iflag){
00519 unsigned int i;
00520 int ord;
00521 Double_t deta;
00522 Double_t chisq = 0.0;
00523 Double_t dfit;
00524
00525 for(i=0; i<TBINCEN.size(); i++){
00526 dfit = 0;
00527 for(ord=0; ord<=5; ord++){
00528 dfit += par[ord] * pow(TBINCEN[i], ord);
00529 }
00530 deta = (dfit - XMEAS[i]) / ERR[i];
00531 chisq += deta * deta;
00532 }
00533
00534 f = chisq;
00535 }
00536
00537 double XtMdcCalib::xtFun(double t, double xtpar[]){
00538 int ord;
00539 double dist = 0.0;
00540 double tmax = xtpar[6];
00541
00542 if(t < tmax ){
00543 for(ord=0; ord<=5; ord++){
00544 dist += xtpar[ord] * pow(t, ord);
00545 }
00546 }else{
00547 for(ord=0; ord<=5; ord++){
00548 dist += xtpar[ord] * pow(tmax, ord);
00549 }
00550 dist += xtpar[7] * (t - tmax);
00551 }
00552
00553 return dist;
00554 }
00555 void XtMdcCalib::fcnXtEdge(Int_t &npar, Double_t *gin, Double_t &f,
00556 Double_t *par, Int_t iflag){
00557 unsigned int i;
00558 Double_t deta;
00559 Double_t chisq = 0.0;
00560 Double_t dfit;
00561
00562 for(i=0; i<TBINCENED.size(); i++){
00563 dfit = par[0] * (TBINCENED[i] - Tmax) + Dmax;
00564 deta = (dfit - XMEASED[i]) / ERRED[i];
00565 chisq += deta * deta;
00566 }
00567
00568 f = chisq;
00569 }
00570