00001 #include "MdcCalibAlg/WrMdcCalib.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 <iomanip>
00011 #include <cstring>
00012
00013 WrMdcCalib::WrMdcCalib(){
00014 for(int lay=0; lay<MdcCalNLayer; lay++){
00015 if(lay < 8){
00016 m_docaMin[lay] = 1.2;
00017 m_docaMax[lay] = 4.8;
00018 } else{
00019 m_docaMin[lay] = 1.6;
00020 m_docaMax[lay] = 6.4;
00021 }
00022 }
00023 }
00024
00025 WrMdcCalib::~WrMdcCalib(){
00026 }
00027
00028 void WrMdcCalib::clear(){
00029 for(int i=0; i<MdcCalTotCell; i++){
00030 delete m_hleft[i];
00031 delete m_hright[i];
00032 }
00033 delete m_hdwxtot;
00034 delete m_hddwx;
00035 delete m_hdwytot;
00036 delete m_hddwy;
00037 delete m_hLrResiSum;
00038 delete m_hLrResiSub;
00039 delete m_fdWire;
00040 delete m_fdResiWire;
00041 MdcCalib::clear();
00042 }
00043
00044 void WrMdcCalib::initialize(TObjArray* hlist, IMdcGeomSvc* mdcGeomSvc,
00045 IMdcCalibFunSvc* mdcFunSvc, IMdcUtilitySvc* mdcUtilitySvc) {
00046 IMessageSvc *msgSvc;
00047 Gaudi::svcLocator()->service("MessageSvc", msgSvc);
00048 MsgStream log(msgSvc, "WrMdcCalib");
00049 log << MSG::INFO << "WrMdcCalib::initialize()" << endreq;
00050
00051 m_hlist = hlist;
00052 m_mdcGeomSvc = mdcGeomSvc;
00053 m_mdcFunSvc = mdcFunSvc;
00054 m_mdcUtilitySvc = mdcUtilitySvc;
00055
00056 MdcCalib::initialize(m_hlist, m_mdcGeomSvc, m_mdcFunSvc, m_mdcUtilitySvc);
00057
00058 int i;
00059 int nwire;
00060 int lay;
00061 int cel;
00062 char hname[200];
00063
00064 m_fdWire = new TFolder("WireCor", "WireCor");
00065 m_hlist->Add(m_fdWire);
00066
00067 m_fdResiWire = new TFolder("ResiWire", "ResiWire");
00068 m_hlist->Add(m_fdResiWire);
00069
00070 nwire = m_mdcGeomSvc -> getWireSize();
00071 for(i=0; i<nwire; i++){
00072 lay = m_mdcGeomSvc -> Wire(i) -> Layer();
00073 cel = m_mdcGeomSvc -> Wire(i) -> Cell();
00074
00075 sprintf(hname, "h%04d_L%02d_%03d_Left", i, lay, cel);
00076 m_hleft[i] = new TH1F(hname, hname, 300, -1.5, 1.5);
00077 m_fdWire->Add(m_hleft[i]);
00078
00079 sprintf(hname, "h%04d_L%02d_%03d_Right", i, lay, cel);
00080 m_hright[i] = new TH1F(hname, hname, 300, -1.5, 1.5);
00081 m_fdWire->Add(m_hright[i]);
00082 }
00083
00084 m_hdwxtot = new TH1F("dwXtot", "", 100, -0.5, 0.5);
00085 m_fdResiWire->Add(m_hdwxtot);
00086
00087 m_hddwx = new TH1F("ddwX", "", 100, -0.5, 0.5);
00088 m_fdResiWire->Add(m_hddwx);
00089
00090 m_hdwytot = new TH1F("dwYtot", "", 100, -0.5, 0.5);
00091 m_fdResiWire->Add(m_hdwytot);
00092
00093 m_hddwy = new TH1F("ddwY", "", 100, -0.5, 0.5);
00094 m_fdResiWire->Add(m_hddwy);
00095
00096 m_hLrResiSum = new TH1F("LrResiSum", "", 200, -0.5, 0.5);
00097 m_fdResiWire->Add(m_hLrResiSum);
00098
00099 m_hLrResiSub = new TH1F("LrResiSub", "", 200, -0.5, 0.5);
00100 m_fdResiWire->Add(m_hLrResiSub);
00101 }
00102
00103 int WrMdcCalib::fillHist(MdcCalEvent* event){
00104 IMessageSvc *msgSvc;
00105 Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
00106 MsgStream log(msgSvc, "WrMdcCalib");
00107 log << MSG::DEBUG << "WrMdcCalib::fillHist()" << endreq;
00108
00109 MdcCalib::fillHist(event);
00110
00111
00112 bool esCutFg = event->getEsCutFlag();
00113 if( ! esCutFg ) return -1;
00114
00115 int i;
00116 int k;
00117 int ntrk;
00118 int nhit;
00119 int stat;
00120
00121 int lay;
00122 int cel;
00123 int wir;
00124 int lr;
00125 double dmeas;
00126 double doca;
00127 double resi;
00128
00129 bool fgHitLay[MdcCalNLayer];
00130
00131 MdcCalRecTrk* rectrk;
00132 MdcCalRecHit* rechit;
00133
00134 ntrk = event->getNTrk();
00135 log << MSG::DEBUG << "number of tracks: " << ntrk << endreq;
00136
00137 for(i=0; i<ntrk; i++){
00138 rectrk = event -> getRecTrk(i);
00139 nhit = rectrk -> getNHits();
00140
00141
00142 double dr = rectrk->getDr();
00143 if(fabs(dr) > m_param.drCut) continue;
00144
00145
00146 double dz = rectrk->getDz();
00147 if(fabs(dz) > m_param.dzCut) continue;
00148
00149 for(lay=0; lay<MdcCalNLayer; lay++) fgHitLay[lay] = false;
00150 for(k=0; k<nhit; k++){
00151 rechit = rectrk -> getRecHit(k);
00152 lay = rechit -> getLayid();
00153 fgHitLay[lay] = true;
00154 }
00155
00156 int nhitlay = 0;
00157 for(lay=0; lay<MdcCalNLayer; lay++) if(fgHitLay[lay]) nhitlay++;
00158 if(nhitlay < m_param.nHitLayCut) continue;
00159
00160 log << MSG::DEBUG << "number of hits: " << nhit << endreq;
00161 for(k=0; k<nhit; k++){
00162 rechit = rectrk -> getRecHit(k);
00163 lay = rechit -> getLayid();
00164 cel = rechit -> getCellid();
00165 wir = m_mdcGeomSvc ->Wire(lay, cel)->Id();
00166 lr = rechit -> getLR();
00167 doca = rechit -> getDocaExc();
00168 dmeas = rechit -> getDmeas();
00169 resi = rechit -> getResiExc();
00170 stat = rechit -> getStat();
00171
00172 if(1 != stat) continue;
00173
00174 if( (fabs(doca) < m_docaMin[lay]) ||
00175 (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 if(wir < (m_mdcGeomSvc->getWireSize()) ){
00189 if( 0 == lr ){
00190 m_hleft[wir] -> Fill(resi);
00191 }else if( 1 == lr ){
00192 m_hright[wir] -> Fill(resi);
00193 }
00194 }else{
00195 std::cout << "wir: " << wir << std::endl;
00196 }
00197 }
00198 }
00199 return 1;
00200 }
00201
00202 int WrMdcCalib::updateConst(MdcCalibConst* calconst){
00203 IMessageSvc *msgSvc;
00204 Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
00205 MsgStream log(msgSvc, "WrMdcCalib");
00206 log << MSG::INFO << "WrMdcCalib::updateConst()" << endreq;
00207
00208 MdcCalib::updateConst(calconst);
00209
00210 int i;
00211 int lay;
00212 int cel;
00213 int nwire = m_mdcGeomSvc -> getWireSize();
00214 double dwphi;
00215 double resiLrSum;
00216 double wpos[6];
00217 double xx;
00218 double yy;
00219 double rr;
00220 double dx;
00221 double dy;
00222 double wphi;
00223
00224 Stat_t entry_l;
00225 double mean_l;
00226 Stat_t entry_r;
00227 double mean_r;
00228 const MdcGeoWire* pWire;
00229
00230 double ddw[MdcCalTotCell][6];
00231 double dwinput[MdcCalTotCell][6];
00232 string strtmp;
00233 ifstream fwpinput(m_param.wpcFile.c_str());
00234 for(i=0; i<7; i++) fwpinput >> strtmp;
00235 for(i=0; i<nwire; i++){
00236 fwpinput >> strtmp >> dwinput[i][0] >> dwinput[i][1] >> dwinput[i][2]
00237 >> dwinput[i][3] >> dwinput[i][4] >> dwinput[i][5];
00238 }
00239 fwpinput.close();
00240
00241 std::cout << "totwire: " << nwire << std::endl;
00242 for(i=0; i<nwire; i++){
00243 pWire = m_mdcGeomSvc -> Wire(i);
00244 lay = pWire -> Layer();
00245 cel = pWire -> Cell();
00246
00247 if(1 == m_param.fgCalib[lay]){
00248 entry_l = m_hleft[i] -> GetEntries();
00249 mean_l = m_hleft[i] -> GetMean();
00250
00251 entry_r = m_hright[i] -> GetEntries();
00252 mean_r = m_hright[i] -> GetMean();
00253
00254 dwphi = 0.5 * (mean_l - mean_r);
00255 } else{
00256 dwphi = 0.0;
00257 }
00258
00259 resiLrSum = 0.5 * (mean_l + mean_r);
00260 m_hLrResiSum->Fill(resiLrSum);
00261 m_hLrResiSub->Fill(dwphi);
00262
00263 wpos[0] = pWire -> Backward().x();
00264 wpos[1] = pWire -> Backward().y();
00265 wpos[3] = pWire -> Forward().x();
00266 wpos[4] = pWire -> Forward().y();
00267
00268 xx = 0.5 * (wpos[0] + wpos[3]);
00269 yy = 0.5 * (wpos[1] + wpos[4]);
00270 rr = sqrt( (xx * xx) + (yy * yy) );
00271
00272 if( yy >= 0 ) wphi = acos(xx / rr);
00273 else wphi = PI2 - acos(xx / rr);
00274
00275 dx = -1.0 * dwphi * sin(wphi);
00276 dy = dwphi * cos(wphi);
00277
00278 ddw[i][0] = dx;
00279 ddw[i][1] = dy;
00280 ddw[i][2] = 0.0;
00281 ddw[i][3] = dx;
00282 ddw[i][4] = dy;
00283 ddw[i][5] = 0.0;
00284 }
00285
00286 ofstream fwpc("WirePosCalib_new.txt");
00287 fwpc << setw(5) << "wireId" << setw(15) << "dx_East(mm)"
00288 << setw(15) << "dy_East(mm)" << setw(15) << "dz_East(mm)"
00289 << setw(15) << "dx_West(mm)" << setw(15) << "dy_West(mm)"
00290 << setw(15) << "dz_West(mm)" << endl;
00291
00292 ofstream fdw("dw.txt");
00293 fdw << setw(5) << "wireId" << setw(15) << "ddx_East(mm)"
00294 << setw(15) << "ddy_East(mm)" << setw(15) << "ddz_East(mm)"
00295 << setw(15) << "ddx_West(mm)" << setw(15) << "ddy_West(mm)"
00296 << setw(15) << "ddz_West(mm)" << endl;
00297
00298 int k;
00299 double dwtot[6];
00300 for(i=0; i<nwire; i++){
00301 for(k=0; k<6; k++) dwtot[k] = dwinput[i][k] + ddw[i][k];
00302 fwpc << setw(5) << i;
00303 for(k=0; k<6; k++) fwpc << setw(15) << dwtot[k];
00304 fwpc << endl;
00305
00306 fdw << setw(5) << i;
00307 for(k=0; k<6; k++) fdw << setw(15) << ddw[i][k];
00308 fdw << endl;
00309 m_hdwxtot->Fill(dwtot[0]);
00310 m_hddwx->Fill(ddw[i][0]);
00311 m_hdwytot->Fill(dwtot[1]);
00312 m_hddwy->Fill(ddw[i][1]);
00313 }
00314 fwpc.close();
00315 fdw.close();
00316 return 1;
00317 }
00318
