#include <MdcPatTuning.h>
Public Member Functions | |
StatusCode | execute () |
StatusCode | execute () |
StatusCode | finalize () |
StatusCode | finalize () |
StatusCode | initialize () |
StatusCode | initialize () |
MdcPatTuning (const std::string &name, ISvcLocator *pSvcLocator) | |
MdcPatTuning (const std::string &name, ISvcLocator *pSvcLocator) | |
Private Member Functions | |
StatusCode | bookNTuple () |
StatusCode | bookNTuple () |
StatusCode | createHitMap (void) |
StatusCode | createHitMap (void) |
int | createSegs (void) |
int | createSegs (void) |
StatusCode | fillInit () |
StatusCode | fillInit () |
void | fillWithStSegs (const RecMdcTrack *tk) |
void | fillWithStSegs (const RecMdcTrack *tk) |
void | iniTrackFinding (void) |
void | iniTrackFinding (void) |
Private Attributes | |
double | Bz |
NTuple::Item< double > | m_1_chisq |
NTuple::Item< double > | m_1_chisq |
NTuple::Item< double > | m_1_ct |
NTuple::Item< double > | m_1_ct |
NTuple::Item< long > | m_1_nHit |
NTuple::Item< long > | m_1_nHit |
NTuple::Item< long > | m_1_nPass |
NTuple::Item< long > | m_1_nPass |
NTuple::Item< long > | m_1_nSlPass |
NTuple::Item< long > | m_1_nSlPass |
NTuple::Item< double > | m_1_phi |
NTuple::Item< double > | m_1_phi |
NTuple::Item< long > | m_1_quality |
NTuple::Item< long > | m_1_quality |
NTuple::Item< long > | m_1_slayer |
NTuple::Item< long > | m_1_slayer |
NTuple::Item< double > | m_1_slope |
NTuple::Item< double > | m_1_slope |
NTuple::Item< long > | m_1_stat |
NTuple::Item< long > | m_1_stat |
NTuple::Item< double > | m_1_z0 |
NTuple::Item< double > | m_1_z0 |
NTuple::Item< double > | m_1c_maxPhiA |
NTuple::Item< double > | m_1c_maxPhiA |
NTuple::Item< double > | m_1c_minPhiA |
NTuple::Item< double > | m_1c_minPhiA |
NTuple::Item< double > | m_1c_phiArg |
NTuple::Item< double > | m_1c_phiArg |
NTuple::Item< double > | m_1c_phiAx |
NTuple::Item< double > | m_1c_phiAx |
bool | m_debug |
bool | m_dropHot |
MdcFlagHold | m_flags |
const MdcDetector * | m_gm |
const MdcDetector * | m_gm |
std::auto_ptr< MdcSegData > | m_hitData |
std::auto_ptr< MdcSegData > | m_hitData |
bool | m_keepBadTdc |
bool | m_keepUnmatch |
int | m_maxMdcDigi |
MdcCalibFunSvc * | m_mdcCalibFunSvc |
MdcCalibFunSvc * | m_mdcCalibFunSvc |
MdcCheckUtil * | m_mdcCheckUtil |
MdcCheckUtil * | m_mdcCheckUtil |
MdcDigiVec | m_mdcDigiVec |
MdcHitCol * | m_mdcHitCol |
MdcHitCol * | m_mdcHitCol |
MdcHitMap * | m_mdcHitMap |
MdcHitMap * | m_mdcHitMap |
bool | m_onlyUnusedHits |
std::string | m_paramFile |
IMagneticFieldSvc * | m_pIMF |
IMagneticFieldSvc * | m_pIMF |
bool | m_poisonHits |
RawDataProviderSvc * | m_rawDataProviderSvc |
RawDataProviderSvc * | m_rawDataProviderSvc |
std::auto_ptr< MdcSegFinder > | m_segFinder |
std::auto_ptr< MdcSegFinder > | m_segFinder |
std::auto_ptr< MdcSegList > | m_segs |
std::auto_ptr< MdcSegList > | m_segs |
std::auto_ptr< MdcTrackListBase > | m_tracks |
std::auto_ptr< MdcTrackListBase > | m_tracks |
NTuple::Tuple * | m_tuple |
NTuple::Tuple * | m_tuple |
NTuple::Tuple * | m_tuple1 |
NTuple::Tuple * | m_tuple1 |
long | t_eventNo |
int | t_nTk |
double | t_recTkNum |
long | t_runNo |
double | t_t0 |
double | t_t0Stat |
double | t_t0Truth |
|
00049 : 00050 Algorithm(name, pSvcLocator) { 00051 declareProperty("maxMdcDigi", m_maxMdcDigi= 0); 00052 declareProperty("keepBadTdc", m_keepBadTdc= 0); 00053 declareProperty("dropHot", m_dropHot= 0); 00054 declareProperty("keepUnmatch", m_keepUnmatch= 0); 00055 00056 declareProperty("paramFile", m_paramFile = "params"); 00057 declareProperty("OnlyUnusedHits",m_onlyUnusedHits = 0); 00058 declareProperty("poisonHits", m_poisonHits= 0); 00059 00060 declareProperty("debug", m_debug= 0); 00061 }
|
|
|
|
|
|
00147 { 00148 MsgStream log(msgSvc(), name()); 00149 StatusCode sc = StatusCode::SUCCESS; 00150 00151 NTuplePtr nt3(ntupleSvc(), "MdcPatTuning/digi"); 00152 if ( nt3 ) { m_tuple = nt3;} 00153 else { 00154 m_tuple = ntupleSvc()->book ("MdcPatTuning/digi", CLID_ColumnWiseTuple, "event"); 00155 if ( m_tuple ) { 00156 sc = m_tuple->addItem ("slayer", m_1_slayer); 00157 sc = m_tuple->addItem ("slope", m_1_slope); 00158 sc = m_tuple->addItem ("phi", m_1_phi); 00159 sc = m_tuple->addItem ("chisq", m_1_chisq); 00160 sc = m_tuple->addItem ("quality", m_1_quality); 00161 sc = m_tuple->addItem ("ct", m_1_ct); 00162 sc = m_tuple->addItem ("z0", m_1_z0); 00163 sc = m_tuple->addItem ("stat", m_1_stat); 00164 sc = m_tuple->addItem ("nSlPass", m_1_nSlPass); 00165 sc = m_tuple->addItem ("nHit", m_1_nHit); 00166 sc = m_tuple->addItem ("nPass", m_1_nPass); 00167 } 00168 } 00169 00170 NTuplePtr nt1(ntupleSvc(), "MdcPatTuning/seg"); 00171 if ( nt1 ) { m_tuple1 = nt1;} 00172 else { 00173 m_tuple1 = ntupleSvc()->book ("MdcPatTuning/seg", CLID_ColumnWiseTuple, "event"); 00174 if ( m_tuple1 ) { 00175 sc = m_tuple->addItem ("phiArg", m_1c_phiArg); 00176 sc = m_tuple->addItem ("phiAx", m_1c_phiAx); 00177 sc = m_tuple->addItem ("maxPhiA", m_1c_maxPhiA); 00178 sc = m_tuple->addItem ("minPhiA", m_1c_minPhiA); 00179 } 00180 } 00181 return StatusCode::SUCCESS; 00182 }
|
|
|
|
00301 { 00302 MsgStream log(msgSvc(), name()); 00303 StatusCode sc; 00304 00305 DataObject *pnode = 0; 00306 sc = eventSvc()->retrieveObject("/Event/Hit",pnode); 00307 if(!sc.isSuccess()) { 00308 pnode = new DataObject; 00309 sc = eventSvc()->registerObject("/Event/Hit",pnode); 00310 if(!sc.isSuccess()) { 00311 log << MSG::FATAL << " Could not register hit branch" <<endreq; 00312 return StatusCode::FAILURE; 00313 } 00314 } 00315 00316 IDataManagerSvc *dataManSvc = dynamic_cast<IDataManagerSvc*> (eventSvc()); 00317 DataObject *hitCol; 00318 eventSvc()->findObject("/Event/Hit/MdcHitCol",hitCol); 00319 if(hitCol!= NULL) { 00320 dataManSvc->clearSubTree("/Event/Hit/MdcHitCol"); 00321 eventSvc()->unregisterObject("/Event/Hit/MdcHitCol"); 00322 } 00323 m_mdcHitCol = new MdcHitCol; 00324 sc = eventSvc()->registerObject("/Event/Hit/MdcHitCol",m_mdcHitCol); 00325 if(!sc.isSuccess()) { 00326 log << MSG::FATAL << " Could not register hit collection" <<endreq; 00327 return StatusCode::FAILURE; 00328 } 00329 00330 DataObject *hitMap; 00331 eventSvc()->findObject("/Event/Hit/MdcHitMap",hitMap); 00332 if(hitMap!= NULL) { 00333 dataManSvc->clearSubTree("/Event/Hit/MdcHitMap"); 00334 eventSvc()->unregisterObject("/Event/Hit/MdcHitMap"); 00335 } 00336 m_mdcHitMap = new MdcHitMap(*m_gm); 00337 sc = eventSvc()->registerObject("/Event/Hit/MdcHitMap",m_mdcHitMap); 00338 if(!sc.isSuccess()) { 00339 log << MSG::FATAL << " Could not register hit map" <<endreq; 00340 return StatusCode::FAILURE; 00341 } 00342 00343 // retrieve Mdc digi vector form RawDataProviderSvc 00344 uint32_t getDigiFlag = 0; 00345 getDigiFlag += m_maxMdcDigi; 00346 if(m_dropHot) getDigiFlag |= MdcRawDataProvider::b_dropHot; 00347 if(m_keepBadTdc) getDigiFlag |= MdcRawDataProvider::b_keepBadTdc; 00348 if(m_keepUnmatch) getDigiFlag |= MdcRawDataProvider::b_keepUnmatch; 00349 m_mdcDigiVec = m_rawDataProviderSvc->getMdcDigiVec(getDigiFlag); 00350 00351 00352 MdcDigiVec::iterator iter = m_mdcDigiVec.begin(); 00353 for (;iter != m_mdcDigiVec.end(); iter++ ) { 00354 const MdcDigi* aDigi = *iter; 00355 MdcHit *hit = new MdcHit(aDigi, m_gm); 00356 hit->setCalibSvc(m_mdcCalibFunSvc); 00357 hit->setCosmicFit(false); 00358 hit->setCountPropTime(true); 00359 m_mdcHitCol->push_back(hit); 00360 m_mdcHitMap->addHit(*hit); 00361 } 00362 00363 return StatusCode::SUCCESS; 00364 00365 }
|
|
|
|
00367 { 00368 //bool useAllAmbig = true; 00369 00370 m_hitData->loadevent(m_mdcHitCol, m_mdcHitMap, t_t0); 00371 00372 int nSegs = m_segFinder->createSegs(m_gm, *m_segs, m_hitData->segUsage(), 00373 m_hitData->hitMap(), t_t0); 00374 00375 if(m_debug) std::cout<<__FILE__<<" "<<__LINE__<<nSegs<< std::endl; 00376 //m_segs->plot(); 00377 00378 return nSegs; 00379 }
|
|
|
|
00106 { 00107 MsgStream log(msgSvc(), name()); 00108 StatusCode sc = StatusCode::SUCCESS; 00109 00110 fillInit(); 00111 00112 00113 // Get RecMdcTrackCol and RecMdcHitCol 00114 SmartDataPtr<RecMdcTrackCol> recMdcTrackCol(eventSvc(), "/Event/Recon/RecMdcTrackCol"); 00115 if (!recMdcTrackCol){ 00116 log << MSG::WARNING<< " Unable to retreve recMdcTrackCol" << endreq; 00117 return StatusCode::SUCCESS; 00118 } 00119 00120 if(StatusCode::SUCCESS != createHitMap()) return StatusCode::FAILURE; 00121 int nSeg = createSegs(); 00122 00123 if(m_debug) std::cout<<__FILE__<<" createSegs find "<<nSeg <<" segs"<< std::endl; 00124 00125 RecMdcTrackCol::iterator it = recMdcTrackCol->begin(); 00126 for(; it != recMdcTrackCol->end(); it++ ) { 00127 fillWithStSegs(*it); 00128 } 00129 00130 00131 m_segs->destroySegs(); 00132 return StatusCode::SUCCESS; 00133 }
|
|
|
|
00184 { 00185 MsgStream log(msgSvc(), name()); 00186 StatusCode sc = StatusCode::SUCCESS; 00187 00188 SmartDataPtr<Event::EventHeader> evtHead(eventSvc(),"/Event/EventHeader"); 00189 if (evtHead) { 00190 t_runNo = evtHead->runNumber(); 00191 t_eventNo = evtHead->eventNumber(); 00192 }else{ 00193 std::cout<<"Could not retrieve event header "<<std::endl; 00194 } 00195 std::cout<< "evtNo:"<<t_eventNo << std::endl;//yzhang debug 00196 00197 //Get event time 00198 t_t0 = -1; 00199 t_t0Stat = -1; 00200 SmartDataPtr<RecEsTimeCol> aevtimeCol(eventSvc(),"/Event/Recon/RecEsTimeCol"); 00201 if (aevtimeCol) { 00202 RecEsTimeCol::iterator iter_evt = aevtimeCol->begin(); 00203 for (;iter_evt!=aevtimeCol->end();iter_evt++){ 00204 t_t0 = (*iter_evt)->getTest(); 00205 t_t0Stat = (*iter_evt)->getStat(); 00206 } 00207 }else{ 00208 log << MSG::WARNING << "Could not find RecEsTimeCol" << endreq; 00209 } 00210 00211 return sc; 00212 }
|
|
|
|
00215 { 00216 00217 //Calc. track passed cell 00218 bool passed[43][288]={false}; 00219 int nSlPass[11]={0,0,0,0,0,0,0,0,0,0,0}; 00220 m_mdcCheckUtil->cellPassed(tk, passed); 00221 int ipass=0; 00222 for(int i=0;i<43;i++){ 00223 for(int j=0; j<288; j++){ 00224 if(passed[i][j]) { 00225 nSlPass[i%4]++; 00226 ipass++; 00227 } 00228 } 00229 } 00230 if(m_debug) std::cout<<__FILE__<<" tk: "<<tk->trackId()<<" npass "<< ipass<< std::endl; 00231 00232 00233 //New TrkRecoTrk and TrkExchangePar 00234 TrkExchangePar par(tk->helix(0),tk->helix(1),tk->helix(2),tk->helix(3),tk->helix(4)); 00235 00236 //Calc. seg info with track 00237 HepVector helix= m_mdcCheckUtil->besPar2PatPar(tk->helix()); 00238 double d0 = helix[0]; 00239 double phi0 = helix[1]; 00240 double kap = helix[2]; 00241 00242 for (int slayer=0; slayer<11; slayer++){ 00243 const MdcSuperLayer *slay = m_gm->SuperLayer(slayer); 00244 MdcSegWorks segStuff;// holds some calculated values to pass to segInfo 00245 00246 double radius = slay->rEnd(); 00247 double rCurl = 99999999.; 00248 if (radius >= rCurl) return; 00249 bool lStraight = false; 00250 if (fabs(kap) < 1.e-9) { lStraight = true; } 00251 00252 double phiArg = kap * radius; 00253 if (lStraight) phiArg = d0 / radius; 00254 if (phiArg < -1.0) phiArg = -1.0; 00255 else if (phiArg > 1.0) phiArg = 1.0; 00256 00257 segStuff.phiArg = phiArg; 00258 segStuff.phiAx.setRad(phi0 + asin(phiArg)); // Approx??!!!! 00259 BesAngle delPhi( fabs(slay->delPhi()) ); 00260 00261 BesAngle maxPhiA = segStuff.phiAx + delPhi + 0.05; // allow a little slop 00262 BesAngle minPhiA = segStuff.phiAx - delPhi - 0.05; 00263 00264 //double maxPhi = maxPhiA; 00265 segStuff.wirLen2inv = 1./ (slay->zEnd() * slay->zEnd()); 00266 00267 MdcSegInfoSterO *info = new MdcSegInfoSterO; 00268 00269 m_1c_phiArg = phiArg; 00270 m_1c_phiAx = segStuff.phiAx; 00271 m_1c_maxPhiA = maxPhiA; 00272 m_1c_minPhiA = minPhiA; 00273 m_tuple1->write(); 00274 00275 MdcSeg* inSeg = (MdcSeg *) m_segs->oneList(slayer)->first(); 00276 while (inSeg != 0) { 00277 if(m_debug) std::cout<<__FILE__<<" "<<__LINE__<< std::endl; 00278 //inSeg->plotSeg(); 00279 int isBad = info->calcStereo(inSeg, par, segStuff, Bz); 00280 m_1_slayer = slayer; 00281 m_1_slope = inSeg->slope(); 00282 m_1_phi = inSeg->phi(); 00283 m_1_chisq = inSeg->chisq(); 00284 m_1_quality = inSeg->quality(); 00285 m_1_ct = info->ct(); 00286 m_1_z0 = info->z0(); 00287 m_1_stat = isBad; 00288 m_1_nSlPass = nSlPass[slayer]; 00289 m_1_nHit = inSeg->nHit(); 00290 int nPass = 0; 00291 for(int i=0; i<inSeg->nHit(); i++){ 00292 int l = inSeg->hit(i)->mdcHit()->layernumber(); 00293 int w = inSeg->hit(i)->mdcHit()->wirenumber(); 00294 if(passed[l][w]) nPass++; 00295 } 00296 m_1_nPass = nPass; 00297 }//seg 00298 }//superlayer 00299 }//fillWithStSegs
|
|
|
|
00136 { 00137 00138 m_segs.reset(0); 00139 m_segFinder.reset(0); 00140 m_hitData.reset(0); 00141 m_tracks.reset(0); 00142 00143 return StatusCode::SUCCESS; 00144 }
|
|
|
|
00065 { 00066 MsgStream log(msgSvc(), name()); 00067 StatusCode sc; 00068 00069 t_nTk = 0; 00070 //Initailize magnetic filed 00071 sc = service ("MagneticFieldSvc",m_pIMF); 00072 if(sc!=StatusCode::SUCCESS) { 00073 log << MSG::ERROR << "Unable to open Magnetic field service"<<endreq; 00074 } 00075 Bz = m_pIMF->getReferField()*1000.; 00076 00077 IMdcCalibFunSvc* imdcCalibSvc; 00078 sc = service ("MdcCalibFunSvc",imdcCalibSvc); 00079 m_mdcCalibFunSvc = dynamic_cast<MdcCalibFunSvc*> (imdcCalibSvc); 00080 if ( sc.isFailure() ){ 00081 log << MSG::FATAL << "Could not load MdcCalibFunSvc!" << endreq; 00082 } 00083 00084 IRawDataProviderSvc* irawDataProviderSvc; 00085 sc = service ("RawDataProviderSvc", irawDataProviderSvc); 00086 m_rawDataProviderSvc = dynamic_cast<RawDataProviderSvc*> (irawDataProviderSvc); 00087 if ( sc.isFailure() ){ 00088 log << MSG::FATAL << "Could not load RawDataProviderSvc!" << endreq; 00089 return StatusCode::FAILURE; 00090 } 00091 00092 if (StatusCode::SUCCESS != bookNTuple()) { 00093 log << MSG::WARNING << " Could not book NTuple" << endreq; 00094 } 00095 00096 m_gm = MdcDetector::instance(false);//do not do sag 00097 00098 m_mdcCheckUtil = new MdcCheckUtil(false); 00099 00100 iniTrackFinding(); 00101 00102 return StatusCode::SUCCESS; 00103 }
|
|
|
|
00382 { 00383 00384 //Initialize hitdata 00385 m_hitData.reset( new MdcSegData( m_poisonHits )); 00386 00387 //Initialize segFinder 00388 if (m_flags.findSegs) { 00389 m_segFinder.reset( new MdcSegFinder(m_flags.segPar.useAllAmbig) ); 00390 } 00391 //Initialize segList 00392 m_segs.reset( new MdcSegList(m_gm->nSuper(), m_flags.segPar) ); 00393 00394 //Initialize trackList 00395 m_tracks.reset( new MdcTrackList(m_flags.tkParTight) ); 00396 00397 00398 }
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|