BOSS6.5.5: WrMdcCalib Class Reference

00013                       {
00014      for(int lay=0; lay<MdcCalNLayer; lay++){
00015           if(lay < 8){
00016                m_docaMin[lay] = 1.2; // mm
00017                m_docaMax[lay] = 4.8; // mm
00018           } else{
00019                m_docaMin[lay] = 1.6; // mm
00020                m_docaMax[lay] = 6.4; // mm
00021           }
00022      }
00023 }

00028                       {
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 }

00102                                           {
00103      IMessageSvc *msgSvc;
00104      Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
00105      MsgStream log(msgSvc, "WrMdcCalib");
00106      log << MSG::DEBUG << "WrMdcCalib::fillHist()" << endreq;
00107 
00108      MdcCalib::fillHist(event);
00109 
00110      // get EsTimeCol
00111      bool esCutFg = event->getEsCutFlag();
00112      if( ! esCutFg ) return -1;
00113 
00114      int i;
00115      int k;
00116      int ntrk;
00117      int nhit;
00118      int stat;
00119 
00120      int lay;
00121      int cel;
00122      int wir;
00123      int lr;
00124      double dmeas;
00125      double doca;
00126      double resi;
00127 
00128      bool fgHitLay[MdcCalNLayer];
00129 
00130      MdcCalRecTrk* rectrk;
00131      MdcCalRecHit* rechit;
00132 
00133      ntrk = event->getNTrk();
00134      log << MSG::DEBUG << "number of tracks: " << ntrk << endreq;
00135 
00136      for(i=0; i<ntrk; i++){
00137           rectrk = event -> getRecTrk(i);
00138           nhit = rectrk -> getNHits();
00139 
00140           // dr cut
00141           double dr = rectrk->getDr();
00142           if(fabs(dr) > m_param.drCut) continue;
00143 
00144           // dz cut
00145           double dz = rectrk->getDz();
00146           if(fabs(dz) > m_param.dzCut) continue;
00147 
00148           for(lay=0; lay<MdcCalNLayer; lay++) fgHitLay[lay] = false;
00149           for(k=0; k<nhit; k++){
00150                rechit = rectrk -> getRecHit(k);
00151                lay = rechit -> getLayid();
00152                fgHitLay[lay] = true;
00153           }
00154 
00155           int nhitlay = 0;
00156           for(lay=0; lay<MdcCalNLayer; lay++) if(fgHitLay[lay]) nhitlay++;
00157           if(nhitlay < m_param.nHitLayCut) continue;
00158 
00159           log << MSG::DEBUG << "number of hits: " << nhit << endreq;
00160           for(k=0; k<nhit; k++){
00161                rechit = rectrk -> getRecHit(k);
00162                lay = rechit -> getLayid();
00163                cel = rechit -> getCellid();
00164                wir = m_mdcGeomSvc ->Wire(lay, cel)->Id();
00165                lr = rechit -> getLR();
00166                doca = rechit -> getDocaExc();
00167                dmeas = rechit -> getDmeas();
00168                resi = rechit -> getResiExc();
00169                stat = rechit -> getStat();
00170 
00171                if(1 != stat) continue;
00172 
00173                if( (fabs(doca) < m_docaMin[lay]) ||
00174                    (fabs(doca) > m_docaMax[lay]) ||
00175                    (fabs(resi) > m_param.maxResi) ){
00176                     continue;
00177                }
00178 
00179                if(0 == lay){
00180                     if( ! fgHitLay[1] ) continue;
00181                } else if(42 == lay){
00182                     if( ! fgHitLay[41] ) continue;
00183                } else{
00184                     if( (!fgHitLay[lay-1]) && (!fgHitLay[lay+1]) ) continue;
00185                }
00186 
00187                if(wir < (m_mdcGeomSvc->getWireSize()) ){
00188                     if( 0 == lr ){
00189                          m_hleft[wir] -> Fill(resi);
00190                     }else if( 1 == lr ){
00191                          m_hright[wir] -> Fill(resi);
00192                     }
00193                }else{
00194                     std::cout << "wir: " << wir << std::endl;
00195                }
00196           }
00197      }
00198      return 1;
00199 }

