00001 #include "MdcCalibAlg/GrXtMdcCalib.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 <cstring>
00013 #include <cmath>
00014
00015 #include "TF1.h"
00016
00017 using namespace std;
00018
00019 double GrXtMdcCalib::DMAX;
00020 double GrXtMdcCalib::TMAX;
00021
00022 GrXtMdcCalib::GrXtMdcCalib(){
00023 m_maxNhit = 5000;
00024 m_fgIni = false;
00025 }
00026
00027 GrXtMdcCalib::~GrXtMdcCalib(){
00028 }
00029
00030 void GrXtMdcCalib::clear(){
00031 cout << "~GrXtMdcCalib" << endl;
00032 for(int lay=0; lay<MdcCalNLayer; lay++){
00033 for(int iEntr=0; iEntr<MdcCalNENTRXT; iEntr++){
00034 for(int iLR=0; iLR<MdcCalLR; iLR++){
00035 delete m_grxt[lay][iEntr][iLR];
00036 }
00037 }
00038 }
00039 delete m_haxis;
00040 delete m_fdXt;
00041
00042 MdcCalib::clear();
00043 }
00044
00045 void GrXtMdcCalib::initialize(TObjArray* hlist, IMdcGeomSvc* mdcGeomSvc,
00046 IMdcCalibFunSvc* mdcFunSvc, IMdcUtilitySvc* mdcUtilitySvc) {
00047 IMessageSvc* msgSvc;
00048 Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
00049 MsgStream log(msgSvc, "GrXtMdcCalib");
00050 log << MSG::INFO << "GrXtMdcCalib::initialize()" << endreq;
00051
00052 m_hlist = hlist;
00053 m_mdcGeomSvc = mdcGeomSvc;
00054 m_mdcFunSvc = mdcFunSvc;
00055 m_mdcUtilitySvc = mdcUtilitySvc;
00056
00057 MdcCalib::initialize(m_hlist, m_mdcGeomSvc, m_mdcFunSvc, m_mdcUtilitySvc);
00058
00059 int lay;
00060 int iLR;
00061 int iEntr;
00062 char hname[200];
00063
00064 m_fdXt = new TFolder("fdXtGr", "fdXtGr");
00065 m_hlist -> Add(m_fdXt);
00066
00067 m_haxis = new TH2F("axis", "", 50, 0, 300, 50, 0, 9);
00068 m_haxis -> SetStats(0);
00069 m_fdXt -> Add(m_haxis);
00070
00071 for(lay=0; lay<MdcCalNLayer; lay++){
00072 for(iEntr=0; iEntr<MdcCalNENTRXT; iEntr++){
00073 for(iLR=0; iLR<MdcCalLR; iLR++){
00074 m_nhit[lay][iEntr][iLR] = 0;
00075
00076 sprintf(hname, "grXt%02d_%02d_lr%01d", lay, iEntr, iLR);
00077 m_grxt[lay][iEntr][iLR] = new TGraph();
00078 m_grxt[lay][iEntr][iLR] -> SetName(hname);
00079 m_grxt[lay][iEntr][iLR] -> SetMarkerStyle(10);
00080 m_grxt[lay][iEntr][iLR] -> SetLineColor(10);
00081 m_fdXt -> Add(m_grxt[lay][iEntr][iLR]);
00082 }
00083 }
00084 }
00085 }
00086
00087 int GrXtMdcCalib::fillHist(MdcCalEvent* event){
00088 IMessageSvc* msgSvc;
00089 Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
00090 MsgStream log(msgSvc, "GrXtMdcCalib");
00091 log << MSG::DEBUG << "GrXtMdcCalib::fillHist()" << endreq;
00092
00093 MdcCalib::fillHist(event);
00094
00095
00096 bool esCutFg = event->getEsCutFlag();
00097 if( ! esCutFg ) return -1;
00098
00099 int i;
00100 int k;
00101 int lay;
00102 int iLR;
00103 int iEntr;
00104
00105 double dr;
00106 double dz;
00107 double doca;
00108 double tdr;
00109 double resi;
00110 double entrance;
00111
00112 int nhitlay;
00113 bool fgHitLay[MdcCalNLayer];
00114 bool fgTrk;
00115
00116 if(! m_fgIni){
00117 for(lay=0; lay<MdcCalNLayer; lay++){
00118 if(lay < 8) m_docaMax[lay] = m_param.maxDocaInner;
00119 else m_docaMax[lay] = m_param.maxDocaOuter;
00120 }
00121 m_fgIni = true;
00122 }
00123
00124 MdcCalRecTrk* rectrk;
00125 MdcCalRecHit* rechit;
00126
00127 int nhit;
00128 int ntrk = event -> getNTrk();
00129 for(i=0; i<ntrk; i++){
00130 fgTrk = true;
00131 rectrk = event->getRecTrk(i);
00132 nhit = rectrk -> getNHits();
00133
00134
00135 dr = rectrk->getDr();
00136 if(fabs(dr) > m_param.drCut) continue;
00137
00138
00139 dz = rectrk->getDz();
00140 if(fabs(dz) > m_param.dzCut) continue;
00141
00142 for(lay=0; lay<MdcCalNLayer; lay++){
00143 fgHitLay[lay] = false;
00144 }
00145
00146 for(k=0; k<nhit; k++){
00147 rechit = rectrk -> getRecHit(k);
00148 lay = rechit -> getLayid();
00149 doca = rechit -> getDocaInc();
00150 resi = rechit -> getResiInc();
00151 fgHitLay[lay] = true;
00152
00153
00154
00155
00156
00157 }
00158 if(! fgTrk) continue;
00159
00160 nhitlay = 0;
00161 for(lay=0; lay<MdcCalNLayer; lay++){
00162 if(fgHitLay[lay]) nhitlay++;
00163 }
00164 if(nhitlay < m_param.nHitLayCut) continue;
00165
00166 for(k=0; k<nhit; k++){
00167 rechit = rectrk -> getRecHit(k);
00168 lay = rechit -> getLayid();
00169 doca = rechit -> getDocaInc();
00170 resi = rechit -> getResiInc();
00171 iLR = rechit -> getLR();
00172 entrance = rechit -> getEntra();
00173 tdr = rechit -> getTdrift();
00174
00175 if( (fabs(doca) > m_docaMax[lay]) ||
00176 (fabs(resi) > m_param.resiCut[lay]) ){
00177 continue;
00178 }
00179
00180 if(0 == lay){
00181 if( ! fgHitLay[1] ) continue;
00182 } else if(42 == lay){
00183 if( ! fgHitLay[41] ) continue;
00184 } else{
00185 if( (!fgHitLay[lay-1]) && (!fgHitLay[lay+1]) ) continue;
00186 }
00187
00188 iEntr = m_mdcFunSvc -> getXtEntrIndex(entrance);
00189
00190 if(iLR < 2){
00191 if(m_nhit[lay][iEntr][iLR] < m_maxNhit){
00192 m_grxt[lay][iEntr][iLR] -> SetPoint(m_nhit[lay][iEntr][iLR],
00193 tdr, fabs(doca));
00194 m_nhit[lay][iEntr][iLR]++;
00195 }
00196 }
00197
00198 if(m_nhit[lay][iEntr][2] < m_maxNhit){
00199 m_grxt[lay][iEntr][2] -> SetPoint(m_nhit[lay][iEntr][2],
00200 tdr, fabs(doca));
00201 m_nhit[lay][iEntr][2]++;
00202 }
00203 }
00204 }
00205 return 1;
00206 }
00207
00208 int GrXtMdcCalib::updateConst(MdcCalibConst* calconst){
00209 IMessageSvc* msgSvc;
00210 Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
00211 MsgStream log(msgSvc, "GrXtMdcCalib");
00212 log << MSG::INFO << "GrXtMdcCalib::updateConst()" << endreq;
00213
00214 MdcCalib::updateConst(calconst);
00215
00216 int ord;
00217 double xtpar[MdcCalNLayer][MdcCalNENTRXT][MdcCalLR][8];
00218 TF1* fxtDr = new TF1("fxtDr", xtfun, 0, 300, 6);
00219 TF1* fxtEd = new TF1("fxtEd", xtedge, 150, 500, 1);
00220 if(1 == m_param.fixXtC0) fxtDr -> FixParameter(0, 0);
00221
00222 for(int lay=0; lay<MdcCalNLayer; lay++){
00223 for(int iEntr=0; iEntr<MdcCalNENTRXT; iEntr++){
00224 for(int lr=0; lr<MdcCalLR; lr++){
00225 m_fgFit[lay][iEntr][lr] = false;
00226 if(0 == m_param.fgCalib[lay]) continue;
00227
00228 if(m_nhit[lay][iEntr][lr] > 1000){
00229 TMAX = calconst -> getXtpar(lay, iEntr, lr, 6);
00230
00231 m_grxt[lay][iEntr][lr] -> Fit("fxtDr", "Q+", "", 0, TMAX);
00232
00233 for(ord=0; ord<6; ord++){
00234 xtpar[lay][iEntr][lr][ord] = fxtDr->GetParameter(ord);
00235 }
00236 xtpar[lay][iEntr][lr][6] = TMAX;
00237
00238 DMAX = 0.0;
00239 for(ord=0; ord<6; ord++) DMAX += xtpar[lay][iEntr][lr][ord] * pow(TMAX, ord);
00240
00241 m_grxt[lay][iEntr][lr] -> Fit("fxtEd", "Q+", "", TMAX, TMAX+300);
00242 xtpar[lay][iEntr][lr][7] = fxtEd->GetParameter(0);
00243 if(xtpar[lay][iEntr][lr][7] < 0.0) xtpar[lay][iEntr][lr][7] = 0.0;
00244 m_fgFit[lay][iEntr][lr] = true;
00245
00246
00247
00248
00249 }
00250
00251 }
00252 }
00253 }
00254
00255 ofstream fxtlog("xtlog");
00256 for(int lay=0; lay<MdcCalNLayer; lay++){
00257 for(int iEntr=0; iEntr<MdcCalNENTRXT; iEntr++){
00258 for(int lr=0; lr<MdcCalLR; lr++){
00259 fxtlog << setw(3) << lay << setw(3) << iEntr << setw(3) << lr;
00260
00261 int fgUpdate = -1;
00262 if(m_fgFit[lay][iEntr][lr]){
00263 fgUpdate = 1;
00264 for(ord=0; ord<8; ord++) calconst->resetXtpar(lay, iEntr, lr, ord, xtpar[lay][iEntr][lr][ord]);
00265 } else{
00266 int iEntrNew = findXtEntr(lay, iEntr, lr);
00267 if(-1 != iEntrNew){
00268 fgUpdate = 2;
00269 for(ord=0; ord<8; ord++){
00270 calconst->resetXtpar(lay, iEntr, lr, ord, xtpar[lay][iEntrNew][lr][ord]);
00271 }
00272 }
00273 }
00274 fxtlog << setw(3) << fgUpdate;
00275 for(ord=0; ord<8; ord++){
00276 double par = calconst -> getXtpar(lay, iEntr, lr, ord);
00277 fxtlog << setw(14) << par;
00278 }
00279 fxtlog << endl;
00280 }
00281 }
00282 }
00283 fxtlog.close();
00284
00285 cout << "Xt update finished. File xtlog was written." << endl;
00286 delete fxtDr;
00287 delete fxtEd;
00288 }
00289
00290 Double_t GrXtMdcCalib::xtfun(Double_t *x, Double_t *par){
00291 Double_t val = 0.0;
00292 for(Int_t ord=0; ord<6; ord++){
00293 val += par[ord] * pow(x[0], ord);
00294 }
00295 return val;
00296 }
00297
00298 Double_t GrXtMdcCalib::xtedge(Double_t *x, Double_t *par){
00299 double val = DMAX + (x[0] - TMAX) * par[0];
00300 return val;
00301 }
00302
00303 int GrXtMdcCalib::findXtEntr(int lay, int iEntr, int lr) const {
00304 int id0 = 8;
00305 int id1 = 9;
00306 int idmax = 17;
00307 int entrId = -1;
00308 if(iEntr <= id0){
00309 int id = -1;
00310 for(int i=iEntr; i<=id0; i++){
00311 if(m_fgFit[lay][i][lr]){
00312 id = i;
00313 break;
00314 }
00315 }
00316 if(-1 != id) entrId = id;
00317 else{
00318 for(int i=iEntr; i>=0; i--){
00319 if(m_fgFit[lay][i][lr]){
00320 id = i;
00321 break;
00322 }
00323 }
00324 entrId = id;
00325 }
00326 } else{
00327 int id = -1;
00328 for(int i=iEntr; i>=id1; i--){
00329 if(m_fgFit[lay][i][lr]){
00330 id = i;
00331 break;
00332 }
00333 }
00334 if(-1 != id) entrId = id;
00335 else{
00336 for(int i=iEntr; i<idmax; i++){
00337 if(m_fgFit[lay][i][lr]){
00338 id = i;
00339 break;
00340 }
00341 }
00342 entrId = id;
00343 }
00344 }
00345 if(-1 == entrId){
00346 cout << "find EntrId error " << "layer " << lay << " iEntr " << iEntr << " lr " << lr << endl;
00347 }
00348
00349 return entrId;
00350 }