00001 #include "MdcCalibAlg/PreT0MdcCalib.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 "GaudiKernel/SmartDataPtr.h"
00010 #include "GaudiKernel/IDataProviderSvc.h"
00011 #include "GaudiKernel/PropertyMgr.h"
00012
00013 #include "EvTimeEvent/RecEsTime.h"
00014 #include "EventModel/Event.h"
00015 #include "MdcRawEvent/MdcDigi.h"
00016 #include "Identifier/Identifier.h"
00017 #include "Identifier/MdcID.h"
00018 #include "TStyle.h"
00019
00020 #include <iostream>
00021 #include <fstream>
00022 #include <iomanip>
00023 #include <cstring>
00024
00025 #include <math.h>
00026 #include "TH2F.h"
00027 #include "TH1F.h"
00028 #include "TFolder.h"
00029 #include "TF1.h"
00030
00031 PreT0MdcCalib::PreT0MdcCalib(){
00032 m_nzbin = 11;
00033 }
00034
00035 PreT0MdcCalib::~PreT0MdcCalib(){
00036 }
00037
00038 void PreT0MdcCalib::clear(){
00039 int lay;
00040 int lr;
00041 int bin;
00042 for(lay=0; lay<MdcCalNLayer; lay++){
00043 for(lr=0; lr<MdcCalLR; lr++){
00044 delete m_hTrec[lay][lr];
00045 for(bin=0; bin<m_nzbin; bin++){
00046 delete m_hTrecZ[lay][lr][bin];
00047 }
00048 }
00049 }
00050 for(lay=0; lay<MdcCalNLayer; lay++){
00051 for(bin=0; bin<2; bin++) delete m_hTrecCosm[lay][bin];
00052 }
00053 delete m_fdTrec;
00054 delete m_fdTrecZ;
00055
00056 MdcCalib::clear();
00057 }
00058
00059 void PreT0MdcCalib::initialize(TObjArray* hlist, IMdcGeomSvc* mdcGeomSvc,
00060 IMdcCalibFunSvc* mdcFunSvc, IMdcUtilitySvc* mdcUtilitySvc) {
00061 IMessageSvc* msgSvc;
00062 Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
00063 MsgStream log(msgSvc, "PreT0MdcCalib");
00064 log << MSG::INFO << "PreT0MdcCalib::initialize()" << endreq;
00065
00066 m_hlist = hlist;
00067 m_mdcGeomSvc = mdcGeomSvc;
00068 m_mdcFunSvc = mdcFunSvc;
00069 m_mdcUtilitySvc = mdcUtilitySvc;
00070
00071
00072
00073 int lay;
00074 int lr;
00075 int bin;
00076 char hname[200];
00077
00078 m_fdTrec = new TFolder("Trec", "Trec");
00079 m_hlist->Add(m_fdTrec);
00080
00081 m_fdTrecZ = new TFolder("TrecZ", "TrecZ");
00082 m_hlist->Add(m_fdTrecZ);
00083
00084 for(lay=0; lay<MdcCalNLayer; lay++){
00085 for(lr=0; lr<MdcCalLR; lr++){
00086 if(0 == lr) sprintf(hname, "Trec%02d_L", lay);
00087 else if(1 == lr) sprintf(hname, "Trec%02d_R", lay);
00088 else sprintf(hname, "Trec%02d_C", lay);
00089
00090 if(lay < 8) m_hTrec[lay][lr] = new TH1F(hname, "", 100, 0, 400);
00091 else m_hTrec[lay][lr] = new TH1F(hname, "", 125, 0, 500);
00092 m_fdTrec -> Add(m_hTrec[lay][lr]);
00093 }
00094 }
00095
00096 for(lay=0; lay<MdcCalNLayer; lay++){
00097 for(int iud=0; iud<2; iud++){
00098 if(0 == iud) sprintf(hname, "TrecCosm%02d_Up", lay);
00099 else sprintf(hname, "TrecCosm%02d_Dw", lay);
00100 if(lay < 8) m_hTrecCosm[lay][iud] = new TH1F(hname, "", 100, 0, 400);
00101 else m_hTrecCosm[lay][iud] = new TH1F(hname, "", 125, 0, 500);
00102 m_fdTrec -> Add(m_hTrecCosm[lay][iud]);
00103 }
00104 }
00105
00106 for(lay=0; lay<MdcCalNLayer; lay++){
00107 for(lr=0; lr<MdcCalLR; lr++){
00108 for(bin=0; bin<m_nzbin; bin++){
00109 if(0 == lr) sprintf(hname, "Trec%02d_z%02d_L", lay, bin);
00110 else if(1 == lr) sprintf(hname, "Trec%02d_z%02d_R", lay, bin);
00111 else sprintf(hname, "Trec%02d_z%02d_C", lay, bin);
00112
00113 if(lay < 8) m_hTrecZ[lay][lr][bin] = new TH1F(hname, "", 100, 0, 400);
00114 else m_hTrecZ[lay][lr][bin] = new TH1F(hname, "", 125, 0, 500);
00115 m_fdTrecZ -> Add(m_hTrecZ[lay][lr][bin]);
00116 }
00117 }
00118 }
00119
00120 double zeast;
00121 double zwest;
00122 for(lay=0; lay<MdcCalNLayer; lay++){
00123 zwest = m_mdcGeomSvc->Wire(lay, 0)->Forward().z();
00124 zeast = m_mdcGeomSvc->Wire(lay, 0)->Backward().z();
00125 m_zwid[lay] = (zeast - zwest) / (double)m_nzbin;
00126
00127 if(lay < 8) m_vp[lay] = 220.0;
00128 else m_vp[lay] = 240.0;
00129
00130 if( 0 == (lay % 2) ){
00131 m_zst[lay] = zwest;
00132 } else{
00133 m_zst[lay] = zeast;
00134 }
00135 }
00136 }
00137
00138 int PreT0MdcCalib::fillHist(MdcCalEvent* event){
00139 IMessageSvc* msgSvc;
00140 Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
00141 MsgStream log(msgSvc, "PreT0MdcCalib");
00142 log << MSG::DEBUG << "PreT0MdcCalib::fillHist()" << endreq;
00143
00144
00145
00146 int i;
00147 int k;
00148 int lay;
00149 int cel;
00150 int lr;
00151 int bin;
00152
00153 double tdc;
00154 double traw;
00155 double trec;
00156 double tdr;
00157 double t0;
00158 double tp = 0.0;
00159 double tof = 0.0;
00160 double zhit;
00161 double zprop;
00162
00163 MdcCalRecTrk* rectrk;
00164 MdcCalRecHit* rechit;
00165
00166 IDataProviderSvc* eventSvc = NULL;
00167 Gaudi::svcLocator()->service("EventDataSvc", eventSvc);
00168
00169
00170 double tes = event->getTes();
00171
00172 bool esCutFg = event->getEsCutFlag();
00173 if( ! esCutFg ) return -1;
00174
00175 int nhit;
00176 int ntrk = event -> getNTrk();
00177 for(i=0; i<ntrk; i++){
00178 rectrk = event->getRecTrk(i);
00179 nhit = rectrk -> getNHits();
00180 for(k=0; k<nhit; k++){
00181 rechit = rectrk -> getRecHit(k);
00182 lay = rechit -> getLayid();
00183 cel = rechit -> getCellid();
00184 lr = rechit -> getLR();
00185 zhit = rechit -> getZhit();
00186 tdc = rechit -> getTdc();
00187 tdr = rechit -> getTdrift();
00188
00189 traw = tdc * MdcCalTdcCnv;
00190
00191
00192 zprop = fabs(zhit - m_zst[lay]);
00193 bin = (int)(zprop / m_zwid[lay]);
00194 tp = zprop / m_vp[lay];
00195 t0 = m_mdcFunSvc -> getT0(lay, cel);
00196 tof = rechit -> getTof();
00197
00198 double adc = rechit -> getQhit();
00199 double tw = m_mdcFunSvc->getTimeWalk(lay, adc);
00200
00201
00202 const MdcGeoWire* pWire = m_mdcGeomSvc -> Wire(lay, cel);
00203 double y1 = pWire -> Backward().y();
00204 double y2 = pWire -> Forward().y();
00205 double ymid = (y1+y2)*0.5;
00206
00207
00208
00209
00210
00211
00212 trec = traw - tp - tof - tes + m_param.timeShift - tw;
00213
00214 m_hTrec[lay][lr] -> Fill(trec);
00215 m_hTrec[lay][2] -> Fill(trec);
00216 if(ymid > 0) m_hTrecCosm[lay][0] -> Fill(trec);
00217 else m_hTrecCosm[lay][1] -> Fill(trec);
00218
00219 if(bin < m_nzbin){
00220 m_hTrecZ[lay][lr][bin] -> Fill(trec);
00221 m_hTrecZ[lay][2][bin] -> Fill(trec);
00222 }
00223 }
00224 }
00225 return 1;
00226 }
00227
00228 int PreT0MdcCalib::updateConst(MdcCalibConst* calconst){
00229 IMessageSvc* msgSvc;
00230 Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
00231 MsgStream log(msgSvc, "PreT0MdcCalib");
00232 log << MSG::DEBUG << "PreT0MdcCalib::updateConst()" << endreq;
00233
00234
00235
00236
00237 int lay;
00238 int wir;
00239 int lr;
00240 double t0Fit[MdcCalNLayer][MdcCalLR];
00241 double t0Cal[MdcCalNLayer][MdcCalLR];
00242 double tmax[MdcCalNLayer][MdcCalLR];
00243 double initT0 = m_param.initT0;
00244
00245 int fitTminFg[MdcCalNLayer][MdcCalLR];
00246 int fitTmaxFg[MdcCalNLayer][MdcCalLR];
00247 double chisq;
00248 double ndf;
00249 double chindfTmin[MdcCalNLayer][MdcCalLR];
00250 double chindfTmax[MdcCalNLayer][MdcCalLR];
00251 char funname[200];
00252
00253
00254 TF1* ftminCosm[MdcCalNLayer][2];
00255 double t0FitCosm[MdcCalNLayer][2];
00256
00257 TF1* ftmin[MdcCalNLayer][MdcCalLR];
00258
00259
00260 for(lay=0; lay<MdcCalNLayer; lay++){
00261 for(lr=0; lr<MdcCalLR; lr++){
00262 fitTminFg[lay][lr] = 0;
00263 chindfTmin[lay][lr] = -1;
00264 sprintf(funname, "ftmin%02d_%d", lay, lr);
00265 ftmin[lay][lr] = new TF1(funname, funTmin, 0, 150, 6);
00266
00267 if(1 == m_param.fgCalib[lay]){
00268 Stat_t nEntryTot = 0;
00269 for(int ibin=1; ibin<=25; ibin++){
00270 Stat_t entry = m_hTrec[lay][lr]->GetBinContent(ibin);
00271 nEntryTot += entry;
00272 }
00273 double c0Ini = (double)nEntryTot / 25.0;
00274 double c1Ini = (m_hTrec[lay][lr]->GetMaximum()) - c0Ini;
00275
00276 ftmin[lay][lr] -> SetParameter(0, c0Ini);
00277 ftmin[lay][lr] -> SetParameter(1, c1Ini);
00278 ftmin[lay][lr] -> SetParameter(2, 0);
00279 ftmin[lay][lr] -> SetParameter(4, initT0);
00280 ftmin[lay][lr] -> SetParameter(5, 2);
00281
00282 m_hTrec[lay][lr] -> Fit(funname, "Q", "",
00283 m_param.tminFitRange[lay][0],
00284 m_param.tminFitRange[lay][1]);
00285 gStyle -> SetOptFit(11);
00286
00287
00288 chisq = ftmin[lay][lr]->GetChisquare();
00289 ndf = ftmin[lay][lr]->GetNDF();
00290 chindfTmin[lay][lr] = chisq / ndf;
00291
00292 if(chindfTmin[lay][lr] < m_param.tminFitChindf){
00293 fitTminFg[lay][lr] = 1;
00294 t0Fit[lay][lr] = ftmin[lay][lr]->GetParameter(4);
00295
00296 t0Fit[lay][lr] += m_param.t0Shift;
00297 t0Cal[lay][lr] = t0Fit[lay][lr] - m_param.timeShift;
00298 }
00299 }
00300
00301 if(0 == fitTminFg[lay][lr]){
00302 wir = m_mdcGeomSvc->Wire(lay, 0)->Id();
00303 t0Cal[lay][lr] = calconst->getT0(wir);
00304 t0Fit[lay][lr] = t0Cal[lay][lr] + m_param.timeShift;
00305 }
00306 }
00307
00308 for(int iud=0; iud<2; iud++){
00309 sprintf(funname, "ftminCosm_%02d_%d", lay, iud);
00310 ftminCosm[lay][iud] = new TF1(funname, funTmin, 0, 150, 6);
00311 ftminCosm[lay][iud] -> SetParameter(0, 0);
00312 ftminCosm[lay][iud] -> SetParameter(4, initT0);
00313 ftminCosm[lay][iud] -> SetParameter(5, 1);
00314 m_hTrecCosm[lay][iud] -> Fit(funname, "Q", "",
00315 m_param.tminFitRange[lay][0],
00316 m_param.tminFitRange[lay][1]);
00317 gStyle -> SetOptFit(11);
00318 t0FitCosm[lay][iud] += m_param.t0Shift;
00319 t0FitCosm[lay][iud] = ftminCosm[lay][iud]->GetParameter(4);
00320 }
00321 }
00322
00323
00324 TF1* ftmax[MdcCalNLayer][MdcCalLR];
00325 for(lay=0; lay<MdcCalNLayer; lay++){
00326 for(lr=0; lr<MdcCalLR; lr++){
00327 fitTmaxFg[lay][lr] = 0;
00328 chindfTmax[lay][lr] = -1;
00329 sprintf(funname, "ftmax%02d_%d", lay, lr);
00330 ftmax[lay][lr] = new TF1(funname, funTmax, 250, 500, 4);
00331
00332 if(1 == m_param.fgCalib[lay]){
00333 ftmax[lay][lr] -> SetParameter(2, m_param.initTm[lay]);
00334 ftmax[lay][lr] -> SetParameter(3, 10);
00335 m_hTrec[lay][lr] -> Fit(funname, "Q+", "",
00336 m_param.tmaxFitRange[lay][0],
00337 m_param.tmaxFitRange[lay][1]);
00338 gStyle -> SetOptFit(11);
00339 chisq = ftmax[lay][lr]->GetChisquare();
00340 ndf = ftmax[lay][lr]->GetNDF();
00341 chindfTmax[lay][lr] = chisq / ndf;
00342 if(chindfTmax[lay][lr] < m_param.tmaxFitChindf){
00343 fitTmaxFg[lay][lr] = 1;
00344 tmax[lay][lr] = ftmax[lay][lr]->GetParameter(2);
00345 }
00346 }
00347
00348 if(0 == fitTmaxFg[lay][lr]){
00349 tmax[lay][lr] = (calconst->getXtpar(lay, 0, lr, 6)) + t0Fit[lay][2];
00350 }
00351 }
00352 }
00353
00354
00355 ofstream ft0("preT0.dat");
00356 for(lay=0; lay<MdcCalNLayer; lay++){
00357 ft0 << setw(5) << lay << setw(3) << fitTminFg[lay][2]
00358 << setw(15) << t0Cal[lay][2] << setw(15) << t0Fit[lay][2]
00359 << setw(15) << chindfTmin[lay][2] << endl;
00360 }
00361 ft0 << endl;
00362 for(lay=0; lay<MdcCalNLayer; lay++){
00363 ft0 << setw(5) << lay
00364 << setw(3) << fitTmaxFg[lay][0] << setw(10) << tmax[lay][0]
00365 << setw(10) << chindfTmax[lay][0]
00366 << setw(3) << fitTmaxFg[lay][1] << setw(10) << tmax[lay][1]
00367 << setw(10) << chindfTmax[lay][1]
00368 << setw(3) << fitTmaxFg[lay][2] << setw(10) << tmax[lay][2]
00369 << setw(10) << chindfTmax[lay][2]
00370 << setw(10) << tmax[lay][0] - t0Fit[lay][2]
00371 << setw(10) << tmax[lay][1] - t0Fit[lay][2]
00372 << setw(10) << tmax[lay][2] - t0Fit[lay][2]
00373 << endl;
00374 }
00375 ft0.close();
00376 cout << "preT0.dat was written." << endl;
00377
00378
00379 ofstream ft0cosm("cosmicT0.dat");
00380 for(lay=0; lay<MdcCalNLayer; lay++){
00381 ft0cosm << setw(5) << lay << setw(15) << t0Fit[lay][2]
00382 << setw(15) << t0FitCosm[lay][0] << setw(15) << t0FitCosm[lay][1] << endl;
00383 }
00384 ft0cosm.close();
00385
00386
00387 int i;
00388 int nwire = m_mdcGeomSvc -> getWireSize();
00389 for(i=0; i<nwire; i++){
00390 lay = m_mdcGeomSvc -> Wire(i) -> Layer();
00391 if(1 == m_param.fgCalib[lay]){
00392 calconst -> resetT0(i, t0Cal[lay][2]);
00393 calconst -> resetDelT0(i, 0.0);
00394 }
00395 }
00396
00397
00398 if(m_param.preT0SetTm){
00399 int iEntr;
00400 double tm;
00401 for(lay=0; lay<MdcCalNLayer; lay++){
00402 if(1 != m_param.fgCalib[lay]) continue;
00403
00404 for(iEntr=0; iEntr<MdcCalNENTRXT; iEntr++){
00405 for(lr=0; lr<MdcCalLR; lr++){
00406 tm = tmax[lay][lr] - t0Fit[lay][2];
00407 if( (tmax[lay][lr] > m_param.tmaxFitRange[lay][0]) &&
00408 (tmax[lay][lr] < m_param.tmaxFitRange[lay][1]) ){
00409 calconst -> resetXtpar(lay, iEntr, lr, 6, tm);
00410 }
00411 }
00412 }
00413 }
00414 }
00415
00416
00417 int bin;
00418 double sdpar = m_param.initSigma;
00419 for(lay=0; lay<MdcCalNLayer; lay++){
00420 for(int iEntr=0; iEntr<MdcCalNENTRSD; iEntr++){
00421 for(lr=0; lr<2; lr++){
00422 for(bin=0; bin<MdcCalSdNBIN; bin++){
00423 calconst -> resetSdpar(lay, iEntr, lr, bin, sdpar);
00424 }
00425 }
00426 }
00427 }
00428
00429
00430 for(lay=0; lay<MdcCalNLayer; lay++){
00431 for(lr=0; lr<MdcCalLR; lr++){
00432 delete ftmin[lay][lr];
00433 delete ftmax[lay][lr];
00434 }
00435 }
00436 return 1;
00437 }
00438 Double_t PreT0MdcCalib::funTmin(Double_t* x, Double_t* par){
00439 Double_t fitval;
00440 fitval = par[0] + par[1]*exp( -par[2]*(x[0]-par[3]) ) /
00441 ( 1 + exp( -(x[0]-par[4])/par[5] ));
00442 return fitval;
00443 }
00444
00445 Double_t PreT0MdcCalib::funTmax(Double_t* x, Double_t* par){
00446 Double_t fitval;
00447 fitval = par[0] + par[1] / (1 + exp((x[0]-par[2])/par[3]));
00448 return fitval;
00449 }
00450