00045                                                              {
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 
00055      MdcCalib::initialize(m_hlist, m_mdcGeomSvc, m_mdcFunSvc);
00056 
00057      int i;
00058      int nwire;
00059      int lay;
00060      int cel;
00061      char hname[200];
00062 
00063      m_fdWire = new TFolder("WireCor", "WireCor");
00064      m_hlist->Add(m_fdWire);
00065 
00066      m_fdResiWire = new TFolder("ResiWire", "ResiWire");
00067      m_hlist->Add(m_fdResiWire);
00068 
00069      nwire = m_mdcGeomSvc -> getWireSize();
00070      for(i=0; i<nwire; i++){
00071           lay = m_mdcGeomSvc -> Wire(i) -> Layer();
00072           cel = m_mdcGeomSvc -> Wire(i) -> Cell();
00073 
00074           sprintf(hname, "h%04d_L%02d_%03d_Left", i, lay, cel);
00075           m_hleft[i] = new TH1F(hname, hname, 300, -1.5, 1.5);
00076           m_fdWire->Add(m_hleft[i]);
00077 
00078           sprintf(hname, "h%04d_L%02d_%03d_Right", i, lay, cel);
00079           m_hright[i] = new TH1F(hname, hname, 300, -1.5, 1.5);
00080           m_fdWire->Add(m_hright[i]);
00081      }
00082 
00083      m_hdwxtot = new TH1F("dwXtot", "", 100, -0.5, 0.5);
00084      m_fdResiWire->Add(m_hdwxtot);
00085 
00086      m_hddwx = new TH1F("ddwX", "", 100, -0.5, 0.5);
00087      m_fdResiWire->Add(m_hddwx);
00088 
00089      m_hdwytot = new TH1F("dwYtot", "", 100, -0.5, 0.5);
00090      m_fdResiWire->Add(m_hdwytot);
00091 
00092      m_hddwy = new TH1F("ddwY", "", 100, -0.5, 0.5);
00093      m_fdResiWire->Add(m_hddwy);
00094 
00095      m_hLrResiSum = new TH1F("LrResiSum", "", 200, -0.5, 0.5);
00096      m_fdResiWire->Add(m_hLrResiSum);
00097 
00098      m_hLrResiSub = new TH1F("LrResiSub", "", 200, -0.5, 0.5);
00099      m_fdResiWire->Add(m_hLrResiSub);
00100 }

00352                                                     {
00353      MdcCalib::read_file(path);
00354 }

00046                                                    {
00047      MdcCalib::setParam(param);
00048      m_param = param;
00049 }

00318                                        {
00319      MdcCalib::settuple(path);
00320      TFile *f1 = new TFile(path.c_str());
00321      std::cout<<path<<std::endl;
00322      //
00323      // m_mdcGeomSvc = mdcGeomSvc;
00324      //m_mdcFunSvc = mdcFunSvc;
00325      int i;
00326      int nwire;
00327      int lay;
00328      int cel;
00329      char hname[200];
00330      m_fdWire=(TFolder *)f1->Get("WireCor;1");
00331      m_fdResiWire=(TFolder *)f1->Get("ResiWire;1");
00332 
00333      nwire = m_mdcGeomSvc -> getWireSize();
00334      for(i=0; i<nwire; i++){
00335           lay = m_mdcGeomSvc -> Wire(i) -> Layer();
00336           cel = m_mdcGeomSvc -> Wire(i) -> Cell();
00337 
00338           sprintf(hname, "h%04d_L%02d_%03d_Left", i, lay, cel);
00339           m_hleft[i]=(TH1F *)m_fdWire->FindObject(hname);
00340           sprintf(hname, "h%04d_L%02d_%03d_Right", i, lay, cel);
00341           m_hright[i] =(TH1F *)m_fdWire->FindObject(hname);
00342      }
00343 
00344      m_hdwxtot = (TH1F *) m_fdResiWire->FindObject("dwXtot");
00345      m_hddwx =(TH1F *) m_fdResiWire->FindObject("ddwX");
00346      m_hdwytot =(TH1F *)m_fdResiWire->FindObject("dwYtot");
00347      m_hddwy =(TH1F *)m_fdResiWire->FindObject("ddwY");
00348      m_hLrResiSum =(TH1F *)m_fdResiWire->FindObject("LrResiSum");
00349      m_hLrResiSub = (TH1F *)m_fdResiWire->FindObject("LrResiSub");
00350 }

