00001 #include "MdcCalibAlg/PreXtMdcCalib.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 <string>
00024 #include <cstring>
00025 #include <cmath>
00026
00027 #include "TF1.h"
00028
00029 using namespace Event;
00030
00031 PreXtMdcCalib::PreXtMdcCalib(){
00032 m_fgIniTm = false;
00033 }
00034
00035 PreXtMdcCalib::~PreXtMdcCalib(){
00036 }
00037
00038 void PreXtMdcCalib::clear(){
00039 for(int lay=0; lay<MdcCalNLayer; lay++){
00040 delete m_grXt[lay];
00041 delete m_htrec[lay];
00042 }
00043 delete m_haxis;
00044 delete m_fdPreXt;
00045 }
00046
00047 void PreXtMdcCalib::initialize(TObjArray* hlist, IMdcGeomSvc* mdcGeomSvc,
00048 IMdcCalibFunSvc* mdcFunSvc, IMdcUtilitySvc* mdcUtilitySvc) {
00049 IMessageSvc* msgSvc;
00050 Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
00051 MsgStream log(msgSvc, "PreXtMdcCalib");
00052 log << MSG::INFO << "PreXtMdcCalib::initialize()" << endreq;
00053
00054 m_hlist = hlist;
00055 m_mdcGeomSvc = mdcGeomSvc;
00056 m_mdcFunSvc = mdcFunSvc;
00057 m_mdcUtilitySvc = mdcUtilitySvc;
00058
00059 int lay;
00060
00061 m_nWire = m_mdcGeomSvc -> getWireSize();
00062 m_nLayer = m_mdcGeomSvc -> getLayerSize();
00063
00064 m_fdPreXt = new TFolder("PreXt", "PreXt");
00065 m_hlist->Add(m_fdPreXt);
00066
00067 m_fdNhit = new TFolder("XtNhit", "XtNhit");
00068 m_hlist->Add(m_fdNhit);
00069
00070 m_haxis = new TH2F("axis", "", 50, 0, 300, 50, 0, 9);
00071 m_haxis -> SetStats(0);
00072 m_fdPreXt -> Add(m_haxis);
00073
00074 char hname[200];
00075 for(lay=0; lay<MdcCalNLayer; lay++){
00076 sprintf(hname, "trec%02d", lay);
00077 m_htrec[lay] = new TH1F(hname, "", 310, -20, 600);
00078 m_fdPreXt -> Add(m_htrec[lay]);
00079 }
00080
00081 m_nhitTot = new TH1F("nhitTot", "", 43, -0.5, 42.5);
00082 m_fdNhit -> Add(m_nhitTot);
00083
00084 for(lay=0; lay<MdcCalNLayer; lay++){
00085 sprintf(hname, "nhitBin%02d", lay);
00086 m_nhitBin[lay] = new TH1F(hname, "", 40, 5.0, 405.0);
00087 m_fdNhit -> Add(m_nhitBin[lay]);
00088 }
00089
00090
00091 int bin;
00092 m_nXtBin = 40;
00093 double twid = 10.0;
00094 for(bin=0; bin<m_nXtBin; bin++) m_tbin[bin] = (double)(bin+1) * twid;
00095
00096 for(lay=0; lay<MdcCalNLayer; lay++){
00097 m_nTot[lay] = 0;
00098 for(bin=0; bin<m_nXtBin; bin++){
00099 m_nEntries[lay][bin] = 0;
00100 }
00101 }
00102 }
00103
00104 int PreXtMdcCalib::fillHist(MdcCalEvent* event){
00105 IMessageSvc* msgSvc;
00106 Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
00107 MsgStream log(msgSvc, "PreXtMdcCalib");
00108 log << MSG::DEBUG << "PreXtMdcCalib::fillHist()" << endreq;
00109
00110 int bin;
00111 int lay;
00112 int digiId;
00113 unsigned fgOverFlow;
00114 double tdc;
00115 double adc;
00116 double traw;
00117 double trec;
00118 Identifier id;
00119
00120 IDataProviderSvc* eventSvc = NULL;
00121 Gaudi::svcLocator()->service("EventDataSvc", eventSvc);
00122
00123 double xtpar[8];
00124 if(! m_fgIniTm){
00125 for(lay=0; lay<MdcCalNLayer; lay++){
00126 m_t0[lay] = m_mdcFunSvc -> getT0(lay, 0);
00127 m_mdcFunSvc->getXtpar(lay, 0, 0, xtpar);
00128 m_tm[lay] = xtpar[6];
00129 }
00130 m_fgIniTm = true;
00131 }
00132
00133
00134 double tes = -9999.0;
00135 int esTimeflag = -1;
00136 SmartDataPtr<RecEsTimeCol> aevtimeCol(eventSvc,"/Event/Recon/RecEsTimeCol");
00137 if( (!aevtimeCol) || (aevtimeCol->size()==0) ){
00138 tes = -9999.0;
00139 esTimeflag = -1;
00140 }else{
00141 RecEsTimeCol::iterator iter_evt = aevtimeCol->begin();
00142 for(; iter_evt!=aevtimeCol->end(); iter_evt++){
00143 tes = (*iter_evt)->getTest();
00144 esTimeflag = (*iter_evt)->getStat();
00145 }
00146 }
00147 bool flagTes = false;
00148 for(int iEs=0; iEs<m_param.nEsFlag; iEs++){
00149 if(esTimeflag == m_param.esFlag[iEs]){
00150 flagTes = true;
00151 break;
00152 }
00153 }
00154 if( (!flagTes) || (tes < m_param.tesMin) || (tes > m_param.tesMax) ) return -1;
00155
00156
00157 SmartDataPtr<MdcDigiCol> mdcDigiCol(eventSvc,"/Event/Digi/MdcDigiCol");
00158 if (!mdcDigiCol) {
00159 log << MSG::FATAL << "Could not find event" << endreq;
00160 }
00161
00162 MdcDigiCol::iterator iter = mdcDigiCol->begin();
00163 digiId = 0;
00164 for(; iter != mdcDigiCol->end(); iter++, digiId++) {
00165 MdcDigi *aDigi = (*iter);
00166 id = (aDigi)->identify();
00167
00168 lay = MdcID::layer(id);
00169
00170 tdc = (aDigi) -> getTimeChannel();
00171 adc = (aDigi) -> getChargeChannel();
00172 fgOverFlow = (aDigi) -> getOverflow();
00173 traw = tdc * MdcCalTdcCnv;
00174 trec = traw - tes - m_t0[lay];
00175
00176
00177
00178
00179 if ( ((fgOverFlow & 1) !=0 ) || ((fgOverFlow & 12) != 0) ||
00180 (aDigi->getTimeChannel() == 0x7FFFFFFF) ||
00181 (aDigi->getChargeChannel() == 0x7FFFFFFF) ) continue;
00182
00183
00184
00185 if(trec > 0){
00186 if(trec < m_tm[lay]){
00187 m_nTot[lay]++;
00188 m_htrec[lay]->Fill(trec);
00189 m_nhitTot->Fill(lay);
00190 }
00191 for(bin=0; bin<m_nXtBin; bin++){
00192 if(trec < m_tbin[bin]){
00193 m_nEntries[lay][bin]++;
00194 m_nhitBin[lay]->Fill(m_tbin[bin]);
00195 }
00196 }
00197 }
00198 }
00199 return 1;
00200 }
00201
00202 int PreXtMdcCalib::updateConst(MdcCalibConst* calconst){
00203 IMessageSvc* msgSvc;
00204 Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
00205 MsgStream log(msgSvc, "PreXtMdcCalib");
00206 log << MSG::DEBUG << "PreXtMdcCalib::updateConst()" << endreq;
00207
00208 int lay;
00209 int bin;
00210 int iLR;
00211 int iEntr;
00212 int ord;
00213
00214 double layRad;
00215 double ncel;
00216 double pi = 3.141592653;
00217
00218 double dm;
00219 double dist[40];
00220 double xtpar[6];
00221 char hname[200];
00222
00223 TF1* funXt = new TF1("funXt", xtfun, 0, 300, 6);
00224 funXt -> FixParameter(0, 0.0);
00225 funXt -> SetParameter(1, 0.03);
00226 funXt -> SetParameter(2, 0.0);
00227 funXt -> SetParameter(3, 0.0);
00228 funXt -> SetParameter(4, 0.0);
00229 funXt -> SetParameter(5, 0.0);
00230
00231 ofstream fxtlog("preXtpar.dat");
00232 for(lay=0; lay<MdcCalNLayer; lay++){
00233 sprintf(hname, "grPreXt%02d", lay);
00234 m_grXt[lay] = new TGraph();
00235 m_grXt[lay] -> SetName(hname);
00236 m_grXt[lay] -> SetMarkerStyle(20);
00237 m_fdPreXt -> Add(m_grXt[lay]);
00238
00239 layRad = m_mdcGeomSvc -> Layer(lay) -> Radius();
00240 ncel = m_mdcGeomSvc -> Layer(lay) -> NCell();
00241 dm = pi * layRad / ncel;
00242
00243 fxtlog << "layer " << lay << endl;
00244 for(bin=0; bin<m_nXtBin; bin++){
00245 dist[bin] = dm * m_nEntries[lay][bin] / m_nTot[lay];
00246 m_grXt[lay] -> SetPoint(bin, m_tbin[bin], dist[bin]);
00247 fxtlog << setw(4) << bin << setw(15) << m_tbin[bin]
00248 << setw(15) << dist[bin] << setw(15) << dm
00249 << setw(10) << m_nEntries[lay][bin]
00250 << setw(10) << m_nTot[lay] << endl;
00251
00252 if(m_tbin[bin] >= m_tm[lay]) break;
00253 }
00254
00255 if(1 == m_param.fgCalib[lay]){
00256 m_grXt[lay] -> Fit(funXt, "Q", "", 0.0, m_tm[lay]);
00257 for(ord=0; ord<6; ord++){
00258 xtpar[ord] = funXt -> GetParameter(ord);
00259 }
00260
00261 for(iEntr=0; iEntr<MdcCalNENTRXT; iEntr++){
00262 for(iLR=0; iLR<MdcCalLR; iLR++){
00263 for(ord=0; ord<6; ord++){
00264 calconst -> resetXtpar(lay, iEntr, iLR, ord, xtpar[ord]);
00265 }
00266 }
00267 }
00268 } else{
00269 for(ord=0; ord<6; ord++) xtpar[ord] = calconst->getXtpar(lay, 0, 0, ord);
00270 }
00271
00272 for(ord=0; ord<6; ord++) fxtlog << setw(14) << xtpar[ord];
00273 fxtlog << setw(10) << m_tm[lay] << " 0" << endl;
00274 }
00275 fxtlog.close();
00276 cout << "preXt.dat was written." << endl;
00277
00278 delete funXt;
00279 return 1;
00280 }
00281
00282 Double_t PreXtMdcCalib::xtfun(Double_t *x, Double_t *par){
00283 Double_t val = 0.0;
00284 for(Int_t ord=0; ord<6; ord++){
00285 val += par[ord] * pow(x[0], ord);
00286 }
00287 return val;
00288 }
00289