00201                                                   {
00202      IMessageSvc *msgSvc;
00203      Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
00204      MsgStream log(msgSvc, "WrMdcCalib");
00205      log << MSG::INFO << "WrMdcCalib::updateConst()" << endreq;
00206 
00207      MdcCalib::updateConst(calconst);
00208 
00209      int i;
00210      int lay;
00211      int cel;
00212      int nwire = m_mdcGeomSvc -> getWireSize();
00213      double dwphi;              // wire derivation in phi direction
00214      double resiLrSum;          // to fill histogram
00215      double wpos[6];
00216      double xx;
00217      double yy;
00218      double rr;
00219      double dx;
00220      double dy;
00221      double wphi;
00222 
00223      Stat_t entry_l;
00224      double mean_l;
00225      Stat_t entry_r;
00226      double mean_r;
00227      const MdcGeoWire* pWire;
00228 
00229      double ddw[MdcCalTotCell][6];
00230      double dwinput[MdcCalTotCell][6];
00231      string strtmp;
00232      ifstream fwpinput(m_param.wpcFile.c_str());
00233      for(i=0; i<7; i++) fwpinput >> strtmp;
00234      for(i=0; i<nwire; i++){
00235           fwpinput >> strtmp >> dwinput[i][0] >> dwinput[i][1] >> dwinput[i][2]
00236                    >> dwinput[i][3] >> dwinput[i][4] >> dwinput[i][5];
00237      }
00238      fwpinput.close();
00239 
00240      std::cout << "totwire: " << nwire << std::endl;
00241      for(i=0; i<nwire; i++){
00242           pWire = m_mdcGeomSvc -> Wire(i);
00243           lay = pWire -> Layer();
00244           cel = pWire -> Cell();
00245 
00246           if(1 == m_param.fgCalib[lay]){
00247                entry_l = m_hleft[i] -> GetEntries();
00248                mean_l = m_hleft[i] -> GetMean();
00249 
00250                entry_r = m_hright[i] -> GetEntries();
00251                mean_r = m_hright[i] -> GetMean();
00252 
00253                dwphi = 0.5 * (mean_l - mean_r);
00254           } else{
00255                dwphi = 0.0;
00256           }
00257 
00258           resiLrSum = 0.5 * (mean_l + mean_r);
00259           m_hLrResiSum->Fill(resiLrSum);
00260           m_hLrResiSub->Fill(dwphi);
00261 
00262           wpos[0] = pWire -> Backward().x(); // east end
00263           wpos[1] = pWire -> Backward().y();
00264           wpos[3] = pWire -> Forward().x(); // west end
00265           wpos[4] = pWire -> Forward().y();
00266 
00267           xx = 0.5 * (wpos[0] + wpos[3]);
00268           yy = 0.5 * (wpos[1] + wpos[4]);
00269           rr = sqrt( (xx * xx) + (yy * yy) );
00270 
00271           if( yy >= 0 )  wphi = acos(xx / rr);
00272           else           wphi = PI2 - acos(xx / rr);
00273 
00274           dx = -1.0 * dwphi * sin(wphi);
00275           dy = dwphi * cos(wphi);
00276 
00277           ddw[i][0] = dx;
00278           ddw[i][1] = dy;
00279           ddw[i][2] = 0.0;
00280           ddw[i][3] = dx;
00281           ddw[i][4] = dy;
00282           ddw[i][5] = 0.0;
00283      }
00284 
00285      ofstream fwpc("WirePosCalib_new.txt");
00286      fwpc << setw(5) << "wireId" << setw(15) << "dx_East(mm)"
00287           << setw(15) << "dy_East(mm)" << setw(15) << "dz_East(mm)"
00288           << setw(15) << "dx_West(mm)" << setw(15) << "dy_West(mm)"
00289           << setw(15) << "dz_West(mm)" << endl;
00290 
00291      ofstream fdw("dw.txt");
00292      fdw << setw(5) << "wireId" << setw(15) << "ddx_East(mm)"
00293          << setw(15) << "ddy_East(mm)" << setw(15) << "ddz_East(mm)"
00294          << setw(15) << "ddx_West(mm)" << setw(15) << "ddy_West(mm)"
00295          << setw(15) << "ddz_West(mm)" << endl;
00296 
00297      int k;
00298      double dwtot[6];
00299      for(i=0; i<nwire; i++){
00300           for(k=0; k<6; k++) dwtot[k] = dwinput[i][k] + ddw[i][k];
00301           fwpc << setw(5) << i;
00302           for(k=0; k<6; k++) fwpc << setw(15) << dwtot[k];
00303           fwpc << endl;
00304 
00305           fdw << setw(5) << i;
00306           for(k=0; k<6; k++) fdw << setw(15) << ddw[i][k];
00307           fdw << endl;
00308           m_hdwxtot->Fill(dwtot[0]);
00309           m_hddwx->Fill(ddw[i][0]);
00310           m_hdwytot->Fill(dwtot[1]);
00311           m_hddwy->Fill(ddw[i][1]);
00312      }
00313      fwpc.close();
00314      fdw.close();
00315      return 1;
00316 }

WrMdcCalib Class Reference

Public Member Functions

Private Attributes

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation


Public Member Functions
void	clear ()
void	clear ()
int	fillHist (MdcCalEvent *event)
int	fillHist (MdcCalEvent *event)
void	initialize (TObjArray hlist, IMdcGeomSvc mdcGeomSvc, IMdcCalibFunSvc *mdcFunSvc)
void	initialize (TObjArray hlist, IMdcGeomSvc mdcGeomSvc, IMdcCalibFunSvc *mdcFunSvc)
virtual void	read_file (std::vector< std::string > path)
virtual void	read_file (std::vector< std::string > path)
void	setParam (MdcCalParams &param)
void	setParam (MdcCalParams &param)
virtual void	settuple (std::string path)
virtual void	settuple (std::string path)
int	updateConst (MdcCalibConst *calconst)
int	updateConst (MdcCalibConst *calconst)
	WrMdcCalib ()
	WrMdcCalib ()
	~WrMdcCalib ()
	~WrMdcCalib ()
Private Attributes
double	m_docaMax [MdcCalNLayer]
double	m_docaMin [MdcCalNLayer]
TFolder *	m_fdResiWire
TFolder *	m_fdResiWire
TFolder *	m_fdWire
TFolder *	m_fdWire
TH1F *	m_hddwx
TH1F *	m_hddwx
TH1F *	m_hddwy
TH1F *	m_hddwy
TH1F *	m_hdwxtot
TH1F *	m_hdwxtot
TH1F *	m_hdwytot
TH1F *	m_hdwytot
TH1F *	m_hleft [MdcCalTotCell]
TH1F *	m_hleft [MdcCalTotCell]
TObjArray *	m_hlist
TObjArray *	m_hlist
TH1F *	m_hLrResiSub
TH1F *	m_hLrResiSub
TH1F *	m_hLrResiSum
TH1F *	m_hLrResiSum
TH1F *	m_hright [MdcCalTotCell]
TH1F *	m_hright [MdcCalTotCell]
IMdcCalibFunSvc *	m_mdcFunSvc
IMdcCalibFunSvc *	m_mdcFunSvc
IMdcGeomSvc *	m_mdcGeomSvc
IMdcGeomSvc *	m_mdcGeomSvc
MdcCalParams	m_param